#!/usr/bin/env python3
"""Reconstruct AEROCAT PFR-Langmuir-Hinshelwood conversion curves.

The smooth outputs reproduce the website model_curves.json basis. They are
model outputs derived from reported equations and constants, not measurements.
"""
from __future__ import annotations
import csv
import math
from pathlib import Path

Q_M3_H = 1.428
CATALYST_MASS_G = 3.0
MOLAR_VOLUME_L_MOL = 24.45

PARAMETERS = {
    "n-hexane": {"mw_g_mol": 86.18, "k": 13.7, "ka": 0.0204},
    "n-pentane": {"mw_g_mol": 72.15, "k": 13.4, "ka": 0.00439},
    "benzene": {"mw_g_mol": 78.11, "k": 6.91, "ka": 0.0752},
    "toluene": {"mw_g_mol": 92.14, "k": 11.5, "ka": 0.0760},
}
C0_PPMV = [1, 2.5, 5, 7.5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100, 105]

def ppmv_to_mg_m3(ppmv: float, mw_g_mol: float) -> float:
    return ppmv * mw_g_mol / MOLAR_VOLUME_L_MOL

def outlet_concentration(c0: float, k: float, ka: float) -> float:
    """Solve ln(c0/c) + Ka(c0-c) = k*Ka*m/Q by bisection."""
    target = k * ka * CATALYST_MASS_G / Q_M3_H
    lo, hi = 1e-12, c0
    for _ in range(200):
        c = (lo + hi) / 2
        residual = math.log(c0 / c) + ka * (c0 - c) - target
        if residual > 0:
            lo = c
        else:
            hi = c
    return (lo + hi) / 2

def conversion_pct(c0_ppmv: float, pollutant: str) -> float:
    p = PARAMETERS[pollutant]
    c0 = ppmv_to_mg_m3(c0_ppmv, p["mw_g_mol"])
    c = outlet_concentration(c0, p["k"], p["ka"])
    return 100 * (1 - c / c0)

def main() -> None:
    out = Path(__file__).with_name("reconstructed_model_curves.csv")
    pollutants = list(PARAMETERS)
    with out.open("w", newline="", encoding="utf-8") as fh:
        writer = csv.writer(fh)
        writer.writerow(["c0_ppmv", *[f"{p}_conversion_pct" for p in pollutants]])
        for c0 in C0_PPMV:
            writer.writerow([c0, *[round(conversion_pct(c0, p), 3) for p in pollutants]])
    print(out)

if __name__ == "__main__":
    main()
