"""
Nitric acid — gas-phase oxidation of NO to NO2 (isothermal)
bdrworks.org · Handbook · Nitric acid, worked model 2

    2 NO + O2 -> 2 NO2          third-order, negative temperature coefficient
    2 NO2 <=> N2O4              fast equilibrium
    NO + NO2 <=> N2O3           fast equilibrium

Rate per unit gas volume (O2 consumed), p in atm, r in kmol/(m3 s):
    r = kp/(R T) * (pNO^2 pO2 - pNO2^2 / K_ox),   log10 kp = 652.1/T - 1.0366
The reverse term makes the rate vanish at equilibrium; K_ox comes from NASA polynomials.
A gas parcel of 1 kmol (feed basis) is followed in time at constant T and P:
    d(xi)/dt = kp * n_tot * (pNO^2 pO2 - pNO2^2/K_ox) / P      xi = kmol O2 reacted
Degree of oxidation: alpha = 2 xi / n_NO,0 (N2O4 and N2O3 do not change it).

Run:  python hno3_oxidation.py      (needs numpy, scipy, matplotlib)
"""
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq
import matplotlib.pyplot as plt

RG = 8.314                       # J/(mol K)
ATM = 1.01325                    # bar per atm

# NASA 7-coefficient polynomials (GRI-Mech 3.0): NO, O2, NO2; [T < 1000 K], [T > 1000 K]
NASA_LO = np.array([
    [4.21847630e+00, -4.63897600e-03, 1.10410220e-05, -9.33613540e-09, 2.80357700e-12, 9.84462300e+03, 2.28084640e+00],   # NO
    [3.78245636e+00, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, -1.06394356e+03, 3.65767573e+00],  # O2
    [3.94403120e+00, -1.58542900e-03, 1.66578120e-05, -2.04754260e-08, 7.83505640e-12, 2.89661790e+03, 6.31199170e+00],   # NO2
])
NASA_HI = np.array([
    [3.26060560e+00, 1.19110430e-03, -4.29170480e-07, 6.94576690e-11, -4.03360990e-15, 9.92097460e+03, 6.36930270e+00],   # NO
    [3.28253784e+00, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1.08845772e+03, 5.45323129e+00],  # O2
    [4.88475420e+00, 2.17239560e-03, -8.28069060e-07, 1.57475100e-10, -1.05108950e-14, 2.31649830e+03, -1.17416950e-01],  # NO2
])


def g_RT(T):
    """Dimensionless standard Gibbs energies G/(RT) at 1 bar."""
    a = NASA_LO if T < 1000.0 else NASA_HI
    h = a[:, 0] + a[:, 1]*T/2 + a[:, 2]*T**2/3 + a[:, 3]*T**3/4 + a[:, 4]*T**4/5 + a[:, 5]/T
    s = a[:, 0]*np.log(T) + a[:, 1]*T + a[:, 2]*T**2/2 + a[:, 3]*T**3/3 + a[:, 4]*T**4/4 + a[:, 6]
    return h - s


def K_ox(T):
    """2 NO + O2 = 2 NO2, K = pNO2^2/(pNO^2 pO2) in 1/atm."""
    g = g_RT(T)
    return np.exp(-(2*g[2] - 2*g[0] - g[1])) * ATM


def kp(T):
    return 10.0**(652.1/T - 1.0366)                  # atm^-2 s^-1


def K4(T):
    return 0.698e-9*np.exp(6866.0/T)                 # N2O4/NO2^2, 1/atm


def K3(T):
    return 65.3e-9*np.exp(4740.0/T)                  # N2O3/(NO NO2), 1/atm


def speciate(NOt, NO2t, n_rest, T, P):
    """Free NO, NO2, N2O4, N2O3 (kmol) from the NO and NO2 totals; P in atm."""
    k4, k3 = K4(T)*P, K3(T)*P
    ntot = n_rest + NOt + NO2t
    for _ in range(50):                               # fixed point on the total moles
        def f(b):
            a = NOt/(1.0 + k3*b/ntot)
            return b + 2.0*k4*b*b/ntot + k3*a*b/ntot - NO2t
        b = brentq(f, 0.0, NO2t) if NO2t > 0 else 0.0
        a = NOt/(1.0 + k3*b/ntot)
        c, d = k4*b*b/ntot, k3*a*b/ntot
        new = n_rest + a + b + c + d
        if abs(new - ntot) < 1e-13:
            break
        ntot = new
    return a, b, c, d, new


def rhs(t, y, T, P, n0):
    xi = y[0]
    NOt, O2, NO2t = n0[0] - 2*xi, n0[1] - xi, 2*xi
    a, b, c, d, nt = speciate(max(NOt, 0.0), NO2t, n0[2] + O2, T, P)
    pNO, pO2, pNO2 = a/nt*P, O2/nt*P, b/nt*P
    return [kp(T)*nt*(pNO**2*pO2 - pNO2**2/K_ox(T))/P]


def oxidise(T, P_bar, n0, t_end=1e4):
    P = P_bar/ATM
    stop = lambda t, y, *args: 2*y[0]/n0[0] - 0.9
    stop.terminal = False
    sol = solve_ivp(rhs, (0, t_end), [0.0], args=(T, P, n0), method="LSODA",
                    rtol=1e-8, atol=1e-12, dense_output=True, events=stop)
    t90 = sol.t_events[0][0] if sol.t_events[0].size else np.nan
    return sol, t90


def alpha_eq(T, P_bar, n0):
    """Equilibrium degree of oxidation (monomers only; dimers are negligible when it matters)."""
    P = P_bar/ATM
    def f(al):
        xi = al*n0[0]/2
        n = np.array([n0[0] - 2*xi, n0[1] - xi, 2*xi]); nt = n.sum() + n0[2]
        p = n/nt*P
        return p[2]**2 - K_ox(T)*p[0]**2*p[1]
    return brentq(f, 1e-12, 1 - 1e-12)


def t_analytic(T, P_bar, n0, x):
    """Closed-form isothermal residence time for -dpNO/dt = 2 kp pNO^2 pO2 at constant total moles
    (no reverse reaction, no dimers, no volume contraction)."""
    P = P_bar/ATM
    a, b, k = n0[0]*P, n0[1]*P, 2*kp(T)
    m = 2*b - a
    return 2/(k*m)*(x/(a*(1 - x)) - np.log((2*b - a*x)/(2*b*(1 - x)))/m)


def nitrous_gas(y_nh3=0.10, S=0.965, secondary=0.14):
    """Dry gas after the burner plus secondary air, per kmol (water condensed out).
    secondary = share of the total air sent to the bleacher and added before absorption."""
    prim = 1.0 - y_nh3
    sec = prim*secondary/(1.0 - secondary)
    NO = y_nh3*S
    O2 = 0.2095*(prim + sec) - y_nh3*(1.25*S + 0.75*(1 - S))
    N2 = 0.7905*(prim + sec) + 0.5*y_nh3*(1 - S)
    n = np.array([NO, O2, N2])
    return n/n.sum()


if __name__ == "__main__":
    n0 = nitrous_gas()
    print("dry nitrous gas: NO %.2f %%, O2 %.2f %%, N2 %.2f %%" % tuple(100*n0))
    fig, ax = plt.subplots(1, 3, figsize=(15, 4.4))
    for Tc, c in [(30, "tab:blue"), (50, "tab:cyan"), (100, "tab:green"), (200, "tab:orange"), (350, "tab:red")]:
        sol, t90 = oxidise(Tc + 273.15, 11.0, n0, 100.0)
        t = np.logspace(-3, 2, 300)
        lab = f"t$_{{90}}$ = {t90:.2g} s" if np.isfinite(t90) else f"equilibrium {100*alpha_eq(Tc + 273.15, 11.0, n0):.0f} %"
        ax[0].semilogx(t, 100*2*sol.sol(t)[0]/n0[0], lw=2, color=c, label=f"{Tc} °C  ({lab})")
    ax[0].axhline(90, color="k", ls=":", lw=1)
    ax[0].set(xlabel="residence time (s)", ylabel="degree of oxidation (%)", title="Isothermal oxidation at 11 bar", ylim=(0, 100))
    ax[0].legend(fontsize=8); ax[0].grid(alpha=0.3, which="both")

    Ts = np.arange(20, 401, 10)
    for P, c in [(1.0, "tab:blue"), (4.5, "tab:orange"), (11.0, "tab:red")]:
        t90 = [oxidise(Tc + 273.15, P, n0)[1] for Tc in Ts]
        ax[1].semilogy(Ts, t90, lw=2, color=c, label=f"{P:g} bar")
        k = np.where(np.isfinite(t90))[0][-1]
        ax[1].plot(Ts[k], t90[k], "o", color=c, ms=4)
    ax[1].text(395, 1.5, "curves end where the\nequilibrium falls below 90 %", fontsize=7, ha="right")
    ax[1].set(xlabel="temperature (°C)", ylabel="time to 90 % oxidation (s)", title="Negative temperature coefficient")
    ax[1].legend(); ax[1].grid(alpha=0.3, which="both")

    Te = np.arange(200, 901, 10)
    for P, c in [(1.0, "tab:blue"), (4.5, "tab:orange"), (11.0, "tab:red")]:
        ax[2].plot(Te, [100*alpha_eq(Tc + 273.15, P, n0) for Tc in Te], lw=2, color=c, label=f"{P:g} bar")
    ax[2].axvspan(840, 950, color="0.85", zorder=0, label="gauze")
    ax[2].set(xlabel="temperature (°C)", ylabel="equilibrium degree of oxidation (%)", title="Equilibrium limit of 2 NO + O$_2$ ⇌ 2 NO$_2$")
    ax[2].legend(); ax[2].grid(alpha=0.3)
    fig.tight_layout(); fig.savefig("hno3_oxidation.png", dpi=150)

    print(f"K_ox: 298 K {K_ox(298.15):.3e}, 600 K {K_ox(600):.3e}, 900 K {K_ox(900):.3e} 1/atm")
    print(f"kp: 298 K {kp(298.15):.2f}, 623 K {kp(623.15):.3f} atm^-2 s^-1;  K4(323 K) = {K4(323.15):.3f}, K3(323 K) = {K3(323.15):.3f} 1/atm")
    for P in [1.0, 4.5, 11.0]:
        row = [oxidise(Tc + 273.15, P, n0)[1] for Tc in [30, 50, 100, 200]]
        print(f"{P:5.1f} bar: t90 = " + ", ".join(f"{t:.3g} s" for t in row) + "  (30, 50, 100, 200 °C)")
    for P in [1.0, 4.5, 11.0]:
        print(f"{P:5.1f} bar, 30 °C: closed-form t90 = {t_analytic(303.15, P, n0, 0.9):.3g} s (constant moles, no dimers)")
    for P in [1.0, 4.5, 11.0]:
        print(f"{P:5.1f} bar: equilibrium oxidation 400 °C {100*alpha_eq(673.15, P, n0):.1f} %, 600 °C {100*alpha_eq(873.15, P, n0):.1f} %, 890 °C {100*alpha_eq(1163.15, P, n0):.2f} %")
    a, b, c, d, nt = speciate(0.1*n0[0], 0.9*n0[0], n0[1] - 0.45*n0[0] + n0[2], 323.15, 11.0/ATM)
    print(f"11 bar, 50 °C, 90 % oxidised: NO {a:.4f}, NO2 {b:.4f}, N2O4 {c:.4f}, N2O3 {d:.4f} kmol; N2O4 carries {100*2*c/(0.9*n0[0]):.0f} % of the NO2")
    plt.show()