import math
import numpy as np

# Fictief eenknoopslichaam; geen voorspeller van echte kern- of huidtemperatuur.
C=210000.; H=8.; kp=40.; tau=120.; grens=60.; basis=37.
# Signed u is extra warmteafvoer t.o.v. de basis; u<0 betekent compensatie.
def sim(h=10.,geregeld=True,eind=10800.):
    n=round(eind/h)
    if not math.isclose(n*h,eind) or not math.isclose(round(3600/h)*h,3600):
        raise ValueError("eind of overgang niet op staprand")
    z=np.zeros(2); zs=[z.copy()]
    for j in range(n):
        w=50. if j*h<3600 else 0.
        def f(z):
            x,u=z
            doel=np.clip(kp*x,-grens,grens) if geregeld else 0.
            return np.array([(w-H*x-u)/C,(doel-u)/tau])
        k1=f(z);k2=f(z+h*k1/2);k3=f(z+h*k2/2);k4=f(z+h*k3)
        z=z+h*(k1+2*k2+2*k3+k4)/6;zs.append(z.copy())
    return np.arange(n+1)*h,np.array(zs)
t,z=sim(); ongeregeld=sim(geregeld=False)
print("50W gedurende60min, x/u na60min",np.round(z[360],6))
print("zonder u x na60min",ongeregeld[1][360,0])
print("x na180min",z[-1,0])
M=np.array([[-H/C,-1/C],[kp/tau,-1/tau]])
print("eigenwaarden per s",np.linalg.eigvals(M))
print("evenwicht bij constant50W x,u",50/(H+kp),kp*50/(H+kp))
print("bij120W en u op60W: x*",(120-60)/H)
print("10MJ/dag gemiddeld W",10e6/86400)
print("stapverschil 10/5s",np.max(np.abs(sim(10)[1]-sim(5)[1][::2])))
assert abs(ongeregeld[1][360,0]-50/H*(1-math.exp(-H/C*3600)))<1e-9
assert np.max(np.abs(sim(10)[1]-sim(5)[1][::2]))<1e-5