# # Numerische Berechnung der Trajektorien des harmonischen Oszillators

# Die Bewegung eines (gedämpften) harmonischen Oszillators wird durch die DGL zweiter Ordnung
# \begin{equation}
# \frac{d^2 x(t)}{dt^2} + 2 \gamma \frac{d x(t)}{dt} + \omega^2 x(t) = 0\tag{1}
# \end{equation}
# beschrieben. 
# Um die DGL numerisch zu lösen, überführen wir sie zunächst in ein System gekoppelter DGLs erster Ordnung:
# \begin{align}
# \frac{d x}{dt} &= f\tag{2}\\
# \frac{d v}{dt} &= g,\tag{3}
# \end{align}
# wobei $v(t) = \dot{x}(t)$ die Geschwindigkeit ist.
# Ihre erste Aufgabe ist es, diese zwei DGL zu vervollständigen (d.h. $f$ und $g$ sind zu bestimmen).
# 
# Im folgenden sollen die DGLs numerisch mithilfe des Eulerverfahrens gelöst werden. 
# Dazu werden die Ableitungen wie folgt diskretisiert, wobei $\Delta t$ der Wert des Zeitschritts ist, den Sie in der Simulation nutzen:
# \begin{align}
# \frac{d x}{dt} &\approx \frac{x(t+\Delta t) - x(t)}{\Delta t} \tag{4}\\
# \frac{d v}{dt} &\approx \frac{v(t+\Delta t) - v(t)}{\Delta t}\,.\tag{5}
# \end{align}
# 
# Setzt man diese Näherungen der Ableitungen in (2) bzw. (3) ein, erhält man eine Approximation für die Änderung von $x$ bzw. $v$ im aktuellen Zeitschritt durch die Beziehungen
# \begin{align}
# x(t+\Delta t) =&\ x(t) + \Delta t \, f\tag{6}\\
# v(t+\Delta t) =&\ v(t) + \Delta t \, g\tag{7} \, .
# \end{align}
# Sie können den Code gleich für den allgemeinen Fall implementieren. 
# Führen Sie aber, wie im Aufgabentext erläutert, zunächst im a)-Teil die erforderlichen Simulationen für den ungedämpften harmonischen Oszillator ($\gamma = 0$) durch, danach im b)-Teil mit dem angegebenen Wert für $\gamma$.\\
# Berechnen und plotten Sie im Verlauf der Simulationen auch die Gesamtenergie als Funktion der Zeit, $E(t) = T(t) + V(t)$, wobei $T(t)$ die Zeitentwicklung der kinetischen Energie und $V(t)$ die potenzielle Energie als Funktion der Zeit ist. 
# (Sie können annehmen, dass die Masse des Oszillators $m=1$ ist.)

# In[ ]:


# Hilfreiche Pakete
import numpy as np
import matplotlib.pyplot as plt
import math

# Konstanten
from scipy.constants import g

# Sinnvolle Konstanten
x0 = 1   
v0 = 0
omega0 = 1
gamma = 0
dt = 1e-2 # Kleiner Zeitschritt
tmax = 20 # Max time for the algorithm to terminate

# Initial values
x_t = [x0]
x_t_exact = [x0]
v_t = [v0]
v_t_exact = [v0]
T_t = [0.5*v0*v0]
V_t = [0.5*omega0*omega0*x0*x0]
E_t = [T_t+V_t]

# exakter Wert der Energie:
E_exact = E_t[0]
Ediff_t = [0]

# Implementation of the Euler-Algorithm
for _ in range(int(tmax // dt)):
    f = ???
    x = x_t[-1] + f * dt # Calc new position
    v = ??? # calc new velocity
    x_t.append(x)
    v_t.append(v)

    x_exact = ???
    v_exact = ???
    x_t_exact.append(x_exact)
    v_t_exact.append(v_exact)

    T = ???
    V = k/2 * x^2 * 
    E = T + V
    T_t.append(T)
    V_t.append(V)
    E_t.append(E)
    Ediff_t.append(E-E_exact)



# TIPP: Zum erstellen mehrerer Plots auf einmal, siehe z.B.:
# https://matplotlib.org/3.1.1/gallery/subplots_axes_and_figures/subplot.html

kwargs = {'c':'b'}
font_kwargs = {'fontsize':14}
times = np.arange(0,len(x_t))*dt
abs_max = max(x_t, key=abs)

fig,ax =  plt.subplots(1,4,figsize=(10,5))

#actual plots
ax[0].plot(times,x_t,**kwargs)
ax[1].plot(x_t,v_t,**kwargs)
ax[2].plot(times,T_t,label = "T")
ax[2].plot(times,V_t,label = "V")
ax[2].plot(times,E_t,label = "E = T+V")
ax[3].plot(times,Ediff_t,label = "Fehler in der Energie")

#style changes
ax[0].set_xlim(0,len(x_t)*dt)
ax[0].set_ylim(-abs_max,abs_max)
ax[0].set_xlabel("t",**font_kwargs)
ax[0].set_ylabel("$x(t)$",**font_kwargs)

ax[1].set_xlabel("$x(t)$",**font_kwargs)
ax[1].set_ylabel("$v(t)$",**font_kwargs)

ax[2].set_xlim(0,len(x_t)*dt)
ax[2].set_ylim(-abs_max,abs_max)
ax[2].set_xlabel("t",**font_kwargs)
ax[2].set_ylabel("Energies",**font_kwargs)

ax[2].set_xlim(0,len(x_t)*dt)
ax[2].set_ylim(-abs_max,abs_max)
ax[2].set_xlabel("t",**font_kwargs)
ax[2].set_ylabel("Differenz",**font_kwargs)

plt.legend()
plt.show()