# reactoren.py — generieke cultuur, niet een werkelijk algenstammodel
import math
import numpy as np
mumax,K,Y=.4,.5,.5 # h^-1, g/L, g droge X per g substraat
# Toestand: X g/L, S g/L, V L, cumulatief ingevoerd S g,
# afgevoerd X g, afgevoerd S g.
def deriv(z,mode):
    X,S,V,I,HX,HS=z
    mu=mumax*max(S,0)/(K+max(S,0))
    if mode=="batch": fin=fout=0; sin=0
    elif mode=="fed": fin=.05; fout=0; sin=20.0
    else: fin=fout=float(mode)*V; sin=10.0
    return np.array([mu*X-fin/V*X,
        fin/V*(sin-S)-mu*X/Y,
        fin-fout,fin*sin,fout*X,fout*S])
def sim(mode,dt,tend):
    z=np.array([.1,10,2,0,0,0],float); rows=[]
    for j in range(round(tend/dt)+1):
        if j%round(.5/dt)==0: rows.append((j*dt,*z))
        k1=deriv(z,mode); k2=deriv(z+dt*k1/2,mode)
        k3=deriv(z+dt*k2/2,mode); k4=deriv(z+dt*k3,mode)
        z+=dt*(k1+2*k2+2*k3+k4)/6
    return np.array(rows)
runs={}
for mode,tend in [("batch",24),("fed",24),(.2,72),(.4,72)]:
    a=sim(mode,.01,tend); b=sim(mode,.005,tend); runs[mode]=a
    # g X-equivalent: reactor + afvoer - ingevoerde opbrengst.
    sluit=a[:,3]*(a[:,1]+Y*a[:,2])+a[:,5]+Y*a[:,6]-Y*a[:,4]
    print("mode",mode,"eind t/X/S/V/I/HX/HS:",a[-1])
    print("max |dt0.01-dt0.005|:",np.max(abs(a[:,1:]-b[:,1:])))
    print("balansafwijking g X-equivalent:",np.max(abs(sluit-sluit[0])))
    assert np.min(a[:,1:4])>=-1e-10
    assert np.max(abs(sluit-sluit[0]))<1e-7
Dcrit=mumax*10/(K+10)
Dopt=mumax*(1-math.sqrt(K/(K+10)))
Sopt=K*Dopt/(mumax-Dopt); Xopt=Y*(10-Sopt)
print(f"Dcrit={Dcrit:.8f} Dopt={Dopt:.8f} Sopt={Sopt:.8f} Xopt={Xopt:.8f}")
print(f"max productiviteit={Dopt*Xopt:.8f} g/(L h)")