import math
import numpy as np

# Eigen klein-signaalmodel; i is een relatief effect, geen insulinedosis.
a=.020; b=.050; c=.040; d=.080; Gbasis=5.
M=np.array([[-a,-b],[c,-d]]) # minuten als tijdeenheid
omega=math.sqrt(b*c-(d-a)**2/4)
def exact(t):
    g=3*math.exp(-(a+d)*t/2)*(math.cos(omega*t)+(d-a)/(2*omega)*math.sin(omega*t))
    i=3*c/omega*math.exp(-(a+d)*t/2)*math.sin(omega*t)
    return np.array([g,i])
def sim(h=.25,eind=120.,tau=6.,invoer=0.,c_model=c):
    n=round(eind/h)
    if not math.isclose(n*h,eind): raise ValueError("niet gehele stapreeks")
    # g (mmol/L), i (relatieve eenheid), gemeten g (mmol/L).
    z=np.array([3.,0.,3.]); zs=[z.copy()]
    def f(z):
        g,i,gm=z
        return np.array([-a*g-b*i+invoer,c_model*g-d*i,(g-gm)/tau])
    for _ in range(n):
        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(); tabel=[]
for h in [2.,1.,.5,.25]:
    ts,zs=sim(h); err=float(np.max(np.abs(zs[:,:2]-np.array([exact(tt) for tt in ts]))))
    tabel.append((h,err));print("RK4 h,maxfout g/i",h,f"{err:.3e}")
print("eigenwaarden per min",np.linalg.eigvals(M))
for tt in [20.,60.,120.]: print("exact t,g,i",tt,np.round(exact(tt),6))
print("sensor op20 min",z[round(20/.25),2],"werkelijk model-g",exact(20)[0])
const=.06; gs=const/(a+b*c/d); is_=c/d*gs
print("constante invoer g*,i*",gs,is_)
laag=sim(c_model=.010)
print("kleinere c, G op60",Gbasis+laag[1][round(60/.25),0])
assert max(q[1] for q in tabel[2:])<1e-5
assert np.all(Gbasis+z[:,0]>0)