# vliegwiel.py - een stalen vliegwiel als buffer voor een wolk over de zonnepanelen
import math

R, dikte, rho = 0.25, 0.05, 7850.0                  # stalen schijf: straal en dikte (m), dichtheid (kg/m³)
m = math.pi * R**2 * dikte * rho
I = 0.5 * m * R**2
w_max, w_min = 6000 * 2 * math.pi / 60, 3000 * 2 * math.pi / 60      # 6000 en 3000 toeren per minuut, in rad/s
eta = 0.90                                          # rendement motor/generator met elektronica, per richting


def E(w):
    return 0.5 * I * w**2


def verlies(w):
    """Stilstandsverlies (W): lagers (evenredig met w) + lucht (met w³); 5 W + 10 W bij 6000 tpm (rekenwaarden)."""
    return 5.0 * (w / w_max) + 10.0 * (w / w_max)**3


def tekort(t):
    """Vermogen (W) dat het vliegwiel moet leveren: een wolk van t = 60 s tot 240 s haalt 500 W weg."""
    return 500.0 if 60.0 <= t < 240.0 else 0.0


def laden(t, w):
    """Na de wolk laadt de zon het vliegwiel weer op met 400 W, tot 6000 tpm."""
    return 400.0 if t >= 240.0 and w < w_max else 0.0


def wolk(dt=0.1, duur=900.0):
    """Euler op I w dw/dt = eta P_laden - P_tekort / eta - P_verlies; geeft tijden en toerentallen (tpm)."""
    w, t, ts, ns = w_max, 0.0, [0.0], [6000.0]
    while t < duur:
        P = eta * laden(t, w) - tekort(t) / eta - verlies(w)
        w = min(w + dt * P / (I * w), w_max)
        t += dt
        ts.append(t)
        ns.append(w * 60 / (2 * math.pi))
    return ts, ns


def leeglopen(uren=24, dt=1.0):
    """Het vliegwiel draait zonder last en verliest alleen door lagers en lucht; per uur het deel van de energie."""
    w, deel = w_max, [1.0]
    for uur in range(uren):
        for _ in range(round(3600 / dt)):
            w -= dt * verlies(w) / (I * w)
        deel.append(E(w) / E(w_max))
    return deel


if __name__ == "__main__":
    print(f"schijf {m:.1f} kg, I = {I:.3f} kg m², vol {E(w_max) / 3600:.0f} Wh, bruikbaar tot 3000 tpm "
          f"{(E(w_max) - E(w_min)) / 3600:.0f} Wh, randsnelheid {w_max * R:.0f} m/s")
    ts, ns = wolk()
    n_laag = min(ns)
    print(f"wolk van 3 minuten (500 W): laagste toerental {n_laag:.0f} tpm, "
          f"{100 * (1 - (n_laag / 6000)**2):.1f} % van de energie gebruikt")
    deel = leeglopen()
    for uur in (1, 2, 5, 10, 24):
        print(f"stilstand, na {uur:2d} h: nog {100 * deel[uur]:5.1f} % van de energie")