# expressie.py — twee gekoppelde toestanden, geen biologisch gekalibreerd model
import numpy as np
# Tijd in h; m en p zijn relatieve hoeveelheden, afzonderlijk gedefinieerd.
def bron(t):
    return 4.0 if 2.0 <= t < 10.0 else 1.0
def f(t,z):
    m,p=z
    return np.array([bron(t)-0.5*m, 2.0*m-0.2*p])
def sim(dt):
    z=np.array([2.0,20.0]); rij=[]
    # Linksconstante bron per tijdstap; stappen liggen precies op schakeltijden.
    for j in range(round(24/dt)+1):
        t=j*dt
        if j%round(.25/dt)==0: rij.append((t,*z))
        s=bron(t)
        def g(a): return np.array([s-.5*a[0],2*a[0]-.2*a[1]])
        k1=g(z); k2=g(z+dt*k1/2); k3=g(z+dt*k2/2); k4=g(z+dt*k3)
        z+=dt*(k1+2*k2+2*k3+k4)/6
    return np.array(rij)
a=sim(.01); b=sim(.005)
print("maximum verschil m/p op kwartierpunten:",np.max(abs(a[:,1:]-b[:,1:])))
for tijd in (0,2,6,10,14,24):
    r=a[round(tijd/.25)]
    print(f"t={tijd:2d} h m={r[1]:.6f} p={r[2]:.6f}")
print("p-piek op kwartierpunten:",a[np.argmax(a[:,2])])