From 01ae4f9084ef82bfffcfe393fc58fce068ef64a3 Mon Sep 17 00:00:00 2001 From: SGCMarkus Date: Mon, 10 Mar 2025 08:26:16 +0100 Subject: [PATCH] rework generate_plot script --- generate_plots.py | 692 +++++++++++++++++++++++++++++++++++++++++ generate_plots_MKGF.py | 667 --------------------------------------- 2 files changed, 692 insertions(+), 667 deletions(-) create mode 100644 generate_plots.py delete mode 100644 generate_plots_MKGF.py diff --git a/generate_plots.py b/generate_plots.py new file mode 100644 index 0000000..4666c1d --- /dev/null +++ b/generate_plots.py @@ -0,0 +1,692 @@ +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) + +# remove any data that has no period +data = data[(data["FitType"] == "sine") | (data["FitType"] == "poly")] + +# remove data with multiple minima/maxima present +filterArray = [] +for ind, row in data.reset_index().iterrows(): + if(len(row["pdcsapPeriodMinima"]) == len(row["pdcsapPeriodMaxima"]) and + len(row["pdcsapPeriodMinima"]) == 1): + filterArray.append(True) + else: + filterArray.append(False) + +filterArray = np.array(filterArray) +data = data[filterArray] + +validStarPeriodMap = [] +starList = set(list(data["StarName"])) + +def allValuesWithin3Std(values: list): + if(not values or len(values) == 1): + return True + + return max(values) - min(values) <= 3*np.std(values) + +def getMeanPeriod(values: list): + return np.mean(values) + +for starName in starList: + periods = list(data[data["StarName"] == starName]["pdcsapPeriod"]) + + validStarPeriodMap.append({"StarName": starName, + "MeanPeriod": getMeanPeriod(periods), + "PeriodWithinStd": allValuesWithin3Std(periods)}) + +validStarPeriodMap = pd.DataFrame(validStarPeriodMap) +periodsCutList = [0.5, 1, 1.5, 2, 5, 10, 15, 20] + +data = pd.merge(data, validStarPeriodMap, on="StarName") +starDB: StarDB = StarDB.getInstance("stars.db") + +def plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, title, filename): + 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) or np.isnan(stdDev)): + stdDev = np.mean(menStd) + if(np.isinf(stdDev) or np.isnan(stdDev)): + stdDev = 0 + axHisto.bar(bincenters[y > 0], y[y > 0], width=0, color='r', yerr=stdDev) + try: + axHisto.set_ylim(0, max(y) + stdDev) + except Exception as e: + print(e) + print(max(y), stdDev) + quit() + axHisto.set_ylabel("Num. flares") + axHisto.set_xlabel("Phase") + axHisto.set_title(title) + 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(filename) + plt.close() + +def plotFlarePhasePeakHistogram(xData, yData, bins, maxY, title, filename): + 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(title) + plt.savefig(filename) + plt.close() + +def plotFlarePeaks(plotdata, filters, labels, colors, maxY, title, filename): + if(isinstance(labels, list) and isinstance(colors, list)): + _, ((axFlarePeaks)) = plt.subplots(nrows=1, ncols=1) + for f, l, c in zip(filters, labels, colors): + axFlarePeaks.scatter(plotdata[f]["PDCSAPNormPhase"], + plotdata[f]["Peak"], + label=l, color=c) + elif(isinstance(labels, str) and isinstance(colors, str)): + _, ((axFlarePeaks)) = plt.subplots(nrows=1, ncols=1) + if(filters is None): + axFlarePeaks.scatter(plotdata[:]["PDCSAPNormPhase"], + plotdata[:]["Peak"], + label=labels, color=colors) + else: + axFlarePeaks.scatter(plotdata[filters]["PDCSAPNormPhase"], + plotdata[filters]["Peak"], + label=labels, color=colors) + else: + return + axFlarePeaks.set_xlim(0, 2) + axFlarePeaks.set_ylim(0.99, maxY) + axFlarePeaks.set_ylabel("Flare peak") + axFlarePeaks.set_xlabel("Phase") + axFlarePeaks.set_title(title) + axFlarePeaks.xaxis.set_major_formatter(tck.FormatStrFormatter('%g $\pi$')) + axFlarePeaks.xaxis.set_major_locator(MaxNLocator(5)) + axFlarePeaks.legend() + plt.savefig(filename) + plt.close() + +def setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, + PDCSAPdataList=None, dataFilter=None): + 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] + if(PDCSAPdataList is not None): + if(dataFilter is not None): + PDCSAPdataList.append(pdcsapbinningData[dataFilter]["PDCSAPNormPhase"]) + else: + PDCSAPdataList.append(pdcsapbinningData[:]["PDCSAPNormPhase"]) + + return xData, yData, PDCSAPdataList + +def generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, histogramTitleArg, histogramFilenameArg, + xData, yData, peak2DHistogramTitleArg, peak2DfilenameArg, + plotdata, filters, flarePlotLabels, flarePlotColors, flarePlotTitleArg, flarePlotFilenameArg): + for bins in binList: + histogramTitle = histogramTitleArg.replace("@bins", str(bins)) + histogramFilename = histogramFilenameArg.replace("@bins", str(bins)) + plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, histogramTitle, histogramFilename) + + for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]: + peak2DHistogramTitle = peak2DHistogramTitleArg.replace("@bins", str(bins)) + peak2Dfilename = peak2DfilenameArg.replace("@bins", str(bins)).replace("@maxY", str(maxY)) + plotFlarePhasePeakHistogram(xData, yData, bins, maxY, peak2DHistogramTitle, peak2Dfilename) + + for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]: + flarePlotTitle = flarePlotTitleArg + flarePlotFilename = flarePlotFilenameArg.replace("@maxY", str(maxY)) + plotFlarePeaks(plotdata, filters, flarePlotLabels, flarePlotColors, maxY, flarePlotTitle, flarePlotFilename) + +# All stars +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) + + 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(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 + else: + if(pdcsapbinningData["SpType"][0][0] == "M"): + color = "red" + elif(pdcsapbinningData["SpType"][0][0] == "K"): + color = "orange" + elif(pdcsapbinningData["SpType"][0][0] == "G"): + color = "yellow" + elif(pdcsapbinningData["SpType"][0][0] == "F"): + color = "greenyellow" + else: + color = "gray" + + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[:]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[:]["Peak"]) + + xData, yData, PDCSAPdataList = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, PDCSAPdataList) + + PDCSAPlabelList.append(f"{starName}") + PDCSAPcolorList.append(color) + + plotdata = pdcsapbinningData + filters = None + flarePlotLabels = f"{starName}" + flarePlotColors = color + generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, f"Flare count in phase of {starName} with @bins bins", f"{locFolder}/{starNameR}-Flarecount-@bins_Bins.png", + xData, yData, f"Flare peak per phase histogram of {starName} with @bins bins", f"{locFolder}/{starNameR}-Flarepeaks-@bins_Bins_maxY-@maxY.png", + plotdata, filters, flarePlotLabels, flarePlotColors, f"Flare peaks per phase of {starName}", f"{locFolder}/{starNameR}-Flarepeaks_maxY-@maxY.png") + + +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]: + # max Flare Peak cut + finalDataMaxFlarePeak = pd.DataFrame() + if("M" in combo): + Mfilter = data["SpType"].str.startswith("M") + Mfilter &= showSourceFilter + finalDataMaxFlarePeak = pd.concat([finalDataMaxFlarePeak, data[Mfilter]], ignore_index=True) + if("K" in combo): + Kfilter = data["SpType"].str.startswith("K") + Kfilter &= showSourceFilter + finalDataMaxFlarePeak = pd.concat([finalDataMaxFlarePeak, data[Kfilter]], ignore_index=True) + if("G" in combo): + Gfilter = data["SpType"].str.startswith("G") + Gfilter &= showSourceFilter + finalDataMaxFlarePeak = pd.concat([finalDataMaxFlarePeak, data[Gfilter]], ignore_index=True) + if("F" in combo): + Ffilter = data["SpType"].str.startswith("F") + Ffilter &= showSourceFilter + finalDataMaxFlarePeak = pd.concat([finalDataMaxFlarePeak, data[Ffilter]], ignore_index=True) + + numStars = len(set(finalDataMaxFlarePeak["StarName"])) + + pdcsapbinningData = [] + locFolder = f"{folderPath}/{''.join(combo)}/maxFlarePeaks/{maxFlarePeak}/" + 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 finalDataMaxFlarePeak.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 = [] + plotFilters = [] + if("M" in combo): + Mfilter = pdcsapbinningData["SpType"] == "M" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Mfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("M Stars") + PDCSAPcolorList.append("red") + plotFilters.append(Mfilter) + if("K" in combo): + Kfilter = pdcsapbinningData["SpType"] == "K" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Kfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("K Stars") + PDCSAPcolorList.append("orange") + plotFilters.append(Kfilter) + if("G" in combo): + Gfilter = pdcsapbinningData["SpType"] == "G" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Gfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("G Stars") + PDCSAPcolorList.append("yellow") + plotFilters.append(Gfilter) + if("F" in combo): + Ffilter = pdcsapbinningData["SpType"] == "F" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Ffilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("F Stars") + PDCSAPcolorList.append("greenyellow") + plotFilters.append(Ffilter) + + xData, yData, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak) + plotdata = pdcsapbinningData + filters = plotFilters + generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, f"Flare count per phase of {', '.join(combo)} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(combo)}_maxFlarePeak_{maxFlarePeak}-Flarecount-@bins_Bins.png", + xData, yData, f"Flare peak per phase histogram of {', '.join(combo)} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(combo)}_maxFlarePeak_{maxFlarePeak}-Flarepeaks-@bins_Bins_maxY-@maxY.png", + plotdata, filters, PDCSAPlabelList, PDCSAPcolorList, f"Flare peaks per phase of {', '.join(combo)} type stars ({numStars} stars)", f"{locFolder}/{''.join(combo)}_maxFlarePeak_{maxFlarePeak}-Flarepeaks_maxY-@maxY.png") + + # all flare peaks + finalDataAllFlarePeaks = pd.DataFrame() + if("M" in combo): + Mfilter = data["SpType"].str.startswith("M") + Mfilter &= showSourceFilter + finalDataAllFlarePeaks = pd.concat([finalDataAllFlarePeaks, data[Mfilter]], ignore_index=True) + if("K" in combo): + Kfilter = data["SpType"].str.startswith("K") + Kfilter &= showSourceFilter + finalDataAllFlarePeaks = pd.concat([finalDataAllFlarePeaks, data[Kfilter]], ignore_index=True) + if("G" in combo): + Gfilter = data["SpType"].str.startswith("G") + Gfilter &= showSourceFilter + finalDataAllFlarePeaks = pd.concat([finalDataAllFlarePeaks, data[Gfilter]], ignore_index=True) + if("F" in combo): + Ffilter = data["SpType"].str.startswith("F") + Ffilter &= showSourceFilter + finalDataAllFlarePeaks = pd.concat([finalDataAllFlarePeaks, data[Ffilter]], ignore_index=True) + + numStars = len(set(finalDataAllFlarePeaks["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 finalDataAllFlarePeaks.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 = [] + plotFilters = [] + if("M" in combo): + Mfilter = pdcsapbinningData["SpType"] == "M" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Mfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Mfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("M Stars") + PDCSAPcolorList.append("red") + plotFilters.append(Mfilter) + if("K" in combo): + Kfilter = pdcsapbinningData["SpType"] == "K" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Kfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Kfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("K Stars") + PDCSAPcolorList.append("orange") + plotFilters.append(Kfilter) + if("G" in combo): + Gfilter = pdcsapbinningData["SpType"] == "G" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Gfilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Gfilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("G Stars") + PDCSAPcolorList.append("yellow") + plotFilters.append(Gfilter) + if("F" in combo): + Ffilter = pdcsapbinningData["SpType"] == "F" + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[Ffilter]["Peak"]) + PDCSAPdataList.append(pdcsapbinningData[Ffilter]["PDCSAPNormPhase"]) + PDCSAPlabelList.append("F Stars") + PDCSAPcolorList.append("greenyellow") + plotFilters.append(Ffilter) + + xData, yData, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak) + plotdata = pdcsapbinningData + filters = plotFilters + generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, f"Flare count per phase of {', '.join(combo)} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(combo)}-Flarecount-@bins_Bins.png", + xData, yData, f"Flare peak per phase histogram of {', '.join(combo)} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(combo)}-Flarepeaks-@bins_Bins_maxY-@maxY.png", + plotdata, filters, PDCSAPlabelList, PDCSAPcolorList, f"Flare peaks per phase of {', '.join(combo)} type stars ({numStars} stars)", f"{locFolder}/{''.join(combo)}-Flarepeaks_maxY-@maxY.png") + + if(comboLength == 1): + spTypes = [f"{combo}0", f"{combo}1", f"{combo}2", f"{combo}3", f"{combo}4", f"{combo}5", f"{combo}6", f"{combo}7", f"{combo}8", f"{combo}9"] + match combo: + case "M": + color = "red" + case "K": + color = "orange" + case "G": + color = "yellow" + case "F": + color = "greenyellow" + + for spTyp in spTypes: + finalDataAccSpTypes = pd.DataFrame() + + typeFilter = data["SpType"].str.startswith(spTyp) + numStars = len(set(data[typeFilter]["StarName"])) + typeFilter &= showSourceFilter + finalDataAccSpTypes = pd.concat([finalDataAccSpTypes, 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 finalDataAccSpTypes.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"]) + PDCSAPlabelList.append(f"{spTyp} Stars") + PDCSAPcolorList.append(color) + xData, yData, PDCSAPdataList = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, PDCSAPdataList, SpTypefilter) + + plotdata = pdcsapbinningData + filters = SpTypefilter + generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, f"Flare count in phase of {spTyp} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(spTyp)}-Flarecount-@bins_Bins.png", + xData, yData, f"Flare peak per phase histogram of {spTyp} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(spTyp)}-Flarepeaks-@bins_Bins_maxY-@maxY.png", + plotdata, filters, PDCSAPlabelList, PDCSAPcolorList, f"Flare peaks per phase of {spTyp} type stars ({numStars} stars)", f"{locFolder}/{''.join(spTyp)}-Flarepeaks_maxY-@maxY.png") + + # per Period + for periodCut in periodsCutList: + finalDataPeriod = pd.DataFrame() + if("M" in combo): + Mfilter = data["SpType"].str.startswith("M") + Mfilter &= showSourceFilter + Mfilter &= data["PeriodWithinStd"] == True + finalDataPeriod = pd.concat([finalDataPeriod, data[Mfilter]], ignore_index=True) + if("K" in combo): + Kfilter = data["SpType"].str.startswith("K") + Kfilter &= showSourceFilter + Kfilter &= data["PeriodWithinStd"] == True + finalDataPeriod = pd.concat([finalDataPeriod, data[Kfilter]], ignore_index=True) + if("G" in combo): + Gfilter = data["SpType"].str.startswith("G") + Gfilter &= showSourceFilter + Gfilter &= data["PeriodWithinStd"] == True + finalDataPeriod = pd.concat([finalDataPeriod, data[Gfilter]], ignore_index=True) + if("F" in combo): + Ffilter = data["SpType"].str.startswith("F") + Ffilter &= showSourceFilter + Ffilter &= data["PeriodWithinStd"] == True + finalDataPeriod = pd.concat([finalDataPeriod, data[Ffilter]], ignore_index=True) + pdcsapbinningDataU = [] + pdcsapbinningDataO = [] + locFolder = f"{folderPath}/{''.join(combo)}/periodCuts/{periodCut}/" + mkdir_p(f"{locFolder}/") + csvFileU = open(f"{locFolder}/under_{periodCut}.csv", "a") + csvFileO = open(f"{locFolder}/over_{periodCut}.csv", "a") + csvFileU.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Mean Period,Normalized Phase of Peak,Peak in Period") + csvFileU.write("\n") + csvFileO.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Mean Period,Normalized Phase of Peak,Peak in Period") + csvFileO.write("\n") + for ind, row in finalDataPeriod.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)) + if(row["MeanPeriod"] <= periodCut): + csvFileU.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{row["MeanPeriod"]},{normPhase},{peak['FlarePeak']}") + csvFileU.write("\n") + pdcsapbinningDataU.append({"SpType": f'{row["SpType"][0]}', + "PDCSAPNormPhase": normPhase, + "Peak": peak["FlarePeak"]}) + else: + csvFileO.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{row["MeanPeriod"]},{normPhase},{peak['FlarePeak']}") + csvFileO.write("\n") + pdcsapbinningDataO.append({"SpType": f'{row["SpType"][0]}', + "PDCSAPNormPhase": normPhase, + "Peak": peak["FlarePeak"]}) + if(peak["FlarePeak"] > 100): + print(row["StarName"], "has over 100 peak") + csvFileU.close() + csvFileO.close() + pdcsapbinningDataU = pd.DataFrame(pdcsapbinningDataU) + pdcsapbinningDataO = pd.DataFrame(pdcsapbinningDataO) + PDCSAPdataListU = []; PDCSAPdataListO = [] + PDCSAPlabelList = [] + PDCSAPcolorList = [] + PDCSAPdataList2dhistPhaseU = []; PDCSAPdataList2dhistPeakU = [] + PDCSAPdataList2dhistPhaseO = []; PDCSAPdataList2dhistPeakO = [] + plotFiltersU = []; plotFiltersO = []; + + if("M" in combo): + PDCSAPlabelList.append("M Stars") + PDCSAPcolorList.append("red") + if("K" in combo): + PDCSAPlabelList.append("K Stars") + PDCSAPcolorList.append("orange") + if("G" in combo): + PDCSAPlabelList.append("G Stars") + PDCSAPcolorList.append("yellow") + if("F" in combo): + PDCSAPlabelList.append("F Stars") + PDCSAPcolorList.append("greenyellow") + + if(len(pdcsapbinningDataU) > 0): + if("M" in combo): + MfilterU = pdcsapbinningDataU["SpType"] == "M" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[MfilterU]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[MfilterU]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[MfilterU]["PDCSAPNormPhase"]) + plotFiltersU.append(MfilterU) + if("K" in combo): + KfilterU = pdcsapbinningDataU["SpType"] == "K" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[KfilterU]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[KfilterU]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[KfilterU]["PDCSAPNormPhase"]) + plotFiltersU.append(KfilterU) + if("G" in combo): + GfilterU = pdcsapbinningDataU["SpType"] == "G" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[GfilterU]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[GfilterU]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[GfilterU]["PDCSAPNormPhase"]) + plotFiltersU.append(GfilterU) + if("F" in combo): + FfilterU = pdcsapbinningDataU["SpType"] == "F" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[FfilterU]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[FfilterU]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[FfilterU]["PDCSAPNormPhase"]) + plotFiltersU.append(FfilterU) + + if(len(pdcsapbinningDataO) > 0): + if("M" in combo): + MfilterO = pdcsapbinningDataO["SpType"] == "M" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[MfilterO]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[MfilterO]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[MfilterO]["PDCSAPNormPhase"]) + plotFiltersO.append(MfilterO) + if("K" in combo): + KfilterO = pdcsapbinningDataO["SpType"] == "K" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[KfilterO]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[KfilterO]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[KfilterO]["PDCSAPNormPhase"]) + plotFiltersO.append(KfilterO) + if("G" in combo): + GfilterO = pdcsapbinningDataO["SpType"] == "G" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[GfilterO]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[GfilterO]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[GfilterO]["PDCSAPNormPhase"]) + plotFiltersO.append(GfilterO) + if("F" in combo): + FfilterO = pdcsapbinningDataO["SpType"] == "F" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[FfilterO]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[FfilterO]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[FfilterO]["PDCSAPNormPhase"]) + plotFiltersO.append(FfilterO) + + if(len(PDCSAPdataList2dhistPhaseU) > 0 and len(PDCSAPdataList2dhistPeakU) > 0): + xDataU, yDataU, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhaseU, PDCSAPdataList2dhistPeakU) + if(len(PDCSAPdataListU) > 0 and len(xDataU) > 0 and len(yDataU) > 0): + generatePlots(PDCSAPdataListU, PDCSAPlabelList, PDCSAPcolorList, f"Flare count per phase of {', '.join(combo)} type stars with @bins bins (Rot. Period under {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarecount-Period_u_{periodCut}-@bins_Bins.png", + xDataU, yDataU, f"Flare peak per phase histogram of {', '.join(combo)} type stars with @bins bins (Rot. Period under {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarepeaks-Period_u_{periodCut}-@bins_Bins_maxY-@maxY.png", + pdcsapbinningDataU, plotFiltersU, PDCSAPlabelList, PDCSAPcolorList, f"Flare peaks per phase of {', '.join(combo)} type stars (Rot. Period under {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarepeaks-Period_u_{periodCut}-maxY_@maxY.png") + + if(len(PDCSAPdataList2dhistPhaseO) > 0 and len(PDCSAPdataList2dhistPeakO) > 0): + xDataO, yDataO, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhaseO, PDCSAPdataList2dhistPeakO) + if(len(PDCSAPdataListO) > 0 and len(xDataO) > 0 and len(yDataO) > 0): + generatePlots(PDCSAPdataListO, PDCSAPlabelList, PDCSAPcolorList, f"Flare count per phase of {', '.join(combo)} type stars with @bins bins (Rot. Period over {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarecount-Period_o_{periodCut}-@bins_Bins.png", + xDataO, yDataO, f"Flare peak per phase histogram of {', '.join(combo)} type stars with @bins bins (Rot. Period over {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarepeaks-Period_o_{periodCut}-@bins_Bins_maxY-@maxY.png", + pdcsapbinningDataO, plotFiltersO, PDCSAPlabelList, PDCSAPcolorList, f"Flare peaks per phase of {', '.join(combo)} type stars (Rot. Period over {periodCut} days)", f"{locFolder}/{''.join(combo)}-Flarepeaks-Period_o_{periodCut}-maxY_@maxY.png") + diff --git a/generate_plots_MKGF.py b/generate_plots_MKGF.py deleted file mode 100644 index 7469edb..0000000 --- a/generate_plots_MKGF.py +++ /dev/null @@ -1,667 +0,0 @@ -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) - -# remove any data that has no period -data = data[(data["FitType"] == "sine") | (data["FitType"] == "poly")] - -# remove data with multiple minima/maxima present -filterArray = [] -for ind, row in data.reset_index().iterrows(): - if(len(row["pdcsapPeriodMinima"]) == len(row["pdcsapPeriodMaxima"]) and - len(row["pdcsapPeriodMinima"]) == 1): - filterArray.append(True) - else: - filterArray.append(False) - -filterArray = np.array(filterArray) -data = data[filterArray] - -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()