photochemistry

photochem.py

import numpy as np
from numpy import pi as π

K_1 = 1e-2
k_2 = 1e5
k_3 = 1e-16
t_d = 24 * 60 * 60

def f(c, t):
    k_1 = K_1 * max(0, np.sin(2 * π * t / t_d))
    return np.array(
        [
            k_1 * c[2] - k_2 * c[0],
            k_1 * c[2] - k_3 * c[1] * c[3],
            k_3 * c[1] * c[3] - k_1 * c[2],
            k_2 * c[0] - k_3 * c[1] * c[3],
        ],
    )


def df(c, t):
    k_1 = K_1 * max(0, np.sin(2 * π * t / t_d))
    return np.array(
        [
            [-k_2, 0,            +k_1,  0         ],
            [0,    -k_3 * c[3],  +k_1, -k_3 * c[1]],
            [0,    +k_3 * c[3],  -k_1, +k_3 * c[1]],
            [+k_2, -k_3 * c[3],  0,    -k_3 * c[1]],
        ],
    )


def solve(c_0, dt, num_steps):
    c = np.zeros((num_steps + 1, len(c_0)))
    c[0] = c_0

    I = np.identity(len(c_0))
    for n in range(num_steps):
        t = n * dt
        A = I - dt * df(c[n], t)
        dc_dt = np.linalg.solve(A, f(c[n], t))
        c[n + 1] = c[n] + dt * dc_dt

    return c