# -*- coding: utf-8 -*-
"""
'Periodensuche für Jahresringe.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


#--------------------------------------------------------------------------
# Aus den EXCEL-Daten wird ein innerhalb von Python verarbeitbarer 
# Dataframe (df) erzeugt.
# Danach werden die Spalten von df mittels der Methode tolist einzeln
# ausgelesen
# zunächst werden nur die Daten aus den Spalten 1 bis 4 (die Radien) 
# eingelesen --> usecols = [1,2,3,4]

daten = "Mappe-Kiefer 2022 - I.xlsx"    # Name der EXCEL-Datendatei
nzeilen = 71                            # Anzahl der einzulesenden Datenzeilen

df = pd.read_excel(daten, skiprows=None, nrows=None,
                   usecols = [1,2,3,4])
y1 = (df["Radius 1"].tolist())
y2 = (df["Radius 2"].tolist())
y3 = (df["Radius 3"].tolist())
y4 = (df["Radius 4"].tolist())

# Aus den eingelesenen Radien werden nun die Radiuszuwächse (Ringdicken)
# berechnet
b=np.array(y1[2:nzeilen+2])
c=np.array(y1[1:nzeilen+1])
y1 = b-c
b=np.array(y2[2:nzeilen+2])
c=np.array(y2[1:nzeilen+1])
y2 = b-c
b=np.array(y3[2:nzeilen+2])
c=np.array(y3[1:nzeilen+1])
y3 = b-c
b=np.array(y4[2:nzeilen+2])
c=np.array(y4[1:nzeilen+1])
y4 = b-c

# Nun werden die Daten aus Spalte 5 (Zeiten) eingelesen
# hierbei werden die ersten 3 (skiprows+1) Zeilen mittels skiprows=2 
# übersprungen
df = pd.read_excel(daten, skiprows=2, nrows=nzeilen,
                   usecols = [5])
t  = (df["zeitliche Mitte [Jahr+Jahresbruchteil]"].tolist())


#--------------------------------------------------------------------------
# zur Reduzierung der "Ausschläge" im Frequenzspektrum
#
mean = statistics.mean(y1)
v1 = np.array(y1)
y1 = v1 - mean 
mean = statistics.mean(y2)
v2 = np.array(y2)
y2 = v2 - mean 
mean = statistics.mean(y3)
v3 = np.array(y3)
y3 = v3 - mean 
mean = statistics.mean(y4)
v4 = np.array(y4)
y4 = v4 - mean 
vmittel = (v1+v2+v3+v4)/4
mean = statistics.mean(vmittel)
ymittel = vmittel - mean


#--------------------------------------------------------------------------
# Fouriertransformation der Testdaten
#
min = 0.02    # kleinste Kreisfrequenz
max = 10      # größte Kreisfrequenz
nout = 10000   # Anzahl der Kreisfrequenzen zwischen min und max
w = np.linspace(min, max, nout)    # w ... Kreisfrequenzen

# Fouriertransformation mittels Lomb-Scargle
# (für ungleichmäßig verteilte Zeitwerte)
# All frequencies in LombScargle are not angular frequencies, but rather 
# frequencies of oscillation (i.e., number of cycles per unit time).
pgram1 = signal.lombscargle(t, y1, w, normalize=True)
pgram2 = signal.lombscargle(t, y2, w, normalize=True)
pgram3 = signal.lombscargle(t, y3, w, normalize=True)
pgram4 = signal.lombscargle(t, y4, w, normalize=True)
pgrammittel= signal.lombscargle(t, ymittel, 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))


# Diagramm 1: Verteilung der Messwerte über der Zeit
ax_t.plot(t, y1, color='green', marker='o', markersize=6, 
          linewidth=0.5, linestyle='dashed', )
ax_t.plot(t, y2, color='red', marker='o', markersize=6, 
          linewidth=0.5, linestyle='dashed', )
ax_t.plot(t, y3, color='yellow', marker='o', markersize=6, 
          linewidth=0.5, linestyle='dashed', )
ax_t.plot(t, y4, color='blue', marker='o', markersize=6, 
          linewidth=0.5, linestyle='dashed', )
ax_t.plot(t, ymittel, color='black', marker='o', markersize=9, 
          linewidth=1, linestyle='dashed', )
ax_t.set_xlabel('Zeit [a]')
ax_t.set_ylabel('Jahresringdicke (Mittelwert bei y=0')


# Diagramm 2: Das normalisierte Periodogramm für die Abtastfrequenzen
ax_w.set_xlim(0.0, max)
ax_w.plot(w, pgram1, color='green')
ax_w.plot(w, pgram2, color='red')
ax_w.plot(w, pgram3, color='yellow')
ax_w.plot(w, pgram4, color='blue')
ax_w.plot(w, pgrammittel, color='black')
ax_w.set_xlabel('Angular frequency [rad/a]')
ax_w.set_ylabel('Normalized amplitude')


# Diagramm 3: Das normalisierte Periodogramm für die Abtastperioden
pmax = (2 * np.pi) / min
ax_p.set_xlim(0.0, pmax)
ax_p.plot(p, pgram1, color='green')
ax_p.plot(p, pgram2, color='red')
ax_p.plot(p, pgram3, color='yellow')
ax_p.plot(p, pgram4, color='blue')
ax_p.plot(p, pgrammittel, color='black',linewidth=3.5)
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(pgrammittel, height=0.05)

ax_p.text(0.7*pmax, 0.15, "Maxima bei:",fontsize=14)
ax_p.text(0.7*pmax, 0.125, "{:.2f} Jahre".format(p[peaks[0]]),
          fontsize=14)
ax_p.text(0.7*pmax, 0.1, "{:.2f} Jahre".format(p[peaks[1]]),
          fontsize=14)
ax_p.text(0.7*pmax, 0.075, "{:.2f} Jahre".format(p[peaks[2]]),
          fontsize=14)
ax_p.text(0.7*pmax, 0.05, "{:.2f} Jahre".format(p[peaks[3]]),
          fontsize=14)
ax_p.text(0.7*pmax, 0.025, "{:.2f} Jahre".format(p[peaks[4]]),
          fontsize=14)

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