vle-thermo 0.16.0

//! Cubic equation of state variants.
//!
//! A cubic EOS expresses pressure as a function of temperature and molar volume
//! using a cubic polynomial in V. The general two-parameter form (Ref (5), Abbott)
//! is: P = RT/(V-b) - a·α(T)/((V+k₁b)(V+k₂b)), where k₁, k₂ are family constants
//! and α(T) is a temperature-dependent function specific to each variant.
//!
//! This module merges all 19 VB6 alpha-function variants (TADiPEDC enum) with
//! 3 additional Pascal models: Schmidt-Wenzel, Patel-Teja, and Patel-Teja USB.
//! The Pascal models are three-parameter EOS that add a third constant (c) to
//! better represent polar and heavy molecules — they require special C-parameter
//! mixing rules defined in the `mixing` module.
//!
//! Numeric discriminants preserve the legacy integer codes used in both VB6
//! `Select Case` and Pascal `case` dispatch. The 3-parameter models are assigned
//! codes 19–21 (continuing from VB6's 0–18 range) since the Pascal codes (6–8)
//! would collide with existing VB6 variants.
//!
//! # References
//! - (5) Abbott — general cubic EOS form with family constants
//! - (4) Da Silva & Báez (1989) — Schmidt-Wenzel, Patel-Teja, Chao-Seader

/// Cubic equation of state model selection.
///
/// Each variant corresponds to a specific α(Tr) function and EOS
/// parameterization. Two-parameter EOS (variants 0–18) use the Abbott
/// generalized form with k₁, k₂ family constants. Three-parameter EOS
/// (variants 19–21) add a c parameter for improved representation of
/// polar/asymmetric molecules.
// Automatically implement common traits for this enum:
//   Debug  — allows printing with {:?} for debugging
//   Clone  — allows creating a deep copy via .clone()
//   Copy   — allows implicit bitwise copies (no need to call .clone())
//   PartialEq, Eq — allows comparing values with == and !=
//   Hash   — allows using this enum as a key in HashMap/HashSet
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
// Store this enum in memory as a 32-bit integer (i32), matching the
// integer variant codes used by the legacy VB6 code (0..21). This also
// makes it safe to pass across the FFI boundary to C/Python.
#[repr(i32)]
// When built with the `python` feature (Maturin/PyO3 wheel build), expose this
// enum as a Python class. `eq` derives Python __eq__/__ne__ from Rust's
// PartialEq; `eq_int` lets Python compare variants against their integer
// codes. The cfg_attr means this annotation is invisible when the crate is
// used purely from Rust via `cargo add vle-thermo`.
#[cfg_attr(feature = "python", pyo3::pyclass(eq, eq_int))]
pub enum CubicEos {
    // --- Two-parameter cubic EOS (from VB6 TADiPEDC) ---
    /// Peng-Robinson (1976). The most widely used cubic EOS for hydrocarbon systems.
    /// α(Tr) = [1 + κ(1 - √Tr)]² where κ = f(ω). Ref (5).
    PR1976 = 0,
    /// Redlich-Kwong (1949). Original temperature-dependent cubic EOS.
    /// α(Tr) = 1/√Tr.
    RK1949 = 1,
    /// Soave-Redlich-Kwong (1972). Improved RK with acentric factor correlation.
    /// α(Tr) = [1 + m(1 - √Tr)]² where m = f(ω).
    RKS1972 = 2,
    /// van der Waals (1870). The first cubic EOS. α(Tr) = 1 (constant).
    VdW1870 = 3,
    /// Peng-Robinson-Lim (1997). Modified PR for improved vapor pressure prediction.
    PRL1997 = 4,
    /// RKS-Lim (1997). Modified RKS for improved vapor pressure prediction.
    RKSL1997 = 5,
    /// RKS-Graboski-Daubert (1978). Modified m(ω) correlation for RKS.
    RKSGD1978 = 6,
    /// Redlich-Prausnitz (1978).
    RP1978 = 7,
    /// Berthelot (1899). Temperature-modified van der Waals.
    Berth1899 = 8,
    /// van der Waals-Adachi (1984). Improved α function for van der Waals family.
    VdWAda1984 = 9,
    /// van der Waals-Valderrama (1989). Modified VdW for polar compounds.
    VdWVald1989 = 10,
    /// RKS-Mathias-Naumann (1980). Polar extension with m, n parameters.
    RKSmn1980 = 11,
    /// RKS-Adachi-Tagawa-Mathias-Naumann (1995). Three-constant α function.
    RKSATmn1995 = 12,
    /// PR-Adachi-Tagawa-Mathias-Naumann-Gasem (1997). Three-constant α for PR family.
    PRATmng1997 = 13,
    /// PR-Mathias-Massih-Naumann (1989). Polar extension with m, n parameters.
    PRMmn1989 = 14,
    /// Peng-Robinson-Stryjek-Vera (1986). Adds component-specific K₁ parameter.
    PRSV1986 = 15,
    /// van der Waals-OL (1998).
    VdWOL1998 = 16,
    /// Redlich-Kwong-OL (1998).
    RKOL1998 = 17,
    /// Peng-Robinson-OL (1998).
    PROL1998 = 18,

    // --- Three-parameter cubic EOS (from Pascal, Ref (4): Da Silva & Báez, 1989) ---
    /// Schmidt-Wenzel (3-parameter). Uses acentric factor ω as the third parameter
    /// to set the c constant. Beta parameter: 0.25989 - 0.02143ω + 0.00337ω².
    /// Ref (4), legacy/pascal/TERMOII.PAS.
    SchmidtWenzel = 19,
    /// Patel-Teja (3-parameter). Uses a fitted Zc correlation as the third parameter.
    /// C-parameter mixing: simple mole-fraction average c_mix = Σ xᵢcᵢ.
    /// Ref (4), legacy/pascal/TERMOII.PAS.
    PatelTeja = 20,
    /// Patel-Teja with Universal Saturation Behavior. Same EOS form as PatelTeja
    /// but uses √B-weighted C-parameter mixing: c_mix = Σ(xᵢ√Bᵢ·cᵢ)/Σ(xᵢ√Bᵢ).
    /// Ref (4), legacy/pascal/TERMOII.PAS.
    PatelTejaUSB = 21,
}

impl CubicEos {
    /// Returns `true` if this is a three-parameter EOS requiring a C-parameter
    /// mixing rule in addition to the standard a, b mixing.
    pub fn is_three_parameter(&self) -> bool {
        matches!(
            self,
            CubicEos::SchmidtWenzel | CubicEos::PatelTeja | CubicEos::PatelTejaUSB
        )
    }
}

/// Vapor-phase model selection.
///
/// In VLE calculations, the vapor phase fugacity can be computed from a cubic EOS,
/// the virial equation (truncated at second coefficient), or the ideal gas assumption.
/// This enum wraps [`CubicEos`] with those additional options.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum VaporModel {
    /// Ideal gas assumption (fugacity coefficient φᵢ = 1). Valid only at low pressures.
    IdealGas,
    /// Second virial equation using Pitzer B⁰/B¹ correlations.
    /// Good for moderate pressures where cubic EOS is overkill.
    Virial,
    /// Cubic EOS — use the associated [`CubicEos`] variant for fugacity calculation.
    Cubic(CubicEos),
}

/// Liquid-phase model selection.
///
/// The liquid phase can be modeled with the same EOS as the vapor (φ-φ approach),
/// with an activity coefficient model (γ-φ approach), or with the Chao-Seader
/// correlation for specific compound classes.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum LiquidModel {
    /// Ideal solution (Raoult's law: γᵢ = 1). Valid for chemically similar components.
    IdealSolution,
    /// Cubic EOS for liquid fugacity (φ-φ approach). Same EOS is used for both phases.
    Cubic(CubicEos),
    /// Activity coefficient model (γ-φ approach). Better for polar/non-ideal liquid mixtures.
    Activity(super::ActivityModel),
    /// Chao-Seader liquid fugacity correlation. Semi-empirical method with special
    /// handling for hydrogen and methane. Ref (4), legacy/pascal/TERMOII.PAS.
    ///
    /// **Legacy path**: `Kᵢ = νᵢ/φ̂ᵢⱽ` with the Grayson–Streed 1963 ν⁰ table
    /// (see [`chao_seader_ln_phi`]) and **no** regular-solution activity
    /// coefficient — exactly what the Pascal program computed. Prefer
    /// [`LiquidModel::GraysonStreed`] for anything new.
    ChaoSeader,
    /// Grayson–Streed (1963) refinery K-value method (Milestone 20):
    /// `Kᵢ = νᵢ·γᵢ/φ̂ᵢⱽ` — the Grayson–Streed νᵢ
    /// ([`regular_solution_ln_nu`]), the Scatchard–Hildebrand regular-solution
    /// γᵢ built from each component's `solubility_param` and `liquid_volume`
    /// (or the `SystemSpec`'s `delta` / `vl` overrides), and the vapor φ̂ᵢⱽ from
    /// whatever `VaporModel` is configured (Redlich–Kwong classically).
    /// Hydrogen and methane are recognised by name
    /// ([`ChaoSeaderSpecies::for_component`]). Both νᵢ and the γ constants are
    /// composition-independent, so the flash caches them per `(T, P)`.
    GraysonStreed,
    /// Braun K10 (BK10) refinery K-value method (Milestone 20): the K-value of
    /// a heavy fraction from its Maxwell–Bonnell vapor pressure,
    /// `Kᵢ = Pᵢᴹᴮ(T; Tb,ᵢ)/P`, corrected by the vapor φ̂ᵢⱽ of the configured
    /// `VaporModel` (ideal gas for the textbook method). Needs `Component::tb`
    /// for every component; the Watson-K correction uses `Component::watson_k`
    /// when it is set. Braun's pressure-correction charts are **not**
    /// implemented — this is the K10 value scaled Raoult-style to the system
    /// pressure, which is what the method's low-pressure (< ~10 bar) validity
    /// range assumes anyway.
    BraunK10,
}

/// Phase identifier used to select liquid or vapor root from the cubic solver.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
#[cfg_attr(feature = "python", pyo3::pyclass(eq, eq_int))]
#[repr(i32)]
pub enum PhaseId {
    Vapor = 0,
    Liquid = 1,
}

// ===========================================================================
// M7.1 / M7.2 — pure-component cubic-EOS algorithms.
//
// This block fills in the family constants, alpha functions, Z-factor,
// fugacity, and departure-property machinery for the two-parameter cubic
// EOS catalogue. M7.1 shipped the *deployable core* (PR1976, RKS1972,
// RK1949, VdW1870 — the four variants Chapter IV uses); M7.2 added the
// remaining twelve two-parameter variants (Berthelot, VdWAda1984,
// RKSGD1978, RKSL1997, RP1978, PRL1997, VdWVald1989, RKSmn1980,
// RKSATmn1995, PRATmng1997, PRMmn1989, PRSV1986), each with an analytical
// dα/dTr.
//
// The remaining CubicEos variants still panic via `unimplemented!` with a
// pointer to their deferred sub-milestone: the OL family (VdWOL1998,
// RKOL1998, PROL1998) is M7.4 because its α is coupled to the reduced
// saturation pressure, and the 3-param Pascal EOS is M7.3. The panic
// message names both the variant and the legacy source line, so a future
// porting session has the receipt it needs to fill in the gap.
//
// Convention used throughout this block (Abbott form, matching VB6):
//
//     P = R·T/(V − b)  −  a·α(T) / (V² + k1·b·V + k2·b²)
//
// with a = OmA · R²·Tc²/Pc and b = OmB · R·Tc/Pc.
// ===========================================================================

use crate::numerics::cubic::{CubicError, solve_real};
use crate::types::Component;
use num_dual::DualNum;
use thiserror::Error;

/// Errors raised by the cubic-EOS pure-component layer.
///
/// `Cubic` passes through any failure from the underlying cubic solver
/// (NaN coefficient, near-zero leading coefficient). `NoRootForPhase` is
/// thrown when the cubic has only one physical root and the caller asked
/// for the wrong phase (e.g., requesting Liquid at supercritical T).
/// `NotImplemented` is thrown for EOS variants whose α has not yet been
/// ported — exactly the variants the deferred sub-milestones cover.
#[derive(Debug, Error, PartialEq)]
pub enum EosError {
    /// The cubic solver itself returned an error (non-finite coefficient
    /// or near-zero leading term).
    #[error("cubic solver failed: {0}")]
    Cubic(#[from] CubicError),

    /// The cubic has no physical root for the requested phase. Usually
    /// means the caller asked for Liquid above the critical point (or
    /// vice versa) — the EOS has only one real root in that region.
    #[error("no real root above B={big_b:.6e} found for phase {phase:?}")]
    NoRootForPhase { phase: PhaseId, big_b: f64 },

    /// The EOS variant's α function has not been ported yet. The variant
    /// in the payload is the exact enum that triggered the panic.
    #[error("EOS variant {0:?} not yet ported — see M7 sub-milestones")]
    NotImplemented(CubicEos),
}

/// Two-parameter cubic-EOS family constants for the Abbott form.
///
/// `k1` and `k2` parameterize the attractive term's denominator
/// `(V² + k1·b·V + k2·b²)`. `om_a` and `om_b` are the Ω_a / Ω_b
/// dimensionless coefficients that multiply `R²·Tc²/Pc` and `R·Tc/Pc`
/// to get the EOS `a` and `b` parameters.
///
/// Constants are stored as `f64` even though VB6 declares them as
/// integers — the discriminant `k1² − 4·k2` is computed in f64 and
/// matching types saves a noisy cast at every call site.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct FamilyConstants {
    /// Coefficient of `b·V` in the attractive-term denominator.
    pub k1: f64,
    /// Coefficient of `b²` in the attractive-term denominator.
    pub k2: f64,
    /// Ω_a — dimensionless coefficient for `a = Ω_a · R²·Tc²/Pc`.
    pub om_a: f64,
    /// Ω_b — dimensionless coefficient for `b = Ω_b · R·Tc/Pc`.
    pub om_b: f64,
}

/// Return the 2-parameter family constants for the given EOS.
///
/// Source: `legacy/vb6/McommonFunctions.bas:273` (`GeneralConstantsEOS`).
/// 3-parameter EOS (Schmidt-Wenzel, Patel-Teja) don't use this table — their
/// Ω_a/Ω_b and the cubic are built in `three_param_aubw`; they return
/// zeros here so callers can still pattern-match without an early panic;
/// the per-component override happens inside their own EOS code path.
pub fn family_constants(eos: CubicEos) -> FamilyConstants {
    use CubicEos::*;
    match eos {
        // PR family: k1=2, k2=-1 → denominator V² + 2bV − b² = (V+(1+√2)b)(V+(1−√2)b).
        PR1976 | RP1978 | PRL1997 | PRATmng1997 | PRSV1986 | PROL1998 | PRMmn1989 => {
            FamilyConstants {
                k1: 2.0,
                k2: -1.0,
                om_a: 0.457235528921382,
                om_b: 0.0777960739038885,
            }
        }
        // RKS family: k1=1, k2=0 → denominator V² + bV = V(V+b).
        RKS1972 | RKSL1997 | RKSGD1978 | RK1949 | VdWVald1989 | RKSATmn1995 | RKOL1998
        | RKSmn1980 => FamilyConstants {
            k1: 1.0,
            k2: 0.0,
            om_a: 0.427480233540341,
            om_b: 0.0866403499649577,
        },
        // VdW family: k1=k2=0 → denominator V². Ω values are the analytical
        // 27/64 and 1/8 (no floating-point rounding involved).
        VdW1870 | Berth1899 | VdWAda1984 | VdWOL1998 => FamilyConstants {
            k1: 0.0,
            k2: 0.0,
            om_a: 27.0 / 64.0,
            om_b: 1.0 / 8.0,
        },
        // 3-parameter EOS: the (Ω_a, Ω_b, cubic) are built per-component in
        // `three_param_aubw`, not from this table. Return zeros so a caller
        // that pattern-matches this function doesn't get an early panic.
        SchmidtWenzel | PatelTeja | PatelTejaUSB => FamilyConstants {
            k1: 0.0,
            k2: 0.0,
            om_a: 0.0,
            om_b: 0.0,
        },
    }
}

// ===========================================================================
// Three-parameter EOS constant helpers (M7.3) — Ref (4): Da Silva & Báez
// (1989), legacy/pascal/TERMOII.PAS. Each helper is a pure ω→constant
// correlation; the α functions keep the convention α(Tr=1)=1 by folding the
// EOS-specific prefactor into Ω_a (see `three_param_aubw`).
// ===========================================================================

/// Patel-Teja α-shape parameter F(ω). TERMOII.PAS:175 (`Numi`).
fn pt_f(w: f64) -> f64 {
    0.452413 + 1.30982 * w - 0.295937 * w * w
}
/// Patel-Teja fitted critical compressibility ξc(ω). TERMOII.PAS:176 (`ZZc`).
fn pt_xi_c(w: f64) -> f64 {
    0.329032 - 0.076799 * w + 0.0211947 * w * w
}
/// Patel-Teja Ω_b(ω). TERMOII.PAS:177 (`Num2i`).
fn pt_om_b(w: f64) -> f64 {
    0.08517138 - 0.02640592 * w + 0.00788769 * w * w
}
/// Patel-Teja Ω_a(ω) = 3ξc² + 3(1−2ξc)Ω_b + Ω_b² + 1 − 3ξc. TERMOII.PAS:178 (`FactorPT`).
fn pt_om_a(w: f64) -> f64 {
    let xc = pt_xi_c(w);
    let ob = pt_om_b(w);
    3.0 * xc * xc + 3.0 * (1.0 - 2.0 * xc) * ob + ob * ob + 1.0 - 3.0 * xc
}

/// Schmidt-Wenzel β(ω) third-parameter coefficient. TERMOII.PAS:167 (`Num2i`).
fn sw_beta(w: f64) -> f64 {
    0.25988221 - 0.02142913 * w + 0.00337143 * w * w
}
/// Schmidt-Wenzel base m₀(ω). TERMOII.PAS:168.
fn sw_m0(w: f64) -> f64 {
    0.465 + 1.347 * w - 0.528 * w * w
}
/// Schmidt-Wenzel m(Tr) — **piecewise** (TERMOII.PAS:169-170). The two
/// branches agree in value at Tr=1 but **not in slope** (a documented kink;
/// see `sw_dm_dtr`).
fn sw_m(w: f64, tr: f64) -> f64 {
    let m0 = sw_m0(w);
    if tr <= 1.0 {
        let g = 5.0 * tr - 3.0 * m0 - 1.0;
        m0 + g * g / 70.0
    } else {
        let g = 4.0 - 3.0 * m0;
        m0 + g * g / 70.0
    }
}
/// dm/dTr for Schmidt-Wenzel — the one-sided derivative of the active branch.
/// Faithful + guarded: the Tr>1 branch has m constant (m'=0); the Tr≤1 branch
/// is quadratic. The slope is discontinuous at Tr=1 — callers at exactly Tr=1
/// get the Tr≤1 (left) derivative, which keeps the entropy departure finite
/// (the legacy returned NaN there; TERMOII.PAS:492).
fn sw_dm_dtr(w: f64, tr: f64) -> f64 {
    if tr <= 1.0 {
        let m0 = sw_m0(w);
        let g = 5.0 * tr - 3.0 * m0 - 1.0;
        2.0 * g * 5.0 / 70.0
    } else {
        0.0
    }
}
/// Schmidt-Wenzel Ω_b(ω) = β / (3(1+βω)). TERMOII.PAS:195.
fn sw_om_b(w: f64) -> f64 {
    let b = sw_beta(w);
    b / (3.0 * (1.0 + b * w))
}
/// Schmidt-Wenzel Ω_a(ω) = [1 − (1−β)/(3(1+βω))]³ — the constant prefactor
/// folded out of α so that α(Tr=1)=1. TERMOII.PAS:173.
fn sw_om_a(w: f64) -> f64 {
    let b = sw_beta(w);
    let inner = 1.0 - (1.0 - b) / (3.0 * (1.0 + b * w));
    inner * inner * inner
}

// ===========================================================================
// OL-family α (M7.4) — Olivera et al. (1998), legacy/vb6/clsQbicsPure.cls:268.
// α = Tr·(1 + SumHk), where SumHk reads the component's *reduced saturation
// pressure* (via comp.sat_model), so these variants are coupled to the
// saturation layer. The OL EOS themselves use the standard 2-parameter
// family constants (VdWOL1998→VdW, RKOL1998→RKS, PROL1998→PR).
// ===========================================================================

/// Per-family SumHk coefficient table `(h[1..10], E)`.
fn ol_coeffs(eos: CubicEos) -> ([f64; 10], f64) {
    use CubicEos::*;
    match eos {
        VdWOL1998 => (
            [
                0.33333333,
                0.35112597,
                0.011287433,
                -0.0038485685,
                0.00064261934,
                -0.000067252383,
                0.0000045962725,
                -0.00000019990875,
                5.0318465e-09,
                -5.5827084e-11,
            ],
            -1.02,
        ),
        RKOL1998 => (
            [
                0.32748,
                0.34376954,
                0.010596403,
                -0.0037538497,
                0.00063257197,
                -0.000066481,
                0.0000045517956,
                -0.00000019796921,
                4.9748592e-09,
                -5.5024228e-11,
            ],
            -1.014,
        ),
        PROL1998 => (
            [
                0.29803582,
                0.015003698,
                -0.0047527103,
                0.0008036716,
                -0.000089548695,
                0.0000068691611,
                -0.00000036067317,
                0.000000012409205,
                -2.5222671e-10,
                2.2955503e-12,
            ],
            -0.00041,
        ),
        _ => unreachable!("ol_coeffs called for non-OL EOS {eos:?}"),
    }
}

/// OL SumHk(Tr) = (1 − L)^E · Σ hₖ·(−L)^k, L = ln(P_sat,r / Tr). Returns `None`
/// if the component's saturation model can't evaluate (e.g. missing coeffs).
fn ol_sumhk(eos: CubicEos, tr: f64, comp: &Component) -> Option<f64> {
    let (h, e_exp) = ol_coeffs(eos);
    let t = tr * comp.tc;
    let pr = crate::saturation::reduced_psat(comp.sat_model, comp, t).ok()?;
    let arg = -(pr / tr).ln(); // −L
    let mut sum = 0.0;
    let mut argk = 1.0; // arg^0
    for hk in h {
        argk *= arg; // arg^k
        sum += hk * argk;
    }
    Some((1.0 + arg).powf(e_exp) * sum)
}

/// OL (SumHk, dSumHk/dTr). dα/dTr is analytical given dP_sat/dT (which is
/// analytical for the Antoine sat model, numerical otherwise).
fn ol_sumhk_and_deriv(eos: CubicEos, tr: f64, comp: &Component) -> Option<(f64, f64)> {
    let (h, e_exp) = ol_coeffs(eos);
    let t = tr * comp.tc;
    let pr = crate::saturation::reduced_psat(comp.sat_model, comp, t).ok()?;
    let dpsat_dt = crate::saturation::d_psat_dt(comp.sat_model, comp, t).ok()?;
    // arg = −L = ln(Tr) − ln(P_sat,r); d(arg)/dTr = 1/Tr − (dP_sat,r/dTr)/P_sat,r,
    // P_sat,r = Psat/Pc, dP_sat,r/dTr = (Tc/Pc)·dPsat/dT.
    let dpr_dtr = (dpsat_dt / comp.pc) * comp.tc;
    let arg = -(pr / tr).ln();
    let darg_dtr = 1.0 / tr - dpr_dtr / pr;
    // P_sum = Σ hₖ argᵏ; dP_sum/darg = Σ k·hₖ arg^(k−1).
    let mut psum = 0.0;
    let mut dpsum = 0.0;
    let mut argk_minus1 = 1.0; // arg^(k−1)
    for (i, hk) in h.iter().enumerate() {
        let k = (i + 1) as f64;
        psum += hk * argk_minus1 * arg;
        dpsum += hk * k * argk_minus1;
        argk_minus1 *= arg;
    }
    let base = (1.0 + arg).powf(e_exp);
    let dbase = e_exp * (1.0 + arg).powf(e_exp - 1.0);
    let s = base * psum;
    let ds_dtr = (dbase * psum + base * dpsum) * darg_dtr;
    Some((s, ds_dtr))
}

/// α(Tr) — the temperature-dependent multiplier on the attractive term.
///
/// # Arguments
/// * `eos` — EOS variant.
/// * `tr` — Reduced temperature T/Tc. **Dimensionless.**
/// * `comp` — Component data. The acentric-factor variants read only
///   `omega`; the polar/fitted variants additionally read `zc`
///   (VdWVald1989), `m_polar`/`n_polar` (RKSmn1980, PRMmn1989),
///   `m_polar`/`n_polar`/`g_polar` (RKSATmn1995, PRATmng1997), or
///   `prsv_k1` (PRSV1986).
///
/// # Returns
/// α evaluated at `tr`, **dimensionless**.
///
/// All 22 variants are implemented. The OL family (saturation-coupled) reads
/// the component's `sat_model`; if that model can't evaluate (e.g. missing
/// coefficients) the OL α returns `NaN` rather than panicking, so a downstream
/// `z_factor` then fails cleanly instead of producing a bogus root.
pub fn alpha(eos: CubicEos, tr: f64, comp: &Component) -> f64 {
    use CubicEos::*;
    let w = comp.omega;
    match eos {
        // ----- Deployable core (M7.1) -----
        VdW1870 => 1.0,
        RK1949 => 1.0 / tr.sqrt(),
        RKS1972 => {
            // α = [1 + m·(1 − √Tr)]² with m = 0.48 + 1.574ω − 0.176ω².
            // VB6: clsQbicsPure.cls:1734
            let m = 0.48 + 1.574 * w - 0.176 * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + m * s).powi(2)
        }
        PR1976 => {
            // α = [1 + κ·(1 − √Tr)]² with κ = 0.37464 + 1.54226ω − 0.26992ω².
            // VB6: clsQbicsPure.cls:1741
            let kappa = 0.37464 + 1.54226 * w - 0.26992 * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + kappa * s).powi(2)
        }
        // ----- M7.2: remaining 2-parameter α variants -----
        // All ported from VB6 `Friend Function Alpha` (clsQbicsPure.cls:1719).

        // Berthelot (1899): a temperature-modified VdW. VB6:1730.
        // α = 1/Tr.
        Berth1899 => 1.0 / tr,

        // van der Waals-Adachi (1984). VB6:1732. α = 10^(m·(1 − Tr)).
        // The VB6 wrote m's linear term in two pieces
        // (0.791981·ω + 0.654505·ω); we fold them into the single
        // coefficient 1.446486·ω — a bit-identical result. (The split
        // is almost certainly a transcription artifact of an intended
        // higher-order term in the original correlation.)
        VdWAda1984 => {
            let m = 0.228165 + 1.446486 * w - 0.648552 * w * w;
            10f64.powf(m * (1.0 - tr))
        }

        // Soave-shaped variants: α = [1 + m(ω)·(1 − √Tr)]², identical in
        // form to RKS1972/PR1976 but with different m(ω) correlations.
        // RKSGD1978 = Graboski-Daubert (VB6:1738).
        RKSGD1978 => {
            let m = 0.48508 + 1.55171 * w - 0.15613 * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + m * s).powi(2)
        }
        // RKSL1997 = Lim's RKS modification, cubic in ω (VB6:1740).
        RKSL1997 => {
            let m = 0.478972559 + 1.576809191 * w - 0.187219516 * w * w + 0.020424946 * w * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + m * s).powi(2)
        }
        // RP1978 = Redlich-Prausnitz, PR family, cubic in ω (VB6:1744).
        RP1978 => {
            let m = 0.379642 + 1.48503 * w - 0.164423 * w * w + 0.016666 * w * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + m * s).powi(2)
        }
        // PRL1997 = Peng-Robinson-Lim, cubic in ω (VB6:1746).
        PRL1997 => {
            let m = 0.378710697 + 1.487972964 * w - 0.166754831 * w * w + 0.017169486 * w * w * w;
            let s = 1.0 - tr.sqrt();
            (1.0 + m * s).powi(2)
        }

        // Mathias-Naumann linear form: α = 1 + (1 − Tr)·(m + n/Tr).
        // VdWVald1989 derives m, n from ω·Zc (VB6:1748); RKSmn1980 reads
        // the component-specific fitted m, n directly (VB6:1753).
        VdWVald1989 => {
            let omegac = w * comp.zc;
            let m = 0.4745 + (2.7349 + 6.0984 * omegac) * omegac;
            let n = 0.0674 + (2.1031 + 3.9512 * omegac) * omegac;
            1.0 + (1.0 - tr) * (m + n / tr)
        }
        RKSmn1980 => {
            let (m, n) = (comp.m_polar, comp.n_polar);
            1.0 + (1.0 - tr) * (m + n / tr)
        }

        // Adachi-Tagawa-Mathias-Naumann exponential form with three fitted
        // constants (m, n, g). The expression is identical for the RKS and
        // PR families — only the family k1/k2 constants differ. VB6:1755/1757.
        // α = exp[ (1 − Tr)·m·|1 − Tr|^(g−1) + n·(1/Tr − 1) ].
        RKSATmn1995 | PRATmng1997 => {
            let (m, n, g) = (comp.m_polar, comp.n_polar, comp.g_polar);
            let u = 1.0 - tr;
            (u * m * u.abs().powf(g - 1.0) + n * (1.0 / tr - 1.0)).exp()
        }

        // PR-Mathias-Massih-Naumann (1989): α = exp[ (1 − Tr)·m + n·(1 − √Tr)² ].
        // m, n are component-specific fitted constants. VB6:1759.
        PRMmn1989 => {
            let (m, n) = (comp.m_polar, comp.n_polar);
            let s = 1.0 - tr.sqrt();
            ((1.0 - tr) * m + n * s * s).exp()
        }

        // Peng-Robinson-Stryjek-Vera (1986): α = [1 + κ·(1 − √Tr)]² with
        // κ = κ₀(ω) + K₁·(1 + √Tr)·(0.7 − Tr). K₁ is the component-specific
        // PRSV parameter (zero recovers the plain PR-like κ₀ form). VB6:1762.
        PRSV1986 => {
            let r = tr.sqrt();
            let kappa0 = 0.378893 + 1.4897153 * w - 0.17131848 * w * w + 0.0196554 * w * w * w;
            let kappa = kappa0 + comp.prsv_k1 * (1.0 + r) * (0.7 - tr);
            let inner = 1.0 + kappa * (1.0 - r);
            inner * inner
        }

        // ----- Deferred: OL family → M7.4 -----
        // The OL-family α is not a function of (Tr, ω): VB6 computes
        // α = Tr·(1 + SumHk) where SumHk depends on the component's reduced
        // saturation pressure (clsQbicsPure.cls:268). It is therefore coupled
        // to the saturation layer and lands with M7.4 alongside the
        // non-Antoine saturation models.
        VdWOL1998 | RKOL1998 | PROL1998 => {
            // OL-family: α = Tr·(1 + SumHk); SumHk reads the reduced saturation
            // pressure via comp.sat_model. NaN if that model can't evaluate
            // (e.g. missing coeffs) — downstream z_factor then errors cleanly.
            ol_sumhk(eos, tr, comp)
                .map(|s| tr * (1.0 + s))
                .unwrap_or(f64::NAN)
        }
        // ----- M7.3: 3-parameter Pascal EOS (Ref (4), TERMOII.PAS) -----
        // α is the pure temperature function (α(Tr=1)=1); the EOS-specific
        // constant prefactors live in `three_param_aubw` (Ω_a/Ω_b).
        SchmidtWenzel => {
            // α = [1 + m(Tr)·(1 − √Tr)]², m piecewise. TERMOII.PAS:171-173.
            let s = 1.0 - tr.sqrt();
            let m = sw_m(w, tr);
            (1.0 + m * s).powi(2)
        }
        PatelTeja | PatelTejaUSB => {
            // Soave-shaped: α = [1 + F·(1 − √Tr)]². PatelTeja and the USB
            // variant share the same pure-component α; they differ only in
            // the mixture C-parameter rule (M8). TERMOII.PAS:179.
            let f = pt_f(w);
            let s = 1.0 - tr.sqrt();
            (1.0 + f * s).powi(2)
        }
    }
}

/// Schmidt-Wenzel `m(ω, Tr)` generic over the scalar type (M12.3).
///
/// Same piecewise form as [`sw_m`] but evaluated with a `DualNum` `tr`, so
/// duals seeded on temperature propagate through the Schmidt-Wenzel α. Branch
/// selection uses `tr.re()` (the value path), matching [`sw_m`]; each piece is
/// smooth in its interior.
fn sw_m_generic<D: DualNum<f64> + Copy>(w: f64, tr: D) -> D {
    let m0 = sw_m0(w);
    if tr.re() <= 1.0 {
        // g = 5·Tr − 3·m0 − 1; m = m0 + g²/70.
        let g = tr * 5.0 - (3.0 * m0 + 1.0);
        g * g * (1.0 / 70.0) + m0
    } else {
        // Tr > 1: g is constant, so m is a constant in Tr (zero slope).
        let g = 4.0 - 3.0 * m0;
        D::from(m0 + g * g / 70.0)
    }
}

/// α(Tr) generic over the scalar type `D: DualNum<f64>` (M12.3).
///
/// This is the temperature-generic sibling of [`alpha`]: with `D = f64` it
/// reproduces `alpha` (to floating-point rounding — one variant, VdWAda1984, is
/// re-expressed as `exp(y·ln 10)` instead of `10^y`, differing only in the last
/// ULP), and with a dual type it carries exact dα/dT along, which is what the
/// T- and P-derivative dual paths (`mixture::d_ln_phi_d_t` / `_d_p`, M12.3) rely
/// on. `tr` is the reduced temperature T/Tc (**dimensionless**); the return is
/// α, **dimensionless**.
///
/// # Coverage
/// Every variant except the OL family is genuinely generic (pure `DualNum`
/// arithmetic). The saturation-coupled OL family (VdWOL1998/RKOL1998/PROL1998)
/// is **value-lifted** from the scalar [`alpha`] with a zero derivative slot —
/// its α depends on the reduced saturation pressure through the saturation
/// layer, which is not part of this generic path. OL is out of scope for the
/// derivative APIs (they gate on classical/2-parameter or dual over the
/// GE-cubic path), so this is a deliberate, documented limitation, not a
/// silent one; do not route OL through a T-derivative and trust the slope.
pub fn alpha_generic<D: DualNum<f64> + Copy>(eos: CubicEos, tr: D, comp: &Component) -> D {
    use CubicEos::*;
    let w = comp.omega;
    // Soave-shaped α = [1 + m·(1 − √Tr)]² helper (RKS/PR/Graboski/Lim/...).
    let soave = |m: f64| -> D {
        let s = -tr.sqrt() + 1.0; // 1 − √Tr
        (s * m + 1.0).powi(2)
    };
    match eos {
        VdW1870 => D::from(1.0),
        RK1949 => tr.sqrt().recip(),
        RKS1972 => soave(0.48 + 1.574 * w - 0.176 * w * w),
        PR1976 => soave(0.37464 + 1.54226 * w - 0.26992 * w * w),
        Berth1899 => tr.recip(),
        VdWAda1984 => {
            // α = 10^(m·(1 − Tr)) = exp(ln(10)·m·(1 − Tr)).
            let m = 0.228165 + 1.446486 * w - 0.648552 * w * w;
            ((-tr + 1.0) * (m * std::f64::consts::LN_10)).exp()
        }
        RKSGD1978 => soave(0.48508 + 1.55171 * w - 0.15613 * w * w),
        RKSL1997 => {
            soave(0.478972559 + 1.576809191 * w - 0.187219516 * w * w + 0.020424946 * w * w * w)
        }
        RP1978 => soave(0.379642 + 1.48503 * w - 0.164423 * w * w + 0.016666 * w * w * w),
        PRL1997 => {
            soave(0.378710697 + 1.487972964 * w - 0.166754831 * w * w + 0.017169486 * w * w * w)
        }
        VdWVald1989 => {
            // α = 1 + (1 − Tr)·(m + n/Tr).
            let omegac = w * comp.zc;
            let m = 0.4745 + (2.7349 + 6.0984 * omegac) * omegac;
            let n = 0.0674 + (2.1031 + 3.9512 * omegac) * omegac;
            (-tr + 1.0) * (tr.recip() * n + m) + 1.0
        }
        RKSmn1980 => {
            let (m, n) = (comp.m_polar, comp.n_polar);
            (-tr + 1.0) * (tr.recip() * n + m) + 1.0
        }
        RKSATmn1995 | PRATmng1997 => {
            // α = exp[(1 − Tr)·m·|1 − Tr|^(g−1) + n·(1/Tr − 1)].
            let (m, n, g) = (comp.m_polar, comp.n_polar, comp.g_polar);
            let u = -tr + 1.0;
            (u * m * u.abs().powf(g - 1.0) + (tr.recip() - 1.0) * n).exp()
        }
        PRMmn1989 => {
            // α = exp[(1 − Tr)·m + n·(1 − √Tr)²].
            let (m, n) = (comp.m_polar, comp.n_polar);
            let s = -tr.sqrt() + 1.0;
            ((-tr + 1.0) * m + s * s * n).exp()
        }
        PRSV1986 => {
            // α = [1 + κ·(1 − √Tr)]², κ = κ₀(ω) + K₁·(1 + √Tr)·(0.7 − Tr).
            let r = tr.sqrt();
            let kappa0 = 0.378893 + 1.4897153 * w - 0.17131848 * w * w + 0.0196554 * w * w * w;
            let kappa = (r + 1.0) * (-tr + 0.7) * comp.prsv_k1 + kappa0;
            let inner = kappa * (-r + 1.0) + 1.0;
            inner * inner
        }
        VdWOL1998 | RKOL1998 | PROL1998 => {
            // Saturation-coupled: value-lifted from the scalar path (see the
            // doc comment). Correct value, zero derivative slot.
            D::from(alpha(eos, tr.re(), comp))
        }
        SchmidtWenzel => {
            let s = -tr.sqrt() + 1.0;
            let m = sw_m_generic(w, tr);
            (m * s + 1.0).powi(2)
        }
        PatelTeja | PatelTejaUSB => soave(pt_f(w)),
    }
}

/// Dimensionless `(A, B, U, W)` for one component, generic in T and P (M12.3).
///
/// The temperature-generic sibling of the `(big_a, big_b, u, w)` computation
/// inside [`EosState::new`]. With `D = f64` it reproduces those fields; with a
/// dual type it carries exact ∂/∂T and/or ∂/∂P along (whichever slot the caller
/// seeded), which is what feeds the mixture T/P-derivative dual paths.
///
/// * `t` — temperature in **K**; `p` — pressure in **kPa absolute**.
///
/// Returns `(A, B, U, W)`, all **dimensionless**.
pub fn eos_dimensionless_generic<D: DualNum<f64> + Copy>(
    eos: CubicEos,
    t: D,
    p: D,
    comp: &Component,
) -> (D, D, D, D) {
    let tr = t / comp.tc;
    let pr = p / comp.pc;
    let a_val = alpha_generic(eos, tr, comp);
    if eos.is_three_parameter() {
        use CubicEos::*;
        let w = comp.omega;
        match eos {
            PatelTeja | PatelTejaUSB => {
                let big_a = a_val * (pt_om_a(w)) * pr / (tr * tr);
                let big_b = pr * pt_om_b(w) / tr;
                let big_c = pr * (1.0 - 3.0 * pt_xi_c(w)) / tr;
                (big_a, big_b, big_b + big_c, -(big_b * big_c))
            }
            SchmidtWenzel => {
                let big_a = a_val * sw_om_a(w) * pr / (tr * tr);
                let big_b = pr * sw_om_b(w) / tr;
                (
                    big_a,
                    big_b,
                    big_b * (1.0 + 3.0 * w),
                    big_b * big_b * (-3.0 * w),
                )
            }
            _ => unreachable!("is_three_parameter but not PT/SW"),
        }
    } else {
        let fc = family_constants(eos);
        let big_a = a_val * fc.om_a * pr / (tr * tr);
        let big_b = pr * fc.om_b / tr;
        (big_a, big_b, big_b * fc.k1, big_b * big_b * fc.k2)
    }
}

/// dα/dTr — analytical derivative of α with respect to reduced temperature.
///
/// Computed analytically per CLAUDE.md "Algorithm Choices": numerical
/// derivatives exist only as test oracles. The closed forms for the four
/// core variants are derived by hand from `alpha` above; see comments
/// for the chain-rule step at each branch.
pub fn d_alpha_d_tr(eos: CubicEos, tr: f64, comp: &Component) -> f64 {
    use CubicEos::*;
    let w = comp.omega;
    match eos {
        // ----- Deployable core (M7.1) -----
        // α = 1 → dα/dTr = 0
        VdW1870 => 0.0,
        // α = Tr^(−1/2) → dα/dTr = −(1/2)·Tr^(−3/2) = −1 / (2·Tr·√Tr)
        RK1949 => -0.5 / (tr * tr.sqrt()),
        RKS1972 => {
            // α = (1 + m·s)²  with s = 1 − √Tr, ds/dTr = −1/(2·√Tr).
            // dα/dTr = 2·(1 + m·s)·m·(−1/(2·√Tr)) = −m·(1 + m·s)/√Tr.
            let m = 0.48 + 1.574 * w - 0.176 * w * w;
            let s = 1.0 - tr.sqrt();
            -m * (1.0 + m * s) / tr.sqrt()
        }
        PR1976 => {
            let kappa = 0.37464 + 1.54226 * w - 0.26992 * w * w;
            let s = 1.0 - tr.sqrt();
            -kappa * (1.0 + kappa * s) / tr.sqrt()
        }
        // ----- M7.2: analytical derivatives for the remaining 2-param variants -----

        // Berthelot: α = 1/Tr → dα/dTr = −1/Tr².
        Berth1899 => -1.0 / (tr * tr),

        // VdWAda: α = 10^(m·(1 − Tr)) = exp(m·(1 − Tr)·ln10) →
        // dα/dTr = −m·ln(10)·α.
        VdWAda1984 => {
            let m = 0.228165 + 1.446486 * w - 0.648552 * w * w;
            let a = 10f64.powf(m * (1.0 - tr));
            -m * std::f64::consts::LN_10 * a
        }

        // Soave form: α = (1 + m·s)², s = 1 − √Tr, ds/dTr = −1/(2√Tr) →
        // dα/dTr = −m·(1 + m·s)/√Tr.
        RKSGD1978 => {
            let m = 0.48508 + 1.55171 * w - 0.15613 * w * w;
            let s = 1.0 - tr.sqrt();
            -m * (1.0 + m * s) / tr.sqrt()
        }
        RKSL1997 => {
            let m = 0.478972559 + 1.576809191 * w - 0.187219516 * w * w + 0.020424946 * w * w * w;
            let s = 1.0 - tr.sqrt();
            -m * (1.0 + m * s) / tr.sqrt()
        }
        RP1978 => {
            let m = 0.379642 + 1.48503 * w - 0.164423 * w * w + 0.016666 * w * w * w;
            let s = 1.0 - tr.sqrt();
            -m * (1.0 + m * s) / tr.sqrt()
        }
        PRL1997 => {
            let m = 0.378710697 + 1.487972964 * w - 0.166754831 * w * w + 0.017169486 * w * w * w;
            let s = 1.0 - tr.sqrt();
            -m * (1.0 + m * s) / tr.sqrt()
        }

        // Mathias-Naumann linear form: α = 1 + (1 − Tr)(m + n/Tr).
        // The (1−Tr)·n/Tr terms collapse exactly: dα/dTr = −m − n/Tr².
        VdWVald1989 => {
            let omegac = w * comp.zc;
            let m = 0.4745 + (2.7349 + 6.0984 * omegac) * omegac;
            let n = 0.0674 + (2.1031 + 3.9512 * omegac) * omegac;
            -m - n / (tr * tr)
        }
        RKSmn1980 => {
            let (m, n) = (comp.m_polar, comp.n_polar);
            -m - n / (tr * tr)
        }

        // ATmn exponential form: α = exp(f) with
        // f = (1 − Tr)·m·|1 − Tr|^(g−1) + n·(1/Tr − 1).
        // Using d/dTr[(1 − Tr)·|1 − Tr|^(g−1)] = −g·|1 − Tr|^(g−1):
        // dα/dTr = α·( −m·g·|1 − Tr|^(g−1) − n/Tr² ).
        RKSATmn1995 | PRATmng1997 => {
            let (m, n, g) = (comp.m_polar, comp.n_polar, comp.g_polar);
            let u = 1.0 - tr;
            let a = (u * m * u.abs().powf(g - 1.0) + n * (1.0 / tr - 1.0)).exp();
            a * (-m * g * u.abs().powf(g - 1.0) - n / (tr * tr))
        }

        // PRMmn: α = exp[(1 − Tr)·m + n·s²], s = 1 − √Tr.
        // dα/dTr = α·( −m − n·s/√Tr ).
        PRMmn1989 => {
            let (m, n) = (comp.m_polar, comp.n_polar);
            let r = tr.sqrt();
            let s = 1.0 - r;
            let a = ((1.0 - tr) * m + n * s * s).exp();
            a * (-m - n * s / r)
        }

        // PRSV: α = inner², inner = 1 + κ·(1 − √Tr),
        // κ = κ₀ + K₁·(1 + √Tr)·(0.7 − Tr). Let r = √Tr.
        // dκ/dTr = K₁·[ (0.7 − Tr)/(2r) − (1 + r) ];
        // d(inner)/dTr = dκ·(1 − r) − κ/(2r);
        // dα/dTr = 2·inner·d(inner)/dTr.
        PRSV1986 => {
            let r = tr.sqrt();
            let kappa0 = 0.378893 + 1.4897153 * w - 0.17131848 * w * w + 0.0196554 * w * w * w;
            let kappa = kappa0 + comp.prsv_k1 * (1.0 + r) * (0.7 - tr);
            let inner = 1.0 + kappa * (1.0 - r);
            let dkappa = comp.prsv_k1 * ((0.7 - tr) / (2.0 * r) - (1.0 + r));
            let dinner = dkappa * (1.0 - r) - kappa / (2.0 * r);
            2.0 * inner * dinner
        }

        // ----- OL family (M7.4, saturation-coupled α) -----
        VdWOL1998 | RKOL1998 | PROL1998 => {
            // α = Tr·(1 + S) → dα/dTr = (1 + S) + Tr·dS/dTr.
            ol_sumhk_and_deriv(eos, tr, comp)
                .map(|(s, ds)| (1.0 + s) + tr * ds)
                .unwrap_or(f64::NAN)
        }
        // ----- M7.3: 3-parameter Pascal EOS (Ref (4), TERMOII.PAS) -----
        SchmidtWenzel => {
            // α = (1 + m·s)², s = 1 − √Tr, m = m(Tr) piecewise.
            // dα/dTr = 2(1 + m·s)(m'·s + m·s'), s' = −1/(2√Tr).
            // m' is the one-sided derivative of the active branch (see
            // `sw_dm_dtr`); the slope is discontinuous at Tr=1 (documented).
            let r = tr.sqrt();
            let s = 1.0 - r;
            let m = sw_m(w, tr);
            let dm = sw_dm_dtr(w, tr);
            let ds = -0.5 / r;
            2.0 * (1.0 + m * s) * (dm * s + m * ds)
        }
        PatelTeja | PatelTejaUSB => {
            // Soave form: α = (1 + F·s)² → dα/dTr = −F(1 + F·s)/√Tr.
            let f = pt_f(w);
            let s = 1.0 - tr.sqrt();
            -f * (1.0 + f * s) / tr.sqrt()
        }
    }
}

// ===========================================================================
// EosState — the per-(EOS, T, P, component) cache (M8.2, PERFORMANCE_PROPOSAL
// §C2). Everything the Z-factor / fugacity / departure machinery needs is
// computed exactly once here: α, dα/dTr, and the dimensionless (A, B, U, W)
// groups of the generalized cubic
//
//     Z³ + (U − B − 1)·Z² + (A + W − U − B·U)·Z − (A·B + W + B·W) = 0
//
// The two-parameter Abbott families are the special case U = k1·B,
// W = k2·B² (verified coefficient-for-coefficient against the previous
// per-family cubic — see tests); the 3-parameter Pascal EOS supply their
// own (U, W). Writing every consumer against (A, B, U, W) is also the
// architecture the M8.3 mixture core builds on (Michelsen–Mollerup style,
// PERFORMANCE_PROPOSAL §B1): mixture code gets written ONCE, and each
// EOS × mixing-rule combination only supplies these four numbers.
// ===========================================================================

/// Cached dimensionless state for one (EOS, T, P, component) point.
///
/// Build it once with [`EosState::new`], then call as many property
/// methods as needed — α and dα/dTr are never recomputed. The free
/// functions ([`z_factor`], [`ln_phi_pure`], …) remain as one-shot
/// conveniences and simply construct a state internally.
#[derive(Debug, Clone, Copy)]
pub struct EosState {
    /// EOS variant this state was built for.
    pub eos: CubicEos,
    /// Temperature in **K**.
    pub t: f64,
    /// Pressure in **kPa absolute**.
    pub p: f64,
    /// Reduced temperature T/Tc. **Dimensionless.**
    pub tr: f64,
    /// α(Tr) — cached once. **Dimensionless.**
    pub alpha: f64,
    /// dα/dTr (analytical) — cached once. **Dimensionless.**
    pub d_alpha_d_tr: f64,
    /// `A = a·α·P/(R·T)²`. **Dimensionless.**
    pub big_a: f64,
    /// `B = b·P/(R·T)`. **Dimensionless.**
    pub big_b: f64,
    /// `U = u·P/(R·T)` — attractive-denominator linear coefficient
    /// (`= k1·B` for 2-parameter families). **Dimensionless.**
    pub u: f64,
    /// `W = w·(P/(R·T))²` — attractive-denominator constant coefficient
    /// (`= k2·B²` for 2-parameter families). **Dimensionless.**
    pub w: f64,
}

impl EosState {
    /// Compute the full cached state for one (EOS, T, P, component).
    ///
    /// # Arguments
    /// * `t` — Temperature in **K**.
    /// * `p` — Pressure in **kPa absolute**.
    /// * `comp` — Component data (critical constants + the parameters the
    ///   chosen α variant reads).
    pub fn new(eos: CubicEos, t: f64, p: f64, comp: &Component) -> Self {
        let tr = t / comp.tc;
        let a_val = alpha(eos, tr, comp);
        let da_val = d_alpha_d_tr(eos, tr, comp);
        let (big_a, big_b, u, w) = if eos.is_three_parameter() {
            three_param_aubw(eos, t, p, comp, a_val)
        } else {
            let fc = family_constants(eos);
            let pr = p / comp.pc;
            let big_a = fc.om_a * a_val * pr / (tr * tr);
            let big_b = fc.om_b * pr / tr;
            (big_a, big_b, fc.k1 * big_b, fc.k2 * big_b * big_b)
        };
        Self {
            eos,
            t,
            p,
            tr,
            alpha: a_val,
            d_alpha_d_tr: da_val,
            big_a,
            big_b,
            u,
            w,
        }
    }

    /// Z = P·V/(R·T) for the requested phase.
    ///
    /// Solves the generalized cubic and selects the appropriate real root:
    /// **liquid** = smallest Z above B, **vapor** = largest Z above B.
    /// Roots ≤ B are unphysical (V < b) and discarded.
    ///
    /// # Errors
    /// `EosError::Cubic` if the solver fails; `NoRootForPhase` if no real
    /// root above B exists for the requested phase.
    pub fn z(&self, phase: PhaseId) -> Result<f64, EosError> {
        let (big_a, big_b, u, w) = (self.big_a, self.big_b, self.u, self.w);
        let a2 = u - big_b - 1.0;
        let a1 = big_a + w - u - big_b * u;
        let a0 = -(big_a * big_b + w + big_b * w);
        let (roots, count) = solve_real(1.0, a2, a1, a0)?;
        select_physical_root(&roots[..count], big_b, phase)
    }

    /// Generalized attractive term `g` such that
    /// `ln φ = Z − 1 − ln(Z − B) − g`. See [`attractive_term_uw`].
    pub fn attractive_term(&self, z: f64) -> f64 {
        attractive_term_uw(z, self.big_a, self.u, self.w)
    }

    /// Pure-component ln(φ) at an already-solved Z (no second cubic solve).
    pub fn ln_phi_at(&self, z: f64) -> f64 {
        z - 1.0 - (z - self.big_b).ln() - self.attractive_term(z)
    }

    /// Pure-component ln(φ) for the requested phase.
    pub fn ln_phi(&self, phase: PhaseId) -> Result<f64, EosError> {
        Ok(self.ln_phi_at(self.z(phase)?))
    }

    /// Departure enthalpy `H^R/(R·T)` at an already-solved Z.
    pub fn h_departure_rt_at(&self, z: f64) -> f64 {
        let g = self.attractive_term(z);
        (z - 1.0) + g * (self.tr * self.d_alpha_d_tr / self.alpha - 1.0)
    }

    /// Departure enthalpy `H^R/(R·T)` for the requested phase.
    pub fn h_departure_rt(&self, phase: PhaseId) -> Result<f64, EosError> {
        Ok(self.h_departure_rt_at(self.z(phase)?))
    }

    /// Departure entropy `S^R/R` at an already-solved Z
    /// (Lewis-Randall identity `S^R/R = H^R/RT − ln φ`).
    pub fn s_departure_r_at(&self, z: f64) -> f64 {
        self.h_departure_rt_at(z) - self.ln_phi_at(z)
    }

    /// Departure entropy `S^R/R` for the requested phase.
    pub fn s_departure_r(&self, phase: PhaseId) -> Result<f64, EosError> {
        Ok(self.s_departure_r_at(self.z(phase)?))
    }
}

/// Compute Z = P·V/(R·T) for the requested phase.
///
/// One-shot convenience over [`EosState`]: builds the cached state and
/// solves. If you need more than one property at the same (T, P), build
/// an [`EosState`] once and reuse it instead.
///
/// # Arguments
/// * `t` — Temperature in **K**.
/// * `p` — Pressure in **kPa absolute**.
/// * `comp` — Component (only Tc, Pc, ω used by the core variants).
/// * `phase` — Which root to return.
///
/// # Returns
/// `Z`, dimensionless. Use `V = Z·R·T/P` to recover the molar volume.
///
/// # Errors
/// `EosError::Cubic` if the underlying solver fails. `NoRootForPhase`
/// if the cubic has no real root above B for the requested phase.
pub fn z_factor(
    eos: CubicEos,
    t: f64,
    p: f64,
    comp: &Component,
    phase: PhaseId,
) -> Result<f64, EosError> {
    EosState::new(eos, t, p, comp).z(phase)
}

/// Select the requested phase's Z root by direct comparison — no
/// filter/collect/sort (PERFORMANCE_PROPOSAL.md §C1). Roots ≤ B are
/// unphysical (V ≤ b, molecular impossibility) and skipped; among the
/// physical roots, liquid = smallest, vapor = largest. `solve_real`
/// returns roots ascending, so a single min/max scan suffices.
fn select_physical_root(roots: &[f64], big_b: f64, phase: PhaseId) -> Result<f64, EosError> {
    let mut selected: Option<f64> = None;
    for &z in roots {
        if z <= big_b {
            continue;
        }
        selected = Some(match (selected, phase) {
            (None, _) => z,
            (Some(cur), PhaseId::Liquid) => cur.min(z),
            (Some(cur), PhaseId::Vapor) => cur.max(z),
        });
    }
    selected.ok_or(EosError::NoRootForPhase { phase, big_b })
}

/// Pure-component fugacity coefficient ln(φ).
///
/// `ln(φ) = Z − 1 − ln(Z − B) − g(Z; A, U, W)`.
///
/// Returns the **natural log** of the fugacity coefficient. To get φ
/// itself, exponentiate. `ln(φ) → 0` as `P → 0` (ideal gas limit).
/// One-shot convenience over [`EosState`].
///
/// # Errors
/// Same as [`z_factor`].
pub fn ln_phi_pure(
    eos: CubicEos,
    t: f64,
    p: f64,
    comp: &Component,
    phase: PhaseId,
) -> Result<f64, EosError> {
    EosState::new(eos, t, p, comp).ln_phi(phase)
}

/// Departure enthalpy in dimensionless form `H^R / (R·T)`.
///
/// `H^R = H_real(T, P) − H_ideal_gas(T)`. Always **negative** for stable
/// liquid and vapor phases at sub-critical conditions (the attractive
/// term lowers the energy below the ideal-gas reference).
///
/// Formula (generalized cubic):
/// ```text
///   H^R/(RT) = (Z − 1) + g · (Tr·dα/dTr/α − 1)
/// ```
/// Derivation: Smith-Van Ness-Abbott §6.4, generalized to the (A, B, U, W)
/// form. Verified against PR/RKS/VdW textbook expressions.
/// One-shot convenience over [`EosState`].
///
/// # Errors
/// Same as [`z_factor`].
pub fn h_departure_rt(
    eos: CubicEos,
    t: f64,
    p: f64,
    comp: &Component,
    phase: PhaseId,
) -> Result<f64, EosError> {
    EosState::new(eos, t, p, comp).h_departure_rt(phase)
}

/// Departure entropy in dimensionless form `S^R / R`.
///
/// `S^R = S_real(T, P) − S_ideal_gas(T, P)`. Uses the identity
/// `S^R/R = H^R/(RT) − G^R/(RT)` and `G^R/(RT) = ln(φ_pure)` for a pure
/// component (Lewis-Randall). One-shot convenience over [`EosState`].
///
/// # Errors
/// Same as [`z_factor`].
pub fn s_departure_r(
    eos: CubicEos,
    t: f64,
    p: f64,
    comp: &Component,
    phase: PhaseId,
) -> Result<f64, EosError> {
    EosState::new(eos, t, p, comp).s_departure_r(phase)
}

// ===========================================================================
// Three-parameter EOS constants (M7.3) — Ref (4): Da Silva & Báez (1989),
// legacy/pascal/TERMOII.PAS. The attractive denominator is V² + uV + w'; in
// dimensionless groups U = uP/(RT), W = w'(P/RT)². The cubic and the
// fugacity/departure algebra are the SAME as the two-parameter case (which is
// the special case U = k1·B, W = k2·B²) — both now flow through `EosState`.
// The (U, W) values were verified to reproduce the legacy Patel-Teja and
// Schmidt-Wenzel cubics coefficient-for-coefficient.
// ===========================================================================

/// Dimensionless `(A, B, U, W)` for a three-parameter cubic EOS.
///
/// `a_val` is the already-computed α(Tr) — passed in by [`EosState::new`]
/// so α is evaluated exactly once per state (M8.2 cache rule).
fn three_param_aubw(
    eos: CubicEos,
    t: f64,
    p: f64,
    comp: &Component,
    a_val: f64,
) -> (f64, f64, f64, f64) {
    use CubicEos::*;
    let tr = t / comp.tc;
    let pr = p / comp.pc;
    let w = comp.omega;
    match eos {
        PatelTeja | PatelTejaUSB => {
            let big_a = pt_om_a(w) * a_val * pr / (tr * tr);
            let big_b = pt_om_b(w) * pr / tr;
            // c-parameter (dimensionless) C = cP/(RT) = (1 − 3ξc)·Pr/Tr.
            // Denominator V² + (b+c)V − bc → U = B + C, W = −B·C.
            let big_c = (1.0 - 3.0 * pt_xi_c(w)) * pr / tr;
            (big_a, big_b, big_b + big_c, -big_b * big_c)
        }
        SchmidtWenzel => {
            let big_a = sw_om_a(w) * a_val * pr / (tr * tr);
            let big_b = sw_om_b(w) * pr / tr;
            // Denominator V² + (1+3ω)bV − 3ω b² → U = (1+3ω)B, W = −3ω B².
            (
                big_a,
                big_b,
                (1.0 + 3.0 * w) * big_b,
                -3.0 * w * big_b * big_b,
            )
        }
        _ => unreachable!("three_param_aubw called for 2-parameter EOS {eos:?}"),
    }
}

/// Generalized attractive term `g(Z; A, U, W)`, defined so that
/// `ln φ = Z − 1 − ln(Z − B) − g`. Equals `(A/B)·F` of the classic
/// two-parameter form, so `ln φ` and the departures share one expression
/// across every EOS family.
///
/// Three branches on the discriminant `U² − 4W` (Müller et al. (9),
/// research paper Eqs 2.31–2.33; Pascal TERMOII.PAS `II`):
///
/// ```text
///   Δ² > 0:  g = (A/Δ)·ln[(2Z + U + Δ)/(2Z + U − Δ)],  Δ = √(U² − 4W)
///   Δ² = 0:  g = 2A/(2Z + U)              (VdW limit; U = W = 0 → A/Z)
///   Δ² < 0:  g = (2A/Δ')·(π/2 − atan[(2Z + U)/Δ']),  Δ' = √(4W − U²)
/// ```
///
/// Each branch is `A·∫_Z^∞ dZ'/(Z'² + U·Z' + W)`, so all three vanish as
/// Z → ∞ (the ideal-gas limit ln φ → 0 holds in every branch — the π/2
/// in the arctan branch is that boundary term).
///
/// The negative branch matters in practice: Schmidt-Wenzel with
/// ω ∈ (−1.94, −0.057) — e.g. hydrogen, ω ≈ −0.216 — puts `U² − 4W`
/// below zero (Ref (4), TERMOII.PAS:355 implements the same arctan
/// form). The `1e-12`-scaled tolerance keeps the degenerate branch from
/// being missed through floating-point noise.
pub fn attractive_term_uw(z: f64, big_a: f64, u: f64, w: f64) -> f64 {
    let disc = u * u - 4.0 * w;
    let scale = (u * u).max(4.0 * w.abs()).max(1e-300);
    if disc.abs() <= 1e-12 * scale {
        // Degenerate: lim_{Δ→0} (A/Δ)·ln[(x+Δ)/(x−Δ)] = 2A/x, x = 2Z+U.
        2.0 * big_a / (2.0 * z + u)
    } else if disc > 0.0 {
        let delta = disc.sqrt();
        (big_a / delta) * ((2.0 * z + u + delta) / (2.0 * z + u - delta)).ln()
    } else {
        // Complex-conjugate denominator roots: the antiderivative is an
        // arctan, and the Z' → ∞ boundary contributes the π/2 (where the
        // real branch's log contributes ln 1 = 0). Positive for A > 0,
        // matching the real branch by analytic continuation.
        let delta = (-disc).sqrt();
        (2.0 * big_a / delta) * (std::f64::consts::FRAC_PI_2 - ((2.0 * z + u) / delta).atan())
    }
}

/// Chao-Seader species selector — hydrogen and methane use distinct
/// coefficient sets from the normal-fluid correlation. Ref (4), TERMOII.PAS:386.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
#[cfg_attr(feature = "python", pyo3::pyclass(eq, eq_int))]
#[repr(i32)]
pub enum ChaoSeaderSpecies {
    /// Normal fluids (the default coefficient set).
    Normal = 0,
    /// Hydrogen-specific coefficients.
    Hydrogen = 1,
    /// Methane-specific coefficients.
    Methane = 2,
}

impl ChaoSeaderSpecies {
    /// Pick the coefficient set for a component from its name — `"hydrogen"`
    /// / `"h2"` and `"methane"` / `"ch4"` (case-insensitive) get their special
    /// sets, everything else is a normal fluid. This is what the
    /// [`LiquidModel::GraysonStreed`] K-value path uses, once per `(T, P)`, so
    /// callers never have to build a species vector by hand.
    pub fn for_component(comp: &Component) -> Self {
        let name = comp.name.trim().to_ascii_lowercase();
        match name.as_str() {
            "hydrogen" | "h2" => ChaoSeaderSpecies::Hydrogen,
            "methane" | "ch4" => ChaoSeaderSpecies::Methane,
            _ => ChaoSeaderSpecies::Normal,
        }
    }
}

/// Which published coefficient set the regular-solution liquid fugacity
/// correlation `log₁₀ ν = ν⁰ + ω·ν¹` uses (Milestone 20).
///
/// Both sets share the functional form and the ν¹ acentric correction; they
/// differ in the ν⁰ table. Grayson & Streed refitted Chao & Seader's ν⁰ to
/// extend the correlation to hydrogen-rich, high-temperature refinery service
/// (up to ~800 °F and 3000 psia), which is why every refinery simulator lists
/// "Grayson–Streed" as its own method.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
#[cfg_attr(feature = "python", pyo3::pyclass(eq, eq_int))]
#[repr(i32)]
pub enum RegularSolutionSet {
    /// Chao & Seader (1961), *AIChE J.* **7**, 598 — the original table.
    ChaoSeader1961 = 0,
    /// Grayson & Streed (1963), *6th World Pet. Congr.* Sect. VII, 233 — the
    /// refit that the legacy Pascal (Ref (4)) and every refinery package use.
    GraysonStreed1963 = 1,
}

/// ν⁰ coefficient table `A0..A9` for one `(set, species)`.
///
/// Sources, each verified against two independent transcriptions before being
/// typed here (a coefficient table is exactly the kind of place a single-source
/// typo hides): Grayson–Streed from Ref (4) `TERMOII.PAS:386` and the COMSOL
/// Liquid & Gas Properties manual Table 2-2; Chao–Seader 1961 from the FOSSEE
/// OpenModelica thermodynamics report (IIT Bombay, 2018) and DWSIM's
/// `ChaoSeader.vb`.
fn nu0_coefficients(set: RegularSolutionSet, species: ChaoSeaderSpecies) -> [f64; 10] {
    match (set, species) {
        (RegularSolutionSet::GraysonStreed1963, ChaoSeaderSpecies::Normal) => [
            2.05135, -2.10899, 0.0, -0.19396, 0.02282, 0.08852, 0.0, -0.00872, -0.00353, 0.00203,
        ],
        (RegularSolutionSet::GraysonStreed1963, ChaoSeaderSpecies::Hydrogen) => [
            1.50709, 2.74283, -0.02110, 0.00011, 0.0, 0.008585, 0.0, 0.0, 0.0, 0.0,
        ],
        (RegularSolutionSet::GraysonStreed1963, ChaoSeaderSpecies::Methane) => [
            1.36822, -1.54831, 0.0, 0.02889, -0.01076, 0.10486, -0.02529, 0.0, 0.0, 0.0,
        ],
        (RegularSolutionSet::ChaoSeader1961, ChaoSeaderSpecies::Normal) => [
            5.75748, -3.01761, -4.98500, 2.02299, 0.0, 0.08427, 0.26667, -0.31138, -0.02655,
            0.02883,
        ],
        (RegularSolutionSet::ChaoSeader1961, ChaoSeaderSpecies::Hydrogen) => [
            1.96718, 1.02972, -0.054009, 0.0005288, 0.0, 0.008585, 0.0, 0.0, 0.0, 0.0,
        ],
        (RegularSolutionSet::ChaoSeader1961, ChaoSeaderSpecies::Methane) => [
            2.43840, -2.24550, -0.34084, 0.00212, -0.00223, 0.10486, -0.03691, 0.0, 0.0, 0.0,
        ],
    }
}

/// Regular-solution pure-liquid fugacity coefficient, returned as **ln(ν)**
/// (ν = fᵢᴸ/(xᵢP) at the mixture's T and P), for either coefficient set.
///
/// The Chao–Seader framework writes `Kᵢ = νᵢ·γᵢ/φ̂ᵢⱽ`, with γᵢ from
/// Scatchard–Hildebrand regular-solution theory and φ̂ᵢⱽ from an EOS
/// (classically Redlich–Kwong). This function is the νᵢ factor only; the
/// [`crate::flash::system`] K-value path assembles the three.
///
/// ```text
///   log₁₀ ν = ν⁰ + ω·ν¹
///   ν⁰ = A0 + A1/Tr + A2·Tr + A3·Tr² + A4·Tr³ + (A5 + A6·Tr + A7·Tr²)·Pr
///        + (A8 + A9·Tr)·Pr² − log₁₀ Pr
///   ν¹ = −4.23893 + 8.65808·Tr − 1.22060/Tr − 3.15224·Tr³ − 0.025·(Pr − 0.6)
/// ```
///
/// # Arguments
/// * `set` — which ν⁰ table ([`RegularSolutionSet`]).
/// * `t` — Temperature in **K**.
/// * `p` — Pressure in **kPa absolute**.
/// * `comp` — Component (uses `tc`, `pc`, `omega`).
/// * `species` — coefficient set selector (normal / hydrogen / methane).
///
/// # Returns
/// ln νᵢ, **dimensionless**.
pub fn regular_solution_ln_nu(
    set: RegularSolutionSet,
    t: f64,
    p: f64,
    comp: &Component,
    species: ChaoSeaderSpecies,
) -> f64 {
    let tr = t / comp.tc;
    let pr = p / comp.pc;
    let a = nu0_coefficients(set, species);
    // ν¹ acentric-correction coefficients A10..A14 (shared by both sets and
    // all species).
    let q = [-4.23893, 8.65808, -1.22060, -3.15224, -0.025];
    let tr2 = tr * tr;
    let nu0 = a[0]
        + a[1] / tr
        + a[2] * tr
        + a[3] * tr2
        + a[4] * tr2 * tr
        + (a[5] + a[6] * tr + a[7] * tr2) * pr
        + (a[8] + a[9] * tr) * pr * pr
        - pr.log10();
    let nu1 = q[0] + q[1] * tr + q[2] / tr + q[3] * tr2 * tr + q[4] * (pr - 0.6);
    // log10(ν) = ν⁰ + ω·ν¹ → ln(ν) = ln(10)·log10(ν).
    (nu0 + comp.omega * nu1) * std::f64::consts::LN_10
}

/// Chao-Seader pure-liquid fugacity coefficient, returned as **ln(ν)** (the
/// natural log of the liquid fugacity coefficient ν = f_L/(xP)).
///
/// Ref (4): Da Silva & Báez (1989), legacy/pascal/TERMOII.PAS:386-405. The
/// legacy returns ν = 10^(ν⁰ + ω·ν¹); we return the natural log to match
/// [`ln_phi_pure`].
///
/// **Coefficient provenance (corrected in Milestone 20):** the table the legacy
/// Pascal carries — and therefore this function — is Grayson & Streed's 1963
/// refit, not Chao & Seader's 1961 original. The name is kept for
/// compatibility; this is exactly
/// [`regular_solution_ln_nu`]`(`[`RegularSolutionSet::GraysonStreed1963`]`, …)`.
/// For the 1961 table call that function with
/// [`RegularSolutionSet::ChaoSeader1961`].
///
/// # Arguments
/// * `t` — Temperature in **K**.
/// * `p` — Pressure in **kPa absolute**.
/// * `comp` — Component (uses `tc`, `pc`, `omega`).
/// * `species` — coefficient set selector.
///
/// # Returns
/// ln of the liquid fugacity coefficient, **dimensionless**.
pub fn chao_seader_ln_phi(t: f64, p: f64, comp: &Component, species: ChaoSeaderSpecies) -> f64 {
    regular_solution_ln_nu(RegularSolutionSet::GraysonStreed1963, t, p, comp, species)
}

// #[cfg(test)] is a conditional compilation attribute. It tells the Rust
// compiler: "only compile the code inside this module when running tests"
// (i.e., when you run `cargo test`). In a normal build (`cargo build`),
// everything inside this block is completely ignored — it won't end up
// in your final binary. This is Rust's built-in way of keeping test code
// next to the code it tests, without bloating the production build.
#[cfg(test)]
mod tests {
    // `use super::*` imports everything from the parent module (the file
    // above this block) into the test module, so we can use CubicEos
    // and its methods directly without fully qualifying them.
    use super::*;

    // #[test] marks this function as a test case. `cargo test` will
    // discover and run every function with this attribute. If the
    // function completes without panicking, the test passes. If any
    // assert! macro fails, it panics and the test is reported as failed.
    #[test]
    fn discriminant_values_match_legacy() {
        // `as i32` casts the enum variant to its integer value (set by
        // #[repr(i32)] above). These assertions verify that our Rust enum
        // values match the exact integer codes used in the legacy VB6
        // program — if someone accidentally reorders the enum or changes
        // a value, this test will catch it.
        assert_eq!(CubicEos::PR1976 as i32, 0);
        assert_eq!(CubicEos::RKS1972 as i32, 2);
        assert_eq!(CubicEos::PRSV1986 as i32, 15);
        assert_eq!(CubicEos::PROL1998 as i32, 18);
        assert_eq!(CubicEos::SchmidtWenzel as i32, 19);
        assert_eq!(CubicEos::PatelTeja as i32, 20);
        assert_eq!(CubicEos::PatelTejaUSB as i32, 21);
    }

    #[test]
    fn three_parameter_detection() {
        // assert!(condition) passes if condition is true, panics if false.
        // The `!` before the call negates the result, so these first two
        // lines check that PR1976 and RKS1972 are NOT three-parameter EOS.
        assert!(!CubicEos::PR1976.is_three_parameter());
        assert!(!CubicEos::RKS1972.is_three_parameter());
        // These three lines verify that the three-parameter EOS variants
        // are correctly identified by the is_three_parameter() method.
        assert!(CubicEos::SchmidtWenzel.is_three_parameter());
        assert!(CubicEos::PatelTeja.is_three_parameter());
        assert!(CubicEos::PatelTejaUSB.is_three_parameter());
    }

    // -----------------------------------------------------------------
    // M7.1 — tests for the deployable core (PR / RKS / RK / VdW).
    // -----------------------------------------------------------------

    fn methane() -> Component {
        // Methane critical properties + acentric factor from NIST.
        Component {
            name: "methane".into(),
            tc: 190.564,
            pc: 4599.0, // kPa
            omega: 0.0115,
            ..Component::default()
        }
    }

    fn n_pentane() -> Component {
        // n-Pentane — a more "typical" hydrocarbon with non-trivial ω.
        Component {
            name: "n-pentane".into(),
            tc: 469.7,
            pc: 3370.0,
            omega: 0.252,
            ..Component::default()
        }
    }

    #[test]
    fn family_constants_match_legacy_table() {
        // PR family — verify the high-precision OmA/OmB constants from
        // legacy/vb6/McommonFunctions.bas:273.
        let fc = family_constants(CubicEos::PR1976);
        assert_eq!(fc.k1, 2.0);
        assert_eq!(fc.k2, -1.0);
        assert!((fc.om_a - 0.457235528921382).abs() < 1e-15);
        assert!((fc.om_b - 0.0777960739038885).abs() < 1e-15);

        // RKS family.
        let fc = family_constants(CubicEos::RKS1972);
        assert_eq!(fc.k1, 1.0);
        assert_eq!(fc.k2, 0.0);
        assert!((fc.om_a - 0.427480233540341).abs() < 1e-15);
        assert!((fc.om_b - 0.0866403499649577).abs() < 1e-15);

        // VdW family — analytical fractions.
        let fc = family_constants(CubicEos::VdW1870);
        assert_eq!(fc.k1, 0.0);
        assert_eq!(fc.k2, 0.0);
        assert_eq!(fc.om_a, 27.0 / 64.0);
        assert_eq!(fc.om_b, 1.0 / 8.0);
    }

    #[test]
    fn alpha_at_tr_one_is_one_for_pr_rks() {
        // For PR and RKS, α(Tr=1) = 1 by construction of the α form.
        // VdW gives α=1 trivially; RK gives α=1/√1=1.
        let c = n_pentane();
        for eos in [
            CubicEos::PR1976,
            CubicEos::RKS1972,
            CubicEos::RK1949,
            CubicEos::VdW1870,
        ] {
            let a = alpha(eos, 1.0, &c);
            assert!(
                (a - 1.0).abs() < 1e-12,
                "{:?}: α(Tr=1) = {} (expected 1.0)",
                eos,
                a
            );
        }
    }

    /// Numerical derivative via central differences — used as test oracle.
    fn d_alpha_numerical(eos: CubicEos, tr: f64, comp: &Component, h: f64) -> f64 {
        (alpha(eos, tr + h, comp) - alpha(eos, tr - h, comp)) / (2.0 * h)
    }

    #[test]
    fn analytical_d_alpha_matches_numerical() {
        // Sweep Tr across a sub/supercritical range and verify the
        // analytical derivative agrees with central differences within
        // 1e-7 relative tolerance — a CLAUDE.md "Algorithm Choices" rule
        // (numerical derivatives are test oracles, not production code).
        let c = n_pentane();
        for eos in [
            CubicEos::PR1976,
            CubicEos::RKS1972,
            CubicEos::RK1949,
            CubicEos::VdW1870,
        ] {
            for tr in [0.5_f64, 0.8, 1.0, 1.2, 2.0] {
                let analytical = d_alpha_d_tr(eos, tr, &c);
                let numerical = d_alpha_numerical(eos, tr, &c, 1e-6);
                let rel = if analytical.abs() < 1e-10 {
                    (analytical - numerical).abs()
                } else {
                    ((analytical - numerical) / analytical).abs()
                };
                assert!(
                    rel < 1e-5,
                    "{:?} Tr={} analytical={} numerical={} rel={}",
                    eos,
                    tr,
                    analytical,
                    numerical,
                    rel
                );
            }
        }
    }

    #[test]
    fn z_factor_methane_supercritical() {
        // Methane at 300 K, 5 MPa (5000 kPa). T > Tc → only one real
        // root — vapor and "liquid" both resolve to the same root.
        let c = methane();
        let z_v = z_factor(CubicEos::PR1976, 300.0, 5000.0, &c, PhaseId::Vapor).unwrap();
        // Z should be < 1 (attractive forces dominate at moderate pressure).
        // For methane at 300 K and 5 MPa, Z ≈ 0.91 with PR.
        assert!(
            z_v > 0.8 && z_v < 1.05,
            "Z(vapor) = {} not in plausible range",
            z_v
        );
    }

    #[test]
    fn z_factor_n_pentane_two_phase() {
        // n-pentane at 400 K (Tr = 0.85), 1500 kPa. Below critical and at
        // moderate pressure → expect both liquid and vapor roots,
        // Z_liquid << Z_vapor.
        let c = n_pentane();
        let z_l = z_factor(CubicEos::PR1976, 400.0, 1500.0, &c, PhaseId::Liquid).unwrap();
        let z_v = z_factor(CubicEos::PR1976, 400.0, 1500.0, &c, PhaseId::Vapor).unwrap();
        assert!(z_l < z_v, "expected Z_liquid={} < Z_vapor={}", z_l, z_v);
        assert!(z_l < 0.1, "liquid Z should be small, got {}", z_l);
        assert!(z_v > 0.5, "vapor Z should be > 0.5, got {}", z_v);
    }

    #[test]
    fn ln_phi_ideal_gas_limit() {
        // At very low pressure, ln(φ) → 0 (any cubic EOS reduces to
        // ideal gas). Use 0.1 kPa as "very low" — should give |ln(φ)|
        // well below 1e-2 for methane at 300 K.
        let c = methane();
        for eos in [
            CubicEos::PR1976,
            CubicEos::RKS1972,
            CubicEos::RK1949,
            CubicEos::VdW1870,
        ] {
            let ln_phi = ln_phi_pure(eos, 300.0, 0.1, &c, PhaseId::Vapor).unwrap();
            assert!(
                ln_phi.abs() < 1e-3,
                "{:?}: ln(φ) at P→0 = {} (expected near 0)",
                eos,
                ln_phi
            );
        }
    }

    // -----------------------------------------------------------------
    // M7.3 — three-parameter EOS (Schmidt-Wenzel, Patel-Teja, PT-USB) and
    // Chao-Seader. Ref (4): Da Silva & Báez (1989), legacy/pascal/TERMOII.PAS.
    // -----------------------------------------------------------------

    const THREE_PARAM: [CubicEos; 3] = [
        CubicEos::SchmidtWenzel,
        CubicEos::PatelTeja,
        CubicEos::PatelTejaUSB,
    ];

    #[test]
    fn three_parameter_alpha_unity_at_critical() {
        // α(Tr=1)=1 for every 3-param EOS — the EOS-specific prefactor is
        // folded into Ω_a, keeping the same convention as the 2-param zoo.
        let c = n_pentane();
        for eos in THREE_PARAM {
            assert!((alpha(eos, 1.0, &c) - 1.0).abs() < 1e-12, "{eos:?}");
        }
    }

    #[test]
    fn three_parameter_d_alpha_matches_numerical() {
        // Analytical dα/dTr vs a central-difference oracle, sampled away from
        // the Schmidt-Wenzel Tr=1 slope kink (the ±h windows stay on one branch).
        let c = n_pentane();
        for eos in THREE_PARAM {
            for tr in [0.5_f64, 0.7, 0.9, 1.2, 1.5, 2.0] {
                let analytical = d_alpha_d_tr(eos, tr, &c);
                let numerical = d_alpha_numerical(eos, tr, &c, 1e-6);
                let rel = if analytical.abs() < 1e-10 {
                    (analytical - numerical).abs()
                } else {
                    ((analytical - numerical) / analytical).abs()
                };
                assert!(
                    rel < 1e-5,
                    "{eos:?} Tr={tr} a={analytical} n={numerical} rel={rel}"
                );
            }
        }
    }

    #[test]
    fn three_parameter_ideal_gas_limit() {
        // As P→0, Z→1 and ln φ→0 for every 3-param EOS.
        let c = n_pentane();
        for eos in THREE_PARAM {
            let z = z_factor(eos, 400.0, 1e-3, &c, PhaseId::Vapor).unwrap();
            assert!((z - 1.0).abs() < 1e-4, "{eos:?} Z={z}");
            let lnphi = ln_phi_pure(eos, 400.0, 1e-3, &c, PhaseId::Vapor).unwrap();
            assert!(lnphi.abs() < 1e-4, "{eos:?} lnphi={lnphi}");
        }
    }

    #[test]
    fn three_parameter_roots_and_fugacity_sane() {
        let c = n_pentane();
        for eos in THREE_PARAM {
            let zv = z_factor(eos, 400.0, 50.0, &c, PhaseId::Vapor).unwrap();
            assert!(zv > 0.0 && zv < 1.1, "{eos:?} Zv={zv}");
            let zl = z_factor(eos, 300.0, 2000.0, &c, PhaseId::Liquid).unwrap();
            assert!(
                zl.is_finite() && zl > 0.0 && zl < zv,
                "{eos:?} Zl={zl} Zv={zv}"
            );
            let lnphi = ln_phi_pure(eos, 400.0, 50.0, &c, PhaseId::Vapor).unwrap();
            assert!(lnphi.is_finite(), "{eos:?} lnphi={lnphi}");
        }
    }

    #[test]
    fn three_parameter_entropy_consistency() {
        // S^R/R = H^R/RT − ln φ (Lewis-Randall) and all finite.
        let c = n_pentane();
        for eos in THREE_PARAM {
            let s = s_departure_r(eos, 400.0, 500.0, &c, PhaseId::Vapor).unwrap();
            let h = h_departure_rt(eos, 400.0, 500.0, &c, PhaseId::Vapor).unwrap();
            let g = ln_phi_pure(eos, 400.0, 500.0, &c, PhaseId::Vapor).unwrap();
            assert!(s.is_finite() && (s - (h - g)).abs() < 1e-9, "{eos:?} s={s}");
        }
    }

    #[test]
    fn schmidt_wenzel_tr1_entropy_finite() {
        // Faithful + guarded: the SW dα/dTr kink at Tr=1 must NOT yield a NaN
        // entropy (the legacy returned NaN there — TERMOII.PAS:492).
        let c = n_pentane();
        let s = s_departure_r(CubicEos::SchmidtWenzel, c.tc, 500.0, &c, PhaseId::Vapor).unwrap();
        assert!(s.is_finite(), "SW entropy at Tr=1 not finite: {s}");
    }

    #[test]
    fn chao_seader_pure_fugacity_sane() {
        // Pair each coefficient set with a representative component (applying,
        // say, the H₂ set to n-pentane is physically meaningless). Each Tr=0.7.
        let hydrogen = Component {
            tc: 33.2,
            pc: 1300.0,
            omega: -0.216,
            ..Component::default()
        };
        let cases = [
            (ChaoSeaderSpecies::Normal, n_pentane()),
            (ChaoSeaderSpecies::Methane, methane()),
            (ChaoSeaderSpecies::Hydrogen, hydrogen),
        ];
        for (species, c) in cases {
            let lnphi = chao_seader_ln_phi(0.7 * c.tc, 500.0, &c, species);
            assert!(lnphi.is_finite(), "{species:?} lnphi not finite: {lnphi}");
            assert!(lnphi.abs() < 50.0, "{species:?} lnphi out of band: {lnphi}");
        }
    }

    // -----------------------------------------------------------------
    // M7.4 — OL-family α (saturation-coupled). Ref (4)/Olivera (1998).
    // -----------------------------------------------------------------

    const OL_FAMILY: [CubicEos; 3] = [CubicEos::VdWOL1998, CubicEos::RKOL1998, CubicEos::PROL1998];

    /// n-pentane with the data the OL α + saturation models need: a reduced
    /// Antoine fit (ln(P/Pc)=a1−a2/(a3+T)), a normal boiling point, and a
    /// liquid molar volume.
    fn pentane_full() -> Component {
        Component {
            name: "n-pentane".into(),
            tc: 469.7,
            pc: 3370.0,
            omega: 0.252,
            tb: 309.2,
            psat_coeffs: vec![6.738, 3165.0, 0.0],
            liquid_volume: 116.0,
            ..Component::default()
        }
    }

    #[test]
    fn ol_alpha_finite_and_positive() {
        let c = pentane_full();
        for eos in OL_FAMILY {
            for tr in [0.6_f64, 0.8, 0.95] {
                let a = alpha(eos, tr, &c);
                assert!(a.is_finite() && a > 0.0, "{eos:?} Tr={tr} α={a}");
            }
        }
    }

    #[test]
    fn ol_d_alpha_matches_numerical() {
        // With the Antoine sat model the OL dα/dTr is fully analytical
        // (analytical dPsat/dT through the chain rule) — match the oracle.
        let c = pentane_full();
        for eos in OL_FAMILY {
            for tr in [0.6_f64, 0.75, 0.9] {
                let analytical = d_alpha_d_tr(eos, tr, &c);
                let numerical = d_alpha_numerical(eos, tr, &c, 1e-6);
                let rel = ((analytical - numerical) / analytical).abs();
                assert!(
                    rel < 1e-4,
                    "{eos:?} Tr={tr} a={analytical} n={numerical} rel={rel}"
                );
            }
        }
    }

    #[test]
    fn ol_z_factor_and_entropy_finite() {
        let c = pentane_full();
        for eos in OL_FAMILY {
            let zv = z_factor(eos, 400.0, 100.0, &c, PhaseId::Vapor).unwrap();
            assert!(zv.is_finite() && zv > 0.0 && zv < 1.1, "{eos:?} Zv={zv}");
            let s = s_departure_r(eos, 400.0, 100.0, &c, PhaseId::Vapor).unwrap();
            assert!(s.is_finite(), "{eos:?} S^R/R={s}");
        }
    }

    #[test]
    fn ol_alpha_nan_without_sat_data() {
        // No Antoine coeffs → reduced_psat fails → α is NaN (downstream
        // z_factor then errors cleanly rather than returning a bogus root).
        let c = n_pentane();
        assert!(alpha(CubicEos::RKOL1998, 0.8, &c).is_nan());
    }

    // -----------------------------------------------------------------
    // M7.2 — tests for the remaining twelve two-parameter α variants.
    // -----------------------------------------------------------------

    /// A synthetic "polar" component carrying every fitted parameter the
    /// M7.2 variants might read (zc, m/n/g, prsv_k1). Critical constants
    /// loosely mimic water so the numbers stay in a plausible range; the
    /// exact values don't matter — the oracle test only checks internal
    /// consistency of α and its analytical derivative.
    fn polar_component() -> Component {
        Component {
            name: "synthetic-polar".into(),
            tc: 647.1,
            pc: 22064.0,
            omega: 0.344,
            zc: 0.229,
            // Fitted polar parameters. g > 1 keeps |1−Tr|^(g−1) smooth at
            // Tr = 1 for the ATmn exponential variants.
            m_polar: 0.45,
            n_polar: 0.12,
            g_polar: 1.5,
            prsv_k1: 0.07,
            ..Component::default()
        }
    }

    /// Every M7.2 variant, plus the synthetic component each one is
    /// exercised against. The OL family is excluded — it is M7.4.
    const M72_VARIANTS: [CubicEos; 12] = [
        CubicEos::Berth1899,
        CubicEos::VdWAda1984,
        CubicEos::RKSGD1978,
        CubicEos::RKSL1997,
        CubicEos::RP1978,
        CubicEos::PRL1997,
        CubicEos::VdWVald1989,
        CubicEos::RKSmn1980,
        CubicEos::RKSATmn1995,
        CubicEos::PRATmng1997,
        CubicEos::PRMmn1989,
        CubicEos::PRSV1986,
    ];

    #[test]
    fn m72_alpha_is_finite_and_positive() {
        // α must be a finite, strictly positive number everywhere — a
        // negative or NaN α would make a·α negative and break the EOS.
        let c = polar_component();
        for eos in M72_VARIANTS {
            for tr in [0.5_f64, 0.7, 0.9, 1.0, 1.3, 2.0] {
                let a = alpha(eos, tr, &c);
                assert!(
                    a.is_finite() && a > 0.0,
                    "{:?} Tr={} gave α={} (expected finite, positive)",
                    eos,
                    tr,
                    a
                );
            }
        }
    }

    #[test]
    fn m72_alpha_at_tr_one_is_one() {
        // By construction, all of these α forms collapse to 1 at Tr = 1
        // (the 1 − √Tr / 1 − Tr / 1/Tr − 1 factors all vanish there). This
        // is the single most important sanity check on the porting: the
        // attractive term must reduce to its critical-point value at Tc.
        let c = polar_component();
        for eos in M72_VARIANTS {
            let a = alpha(eos, 1.0, &c);
            assert!(
                (a - 1.0).abs() < 1e-12,
                "{:?}: α(Tr=1) = {} (expected 1.0)",
                eos,
                a
            );
        }
    }

    #[test]
    fn m72_analytical_d_alpha_matches_numerical() {
        // The CLAUDE.md "Algorithm Choices" rule: analytical dα/dTr is the
        // production path, the central-difference value is only an oracle.
        // Tr = 1.0 is skipped for the abs-power ATmn variants — even with
        // g = 1.5 the |1−Tr|^(g−1) factor is only C¹ there, so a symmetric
        // difference would clip the corner; every other Tr is smooth.
        let c = polar_component();
        for eos in M72_VARIANTS {
            for tr in [0.55_f64, 0.7, 0.85, 1.1, 1.4] {
                let analytical = d_alpha_d_tr(eos, tr, &c);
                let numerical = d_alpha_numerical(eos, tr, &c, 1e-6);
                let rel = if analytical.abs() < 1e-8 {
                    (analytical - numerical).abs()
                } else {
                    ((analytical - numerical) / analytical).abs()
                };
                assert!(
                    rel < 1e-5,
                    "{:?} Tr={} analytical={} numerical={} rel={}",
                    eos,
                    tr,
                    analytical,
                    numerical,
                    rel
                );
            }
        }
    }

    #[test]
    fn prsv_k1_recovers_kappa0_when_zero() {
        // With K₁ = 0 the PRSV κ collapses to κ₀(ω), so α(PRSV) must match
        // the bracketed [1 + κ₀(1−√Tr)]² form exactly. This pins the K₁
        // contribution as a pure additive correction.
        let mut c = polar_component();
        c.prsv_k1 = 0.0;
        let w = c.omega;
        let kappa0 = 0.378893 + 1.4897153 * w - 0.17131848 * w * w + 0.0196554 * w * w * w;
        for tr in [0.6_f64, 0.8, 1.2] {
            let s = 1.0 - tr.sqrt();
            let expected = (1.0 + kappa0 * s).powi(2);
            let got = alpha(CubicEos::PRSV1986, tr, &c);
            assert!(
                (got - expected).abs() < 1e-12,
                "PRSV K₁=0 Tr={}: got {} expected {}",
                tr,
                got,
                expected
            );
        }
    }

    // -----------------------------------------------------------------
    // M8.2 — EosState cache + generalized attractive term.
    // -----------------------------------------------------------------

    #[test]
    fn eos_state_reuse_matches_one_shot_functions() {
        // One EosState reused across all four properties must reproduce
        // the one-shot free functions exactly, for a 2-param and a
        // 3-param EOS (both flow through the same generalized path now).
        let c = n_pentane();
        for eos in [CubicEos::PR1976, CubicEos::RKS1972, CubicEos::PatelTeja] {
            let st = EosState::new(eos, 400.0, 1500.0, &c);
            for phase in [PhaseId::Vapor, PhaseId::Liquid] {
                assert_eq!(
                    st.z(phase).unwrap(),
                    z_factor(eos, 400.0, 1500.0, &c, phase).unwrap(),
                    "{eos:?} {phase:?} Z"
                );
                assert_eq!(
                    st.ln_phi(phase).unwrap(),
                    ln_phi_pure(eos, 400.0, 1500.0, &c, phase).unwrap(),
                    "{eos:?} {phase:?} ln_phi"
                );
                assert_eq!(
                    st.h_departure_rt(phase).unwrap(),
                    h_departure_rt(eos, 400.0, 1500.0, &c, phase).unwrap(),
                    "{eos:?} {phase:?} H^R"
                );
                assert_eq!(
                    st.s_departure_r(phase).unwrap(),
                    s_departure_r(eos, 400.0, 1500.0, &c, phase).unwrap(),
                    "{eos:?} {phase:?} S^R"
                );
            }
        }
    }

    #[test]
    fn generalized_uw_reproduces_two_param_families() {
        // U = k1·B, W = k2·B² must give the same attractive term as the
        // old per-family (A/B)·F closed form. Check PR (log branch) and
        // VdW (degenerate branch) hand-built references.
        let (big_a, big_b, z) = (0.5_f64, 0.05_f64, 0.8_f64);
        // PR: k1=2, k2=−1 → disc = 8B², Δ = 2√2·B.
        let sd = 8.0_f64.sqrt();
        let f_pr = (1.0 / sd)
            * ((2.0 * z + big_b * (2.0 + sd)) / (2.0 * z + big_b * (2.0 - sd))).ln()
            * (big_a / big_b);
        let g_pr = attractive_term_uw(z, big_a, 2.0 * big_b, -big_b * big_b);
        assert!((f_pr - g_pr).abs() < 1e-14, "PR: {f_pr} vs {g_pr}");
        // VdW: U = W = 0 → g = A/Z.
        let g_vdw = attractive_term_uw(z, big_a, 0.0, 0.0);
        assert!((g_vdw - big_a / z).abs() < 1e-14, "VdW: {g_vdw}");
    }

    #[test]
    fn attractive_term_arctan_branch_hydrogen_like_sw() {
        // Schmidt-Wenzel with hydrogen's ω ≈ −0.216 puts U² − 4W < 0
        // (the 9ω² + 18ω + 1 quadratic is negative there), exercising the
        // arctan branch. Ref (4), TERMOII.PAS:355. It must be finite and
        // → 0 as Z → ∞ (ideal-gas limit — the π/2 boundary term).
        let hydrogen = Component {
            name: "hydrogen".into(),
            tc: 33.2,
            pc: 1300.0,
            omega: -0.216,
            ..Component::default()
        };
        let st = EosState::new(CubicEos::SchmidtWenzel, 40.0, 500.0, &hydrogen);
        assert!(
            st.u * st.u - 4.0 * st.w < 0.0,
            "test premise: SW hydrogen must hit the negative-discriminant branch"
        );
        let g = st.attractive_term(0.9);
        assert!(g.is_finite() && g > 0.0, "g = {g}");
        // Ideal-gas limit: g(Z→∞) → 0.
        assert!(st.attractive_term(1e9).abs() < 1e-6);
        // And the full ln φ machinery stays finite through this branch.
        let lnphi = st.ln_phi(PhaseId::Vapor).unwrap();
        assert!(lnphi.is_finite(), "ln φ = {lnphi}");
    }

    #[test]
    fn m72_z_factor_and_ln_phi_work() {
        // The full Z-factor / ln(φ) machinery must run end-to-end for the
        // new variants (they were panicking before M7.2). Use the polar
        // component below its critical point.
        let c = polar_component();
        for eos in M72_VARIANTS {
            let z = z_factor(eos, 500.0, 2000.0, &c, PhaseId::Vapor).unwrap();
            assert!(z > 0.0 && z < 1.2, "{:?}: Z={} out of range", eos, z);
            let ln_phi = ln_phi_pure(eos, 500.0, 2000.0, &c, PhaseId::Vapor).unwrap();
            assert!(ln_phi.is_finite(), "{:?}: ln(φ)={} not finite", eos, ln_phi);
        }
    }

    // ── M12.3: the generic α / dimensionless-param path must match scalar ──

    /// Every non-OL EOS variant, paired with a component that carries the
    /// parameters it reads. OL is excluded (value-lifted, tested separately).
    fn non_ol_variants_with_comp() -> Vec<(CubicEos, Component)> {
        let mut v: Vec<(CubicEos, Component)> = vec![
            (CubicEos::VdW1870, n_pentane()),
            (CubicEos::RK1949, n_pentane()),
            (CubicEos::RKS1972, n_pentane()),
            (CubicEos::PR1976, n_pentane()),
            (CubicEos::SchmidtWenzel, n_pentane()),
            (CubicEos::PatelTeja, n_pentane()),
            (CubicEos::PatelTejaUSB, n_pentane()),
        ];
        for eos in M72_VARIANTS {
            v.push((eos, polar_component()));
        }
        v
    }

    #[test]
    fn alpha_generic_f64_matches_scalar_alpha() {
        // alpha_generic::<f64> must reproduce the scalar `alpha` across every
        // variant and a Tr sweep spanning sub- and supercritical (the
        // Schmidt-Wenzel piecewise breakpoint at Tr = 1 included). VdWAda1984
        // is re-expressed exp(y·ln10) vs 10^y, so allow a tiny rel tolerance.
        for (eos, c) in non_ol_variants_with_comp() {
            for &tr in &[0.5, 0.7, 0.95, 1.0, 1.05, 1.5, 2.0] {
                let scalar = alpha(eos, tr, &c);
                let generic = alpha_generic::<f64>(eos, tr, &c);
                let tol = 1e-12 * scalar.abs().max(1.0);
                assert!(
                    (scalar - generic).abs() <= tol,
                    "{eos:?} Tr={tr}: scalar={scalar} generic={generic}"
                );
            }
        }
    }

    #[test]
    fn alpha_generic_dual_matches_analytic_d_alpha() {
        // The whole point of the generic path: a first-order dual seeded on Tr
        // must reproduce the hand-analytic d_alpha_d_tr. This is the α-side
        // guarantee underpinning the T-derivative dual sweeps.
        use num_dual::Dual64;
        for (eos, c) in non_ol_variants_with_comp() {
            for &tr in &[0.6, 0.85, 1.2, 1.8] {
                let d = alpha_generic(eos, Dual64::new(tr, 1.0), &c);
                let analytic = d_alpha_d_tr(eos, tr, &c);
                let tol = 1e-9 * analytic.abs().max(1.0);
                assert!(
                    (d.eps - analytic).abs() <= tol,
                    "{eos:?} Tr={tr}: dual dα/dTr={} analytic={analytic}",
                    d.eps
                );
            }
        }
    }

    #[test]
    fn eos_dimensionless_generic_f64_matches_eos_state() {
        // (A, B, U, W) from the generic path must equal the EosState fields
        // for every variant, so the mixture T/P-derivative path is anchored to
        // the same numbers the scalar engine uses.
        let (t, p) = (360.0, 2500.0);
        for (eos, c) in non_ol_variants_with_comp() {
            let st = EosState::new(eos, t, p, &c);
            let (a, b, u, w) = eos_dimensionless_generic::<f64>(eos, t, p, &c);
            let tol = 1e-10;
            assert!(
                (a - st.big_a).abs() <= tol * st.big_a.abs().max(1.0),
                "{eos:?} A"
            );
            assert!(
                (b - st.big_b).abs() <= tol * st.big_b.abs().max(1.0),
                "{eos:?} B"
            );
            assert!((u - st.u).abs() <= tol * st.u.abs().max(1.0), "{eos:?} U");
            assert!((w - st.w).abs() <= tol * st.w.abs().max(1.0), "{eos:?} W");
        }
    }

    #[test]
    fn grayson_streed_set_is_the_legacy_chao_seader_table_and_1961_differs() {
        // M20: `chao_seader_ln_phi` carries the Grayson-Streed 1963 table.
        let c = n_pentane();
        for species in [
            ChaoSeaderSpecies::Normal,
            ChaoSeaderSpecies::Hydrogen,
            ChaoSeaderSpecies::Methane,
        ] {
            let legacy = chao_seader_ln_phi(400.0, 1500.0, &c, species);
            let gs = regular_solution_ln_nu(
                RegularSolutionSet::GraysonStreed1963,
                400.0,
                1500.0,
                &c,
                species,
            );
            let cs = regular_solution_ln_nu(
                RegularSolutionSet::ChaoSeader1961,
                400.0,
                1500.0,
                &c,
                species,
            );
            assert_eq!(legacy, gs);
            assert!((cs - gs).abs() > 1e-3, "{species:?}: {cs} vs {gs}");
            // For a normal fluid at moderate conditions the two tables agree to
            // within a few tenths in log10 ν — they were fitted to the same
            // data. (The H2/CH4 sets applied to pentane's Tc/Pc are not
            // physically meaningful and are only checked for being distinct.)
            if species == ChaoSeaderSpecies::Normal {
                assert!(
                    (cs - gs).abs() < 0.5 * std::f64::consts::LN_10,
                    "{species:?}: {cs} vs {gs}"
                );
            }
        }
    }

    #[test]
    fn species_selection_by_name() {
        let mk = |n: &str| Component {
            name: n.into(),
            ..Component::default()
        };
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("Hydrogen")),
            ChaoSeaderSpecies::Hydrogen
        );
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("H2")),
            ChaoSeaderSpecies::Hydrogen
        );
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("methane")),
            ChaoSeaderSpecies::Methane
        );
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("CH4")),
            ChaoSeaderSpecies::Methane
        );
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("n-decane")),
            ChaoSeaderSpecies::Normal
        );
        assert_eq!(
            ChaoSeaderSpecies::for_component(&mk("PC-12")),
            ChaoSeaderSpecies::Normal
        );
    }
}