import numpy as np
from scipy.linalg import expm

C_a,C_b,H,G = 80000.0,1100000.0,20.2,30.0
A=np.array([[-(H+G)/C_a,G/C_a],[G/C_b,-G/C_b]])
B=np.array([1/C_a,0.0])
evenwicht=np.linalg.solve(-A,B*600)+4
t=np.linspace(0,24*3600,145)
temp=np.array([evenwicht+expm(A*s)@(np.array([12.,12.])-evenwicht) for s in t])
polen=np.linalg.eigvals(A)
print("tijdschalen h",sorted(-1/polen/3600))
print("na24h",*(f"{v:.6f}" for v in temp[-1]))
print("evenwicht",evenwicht)

def stap2(t,zeta,wn=1.0):
    t=np.asarray(t)
    if zeta<1:
        wd=wn*np.sqrt(1-zeta*zeta)
        return 1-np.exp(-zeta*wn*t)*(np.cos(wd*t)+zeta*wn/wd*np.sin(wd*t))
    if zeta==1: return 1-(1+wn*t)*np.exp(-wn*t)
    p1=-wn*(zeta-np.sqrt(zeta*zeta-1))
    p2=-wn*(zeta+np.sqrt(zeta*zeta-1))
    return 1+(p2*np.exp(p1*t)-p1*np.exp(p2*t))/(p1-p2)

zeta=0.4
mp=np.exp(-np.pi*zeta/np.sqrt(1-zeta*zeta))
tp=np.pi/np.sqrt(1-zeta*zeta)
print("zeta0.4 overshootpct",100*mp,"piekbij",tp)
assert abs(temp[-1,0]-28.375973)<1e-6
# Identificatie zonder ruis: bekende15K-evenwichtsverhoging.
ti=np.array([1,2,3,4,5])*3600
yi=15*(1-np.exp(-ti/10800))
helling=np.polyfit(ti,np.log(1-yi/15),1)[0]
print("geidentificeerde tau",-1/helling,"s")