374 lines
18 KiB
Python
374 lines
18 KiB
Python
import numpy as np
|
|
import pandas as pd
|
|
import itertools
|
|
from main.astrodatagui.db.StarsDB import StarDB
|
|
import matplotlib.ticker as tck
|
|
from matplotlib.pyplot import MaxNLocator
|
|
import matplotlib.pyplot as plt
|
|
from datetime import datetime
|
|
from errno import EEXIST
|
|
from os import makedirs, path
|
|
import shutil
|
|
|
|
def normalizePhase(phase, phaseMin = None, phaseMax = None):
|
|
if(phaseMin is None):
|
|
phaseMin = np.abs(np.min(phase))
|
|
if(phaseMax is None):
|
|
phaseMax = np.abs(np.max(phase))
|
|
return (phase + phaseMin) / (phaseMin + phaseMax) * (2)
|
|
|
|
def mkdir_p(mypath):
|
|
'''Creates a directory. equivalent to using mkdir -p on the command line'''
|
|
|
|
try:
|
|
makedirs(mypath)
|
|
except OSError as exc: # Python >2.5
|
|
if exc.errno == EEXIST and path.isdir(mypath):
|
|
pass
|
|
else: raise
|
|
|
|
fileName = "data.cff"
|
|
data = pd.read_pickle(fileName)
|
|
|
|
binList = [10, 20, 30]
|
|
spType = ["M", "K", "G", "F"]
|
|
|
|
useKepler = True
|
|
useK2 = True
|
|
useTESS = True
|
|
|
|
showSourceFilter = np.full(len(data), False)
|
|
if(useKepler):
|
|
showKepler = data["Source"] == "Kepler"
|
|
showSourceFilter |= showKepler
|
|
if(useK2):
|
|
showK2 = data["Source"] == "K2"
|
|
showSourceFilter |= showK2
|
|
if(useTESS):
|
|
showTESS = data["Source"] == "TESS"
|
|
showSourceFilter |= showTESS
|
|
|
|
current = datetime.now()
|
|
date = f"{current.year}-{current.month}-{current.day}"
|
|
time = f"{current.hour}-{current.minute}-{current.second}"
|
|
folderPath = f"../{date}/"
|
|
#folderPath = f"G:/Meine Ablage/Masterthesis/{date}/"
|
|
mkdir_p(folderPath)
|
|
|
|
for comboLength in range(1, len(spType) + 1):
|
|
for combo in itertools.combinations(spType, comboLength):
|
|
finalData = pd.DataFrame()
|
|
if("M" in combo):
|
|
Mfilter = data["SpType"].str.startswith("M")
|
|
Mfilter &= showSourceFilter
|
|
finalData = pd.concat([finalData, data[Mfilter]], ignore_index=True)
|
|
if("K" in combo):
|
|
Kfilter = data["SpType"].str.startswith("K")
|
|
Kfilter &= showSourceFilter
|
|
finalData = pd.concat([finalData, data[Kfilter]], ignore_index=True)
|
|
if("G" in combo):
|
|
Gfilter = data["SpType"].str.startswith("G")
|
|
Gfilter &= showSourceFilter
|
|
finalData = pd.concat([finalData, data[Gfilter]], ignore_index=True)
|
|
if("F" in combo):
|
|
Ffilter = data["SpType"].str.startswith("F")
|
|
Ffilter &= showSourceFilter
|
|
finalData = pd.concat([finalData, data[Ffilter]], ignore_index=True)
|
|
|
|
|
|
|
|
numStars = len(set(finalData["StarName"]))
|
|
|
|
pdcsapbinningData = []
|
|
locFolder = f"{folderPath}/{''.join(combo)}/"
|
|
mkdir_p(f"{locFolder}/")
|
|
csvFile = open(f"{locFolder}/{''.join(combo)}.csv", "a")
|
|
csvFile.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Normalized Phase of Peak,Peak in Period")
|
|
csvFile.write("\n")
|
|
for ind, row in finalData.reset_index().iterrows():
|
|
PDCSAPminOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][0]
|
|
PDCSAPmaxOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][1]
|
|
if(len(row["pdcsapFoldedPeaksPhasePair"]) > 0):
|
|
pdcsapVals = pd.DataFrame(row["pdcsapFoldedPeaksPhasePair"])
|
|
for td, peak, pv in zip(pdcsapVals["Phase"], pdcsapVals["Peak"], row["pdcsapPeaks"]):
|
|
normPhase = normalizePhase(td.value, np.abs(PDCSAPminOrigPhase), np.abs(PDCSAPmaxOrigPhase))
|
|
csvFile.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{normPhase},{peak['FlarePeak']}")
|
|
csvFile.write("\n")
|
|
pdcsapbinningData.append({"SpType": row["SpType"][0],
|
|
"PDCSAPNormPhase": normPhase,
|
|
"Peak": peak["FlarePeak"]})
|
|
if(peak["FlarePeak"] > 100):
|
|
print(row["StarName"], "has over 100 peak")
|
|
|
|
csvFile.close()
|
|
pdcsapbinningData = pd.DataFrame(pdcsapbinningData)
|
|
PDCSAPdataList = []
|
|
PDCSAPlabelList = []
|
|
PDCSAPcolorList = []
|
|
PDCSAPdataList2dhistPhase = []
|
|
PDCSAPdataList2dhistPeak = []
|
|
PDCSAPlabelList2dhist = []
|
|
PDCSAPcolorList2dhist = []
|
|
if("M" in combo):
|
|
Mfilter = pdcsapbinningData["SpType"] == "M"
|
|
PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"])
|
|
PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Mfilter]["Peak"])
|
|
PDCSAPlabelList2dhist.append("M Stars")
|
|
PDCSAPcolorList2dhist.append("red")
|
|
|
|
PDCSAPdataList.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"])
|
|
PDCSAPlabelList.append("M Stars")
|
|
PDCSAPcolorList.append("red")
|
|
if("K" in combo):
|
|
Kfilter = pdcsapbinningData["SpType"] == "K"
|
|
PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"])
|
|
PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Kfilter]["Peak"])
|
|
PDCSAPlabelList2dhist.append("K Stars")
|
|
PDCSAPcolorList2dhist.append("orange")
|
|
|
|
PDCSAPdataList.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"])
|
|
PDCSAPlabelList.append("K Stars")
|
|
PDCSAPcolorList.append("orange")
|
|
if("G" in combo):
|
|
Gfilter = pdcsapbinningData["SpType"] == "G"
|
|
PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"])
|
|
PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Gfilter]["Peak"])
|
|
PDCSAPlabelList2dhist.append("G Stars")
|
|
PDCSAPcolorList2dhist.append("yellow")
|
|
|
|
PDCSAPdataList.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"])
|
|
PDCSAPlabelList.append("G Stars")
|
|
PDCSAPcolorList.append("yellow")
|
|
if("F" in combo):
|
|
Ffilter = pdcsapbinningData["SpType"] == "F"
|
|
PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"])
|
|
PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Ffilter]["Peak"])
|
|
PDCSAPlabelList2dhist.append("F Stars")
|
|
PDCSAPcolorList2dhist.append("greenyellow")
|
|
|
|
PDCSAPdataList.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"])
|
|
PDCSAPlabelList.append("F Stars")
|
|
PDCSAPcolorList.append("greenyellow")
|
|
|
|
xData = pd.DataFrame()
|
|
yData = pd.DataFrame()
|
|
for aX, aY in zip(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak):
|
|
xData = pd.concat([xData, aX], ignore_index=True)
|
|
yData = pd.concat([yData, aY], ignore_index=True)
|
|
|
|
xData = np.asarray(xData.values)[:,0]
|
|
yData = np.asarray(yData.values)[:,0]
|
|
|
|
for bins in binList:
|
|
figHisto, ((axHisto)) = plt.subplots(nrows=1, ncols=1)
|
|
y, binEdges, _ = axHisto.hist(PDCSAPdataList, bins,
|
|
label=PDCSAPlabelList,
|
|
color=PDCSAPcolorList,
|
|
stacked=True,
|
|
range=[0, 2])
|
|
bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
|
|
if(isinstance(y[0], np.ndarray)):
|
|
y = y[-1]
|
|
n_i = y
|
|
m_i = bincenters * np.pi
|
|
N = np.sum(n_i)
|
|
mean = np.sum(n_i * m_i)/N
|
|
stdDev = np.sqrt(np.sum(((n_i - mean)**2)) / (N-1))
|
|
menStd = np.sqrt(y)
|
|
axHisto.bar(bincenters[y > 0], y[y > 0], width=0, color='r', yerr=stdDev)
|
|
axHisto.set_ylim(0, max(y) + stdDev)
|
|
axHisto.set_ylabel("Num. flares")
|
|
axHisto.set_xlabel("Phase")
|
|
axHisto.set_title(f"Flare count per phase of {', '.join(combo)} type stars with {bins} bins ({numStars} stars)")
|
|
axHisto.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axHisto.xaxis.set_major_locator(MaxNLocator(5))
|
|
axHisto.legend()
|
|
|
|
axHistoPhase = axHisto.twinx()
|
|
secAxisXdata = np.linspace(0, 2, num=10000)
|
|
secAxisYdata = np.cos(secAxisXdata*np.pi) + 1
|
|
axHistoPhase.plot(secAxisXdata, secAxisYdata)
|
|
axHistoPhase.set_ylim(0, 7)
|
|
|
|
plt.savefig(f"{locFolder}/{''.join(combo)}-Flarecount-{bins}_Bins.png")
|
|
plt.close()
|
|
|
|
for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
|
|
figFlarepeakHist, ((axFlarepeakHist)) = plt.subplots(nrows=1, ncols=1)
|
|
axFlarepeakHist.set_ylabel("Flare peak")
|
|
axFlarepeakHist.set_xlabel("Phase")
|
|
axFlarepeakHist.set_title(f"Flare peak per phase histogram of {', '.join(combo)} type stars with {bins} bins ({numStars} stars)")
|
|
H, xedges, yedges = np.histogram2d(xData, yData, bins=bins, range=[[0, 2], [0.99, maxY]])
|
|
cmax = 11
|
|
H_clipped = np.clip(H, None, cmax)
|
|
im = axFlarepeakHist.imshow(H_clipped.T, origin='lower', interpolation='nearest',
|
|
extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]],
|
|
aspect='auto', cmap='viridis')
|
|
figFlarepeakHist.colorbar(im, label='Counts', ax=axFlarepeakHist)
|
|
axFlarepeakHist.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axFlarepeakHist.xaxis.set_major_locator(MaxNLocator(5))
|
|
plt.savefig(f"{locFolder}/{''.join(combo)}-Flarepeaks-{bins}_Bins_maxY-{maxY}.png")
|
|
plt.close()
|
|
|
|
for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
|
|
figFlarePeaks, ((axFlarePeaks)) = plt.subplots(nrows=1, ncols=1)
|
|
axFlarePeaks.set_title(f"Flare peaks per phase of {', '.join(combo)} type stars ({numStars} stars)")
|
|
if("M" in combo):
|
|
axFlarePeaks.scatter(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"],
|
|
pdcsapbinningData[Mfilter]["Peak"],
|
|
label="M Stars", color="red")
|
|
if("K" in combo):
|
|
axFlarePeaks.scatter(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"],
|
|
pdcsapbinningData[Kfilter]["Peak"],
|
|
label="K Stars", color="orange")
|
|
if("G" in combo):
|
|
axFlarePeaks.scatter(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"],
|
|
pdcsapbinningData[Gfilter]["Peak"],
|
|
label="G Stars", color="yellow")
|
|
if("F" in combo):
|
|
axFlarePeaks.scatter(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"],
|
|
pdcsapbinningData[Ffilter]["Peak"],
|
|
label="F Stars", color="greenyellow")
|
|
axFlarePeaks.set_xlim(0, 2)
|
|
axFlarePeaks.set_ylim(0.99, maxY)
|
|
axFlarePeaks.set_ylabel("Flare peak")
|
|
axFlarePeaks.set_xlabel("Phase")
|
|
axFlarePeaks.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axFlarePeaks.xaxis.set_major_locator(MaxNLocator(5))
|
|
axFlarePeaks.legend()
|
|
|
|
plt.savefig(f"{locFolder}/{''.join(combo)}-Flarepeaks_maxY-{maxY}.png")
|
|
plt.close()
|
|
|
|
for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyellow"]):
|
|
spTypes = [f"{sT}0", f"{sT}1", f"{sT}2", f"{sT}3", f"{sT}4", f"{sT}5", f"{sT}6", f"{sT}7", f"{sT}8", f"{sT}9"]
|
|
for spTyp in spTypes:
|
|
finalData = pd.DataFrame()
|
|
|
|
Mfilter = data["SpType"].str.startswith(spTyp)
|
|
numStars = len(set(data[Mfilter]["StarName"]))
|
|
Mfilter &= showSourceFilter
|
|
finalData = pd.concat([finalData, data[Mfilter]], ignore_index=True)
|
|
|
|
pdcsapbinningData = []
|
|
locFolder = f"{folderPath}/{spTyp}/"
|
|
mkdir_p(f"{locFolder}/")
|
|
csvFile = open(f"{locFolder}/{spTyp}.csv", "a")
|
|
csvFile.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Normalized Phase of Peak,Peak in Period")
|
|
csvFile.write("\n")
|
|
for ind, row in finalData.reset_index().iterrows():
|
|
PDCSAPminOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][0]
|
|
PDCSAPmaxOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][1]
|
|
if(len(row["pdcsapFoldedPeaksPhasePair"]) > 0):
|
|
pdcsapVals = pd.DataFrame(row["pdcsapFoldedPeaksPhasePair"])
|
|
for td, peak, pv in zip(pdcsapVals["Phase"], pdcsapVals["Peak"], row["pdcsapPeaks"]):
|
|
normPhase = normalizePhase(td.value, np.abs(PDCSAPminOrigPhase), np.abs(PDCSAPmaxOrigPhase))
|
|
csvFile.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{normPhase},{peak['FlarePeak']}")
|
|
csvFile.write("\n")
|
|
pdcsapbinningData.append({"SpType": f'{row["SpType"][0]}{row["SpType"][1]}',
|
|
"PDCSAPNormPhase": normPhase,
|
|
"Peak": peak["FlarePeak"]})
|
|
if(peak["FlarePeak"] > 100):
|
|
print(row["StarName"], "has over 100 peak")
|
|
csvFile.close()
|
|
pdcsapbinningData = pd.DataFrame(pdcsapbinningData)
|
|
PDCSAPdataList = []
|
|
PDCSAPlabelList = []
|
|
PDCSAPcolorList = []
|
|
PDCSAPdataList2dhistPhase = []
|
|
PDCSAPdataList2dhistPeak = []
|
|
PDCSAPlabelList2dhist = []
|
|
PDCSAPcolorList2dhist = []
|
|
try:
|
|
SpTypefilter = pdcsapbinningData["SpType"] == spTyp
|
|
except:
|
|
shutil.rmtree(locFolder)
|
|
continue
|
|
|
|
PDCSAPdataList2dhistPhase.append(pdcsapbinningData[SpTypefilter]["PDCSAPNormPhase"])
|
|
PDCSAPdataList2dhistPeak.append(pdcsapbinningData[SpTypefilter]["Peak"])
|
|
PDCSAPlabelList2dhist.append(f"{spTyp} Stars")
|
|
PDCSAPcolorList2dhist.append(color)
|
|
|
|
xData = pd.DataFrame()
|
|
yData = pd.DataFrame()
|
|
for aX, aY in zip(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak):
|
|
xData = pd.concat([xData, aX], ignore_index=True)
|
|
yData = pd.concat([yData, aY], ignore_index=True)
|
|
|
|
xData = np.asarray(xData.values)[:,0]
|
|
yData = np.asarray(yData.values)[:,0]
|
|
PDCSAPdataList.append(pdcsapbinningData[SpTypefilter]["PDCSAPNormPhase"])
|
|
|
|
PDCSAPlabelList.append(f"{spTyp} Stars")
|
|
PDCSAPcolorList.append(color)
|
|
|
|
for bins in binList:
|
|
figHisto, ((axHisto)) = plt.subplots(nrows=1, ncols=1)
|
|
y, binEdges, _ = axHisto.hist(PDCSAPdataList, bins,
|
|
label=PDCSAPlabelList,
|
|
color=PDCSAPcolorList,
|
|
stacked=True,
|
|
range=[0, 2])
|
|
bincenters = 0.5*(binEdges[1:]+binEdges[:-1])
|
|
if(isinstance(y[0], np.ndarray)):
|
|
y = y[-1]
|
|
n_i = y
|
|
m_i = bincenters * np.pi
|
|
N = np.sum(n_i)
|
|
mean = np.sum(n_i * m_i)/N
|
|
stdDev = np.sqrt(np.sum(((n_i - mean)**2)) / (N-1))
|
|
menStd = np.sqrt(y)
|
|
axHisto.bar(bincenters[y > 0], y[y > 0], width=0, color='r', yerr=stdDev)
|
|
if(~np.isnan(stdDev) & ~np.isinf(stdDev)):
|
|
axHisto.set_ylim(0, max(y[y > 0 & ~np.isnan(y) & ~np.isinf(y)] + stdDev))
|
|
axHisto.set_ylabel("Num. flares")
|
|
axHisto.set_xlabel("Phase")
|
|
axHisto.set_title(f"Flare count in phase of {spTyp} type stars with {bins} bins ({numStars} stars)")
|
|
axHisto.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axHisto.xaxis.set_major_locator(MaxNLocator(5))
|
|
axHisto.legend()
|
|
|
|
axHistoPhase = axHisto.twinx()
|
|
secAxisXdata = np.linspace(0, 2, num=10000)
|
|
secAxisYdata = np.cos(secAxisXdata*np.pi) + 1
|
|
axHistoPhase.plot(secAxisXdata, secAxisYdata)
|
|
axHistoPhase.set_ylim(0, 7)
|
|
|
|
plt.savefig(f"{locFolder}/{''.join(spTyp)}-Flarecount-{bins}_Bins.png")
|
|
plt.close()
|
|
|
|
for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
|
|
figFlarepeakHist, ((axFlarepeakHist)) = plt.subplots(nrows=1, ncols=1)
|
|
axFlarepeakHist.set_ylabel("Flare peak")
|
|
axFlarepeakHist.set_xlabel("Phase")
|
|
H, xedges, yedges = np.histogram2d(xData, yData, bins=bins, range=[[0, 2], [min(yData), maxY]])
|
|
cmax = 11
|
|
H_clipped = np.clip(H, None, cmax)
|
|
im = axFlarepeakHist.imshow(H_clipped.T, origin='lower', interpolation='nearest',
|
|
extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]],
|
|
aspect='auto', cmap='viridis')
|
|
figFlarepeakHist.colorbar(im, label='Counts', ax=axFlarepeakHist)
|
|
axFlarepeakHist.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axFlarepeakHist.xaxis.set_major_locator(MaxNLocator(5))
|
|
axFlarepeakHist.set_title(f"Flare peak per phase histogram of {spTyp} type stars with {bins} bins ({numStars} stars)")
|
|
plt.savefig(f"{locFolder}/{''.join(spTyp)}-Flarepeaks-{bins}_Bins_maxY-{maxY}.png")
|
|
plt.close()
|
|
|
|
for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]:
|
|
figFlarePeaks, ((axFlarePeaks)) = plt.subplots(nrows=1, ncols=1)
|
|
axFlarePeaks.scatter(pdcsapbinningData[SpTypefilter]["PDCSAPNormPhase"],
|
|
pdcsapbinningData[SpTypefilter]["Peak"],
|
|
label=f"{spTyp} Stars", color=color)
|
|
axFlarePeaks.set_xlim(0, 2)
|
|
axFlarePeaks.set_ylim(0.99, maxY)
|
|
axFlarePeaks.set_ylabel("Flare peak")
|
|
axFlarePeaks.set_xlabel("Phase")
|
|
axFlarePeaks.set_title(f"Flare peaks per phase of {spTyp} type stars ({numStars} stars)")
|
|
axFlarePeaks.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$'))
|
|
axFlarePeaks.xaxis.set_major_locator(MaxNLocator(5))
|
|
axFlarePeaks.legend()
|
|
|
|
plt.savefig(f"{locFolder}/{''.join(spTyp)}-Flarepeaks_maxY-{maxY}.png")
|
|
plt.close()
|