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 = "datav3.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()