"""
Nitric acid — ammonia burner (catalytic oxidation of NH3 on Pt–Rh gauze)
bdrworks.org · Handbook · Nitric acid, worked model 1

Adiabatic energy balance of the gauze reactor:
    4 NH3 + 5 O2 -> 4 NO + 6 H2O   (selectivity S)
    4 NH3 + 3 O2 -> 2 N2 + 6 H2O   (the rest)
The gauze temperature follows from H_in(T_mix) = H_out(T_gauze); the oxygen balance
shows how much secondary air the absorption section needs.

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

# ---------------------------------------------------------------- data
SPECIES = ["NH3", "O2", "N2", "NO", "H2O"]
RG = 8.314                                   # J/(mol K)
Y_O2_AIR = 0.2095                            # O2 in dry air (Ar lumped with N2)

# NASA 7-coefficient polynomials (GRI-Mech 3.0 thermodynamic data): [T < 1000 K], [T > 1000 K]
NASA_LO = np.array([
    [4.28602740e+00, -4.66052300e-03, 2.17185130e-05, -2.28088870e-08, 8.26380460e-12, -6.74172850e+03, -6.25372770e-01],   # NH3
    [3.78245636e+00, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, -1.06394356e+03, 3.65767573e+00],   # O2
    [3.29867700e+00, 1.40824040e-03, -3.96322200e-06, 5.64151500e-09, -2.44485400e-12, -1.02089990e+03, 3.95037200e+00],   # N2
    [4.21847630e+00, -4.63897600e-03, 1.10410220e-05, -9.33613540e-09, 2.80357700e-12, 9.84462300e+03, 2.28084640e+00],    # NO
    [4.19864056e+00, -2.03643410e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, -3.02937267e+04, -8.49032208e-01],  # H2O
])
NASA_HI = np.array([
    [2.63445210e+00, 5.66625600e-03, -1.72786760e-06, 2.38671610e-10, -1.25787860e-14, -6.54469580e+03, 6.56629280e+00],   # NH3
    [3.28253784e+00, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1.08845772e+03, 5.45323129e+00],   # O2
    [2.92664000e+00, 1.48797680e-03, -5.68476000e-07, 1.00970380e-10, -6.75335100e-15, -9.22797700e+02, 5.98052800e+00],   # N2
    [3.26060560e+00, 1.19110430e-03, -4.29170480e-07, 6.94576690e-11, -4.03360990e-15, 9.92097460e+03, 6.36930270e+00],    # NO
    [3.03399249e+00, 2.17691804e-03, -1.64072518e-07, -9.70419870e-11, 1.68200992e-14, -3.00042971e+04, 4.96677010e+00],   # H2O
])


def h(T):
    """Molar enthalpies including heats of formation, J/mol."""
    a = NASA_LO if T < 1000.0 else NASA_HI
    return RG * T * (a[:, 0] + a[:, 1]*T/2 + a[:, 2]*T**2/3 + a[:, 3]*T**3/4 + a[:, 4]*T**4/5 + a[:, 5]/T)


def feed(y_nh3, hum=0.0):
    """1 mol of NH3-air mixture; hum = mol H2O per mol dry air."""
    air = 1.0 - y_nh3
    dry = air / (1.0 + hum)
    return np.array([y_nh3, Y_O2_AIR * dry, (1 - Y_O2_AIR) * dry, 0.0, air - dry])


def burn(n0, S):
    """Outlet moles for complete NH3 conversion with NO selectivity S."""
    a = n0[0]
    return n0 + np.array([-a, -a * (1.25 * S + 0.75 * (1 - S)), 0.5 * a * (1 - S), a * S, 1.5 * a])


def gauze_temperature(y_nh3, T_mix, S=0.96, hum=0.0):
    n0 = feed(y_nh3, hum)
    n1 = burn(n0, S)
    H_in = n0 @ h(T_mix)
    return brentq(lambda T: n1 @ h(T) - H_in, T_mix, 2500.0), n1


if __name__ == "__main__":
    ys = np.linspace(0.06, 0.14, 41)
    fig, ax = plt.subplots(1, 2, figsize=(12, 4.4))
    for tmix, c in [(100, "tab:blue"), (200, "tab:orange"), (300, "tab:red")]:
        T = [gauze_temperature(y, tmix + 273.15)[0] - 273.15 for y in ys]
        ax[0].plot(100 * ys, T, lw=2, color=c, label=f"energy balance, mixture at {tmix} °C")
        ax[0].plot(100 * ys, tmix + 68 * 100 * ys, "--", color=c, lw=1, label="rule of thumb t$_{mix}$ + 68·(% NH$_3$)" if tmix == 100 else None)
    ax[0].axhspan(840, 950, color="0.85", zorder=0, label="usual gauze temperatures")
    ax[0].plot([10.0, 10.7], [890, 900], "ks", ms=6, label="plant data: 10 % at 220 °C, 10.7 % at 175 °C")
    ax[0].axvline(11.0, color="k", ls=":", lw=1); ax[0].axvline(13.5, color="k", ls=":", lw=1)
    ax[0].text(10.9, 530, "explosion-safety limit,\nhigh pressure", fontsize=7, ha="right"); ax[0].text(13.4, 530, "low\npressure", fontsize=7, ha="right")
    ax[0].set(xlabel="NH$_3$ in NH$_3$–air mixture (mol-%)", ylabel="gauze temperature (°C)", title="Adiabatic gauze temperature, S = 96 %")
    ax[0].legend(fontsize=7, loc="upper left"); ax[0].grid(alpha=0.3)
    # oxygen balance: O2 left after the burner vs O2 still needed for 4 NO + 3 O2 + 2 H2O -> 4 HNO3
    for S in [0.93, 0.96, 0.98]:
        left, need = [], []
        for y in ys:
            n1 = burn(feed(y), S)
            left.append(n1[1] / y); need.append(0.75 * n1[3] / y)
        ax[1].plot(100 * ys, np.array(left) - np.array(need), lw=2, label=f"S = {100*S:.0f} %")
    ax[1].axhline(0, color="k", lw=1)
    ax[1].set(xlabel="NH$_3$ in NH$_3$–air mixture (mol-%)", ylabel="O$_2$ surplus after burner (mol per mol NH$_3$)",
              title="Oxygen left for NO oxidation and absorption")
    ax[1].legend(); ax[1].grid(alpha=0.3)
    fig.tight_layout(); fig.savefig("hno3_burner.png", dpi=150)
    for tmix in [100, 200, 300]:
        T, n1 = gauze_temperature(0.10, tmix + 273.15)
        y = n1 / n1.sum()
        print(f"10 % NH3, mixture {tmix} °C: gauze {T-273.15:.0f} °C (rule {tmix+680:.0f} °C); outlet NO {100*y[3]:.2f} %, O2 {100*y[1]:.2f} %, H2O {100*y[4]:.2f} %")
    for y, tmix, S, T_plant in [(0.10, 220, 0.965, 890), (0.107, 175, 0.96, 900)]:
        T, n1 = gauze_temperature(y, tmix + 273.15, S)
        print(f"plant check: {100*y:.1f} % NH3 at {tmix} °C, S = {S}: model {T-273.15:.0f} °C, plant ca. {T_plant} °C")
    for S in [0.93, 0.96, 0.98]:
        n1 = burn(feed(0.10), S)
        print(f"S = {S:.2f}: O2 surplus after burner = {(n1[1] - 0.75*n1[3])/0.10:+.3f} mol O2 per mol NH3")
    plt.show()