# -*- coding: utf-8 -*-
"""
'Periodensuche für Sonnenflecke-Relativzahlen.py'
Created on Wed Apr  6 10:57:21 2022
@author: fischer
"""

#--------------------------------------------------------------------------
# Zugriff auf vorgefertigte Programme
import numpy as np               # Rechnen mit Listen
import matplotlib.pyplot as plt  # Erstellung von Grafiken
import scipy.signal as signal    # hier zur Fouriertransformation
import statistics                # Statistik, hier: Mittelwertbildung
import pandas as pd              # Einlesen von EXCEL-Daten
from scipy.signal import find_peaks  # Suche nach Maxima in einer Liste

#--------------------------------------------------------------------------
# Die folgenden Testdaten - hier Sonnenflecken-Relativzahlen, die durch Schüler
# ermittelt wurden - können genutzt werden, um das Programm zur Periodensuche
# kennenzulernen.
t = [1995.93,1996.43,1996.93,1997.43,1997.93,1998.43,1998.93,1999.43,1999.93,
     2000.43,2000.93,2001.43,2001.93,2002.43,2002.93,2003.43,2003.93,
     2004.43,2004.93,2005.43,2005.93,2006.43,2006.93,2007.43,2007.93,
     2008.43,2008.93,2009.43,2009.93,2010.43,2010.93,2011.43,2011.93,
     2012.43,2012.93,2013.43,2013.93,2014.43,2014.93,2015.43,2015.93,
     2016.43,2016.93,2017.43,2017.93,2018.43,2018.93,2019.43,2019.93,
     2020.43,2020.93,2021.43]

y = [49.5,6.5,13,16,31,28,49.5,61,111,59,81.5,66.5,108.5,63.5,75,45.5,58.5,
     51,19.5,52,41.5,10,17.5,5,10.5,9.5,0,10,24.5,6.5,23,49,60,60.5,25,44.5,
     51,41.5,104,26,29.5,28.5,24.5,1.5,0,0,0,0,0,0,18.5,10.5] 

#--------------------------------------------------------------------------
# zur Reduzierung der "Ausschläge" im Frequenzspektrum
mean = statistics.mean(y)
v = np.array(y)
y = v - mean 

#--------------------------------------------------------------------------
# Fouriertransformation der Testdaten
min = 0.05    # kleinste Kreisfrequenz
max = 10      # größte Kreisfrequenz
nout = 1000   # Anzahl der Kreisfrequenzen zwischen min und max
w = np.linspace(min, max, nout)    # w ... Kreisfrequenzen
pgram = signal.lombscargle(t, y, w, normalize=True)
# p ... Umwandlung in Perioden
p = (2 * np.pi) / w

#--------------------------------------------------------------------------
# Vorbereitung der Diagrammausgabe (Größe und Teilbereiche)
fig, (ax_t, ax_w, ax_p) = plt.subplots(3, 1, #constrained_layout=True, 
     figsize=(16,15))

# Diagrammausgabe 1: Verteilung der Messwerte über der Zeit
ax_t.plot(t, y, color='green', marker='o', markersize=12, 
          linewidth=0.5, linestyle='dashed', )
ax_t.set_xlabel('Zeit [a]')
ax_t.set_ylabel('Sonnenfleckenrelativzahl (Mittelwert bei y=0')

# Diagrammausgabe 2: Das normalisierte Periodogramm für die Abtastfrequenzen
ax_w.set_xlim(0.0, max)
ax_w.plot(w, pgram, color='green')
ax_w.set_xlabel('Angular frequency [rad/a]')
ax_w.set_ylabel('Normalized amplitude')

# Diagrammausgabe3: Das normalisierte Periodogramm für die Abtastperioden
pmax = (2 * np.pi) / min
ax_p.set_xlim(0.0, pmax)
ax_p.plot(p, pgram, color='green')
ax_p.set_xlabel('Periode [a]')
ax_p.set_ylabel('Normalized amplitude')
ax_p.grid(which='minor', alpha=0.1)
ax_p.grid(which='major', alpha=0.5)

#--------------------------------------------------------------------------
# Suche nach den Maxima in der Verteilung der Periodendauern
peaks, _ = find_peaks(pgram, height=0.06)
print(p[peaks])

ax_p.text(25.0, 0.45, "Maxima bei:", fontsize=20)
ax_p.text(35, 0.35, "{:.2f} Jahre".format(p[peaks[0]]),
          fontsize=20)
ax_p.text(35, 0.25, "{:.2f} Jahre".format(p[peaks[1]]),
          fontsize=20)

# Ausgabe der Diagramme in eine png-Datei und auf dem Bildschirm
plt.savefig("Sonnenaktivität.png", dpi=100, bbox_inches='tight', pad_inches = 0)
plt.show()