# kas.py - eerste simulatie van de meetkas op de zonnige aprildag (klok in wintertijd, UTC+1)
import math
import numpy as np

# uurwaarden van 0:00 tot 24:00 (hoofdstuk 4 van Fundamental)
KAS = [9.5, 9.0, 8.6, 8.2, 7.9, 7.7, 7.7, 8.2, 11.0, 15.5, 20.5, 25.0, 28.5, 31.0, 32.0, 30.5, 27.5, 23.5,
       19.5, 16.0, 13.5, 12.0, 11.0, 10.3, 9.8]
BUITEN = [8.5, 8.0, 7.6, 7.2, 6.9, 6.7, 6.6, 7.2, 8.6, 10.4, 12.3, 14.0, 15.4, 16.5, 17.3, 17.8, 17.9, 17.3,
          16.0, 14.3, 12.6, 11.3, 10.3, 9.5, 8.8]
UA = 16.7                                # warmteverlies van de meetkas (W/K)


def zon(t, dag=110, lat=50.8, lon=4.35):
    """Instraling op een horizontaal vlak bij heldere hemel (W/m²); t = kloktijd in wintertijd (h)."""
    d = math.radians(23.44 * math.sin(math.radians(360 / 365 * (284 + dag))))            # declinatie
    B = math.radians(360 / 365 * (dag - 81))
    tv = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)                 # tijdvereffening (min)
    w = math.radians(15 * (t - 1 + lon / 15 + tv / 60 - 12))                            # uurhoek
    s = (math.sin(math.radians(lat)) * math.sin(d)
         + math.cos(math.radians(lat)) * math.cos(d) * math.cos(w))                     # sinus van de zonnehoogte
    if s <= 0:
        return 0.0
    direct = 1361 * 0.7 ** ((1 / max(s, 0.02)) ** 0.678)                                # loodrecht op de stralen
    return direct * s + 0.12 * direct * s + 20 * math.sqrt(s)                           # direct + diffuus


def simuleer(C, a, dt=60.0, T0=KAS[0], P_verw=0.0):
    """Euler voor C dT/dt = a G(t) + P - UA (T - T_buiten(t)), een dag lang."""
    n = round(24 * 3600 / dt)
    t = np.arange(n + 1) * dt / 3600                          # tijd in uren
    T = np.empty(n + 1)
    T[0] = T0
    for k in range(n):
        T_buiten = np.interp(t[k], range(25), BUITEN)
        P = a * zon(t[k]) + P_verw - UA * (T[k] - T_buiten)    # netto warmtestroom (W)
        T[k + 1] = T[k] + dt * P / C                           # Euler-stap
    return t, T


def rms(T, dt=60.0):
    """Wortel uit het gemiddelde kwadraat van (model - meting) op de hele uren."""
    return math.sqrt(np.mean((T[::round(3600 / dt)] - np.array(KAS)) ** 2))


for naam, C, a in (("alle massa, heel het dak ", 248e3, 0.58), ("snelle massa, heel het dak", 30e3, 0.58),
                   ("snelle massa, a = 0,30   ", 30e3, 0.30)):
    t, T = simuleer(C, a)
    print(f"{naam}: max {T.max():5.1f} °C om {t[T.argmax()]:5.2f} h, min {T.min():4.1f} °C, rms {rms(T):4.2f} °C")

print("\nrms (°C)   a = " + "  ".join(f"{a:.2f}" for a in (0.20, 0.25, 0.30, 0.35, 0.40)))
for C in (20e3, 30e3, 40e3, 60e3, 100e3, 150e3, 250e3):
    fouten = [rms(simuleer(C, a, dt=120.0)[1], dt=120.0) for a in (0.20, 0.25, 0.30, 0.35, 0.40)]
    print(f"C = {C / 1000:3.0f} kJ/K  " + "  ".join(f"{f:4.2f}" for f in fouten))