"""
Methanol synthesis — chemical equilibrium of the CO2-hydrogenation and shift reactions
bdrworks.org · Handbook · Methanol, supporting script

Equilibrium CH3OH content and carbon conversion of a synthesis gas as a function of
temperature and pressure (ideal gas, Graaf equilibrium constants).

Run:  python meoh_equilibrium.py
"""
import numpy as np
import matplotlib.pyplot as plt

# ---------------------------------------------------------------- data
SPECIES = ["CO", "CO2", "H2", "CH3OH", "H2O", "CH4", "N2"]
#                 CO   CO2   H2  MeOH  H2O  CH4  N2
NU = np.array([[0.0, -1.0, -3.0, 1.0, 1.0, 0.0, 0.0],     # 1: CO2 + 3H2 <=> CH3OH + H2O
               [1.0, -1.0, -1.0, 0.0, 1.0, 0.0, 0.0]])    # 2: CO2 + H2  <=> CO + H2O (RWGS)
RG = 8.314                                                 # J/(mol K)

# ideal-gas cp = A + B T + C T^2 + D T^3, J/(mol K)  (Reid, Prausnitz & Poling)
CP = np.array([[30.87, -1.285e-2, 2.789e-5, -1.272e-8],   # CO
               [19.80, 7.344e-2, -5.602e-5, 1.715e-8],    # CO2
               [27.14, 9.274e-3, -1.381e-5, 7.645e-9],    # H2
               [21.15, 7.092e-2, 2.587e-5, -2.852e-8],    # CH3OH
               [32.24, 1.924e-3, 1.055e-5, -3.596e-9],    # H2O
               [19.25, 5.213e-2, 1.197e-5, -1.132e-8],    # CH4
               [31.15, -1.357e-2, 2.680e-5, -1.168e-8]])  # N2
HF = np.array([-110.53, -393.51, 0.0, -200.94, -241.83, -74.87, 0.0]) * 1e3   # J/mol, 298.15 K


def cp(T):
    """Heat capacities of all species, J/(mol K)."""
    return CP @ np.array([1.0, T, T**2, T**3])


def h(T):
    """Molar enthalpies H_f(298) + integral of cp, J/mol."""
    T0 = 298.15
    dT = np.array([T - T0, (T**2 - T0**2) / 2, (T**3 - T0**3) / 3, (T**4 - T0**4) / 4])
    return HF + CP @ dT


def dH_rxn(T):
    """Heats of reaction 1 and 2, J/mol."""
    return NU @ h(T)


def K_eq(T):
    """Equilibrium constants (Graaf et al.), partial pressures in bar.
    K1 for CO2 + 3H2 = CH3OH + H2O (bar^-2);  K3 for CO + H2O = CO2 + H2 (-)."""
    K1 = 10.0 ** (3066.0 / T - 10.592)
    K3 = 10.0 ** (2073.0 / T - 2.029)
    return K1, K3


def rates(T, p):
    """Vanden Bussche & Froment (1996) rates, mol/(kg_cat s); p = partial pressures in bar."""
    pCO, pCO2, pH2, pM, pW = [max(x, 1e-10) for x in p[:5]]
    k = lambda A, B: A * np.exp(B / (RG * T))
    k1, k2, k3, k4, k5 = k(1.07, 36696), k(3453.38, 0), k(0.499, 17197), k(6.62e-11, 124119), k(1.22e10, -94765)
    K1, K3 = K_eq(T)
    den = 1.0 + k2 * pW / pH2 + k3 * np.sqrt(pH2) + k4 * pW
    r1 = k1 * pCO2 * pH2 * (1.0 - pW * pM / (K1 * pH2**3 * pCO2)) / den**3
    r2 = k5 * pCO2 * (1.0 - K3 * pW * pCO / (pCO2 * pH2)) / den
    return np.array([r1, r2])


def equilibrium(T, P, F):
    """Chemical equilibrium of reactions 1 and 2 from flows F (any unit); returns outlet flows."""
    from scipy.optimize import fsolve
    K1, K3 = K_eq(T)
    def res(x):
        n = F + NU.T @ x
        n = np.maximum(n, 1e-12)
        p = P * n / n.sum()
        return [np.log(p[3] * p[4] / (p[1] * p[2]**3)) - np.log(K1),
                np.log(p[0] * p[4] / (p[1] * p[2])) + np.log(K3)]
    x = fsolve(res, [0.1 * F.sum() * 0.05, 0.0], xtol=1e-12)
    return F + NU.T @ x


if __name__ == "__main__":
    Y0 = np.array([0.10, 0.05, 0.70, 0.005, 0.001, 0.12, 0.024]); Y0 = Y0 / Y0.sum()
    Ts = np.linspace(470, 620, 61)
    fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
    for P in [50, 75, 100]:
        n = [equilibrium(t, P, Y0) for t in Ts]
        ax[0].plot(Ts, [100 * x[3] / x.sum() for x in n], lw=2, label=f"{P} bar")
        ax[1].plot(Ts, [100 * (1 - (x[0] + x[1]) / (Y0[0] + Y0[1])) for x in n], lw=2, label=f"{P} bar")
        e = equilibrium(523.0, P, Y0)
        print(f"{P} bar, 523 K: CH3OH {100*e[3]/e.sum():.2f} %, carbon conversion {100*(1-(e[0]+e[1])/(Y0[0]+Y0[1])):.1f} %")
    ax[0].set(xlabel="temperature (K)", ylabel="CH$_3$OH at equilibrium (mol-%)", title="Equilibrium methanol")
    ax[1].set(xlabel="temperature (K)", ylabel="CO + CO$_2$ converted (%)", title="Equilibrium carbon conversion")
    for a in ax:
        a.grid(alpha=0.3); a.legend()
    fig.tight_layout(); fig.savefig("meoh_equilibrium.png", dpi=150); plt.show()