import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.special import ellipj
# =============================================================================
# 1. PHYSIKALISCHE PARAMETER & INITIALISIERUNG
# =============================================================================
# Traegheitsmomente des asymmetrischen Kreisels
Ix, Iy, Iz = 1.0, 2.0, 3.0
# Anfangsbedingungen (w0_y = 0 fuer den exakten Phasen-Match bei t=0)
w0 = [1.0, 0.0, 0.2]
# Langzeittest: 2 Stunden = 120 Minuten = 7200 Sekunden
t_start = 0
t_end = 7200
t_span = (t_start, t_end)
# Hochaufloesendes Zeitgitter fuer die Auswertung
t_eval = np.linspace(t_start, t_end, 150000)
# =============================================================================
# WEG A: NUMERISCHER SOLVER (Runge-Kutta RK45)
# =============================================================================
def euler_equations(t, w):
wx, wy, wz = w
return [
((Iy - Iz) / Ix) * wy * wz,
((Iz - Ix) / Iy) * wz * wx,
((Ix - Iy) / Iz) * wx * wy
]
print("Starte 2-Stunden-Numerik-Simulation (RK45)...")
sol_num = solve_ivp(euler_equations, t_span, w0, t_eval=t_eval, method='RK45', rtol=1e-8, atol=1e-10)
print("Numerische Integration abgeschlossen.")
# =============================================================================
# WEG B: ANALYTISCHER WEG (Exakte Theorie via elliptische Funktionen)
# =============================================================================
print("Berechne analytische Referenzloesung...")
E_start = 0.5 * (Ix*w0[0]**2 + Iy*w0[1]**2 + Iz*w0[2]**2)
L2 = (Ix*w0[0])**2 + (Iy*w0[1])**2 + (Iz*w0[2])**2
# Exakte Amplitudenmaxima fuer wy
wy_max = np.sqrt((L2 - 2*E_start*Ix) / (Iy * (Iy - Ix)))
# Modul k fuer die Jacobi-Elliptik
k_sq = ((Iy - Ix) * (Iz * 2 * E_start - L2)) / ((Iz - Iy) * (L2 - Ix * 2 * E_start))
# Exakte zeitliche Frequenzskalierung tau
tau_scale = np.sqrt(((Iz - Iy) * (2 * E_start * Iz - L2)) / (Ix * Iy * Iz))
# Berechnung der Jacobischen Amplitudenfunktionen (sn, cn, dn)
sn, cn, dn, _ = ellipj(t_eval * tau_scale, k_sq)
wy_analytisch = wy_max * sn
# =============================================================================
# AUSWERTUNG DER ENERGIEERHALTUNG
# =============================================================================
E_num = 0.5 * (Ix*sol_num.y[0]**2 + Iy*sol_num.y[1]**2 + Iz*sol_num.y[2]**2)
delta_E = E_num - E_start
# =============================================================================
# GRAFISCHE AUSWERTUNG (Zwei separate Plots)
# =============================================================================
print("Generiere Grafiken...")
plt.figure(figsize=(12, 10))
# --- PLOT 1: Das Finale nach 2 Stunden (Die letzten 30 Sekunden) ---
plt.subplot(2, 1, 1)
finale_indices = np.where(sol_num.t >= (t_end - 30))[0]
plt.plot(sol_num.t[finale_indices], sol_num.y[1][finale_indices], 'r-', linewidth=2, label='Numerisch (RK45)')
plt.plot(t_eval[finale_indices], wy_analytisch[finale_indices], 'b--', linewidth=2, label='Analytisch (Exakt)')
plt.title('Euler-Kreisel Langzeittest: Phasen-Vergleich im Finale (Nach 2 Stunden)')
plt.xlabel('Zeit (s)')
plt.ylabel('Winkelgeschwindigkeit wy (rad/s)')
plt.grid(True)
plt.legend(loc='upper right')
# --- PLOT 2: Der Energie-Fehlertrend (Hamilton-Drift) ---
plt.subplot(2, 1, 2)
plt.plot(sol_num.t, delta_E, color='darkred', linewidth=1.5)
plt.title('Numerischer Energiefehler (Drift vom Hamilton-Sollwert)')
plt.xlabel('Zeit (s)')
plt.ylabel('Absolute Energieabweichung Delta E')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.grid(True)
plt.tight_layout()
plt.show()
print("Fertig!")
