# compost.py — gekoppeld onderwijsmodel, geen hygiënevalidatie
import math
import numpy as np
R, Ea = 8.31446261815324, 50000.0

def k_dag(tc):
    arr = 0.08*math.exp(-Ea/R*(1/(tc+273.15)-1/293.15))
    # Fenomenologische hogetemperatuurremming, NIET uit alleen Arrhenius.
    return arr/(1+math.exp((tc-55)/3))

def vat(hverlies, dt=60.0, dagen=10.0):
    c, q, buiten = 180000.0, 12e6, -2.0
    massa, tc = 8.0, 8.0
    rij=[]
    for i in range(round(dagen*86400/dt)+1):
        if i % max(1,round(1800/dt)) == 0:
            rij.append((i*dt/86400,tc,massa))
        # Exacte afname bij tijdens deze stap constant veronderstelde T.
        dm=massa*(-math.expm1(-k_dag(tc)*dt/86400))
        p=q*dm/dt
        tc += dt*(p-hverlies*(tc-buiten))/c
        massa -= dm
    return np.array(rij)

goed, lek = vat(1.0), vat(4.0)
half = vat(1.0,30.0)
for naam,rij in (("H=1 W/K",goed),("H=4 W/K",lek)):
    top=rij[rij[:,1].argmax()]
    print(f"{naam}: top {top[1]:.3f} °C op dag {top[0]:.3f}; "
          f"dag10 {rij[-1,1]:.3f} °C, substraat {rij[-1,2]:.3f} kg")
print(f"tijdstapcontrole top: {goed[:,1].max():.3f} versus {half[:,1].max():.3f} °C")
print(f"max temperatuurverschil op dezelfde halfuursstippen: "
      f"{np.max(np.abs(goed[:,1]-half[:,1])):.5f} K")

print(f"totale vrijgekomen warmte H=1: {12*(8-goed[-1,2]):.3f} MJ")

# Synthetische chemische concentratiemetingen voor eerste-ordefit.
t=np.arange(0.0,1801.0,150.0)
rng=np.random.default_rng(25)
c=0.100*np.exp(-0.0012*t)*np.exp(rng.normal(0,0.015,len(t)))
b,a=np.polyfit(t,np.log(c),1)
fit=np.exp(a+b*t)
print(f"eerste-ordefit k={-b:.7f} 1/s, c0={math.exp(a):.6f} mol/L")
print(f"halveringstijd={math.log(2)/(-b):.2f} s")
print(f"max relatief residu={100*np.max(np.abs((c-fit)/fit)):.2f}%")