# pakket.py - een pakket van 2,0 kg aan een parachute, uit de drone op 40 m (hoofdstuk 4 en 8)
import math
import numpy as np

g, rho, m = 9.81, 1.2, 2.0
k_doos = 0.5 * rho * 1.05 * 0.30**2                 # doos van 30 cm, plat vlak naar voren
D, cd_par = 1.5, 0.8                                # ronde parachute: diameter (m), cd op de nominale oppervlakte
k_par = 0.5 * rho * cd_par * math.pi / 4 * D**2


def open_deel(t, t_start=0.5, t_vol=1.5):
    """Welk deel van de parachute open is: niets tot t_start, dan lineair tot volledig open op t_vol."""
    return min(1.0, max(0.0, (t - t_start) / (t_vol - t_start)))


def val(h0=40.0, vx0=0.0, wind=0.0, parachute=True, dt=0.001):
    """Heun op s = [x, y, vx, vy]; geeft de tijden, de toestanden en de grootste kracht van de lucht."""
    def f(s, t):
        k = k_doos + (k_par * open_deel(t) if parachute else 0.0)
        rx, ry = s[2] - wind, s[3]                   # snelheid tegenover de lucht
        v = math.hypot(rx, ry)
        return np.array([s[2], s[3], -k / m * v * rx, -g - k / m * v * ry]), k * v * v

    s, t, F_max = np.array([0.0, h0, vx0, 0.0]), 0.0, 0.0
    ts, ss = [t], [s]
    while s[1] > 0:
        k1, F = f(s, t)
        k2, _ = f(s + dt * k1, t + dt)
        s, t = s + dt / 2 * (k1 + k2), t + dt
        F_max = max(F_max, F)
        ts.append(t)
        ss.append(s)
    return np.array(ts), np.array(ss), F_max


if __name__ == "__main__":
    ve = math.sqrt(m * g / (k_doos + k_par))
    print(f"eindsnelheid met open parachute: {ve:.2f} m/s")
    for naam, kw in (("zonder parachute", dict(parachute=False)), ("met parachute", {}),
                     ("wind 3 m/s", dict(wind=3.0)), ("wind 5 m/s", dict(wind=5.0)),
                     ("drone 8 m/s vooruit", dict(vx0=8.0))):
        ts, ss, F_max = val(**kw)
        x, y, vx, vy = ss[-1]
        print(f"{naam:20s} {ts[-1]:5.2f} s  landt {x:5.1f} m verder  v = {math.hypot(vx, vy):5.2f} m/s"
              f"  grootste luchtkracht {F_max:5.1f} N")