# aanloop.py - de rotor versnelt door de eerste eigenfrequentie van de mast
import math

k, m, zeta = 9839.0, 233.5, 0.005          # mast als veer (N/m), meebewegende massa (kg), demping (rotor stil)
c = 2 * zeta * math.sqrt(k * m)           # dempingsconstante (N s/m)
U = 0.05                                  # onbalans van de rotor (kg m): de kracht is U * W**2
w0 = math.sqrt(k / m)


def afgeleide(t, s, T_op, W_eind):
    """s = (x, v, hoek van de rotor); de rotor versnelt lineair tot W_eind in T_op seconden."""
    x, v, hoek = s
    W = W_eind * min(t / T_op, 1.0)
    return (v, (U * W**2 * math.cos(hoek) - c * v - k * x) / m, W)


def aanloop(T_op, n_eind=100, duur=150.0, h=0.005):
    """Heun (hoofdstuk 4): Euler als voorspeller, het gemiddelde van twee hellingen als verbeteraar."""
    W_eind = 2 * math.pi * n_eind / 60
    s, t, spoor = (0.0, 0.0, 0.0), 0.0, []
    while t < duur:
        a = afgeleide(t, s, T_op, W_eind)
        p = tuple(si + h * ai for si, ai in zip(s, a))
        b = afgeleide(t + h, p, T_op, W_eind)
        s = tuple(si + h / 2 * (ai + bi) for si, ai, bi in zip(s, a, b))
        t += h
        spoor.append((t, s[0]))
    return spoor


if __name__ == "__main__":
    print(f"eigenfrequentie {w0 / (2 * math.pi):.3f} Hz = {60 * w0 / (2 * math.pi):.0f} tpm, tijdconstante {1 / (zeta * w0):.0f} s")
    for T_op in (4, 30, 120):
        spoor = aanloop(T_op)
        tmax, xmax = max(spoor, key=lambda p: abs(p[1]))
        print(f"van 0 naar 100 tpm in {T_op:3d} s: grootste uitwijking {1000 * abs(xmax):4.1f} mm (na {tmax:5.1f} s)")
    F = U * w0**2
    print(f"blijft de rotor op {60 * w0 / (2 * math.pi):.0f} tpm draaien: {1000 * F / (c * w0):.1f} mm (evenwichtstoestand)")