"""worked_examples.py
=====================================================================
Ground-truth numerics for notes/examples.tex ("Worked examples").

Every number quoted in the LaTeX chapter is printed by this script.
Nothing is estimated by hand.

The verified routines of the lecture package are reused throughout:

    eosstyle.py    R, Rb, SPECIES, cubic_params, Zroots, P_of_V
    figB_common.py ln_phi, Psat_solve, maxwell_equal_area_check

Run with:   python3 /home/claude/eos_lecture/code/worked_examples.py
=====================================================================
"""
import sys

sys.path.insert(0, "/home/claude/eos_lecture/code")

import numpy as np
from scipy.optimize import brentq
from scipy.integrate import quad

from eosstyle import R, Rb, SPECIES, cubic_params, Zroots, P_of_V
from figB_common import ln_phi, Psat_solve, maxwell_equal_area_check

S2 = np.sqrt(2.0)


def head(n, title):
    """Labelled block header, one per worked example."""
    print("\n" + "=" * 78)
    print(f"EXAMPLE {n}: {title}")
    print("=" * 78)


def kappa_pr(w):
    return 0.37464 + 1.54226 * w - 0.26992 * w ** 2


def m_srk(w):
    return 0.480 + 1.574 * w - 0.176 * w ** 2


def psat_fine(model, sp, T, Plo, Phi_, n=4000):
    """Isofugacity saturation pressure found on a fine linear pressure grid.

    Psat_solve() in figB_common uses a coarse geometric grid, which can miss
    the very narrow three-root window of the van der Waals equation close to
    the critical point.  This routine is used only where that happens; it
    solves exactly the same equation, ln(phi_V) = ln(phi_L).
    """
    def f(P):
        r = Zroots(model, sp, T, P)
        if len(r) < 2:
            return np.nan
        return ln_phi(model, sp, T, P, r[-1]) - ln_phi(model, sp, T, P, r[0])

    Ps = np.linspace(Plo, Phi_, n)
    v = np.array([f(P) for P in Ps])
    ok = np.isfinite(v)
    Pf, vf = Ps[ok], v[ok]
    sc = np.where(np.diff(np.sign(vf)) != 0)[0]
    if len(sc) == 0:
        return np.nan
    i = sc[0]
    return brentq(f, Pf[i], Pf[i + 1], xtol=1e-12)


# =====================================================================
# EXAMPLE 1  Compressibility factor and its engineering consequence
# =====================================================================
head(1, "Compressibility factor from an EOS and its engineering consequence")

sp, T, P = "methane", 300.0, 100.0            # K, bar
d = SPECIES[sp]
M_CH4 = 16.043                                 # g/mol, molar mass
V_VESSEL = 50.0                                # m3
NDOT = 100.0                                   # kmol/h compressor feed

print(f"Species        : {sp}  Tc={d['Tc']} K  Pc={d['Pc']} bar  w={d['w']}, "
      f"M = {M_CH4} g/mol")
print(f"State          : T = {T:.1f} K, P = {P:.1f} bar")
print(f"Reduced        : Tr = {T/d['Tc']:.4f}, Pr = {P/d['Pc']:.4f}")

# --- ideal gas -------------------------------------------------------
V_ig = Rb * T / P                              # cm3/mol
print(f"\nIdeal gas      : Z = 1, V = RT/P = {Rb:.2f}*{T:.0f}/{P:.0f}"
      f" = {V_ig:.4f} cm3/mol")

# --- Peng-Robinson ---------------------------------------------------
a_pr, b_pr, _, _ = cubic_params("pr", d["Tc"], d["Pc"], d["w"], T)
kap1 = kappa_pr(d["w"])
al1 = (1 + kap1 * (1 - np.sqrt(T / d["Tc"]))) ** 2
A_pr = a_pr * P / (Rb * T) ** 2
B_pr = b_pr * P / (Rb * T)
roots = Zroots("pr", sp, T, P)
Z_pr = roots[-1]
V_pr = Z_pr * Rb * T / P
print(f"PR parameters  : kappa = {kap1:.6f}, alpha = {al1:.6f}")
print(f"                 a(T) = {a_pr:.4f} cm6 bar/mol2, b = {b_pr:.4f} cm3/mol")
print(f"                 A = {A_pr:.6f}, B = {B_pr:.6f}")
print(f"                 cubic coefficients: Z^3 + ({-(1-B_pr):.6f})Z^2 + "
      f"({A_pr-2*B_pr-3*B_pr**2:.6f})Z + ({-(A_pr*B_pr-B_pr**2-B_pr**3):.6f})")
print(f"                 real roots of the cubic: {np.round(roots, 6)}")
print(f"PR             : Z = {Z_pr:.6f}, V = ZRT/P = {V_pr:.4f} cm3/mol")
print(f"                 difference V_ig - V_PR = {V_ig - V_pr:.4f} cm3/mol")

# --- consequence 1: inventory in a fixed vessel ----------------------
n_ig = V_VESSEL * 1e6 / V_ig                   # mol
n_pr = V_VESSEL * 1e6 / V_pr
m_ig = n_ig * M_CH4 / 1000.0                   # kg
m_pr = n_pr * M_CH4 / 1000.0
print(f"\nVessel {V_VESSEL:.0f} m3   : n_ig = V/V_ig = {n_ig/1000:.4f} kmol, "
      f"n_PR = {n_pr/1000:.4f} kmol")
print(f"                 m_ig = {m_ig:.4f} kg, m_PR = {m_pr:.4f} kg, "
      f"shortfall = {m_pr - m_ig:.4f} kg")
print(f"                 m_ig/m_PR - 1 = Z - 1 = {(Z_pr-1)*100:.4f} %, "
      f"i.e. the ideal gas law understates the inventory by "
      f"{(1-Z_pr)*100:.4f} %")
print(f"                 the ideal-gas error in Z itself is 1 - Z = "
      f"{1-Z_pr:.6f}")

# --- consequence 2: compressor volumetric flow -----------------------
Q_ig = NDOT * 1000.0 * V_ig / 1e6              # m3/h
Q_pr = NDOT * 1000.0 * V_pr / 1e6
print(f"Flow {NDOT:.0f} kmol/h: Q_ig = {Q_ig:.4f} m3/h, Q_PR = {Q_pr:.4f} m3/h")
print(f"                 Q_ig/Q_PR - 1 = 1/Z - 1 = {(1/Z_pr - 1)*100:.4f} %, "
      f"so a machine sized on the ideal gas law is oversized by that margin")

# --- checks ----------------------------------------------------------
print(f"\n[check] back-substitution P(V_PR) = {P_of_V('pr', sp, T, V_pr):.6f} bar "
      f"against the specified {P:.1f} bar")
B2_pr = b_pr - a_pr / (Rb * T)                 # cm3/mol
Z_virial = 1 + B2_pr * P / (Rb * T)
print(f"[check] PR second virial coefficient B2 = b - a/RT = "
      f"{B2_pr:.4f} cm3/mol; two-term virial Z = 1 + B2 P/RT = "
      f"{Z_virial:.6f}")
print(f"[check] the virial truncation has the right sign and magnitude but is "
      f"{abs(Z_virial-Z_pr):.4f} low at Pr = {P/d['Pc']:.4f}, as expected of a "
      f"density expansion truncated after one term")


# =====================================================================
# EXAMPLE 2  van der Waals parameters and the cubic solved by hand
# =====================================================================
head(2, "van der Waals parameters from critical constants; the cubic in V")

sp2, T2_, P2_ = "CO2", 300.0, 60.0
d2 = SPECIES[sp2]
VC_CO2 = 94.0                                  # cm3/mol, literature critical volume

a_vdw = 27 * Rb ** 2 * d2["Tc"] ** 2 / (64 * d2["Pc"])
b_vdw = Rb * d2["Tc"] / (8 * d2["Pc"])
print(f"Species        : {sp2}  Tc={d2['Tc']} K  Pc={d2['Pc']} bar")
print(f"State          : T = {T2_:.1f} K, P = {P2_:.1f} bar, "
      f"Tr = {T2_/d2['Tc']:.4f}, Pr = {P2_/d2['Pc']:.4f}")
print(f"\na = 27 R^2 Tc^2 / (64 Pc) = 27*{Rb:.2f}^2*{d2['Tc']}^2/(64*{d2['Pc']})"
      f" = {a_vdw:.4f} cm6 bar/mol2")
print(f"b = R Tc / (8 Pc)         = {Rb:.2f}*{d2['Tc']}/(8*{d2['Pc']})"
      f" = {b_vdw:.4f} cm3/mol")

# cubic in V:  V^3 - (b + RT/P) V^2 + (a/P) V - ab/P = 0
c2_ = -(b_vdw + Rb * T2_ / P2_)
c1_ = a_vdw / P2_
c0_ = -a_vdw * b_vdw / P2_
print(f"\nRT/P = {Rb*T2_/P2_:.4f} cm3/mol")
print(f"Cubic in V     : V^3 + ({c2_:.4f}) V^2 + ({c1_:.4f}) V "
      f"+ ({c0_:.4f}) = 0")
Vr_all = np.roots([1.0, c2_, c1_, c0_])
Vreal = np.sort(np.real(Vr_all[np.abs(np.imag(Vr_all)) < 1e-9]))
Vcplx = Vr_all[np.abs(np.imag(Vr_all)) >= 1e-9]
print(f"Real roots     : {np.round(Vreal, 4)} cm3/mol   "
      f"({len(Vreal)} real, {len(Vcplx)} complex conjugate)")
if len(Vcplx):
    print(f"Complex pair   : {Vcplx[0].real:.4f} +/- {abs(Vcplx[0].imag):.4f} i "
          f"cm3/mol (no physical meaning)")

Zs2 = Vreal * P2_ / (Rb * T2_)
for Vv, Zz in zip(Vreal, Zs2):
    lp = ln_phi("vdw", sp2, T2_, P2_, Zz)
    hh = 1e-4 * Vv
    dpdv = (P_of_V("vdw", sp2, T2_, Vv + hh)
            - P_of_V("vdw", sp2, T2_, Vv - hh)) / (2 * hh)
    print(f"   V = {Vv:10.4f} cm3/mol  Z = PV/RT = {Zz:.6f}  "
          f"ln phi = {lp:.6f}  phi = {np.exp(lp):.6f}  "
          f"(dP/dV)_T = {dpdv:.6e} bar mol/cm3")
print(f"Selection      : (i) discard roots with V < b = {b_vdw:.4f} cm3/mol; "
      f"(ii) discard any root with (dP/dV)_T > 0 (mechanically unstable); "
      f"(iii) of the survivors keep the one of lowest G = RT ln phi.")
print(f"Selected root  : V = {Vreal[-1]:.4f} cm3/mol, Z = {Zs2[-1]:.6f}, "
      f"phi = {np.exp(ln_phi('vdw', sp2, T2_, P2_, Zs2[-1])):.6f} "
      f"(single real root, so the selection is forced)")

# saturation pressure of the vdW model at this temperature
Psat_vdw = psat_fine("vdw", sp2, T2_, 60.0, 73.7)
print(f"\nvdW Psat({T2_:.0f} K) = {Psat_vdw:.4f} bar (fine-grid isofugacity), "
      f"so P = {P2_:.1f} bar lies below it and the fluid is a vapour")
print(f"For comparison : PR Psat({T2_:.0f} K) = "
      f"{Psat_solve('pr', sp2, T2_):.4f} bar")
Z2_pr = Zroots("pr", sp2, T2_, P2_)[-1]
print(f"For comparison : PR at the same state gives Z = {Z2_pr:.6f}, "
      f"V = {Z2_pr*Rb*T2_/P2_:.4f} cm3/mol, i.e. the vdW molar volume is high "
      f"by {(Vreal[-1]/(Z2_pr*Rb*T2_/P2_) - 1)*100:.4f} %")

# the three-root case at 70 bar, same temperature
P2b = 70.0
r2b = Zroots("vdw", sp2, T2_, P2b)
V2b = r2b * Rb * T2_ / P2b
lp2b = [ln_phi("vdw", sp2, T2_, P2b, z) for z in r2b]
print(f"\nAt {P2b:.0f} bar and the same temperature the same cubic has "
      f"{len(r2b)} real roots, so the criterion has work to do:")
for nm, z, v, lp in zip(["liquid", "unstable", "vapour"], r2b, V2b, lp2b):
    print(f"   {nm:<10} Z = {z:.6f}  V = {v:9.4f} cm3/mol  "
          f"ln phi = {lp:.8f}")
print(f"   G_V - G_L = RT[ln phi_V - ln phi_L] = "
      f"{R*T2_*(lp2b[-1]-lp2b[0]):.4f} J/mol -> the "
      f"{'liquid' if lp2b[-1] > lp2b[0] else 'vapour'} root is stable, "
      f"consistent with P = {P2b:.0f} bar > Psat = {Psat_vdw:.4f} bar")

# vdW Zc versus the measured value
Zc_vdw = 3.0 / 8.0
Zc_meas = d2["Pc"] * VC_CO2 / (Rb * d2["Tc"])
print(f"\nvdW Zc = 3/8 = {Zc_vdw:.4f}; measured CO2 Zc = Pc Vc/(R Tc) = "
      f"{d2['Pc']}*{VC_CO2}/({Rb:.2f}*{d2['Tc']}) = {Zc_meas:.4f} "
      f"(Vc = {VC_CO2:.1f} cm3/mol, literature)")
print(f"vdW overpredicts Zc by {(Zc_vdw/Zc_meas - 1)*100:.4f} %; "
      f"Vc(vdW) = 3b = {3*b_vdw:.4f} cm3/mol against {VC_CO2:.1f} cm3/mol, "
      f"the same {(3*b_vdw/VC_CO2 - 1)*100:.4f} % error")
print(f"[check] the vdW critical point is reproduced exactly by its own "
      f"constants: Z at (Tc, Pc) = "
      f"{Zroots('vdw', sp2, d2['Tc'], d2['Pc'])[-1]:.6f} vs 3/8 = 0.375000")


# =====================================================================
# EXAMPLE 3  Three real roots and root selection
# =====================================================================
head(3, "Three real roots and root selection (the commonest practical error)")

sp3, T3, P3 = "n-butane", 350.0, 10.0
d3 = SPECIES[sp3]
a3, b3, _, _ = cubic_params("pr", d3["Tc"], d3["Pc"], d3["w"], T3)
A3 = a3 * P3 / (Rb * T3) ** 2
B3 = b3 * P3 / (Rb * T3)
k3 = kappa_pr(d3["w"])
al3 = (1 + k3 * (1 - np.sqrt(T3 / d3["Tc"]))) ** 2
print(f"Species        : {sp3}  Tc={d3['Tc']} K  Pc={d3['Pc']} bar  w={d3['w']}")
print(f"State          : T = {T3:.1f} K, P = {P3:.1f} bar, "
      f"Tr = {T3/d3['Tc']:.4f}, Pr = {P3/d3['Pc']:.4f}")
print(f"PR             : kappa = {k3:.6f}, alpha = {al3:.6f}")
print(f"                 a(T) = {a3:.4f} cm6 bar/mol2, b = {b3:.4f} cm3/mol")
print(f"                 A = {A3:.6f}, B = {B3:.6f}")
print(f"Cubic          : Z^3 - (1-B)Z^2 + (A-2B-3B^2)Z - (AB-B^2-B^3) = 0")
print(f"                 Z^3 + ({-(1-B3):.6f})Z^2 + ({A3-2*B3-3*B3**2:.6f})Z "
      f"+ ({-(A3*B3-B3**2-B3**3):.6f}) = 0")

r3 = Zroots("pr", sp3, T3, P3)
print(f"Real roots     : {np.round(r3, 6)}  ({len(r3)} roots)")
V3 = r3 * Rb * T3 / P3
names = ["liquid", "unstable (middle)", "vapour"]
lp3 = [ln_phi("pr", sp3, T3, P3, z) for z in r3]
dPdV3 = []
for v in V3:
    h = 1e-4 * v
    dPdV3.append((P_of_V("pr", sp3, T3, v + h)
                  - P_of_V("pr", sp3, T3, v - h)) / (2 * h))
print(f"\n{'branch':<20}{'Z':>11}{'V/cm3 mol-1':>14}{'ln phi':>11}{'phi':>11}"
      f"{'(dP/dV)_T':>15}")
for nm, z, v, lp, dp in zip(names, r3, V3, lp3, dPdV3):
    print(f"{nm:<20}{z:>11.6f}{v:>14.4f}{lp:>11.6f}{np.exp(lp):>11.6f}"
          f"{dp:>15.6e}")

dG3 = R * T3 * (lp3[-1] - lp3[0])
print(f"\nG_vap - G_liq = RT[ln phi_V - ln phi_L] = "
      f"{R:.3f}*{T3:.0f}*({lp3[-1]:.6f} - ({lp3[0]:.6f})) = {dG3:.4f} J/mol")
print(f"Stable phase   : {'liquid' if dG3 > 0 else 'vapour'} "
      f"(the lower molar Gibbs energy)")
Psat3 = Psat_solve("pr", sp3, T3)
print(f"[check] PR Psat({T3:.0f} K) = {Psat3:.4f} bar; the state pressure "
      f"{P3:.1f} bar exceeds it, so the fluid must be a subcooled liquid, "
      f"which is what the Gibbs test returned")
print(f"[check] taking the vapour root by mistake would give "
      f"V_V/V_L = {V3[-1]/V3[0]:.4f}, a factor {V3[-1]/V3[0]:.2f} error in "
      f"molar volume and hence in density, mass inventory and holdup")
print(f"[check] mass in a 50 m3 drum: liquid root "
      f"{50e6/V3[0]*58.122/1e6:.4f} t, vapour root "
      f"{50e6/V3[-1]*58.122/1e6:.4f} t (M = 58.122 g/mol)")


# =====================================================================
# EXAMPLE 4  RK, SRK and PR compared on one state point
# =====================================================================
head(4, "RK, SRK and PR compared on one state point")

sp4, T4, P4 = "propane", 350.0, 15.0
d4 = SPECIES[sp4]
print(f"Species        : {sp4}  Tc={d4['Tc']} K  Pc={d4['Pc']} bar  w={d4['w']}")
print(f"State          : T = {T4:.1f} K, P = {P4:.1f} bar, "
      f"Tr = {T4/d4['Tc']:.4f}, Pr = {P4/d4['Pc']:.4f}")
print(f"PR Psat({T4:.0f} K) = {Psat_solve('pr', sp4, T4):.4f} bar, so the state "
      f"is a superheated vapour and one real root is expected")
print(f"\n{'model':<6}{'a(T)/cm6 bar mol-2':>22}{'b/cm3 mol-1':>14}"
      f"{'A':>11}{'B':>11}{'Z':>11}{'V/cm3 mol-1':>14}{'phi':>10}{'roots':>7}")
ex4 = {}
for m in ("rk", "srk", "pr"):
    am, bm, _, _ = cubic_params(m, d4["Tc"], d4["Pc"], d4["w"], T4)
    Am = am * P4 / (Rb * T4) ** 2
    Bm = bm * P4 / (Rb * T4)
    rr = Zroots(m, sp4, T4, P4)
    Zm = rr[-1]
    Vm = Zm * Rb * T4 / P4
    phim = np.exp(ln_phi(m, sp4, T4, P4, Zm))
    ex4[m] = (am, bm, Am, Bm, Zm, Vm, phim, len(rr))
    print(f"{m.upper():<6}{am:>22.1f}{bm:>14.4f}{Am:>11.6f}{Bm:>11.6f}"
          f"{Zm:>11.6f}{Vm:>14.4f}{phim:>10.6f}{len(rr):>7d}")

print(f"\nZ(ideal gas) = 1, V(ideal gas) = {Rb*T4/P4:.4f} cm3/mol")
for m in ("rk", "srk", "pr"):
    print(f"  {m.upper():<4} Z - 1 = {ex4[m][4]-1:+.6f} "
          f"({(ex4[m][4]-1)*100:+.4f} %), "
          f"V_ig - V = {Rb*T4/P4 - ex4[m][5]:.4f} cm3/mol")
print(f"alpha values   : RK {(T4/d4['Tc'])**-0.5:.6f}, "
      f"SRK {(1+m_srk(d4['w'])*(1-np.sqrt(T4/d4['Tc'])))**2:.6f}, "
      f"PR {(1+kappa_pr(d4['w'])*(1-np.sqrt(T4/d4['Tc'])))**2:.6f}")
print(f"Spread in Z (RK to PR): {abs(ex4['rk'][4]-ex4['pr'][4]):.6f}, "
      f"{abs(ex4['rk'][4]-ex4['pr'][4])/ex4['pr'][4]*100:.4f} % of the PR value")
print(f"Spread in phi (RK to PR): {abs(ex4['rk'][6]-ex4['pr'][6]):.6f}, "
      f"{abs(ex4['rk'][6]-ex4['pr'][6])/ex4['pr'][6]*100:.4f} %")
print(f"Ratio b(PR)/b(SRK) = {ex4['pr'][1]/ex4['srk'][1]:.6f} "
      f"= 0.07780/0.08664 = {0.07780/0.08664:.6f}")

print()
for m in ("rk", "srk", "pr"):
    def integ(Pv, m=m):
        if Pv < 1e-12:
            am, bm, _, _ = cubic_params(m, d4["Tc"], d4["Pc"], d4["w"], T4)
            return (bm - am / (Rb * T4)) / (Rb * T4)
        return (Zroots(m, sp4, T4, Pv)[-1] - 1.0) / Pv
    val, _ = quad(integ, 0.0, P4, limit=200)
    ana = ln_phi(m, sp4, T4, P4, ex4[m][4])
    print(f"[check] {m.upper():<4} ln phi analytic = {ana: .8f}, "
          f"by integration of (Z-1)dP/P = {val: .8f}, diff = "
          f"{abs(val-ana):.2e}")


# =====================================================================
# EXAMPLE 5  The alpha function and its sensitivity
# =====================================================================
head(5, "The alpha function and its sensitivity to the acentric factor")

sp5 = "propane"
d5 = SPECIES[sp5]
TRS = [0.6, 0.8, 1.0, 1.3]
ac_rk = 0.42748 * Rb ** 2 * d5["Tc"] ** 2 / d5["Pc"]
ac_pr = 0.45724 * Rb ** 2 * d5["Tc"] ** 2 / d5["Pc"]
ms5 = m_srk(d5["w"])
kp5 = kappa_pr(d5["w"])
print(f"Species        : {sp5}  Tc={d5['Tc']} K  Pc={d5['Pc']} bar  w={d5['w']}")
print(f"ac(RK and SRK) = 0.42748 R^2 Tc^2/Pc = {ac_rk:.1f} cm6 bar/mol2")
print(f"ac(PR)         = 0.45724 R^2 Tc^2/Pc = {ac_pr:.1f} cm6 bar/mol2")
print(f"m(SRK)   = 0.480 + 1.574 w - 0.176 w^2   = {ms5:.6f}")
print(f"kappa(PR) = 0.37464 + 1.54226 w - 0.26992 w^2 = {kp5:.6f}")
print(f"alpha_RK = Tr^-1/2, alpha_SRK = [1+m(1-sqrt(Tr))]^2, "
      f"alpha_PR = [1+kappa(1-sqrt(Tr))]^2")

print(f"\n{'Tr':>6}{'T/K':>10}"
      f"{'al_RK':>10}{'al_SRK':>10}{'al_PR':>10}"
      f"{'a_RK/1e6':>12}{'a_SRK/1e6':>12}{'a_PR/1e6':>12}")
for tr in TRS:
    Tx = tr * d5["Tc"]
    al_rk = tr ** -0.5
    al_srk = (1 + ms5 * (1 - np.sqrt(tr))) ** 2
    al_pr = (1 + kp5 * (1 - np.sqrt(tr))) ** 2
    a_rk = cubic_params("rk", d5["Tc"], d5["Pc"], d5["w"], Tx)[0]
    a_srk = cubic_params("srk", d5["Tc"], d5["Pc"], d5["w"], Tx)[0]
    a_pr = cubic_params("pr", d5["Tc"], d5["Pc"], d5["w"], Tx)[0]
    print(f"{tr:>6.1f}{Tx:>10.2f}{al_rk:>10.4f}{al_srk:>10.4f}{al_pr:>10.4f}"
          f"{a_rk/1e6:>12.4f}{a_srk/1e6:>12.4f}{a_pr/1e6:>12.4f}")
    assert abs(a_rk - ac_rk * al_rk) < 1e-6 * ac_rk
    assert abs(a_srk - ac_rk * al_srk) < 1e-6 * ac_rk
    assert abs(a_pr - ac_pr * al_pr) < 1e-6 * ac_pr
print("[check] a(T) = ac alpha reproduced to better than 1e-6 relative for "
      "all three models at all four reduced temperatures")
print(f"[check] the RK and SRK alpha functions cross at Tr = 1 by "
      f"construction: alpha_RK(1) = {1.0:.4f}, alpha_SRK(1) = "
      f"{(1+ms5*(1-1))**2:.4f}, alpha_PR(1) = {(1+kp5*(1-1))**2:.4f}")

# --- sensitivity of Psat to a 5 % error in omega ---------------------
w0 = d5["w"]
w1 = 1.05 * w0
SPECIES["propane_w5"] = dict(Tc=d5["Tc"], Pc=d5["Pc"], w=w1)
print(f"\nAcentric factor: w = {w0:.4f}; a +5 % error gives w' = {w1:.4f}")
print(f"kappa(w) = {kp5:.6f} -> kappa(w') = {kappa_pr(w1):.6f} "
      f"({kappa_pr(w1)/kp5*100 - 100:+.4f} %)")
print(f"\n{'Tr':>6}{'T/K':>10}{'Psat(w)/bar':>14}{'Psat(w+5%)/bar':>17}"
      f"{'change %':>11}{'amplification':>15}")
sens = {}
for tr in [0.6, 0.7, 0.8, 0.9, 0.95]:
    Tx = tr * d5["Tc"]
    p0 = Psat_solve("pr", sp5, Tx)
    p1 = Psat_solve("pr", "propane_w5", Tx)
    ch = (p1 - p0) / p0 * 100
    sens[tr] = (Tx, p0, p1, ch)
    print(f"{tr:>6.2f}{Tx:>10.2f}{p0:>14.6f}{p1:>17.6f}{ch:>11.4f}"
          f"{abs(ch)/5.0:>15.5f}")
print(f"Range of the saturation pressure over the tabulated interval: "
      f"Psat(Tr=0.95)/Psat(Tr=0.60) = {sens[0.95][1]/sens[0.6][1]:.4f}")
print(f"[check] the largest amplification |dPsat/Psat|/|dw/w| is "
      f"{max(abs(v[3]) for v in sens.values())/5.0:.4f} at the lowest reduced "
      f"temperature tabulated")
print(f"[check] alpha(Tr = 1) = 1 whatever w, so a(Tc) and hence Psat(Tc) = Pc "
      f"are independent of w and the sensitivity must decay to zero as "
      f"Tr -> 1; the tabulated amplification falls monotonically from "
      f"{abs(sens[0.6][3])/5.0:.4f} at Tr = 0.60 to "
      f"{abs(sens[0.95][3])/5.0:.4f} at Tr = 0.95")
print(f"[check] a(Tc) is the same for both acentric factors: "
      f"{cubic_params('pr', d5['Tc'], d5['Pc'], w0, d5['Tc'])[0]:.4f} and "
      f"{cubic_params('pr', d5['Tc'], d5['Pc'], w1, d5['Tc'])[0]:.4f} "
      f"cm6 bar/mol2")


# =====================================================================
# EXAMPLE 6  Fugacity coefficient from the PR EOS, verified by quadrature
# =====================================================================
head(6, "Fugacity coefficient from the PR EOS, verified by quadrature")

sp6, T6, P6 = "CO2", 320.0, 70.0
d6 = SPECIES[sp6]
a6, b6, _, _ = cubic_params("pr", d6["Tc"], d6["Pc"], d6["w"], T6)
k6 = kappa_pr(d6["w"])
al6 = (1 + k6 * (1 - np.sqrt(T6 / d6["Tc"]))) ** 2
A6 = a6 * P6 / (Rb * T6) ** 2
B6 = b6 * P6 / (Rb * T6)
r6 = Zroots("pr", sp6, T6, P6)
Z6 = r6[-1]
V6 = Z6 * Rb * T6 / P6
print(f"Species        : {sp6}  Tc={d6['Tc']} K  Pc={d6['Pc']} bar  w={d6['w']}")
print(f"State          : T = {T6:.1f} K, P = {P6:.1f} bar, "
      f"Tr = {T6/d6['Tc']:.4f}, Pr = {P6/d6['Pc']:.4f}; supercritical, "
      f"{len(r6)} real root")
print(f"kappa = {k6:.6f}, alpha = {al6:.6f}")
print(f"a(T) = {a6:.4f} cm6 bar/mol2, b = {b6:.4f} cm3/mol")
print(f"A = a P/(RT)^2 = {A6:.6f}, B = b P/(RT) = {B6:.6f}")
print(f"Z = {Z6:.6f}, V = {V6:.4f} cm3/mol")

t1 = Z6 - 1.0
t2 = -np.log(Z6 - B6)
targ = (Z6 + (1 + S2) * B6) / (Z6 + (1 - S2) * B6)
t3 = -A6 / (2 * S2 * B6) * np.log(targ)
lnphi6 = t1 + t2 + t3
print(f"\nA/(2 sqrt2 B)                 = {A6/(2*S2*B6):.8f}")
print(f"log argument                  = {targ:.8f}, ln = {np.log(targ):.8f}")
print(f"term 1  Z - 1                 = {t1:.8f}")
print(f"term 2  -ln(Z - B)            = {t2:.8f}")
print(f"term 3  -A/(2 sqrt2 B) ln(..) = {t3:.8f}")
print(f"ln phi (analytic)             = {lnphi6:.8f}")
print(f"phi    (analytic)             = {np.exp(lnphi6):.8f}")
print(f"fugacity f = phi P              = {np.exp(lnphi6)*P6:.4f} bar")
print(f"[cross-check with figB_common.ln_phi] "
      f"{ln_phi('pr', sp6, T6, P6, Z6):.8f}")


def integ6(Pv):
    """(Z - 1)/P along the 320 K PR isotherm; the P -> 0 limit is
    (b - a/RT)/RT, the second virial coefficient divided by RT."""
    if Pv < 1e-12:
        return (b6 - a6 / (Rb * T6)) / (Rb * T6)
    return (Zroots("pr", sp6, T6, Pv)[-1] - 1.0) / Pv


val6, err6 = quad(integ6, 0.0, P6, limit=400)
print(f"\nln phi = int_0^P (Z-1) dP/P   = {val6:.8f} "
      f"(quadrature error estimate {err6:.1e})")
print(f"phi    (quadrature)           = {np.exp(val6):.8f}")
print(f"absolute difference in ln phi = {abs(val6-lnphi6):.3e}")
print(f"relative difference in phi    = "
      f"{abs(np.exp(val6)-np.exp(lnphi6))/np.exp(lnphi6):.3e}")
print(f"integrand (Z-1)/P per bar     : "
      + ", ".join(f"P={p:.1f}: {integ6(p):.6f}" for p in [0.0, 17.5, 35.0, 52.5, 70.0]))
print(f"[check] B2 = b - a/RT = {b6 - a6/(Rb*T6):.4f} cm3/mol, so the P -> 0 "
      f"integrand is B2/RT = {(b6-a6/(Rb*T6))/(Rb*T6):.6f} bar^-1")


# =====================================================================
# EXAMPLE 7  Vapour pressure by the isofugacity condition
# =====================================================================
head(7, "Vapour pressure of propane at 300 K by isofugacity (PR)")

sp7, T7 = "propane", 300.0
M_C3 = 44.097                                   # g/mol
d7 = SPECIES[sp7]
Tr7 = T7 / d7["Tc"]
P0 = d7["Pc"] * np.exp(5.373 * (1 + d7["w"]) * (1 - d7["Tc"] / T7))
print(f"Species        : {sp7}  Tc={d7['Tc']} K  Pc={d7['Pc']} bar  w={d7['w']}")
print(f"T = {T7:.1f} K, Tr = {Tr7:.6f}")
print(f"Wilson initial guess P0 = Pc exp[5.373(1+w)(1 - Tc/T)] with "
      f"1 + w = {1+d7['w']:.3f}: P0 = {P0:.6f} bar")
print(f"\nSuccessive substitution  P_(k+1) = P_k phi_L/phi_V")
print(f"{'k':>3}{'P_k/bar':>12}{'Z_L':>11}{'Z_V':>11}"
      f"{'phi_L':>11}{'phi_V':>11}{'phi_L/phi_V':>14}")
Pk = P0
Pseq = []
for k in range(6):
    rr = Zroots("pr", sp7, T7, Pk)
    ZL, ZV = rr[0], rr[-1]
    pL = np.exp(ln_phi("pr", sp7, T7, Pk, ZL))
    pV = np.exp(ln_phi("pr", sp7, T7, Pk, ZV))
    ratio = pL / pV
    Pseq.append(Pk)
    print(f"{k:>3}{Pk:>12.6f}{ZL:>11.6f}{ZV:>11.6f}"
          f"{pL:>11.6f}{pV:>11.6f}{ratio:>14.8f}")
    Pk = Pk * ratio

Psat7 = Psat_solve("pr", sp7, T7)
print(f"\nConverged (Brent root of ln phi_V - ln phi_L): "
      f"Psat = {Psat7:.6f} bar")
print(f"Wilson initial estimate was high by "
      f"{(P0/Psat7 - 1)*100:.4f} % of the converged value")
errs = np.abs(np.array(Pseq) - Psat7) / Psat7
print(f"Error of each iterate relative to the converged value: "
      + ", ".join(f"{e:.3e}" for e in errs))
print(f"Successive error ratios (linear convergence): "
      + ", ".join(f"{errs[i+1]/errs[i]:.4f}" for i in range(len(errs) - 1)))
rr7 = Zroots("pr", sp7, T7, Psat7)
ZL7, ZV7 = rr7[0], rr7[-1]
VL7 = ZL7 * Rb * T7 / Psat7
VV7 = ZV7 * Rb * T7 / Psat7
print(f"Z_L = {ZL7:.6f}, Z_V = {ZV7:.6f}")
print(f"V_L = {VL7:.4f} cm3/mol, V_V = {VV7:.4f} cm3/mol, "
      f"V_V - V_L = {VV7-VL7:.4f} cm3/mol")
print(f"saturated liquid density = M/V_L = {M_C3/VL7*1000:.4f} kg/m3 "
      f"(M = {M_C3} g/mol)")
print(f"residual |ln phi_V - ln phi_L| at convergence = "
      f"{abs(ln_phi('pr', sp7, T7, Psat7, ZV7) - ln_phi('pr', sp7, T7, Psat7, ZL7)):.3e}")

# Clapeyron
dT = 0.5
Pm = Psat_solve("pr", sp7, T7 - dT)
Pp = Psat_solve("pr", sp7, T7 + dT)
dPdT = (Pp - Pm) / (2 * dT)
dHvap = T7 * (VV7 - VL7) * dPdT * 0.1     # K * cm3/mol * bar/K -> 0.1 J/mol
print(f"\nClapeyron      : Psat({T7-dT:.1f} K) = {Pm:.6f} bar, "
      f"Psat({T7+dT:.1f} K) = {Pp:.6f} bar")
print(f"                 dP/dT = ({Pp:.6f} - {Pm:.6f})/{2*dT:.1f} = "
      f"{dPdT:.6f} bar/K")
print(f"dHvap = T (V_V - V_L) dP/dT = {T7:.0f} * {VV7-VL7:.4f} * "
      f"{dPdT:.6f} * 0.1 = {dHvap:.2f} J/mol = {dHvap/1000:.4f} kJ/mol")
print(f"                 (1 cm3 bar = 0.1 J)")

print()
area7 = maxwell_equal_area_check("pr", sp7, T7, Psat7)
print(f"[check] Maxwell equal-area residual {area7:.3e} bar cm3/mol, "
      f"{abs(area7)/(Psat7*(VV7-VL7))*100:.3e} % of P(V_V - V_L)")

Pc7_Pa = d7["Pc"] * 1e5
ac7 = 0.45724 * R ** 2 * d7["Tc"] ** 2 / Pc7_Pa
kk7 = kappa_pr(d7["w"])
s7 = 1 + kk7 * (1 - np.sqrt(Tr7))
aa7 = ac7 * s7 ** 2
bb7 = 0.07780 * R * d7["Tc"] / Pc7_Pa
dadT7 = -ac7 * kk7 * s7 / np.sqrt(T7 * d7["Tc"])
BB7 = bb7 * Psat7 * 1e5 / (R * T7)


def hres_pr(Z, B, a, dadT, b, T):
    return (R * T * (Z - 1)
            + (T * dadT - a) / (2 * S2 * b)
            * np.log((Z + (1 + S2) * B) / (Z + (1 - S2) * B)))


hV = hres_pr(ZV7, BB7, aa7, dadT7, bb7, T7)
hL = hres_pr(ZL7, BB7, aa7, dadT7, bb7, T7)
print(f"[check] dHvap from the PR residual enthalpies: H_res(V) = {hV:.2f} "
      f"J/mol, H_res(L) = {hL:.2f} J/mol, difference = {hV-hL:.2f} J/mol")
print(f"[check] Clapeyron and departure routes differ by "
      f"{abs(hV-hL-dHvap)/dHvap*100:.4f} %")


# =====================================================================
# EXAMPLE 8  Binary VLE, K-values and one Rachford-Rice iteration
# =====================================================================
head(8, "Binary VLE: methane + n-butane, K-values and Rachford-Rice")

T8, P8 = 350.0, 20.0
comp = ["methane", "n-butane"]
z8 = np.array([0.5, 0.5])
KIJ = 0.0
Tc8 = np.array([SPECIES[c]["Tc"] for c in comp])
Pc8 = np.array([SPECIES[c]["Pc"] for c in comp])
w8 = np.array([SPECIES[c]["w"] for c in comp])
print(f"Components     : {comp[0]} (Tc={Tc8[0]} K, Pc={Pc8[0]} bar, w={w8[0]}), "
      f"{comp[1]} (Tc={Tc8[1]} K, Pc={Pc8[1]} bar, w={w8[1]})")
print(f"Feed           : z = [{z8[0]}, {z8[1]}], T = {T8:.1f} K, "
      f"P = {P8:.1f} bar, k_ij = {KIJ:.1f} (assumed)")

a_i = np.array([cubic_params("pr", Tc8[i], Pc8[i], w8[i], T8)[0] for i in range(2)])
b_i = np.array([cubic_params("pr", Tc8[i], Pc8[i], w8[i], T8)[1] for i in range(2)])


def aij_of(kij):
    return np.sqrt(np.outer(a_i, a_i)) * (1 - kij)


AIJ = aij_of(KIJ)
print(f"\nPure PR params : a1 = {a_i[0]:.1f}, a2 = {a_i[1]:.1f} cm6 bar/mol2")
print(f"                 b1 = {b_i[0]:.4f}, b2 = {b_i[1]:.4f} cm3/mol")
print(f"                 a12 = sqrt(a1 a2)(1 - k12) = {AIJ[0,1]:.1f} "
      f"cm6 bar/mol2")


def mix_props(y, aij=None):
    aij = AIJ if aij is None else aij
    return float(y @ aij @ y), float(y @ b_i)


def ln_phi_i(y, T, P, want, aij=None):
    """Partial fugacity coefficients of both components in phase 'want'
    from PR with van der Waals one-fluid mixing rules."""
    aij = AIJ if aij is None else aij
    am, bm = mix_props(y, aij)
    A = am * P / (Rb * T) ** 2
    B = bm * P / (Rb * T)
    c = [1.0, -(1 - B), A - 2 * B - 3 * B ** 2, -(A * B - B ** 2 - B ** 3)]
    r = np.roots(c)
    r = np.sort(np.real(r[np.abs(np.imag(r)) < 1e-9]))
    r = r[r > B]
    Z = r[0] if want == "L" else r[-1]
    L = np.log((Z + (1 + S2) * B) / (Z + (1 - S2) * B))
    out = np.empty(2)
    for i in range(2):
        sm = float(aij[i] @ y)
        out[i] = (b_i[i] / bm * (Z - 1) - np.log(Z - B)
                  - A / (2 * S2 * B) * (2 * sm / am - b_i[i] / bm) * L)
    return out, Z, A, B, am, bm


def RR(psi, K, z):
    return float(np.sum(z * (K - 1) / (1 + psi * (K - 1))))


def dRR(psi, K, z):
    return float(-np.sum(z * (K - 1) ** 2 / (1 + psi * (K - 1)) ** 2))


def flash(P, T, z, aij=None, tol=1e-12, nmax=500):
    """Successive-substitution two-phase flash; returns (psi, K, x, y)."""
    K = Pc8 / P * np.exp(5.373 * (1 + w8) * (1 - Tc8 / T))
    psi = 0.5
    for _ in range(nmax):
        if RR(0.0, K, z) <= 0 or RR(1.0, K, z) >= 0:
            return np.nan, K, None, None      # single phase
        psi = brentq(RR, 0.0, 1.0, args=(K, z), xtol=1e-15)
        x = z / (1 + psi * (K - 1))
        y = K * x
        lL = ln_phi_i(x, T, P, "L", aij)[0]
        lV = ln_phi_i(y, T, P, "V", aij)[0]
        Kn = np.exp(lL - lV)
        if np.max(np.abs(Kn / K - 1)) < tol:
            K = Kn
            break
        K = Kn
    psi = brentq(RR, 0.0, 1.0, args=(K, z), xtol=1e-15)
    x = z / (1 + psi * (K - 1))
    return psi, K, x, K * x


# --- step 1: Wilson initialisation -----------------------------------
K_w = Pc8 / P8 * np.exp(5.373 * (1 + w8) * (1 - Tc8 / T8))
print(f"\nWilson initialisation K_i = (Pc_i/P) exp[5.373(1+w_i)(1 - Tc_i/T)]:")
print(f"  K1 = {K_w[0]:.6f}, K2 = {K_w[1]:.6f}")
psi_w = brentq(RR, 0.0, 1.0, args=(K_w, z8), xtol=1e-15)
x_w = z8 / (1 + psi_w * (K_w - 1))
y_w = K_w * x_w
print(f"  Rachford-Rice with these K: psi = {psi_w:.6f}")
print(f"  x = [{x_w[0]:.6f}, {x_w[1]:.6f}],  y = [{y_w[0]:.6f}, {y_w[1]:.6f}]")

# --- step 2: PR fugacity coefficients and K-values --------------------
lpL, ZL8, AL8, BL8, amL, bmL = ln_phi_i(x_w, T8, P8, "L")
lpV, ZV8, AV8, BV8, amV, bmV = ln_phi_i(y_w, T8, P8, "V")
print(f"\nLiquid trial phase : a_mix = {amL:.1f} cm6 bar/mol2, "
      f"b_mix = {bmL:.4f} cm3/mol")
print(f"                     A = {AL8:.6f}, B = {BL8:.6f}, Z_L = {ZL8:.6f}, "
      f"V_L = {ZL8*Rb*T8/P8:.4f} cm3/mol")
print(f"Vapour trial phase : a_mix = {amV:.1f} cm6 bar/mol2, "
      f"b_mix = {bmV:.4f} cm3/mol")
print(f"                     A = {AV8:.6f}, B = {BV8:.6f}, Z_V = {ZV8:.6f}, "
      f"V_V = {ZV8*Rb*T8/P8:.4f} cm3/mol")
print(f"\nphi_hat_i^L = [{np.exp(lpL[0]):.6f}, {np.exp(lpL[1]):.6f}]")
print(f"phi_hat_i^V = [{np.exp(lpV[0]):.6f}, {np.exp(lpV[1]):.6f}]")
K_pr = np.exp(lpL - lpV)
print(f"K_i = phi_hat_i^L/phi_hat_i^V = [{K_pr[0]:.6f}, {K_pr[1]:.6f}]")

# --- step 3: one Rachford-Rice Newton iteration -----------------------
F0 = RR(psi_w, K_pr, z8)
dF0 = dRR(psi_w, K_pr, z8)
psi_new = psi_w - F0 / dF0
print(f"\nRachford-Rice residual F(psi) = sum_i z_i (K_i - 1)/[1 + psi(K_i - 1)]")
print(f"  at psi = {psi_w:.6f}: term 1 = "
      f"{z8[0]*(K_pr[0]-1)/(1+psi_w*(K_pr[0]-1)):.6f}, term 2 = "
      f"{z8[1]*(K_pr[1]-1)/(1+psi_w*(K_pr[1]-1)):.6f}")
print(f"  F = {F0:.6f},  dF/dpsi = {dF0:.6f}")
print(f"  Newton step psi_new = {psi_w:.6f} - ({F0:.6f})/({dF0:.6f}) = "
      f"{psi_new:.6f}")
print(f"  bracketing test: F(0) = {RR(0.0, K_pr, z8):.6f}, "
      f"F(1) = {RR(1.0, K_pr, z8):.6f}")
print(f"  both positive and psi_new > 1, so no root of F lies in [0,1]: "
      f"the feed at {P8:.0f} bar is a single superheated vapour phase")

# --- checks -----------------------------------------------------------
for i, c in enumerate(comp):
    yy = np.zeros(2)
    yy[i] = 1.0
    lp1 = ln_phi_i(yy, T8, P8, "V")[0]
    lp2 = ln_phi("pr", c, T8, P8, Zroots("pr", c, T8, P8)[-1])
    print(f"[check] pure-component limit for {c}: mixture phi_hat = "
          f"{np.exp(lp1[i]):.8f}, pure PR phi = {np.exp(lp2):.8f}, "
          f"diff = {abs(np.exp(lp1[i])-np.exp(lp2)):.2e}")


def dew_res(P):
    """sum_i y_i/K_i - 1 at fixed vapour composition y = z."""
    K = Pc8 / P * np.exp(5.373 * (1 + w8) * (1 - Tc8 / T8))
    x = z8 / K
    for _ in range(500):
        x = x / np.sum(x)
        lL = ln_phi_i(x, T8, P, "L")[0]
        lV = ln_phi_i(z8, T8, P, "V")[0]
        Kn = np.exp(lL - lV)
        xn = z8 / Kn
        if np.max(np.abs(xn / np.sum(xn) - x)) < 1e-13:
            K = Kn
            break
        x, K = xn, Kn
    return float(np.sum(z8 / K) - 1.0)


P_dew = brentq(dew_res, 20.0, 30.0, xtol=1e-10)
print(f"[check] dew-point pressure of the equimolar feed at {T8:.0f} K = "
      f"{P_dew:.4f} bar; the specified {P8:.0f} bar lies below it, "
      f"confirming a single vapour phase")

P_two = 40.0
psi2, K2f, x2, y2 = flash(P_two, T8, z8)
print(f"[check] the same procedure at {P_two:.0f} bar, inside the two-phase "
      f"region, converges to psi = {psi2:.6f}, K = [{K2f[0]:.6f}, "
      f"{K2f[1]:.6f}], x = [{x2[0]:.6f}, {x2[1]:.6f}], "
      f"y = [{y2[0]:.6f}, {y2[1]:.6f}]")
lLc = ln_phi_i(x2, T8, P_two, "L")[0]
lVc = ln_phi_i(y2, T8, P_two, "V")[0]
print(f"[check] isofugacity residual max|x_i phi_i^L - y_i phi_i^V| = "
      f"{np.max(np.abs(x2*np.exp(lLc) - y2*np.exp(lVc))):.3e}; "
      f"material balance residual = "
      f"{np.max(np.abs(z8 - ((1-psi2)*x2 + psi2*y2))):.3e}")

# effect of a non-zero k12
print()
for kij in (0.02, 0.05):
    aij = aij_of(kij)
    lpLk = ln_phi_i(x_w, T8, P8, "L", aij)[0]
    lpVk = ln_phi_i(y_w, T8, P8, "V", aij)[0]
    Kk = np.exp(lpLk - lpVk)
    pk, Kf, _, _ = flash(P_two, T8, z8, aij)
    print(f"[check] k12 = {kij:.2f}: first-iterate K = [{Kk[0]:.6f}, "
          f"{Kk[1]:.6f}] ({(Kk[0]/K_pr[0]-1)*100:+.4f} %, "
          f"{(Kk[1]/K_pr[1]-1)*100:+.4f} % against k12 = 0); "
          f"converged psi at {P_two:.0f} bar = {pk:.6f} "
          f"({(pk/psi2-1)*100:+.4f} %)")


# =====================================================================
# EXAMPLE 9  Departure functions and a compressor duty
# =====================================================================
head(9, "PR departure functions for CO2 and an isentropic compressor duty")

sp9 = "CO2"
d9 = SPECIES[sp9]
T9, P9 = 350.0, 80.0
CPIG = 42.0            # J/mol/K, assumed constant ideal-gas heat capacity
NDOT9 = 100.0          # kmol/h
P_IN, P_OUT, T_IN = 10.0, 80.0, 350.0

Pc9 = d9["Pc"] * 1e5
ac9 = 0.45724 * R ** 2 * d9["Tc"] ** 2 / Pc9      # Pa m6 / mol2
b9 = 0.07780 * R * d9["Tc"] / Pc9                 # m3/mol
kap9 = kappa_pr(d9["w"])
print(f"Species        : {sp9}  Tc={d9['Tc']} K  Pc={d9['Pc']} bar  w={d9['w']}")
print(f"SI parameters  : ac = {ac9:.6e} Pa m6/mol2, b = {b9:.6e} m3/mol, "
      f"kappa = {kap9:.6f}")


def pr_ab(T):
    """PR a(T), da/dT and b in SI units for CO2."""
    Tr = T / d9["Tc"]
    s = 1 + kap9 * (1 - np.sqrt(Tr))
    a = ac9 * s ** 2
    dadT = -ac9 * kap9 * s / np.sqrt(T * d9["Tc"])
    return a, dadT, b9


def pr_state(T, P_bar):
    a, dadT, b = pr_ab(T)
    Pa = P_bar * 1e5
    A = a * Pa / (R * T) ** 2
    B = b * Pa / (R * T)
    Z = Zroots("pr", sp9, T, P_bar)[-1]
    L = np.log((Z + (1 + S2) * B) / (Z + (1 - S2) * B))
    return Z, A, B, L, a, dadT, b


def h_res(T, P_bar):
    Z, A, B, L, a, dadT, b = pr_state(T, P_bar)
    return R * T * (Z - 1) + (T * dadT - a) / (2 * S2 * b) * L


def s_res(T, P_bar):
    Z, A, B, L, a, dadT, b = pr_state(T, P_bar)
    return R * np.log(Z - B) + dadT / (2 * S2 * b) * L


Z9, A9, B9, L9, a9, dadT9, _ = pr_state(T9, P9)
print(f"\nState          : T = {T9:.1f} K, P = {P9:.1f} bar, "
      f"Tr = {T9/d9['Tc']:.4f}, Pr = {P9/d9['Pc']:.4f}")
print(f"a(T) = {a9:.6e} Pa m6/mol2, da/dT = {dadT9:.6e} Pa m6/mol2/K")
print(f"A = {A9:.6f}, B = {B9:.6f}, Z = {Z9:.6f}, "
      f"V = {Z9*R*T9/(P9*1e5)*1e6:.4f} cm3/mol")
print(f"L = ln[(Z+(1+sqrt2)B)/(Z+(1-sqrt2)B)] = {L9:.6f}")
print(f"(T da/dT - a)/(2 sqrt2 b) = "
      f"{(T9*dadT9-a9)/(2*S2*b9):.4f} J/mol")
print(f"(da/dT)/(2 sqrt2 b)       = {dadT9/(2*S2*b9):.6f} J/mol/K")
H9 = h_res(T9, P9)
S9 = s_res(T9, P9)
lnp9 = ln_phi("pr", sp9, T9, P9, Z9)
G9 = R * T9 * lnp9
print(f"\nH - H^ig = RT(Z-1) + (T da/dT - a)/(2 sqrt2 b) L = "
      f"{H9:.4f} J/mol = {H9/1000:.4f} kJ/mol")
print(f"S - S^ig = R ln(Z-B) + (da/dT)/(2 sqrt2 b) L    = "
      f"{S9:.6f} J/mol/K")
print(f"G - G^ig = RT ln phi = {G9:.4f} J/mol, phi = {np.exp(lnp9):.6f}")
print(f"[check] H_res - T S_res = {H9 - T9*S9:.4f} J/mol against "
      f"G_res = {G9:.4f} J/mol, difference {abs(H9-T9*S9-G9):.2e} J/mol")


def s_res_num(T, P_bar, h=0.05):
    """-d(G_res)/dT at constant P, by central difference on RT ln phi."""
    def g(TT):
        return R * TT * ln_phi("pr", sp9, TT, P_bar,
                               Zroots("pr", sp9, TT, P_bar)[-1])
    return -(g(T + h) - g(T - h)) / (2 * h)


S9n = s_res_num(T9, P9)
H9n = G9 + T9 * S9n
print(f"[check] S_res by numerical differentiation of G_res: {S9n:.4f} "
      f"J/mol/K (analytic {S9:.4f}, difference {abs(S9n-S9):.2e})")
print(f"[check] H_res = G_res + T S_res = {H9n:.4f} J/mol "
      f"(analytic {H9:.4f}, difference {abs(H9n-H9):.4f})")

# --- isentropic compression ------------------------------------------
print(f"\nCompressor     : {NDOT9:.0f} kmol/h CO2, {P_IN:.0f} bar to "
      f"{P_OUT:.0f} bar, inlet T = {T_IN:.0f} K")
print(f"Assumption     : constant ideal-gas heat capacity "
      f"Cp^ig = {CPIG:.1f} J/mol/K over the whole path")
s_in = s_res(T_IN, P_IN)
h_in = h_res(T_IN, P_IN)
print(f"Inlet          : Z = {Zroots('pr', sp9, T_IN, P_IN)[-1]:.6f}, "
      f"H_res = {h_in:.4f} J/mol, S_res = {s_in:.4f} J/mol/K")


def ds_total(T2):
    return (CPIG * np.log(T2 / T_IN) - R * np.log(P_OUT / P_IN)
            + s_res(T2, P_OUT) - s_in)


T2s = brentq(ds_total, T_IN, 1200.0, xtol=1e-10)
s_out = s_res(T2s, P_OUT)
h_out = h_res(T2s, P_OUT)
print(f"Solve Cp ln(T2/T1) - R ln(P2/P1) + S_res(T2,P2) - S_res(T1,P1) = 0")
print(f"  Cp ln(T2/T1) = {CPIG*np.log(T2s/T_IN):.4f}, "
      f"-R ln(P2/P1) = {-R*np.log(P_OUT/P_IN):.4f}, "
      f"S_res(out) - S_res(in) = {s_out-s_in:.4f} J/mol/K")
print(f"  -> T2s = {T2s:.4f} K, Z_out = "
      f"{Zroots('pr', sp9, T2s, P_OUT)[-1]:.6f}, "
      f"H_res(out) = {h_out:.4f} J/mol, S_res(out) = {s_out:.4f} J/mol/K")
print(f"  entropy balance residual at the solution = {ds_total(T2s):.2e} J/mol/K")

Ws = CPIG * (T2s - T_IN) + h_out - h_in
ndot_mol_s = NDOT9 * 1000.0 / 3600.0
print(f"In kJ/mol: H_res(in) = {h_in/1000:.4f}, H_res(out) = {h_out/1000:.4f}")
print(f"\nW_s = Cp(T2s - T1) + H_res(out) - H_res(in)")
print(f"    = {CPIG:.1f}*({T2s:.4f} - {T_IN:.1f}) + ({h_out:.4f}) - "
      f"({h_in:.4f}) = {Ws:.4f} J/mol")
print(f"Duty = W_s ndot = {Ws:.2f} J/mol * {ndot_mol_s:.4f} mol/s = "
      f"{Ws*ndot_mol_s/1000:.4f} kW")

T2_ig = T_IN * (P_OUT / P_IN) ** (R / CPIG)
W_ig = CPIG * (T2_ig - T_IN)
print(f"\nIdeal gas only : R/Cp = {R/CPIG:.6f}, "
      f"T2 = T1 (P2/P1)^(R/Cp) = {T_IN:.0f}*{P_OUT/P_IN:.0f}^"
      f"{R/CPIG:.6f} = {T2_ig:.4f} K")
print(f"                 W = Cp(T2 - T1) = {W_ig:.4f} J/mol, duty = "
      f"{W_ig*ndot_mol_s/1000:.4f} kW")
print(f"Comparison     : W_s - W_ig = {Ws - W_ig:+.4f} J/mol, "
      f"{(Ws-W_ig)*ndot_mol_s/1000:+.4f} kW, "
      f"{(Ws/W_ig - 1)*100:+.4f} % of the ideal-gas duty")
print(f"                 discharge temperature error of the ideal-gas route: "
      f"{T2_ig - T2s:+.4f} K")


# =====================================================================
# EXAMPLE 10  Corresponding states and the acentric factor in practice
# =====================================================================
head(10, "Corresponding states: omega from PR saturation; Z at Tr=1.2, Pr=1.5")

sp10 = "n-butane"
d10 = SPECIES[sp10]
T07 = 0.7 * d10["Tc"]
Ps07 = Psat_solve("pr", sp10, T07)
Pr07 = Ps07 / d10["Pc"]
w_rec = -1.0 - np.log10(Pr07)
print(f"Species        : {sp10}  Tc={d10['Tc']} K  Pc={d10['Pc']} bar, "
      f"tabulated w = {d10['w']:.3f}")
print(f"Tr = 0.7 gives T = 0.7*{d10['Tc']} = {T07:.4f} K")
print(f"PR Psat = {Ps07:.6f} bar, Pr_sat = {Ps07:.6f}/{d10['Pc']} = "
      f"{Pr07:.6f}")
print(f"w = -1 - log10(Pr_sat) = -1 - ({np.log10(Pr07):.6f}) = {w_rec:.6f}")
print(f"deviation from the tabulated value: {w_rec - d10['w']:+.6f} "
      f"({(w_rec/d10['w'] - 1)*100:+.4f} %)")

# --- Lee-Kesler three-parameter corresponding states -----------------
LK0 = dict(b1=0.1181193, b2=0.265728, b3=0.154790, b4=0.030323,
           c1=0.0236744, c2=0.0186984, c3=0.0, c4=0.042724,
           d1=0.155488e-4, d2=0.623689e-4, beta=0.65392, gam=0.060167)
LKR = dict(b1=0.2026579, b2=0.331511, b3=0.027655, b4=0.203488,
           c1=0.0313385, c2=0.0503618, c3=0.016901, c4=0.041577,
           d1=0.48736e-4, d2=0.0740336e-4, beta=1.226, gam=0.03754)
W_REF = 0.3978


def lk_Z(Tr, Pr, cst):
    """Lee-Kesler modified BWR solved for the ideal reduced volume."""
    B = cst["b1"] - cst["b2"] / Tr - cst["b3"] / Tr ** 2 - cst["b4"] / Tr ** 3
    C = cst["c1"] - cst["c2"] / Tr + cst["c3"] / Tr ** 3
    D = cst["d1"] + cst["d2"] / Tr

    def f(Vr):
        Z = (1 + B / Vr + C / Vr ** 2 + D / Vr ** 5
             + cst["c4"] / (Tr ** 3 * Vr ** 2)
             * (cst["beta"] + cst["gam"] / Vr ** 2)
             * np.exp(-cst["gam"] / Vr ** 2))
        return Z - Pr * Vr / Tr

    lo, hi = 0.05, 1.0e5
    Vr = brentq(f, lo, hi, xtol=1e-14, rtol=1e-14)
    return Pr * Vr / Tr, Vr


TR10, PR10 = 1.2, 1.5
Z0, Vr0 = lk_Z(TR10, PR10, LK0)
Zref, VrR = lk_Z(TR10, PR10, LKR)
Z1 = (Zref - Z0) / W_REF
Z_cs = Z0 + d10["w"] * Z1
print(f"\nThree-parameter corresponding states (Lee-Kesler) at "
      f"Tr = {TR10:.1f}, Pr = {PR10:.1f}:")
print(f"  simple fluid    Z0 = {Z0:.6f} (reduced volume Vr = {Vr0:.6f})")
print(f"  reference fluid Zr = {Zref:.6f} (n-octane, w_ref = {W_REF})")
print(f"  Z1 = (Zr - Z0)/w_ref = ({Zref:.6f} - {Z0:.6f})/{W_REF} = {Z1:.6f}")
print(f"  Z = Z0 + w Z1 = {Z0:.6f} + {d10['w']:.3f}*{Z1:.6f} = {Z_cs:.6f}")

T10 = TR10 * d10["Tc"]
P10 = PR10 * d10["Pc"]
Z_pr10 = Zroots("pr", sp10, T10, P10)[-1]
print(f"\nDirect PR at T = {T10:.4f} K, P = {P10:.4f} bar: Z = {Z_pr10:.6f}")
print(f"V(PR) = {Z_pr10*Rb*T10/P10:.4f} cm3/mol, "
      f"V(LK, three-parameter) = {Z_cs*Rb*T10/P10:.4f} cm3/mol, "
      f"V(LK, two-parameter) = {Z0*Rb*T10/P10:.4f} cm3/mol")
print(f"Difference     : Z(PR) - Z(LK) = {Z_pr10 - Z_cs:+.6f} "
      f"({(Z_pr10/Z_cs - 1)*100:+.4f} % of the Lee-Kesler value)")
print(f"Two-parameter estimate (w set to zero): Z = Z0 = {Z0:.6f}, "
      f"{(Z0/Z_cs - 1)*100:+.4f} % against the three-parameter value; the "
      f"acentric correction is worth {d10['w']*Z1:+.6f} in Z")

Zc0, _ = lk_Z(1.0, 1.0, LK0)
Zlow, _ = lk_Z(2.0, 1e-4, LK0)
B_lk = (LK0["b1"] - LK0["b2"] / 2.0 - LK0["b3"] / 4.0 - LK0["b4"] / 8.0)
B_num = (Zlow - 1.0) * 2.0 / 1e-4
print(f"\n[check] Lee-Kesler simple fluid at Tr = 2, Pr = 1e-4: Z0 = "
      f"{Zlow:.8f}, correctly approaching unity; the implied second virial "
      f"coefficient (Z0 - 1)Tr/Pr = {B_num:.6f} against the analytic "
      f"B(Tr = 2) = {B_lk:.6f} of the same correlation")
print(f"[check] Lee-Kesler simple fluid at its own critical point "
      f"(Tr = Pr = 1): Z0 = {Zc0:.6f} against the correlation's stated "
      f"Zc = 0.2901; the residual is the accuracy attainable at a triple "
      f"root")
Zc_num = np.mean(np.real(np.roots(
    [1.0, -(1 - 0.07780), 0.45724 - 2 * 0.07780 - 3 * 0.07780 ** 2,
     -(0.45724 * 0.07780 - 0.07780 ** 2 - 0.07780 ** 3)])))
print(f"[check] PR at its own critical point has a triple root whose mean is "
      f"(1 - B)/3 = {Zc_num:.6f} against the design value 0.30740")
B0ab = 0.083 - 0.422 / TR10 ** 1.6
B1ab = 0.139 - 0.172 / TR10 ** 4.2
Z_ab = 1 + (B0ab + d10["w"] * B1ab) * PR10 / TR10
print(f"[check] the Pitzer two-term virial correlation "
      f"(B0 = {B0ab:.6f}, B1 = {B1ab:.6f}) would give Z = {Z_ab:.6f}, but its "
      f"validity criterion Pr < Tr/2 = {TR10/2:.2f} is violated at "
      f"Pr = {PR10:.1f}, which is precisely why a full corresponding-states "
      f"correlation is needed here")
w_c3 = -1 - np.log10(Psat_solve("pr", "propane",
                                0.7 * SPECIES["propane"]["Tc"])
                     / SPECIES["propane"]["Pc"])
print(f"[check] the same recovery for propane returns w = {w_c3:.6f} against "
      f"a tabulated {SPECIES['propane']['w']:.3f}, so the small offset is a "
      f"property of the kappa correlation and not of one substance")

print("\n" + "=" * 78)
print("All ten examples completed.")
print("=" * 78)
