# waterraket.py - een PET-fles van 1,5 L met water en lucht onder druk, verticaal gelanceerd
import math

g, rho, p_at, gamma = 9.81, 1000.0, 1.013e5, 1.4
V_fles, d_mond, m_leeg = 1.5e-3, 0.0215, 0.15      # m³, m, kg
A_mond = math.pi / 4 * d_mond**2
k = 0.5 * 1.2 * 0.4 * math.pi / 4 * 0.09**2          # luchtweerstand van de fles (Ø 9 cm, cd 0,4)


def lancering(V_water, p_over=4.0e5, dt=1e-4):
    """Geeft (brandtijd, snelheid bij brandsluiting, hoogste punt). De luchtstoot na het water verwaarlozen we."""
    V_lucht0, p0 = V_fles - V_water, p_over + p_at
    V_w, v, h, t, t_brand, v_brand = V_water, 0.0, 0.0, 0.0, None, None
    while v >= 0.0:
        m = m_leeg + rho * V_w
        F = 0.0
        if V_w > 0:
            p = p0 * (V_lucht0 / (V_fles - V_w))**gamma  # lucht zet adiabatisch uit
            if p > p_at:
                u = math.sqrt(2 * (p - p_at) / rho)      # uitstroomsnelheid (Bernoulli)
                F = rho * A_mond * u * u                 # stuwkracht = massadebiet x uitstroomsnelheid
                V_w = max(V_w - A_mond * u * dt, 0.0)
            else:
                V_w = 0.0                                # de druk is op: resterend water blijft zitten
        elif t_brand is None:
            t_brand, v_brand = t, v
        a = F / m - g - k / m * v * abs(v)
        v, h, t = v + a * dt, h + v * dt, t + dt
    return t_brand, v_brand, h


print("water (L)  brandtijd (s)  v einde stuwing (m/s)  hoogste punt (m)")
for V in (0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0):
    tb, vb, h = lancering(V * 1e-3)
    print(f"{V:6.1f}  {tb:12.3f}  {vb:18.1f}  {h:16.1f}")