# -*- coding: utf-8 -*-
"""
Created on Fri Mar 27 13:40:49 2026

@author: J. Grabow (Feinmechaniker)
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parameter
m = 175.0
cw = 0.78
A = 1.0
h0 = 41000.0
g = 9.81
rho0 = 1.225
gamma = 1.4      # Adiabatenexponent Luft
R_spec = 287.05  # Spezifische Gaskonstante Luft J/(kg*K)

def get_atmo_params(h):
    """Berechnet Luftdichte und Schallgeschwindigkeit nach ISA"""
    # Vereinfachtes Modell der Temperaturschichten
    if h > 20000: # Stratosphäre Teil 2
        T = 216.65 + (h - 20000) * 0.001 
    elif h > 11000: # Stratosphäre Teil 1 (Isotherm)
        T = 216.65
    else: # Troposphäre
        T = 288.15 - h * 0.0065
    
    rho = rho0 * np.exp(-g * h / (R_spec * T)) # Barometrische Höhenformel (lokal)
    cs = np.sqrt(gamma * R_spec * T)
    return rho, cs

def ode_fall(t, y):
    h, v = y
    if h <= 0: return [0, 0]
    rho, _ = get_atmo_params(h)
    drag = 0.5 * rho * v**2 * cw * A * np.sign(v)
    dvdt = -g - (drag / m)
    return [v, dvdt]

def hit_ground(t, y): return y[0]
hit_ground.terminal = True

sol = solve_ivp(ode_fall, (0, 600), [h0, 0.0], events=hit_ground, max_step=0.5)

# Daten für Plots
t = sol.t
h = sol.y[0]
v_ms = np.abs(sol.y[1])
rho_list, cs_list = np.vectorize(get_atmo_params)(h)
mach_no = v_ms / cs_list

# Visualisierung
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

# Plot 1: Höhe über Zeit
ax1.plot(t, h / 1000, color='green', lw=2)
ax1.set_title("Fallhöhe über die Zeit")
ax1.set_ylabel("Höhe [km]")
ax1.grid(True, linestyle='--')

# Plot 2: Geschwindigkeit & Schallmauer
ax2.plot(t, v_ms, label='Fallgeschwindigkeit [m/s]', lw=2)
ax2.plot(t, cs_list, label='Lokale Schallgeschwindigkeit $c_s$', linestyle='--', color='orange')
ax2.fill_between(t, v_ms, cs_list, where=(v_ms > cs_list), color='red', alpha=0.2, label='Überschallbereich')
ax2.set_title("Geschwindigkeit und Schallgrenze")
ax2.set_xlabel("Zeit [s]")
ax2.set_ylabel("Geschwindigkeit [m/s]")
ax2.legend()
ax2.grid(True, linestyle='--')

plt.tight_layout()
plt.show()

print(f"Maximale Mach-Zahl: {np.max(mach_no):.2f}")
print(f"Zeit im freien Fall bis Aufprall: {t[-1]:.1f} s")