import math
import numpy as np

# g stikstof in drie fictieve pools; geen reactorontwerp of zuurstofgarantie.
mu=.80; kr=1.; m=.10; k=.20; begin=np.array([2.,1.,9.,0.])
def flux(z,licht=1.,E=0.):
    A,D,R,H=z
    U=mu*licht*R/(kr+R)*A
    sterfte=m*A; retour=k*D; oogst=E*A
    return np.array([U-sterfte-oogst,sterfte-retour,retour-U,oogst])
def sim(h=.0025,eind=60.,scenario="licht"):
    n=round(eind/h)
    if not math.isclose(n*h,eind): raise ValueError("niet gehele stapreeks")
    # Alle scenario-overgangen vallen op een staprand.
    if any(not math.isclose(round(v/h)*h,v) for v in [10.,13.,20.]):
        raise ValueError("overgang niet op staprand")
    z=begin.copy(); zs=[z.copy()]
    for j in range(n):
        t=j*h
        licht=0. if scenario=="donker" and 10<=t<13 else 1.
        E=.20 if scenario=="oogst" and 10<=t<20 else 0.
        # Instellingen blijven constant over deze stap; ook bij RK4-eindpunt.
        k1=flux(z,licht,E); k2=flux(z+h*k1/2,licht,E)
        k3=flux(z+h*k2/2,licht,E); k4=flux(z+h*k3,licht,E)
        z=z+h*(k1+2*k2+2*k3+k4)/6
        if min(z)<-1e-10 or not np.all(np.isfinite(z)):
            raise ValueError("niet-fysische pool; geen stilzwijgend afkappen")
        zs.append(z.copy())
    return np.arange(n+1)*h,np.array(zs)
Rstar=m*kr/(mu-m); Astar=(sum(begin)-Rstar)/(1+m/k); Dstar=m*Astar/k
b=Astar*mu*kr/(kr+Rstar)**2
J=np.array([[0,b],[-k-m,-k-b]])
print("evenwicht A,D,R (g N)",Astar,Dstar,Rstar)
print("gereduceerde eigenwaarden per dag",np.linalg.eigvals(J))
scenarios={}
for s in ["licht","donker","oogst"]:
    t,z=sim(scenario=s); scenarios[s]=(t,z)
    drift=np.max(np.abs(z.sum(axis=1)-sum(begin)))
    print(s,"dag60 A,D,R,H",np.round(z[-1],6),"balansfout",f"{drift:.2e}")
z=scenarios["donker"][1]; i10=round(10/.0025); i13=round(13/.0025)
print("A13/A10 donker",z[i13,0]/z[i10,0],"exact",math.exp(-m*3))
assert all(np.max(np.abs(z.sum(axis=1)-12))<1e-10 for t,z in scenarios.values())
assert abs(z[i13,0]/z[i10,0]-math.exp(-.3))<1e-10