From a901ba9c45ce944ee8c845b78723e67f94b6515e Mon Sep 17 00:00:00 2001 From: SGCMarkus Date: Wed, 4 Dec 2024 14:05:26 +0100 Subject: [PATCH] generate_plots_MKGF: generate plots individually per star too --- generate_plots_MKGF.py | 149 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 143 insertions(+), 6 deletions(-) diff --git a/generate_plots_MKGF.py b/generate_plots_MKGF.py index 808c7f2..47c8a0a 100644 --- a/generate_plots_MKGF.py +++ b/generate_plots_MKGF.py @@ -27,7 +27,7 @@ def mkdir_p(mypath): pass else: raise -fileName = "datav3.cff" +fileName = "datav4.cff" data = pd.read_pickle(fileName) binList = [10, 20, 30] @@ -245,10 +245,10 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel 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) + 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}/" @@ -342,7 +342,7 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel 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]]) + 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', @@ -371,3 +371,140 @@ for sT, color in zip(["M", "K", "G", "F"], ["red", "orange", "yellow", "greenyel 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()