# time constant and gain of a first-order process
tau = 5.0        # time constant, min
K = 1.0          # steady-state gain
print(tau, K)    # print shows the values

t = 10.0                         # min
ratio = t / tau                  # how many time constants?
print(f"t/tau = {ratio}")
print(f"t/tau = {ratio:.2f}")    # two decimal places
print(f"tau squared = {tau**2}") # ** is a power

times = [0, 5, 10, 15]           # a list of times, min
for x in times:                  # x takes each value in turn
    print(x, x / tau)            # indented = inside the loop
# times / tau   # <- uncomment: TypeError. Lists can't do math.

import math                          # standard-library math

def step_resp(t, tau, K=1.0):
    # unit-step response of a first-order system
    return K * (1 - math.exp(-t / tau))

print(step_resp(5.0, tau=5.0))       # t = tau  ->  63.2 %

# step_resp([0, 5, 10], tau=5.0)   # <- uncomment: TypeError
for x in [0, 5, 10]:               # workable, but clumsy
    print(x, round(step_resp(x, 5.0), 3))

import numpy as np

t = np.array([0, 5, 10, 15])   # array of times, min
print(t / 5)                   # math applies to every element
print(t.shape)                 # (4,) -> four numbers

print(np.linspace(0, 25, 6))    # 6 points from 0 to 25
print(np.arange(0, 25, 5))      # step 5, stops before 25
t = np.linspace(0, 25, 251)     # fine grid: dt = 0.1 min
print(len(t), t[:3])            # count, first three values

tau = 5.0
y = 1 - np.exp(-t / tau)           # step response, all 251 points
print(np.round(y[:5], 4))          # peek at the first five
print(f"y at t = 25: {y[-1]:.4f}") # index -1 is the last element

t5 = np.linspace(0, 25, 6)         # 0, 5, 10, ...
y5 = 1 - np.exp(-t5 / 5)
print(y5[1])                       # t = 5 = tau: 63.2 %
print(y5[1:4])                     # elements 1, 2, 3
print(y5[::2])                     # every second element

mask = y >= 0.98                   # True once we are within 2 %
i = np.argmax(mask)                # first index where True
print(f"98% at t = {t[i]:.1f} min  (4 tau = {4*tau} min)")

A = np.array([[0, 1],
              [-2, -3]])           # 2 x 2 matrix
x = np.array([1, 0])
print(A @ x)                       # matrix-vector product
print(np.linalg.eigvals(A))        # poles at -1 and -2

import matplotlib.pyplot as plt

plt.plot(t, y)               # y = 1 - exp(-t/5) from before
plt.show()                   # display the figure

plt.rcParams.update({"font.size": 15})   # bigger text everywhere

plt.plot(t, y)
plt.xlabel("t (min)")
plt.ylabel("y (fraction of final value)")
plt.title("First-order step response, tau = 5 min")
plt.grid(True)
plt.show()

fig, ax = plt.subplots(figsize=(7, 4.5))
for tau_i in [2, 5, 10]:                    # three time constants
    ax.plot(t, 1 - np.exp(-t / tau_i), label=f"tau = {tau_i}")
ax.set_xlabel("t (min)"); ax.set_ylabel("y")
ax.grid(True); ax.legend()
plt.show()

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(t, 1 - np.exp(-t / 5))
ax.axhline(0.632, ls="--", color="k")      # 63.2 % level
ax.axvline(5, ls="--", color="k")          # t = tau
ax.set_xlabel("t (min)"); ax.set_ylabel("y")
ax.grid(True)
plt.show()

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(7, 5))
ax1.step(t, np.ones_like(t), where="post")     # unit step input
ax1.set_ylabel("x (input)")
ax2.plot(t, 1 - np.exp(-t / 5))                # response
ax2.set_ylabel("y (output)"); ax2.set_xlabel("t (min)")
plt.show()

import sympy as sp

s, t = sp.symbols("s t")        # symbols, not numbers
G = 1 / (5*s + 1)               # a transfer function
G                               # last line is displayed

den = s**2 + 3*s + 2
print(sp.factor(den))                       # (s + 1)(s + 2)
print(sp.expand((s + 1) * (s + 2)))         # back again
print(sp.solve(den, s))                     # roots: the poles

y_t = 1 - sp.exp(-t / 5)                    # step response
print(sp.diff(y_t, t))                      # slope of the curve
print(sp.diff(y_t, t).subs(t, 0))           # initial slope = 1/tau
print(sp.integrate(y_t, (t, 0, 5)))         # area from 0 to 5

f = sp.exp(-t / 5)
F = sp.laplace_transform(f, t, s, noconds=True)
F                                           # 1 / (s + 1/5)

Y = 1 / (s * (5*s + 1))            # unit step through 1/(5s+1)
print(sp.apart(Y, s))              # partial fractions
sp.inverse_laplace_transform(Y, s, t)

C = sp.Function("C")
ode = sp.Eq(5 * C(t).diff(t) + C(t), 2)
sol = sp.dsolve(ode, C(t), ics={C(0): 3})
sol                                # C(t) = 2 + exp(-t/5)

C_exact = sp.lambdify(t, sol.rhs, "numpy")   # symbolic -> function
tt = np.array([0, 5, 10, 25])
print(C_exact(tt))                           # 3.0, 2.37, 2.14, 2.007

from scipy.integrate import solve_ivp

def mixing(t, C):                      # C is a list of states
    return [(2 - C[0]) / 5]            # dC/dt

tt = np.linspace(0, 25, 251)
sol_n = solve_ivp(mixing, [0, 25], [3.0], t_eval=tt)
print(sol_n.y[0][:3])                  # first three values

C_num = sol_n.y[0]
C_ex = 2 + np.exp(-tt / 5)             # from SymPy earlier
err = np.abs(C_num - C_ex).max()
print(f"max |numerical - exact| = {err:.1e} g/L")

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(tt, C_ex, label="exact")
ax.plot(tt[::15], C_num[::15], "o", label="solve_ivp")
ax.set_xlabel("t (min)"); ax.set_ylabel("C_A (g/L)")
ax.grid(True); ax.legend()
plt.show()

def first_order(t, y, tau, K, x):
    return [(K * x - y[0]) / tau]      # tau y' + y = K x

r = solve_ivp(first_order, [0, 25], [0.0], t_eval=tt,
              args=(5.0, 2.0, 1.0))    # tau=5, K=2, step of 1
print(f"final value: {r.y[0][-1]:.3f}  (K*x = 2, after 5 tau)")

from scipy.optimize import curve_fit

rng = np.random.default_rng(0)               # repeatable noise
y_meas = 2*(1 - np.exp(-tt/5)) + rng.normal(0, 0.05, tt.size)

def model(t, K, tau):
    return K * (1 - np.exp(-t / tau))

(K_fit, tau_fit), _ = curve_fit(model, tt, y_meas, p0=[1, 1])
print(f"K = {K_fit:.2f} (true 2),  tau = {tau_fit:.2f} (true 5)")

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(tt, y_meas, ".", alpha=0.5, label="data")
ax.plot(tt, model(tt, K_fit, tau_fit), lw=2, label="fit")
ax.set_xlabel("t (min)"); ax.set_ylabel("y")
ax.grid(True); ax.legend()
plt.show()

# pip install control   (Colab only)
import control as ct

G = ct.tf([1], [5, 1])       # G(s) = 1 / (5 s + 1)
print(G)

print("poles:", ct.poles(G))          # -0.2  ->  tau = 5
print("zeros:", ct.zeros(G))          # none
print("gain :", ct.dcgain(G))         # G(0) = 1

t = np.linspace(0, 25, 251)
resp = ct.step_response(G, t)
print(resp.outputs[:3])                       # first three values
print(np.abs(resp.outputs - (1 - np.exp(-t/5))).max())

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(resp.time, resp.outputs)
ax.set_xlabel("t (min)"); ax.set_ylabel("y")
ax.grid(True)
plt.show()

G1 = ct.tf([1], [2, 1])         # tank 1, tau = 2
G2 = ct.tf([1], [3, 1])         # tank 2, tau = 3
G12 = G1 * G2                   # series connection
print(G12)
print("poles:", ct.poles(G12))  # -0.5 and -0.33

fig, ax = plt.subplots(figsize=(7, 4.5))
for zeta in [0.3, 1.0, 2.0]:
    G2nd = ct.tf([1], [1, 2*zeta, 1])            # tau = 1
    r = ct.step_response(G2nd, np.linspace(0, 12, 400))
    ax.plot(r.time, r.outputs, label=f"zeta = {zeta}")
ax.set_xlabel("t"); ax.set_ylabel("y"); ax.grid(True); ax.legend()
plt.show()

t = np.linspace(0, 30, 601)
u = np.sin(1.0 * t)                          # input, 1 rad/min
out = ct.forced_response(G, t, u)            # G = 1/(5s+1)

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(t, u, "--", label="input")
ax.plot(t, out.outputs, label="output")
ax.set_xlabel("t (min)"); ax.grid(True); ax.legend()
plt.show()

tau, K = 5.0, 2.0
t = np.linspace(0, 25, 251)

exact = K * (1 - np.exp(-t / tau))
ode = solve_ivp(lambda t, y: [(K - y[0]) / tau],
                [0, 25], [0.0], t_eval=t).y[0]
tf_out = ct.step_response(ct.tf([K], [tau, 1]), t).outputs

# largest differences: small numbers mean the routes agree
print(np.abs(exact - ode).max(), np.abs(exact - tf_out).max())