"""Reimplementation of the conserving cyanide model from Our_Explanation.html.

State: M1, M2 (mobile, two exchange contributions), B (retained product).
  dM_i/dt = Q w_i c / tau_i - M_i/tau_i - (M_i/M) r
  r = Q a (M/Q)^n ; dB/dt = r ; r = 0 if M<=0 or B >= iron cap.
Calibration: choose a so that 24.75 h at c=1 + 71 d at c=0 gives B = 115.34.
"""
import numpy as np, numba, sys, time

W1, W2 = 0.18658, 0.81342
T1, T2 = 7.2, 96.0
BMAX_REF = 11980.20
LAB_TARGET = 115.34
# 10 g HCN/m^3 at 11 C, relative to lab reference 2% by volume (0.02 atm)
C_ROOM = (10.0 / 27.0253) * 8.314462 * 284.15 / 101325.0 / 0.02
DT = 0.01

print(f"C_ROOM = {C_ROOM:.4f}", flush=True)

@numba.njit
def run(a, n, Q, iron_scale, durs, cs):
    M1 = 0.0
    M2 = 0.0
    B = 0.0
    Bmax = BMAX_REF * iron_scale
    Qw1 = Q * W1
    Qw2 = Q * W2
    inv1 = 1.0 / T1
    inv2 = 1.0 / T2
    for p in range(durs.shape[0]):
        steps = int(durs[p] / DT + 0.5)
        c = cs[p]
        s1 = Qw1 * c * inv1
        s2 = Qw2 * c * inv2
        for _ in range(steps):
            M = M1 + M2
            if M > 1e-12 and B < Bmax and a > 0.0:
                r = Q * a * (M / Q) ** n
                f1 = M1 / M
                f2 = M2 / M
            else:
                r = 0.0
                f1 = 0.0
                f2 = 0.0
            M1 += DT * (s1 - M1 * inv1 - f1 * r)
            M2 += DT * (s2 - M2 * inv2 - f2 * r)
            if M1 < 0.0:
                M1 = 0.0
            if M2 < 0.0:
                M2 = 0.0
            B += DT * r
    return B

def calibrate(n, Q):
    durs = np.array([24.75, 1704.0])
    cs = np.array([1.0, 0.0])
    lo, hi = 0.0, 2.0
    while run(hi, n, Q, 1.0, durs, cs) < LAB_TARGET:
        hi *= 2.0
        if hi > 1e6:
            raise RuntimeError("calibration failed")
    for _ in range(22):
        mid = 0.5 * (lo + hi)
        if run(mid, n, Q, 1.0, durs, cs) < LAB_TARGET:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)

def chamber_sched(pulse_h):
    durs, cs = [], []
    for _ in range(400):
        durs += [pulse_h, 24.0 - pulse_h]
        cs += [C_ROOM, 0.0]
    durs += [1704.0]
    cs += [0.0]
    return np.array(durs), np.array(cs)

def delousing_sched():
    durs, cs = [], []
    t = 0.0
    for k in range(400):
        start = k * 16.2
        if start > t:
            durs += [start - t]
            cs += [0.0]
        durs += [6.0]
        cs += [C_ROOM]
        t = start + 6.0
    durs += [6480.0 - t, 1704.0]
    cs += [0.0, 0.0]
    return np.array(durs), np.array(cs)

def scenario(n, Q, pulse_h, ch_sched, del_sched):
    a = calibrate(n, Q)
    B_ch = run(a, n, Q, 1.0, *ch_sched)
    B_del = run(a, n, Q, 1.0, *del_sched)
    return a, B_ch, 0.85 * B_del, 0.90 * B_del

if __name__ == "__main__":
    mode = sys.argv[1] if len(sys.argv) > 1 else "validate"
    # warm up numba
    run(0.1, 2.0, 500.0, 1.0, np.array([1.0]), np.array([1.0]))

    ch_sched = chamber_sched(0.24)
    del_sched = delousing_sched()

    if mode == "validate":
        # (n, Q, expected_a, expected_ch, expected_R12, expected_R13)
        cases = [
            (2.10732997, 500, 0.1669961925, 8.18, 2884, 3053),
            (2.02, 500, None, 11.02, 2953, 3127),
            (1.0, 500, None, 192.99, 4101, 4342),
            (2.0, 500, None, 11.79, 2969, 3144),
            (1.76394234, 500, None, 25.45, 3174, 3361),
            (1.11636664, 500, None, 148.09, 3914, 4144),
            (0.46932662, 500, None, 520.46, 5536, 5861),
            (0.41315569, 500, None, 551.93, 5787, 6127),
            (1.93, 500, None, 14.89, None, None),
            (2.11, 500, None, 8.11, None, None),
            (2.10732997, 200, None, 32.40, None, None),
            (2.10732997, 2000, None, 4.54, None, None),
        ]
        print(f"{'n':>12} {'Q':>6} | {'a':>12} {'a_exp':>12} | {'ch':>9} {'ch_exp':>7} | {'R12':>7} {'R12exp':>6} | {'R13':>7} {'R13exp':>6}")
        t0 = time.time()
        for n, Q, a_e, ch_e, r12_e, r13_e in cases:
            a, B_ch, R12, R13 = scenario(n, Q, 0.24, ch_sched, del_sched)
            ae = f"{a_e:12.6f}" if a_e else " " * 12
            r12e = f"{r12_e:6.0f}" if r12_e else " " * 6
            r13e = f"{r13_e:6.0f}" if r13_e else " " * 6
            print(f"{n:12.8f} {Q:6.0f} | {a:12.6f} {ae} | {B_ch:9.2f} {ch_e:7.2f} | {R12:7.0f} {r12e} | {R13:7.0f} {r13e}", flush=True)
        print(f"elapsed {time.time()-t0:.1f}s", flush=True)

    elif mode == "sweep":
        import csv
        ns = np.linspace(0.4, 2.2, 19)
        Qs = [200, 350, 500, 750, 1000, 1500, 2000]
        widths = [0.24, 0.5, 1.2]
        ch_scheds = {w: chamber_sched(w) for w in widths}
        out = "/home/hatch/workspace/user/files/rudolf_model_sweep/sweep.csv"
        t0 = time.time()
        with open(out, "w", newline="") as f:
            wr = csv.writer(f)
            wr.writerow(["n", "Q", "a"] + [f"ch_{w}" for w in widths] + ["del_unscaled", "R12", "R13"])
            total = len(ns) * len(Qs)
            done = 0
            for n in ns:
                for Q in Qs:
                    a = calibrate(n, Q)
                    chs = [run(a, n, Q, 1.0, *ch_scheds[w]) for w in widths]
                    B_del = run(a, n, Q, 1.0, *del_sched)
                    wr.writerow([f"{n:.6f}", Q, f"{a:.8f}"] +
                                [f"{c:.2f}" for c in chs] +
                                [f"{B_del:.2f}", f"{0.85*B_del:.2f}", f"{0.90*B_del:.2f}"])
                    f.flush()
                    done += 1
                    if done % 20 == 0:
                        print(f"{done}/{total}  {time.time()-t0:.0f}s", flush=True)
        print(f"sweep done in {time.time()-t0:.0f}s -> {out}", flush=True)