diff --git a/datav5.cff b/datav5.cff new file mode 100644 index 0000000..dd0c339 Binary files /dev/null and b/datav5.cff differ diff --git a/generate_plots.py b/generate_plots.py index 41b8ea8..9e63b9b 100644 --- a/generate_plots.py +++ b/generate_plots.py @@ -27,7 +27,7 @@ def mkdir_p(mypath): pass else: raise -fileName = "datav4.cff" +fileName = "datav5.cff" data = pd.read_pickle(fileName) binList = [10, 20, 30] @@ -56,19 +56,8 @@ folderPath = f"../{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] +#data = data[(data["FitType"] == "sine") | (data["FitType"] == "poly")] +data = data[(data["FitType"] == "sine")] validStarPeriodMap = [] starList = set(list(data["StarName"])) @@ -95,7 +84,7 @@ 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): +def plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, title, filename, foldedFits): figHisto, ((axHisto)) = plt.subplots(nrows=1, ncols=1) y, binEdges, _ = axHisto.hist(PDCSAPdataList, bins, label=PDCSAPlabelList, @@ -130,10 +119,17 @@ def plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, ti 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) + if(foldedFits is None): + secAxisXdata = np.linspace(0, 2, num=10000) + secAxisYdata = np.cos(secAxisXdata*np.pi) + 1 + axHistoPhase.plot(secAxisXdata, secAxisYdata, color="blue") + axHistoPhase.set_ylim(0, 7) + else: + for fit in foldedFits: + axHistoPhase.plot(fit[0], fit[1], color="blue") + fitCol = [fit[1] for fit in foldedFits] + fitCol = np.array(list(itertools.chain.from_iterable(fitCol))) + axHistoPhase.set_ylim(np.min(fitCol), np.max(fitCol)*1.2) plt.savefig(filename) plt.close() @@ -205,11 +201,12 @@ def setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, def generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, histogramTitleArg, histogramFilenameArg, xData, yData, peak2DHistogramTitleArg, peak2DfilenameArg, - plotdata, filters, flarePlotLabels, flarePlotColors, flarePlotTitleArg, flarePlotFilenameArg): + plotdata, filters, flarePlotLabels, flarePlotColors, flarePlotTitleArg, flarePlotFilenameArg, + foldedFits=None): for bins in binList: histogramTitle = histogramTitleArg.replace("@bins", str(bins)) histogramFilename = histogramFilenameArg.replace("@bins", str(bins)) - plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, histogramTitle, histogramFilename) + plotBinsHistogram(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, bins, histogramTitle, histogramFilename, foldedFits) for maxY in [1.05, 1.1, 1.2, 1.5, 2, 2.5, 3, 5, max(yData)]: peak2DHistogramTitle = peak2DHistogramTitleArg.replace("@bins", str(bins)) @@ -230,10 +227,13 @@ for starName in starDB.getAllStars(): finalData = pd.concat([finalData, data[nameFilter]], ignore_index=True) pdcsapbinningData = [] + pdcsapbinningDataSpotModDiffPeriod = [] + foldedFits = [] + foldedPeriodFits = [] 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("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Spot Modulation,Normalized Phase of Peak,Peak in Period") csvFile.write("\n") for ind, row in finalData.reset_index().iterrows(): PDCSAPminOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][0] @@ -243,24 +243,34 @@ for starName in starDB.getAllStars(): 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(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{row['pdcsapSpotModulation']},{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") + foldedFits.append([normalizePhase(row["pdcsapFoldedFitPhase"]), row["pdcsapFoldedFit"]]) + if(row["pdcsapPeriodFoldedPhase"] is not None): + periodPdcsapVals = pd.DataFrame(row["pdcsapPeriodFoldedPeaksPhasePair"]) + for td, peak in zip(periodPdcsapVals["Phase"], periodPdcsapVals["Peak"]): + normPhasePeriod = normalizePhase(td.value, np.abs(row["pdcsapPeriodFoldedFitPhaseStarEnd"][0]), np.abs(row["pdcsapPeriodFoldedFitPhaseStarEnd"][1])) + pdcsapbinningDataSpotModDiffPeriod.append({"SpType": f'{row["SpType"][0:2] if len(row["SpType"]) > 1 else row["SpType"][0]}', + "PDCSAPNormPhasePeriod": normPhasePeriod, + "PeakPeriod": peak["FlarePeak"]}) + if(peak["FlarePeak"] > 100): + print(row["StarName"], "has over 100 peak") + foldedPeriodFits.append([normalizePhase(row["pdcsapPeriodFoldedFitPhase"]), row["pdcsapPeriodFoldedFit"]]) csvFile.close() pdcsapbinningData = pd.DataFrame(pdcsapbinningData) + pdcsapbinningDataSpotModDiffPeriod = pd.DataFrame(pdcsapbinningDataSpotModDiffPeriod) if len(pdcsapbinningDataSpotModDiffPeriod) > 0 else None PDCSAPdataList = [] PDCSAPlabelList = [] PDCSAPcolorList = [] PDCSAPdataList2dhistPhase = [] PDCSAPdataList2dhistPeak = [] - if(len(pdcsapbinningData) < 1): - continue - else: + if(len(pdcsapbinningData) > 0): if(pdcsapbinningData["SpType"][0][0] == "M"): color = "red" elif(pdcsapbinningData["SpType"][0][0] == "K"): @@ -272,21 +282,49 @@ for starName in starDB.getAllStars(): else: color = "gray" - PDCSAPdataList2dhistPhase.append(pdcsapbinningData[:]["PDCSAPNormPhase"]) - PDCSAPdataList2dhistPeak.append(pdcsapbinningData[:]["Peak"]) + PDCSAPdataList2dhistPhase.append(pdcsapbinningData[:]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningData[:]["Peak"]) - xData, yData, PDCSAPdataList = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, PDCSAPdataList) + xData, yData, PDCSAPdataList = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, PDCSAPdataList) - PDCSAPlabelList.append(f"{starName}") - PDCSAPcolorList.append(color) + 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") + 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", + pdcsapbinningData, None, f"{starName}", color, f"Flare peaks per phase of {starName}", f"{locFolder}/{starNameR}-Flarepeaks_maxY-@maxY.png", + foldedFits) + + if(pdcsapbinningDataSpotModDiffPeriod is not None): + PDCSAPdataList = [] + PDCSAPlabelList = [] + PDCSAPcolorList = [] + PDCSAPdataList2dhistPhase = [] + PDCSAPdataList2dhistPeak = [] + + 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(pdcsapbinningDataSpotModDiffPeriod[:]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeak.append(pdcsapbinningDataSpotModDiffPeriod[:]["Peak"]) + + xData, yData, PDCSAPdataList = setupxyDataAndDataList(PDCSAPdataList2dhistPhase, PDCSAPdataList2dhistPeak, PDCSAPdataList) + + PDCSAPlabelList.append(f"{starName}") + PDCSAPcolorList.append(color) + + generatePlots(PDCSAPdataList, PDCSAPlabelList, PDCSAPcolorList, f"Flare count in phase of {starName} with @bins bins", f"{locFolder}/{starNameR}-Flarecount-@bins_Bins_Period.png", + xData, yData, f"Flare peak per phase histogram of {starName} with @bins bins", f"{locFolder}/{starNameR}-Flarepeaks-@bins_Bins_maxY-@maxY_Period.png", + pdcsapbinningDataSpotModDiffPeriod, None, f"{starName}", color, f"Flare peaks per phase of {starName}", f"{locFolder}/{starNameR}-Flarepeaks_maxY-@maxY_Period.png", + foldedPeriodFits) for comboLength in range(1, len(spType) + 1): @@ -313,12 +351,18 @@ for comboLength in range(1, len(spType) + 1): 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") + pdcsapbinningDataU = [] + pdcsapbinningDataO = [] + locFolderU = f"{folderPath}/{''.join(combo)}/maxFlarePeaks/{maxFlarePeak}/" + locFolderO = f"{folderPath}/{''.join(combo)}/minFlarePeaks/{maxFlarePeak}/" + mkdir_p(f"{locFolderU}/") + mkdir_p(f"{locFolderO}/") + csvFileU = open(f"{locFolderU}/{''.join(combo)}_maxFlarePeak_{maxFlarePeak}.csv", "a") + csvFileU.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Normalized Phase of Peak,Peak in Period") + csvFileU.write("\n") + csvFileO = open(f"{locFolderO}/{''.join(combo)}_minFlarePeak_{maxFlarePeak}.csv", "a") + csvFileO.write("Star Name,Spectral Type,Source,File,Flare Time,Flare Peak,Period,Normalized Phase of Peak,Peak in Period") + csvFileO.write("\n") for ind, row in finalDataMaxFlarePeak.reset_index().iterrows(): PDCSAPminOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][0] PDCSAPmaxOrigPhase = row["pdcsapFoldedFitPhaseStarEnd"][1] @@ -327,66 +371,117 @@ for comboLength in range(1, len(spType) + 1): 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], + csvFileU.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{normPhase},{peak['FlarePeak']}") + csvFileU.write("\n") + pdcsapbinningDataU.append({"SpType": row["SpType"][0], + "PDCSAPNormPhase": normPhase, + "Peak": peak["FlarePeak"]}) + if(peak["FlarePeak"] > 100): + print(row["StarName"], "has over 100 peak") + else: + normPhase = normalizePhase(td.value, np.abs(PDCSAPminOrigPhase), np.abs(PDCSAPmaxOrigPhase)) + csvFileO.write(f"{row['StarName']},{row['SpType']},{row['Source']},{row['FilePath']},{pv['FlarePeakTime']},{pv['FlarePeak']},{row['pdcsapPeriod']},{normPhase},{peak['FlarePeak']}") + csvFileO.write("\n") + pdcsapbinningDataO.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) + csvFileU.close() + csvFileO.close() + if(len(pdcsapbinningDataU) > 0): + pdcsapbinningDataU = pd.DataFrame(pdcsapbinningDataU) + PDCSAPdataListU = [] + PDCSAPlabelListU = [] + PDCSAPcolorListU = [] + PDCSAPdataList2dhistPhaseU = [] + PDCSAPdataList2dhistPeakU = [] + plotFiltersU = [] + if("M" in combo): + Mfilter = pdcsapbinningDataU["SpType"] == "M" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[Mfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[Mfilter]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[Mfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListU.append("M Stars") + PDCSAPcolorListU.append("red") + plotFiltersU.append(Mfilter) + if("K" in combo): + Kfilter = pdcsapbinningDataU["SpType"] == "K" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[Kfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[Kfilter]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[Kfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListU.append("K Stars") + PDCSAPcolorListU.append("orange") + plotFiltersU.append(Kfilter) + if("G" in combo): + Gfilter = pdcsapbinningDataU["SpType"] == "G" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[Gfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[Gfilter]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[Gfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListU.append("G Stars") + PDCSAPcolorListU.append("yellow") + plotFiltersU.append(Gfilter) + if("F" in combo): + Ffilter = pdcsapbinningDataU["SpType"] == "F" + PDCSAPdataList2dhistPhaseU.append(pdcsapbinningDataU[Ffilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakU.append(pdcsapbinningDataU[Ffilter]["Peak"]) + PDCSAPdataListU.append(pdcsapbinningDataU[Ffilter]["PDCSAPNormPhase"]) + PDCSAPlabelListU.append("F Stars") + PDCSAPcolorListU.append("greenyellow") + plotFiltersU.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") + xData, yData, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhaseU, PDCSAPdataList2dhistPeakU) + generatePlots(PDCSAPdataListU, PDCSAPlabelListU, PDCSAPcolorListU, 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", + pdcsapbinningDataU, plotFiltersU, PDCSAPlabelListU, PDCSAPcolorListU, f"Flare peaks per phase of {', '.join(combo)} type stars ({numStars} stars)", f"{locFolder}/{''.join(combo)}_maxFlarePeak_{maxFlarePeak}-Flarepeaks_maxY-@maxY.png") + + if(len(pdcsapbinningDataO) > 0): + pdcsapbinningDataO = pd.DataFrame(pdcsapbinningDataO) + PDCSAPdataListO = [] + PDCSAPlabelListO = [] + PDCSAPcolorListO = [] + PDCSAPdataList2dhistPhaseO = [] + PDCSAPdataList2dhistPeakO = [] + plotFiltersO = [] + if("M" in combo): + Mfilter = pdcsapbinningDataO["SpType"] == "M" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[Mfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[Mfilter]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[Mfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListO.append("M Stars") + PDCSAPcolorListO.append("red") + plotFiltersO.append(Mfilter) + if("K" in combo): + Kfilter = pdcsapbinningDataO["SpType"] == "K" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[Kfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[Kfilter]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[Kfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListO.append("K Stars") + PDCSAPcolorListO.append("orange") + plotFiltersO.append(Kfilter) + if("G" in combo): + Gfilter = pdcsapbinningDataO["SpType"] == "G" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[Gfilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[Gfilter]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[Gfilter]["PDCSAPNormPhase"]) + PDCSAPlabelListO.append("G Stars") + PDCSAPcolorListO.append("yellow") + plotFiltersO.append(Gfilter) + if("F" in combo): + Ffilter = pdcsapbinningDataO["SpType"] == "F" + PDCSAPdataList2dhistPhaseO.append(pdcsapbinningDataO[Ffilter]["PDCSAPNormPhase"]) + PDCSAPdataList2dhistPeakO.append(pdcsapbinningDataO[Ffilter]["Peak"]) + PDCSAPdataListO.append(pdcsapbinningDataO[Ffilter]["PDCSAPNormPhase"]) + PDCSAPlabelListO.append("F Stars") + PDCSAPcolorListO.append("greenyellow") + plotFiltersO.append(Ffilter) + + xData, yData, _ = setupxyDataAndDataList(PDCSAPdataList2dhistPhaseO, PDCSAPdataList2dhistPeakO) + generatePlots(PDCSAPdataListO, PDCSAPlabelListO, PDCSAPcolorListO, f"Flare count per phase of {', '.join(combo)} type stars with @bins bins ({numStars} stars)", f"{locFolder}/{''.join(combo)}_minFlarePeak_{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)}_minFlarePeak_{maxFlarePeak}-Flarepeaks-@bins_Bins_maxY-@maxY.png", + pdcsapbinningDataO, plotFiltersO, PDCSAPlabelListO, PDCSAPcolorListO, f"Flare peaks per phase of {', '.join(combo)} type stars ({numStars} stars)", f"{locFolder}/{''.join(combo)}_minFlarePeak_{maxFlarePeak}-Flarepeaks_maxY-@maxY.png") # all flare peaks finalDataAllFlarePeaks = pd.DataFrame() if("M" in combo): diff --git a/main/astrodatagui/AstrodataGUI.py b/main/astrodatagui/AstrodataGUI.py index fe01efc..1ddb56a 100644 --- a/main/astrodatagui/AstrodataGUI.py +++ b/main/astrodatagui/AstrodataGUI.py @@ -39,10 +39,12 @@ class AstrodataGUI(QtWidgets.QMainWindow): self.cbPlotRemoveOutliersEnable.stateChanged.connect(self.plotOptionsCBChecked) self.cbPlotRemoveNansEnable.stateChanged.connect(self.plotOptionsCBChecked) self.cbPlotBinEnable.stateChanged.connect(self.plotOptionsCBChecked) - self.cbPlotFoldOptimize.stateChanged.connect(self.plotOptionsCBChecked) self.cbPlotFlattenPlotEnable.stateChanged.connect(self.plotOptionsExclusiveCBChecked) self.cbPlotFoldEnable.stateChanged.connect(self.plotOptionsExclusiveCBChecked) self.cbPlotPeriodogramEnable.stateChanged.connect(self.plotOptionsExclusiveCBChecked) + self.cbPlotFoldShowSpotModulation.stateChanged.connect(self.plotOptionsCBChecked) + + self.gbPlotFoldOptimize.toggled.connect(self.plotOptionsCBChecked) self.comboPlotFluxType.currentTextChanged.connect(self.plotOptionsComboBoxTextChanged) self.comboPlotNormalizeUnit.currentTextChanged.connect(self.plotOptionsComboBoxTextChanged) @@ -252,7 +254,10 @@ class AstrodataGUI(QtWidgets.QMainWindow): self.flaredetectorPreview.setFoldedFitType(self.foldedFitType) def updatePeriods(self, periods: list): - self.edPlotFoldPeriod.setText(str(periods[0].value)) + try: + self.edPlotFoldPeriod.setText(str(periods[0].value)) + except: + self.edPlotFoldPeriod.setText(str(periods[0])) def updateEpochPeriod(self, epoch: float): self.edPlotFoldEpochTime.setText(str(epoch)) @@ -282,7 +287,9 @@ class AstrodataGUI(QtWidgets.QMainWindow): case self.cbPlotBinEnable: self.updateFlaredetectionWidgetBin() - case self.cbPlotFoldOptimize: + case self.gbPlotFoldOptimize: + self.updateFlaredetectionWidgetFold() + case self.cbPlotFoldShowSpotModulation: self.updateFlaredetectionWidgetFold() def plotOptionsExclusiveCBChecked(self, state): @@ -399,14 +406,16 @@ class AstrodataGUI(QtWidgets.QMainWindow): def updateFlaredetectionWidgetFold(self): foldEnabled = self.cbPlotFoldEnable.isChecked() - optimizeEnabled = self.cbPlotFoldOptimize.isChecked() + optimizeEnabled = self.gbPlotFoldOptimize.isChecked() + showSpotModulationEnabled = self.cbPlotFoldShowSpotModulation.isChecked() + try: period = float(self.edPlotFoldPeriod.text()) epoch = float(self.edPlotFoldEpochTime.text()) except: print(f"Invalid period/epoch detected, aborting") return - self.flaredetectorPreview.setFoldState(foldEnabled, period, epoch, optimizeEnabled) + self.flaredetectorPreview.setFoldState(foldEnabled, period, epoch, optimizeEnabled, showSpotModulationEnabled) def updateFlaredetectionWidgetPeriodogram(self): periodogramEnabled = self.cbPlotPeriodogramEnable.isChecked() diff --git a/main/astrodatagui/CalcAllFlaresThread.py b/main/astrodatagui/CalcAllFlaresThread.py index c48ae7c..52e5c14 100644 --- a/main/astrodatagui/CalcAllFlaresThread.py +++ b/main/astrodatagui/CalcAllFlaresThread.py @@ -89,47 +89,70 @@ def getFlareCount(filesDict): plt.savefig(f"{starFolder}/{starNameR}_{source}-{sequence}-periodogram-marked_max.png") plt.close() - pdcsapEpochTime = getEpochTime(lc) - pdcsapFoldedLC = lc.fold(period=pdcsapPeakPeriod, epoch_time=pdcsapEpochTime) - pdcsapPhase, pdcsapSineFit, pdcsapFitType = getFoldedBestFit(pdcsapFoldedLC, fitType=filesDict["FitType"]) - pdcsapMinima, pdcsapMaxima = getFoldedFitPeakValley(pdcsapSineFit) - pdcsapminPhasesBounds, pdcsapmaxPhasesBounds = getPhaseRangesNearPeak((pdcsapMinima, pdcsapMaxima), pdcsapPhase, returnPhaseValue=True) - pdcsapFoldedPeaks = [] - pdcsapFoldedPeaksPhasePair = [] - - pdcsapFoldedLC.scatter() + #pdcsapEpochTime = getEpochTime(lc) + #pdcsapFoldedLC = lc.fold(period=pdcsapPeakPeriod, epoch_time=pdcsapEpochTime) + #pdcsapPhase, pdcsapSineFit, pdcsapFitType = getFoldedBestFit(pdcsapFoldedLC, fitType=filesDict["FitType"]) + optimizedFit = getOptimizedFold(lc.normalize(), filesDict["FitType"]) + #pdcsapMinima, pdcsapMaxima = getFoldedFitPeakValley(pdcsapSineFit) + #pdcsapminPhasesBounds, pdcsapmaxPhasesBounds = getPhaseRangesNearPeak((pdcsapMinima, pdcsapMaxima), pdcsapPhase, returnPhaseValue=True) + pdcsapFoldedPeaks = []; pdcsapPeriodFoldedPeaks = [] + pdcsapFoldedPeaksPhasePair = []; pdcsapPeriodFoldedPeaksPhasePair = [] + optimizedFit["foldedLC"].scatter() plt.title(f"{starName} - folded lightcurve") - plt.plot(pdcsapPhase, pdcsapSineFit, color="blue", label=f"{pdcsapFitType}-fit") + plt.plot(optimizedFit["phase"], optimizedFit["fit"], color="blue", label=f"{optimizedFit['fitType']}-fit") for peak in pdcsapPeaks: - cycle, foldedIndex = convertStarndardIndexToFoldedIndex(pdcsapFoldedLC, peak["StandardIndex"]) + cycle, foldedIndex = convertStarndardIndexToFoldedIndex(optimizedFit["foldedLC"], peak["StandardIndex"]) pdcsapFoldedPeaks.append({"Cycle: ": cycle, "Index": foldedIndex}) - pdcsapFoldedPeaksPhasePair.append({"Phase": pdcsapFoldedLC.phase[pdcsapFoldedLC.cycle == cycle][foldedIndex], "Peak": peak}) - plt.plot(pdcsapFoldedLC.phase[pdcsapFoldedLC.cycle == cycle][foldedIndex].value, - pdcsapFoldedLC.flux[pdcsapFoldedLC.cycle == cycle][foldedIndex], "x", color="red") + pdcsapFoldedPeaksPhasePair.append({"Phase": optimizedFit["foldedLC"].phase[optimizedFit["foldedLC"].cycle == cycle][foldedIndex], "Peak": peak}) + plt.plot(optimizedFit["foldedLC"].phase[optimizedFit["foldedLC"].cycle == cycle][foldedIndex].value, + optimizedFit["foldedLC"].flux[optimizedFit["foldedLC"].cycle == cycle][foldedIndex], "x", color="red") plt.plot([], [], "x", color="red", label="Flare peaks") plt.legend() plt.savefig(f"{starFolder}/{starNameR}_{source}-{sequence}-foldedLC-marked_fit_flares.png") plt.close() + + if(optimizedFit["periodFoldedLC"] is not None): + optimizedFit["periodFoldedLC"].scatter() + plt.title(f"{starName} - folded lightcurve") + plt.plot(optimizedFit["periodFoldedPhase"], optimizedFit["periodFoldedFit"], color="blue", label=f"{optimizedFit['fitType']}-fit") + for peak in pdcsapPeaks: + cycle, foldedIndex = convertStarndardIndexToFoldedIndex(optimizedFit["periodFoldedLC"], peak["StandardIndex"]) + pdcsapPeriodFoldedPeaks.append({"Cycle: ": cycle, "Index": foldedIndex}) + pdcsapPeriodFoldedPeaksPhasePair.append({"Phase": optimizedFit["periodFoldedLC"].phase[optimizedFit["periodFoldedLC"].cycle == cycle][foldedIndex], "Peak": peak}) + plt.plot(optimizedFit["periodFoldedLC"].phase[optimizedFit["periodFoldedLC"].cycle == cycle][foldedIndex].value, + optimizedFit["periodFoldedLC"].flux[optimizedFit["periodFoldedLC"].cycle == cycle][foldedIndex], "x", color="red") + plt.plot([], [], "x", color="red", label="Flare peaks") + plt.legend() + plt.savefig(f"{starFolder}/{starNameR}_{source}-{sequence}-periodFoldedLC-marked_fit_flares.png") + plt.close() filesDict["pdcsapPeaks"] = pdcsapPeaks filesDict["pdcsapPeaksCount"] = len(pdcsapPeaks) filesDict["pdcsapFits"] = pdcsapFits filesDict["pdcsapValidTimespans"] = pdcsapTds filesDict["pdcsapValidSeconds"] = pdcsapValSec - filesDict["pdcsapFoldedCycle"] = sapFoldedLC.cycle - #filesDict["pdcsapFoldedPhase"] = sapFoldedLC.phase.value - #filesDict["pdcsapFoldedFitPhase"] = pdcsapPhase - filesDict["pdcsapFoldedFitPhaseStarEnd"] = (pdcsapFoldedLC.phase.value[0], pdcsapFoldedLC.phase.value[-1]) + filesDict["pdcsapFoldedEpoch"] = optimizedFit["epoch"] + filesDict["pdcsapFoldedCycle"] = optimizedFit["foldedLC"].cycle + filesDict["pdcsapFoldedPhase"] = optimizedFit["foldedLC"].phase.value + filesDict["pdcsapFoldedFitPhase"] = optimizedFit["phase"] + filesDict["pdcsapFoldedFit"] = optimizedFit["fit"] + filesDict["pdcsapFoldedFitPhaseStarEnd"] = (optimizedFit["foldedLC"].phase.value[0], optimizedFit["foldedLC"].phase.value[-1]) filesDict["pdcsapFoldedPeaksPhasePair"] = pdcsapFoldedPeaksPhasePair - filesDict["pdcsapPeriod"] = pdcsapPeakPeriod.value - filesDict["pdcsapPeriodMinima"] = pdcsapMinima - filesDict["pdcsapPeriodMinimaBoundaries"] = pdcsapminPhasesBounds - filesDict["pdcsapPeriodMaxima"] = pdcsapMaxima - filesDict["pdcsapPeriodMaximaBoundaries"] = pdcsapmaxPhasesBounds + filesDict["pdcsapPeriod"] = optimizedFit["Period"] + filesDict["pdcsapSpotModulation"] = optimizedFit["SpotModulation"] + filesDict['FitType'] = optimizedFit["fitType"] + + filesDict["pdcsapPeriodFoldedCycle"] = optimizedFit["periodFoldedLC"].cycle if optimizedFit["periodFoldedLC"] is not None else None + filesDict["pdcsapPeriodFoldedPhase"] = optimizedFit["periodFoldedLC"].phase.value if optimizedFit["periodFoldedLC"] is not None else None + filesDict["pdcsapPeriodFoldedFitPhase"] = optimizedFit["periodFoldedPhase"] + filesDict["pdcsapPeriodFoldedFit"] = optimizedFit["periodFoldedFit"] + filesDict["pdcsapPeriodFoldedFitPhaseStarEnd"] = (optimizedFit["periodFoldedLC"].phase.value[0], optimizedFit["periodFoldedLC"].phase.value[-1]) if optimizedFit["periodFoldedLC"] is not None else (None, None) + filesDict["pdcsapPeriodFoldedPeaksPhasePair"] = pdcsapPeriodFoldedPeaksPhasePair + filesDict['periodFitType'] = optimizedFit["periodFitType"] csvFile = open(f"{starFolder}/{starNameR}_{source}-{sequence}.csv", "a") csvFile.write("StarName,Spectral Type,Rotational Velocity,Rotenional Velocity Unit,Distance,Distance Unit,Source,Sequence,File Path,Initial folded Fit Type,Used folded Fit Type") - csvFile.write(f"{filesDict['StarName']},{filesDict['SpType']},{filesDict['RotVel']},{filesDict['RotVelUnit']},{filesDict['Distance']},{filesDict['DistanceUnit']},{filesDict['Source']},{filesDict['Sequence']},{filesDict['FilePath']},{filesDict['FitType']},{pdcsapFitType}") + csvFile.write(f"{filesDict['StarName']},{filesDict['SpType']},{filesDict['RotVel']},{filesDict['RotVelUnit']},{filesDict['Distance']},{filesDict['DistanceUnit']},{filesDict['Source']},{filesDict['Sequence']},{filesDict['FilePath']},{filesDict['FitType']},{optimizedFit['fitType']}") csvFile.close() del lc print(f"Finished {filesDict['StarName']}, {filesDict['Sequence']}") diff --git a/main/astrodatagui/astrodatagui.ui b/main/astrodatagui/astrodatagui.ui index fc7ae16..7561dba 100644 --- a/main/astrodatagui/astrodatagui.ui +++ b/main/astrodatagui/astrodatagui.ui @@ -174,6 +174,21 @@ Sequences/Target Table ID + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + @@ -255,12 +270,39 @@ Plot options + + 0 + + + 9 + + + 9 + + + 9 + + + 9 + + + 9 + + + 9 + + + 9 + + + 9 + @@ -346,6 +388,21 @@ Normalize + + 9 + + + 9 + + + 9 + + + 9 + + + 6 + @@ -412,6 +469,21 @@ Remove Outliers + + 9 + + + 9 + + + 9 + + + 9 + + + 6 + @@ -464,6 +536,21 @@ Remove nans/infs + + 9 + + + 9 + + + 9 + + + 9 + + + 6 + @@ -522,6 +609,21 @@ false + + 9 + + + 0 + + + 9 + + + 9 + + + 6 + @@ -625,6 +727,21 @@ Flatten + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + @@ -706,34 +823,31 @@ false + + 9 + + + 0 + + + 9 + + + 0 + + + 6 + - - - - - 0 - 0 - - - - - 60 - 0 - - - - - 60 - 16777215 - - + + - 1 + Enable - + @@ -758,7 +872,7 @@ - + @@ -783,7 +897,39 @@ - + + + + + + + + + + + + 0 + 0 + + + + + 60 + 0 + + + + + 60 + 16777215 + + + + 1 + + + + @@ -808,32 +954,26 @@ - - - - + + + + Optimize - - - - - - Enable + + true - - - - - - - - - - - - - Optimize: + + false + + + + + Show Spot Modulation + + + + @@ -870,6 +1010,21 @@ Periodogram + + 9 + + + 9 + + + 9 + + + 9 + + + 6 + @@ -980,6 +1135,21 @@ Star infos + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + diff --git a/main/astrodatagui/ui/FlaredetectorWidget.py b/main/astrodatagui/ui/FlaredetectorWidget.py index d924cf5..1c33b79 100644 --- a/main/astrodatagui/ui/FlaredetectorWidget.py +++ b/main/astrodatagui/ui/FlaredetectorWidget.py @@ -51,14 +51,14 @@ class FlaredetectorWidget(QtWidgets.QWidget): self.periodogramInfoGroupBox.setSizePolicy(sp) self.mainLayout.addWidget(self.periodogramInfoGroupBox) - self.fluxType = "sap_flux" - self.fluxErrType = "sap_flux_err" + self.fluxType = "pdcsap_flux" + self.fluxErrType = "pdcsap_flux_err" self.NormalizeState = {"Enabled": False, "Scale": "unscaled"} self.RemoveOutliersState = {"Enabled": False, "Sigma": 5.0} self.RemoveNansState = {"Enabled": False} self.BinState = {"Enabled": False, "Size": None} self.FlattenState = {"Enabled": False, "WindowLength": 101, "PolynomialOrder": 2} - self.FoldState = {"Enabled": False, "Period": 1, "EpochTime": 0} + self.FoldState = {"Enabled": False, "Period": 1, "EpochTime": 0, "Optimize": False, "spotModulation": False} self.PeriodogramState = {"Enabled": False, "Method": "lombscargle", "View": "frequency"} self.ShowQualityState = {"Enabled": False} @@ -156,9 +156,10 @@ class FlaredetectorWidget(QtWidgets.QWidget): self.updateFit() self.updatePlot() - def setFoldState(self, enabled: bool, period: float, epoch: float, optimize: bool): + def setFoldState(self, enabled: bool, period: float, epoch: float, optimize: bool, spotModulation: bool): self.FoldState["Enabled"] = enabled self.FoldState["Optimize"] = optimize + self.FoldState["SpotModulation"] = spotModulation if(period < 0): print(f"Negative period detected, setting default 1") period = 1 @@ -248,7 +249,24 @@ class FlaredetectorWidget(QtWidgets.QWidget): label += " - flattened" if(self.FoldState["Enabled"]): if(self.FoldState["Optimize"]): - pass + optimizedFold = getOptimizedFold(lc.normalize(), self.foldedFitType) + if(self.FoldState["SpotModulation"]): + lc = optimizedFold["foldedLC"] + self.periodsCalculated.emit([optimizedFold["SpotModulation"]]) + foldOptimizePhase = optimizedFold["phase"] + foldOptimizeFit = optimizedFold["fit"] + elif(not self.FoldState["SpotModulation"] and optimizedFold["periodFoldedLC"] is not None): + lc = optimizedFold["periodFoldedLC"] + self.periodsCalculated.emit([optimizedFold["Period"]]) + foldOptimizePhase = optimizedFold["periodFoldedPhase"] + foldOptimizeFit = optimizedFold["periodFoldedFit"] + else: + lc = optimizedFold["foldedLC"] + self.periodsCalculated.emit([optimizedFold["SpotModulation"]]) + foldOptimizePhase = optimizedFold["phase"] + foldOptimizeFit = optimizedFold["fit"] + print(self.FoldState["SpotModulation"]) + self.epochCalculated.emit(optimizedFold["epoch"]) else: lc = lc.fold(period=self.FoldState["Period"], epoch_time=self.FoldState["EpochTime"]) @@ -275,15 +293,19 @@ class FlaredetectorWidget(QtWidgets.QWidget): cycle, foldedIndex = convertStarndardIndexToFoldedIndex(lc, peak["StandardIndex"]) self.figureAxis.plot(lc.phase[lc.cycle == cycle][foldedIndex].value, lc.flux[lc.cycle == cycle][foldedIndex], "x", color="red") - try: - phase, sineFit, _ = getFoldedBestFit(lc, fitType=self.foldedFitType) - self.figureAxis.plot(phase, sineFit, color="red") - print(phase, sineFit) - minPhasesBoundsIndices, maxPhasesBoundsIndices = getPhaseRangesNearPeak(getFoldedFitPeakValley(sineFit), phase) - plotPhaseRangesNearPeak((minPhasesBoundsIndices, maxPhasesBoundsIndices), phase, ax=self.figureAxis) - except Exception as e: - print("Failed to get fit") - print(e) + if(self.FoldState["Optimize"]): + self.figureAxis.plot(foldOptimizePhase, foldOptimizeFit, color="red") + + else: + try: + phase, sineFit, _ = getFoldedBestFit(lc, fitType=self.foldedFitType) + self.figureAxis.plot(phase, sineFit, color="red") + print(phase, sineFit) + minPhasesBoundsIndices, maxPhasesBoundsIndices = getPhaseRangesNearPeak(getFoldedFitPeakValley(sineFit), phase) + plotPhaseRangesNearPeak((minPhasesBoundsIndices, maxPhasesBoundsIndices), phase, ax=self.figureAxis) + except Exception as e: + print("Failed to get fit") + print(e) elif(self.PeriodogramState["Enabled"]): lc.plot(label=label, ax=self.figureAxis, view=self.PeriodogramState["View"]) if(self.PeriodogramState["View"] == "period"): diff --git a/main/flaredetector/util.py b/main/flaredetector/util.py index 075f385..45f2128 100644 --- a/main/flaredetector/util.py +++ b/main/flaredetector/util.py @@ -1,7 +1,8 @@ import numpy as np from numpy.polynomial.polynomial import Polynomial -from scipy.signal import find_peaks, argrelextrema +from scipy.signal import find_peaks, argrelextrema, periodogram from scipy.optimize import curve_fit +import itertools def findMaxIndices(lc, num=3, distance=100, height=(None, None), sortByHighest=False): if(hasattr(lc, "power")): @@ -392,5 +393,88 @@ def plotPhaseRangesNearPeak(indices, phase, ax=None): def getEpochTime(lc): return lc.time.value[argrelextrema(np.asarray(lc.flux.value), np.less, order=500)[0]][0] -def getOptimizedFold(normalizedLC, period): - pass +def getOptimizedFold(normalizedLC, preferedFoldedFitType): + isValid = False + + periodogramLS = normalizedLC.to_periodogram(method="lombscargle") + periodogramBLS = normalizedLC.to_periodogram(method="boxleastsquares") + + lsPeriods = periodogramLS.period[findMaxIndices(periodogramLS, num=4, distance=100, sortByHighest=True)] + blsPeriods = periodogramBLS.period[findMaxIndices(periodogramBLS, num=4, distance=100, sortByHighest=True)] + + period = -100 + spotModulation = -100 + for lsP, blsP in itertools.product(lsPeriods, blsPeriods): + if(abs(lsP.value - blsP.value) < max(lsP.value, blsP.value)*0.05): + period = np.average([lsP.value, blsP.value]) + print("Period found: ", period) + break + else: + print("No matching LS/BLS peak -> use highest LS") + period = lsPeriods[0].value + + spotModulation = period + epoch = getEpochTime(normalizedLC) + tries = 0 + periodFoldedLC = None + periodFoldedPhase = None + periodFoldedFit = None + periodFitType = None + while(tries < 30): + tries += 1 + foldedLC = normalizedLC.fold(period=spotModulation, epoch_time=epoch) + phase, fit, fitType = getFoldedBestFit(foldedLC, fitType=preferedFoldedFitType) + phaseLength = abs(phase[0]) + abs(phase[-1]) + fitMaximaArgs = argrelextrema(np.asarray(fit), np.greater)[0] + spotModulationBefore = spotModulation + if(len(fitMaximaArgs) > 1): + print("fitMaximaArgs > 1: ", fitMaximaArgs) + if(len(fitMaximaArgs) == 2 and fitMaximaArgs[0] > len(phase)*0.1 and fitMaximaArgs[1] < len(phase)*0.1): + halfPeriod = period/2 + for lsP, blsP in zip(lsPeriods, blsPeriods): + if(abs(lsP.value - halfPeriod) < period*0.05): + spotModulation = lsP.value + print("lsP.value", lsP.value) + break + elif(abs(blsP.value - halfPeriod) < period*0.05): + spotModulation = blsP.value + print("blsP.value", blsP.value) + break + else: + print("2 fit peaks, but no spot modulation") + + if(spotModulationBefore != spotModulation): + periodFoldedLC = foldedLC + periodFoldedPhase = phase + periodFoldedFit = fit + periodFitType = fitType + foldedLC = normalizedLC.fold(period=spotModulation, epoch_time=epoch) + phase, fit, fitType = getFoldedBestFit(foldedLC, fitType=preferedFoldedFitType) + phaseLength = abs(phase[0]) + abs(phase[-1]) + + fitMinimaArgs = argrelextrema(np.asarray(fit), np.less)[0] + fitMinima = phase[0] + for args in fitMinimaArgs: + if(abs(phase[args]) < abs(fitMinima)): + fitMinima = phase[args] + + if(abs(fitMinima) < phaseLength*0.01): + isValid = True + break; + else: + epoch += fitMinima + + return {"Period": period, + "SpotModulation": spotModulation, + "lsPeriods": lsPeriods, + "blsPeriods": blsPeriods, + "foldedLC": foldedLC, + "phase": phase, + "fit": fit, + "periodFoldedLC": periodFoldedLC, + "periodFoldedPhase": periodFoldedPhase, + "periodFoldedFit": periodFoldedFit, + "periodFitType": periodFitType, + "fitType": fitType, + "epoch": epoch, + "isValid": isValid} diff --git a/show_unknown_sptype_stars.py b/show_unknown_sptype_stars.py index 8b8e29d..27e9ce0 100644 --- a/show_unknown_sptype_stars.py +++ b/show_unknown_sptype_stars.py @@ -58,11 +58,11 @@ nonKicStars = pd.DataFrame(nonKicStars) print(kicNames) #print(nonKicStars) -KeplerM = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/M-type-superflares.csv") -KeplerK = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/K-type-superflares.csv") -KeplerG = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/G-type-superflares.csv") -KeplerF = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/F-type-superflares.csv") -KeplerA = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/A-type-superflares.csv") +#KeplerM = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/M-type-superflares.csv") +#KeplerK = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/K-type-superflares.csv") +#KeplerG = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/G-type-superflares.csv") +#KeplerF = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/F-type-superflares.csv") +#KeplerA = pd.read_csv("D:/Masterthesis/ressources/akthukair_AFD/Results/A-type-superflares.csv") #kicBasedSP = [] #for fullKIC in kicNames["kicName"].values: