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 = "datav4.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): for maxFlarePeak in [1.01, 1.05, 1.1, 1.25, 1.5]: 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)}_maxFlarePeak_{maxFlarePeak}.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"]): if(peak["FlarePeak"] <= maxFlarePeak): 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() if(len(pdcsapbinningData) < 1): continue 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) if(np.isinf(stdDev)): stdDev = np.mean(menStd) 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)}_maxFlarePeak_{maxFlarePeak}-Flarecount-{bins}_Bins.png") plt.close() 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() typeFilter = data["SpType"].str.startswith(spTyp) numStars = len(set(data[typeFilter]["StarName"])) typeFilter &= showSourceFilter finalData = pd.concat([finalData, data[typeFilter]], 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], [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)) 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() starDB: StarDB = StarDB.getInstance("stars.db") for starName in starDB.getAllStars(): finalData = pd.DataFrame() starNameR = starName.replace('*', '_star_') nameFilter = data["StarName"] == starName nameFilter &= showSourceFilter finalData = pd.concat([finalData, data[nameFilter]], ignore_index=True) color = "" pdcsapbinningData = [] locFolder = f"{folderPath}/stars/{starNameR}/" mkdir_p(f"{locFolder}/") csvFile = open(f"{locFolder}/{starNameR}.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(row["SpType"][0] == "M"): color = "red" elif(row["SpType"][0] == "K"): color = "orange" elif(row["SpType"][0] == "G"): color = "yellow" elif(row["SpType"][0] == "F"): color = "greenyellow" else: color = "gray" 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:2] if len(row["SpType"]) > 1 else 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 = [] if(len(pdcsapbinningData) < 1): continue PDCSAPdataList2dhistPhase.append(pdcsapbinningData[:]["PDCSAPNormPhase"]) PDCSAPdataList2dhistPeak.append(pdcsapbinningData[:]["Peak"]) 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[:]["PDCSAPNormPhase"]) PDCSAPlabelList.append(f"{starName}") 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 {starName} with {bins} bins") 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}/{starNameR}-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], [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)) axFlarepeakHist.set_title(f"Flare peak per phase histogram of {starName} with {bins} bins") plt.savefig(f"{locFolder}/{starNameR}-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[:]["PDCSAPNormPhase"], pdcsapbinningData[:]["Peak"], label=f"{starName}", 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 {starName}") axFlarePeaks.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$')) axFlarePeaks.xaxis.set_major_locator(MaxNLocator(5)) axFlarePeaks.legend() plt.savefig(f"{locFolder}/{starNameR}-Flarepeaks_maxY-{maxY}.png") plt.close()