import math
import numpy as np

# Afgeronde leswaarden uit e04, alleen voor het testreservoir.
r=.70; K=770000.; N0=1900.; alarm=500000.
def exact(t,N=N0):
    return K/(1+(K/N-1)*math.exp(-r*t))
def euler_log(h,eind=10.,N=N0):
    n=round(eind/h)
    if not math.isclose(n*h,eind): raise ValueError("niet gehele stapreeks")
    ys=[N]
    for _ in range(n):
        N=N+h*r*N*(1-N/K); ys.append(N)
    return np.arange(n+1)*h,np.array(ys)
t_alarm=math.log((K/N0-1)*alarm/(K-alarm))/r
t_half=math.log(K/N0-1)/r
print("exact dag 10",round(exact(10),3),"cellen/mL")
print("alarmtijd",round(t_alarm,6),"dag; buigpunt",round(t_half,6))
rows=[]
for h in [1.,.5,.25,.125,.0625]:
    t,n=euler_log(h); err=n[-1]-exact(10)
    rows.append((h,n[-1],err))
    print("h,N(10),fout",h,round(n[-1],3),round(err,3))
print("groei bij K/2",r*K/4,"cellen/(mL dag)")
for E in [.25,.35,.50,.80]:
    ns=max(0.,K*(1-E/r)); print("oogst E,N*,Y",E,round(ns,3),round(E*ns,3))
assert abs(exact(0)-N0)<1e-10
assert all(rows[j+1][2]>rows[j][2] for j in range(len(rows)-1))