# module_thermisch.py — twee knopen, SI, modeltijd sinds 00:00 UTC+1
import numpy as np
CA, CB = 80000.0, 1100000.0             # J/K: snelle massa en trage buffer
H, G = 20.2, 30.0                      # W/K: buitenverlies en interne koppeling

def simuleer(dt=10.0, uren=24, vermogen=600.0, buiten=4.0, start=12.0):
    a = b = start
    rij = []
    for k in range(round(uren * 3600 / dt) + 1):
        rij.append((k * dt / 3600, a, b))
        qa = vermogen - H * (a - buiten) - G * (a - b)
        qb = G * (a - b)
        a, b = a + dt * qa / CA, b + dt * qb / CB
    return np.asarray(rij)

A = np.array([[-(H + G)/CA, G/CA], [G/CB, -G/CB]])
polen = np.linalg.eigvals(A)
tau = sorted(-1 / polen / 3600)
rij = simuleer()
for uur in (0, 1, 4, 8, 24):
    k = round(uur * 3600 / 10)
    print(f"{uur:2d} h: lucht {rij[k,1]:.2f} °C, buffer {rij[k,2]:.2f} °C")
print("tijdconstanten (h):", ", ".join(f"{x:.3f}" for x in tau))
fijn = simuleer(dt=5)
print(f"verschil laatste luchtwaarde 10/5 s: {abs(rij[-1,1]-fijn[-1,1]):.5f} K")