vle-thermo 0.16.0

//! Python bindings for the VLE engine.
//!
//! Builds a Python extension module `vle._engine` when this crate is
//! compiled with the `python` feature. End-users import it as
//! `from vle import _engine` (or, more commonly, via the higher-level
//! `vle` Python wrapper that re-exports the pieces they need).
//!
//! ## What lives here in Milestone 5
//!
//! This is the first PyO3 module the crate ships. The Milestone 5 goal is
//! to prove the boundary end-to-end (wheel builds, installs, imports,
//! tests pass on every platform in CI). To that end this module exposes
//! a small but useful surface:
//!
//! - [`version()`] — a free function returning the crate's semver string.
//!   Lets calling Python code introspect which version of the engine is
//!   loaded.
//! - Four enum classes — [`CubicEos`], [`ActivityModel`], [`MixingRule`],
//!   [`SatPressureModel`]. These are the model-selection types Python
//!   callers will pass to every subsequent calculation in M6+. Exposing
//!   them now means every later milestone can wire new functions in
//!   without first having to set up enum bindings.
//!
//! ## Adding more bindings in M6+
//!
//! Per CLAUDE.md's "PyO3 Bindings Rule (M5+)", every milestone after M5
//! that adds public Rust functionality must also expose it via PyO3 in
//! the same commit series. The convention:
//!
//! - **Free functions**: `#[pyfunction]` in this file or a co-located
//!   module file (e.g., `engine/src/py_bindings/numerics.rs`); register
//!   in the `#[pymodule]` block below via `wrap_pyfunction!`.
//! - **Types**: add `#[cfg_attr(feature = "python", pyo3::pyclass(...))]`
//!   to the Rust struct/enum definition where it already lives, then
//!   register here via `m.add_class::<T>()`.
//! - **Tests**: every new binding gets at least one round-trip test in
//!   `python/tests/test_engine.py` (or a sibling file).
//!
//! The CI matrix exercises every wheel via `pytest`, so a missing binding
//! is a hard failure — not a code-review oversight.

use std::cell::RefCell;

use pyo3::exceptions::{PyRuntimeError, PyValueError};
use pyo3::prelude::*;

use crate::activity::{
    ActivityModel, excess_enthalpy as excess_enthalpy_rs, excess_entropy as excess_entropy_rs,
    excess_gibbs as excess_gibbs_rs, ln_gamma as ln_gamma_rs,
};
use crate::eos::{
    ChaoSeaderSpecies, CubicEos, EosError, PhaseId, alpha as eos_alpha_rs,
    chao_seader_ln_phi as chao_seader_ln_phi_rs, d_alpha_d_tr as eos_d_alpha_rs, family_constants,
    h_departure_rt, ln_phi_pure, s_departure_r, z_factor,
};
use crate::liquid_volume::{VolumeModel, liquid_molar_volume as liquid_molar_volume_rs};
use crate::saturation::{
    SatError, SatPressureModel, boiling_temperature as boiling_temperature_rs,
    d_psat_dt as d_psat_dt_rs, d_psat_dt_antoine, poynting_factor as poynting_factor_rs,
    psat as psat_rs, psat_antoine, psat_maxwell as psat_maxwell_rs,
    reduced_psat as reduced_psat_rs,
};
use crate::types::Component;
use crate::virial::{
    b_mix as virial_b_mix, h_departure_rt_virial, ln_phi_mix_virial, ln_phi_pure_virial, pitzer_b,
    pitzer_b0, pitzer_b1, pitzer_d_b_d_t, s_departure_r_virial, z_factor_virial,
};

/// Return the engine crate's version string (matches `Cargo.toml`).
///
/// Useful for diagnostics: `print(vle._engine.version())` confirms which
/// wheel is loaded when troubleshooting installs.
#[pyfunction]
fn version() -> &'static str {
    env!("CARGO_PKG_VERSION")
}

/// Return the raw TOML text of the canonical unit catalog.
///
/// Same bytes that seed `UnitRegistry::with_vle_defaults()` on the Rust
/// side. The Python `vle.units` module parses this to register gauge units
/// (barg/psig/kPag), `lbmol`, and `BTU/lbmol` with `pint`, so the two
/// halves of the stack share a single source of truth for unit constants.
#[pyfunction]
fn default_units_toml() -> &'static str {
    vle_units::default_units_toml()
}

// ──────────────────────────────────────────────────────────────────────
// M6.1 — numerics bindings.
//
// All four solvers + the small utility surface are exposed under flat
// names (no `numerics.` prefix on the Python side) so user code reads
// like `from vle._engine import brent, solve_cubic, halley`. The
// higher-level `vle` Python wrapper can re-export these into more
// pythonic submodules later (`vle.numerics.brent` etc.) without
// changing the binding shape.
//
// The three callback-taking functions (brent, illinois, halley) follow
// the same pattern: cache any Python error raised inside the callback
// in a RefCell, return NaN to the Rust solver (which handles NaN
// gracefully via its `NanEvaluation` / `NonFiniteEvaluation` error
// variants), then re-raise the original Python error after the solver
// returns. This preserves Python tracebacks across the FFI boundary
// instead of collapsing them into a generic Rust error message.
// ──────────────────────────────────────────────────────────────────────

/// Solve `a·x³ + b·x² + c·x + d = 0` over the reals.
///
/// Returns a list of every real root in ascending order — 1, 2, or 3
/// entries depending on the discriminant. Raises `ValueError` if `a`
/// is near zero (not actually a cubic) or any coefficient is non-finite.
#[pyfunction]
fn solve_cubic(a: f64, b: f64, c: f64, d: f64) -> PyResult<Vec<f64>> {
    // The Rust core returns a stack array + live count (allocation-free hot
    // path, M8.2); the FFI boundary converts to a Python list, so the
    // Python-visible API is unchanged.
    crate::numerics::cubic::solve_real(a, b, c, d)
        .map(|(roots, count)| roots[..count].to_vec())
        .map_err(|e| PyValueError::new_err(e.to_string()))
}

/// Brent's method scalar root finder on a bracketed interval.
///
/// `f(a)` and `f(b)` must have opposite signs. Returns the root in
/// `[a, b]`. Raises `RuntimeError` if the bracket is invalid, the
/// iteration limit is hit, or `f` returns NaN. Any exception raised
/// inside `f` is re-raised through the Rust solver.
#[pyfunction]
#[pyo3(signature = (f, a, b, xtol = 1e-9, max_iter = 100))]
fn brent(py: Python<'_>, f: PyObject, a: f64, b: f64, xtol: f64, max_iter: usize) -> PyResult<f64> {
    let err_cache: RefCell<Option<PyErr>> = RefCell::new(None);
    let result = crate::numerics::root_finding::brent(
        |x| call_scalar_callback(py, &f, x, &err_cache),
        a,
        b,
        xtol,
        max_iter,
    );
    if let Some(e) = err_cache.into_inner() {
        return Err(e);
    }
    result.map_err(|e| PyRuntimeError::new_err(e.to_string()))
}

/// Illinois algorithm (modified Regula Falsi) scalar root finder.
///
/// Lighter-weight alternative to [`brent`]: predictable per-iteration
/// cost, slightly slower convergence. Same bracketing requirement and
/// error semantics.
#[pyfunction]
#[pyo3(signature = (f, a, b, xtol = 1e-9, max_iter = 100))]
fn illinois(
    py: Python<'_>,
    f: PyObject,
    a: f64,
    b: f64,
    xtol: f64,
    max_iter: usize,
) -> PyResult<f64> {
    let err_cache: RefCell<Option<PyErr>> = RefCell::new(None);
    let result = crate::numerics::root_finding::illinois(
        |x| call_scalar_callback(py, &f, x, &err_cache),
        a,
        b,
        xtol,
        max_iter,
    );
    if let Some(e) = err_cache.into_inner() {
        return Err(e);
    }
    result.map_err(|e| PyRuntimeError::new_err(e.to_string()))
}

/// Halley's method scalar root finder (cubic convergence).
///
/// `f_and_derivs(x)` must return a 3-tuple `(f, f', f'')`. Not
/// bracketed — caller is responsible for a good initial guess `x0`.
/// Raises `RuntimeError` on non-convergence, singular step, or
/// non-finite evaluation.
#[pyfunction]
#[pyo3(signature = (f_and_derivs, x0, xtol = 1e-12, max_iter = 50))]
fn halley(
    py: Python<'_>,
    f_and_derivs: PyObject,
    x0: f64,
    xtol: f64,
    max_iter: usize,
) -> PyResult<f64> {
    let err_cache: RefCell<Option<PyErr>> = RefCell::new(None);
    let result = crate::numerics::halley::halley(
        |x| match f_and_derivs
            .call1(py, (x,))
            .and_then(|r| r.extract::<(f64, f64, f64)>(py))
        {
            Ok(triple) => triple,
            Err(e) => {
                if err_cache.borrow().is_none() {
                    *err_cache.borrow_mut() = Some(e);
                }
                (f64::NAN, f64::NAN, f64::NAN)
            }
        },
        x0,
        xtol,
        max_iter,
    );
    if let Some(e) = err_cache.into_inner() {
        return Err(e);
    }
    result.map_err(|e| PyRuntimeError::new_err(e.to_string()))
}

/// Broyden's "good" quasi-Newton method for `F(x) = 0` with `F: ℝⁿ → ℝⁿ`.
///
/// `f(x)` must accept a sequence of `n` floats and return a sequence of
/// `n` floats. `x0` sets the problem size. Tolerances and refresh
/// cadence default to engineering values (xtol=ftol=1e-8, max_iter=100,
/// refresh_every=5, fd_step=1e-7); tighten via the keyword args when
/// you need it.
///
/// Raises `ValueError` for input-shape problems and `RuntimeError` for
/// convergence/singular-Jacobian failures. Python exceptions raised
/// inside `f` re-raise verbatim through the Rust solver.
#[pyfunction]
#[pyo3(signature = (f, x0, xtol = 1e-8, ftol = 1e-8, max_iter = 100, refresh_every = 5, fd_step = 1e-7))]
#[allow(clippy::too_many_arguments)]
fn broyden(
    py: Python<'_>,
    f: PyObject,
    x0: Vec<f64>,
    xtol: f64,
    ftol: f64,
    max_iter: usize,
    refresh_every: usize,
    fd_step: f64,
) -> PyResult<Vec<f64>> {
    let cfg = crate::numerics::broyden::BroydenConfig {
        xtol,
        ftol,
        max_iter,
        refresh_every,
        fd_step,
    };
    // Cache the first Python exception so the Rust solver can keep
    // running on a "NaN cycle" of empty residuals, then we re-raise
    // the original traceback after the solver returns. Matches the
    // pattern used by brent / illinois / halley above.
    let err_cache: RefCell<Option<PyErr>> = RefCell::new(None);
    let result = crate::numerics::broyden::broyden(
        |x| match f
            .call1(py, (x.to_vec(),))
            .and_then(|r| r.extract::<Vec<f64>>(py))
        {
            Ok(v) => v,
            Err(e) => {
                if err_cache.borrow().is_none() {
                    *err_cache.borrow_mut() = Some(e);
                }
                // Return a NaN vector of the right length so the Rust
                // solver hits NonFiniteEvaluation on the very next
                // check and exits cleanly.
                vec![f64::NAN; x.len()]
            }
        },
        &x0,
        cfg,
    );
    if let Some(e) = err_cache.into_inner() {
        return Err(e);
    }
    result.map_err(|e| match e {
        crate::numerics::broyden::BroydenError::DimensionMismatch { .. } => {
            PyValueError::new_err(e.to_string())
        }
        _ => PyRuntimeError::new_err(e.to_string()),
    })
}

/// Helper for the scalar-callback solvers: call into Python, extract
/// f64, cache the first error and return NaN so the Rust solver fails
/// fast with `NanEvaluation`.
fn call_scalar_callback(
    py: Python<'_>,
    f: &PyObject,
    x: f64,
    err_cache: &RefCell<Option<PyErr>>,
) -> f64 {
    // Short-circuit subsequent calls once an error is pending — saves
    // wasted Python-call work while we wait for the Rust loop to exit.
    if err_cache.borrow().is_some() {
        return f64::NAN;
    }
    match f.call1(py, (x,)).and_then(|r| r.extract::<f64>(py)) {
        Ok(v) => v,
        Err(e) => {
            *err_cache.borrow_mut() = Some(e);
            f64::NAN
        }
    }
}

// ---- Utility functions ----------------------------------------------

/// Compute `1 - sum(xs)`. Used as the Rachford-Rice / bubble-point
/// composition-closure residual.
#[pyfunction]
fn sum_frac_residual(xs: Vec<f64>) -> f64 {
    crate::numerics::utils::sum_frac_residual(&xs)
}

/// L1 (sum-of-absolutes) norm of a vector.
#[pyfunction]
fn norm_l1(xs: Vec<f64>) -> f64 {
    crate::numerics::utils::norm_l1(&xs)
}

/// L2 (Euclidean) norm of a vector.
#[pyfunction]
fn norm_l2(xs: Vec<f64>) -> f64 {
    crate::numerics::utils::norm_l2(&xs)
}

/// L∞ (max-of-absolutes) norm of a vector.
#[pyfunction]
fn norm_linf(xs: Vec<f64>) -> f64 {
    crate::numerics::utils::norm_linf(&xs)
}

// ──────────────────────────────────────────────────────────────────────
// M7.1 — pure-component EOS bindings.
//
// Cubic EOS expressions take scalar floats for the component data
// (tc, pc, omega) rather than a `Component` pyclass — same convention
// as the M6.1 numerics bindings. The higher-level `vle` Python wrapper
// can pull these into a class-shaped API once we have a richer property
// database; in M7 keeping the C ABI scalar-only keeps the FFI surface
// minimal and the test matrix tractable.
// ──────────────────────────────────────────────────────────────────────

/// Build a minimal Component for the cubic-EOS layer. Internal helper —
/// not exposed to Python. Only the four fields used by the core PR/RKS/
/// RK/VdW path are populated; everything else stays at its `Default`.
fn comp_for_eos(tc: f64, pc: f64, omega: f64) -> Component {
    Component {
        tc,
        pc,
        omega,
        ..Component::default()
    }
}

/// Build a Component carrying the parameters the M7.2 α functions read.
/// Internal helper for the `_ex` α bindings. Tc/Pc are left at 1.0 because
/// α(Tr) is a pure function of the *reduced* temperature and these polar
/// parameters — it never touches the dimensional critical constants.
fn comp_for_alpha(omega: f64, zc: f64, m: f64, n: f64, g: f64, prsv_k1: f64) -> Component {
    Component {
        tc: 1.0,
        pc: 1.0,
        omega,
        zc,
        m_polar: m,
        n_polar: n,
        g_polar: g,
        prsv_k1,
        ..Component::default()
    }
}

fn phase_from_str(s: &str) -> PyResult<PhaseId> {
    match s.to_ascii_lowercase().as_str() {
        "vapor" | "v" | "gas" => Ok(PhaseId::Vapor),
        "liquid" | "l" => Ok(PhaseId::Liquid),
        other => Err(PyValueError::new_err(format!(
            "phase must be 'vapor' or 'liquid' (got {other:?})"
        ))),
    }
}

/// Map an `EosError` to a Python exception. `NotImplemented` becomes a
/// `NotImplementedError`; the rest become `RuntimeError` so the
/// originating message survives the FFI hop.
fn map_eos_err(e: EosError) -> PyErr {
    match e {
        EosError::NotImplemented(_) => {
            pyo3::exceptions::PyNotImplementedError::new_err(e.to_string())
        }
        _ => PyRuntimeError::new_err(e.to_string()),
    }
}

/// α(Tr) for the requested cubic EOS.
///
/// This is the acentric-factor-only entry point: it covers every variant
/// whose α depends solely on ω (PR1976, RKS1972, RK1949, VdW1870,
/// Berth1899, VdWAda1984, RKSGD1978, RKSL1997, RP1978, PRL1997). For the
/// variants that additionally read `Zc` or a fitted polar parameter
/// (VdWVald1989, RKSmn1980, RKSATmn1995, PRATmng1997, PRMmn1989, PRSV1986)
/// use [`eos_alpha_ex`], which accepts them.
///
/// Panics (surfaced as a Python `PanicException`) for the OL family
/// (M7.4 — saturation-coupled) and the 3-param Pascal EOS (M7.3).
#[pyfunction]
fn eos_alpha(eos: CubicEos, tr: f64, omega: f64) -> f64 {
    eos_alpha_rs(eos, tr, &comp_for_eos(1.0, 1.0, omega))
}

/// dα/dTr for the requested cubic EOS (analytical). ω-only entry point;
/// see [`eos_d_alpha_d_tr_ex`] for the parameterized variants.
#[pyfunction]
fn eos_d_alpha_d_tr(eos: CubicEos, tr: f64, omega: f64) -> f64 {
    eos_d_alpha_rs(eos, tr, &comp_for_eos(1.0, 1.0, omega))
}

/// α(Tr) for the requested cubic EOS, with the full parameter set.
///
/// Extends [`eos_alpha`] with the fields the M7.2 polar/fitted variants
/// read: `zc` (VdWVald1989), `m`/`n` (RKSmn1980, PRMmn1989), `m`/`n`/`g`
/// (RKSATmn1995, PRATmng1997), and `prsv_k1` (PRSV1986). The ω-only
/// variants ignore the extra arguments, so this is a strict superset of
/// `eos_alpha` — passing the defaults reproduces it exactly.
///
/// All extra parameters default to 0, so callers only specify what their
/// chosen variant needs.
#[pyfunction]
#[pyo3(signature = (eos, tr, omega, zc=0.0, m=0.0, n=0.0, g=0.0, prsv_k1=0.0))]
#[allow(clippy::too_many_arguments)]
fn eos_alpha_ex(
    eos: CubicEos,
    tr: f64,
    omega: f64,
    zc: f64,
    m: f64,
    n: f64,
    g: f64,
    prsv_k1: f64,
) -> f64 {
    eos_alpha_rs(eos, tr, &comp_for_alpha(omega, zc, m, n, g, prsv_k1))
}

/// dα/dTr for the requested cubic EOS (analytical), with the full
/// parameter set. The `_ex` counterpart of [`eos_d_alpha_d_tr`] — see
/// [`eos_alpha_ex`] for the parameter meanings.
#[pyfunction]
#[pyo3(signature = (eos, tr, omega, zc=0.0, m=0.0, n=0.0, g=0.0, prsv_k1=0.0))]
#[allow(clippy::too_many_arguments)]
fn eos_d_alpha_d_tr_ex(
    eos: CubicEos,
    tr: f64,
    omega: f64,
    zc: f64,
    m: f64,
    n: f64,
    g: f64,
    prsv_k1: f64,
) -> f64 {
    eos_d_alpha_rs(eos, tr, &comp_for_alpha(omega, zc, m, n, g, prsv_k1))
}

/// Family constants (k1, k2, OmA, OmB) for the EOS, as a 4-tuple.
#[pyfunction]
fn eos_family_constants(eos: CubicEos) -> (f64, f64, f64, f64) {
    let fc = family_constants(eos);
    (fc.k1, fc.k2, fc.om_a, fc.om_b)
}

/// Z = P·V/(R·T) for the requested phase.
///
/// `t` in K, `p` in kPa absolute. `phase` is "vapor" or "liquid".
#[pyfunction]
fn eos_z_factor(
    eos: CubicEos,
    t: f64,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    phase: &str,
) -> PyResult<f64> {
    let comp = comp_for_eos(tc, pc, omega);
    let phase = phase_from_str(phase)?;
    z_factor(eos, t, p, &comp, phase).map_err(map_eos_err)
}

/// Pure-component ln(φ) for the requested phase.
#[pyfunction]
fn eos_ln_phi_pure(
    eos: CubicEos,
    t: f64,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    phase: &str,
) -> PyResult<f64> {
    let comp = comp_for_eos(tc, pc, omega);
    let phase = phase_from_str(phase)?;
    ln_phi_pure(eos, t, p, &comp, phase).map_err(map_eos_err)
}

/// Departure enthalpy H^R/(R·T), dimensionless.
#[pyfunction]
fn eos_h_departure_rt(
    eos: CubicEos,
    t: f64,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    phase: &str,
) -> PyResult<f64> {
    let comp = comp_for_eos(tc, pc, omega);
    let phase = phase_from_str(phase)?;
    h_departure_rt(eos, t, p, &comp, phase).map_err(map_eos_err)
}

/// Departure entropy S^R/R, dimensionless.
#[pyfunction]
fn eos_s_departure_r(
    eos: CubicEos,
    t: f64,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    phase: &str,
) -> PyResult<f64> {
    let comp = comp_for_eos(tc, pc, omega);
    let phase = phase_from_str(phase)?;
    s_departure_r(eos, t, p, &comp, phase).map_err(map_eos_err)
}

/// Chao-Seader pure-liquid fugacity coefficient as **ln(ν)** (M7.3).
///
/// `species` selects the coefficient set (`ChaoSeaderSpecies.Normal`,
/// `.Hydrogen`, or `.Methane`). T in **K**, P in **kPa absolute**.
/// Ref (4): Da Silva & Báez (1989), legacy/pascal/TERMOII.PAS.
#[pyfunction]
fn chao_seader_ln_phi(
    t: f64,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    species: ChaoSeaderSpecies,
) -> f64 {
    let comp = comp_for_eos(tc, pc, omega);
    chao_seader_ln_phi_rs(t, p, &comp, species)
}

// ──────────────────────────────────────────────────────────────────────
// M7.5 — Antoine saturation bindings.
// ──────────────────────────────────────────────────────────────────────

fn map_sat_err(e: SatError) -> PyErr {
    match e {
        SatError::NotImplemented(_) => {
            pyo3::exceptions::PyNotImplementedError::new_err(e.to_string())
        }
        SatError::BadCoefficients { .. } => PyValueError::new_err(e.to_string()),
        SatError::OutOfRange(_) => PyValueError::new_err(e.to_string()),
        SatError::Maxwell(_) => PyRuntimeError::new_err(e.to_string()),
    }
}

/// Build a Component carrying the fields the saturation layer reads.
fn comp_for_sat(tc: f64, pc: f64, omega: f64, tb: f64, coeffs: Vec<f64>) -> Component {
    Component {
        tc,
        pc,
        omega,
        tb,
        psat_coeffs: coeffs,
        ..Component::default()
    }
}

/// Antoine saturation pressure `ln(Psat/Pc) = a1 − a2/(a3 + T)`.
///
/// Returns Psat in **kPa**. `coeffs` is `[a1, a2, a3]`.
#[pyfunction]
fn antoine_psat(t: f64, pc: f64, coeffs: Vec<f64>) -> PyResult<f64> {
    let comp = Component {
        pc,
        psat_coeffs: coeffs,
        ..Component::default()
    };
    psat_antoine(&comp, t).map_err(map_sat_err)
}

/// Analytical dPsat/dT for the Antoine form. Returns kPa/K.
#[pyfunction]
fn antoine_d_psat_dt(t: f64, pc: f64, coeffs: Vec<f64>) -> PyResult<f64> {
    let comp = Component {
        pc,
        psat_coeffs: coeffs,
        ..Component::default()
    };
    d_psat_dt_antoine(&comp, t).map_err(map_sat_err)
}

// ──────────────────────────────────────────────────────────────────────
// M7.4 — advanced saturation models + OL-family α.
// ──────────────────────────────────────────────────────────────────────

/// Saturation pressure for any non-Maxwell model. Returns Psat in **kPa**.
/// `coeffs` are the Antoine `[a1,a2,a3]` (or DIPPR `[c0..c4]` for `Polynomial`);
/// the corresponding-states models (Riedel/Müller/RPM) read `tc`, `pc`, `omega`, `tb`.
#[pyfunction]
#[pyo3(signature = (model, t, tc, pc, omega=0.0, tb=0.0, coeffs=vec![]))]
fn sat_psat(
    model: SatPressureModel,
    t: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    tb: f64,
    coeffs: Vec<f64>,
) -> PyResult<f64> {
    psat_rs(model, &comp_for_sat(tc, pc, omega, tb, coeffs), t).map_err(map_sat_err)
}

/// `dPsat/dT` (kPa/K) for any non-Maxwell model — analytical for Antoine,
/// central-difference otherwise.
#[pyfunction]
#[pyo3(signature = (model, t, tc, pc, omega=0.0, tb=0.0, coeffs=vec![]))]
fn sat_d_psat_dt(
    model: SatPressureModel,
    t: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    tb: f64,
    coeffs: Vec<f64>,
) -> PyResult<f64> {
    d_psat_dt_rs(model, &comp_for_sat(tc, pc, omega, tb, coeffs), t).map_err(map_sat_err)
}

/// Reduced saturation pressure `Psat/Pc` (dimensionless) for a non-Maxwell model.
#[pyfunction]
#[pyo3(signature = (model, t, tc, pc, omega=0.0, tb=0.0, coeffs=vec![]))]
fn sat_reduced_psat(
    model: SatPressureModel,
    t: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    tb: f64,
    coeffs: Vec<f64>,
) -> PyResult<f64> {
    reduced_psat_rs(model, &comp_for_sat(tc, pc, omega, tb, coeffs), t).map_err(map_sat_err)
}

/// Maxwell equal-area saturation pressure for a cubic EOS at temperature `t`.
/// Returns Psat in **kPa**. `coeffs` seed the initial guess (Antoine if length 3).
#[pyfunction]
#[pyo3(signature = (eos, t, tc, pc, omega=0.0, coeffs=vec![]))]
fn sat_maxwell(
    eos: CubicEos,
    t: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    coeffs: Vec<f64>,
) -> PyResult<f64> {
    psat_maxwell_rs(eos, &comp_for_sat(tc, pc, omega, 0.0, coeffs), t).map_err(map_sat_err)
}

/// Boiling temperature (**K**): invert `Psat(T)=P` for the given model.
#[pyfunction]
#[pyo3(signature = (model, p, tc, pc, omega=0.0, tb=0.0, coeffs=vec![]))]
fn boiling_temperature(
    model: SatPressureModel,
    p: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    tb: f64,
    coeffs: Vec<f64>,
) -> PyResult<f64> {
    boiling_temperature_rs(model, &comp_for_sat(tc, pc, omega, tb, coeffs), p).map_err(map_sat_err)
}

/// Poynting factor `exp[V_L·(P − Psat)/(R·T)]` (dimensionless). `liquid_volume`
/// in **cm³/mol**; `p`, `psat` in **kPa**; `t` in **K**.
#[pyfunction]
fn poynting_factor(p: f64, psat: f64, t: f64, liquid_volume: f64) -> f64 {
    let comp = Component {
        liquid_volume,
        ..Component::default()
    };
    poynting_factor_rs(&comp, p, psat, t)
}

/// OL-family α(Tr) — needs the component's saturation data. `sat_model` selects
/// the reduced-saturation-pressure source; pass Antoine `coeffs=[a1,a2,a3]` (or
/// `tb` for the corresponding-states models).
#[pyfunction]
#[pyo3(signature = (eos, tr, tc, pc, omega, sat_model, tb=0.0, coeffs=vec![]))]
#[allow(clippy::too_many_arguments)]
fn eos_alpha_ol(
    eos: CubicEos,
    tr: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    sat_model: SatPressureModel,
    tb: f64,
    coeffs: Vec<f64>,
) -> f64 {
    let mut comp = comp_for_sat(tc, pc, omega, tb, coeffs);
    comp.sat_model = sat_model;
    eos_alpha_rs(eos, tr, &comp)
}

/// Analytical dα/dTr for the OL family (companion to [`eos_alpha_ol`]).
#[pyfunction]
#[pyo3(signature = (eos, tr, tc, pc, omega, sat_model, tb=0.0, coeffs=vec![]))]
#[allow(clippy::too_many_arguments)]
fn eos_d_alpha_d_tr_ol(
    eos: CubicEos,
    tr: f64,
    tc: f64,
    pc: f64,
    omega: f64,
    sat_model: SatPressureModel,
    tb: f64,
    coeffs: Vec<f64>,
) -> f64 {
    let mut comp = comp_for_sat(tc, pc, omega, tb, coeffs);
    comp.sat_model = sat_model;
    eos_d_alpha_rs(eos, tr, &comp)
}

// ──────────────────────────────────────────────────────────────────────
// M7.6 — Virial bindings.
// ──────────────────────────────────────────────────────────────────────

/// Pitzer B⁰(Tr) — simple-fluid reduced second virial coefficient.
#[pyfunction]
fn virial_pitzer_b0(tr: f64) -> f64 {
    pitzer_b0(tr)
}

/// Pitzer B¹(Tr) — acentric correction to the reduced second virial coefficient.
#[pyfunction]
fn virial_pitzer_b1(tr: f64) -> f64 {
    pitzer_b1(tr)
}

/// Second virial coefficient B(T) for a pure component. Returns cm³/mol.
#[pyfunction]
fn virial_b_pure(tc: f64, pc: f64, omega: f64, t: f64) -> f64 {
    pitzer_b(&comp_for_eos(tc, pc, omega), t)
}

/// dB/dT for a pure component. Returns cm³/(mol·K).
#[pyfunction]
fn virial_d_b_d_t_pure(tc: f64, pc: f64, omega: f64, t: f64) -> f64 {
    pitzer_d_b_d_t(&comp_for_eos(tc, pc, omega), t)
}

/// Z = 1 + B·P/(R·T) — pure-component truncated virial.
#[pyfunction]
fn virial_z(tc: f64, pc: f64, omega: f64, t: f64, p: f64) -> f64 {
    z_factor_virial(&comp_for_eos(tc, pc, omega), t, p)
}

/// ln(φ) from the truncated virial — pure.
#[pyfunction]
fn virial_ln_phi(tc: f64, pc: f64, omega: f64, t: f64, p: f64) -> f64 {
    ln_phi_pure_virial(&comp_for_eos(tc, pc, omega), t, p)
}

/// Departure enthalpy H^R/(R·T) from the truncated virial — pure.
#[pyfunction]
fn virial_h_dep_rt(tc: f64, pc: f64, omega: f64, t: f64, p: f64) -> f64 {
    h_departure_rt_virial(&comp_for_eos(tc, pc, omega), t, p)
}

/// Departure entropy S^R/R from the truncated virial — pure.
#[pyfunction]
fn virial_s_dep_r(tc: f64, pc: f64, omega: f64, t: f64, p: f64) -> f64 {
    s_departure_r_virial(&comp_for_eos(tc, pc, omega), t, p)
}

/// Mixture B(T, x) — quadratic Lewis-Randall mixing with Pitzer cross-terms.
///
/// Parallel input vectors `tcs`, `pcs`, `omegas`, `mole_fractions` all
/// must have the same length. Returns B_mix in cm³/mol.
#[pyfunction]
fn virial_b_mix_py(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    mole_fractions: Vec<f64>,
    t: f64,
) -> PyResult<f64> {
    let n = tcs.len();
    if pcs.len() != n || omegas.len() != n || mole_fractions.len() != n {
        return Err(PyValueError::new_err(
            "tcs, pcs, omegas, mole_fractions must all have the same length",
        ));
    }
    let comps: Vec<Component> = (0..n)
        .map(|i| comp_for_eos(tcs[i], pcs[i], omegas[i]))
        .collect();
    virial_b_mix(&comps, &mole_fractions, t).map_err(|e| PyRuntimeError::new_err(e.to_string()))
}

/// Mixture partial fugacity coefficients ln(φᵢ) from the truncated virial.
///
/// Returns one ln(φᵢ) per component in input order.
#[pyfunction]
fn virial_ln_phi_mix(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    mole_fractions: Vec<f64>,
    t: f64,
    p: f64,
) -> PyResult<Vec<f64>> {
    let n = tcs.len();
    if pcs.len() != n || omegas.len() != n || mole_fractions.len() != n {
        return Err(PyValueError::new_err(
            "tcs, pcs, omegas, mole_fractions must all have the same length",
        ));
    }
    let comps: Vec<Component> = (0..n)
        .map(|i| comp_for_eos(tcs[i], pcs[i], omegas[i]))
        .collect();
    ln_phi_mix_virial(&comps, &mole_fractions, t, p)
        .map_err(|e| PyRuntimeError::new_err(e.to_string()))
}

// ──────────────────────────────────────────────────────────────────────
// M8.1 — liquid molar volume + activity-coefficient models.
// ──────────────────────────────────────────────────────────────────────

/// Saturated-liquid molar volume Vᴸ for one component.
///
/// `model` is [`VolumeModel::Rackett`] (uses `zra`) or
/// [`VolumeModel::Thomson`] (uses `vstar` = V* in cm³/mol and `omega_srk`).
/// `tc` in K, `pc` in kPa, `t` in K. Returns cm³/mol.
#[pyfunction]
#[pyo3(signature = (model, tc, pc, t, zra=0.0, vstar=0.0, omega_srk=0.0))]
fn liquid_molar_volume(
    model: VolumeModel,
    tc: f64,
    pc: f64,
    t: f64,
    zra: f64,
    vstar: f64,
    omega_srk: f64,
) -> f64 {
    let comp = Component {
        tc,
        pc,
        zra,
        liquid_volume: vstar,
        omega_srk,
        ..Component::default()
    };
    liquid_molar_volume_rs(model, &comp, t)
}

/// ln(γᵢ) for component `i`. See the Rust [`crate::activity::ln_gamma`] for the
/// per-model `aij` convention; `alpha` is the N×N symmetric NRTL non-randomness
/// matrix (empty for every other model); `vl` (cm³/mol) and `delta`
/// ((cal/cm³)^0.5) may be empty lists for models that don't use them. `t` in K.
#[pyfunction]
#[pyo3(signature = (model, i, x, aij, alpha=vec![], vl=vec![], delta=vec![], t=298.15))]
#[allow(clippy::too_many_arguments)]
fn activity_ln_gamma(
    model: ActivityModel,
    i: usize,
    x: Vec<f64>,
    aij: Vec<Vec<f64>>,
    alpha: Vec<Vec<f64>>,
    vl: Vec<f64>,
    delta: Vec<f64>,
    t: f64,
) -> f64 {
    ln_gamma_rs(model, i, &x, &aij, &alpha, &vl, &delta, t)
}

/// Excess Gibbs energy Gᴱ = RT Σ xᵢ ln γᵢ, in kJ/kmol. Args as in
/// [`activity_ln_gamma`] (minus the component index).
#[pyfunction]
#[pyo3(signature = (model, x, aij, alpha=vec![], vl=vec![], delta=vec![], t=298.15))]
fn activity_excess_gibbs(
    model: ActivityModel,
    x: Vec<f64>,
    aij: Vec<Vec<f64>>,
    alpha: Vec<Vec<f64>>,
    vl: Vec<f64>,
    delta: Vec<f64>,
    t: f64,
) -> f64 {
    excess_gibbs_rs(model, &x, &aij, &alpha, &vl, &delta, t)
}

/// Excess enthalpy Hᴱ (analytical), in kJ/kmol. Args as in [`activity_excess_gibbs`].
#[pyfunction]
#[pyo3(signature = (model, x, aij, alpha=vec![], vl=vec![], delta=vec![], t=298.15))]
fn activity_excess_enthalpy(
    model: ActivityModel,
    x: Vec<f64>,
    aij: Vec<Vec<f64>>,
    alpha: Vec<Vec<f64>>,
    vl: Vec<f64>,
    delta: Vec<f64>,
    t: f64,
) -> f64 {
    excess_enthalpy_rs(model, &x, &aij, &alpha, &vl, &delta, t)
}

/// Excess entropy Sᴱ = (Hᴱ − Gᴱ)/T, in kJ/(kmol·K). Args as in
/// [`activity_excess_gibbs`].
#[pyfunction]
#[pyo3(signature = (model, x, aij, alpha=vec![], vl=vec![], delta=vec![], t=298.15))]
fn activity_excess_entropy(
    model: ActivityModel,
    x: Vec<f64>,
    aij: Vec<Vec<f64>>,
    alpha: Vec<Vec<f64>>,
    vl: Vec<f64>,
    delta: Vec<f64>,
    t: f64,
) -> f64 {
    excess_entropy_rs(model, &x, &aij, &alpha, &vl, &delta, t)
}

// ──────────────────────────────────────────────────────────────────────
// M8.3 / M8.4 — mixture core: mixing rules, multicomponent fugacity, the
// exact composition-derivative Jacobian, and mixture energy properties.
//
// The scalar-array FFI convention (tcs/pcs/omegas + parallel lists)
// continues from the multicomponent virial bindings. The mixing rule and
// (optional) activity model are passed as the already-exposed enum
// classes. GE-based rules (Wong-Sandler, Huron-Vidal, MHV1/2) require the
// `ge_model` plus its `ge_aij`/`ge_vl`/`ge_delta` arrays; classical rules
// leave `ge_model=None`.
// ──────────────────────────────────────────────────────────────────────

use crate::mixture::{
    GeSpec, MixError, MixtureSpec, chao_seader_ln_phi_mix as chao_seader_ln_phi_mix_rs,
    d_ln_phi_d_n as d_ln_phi_d_n_rs, ln_phi_mix as ln_phi_mix_rs, z_mix as z_mix_rs,
};

/// Build the component vector for a mixture binding from parallel scalar
/// arrays. `cp` rows (may be empty) carry the ideal-Cp/R coefficients used
/// only by the energy bindings; everything else stays at `Default`.
fn mix_components(
    tcs: &[f64],
    pcs: &[f64],
    omegas: &[f64],
    cp: &[Vec<f64>],
) -> PyResult<Vec<Component>> {
    let n = tcs.len();
    if pcs.len() != n || omegas.len() != n {
        return Err(PyValueError::new_err(
            "tcs, pcs, omegas must have the same length",
        ));
    }
    (0..n)
        .map(|i| {
            let mut c = comp_for_eos(tcs[i], pcs[i], omegas[i]);
            if let Some(row) = cp.get(i) {
                if row.len() == 5 {
                    c.cp_coeffs = [row[0], row[1], row[2], row[3], row[4]];
                } else if !row.is_empty() {
                    return Err(PyValueError::new_err(
                        "each cp_coeffs row must have exactly 5 entries",
                    ));
                }
            }
            Ok(c)
        })
        .collect()
}

/// Map a [`MixError`] to a Python exception (ValueError for input/pairing
/// problems, RuntimeError for numerical failures).
fn map_mix_err(e: MixError) -> PyErr {
    match e {
        MixError::Dimension(_) | MixError::Unsupported(_) => PyValueError::new_err(e.to_string()),
        _ => PyRuntimeError::new_err(e.to_string()),
    }
}

/// Assemble a [`GeSpec`] from the optional activity model + its arrays.
/// Returns `None` when `ge_model` is `None` (classical rules).
fn ge_spec<'a>(
    ge_model: Option<ActivityModel>,
    ge_aij: &'a [Vec<f64>],
    ge_vl: &'a [f64],
    ge_delta: &'a [f64],
) -> Option<GeSpec<'a>> {
    ge_model.map(|model| GeSpec {
        model,
        aij: ge_aij,
        alpha: &[],
        vl: ge_vl,
        delta: ge_delta,
    })
}

/// Mixture compressibility factor Z for the requested phase.
///
/// `eos`/`rule` are the model-selection enums; `tcs`/`pcs`/`omegas`/`x` are
/// parallel per-component arrays; `kij` is the N×N interaction matrix (pass
/// `[]` for all-zero). GE-based rules also need `ge_model` + `ge_aij`/
/// `ge_vl`/`ge_delta`. `t` in K, `p` in kPa abs, `phase` = "vapor"/"liquid".
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, x, kij, t, p, phase,
    ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_z(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<f64> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    z_mix_rs(&spec, t, p, &x, phase).map_err(map_mix_err)
}

/// Partial fugacity coefficients ln φ̂ᵢ for every component (a list).
/// Arguments as in [`mixture_z`].
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, x, kij, t, p, phase,
    ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_ln_phi(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<Vec<f64>> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    ln_phi_mix_rs(&spec, t, p, &x, phase).map_err(map_mix_err)
}

/// Exact composition Jacobian ∂ln φ̂ᵢ/∂nⱼ, returned as an N×N list of lists.
/// Analytic for classical mixing + 2-parameter EOS; dual-number AD
/// otherwise (never finite differences). Arguments as in [`mixture_z`].
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, x, kij, t, p, phase,
    ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_d_ln_phi_d_n(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<Vec<Vec<f64>>> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    d_ln_phi_d_n_rs(&spec, t, p, &x, phase).map_err(map_mix_err)
}

/// Mixture residual enthalpy H^R/(R·T) (dimensionless) for the phase.
/// Arguments as in [`mixture_z`].
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, x, kij, t, p, phase,
    ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_h_departure_rt(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<f64> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    crate::energy::h_departure_rt_mix(&spec, t, p, &x, phase).map_err(map_mix_err)
}

/// Mixture residual entropy S^R/R (dimensionless) for the phase.
/// Arguments as in [`mixture_z`].
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, x, kij, t, p, phase,
    ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_s_departure_r(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<f64> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    crate::energy::s_departure_r_mix(&spec, t, p, &x, phase).map_err(map_mix_err)
}

/// Total molar enthalpy and entropy of one phase (ideal + residual),
/// returned as an `(H, S)` tuple — H in kJ/kmol, S in kJ/(kmol·K).
///
/// `cp_coeffs` is an N×5 list of ideal-Cp/R polynomial coefficients (one
/// row per component). `t_ref`/`p_ref` set the ideal-gas reference; the
/// per-component Hᵢ°/Sᵢ° default to 0. Other arguments as in [`mixture_z`].
#[pyfunction]
#[pyo3(signature = (eos, rule, tcs, pcs, omegas, cp_coeffs, x, kij, t, p, phase,
    t_ref=298.15, p_ref=101.325, ge_model=None, ge_aij=vec![], ge_vl=vec![], ge_delta=vec![]))]
#[allow(clippy::too_many_arguments)]
fn mixture_phase_enthalpy_entropy(
    eos: CubicEos,
    rule: crate::mixing::MixingRule,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    cp_coeffs: Vec<Vec<f64>>,
    x: Vec<f64>,
    kij: Vec<Vec<f64>>,
    t: f64,
    p: f64,
    phase: &str,
    t_ref: f64,
    p_ref: f64,
    ge_model: Option<ActivityModel>,
    ge_aij: Vec<Vec<f64>>,
    ge_vl: Vec<f64>,
    ge_delta: Vec<f64>,
) -> PyResult<(f64, f64)> {
    let comps = mix_components(&tcs, &pcs, &omegas, &cp_coeffs)?;
    let phase = phase_from_str(phase)?;
    let ge = ge_spec(ge_model, &ge_aij, &ge_vl, &ge_delta);
    let spec = MixtureSpec {
        eos,
        rule,
        components: &comps,
        kij: &kij,
        ge,
    };
    crate::energy::phase_enthalpy_entropy(&spec, t, p, &x, phase, t_ref, p_ref, &[], &[])
        .map_err(map_mix_err)
}

/// Ideal-gas mixture enthalpy relative to `t_ref`, in kJ/kmol.
/// `cp_coeffs` is N×5 (ideal-Cp/R coefficients per component).
#[pyfunction]
#[pyo3(signature = (tcs, pcs, omegas, cp_coeffs, x, t, t_ref=298.15))]
fn mixture_ideal_enthalpy(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    cp_coeffs: Vec<Vec<f64>>,
    x: Vec<f64>,
    t: f64,
    t_ref: f64,
) -> PyResult<f64> {
    let comps = mix_components(&tcs, &pcs, &omegas, &cp_coeffs)?;
    Ok(crate::energy::ideal_enthalpy_mix(&comps, &x, t, t_ref, &[]))
}

/// Ideal-gas mixture entropy relative to `(t_ref, p_ref)` including the
/// `−R Σ xᵢ ln xᵢ` mixing term, in kJ/(kmol·K). `cp_coeffs` is N×5.
#[pyfunction]
#[pyo3(signature = (tcs, pcs, omegas, cp_coeffs, x, t, p, t_ref=298.15, p_ref=101.325))]
#[allow(clippy::too_many_arguments)]
fn mixture_ideal_entropy(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    cp_coeffs: Vec<Vec<f64>>,
    x: Vec<f64>,
    t: f64,
    p: f64,
    t_ref: f64,
    p_ref: f64,
) -> PyResult<f64> {
    let comps = mix_components(&tcs, &pcs, &omegas, &cp_coeffs)?;
    Ok(crate::energy::ideal_entropy_mix(
        &comps,
        &x,
        t,
        p,
        t_ref,
        p_ref,
        &[],
    ))
}

/// Chao-Seader multicomponent liquid ln νᵢ (a list). `species` is a list of
/// [`ChaoSeaderSpecies`] (one per component). `t` in K, `p` in kPa abs.
#[pyfunction]
fn mixture_chao_seader_ln_phi(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    species: Vec<ChaoSeaderSpecies>,
    t: f64,
    p: f64,
) -> PyResult<Vec<f64>> {
    let comps = mix_components(&tcs, &pcs, &omegas, &[])?;
    chao_seader_ln_phi_mix_rs(&comps, &species, t, p).map_err(map_mix_err)
}

// ──────────────────────────────────────────────────────────────────────
// M9 — flash, bubble/dew, and stability bindings.
//
// The `SystemSpec` combines a vapor and a liquid model; Python selects each
// with a small string tag plus an optional enum:
//   vapor_kind ∈ {"ideal", "virial", "cubic"}  (+ vapor_eos for "cubic")
//   liquid_kind ∈ {"ideal", "cubic", "activity", "chao_seader"}
//       (+ liquid_eos for "cubic", + liquid_activity for "activity")
// Component data is the usual parallel arrays, extended with per-component
// `psat_coeffs` (reduced-Antoine, for the γ-φ path) and `vl` (liquid molar
// volume, for Poynting / Wilson / Scatchard). Results come back as tuples.
// ──────────────────────────────────────────────────────────────────────

use crate::eos::{LiquidModel, VaporModel};
use crate::flash::adiabatic::flash_adiabatic;
use crate::flash::aij_regression::{AijBubblePoint, fit_aij};
use crate::flash::bubble::{bubble_pressure, bubble_temperature};
use crate::flash::critical::critical_point;
use crate::flash::dew::{dew_pressure, dew_temperature};
use crate::flash::envelope::trace_envelope;
use crate::flash::isothermal::{flash_isothermal, rachford_rice as rachford_rice_rs};
use crate::flash::kij_regression::{BubblePoint, fit_kij};
use crate::flash::stability::{Stability, stability_analysis};
use crate::flash::{FlashError, SystemSpec, k_values as flash_k_values};
use crate::mixing::MixingRule;

/// Map a [`FlashError`] to a Python exception. Input-condition problems
/// (bad shapes, unsupported model combos, no two-phase root) become
/// `ValueError`; numerical failures become `RuntimeError`.
fn map_flash_err(e: FlashError) -> PyErr {
    match e {
        // Bad-argument errors surface as Python `ValueError`; genuine runtime
        // failures (non-convergence, thermodynamic breakdown) as `RuntimeError`.
        FlashError::Dimension(_)
        | FlashError::InvalidInput(_)
        | FlashError::Unsupported(_)
        | FlashError::NoRachfordRiceRoot { .. } => PyValueError::new_err(e.to_string()),
        _ => PyRuntimeError::new_err(e.to_string()),
    }
}

/// Build the flash `Component` vector including the saturation coefficients
/// and liquid molar volumes the γ-φ path needs.
fn flash_components(
    tcs: &[f64],
    pcs: &[f64],
    omegas: &[f64],
    psat_coeffs: &[Vec<f64>],
    vl: &[f64],
) -> PyResult<Vec<Component>> {
    let n = tcs.len();
    if pcs.len() != n || omegas.len() != n {
        return Err(PyValueError::new_err(
            "tcs, pcs, omegas must have the same length",
        ));
    }
    Ok((0..n)
        .map(|i| Component {
            tc: tcs[i],
            pc: pcs[i],
            omega: omegas[i],
            psat_coeffs: psat_coeffs.get(i).cloned().unwrap_or_default(),
            liquid_volume: vl.get(i).copied().unwrap_or(0.0),
            ..Component::default()
        })
        .collect())
}

/// Resolve the vapor model from a string tag + optional EOS.
fn vapor_model(kind: &str, eos: Option<CubicEos>) -> PyResult<VaporModel> {
    match kind.to_ascii_lowercase().as_str() {
        "ideal" | "idealgas" | "ideal_gas" => Ok(VaporModel::IdealGas),
        "virial" => Ok(VaporModel::Virial),
        "cubic" | "eos" => eos
            .map(VaporModel::Cubic)
            .ok_or_else(|| PyValueError::new_err("vapor_kind='cubic' needs vapor_eos")),
        other => Err(PyValueError::new_err(format!(
            "vapor_kind must be 'ideal', 'virial', or 'cubic' (got {other:?})"
        ))),
    }
}

/// Resolve the liquid model from a string tag + optional EOS / activity model.
fn liquid_model(
    kind: &str,
    eos: Option<CubicEos>,
    activity: Option<ActivityModel>,
) -> PyResult<LiquidModel> {
    match kind.to_ascii_lowercase().as_str() {
        "ideal" | "idealsolution" | "ideal_solution" => Ok(LiquidModel::IdealSolution),
        "cubic" | "eos" => eos
            .map(LiquidModel::Cubic)
            .ok_or_else(|| PyValueError::new_err("liquid_kind='cubic' needs liquid_eos")),
        "activity" | "gamma" => activity
            .map(LiquidModel::Activity)
            .ok_or_else(|| PyValueError::new_err("liquid_kind='activity' needs liquid_activity")),
        "chao_seader" | "chaoseader" => Ok(LiquidModel::ChaoSeader),
        "grayson_streed" | "graysonstreed" | "gs" => Ok(LiquidModel::GraysonStreed),
        "bk10" | "braun_k10" | "braunk10" => Ok(LiquidModel::BraunK10),
        other => Err(PyValueError::new_err(format!(
            "liquid_kind must be 'ideal', 'cubic', 'activity', 'chao_seader', \
             'grayson_streed', or 'bk10' (got {other:?})"
        ))),
    }
}

/// Solve the scalar Rachford-Rice equation for the vapor fraction β at fixed
/// K-values (§F). Raises `ValueError` if there is no interior root (the
/// mixture cannot be two-phase at these K), `RuntimeError` on non-convergence.
#[pyfunction]
#[pyo3(signature = (z, k, tol=1e-12, max_iter=200))]
fn rachford_rice(z: Vec<f64>, k: Vec<f64>, tol: f64, max_iter: usize) -> PyResult<f64> {
    rachford_rice_rs(&z, &k, tol, max_iter).map_err(map_flash_err)
}

/// A bundle of the model-selection arguments shared by every flash binding,
/// kept as one macro to avoid repeating the long signature. Returns the
/// assembled `(components, vapor, liquid)` for use inside a binding body.
#[allow(clippy::too_many_arguments)]
fn build_flash_pieces(
    tcs: &[f64],
    pcs: &[f64],
    omegas: &[f64],
    psat_coeffs: &[Vec<f64>],
    vl: &[f64],
    vapor_kind: &str,
    vapor_eos: Option<CubicEos>,
    liquid_kind: &str,
    liquid_eos: Option<CubicEos>,
    liquid_activity: Option<ActivityModel>,
) -> PyResult<(Vec<Component>, VaporModel, LiquidModel)> {
    let comps = flash_components(tcs, pcs, omegas, psat_coeffs, vl)?;
    let vapor = vapor_model(vapor_kind, vapor_eos)?;
    let liquid = liquid_model(liquid_kind, liquid_eos, liquid_activity)?;
    Ok((comps, vapor, liquid))
}

/// Isothermal (PT) flash. Returns `(beta, x, y, k, iterations, two_phase)`.
///
/// See the module tag docs for the `*_kind` / `*_eos` / `liquid_activity`
/// model selectors. `kij`/`aij`/`vl` may be empty; `t` in K, `p` in kPa abs.
#[pyfunction]
#[pyo3(signature = (tcs, pcs, omegas, z, t, p,
    vapor_kind, liquid_kind, vapor_eos=None, liquid_eos=None, liquid_activity=None,
    mixing_rule=MixingRule::Classical, kij=vec![], aij=vec![], vl=vec![],
    psat_coeffs=vec![], ge_model=None, tol=1e-10, max_iter=200))]
#[allow(clippy::too_many_arguments)]
fn flash_pt(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    z: Vec<f64>,
    t: f64,
    p: f64,
    vapor_kind: &str,
    liquid_kind: &str,
    vapor_eos: Option<CubicEos>,
    liquid_eos: Option<CubicEos>,
    liquid_activity: Option<ActivityModel>,
    mixing_rule: MixingRule,
    kij: Vec<Vec<f64>>,
    aij: Vec<Vec<f64>>,
    vl: Vec<f64>,
    psat_coeffs: Vec<Vec<f64>>,
    ge_model: Option<ActivityModel>,
    tol: f64,
    max_iter: usize,
) -> PyResult<(f64, Vec<f64>, Vec<f64>, Vec<f64>, usize, bool)> {
    let (comps, vapor, liquid) = build_flash_pieces(
        &tcs,
        &pcs,
        &omegas,
        &psat_coeffs,
        &vl,
        vapor_kind,
        vapor_eos,
        liquid_kind,
        liquid_eos,
        liquid_activity,
    )?;
    let spec = SystemSpec {
        components: &comps,
        vapor,
        liquid,
        mixing_rule,
        kij: &kij,
        aij: &aij,
        alpha: &[],
        vl: &vl,
        delta: &[],
        sat_models: &[],
        ge_model,
    };
    let r = flash_isothermal(&spec, t, p, &z, tol, max_iter).map_err(map_flash_err)?;
    Ok((r.beta, r.x, r.y, r.k, r.iterations, r.two_phase))
}

/// K-values `Kᵢ = yᵢ/xᵢ` at a trial `(t, p, x, y)` (a list). Model selectors
/// as in [`flash_pt`].
#[pyfunction]
#[pyo3(signature = (tcs, pcs, omegas, x, y, t, p,
    vapor_kind, liquid_kind, vapor_eos=None, liquid_eos=None, liquid_activity=None,
    mixing_rule=MixingRule::Classical, kij=vec![], aij=vec![], vl=vec![],
    psat_coeffs=vec![], ge_model=None))]
#[allow(clippy::too_many_arguments)]
fn flash_k_values_py(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    x: Vec<f64>,
    y: Vec<f64>,
    t: f64,
    p: f64,
    vapor_kind: &str,
    liquid_kind: &str,
    vapor_eos: Option<CubicEos>,
    liquid_eos: Option<CubicEos>,
    liquid_activity: Option<ActivityModel>,
    mixing_rule: MixingRule,
    kij: Vec<Vec<f64>>,
    aij: Vec<Vec<f64>>,
    vl: Vec<f64>,
    psat_coeffs: Vec<Vec<f64>>,
    ge_model: Option<ActivityModel>,
) -> PyResult<Vec<f64>> {
    let (comps, vapor, liquid) = build_flash_pieces(
        &tcs,
        &pcs,
        &omegas,
        &psat_coeffs,
        &vl,
        vapor_kind,
        vapor_eos,
        liquid_kind,
        liquid_eos,
        liquid_activity,
    )?;
    let spec = SystemSpec {
        components: &comps,
        vapor,
        liquid,
        mixing_rule,
        kij: &kij,
        aij: &aij,
        alpha: &[],
        vl: &vl,
        delta: &[],
        sat_models: &[],
        ge_model,
    };
    flash_k_values(&spec, t, p, &x, &y).map_err(map_flash_err)
}

/// Tangent-plane-distance stability analysis (cubic/φ-φ only). Returns
/// `(is_stable, trial_k, tpd)`; on a stable feed `trial_k` is empty and
/// `tpd = 0`.
#[pyfunction]
#[pyo3(signature = (tcs, pcs, omegas, z, t, p, eos,
    mixing_rule=MixingRule::Classical, kij=vec![], max_iter=100))]
#[allow(clippy::too_many_arguments)]
fn flash_stability(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    z: Vec<f64>,
    t: f64,
    p: f64,
    eos: CubicEos,
    mixing_rule: MixingRule,
    kij: Vec<Vec<f64>>,
    max_iter: usize,
) -> PyResult<(bool, Vec<f64>, f64)> {
    let comps = flash_components(&tcs, &pcs, &omegas, &[], &[])?;
    let spec = SystemSpec {
        components: &comps,
        vapor: VaporModel::Cubic(eos),
        liquid: LiquidModel::Cubic(eos),
        mixing_rule,
        kij: &kij,
        aij: &[],
        alpha: &[],
        vl: &[],
        delta: &[],
        sat_models: &[],
        ge_model: None,
    };
    match stability_analysis(&spec, t, p, &z, max_iter).map_err(map_flash_err)? {
        Stability::Stable => Ok((true, vec![], 0.0)),
        Stability::Unstable { trial_k, tpd } => Ok((false, trial_k, tpd)),
    }
}

/// Shared body for the four bubble/dew bindings: assemble the spec and call
/// `solver`, returning `(value, incipient, k)`.
#[allow(clippy::too_many_arguments)]
fn saturation_binding(
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    comp: Vec<f64>,
    fixed: f64,
    vapor_kind: &str,
    liquid_kind: &str,
    vapor_eos: Option<CubicEos>,
    liquid_eos: Option<CubicEos>,
    liquid_activity: Option<ActivityModel>,
    mixing_rule: MixingRule,
    kij: Vec<Vec<f64>>,
    aij: Vec<Vec<f64>>,
    vl: Vec<f64>,
    psat_coeffs: Vec<Vec<f64>>,
    ge_model: Option<ActivityModel>,
    tol: f64,
    max_iter: usize,
    solver: impl Fn(
        &SystemSpec,
        f64,
        &[f64],
        f64,
        usize,
    ) -> Result<crate::flash::bubble::SaturationResult, FlashError>,
) -> PyResult<(f64, Vec<f64>, Vec<f64>)> {
    let (comps, vapor, liquid) = build_flash_pieces(
        &tcs,
        &pcs,
        &omegas,
        &psat_coeffs,
        &vl,
        vapor_kind,
        vapor_eos,
        liquid_kind,
        liquid_eos,
        liquid_activity,
    )?;
    let spec = SystemSpec {
        components: &comps,
        vapor,
        liquid,
        mixing_rule,
        kij: &kij,
        aij: &aij,
        alpha: &[],
        vl: &vl,
        delta: &[],
        sat_models: &[],
        ge_model,
    };
    let r = solver(&spec, fixed, &comp, tol, max_iter).map_err(map_flash_err)?;
    Ok((r.value, r.incipient, r.k))
}

macro_rules! saturation_pyfn {
    ($name:ident, $solver:path, $fixed_doc:literal, $comp_doc:literal) => {
        #[doc = concat!("Saturation point — ", $fixed_doc, ". Returns `(value, incipient, k)` where ", $comp_doc, ".")]
        #[pyfunction]
        #[pyo3(signature = (tcs, pcs, omegas, comp, fixed, vapor_kind, liquid_kind,
            vapor_eos=None, liquid_eos=None, liquid_activity=None,
            mixing_rule=MixingRule::Classical, kij=vec![], aij=vec![], vl=vec![],
            psat_coeffs=vec![], ge_model=None, tol=1e-9, max_iter=200))]
        #[allow(clippy::too_many_arguments)]
        fn $name(
            tcs: Vec<f64>,
            pcs: Vec<f64>,
            omegas: Vec<f64>,
            comp: Vec<f64>,
            fixed: f64,
            vapor_kind: &str,
            liquid_kind: &str,
            vapor_eos: Option<CubicEos>,
            liquid_eos: Option<CubicEos>,
            liquid_activity: Option<ActivityModel>,
            mixing_rule: MixingRule,
            kij: Vec<Vec<f64>>,
            aij: Vec<Vec<f64>>,
            vl: Vec<f64>,
            psat_coeffs: Vec<Vec<f64>>,
            ge_model: Option<ActivityModel>,
            tol: f64,
            max_iter: usize,
        ) -> PyResult<(f64, Vec<f64>, Vec<f64>)> {
            saturation_binding(
                tcs, pcs, omegas, comp, fixed, vapor_kind, liquid_kind, vapor_eos, liquid_eos,
                liquid_activity, mixing_rule, kij, aij, vl, psat_coeffs, ge_model, tol, max_iter,
                $solver,
            )
        }
    };
}

saturation_pyfn!(
    bubble_pressure_py,
    bubble_pressure,
    "bubble pressure given `(fixed=T [K], comp=x)`",
    "`value` = P in kPa and `incipient` = the vapor y"
);
saturation_pyfn!(
    bubble_temperature_py,
    bubble_temperature,
    "bubble temperature given `(fixed=P [kPa], comp=x)`",
    "`value` = T in K and `incipient` = the vapor y"
);
saturation_pyfn!(
    dew_pressure_py,
    dew_pressure,
    "dew pressure given `(fixed=T [K], comp=y)`",
    "`value` = P in kPa and `incipient` = the liquid x"
);
saturation_pyfn!(
    dew_temperature_py,
    dew_temperature,
    "dew temperature given `(fixed=P [kPa], comp=y)`",
    "`value` = T in K and `incipient` = the liquid x"
);

/// Mixture critical point (Heidemann §G) for a two-parameter cubic EOS with
/// classical mixing. Returns `(Tc [K], Pc [kPa], Vc [m³/kmol])`.
///
/// `t_init` is an initial temperature guess (pass 0.0 to use the
/// mole-fraction-average Tc).
#[pyfunction]
#[pyo3(signature = (eos, tcs, pcs, omegas, z, t_init=0.0, kij=vec![], max_iter=200))]
#[allow(clippy::too_many_arguments)]
fn critical_point_py(
    eos: CubicEos,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    z: Vec<f64>,
    t_init: f64,
    kij: Vec<Vec<f64>>,
    max_iter: usize,
) -> PyResult<(f64, f64, f64)> {
    let comps = flash_components(&tcs, &pcs, &omegas, &[], &[])?;
    let spec = SystemSpec {
        components: &comps,
        vapor: VaporModel::Cubic(eos),
        liquid: LiquidModel::Cubic(eos),
        mixing_rule: MixingRule::Classical,
        kij: &kij,
        aij: &[],
        alpha: &[],
        vl: &[],
        delta: &[],
        sat_models: &[],
        ge_model: None,
    };
    let cp = critical_point(&spec, &z, t_init, max_iter).map_err(map_flash_err)?;
    Ok((cp.tc, cp.pc, cp.vc))
}

/// Adiabatic (PH) flash for a φ-φ cubic system. Returns
/// `(t, beta, x, y, enthalpy)` — the flash temperature (K), vapor fraction,
/// phase compositions, and the stream enthalpy (kJ/kmol).
#[pyfunction]
#[pyo3(signature = (eos, tcs, pcs, omegas, cp_coeffs, z, p, h_feed, t_lo, t_hi,
    t_ref=298.15, p_ref=101.325, kij=vec![], tol=1e-4, max_iter=200))]
#[allow(clippy::too_many_arguments)]
fn flash_adiabatic_py(
    eos: CubicEos,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    cp_coeffs: Vec<Vec<f64>>,
    z: Vec<f64>,
    p: f64,
    h_feed: f64,
    t_lo: f64,
    t_hi: f64,
    t_ref: f64,
    p_ref: f64,
    kij: Vec<Vec<f64>>,
    tol: f64,
    max_iter: usize,
) -> PyResult<(f64, f64, Vec<f64>, Vec<f64>, f64)> {
    let comps = mix_components(&tcs, &pcs, &omegas, &cp_coeffs)?;
    let spec = SystemSpec {
        components: &comps,
        vapor: VaporModel::Cubic(eos),
        liquid: LiquidModel::Cubic(eos),
        mixing_rule: MixingRule::Classical,
        kij: &kij,
        aij: &[],
        alpha: &[],
        vl: &[],
        delta: &[],
        sat_models: &[],
        ge_model: None,
    };
    let r = flash_adiabatic(
        &spec, p, &z, h_feed, t_ref, p_ref, t_lo, t_hi, tol, max_iter,
    )
    .map_err(map_flash_err)?;
    Ok((r.t, r.flash.beta, r.flash.x, r.flash.y, r.enthalpy))
}

/// Fit a binary EOS `k₁₂` to bubble-pressure data. `data` is a list of
/// `(T [K], x1, P_exp [kPa])` triples. Returns `(kij, sse, rmse)`.
#[pyfunction]
#[pyo3(signature = (eos, tcs, pcs, omegas, psat_coeffs, data, k_lo=-0.1, k_hi=0.3, tol=1e-6, max_iter=100))]
#[allow(clippy::too_many_arguments)]
fn fit_kij_py(
    eos: CubicEos,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    psat_coeffs: Vec<Vec<f64>>,
    data: Vec<(f64, f64, f64)>,
    k_lo: f64,
    k_hi: f64,
    tol: f64,
    max_iter: usize,
) -> PyResult<(f64, f64, f64)> {
    let comps = flash_components(&tcs, &pcs, &omegas, &psat_coeffs, &[])?;
    let pts: Vec<BubblePoint> = data
        .into_iter()
        .map(|(t, x1, p_exp)| BubblePoint { t, x1, p_exp })
        .collect();
    let fit = fit_kij(eos, &comps, &pts, k_lo, k_hi, tol, max_iter).map_err(map_flash_err)?;
    Ok((fit.kij, fit.sse, fit.rmse))
}

/// Fit binary activity-model parameters `(A₁₂, A₂₁)` to bubble-pressure data
/// by Levenberg–Marquardt. `data` is a list of `(T [K], x1, P_exp [kPa])`.
/// `alpha` is the fixed 2×2 NRTL non-randomness matrix (pass `[]` for the
/// non-NRTL models; α is held constant, only the energies are fitted).
/// Returns `(a12, a21, sse, rmse, iterations)`.
#[pyfunction]
#[pyo3(signature = (model, tcs, pcs, omegas, psat_coeffs, data, a12_0, a21_0,
    alpha=vec![], vl=vec![], tol=1e-10, max_iter=100))]
#[allow(clippy::too_many_arguments)]
fn fit_aij_py(
    model: ActivityModel,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    psat_coeffs: Vec<Vec<f64>>,
    data: Vec<(f64, f64, f64)>,
    a12_0: f64,
    a21_0: f64,
    alpha: Vec<Vec<f64>>,
    vl: Vec<f64>,
    tol: f64,
    max_iter: usize,
) -> PyResult<(f64, f64, f64, f64, usize)> {
    let comps = flash_components(&tcs, &pcs, &omegas, &psat_coeffs, &vl)?;
    let pts: Vec<AijBubblePoint> = data
        .into_iter()
        .map(|(t, x1, p_exp)| AijBubblePoint { t, x1, p_exp })
        .collect();
    let fit = fit_aij(
        model, &comps, &alpha, &vl, &pts, a12_0, a21_0, tol, max_iter,
    )
    .map_err(map_flash_err)?;
    Ok((fit.a12, fit.a21, fit.sse, fit.rmse, fit.iterations))
}

/// Trace the phase envelope of a φ-φ cubic mixture at composition `z`
/// (Michelsen continuation, §K). Returns a list of `(T [K], P [kPa])`
/// points from the low-pressure seed up through the critical region.
#[pyfunction]
#[pyo3(signature = (eos, tcs, pcs, omegas, z, p_start=100.0, kij=vec![], max_points=60))]
#[allow(clippy::too_many_arguments)]
fn trace_envelope_py(
    eos: CubicEos,
    tcs: Vec<f64>,
    pcs: Vec<f64>,
    omegas: Vec<f64>,
    z: Vec<f64>,
    p_start: f64,
    kij: Vec<Vec<f64>>,
    max_points: usize,
) -> PyResult<Vec<(f64, f64)>> {
    let comps = flash_components(&tcs, &pcs, &omegas, &[], &[])?;
    let spec = SystemSpec {
        components: &comps,
        vapor: VaporModel::Cubic(eos),
        liquid: LiquidModel::Cubic(eos),
        mixing_rule: MixingRule::Classical,
        kij: &kij,
        aij: &[],
        alpha: &[],
        vl: &[],
        delta: &[],
        sat_models: &[],
        ge_model: None,
    };
    let pts = trace_envelope(&spec, &z, p_start, max_points).map_err(map_flash_err)?;
    Ok(pts.into_iter().map(|p| (p.t, p.p)).collect())
}

// ──────────────────────────────────────────────────────────────────────
// M12.2 — Rust-side component database bindings (gap G3).
//
// These expose `vle_thermo::db::{component, available}` so the Rust DB is
// testable through the wheel (the M5+ same-commit binding rule). The wheel's
// *primary* component loader stays the pure-Python `vle.components` — no churn
// there — so these exist for parity testing and for callers that want the
// crate's own catalog. Gated on `component-db`, which the `python` feature
// turns on transitively (see engine/Cargo.toml).
//
// A looked-up component is returned as a plain `dict` (rather than a bespoke
// pyclass): it carries every field the JSON supplies, reads naturally on the
// Python side (`db_component("toluene")["tc"]`), and keeps the FFI surface
// minimal — the same "no `Component` pyclass" convention as the M7 scalar
// bindings above.
// ──────────────────────────────────────────────────────────────────────

/// Look up a bundled component by name (case-insensitive), or `None` if unknown.
///
/// Returns a dict of canonical-unit properties: `name`, `tc`/`tb` in **K**,
/// `pc` in **kPa** (absolute), `vc`/`liquid_volume` in **cm³/mol**, `mw` in
/// **g/mol**, `omega`/`zc` dimensionless, `psat_coeffs` the reduced-Antoine
/// coefficient list `ln(Psat/Pc) = a₁ − a₂/(a₃ + T)` (T in **K**), and
/// `cp_coeffs` the dimensionless ideal-gas Cp°/R polynomial `Σₖ aₖ·Tᵏ`
/// (T in **K**). Name normalization matches `vle.components.get`.
#[cfg(feature = "component-db")]
#[pyfunction]
fn db_component<'py>(
    py: Python<'py>,
    name: &str,
) -> PyResult<Option<Bound<'py, pyo3::types::PyDict>>> {
    let Some(c) = crate::db::component(name) else {
        return Ok(None);
    };
    let d = pyo3::types::PyDict::new_bound(py);
    d.set_item("name", c.name)?;
    d.set_item("tc", c.tc)?;
    d.set_item("pc", c.pc)?;
    d.set_item("vc", c.vc)?;
    d.set_item("zc", c.zc)?;
    d.set_item("omega", c.omega)?;
    d.set_item("tb", c.tb)?;
    d.set_item("mw", c.mw)?;
    d.set_item("cp_coeffs", c.cp_coeffs.to_vec())?;
    d.set_item("psat_coeffs", c.psat_coeffs)?;
    d.set_item("liquid_volume", c.liquid_volume)?;
    Ok(Some(d))
}

/// Names of all bundled components, sorted. Mirrors `vle.components.available()`.
#[cfg(feature = "component-db")]
#[pyfunction]
fn db_available() -> Vec<String> {
    crate::db::available()
}

/// PyO3 module entry point.
///
/// Maturin builds this into `vle/_engine.<platform>.<ext>` and Python
/// imports it as `vle._engine`. The function name (`_engine`) must match
/// the `module-name` set in `python/pyproject.toml`'s `[tool.maturin]`.
#[pymodule]
fn _engine(_py: Python<'_>, m: &Bound<'_, PyModule>) -> PyResult<()> {
    m.add_function(wrap_pyfunction!(version, m)?)?;
    m.add_function(wrap_pyfunction!(default_units_toml, m)?)?;

    // M6.1 numerics — solvers + utility functions.
    m.add_function(wrap_pyfunction!(solve_cubic, m)?)?;
    m.add_function(wrap_pyfunction!(brent, m)?)?;
    m.add_function(wrap_pyfunction!(illinois, m)?)?;
    m.add_function(wrap_pyfunction!(halley, m)?)?;
    m.add_function(wrap_pyfunction!(broyden, m)?)?;
    m.add_function(wrap_pyfunction!(sum_frac_residual, m)?)?;
    m.add_function(wrap_pyfunction!(norm_l1, m)?)?;
    m.add_function(wrap_pyfunction!(norm_l2, m)?)?;
    m.add_function(wrap_pyfunction!(norm_linf, m)?)?;

    // The four model-selection enums. Each is `#[pyclass(eq, eq_int)]`
    // at its definition site; here we just register the class with the
    // module so `vle._engine.CubicEos` etc. resolve from Python.
    m.add_class::<crate::eos::CubicEos>()?;
    m.add_class::<crate::activity::ActivityModel>()?;
    m.add_class::<crate::mixing::MixingRule>()?;
    m.add_class::<crate::saturation::SatPressureModel>()?;
    m.add_class::<crate::eos::PhaseId>()?;
    m.add_class::<crate::eos::ChaoSeaderSpecies>()?;

    // M19 petroleum characterization — the `Assay` class plus the `petro_*`
    // correlation functions (see engine/src/py_petroleum.rs).
    crate::py_petroleum::register(m)?;
    crate::py_refinery::register(m)?;

    // M7.1 cubic-EOS bindings.
    m.add_function(wrap_pyfunction!(eos_alpha, m)?)?;
    m.add_function(wrap_pyfunction!(eos_d_alpha_d_tr, m)?)?;
    // M7.2 — α variants taking the full polar/fitted parameter set.
    m.add_function(wrap_pyfunction!(eos_alpha_ex, m)?)?;
    m.add_function(wrap_pyfunction!(eos_d_alpha_d_tr_ex, m)?)?;
    m.add_function(wrap_pyfunction!(eos_family_constants, m)?)?;
    m.add_function(wrap_pyfunction!(eos_z_factor, m)?)?;
    m.add_function(wrap_pyfunction!(eos_ln_phi_pure, m)?)?;
    m.add_function(wrap_pyfunction!(eos_h_departure_rt, m)?)?;
    m.add_function(wrap_pyfunction!(eos_s_departure_r, m)?)?;
    m.add_function(wrap_pyfunction!(chao_seader_ln_phi, m)?)?;

    // M7.5 Antoine saturation bindings.
    m.add_function(wrap_pyfunction!(antoine_psat, m)?)?;
    m.add_function(wrap_pyfunction!(antoine_d_psat_dt, m)?)?;

    // M7.4 advanced saturation models + OL-family α.
    m.add_function(wrap_pyfunction!(sat_psat, m)?)?;
    m.add_function(wrap_pyfunction!(sat_d_psat_dt, m)?)?;
    m.add_function(wrap_pyfunction!(sat_reduced_psat, m)?)?;
    m.add_function(wrap_pyfunction!(sat_maxwell, m)?)?;
    m.add_function(wrap_pyfunction!(boiling_temperature, m)?)?;
    m.add_function(wrap_pyfunction!(poynting_factor, m)?)?;
    m.add_function(wrap_pyfunction!(eos_alpha_ol, m)?)?;
    m.add_function(wrap_pyfunction!(eos_d_alpha_d_tr_ol, m)?)?;

    // M7.6 Virial bindings.
    m.add_function(wrap_pyfunction!(virial_pitzer_b0, m)?)?;
    m.add_function(wrap_pyfunction!(virial_pitzer_b1, m)?)?;
    m.add_function(wrap_pyfunction!(virial_b_pure, m)?)?;
    m.add_function(wrap_pyfunction!(virial_d_b_d_t_pure, m)?)?;
    m.add_function(wrap_pyfunction!(virial_z, m)?)?;
    m.add_function(wrap_pyfunction!(virial_ln_phi, m)?)?;
    m.add_function(wrap_pyfunction!(virial_h_dep_rt, m)?)?;
    m.add_function(wrap_pyfunction!(virial_s_dep_r, m)?)?;
    m.add_function(wrap_pyfunction!(virial_b_mix_py, m)?)?;
    m.add_function(wrap_pyfunction!(virial_ln_phi_mix, m)?)?;

    // M8.1 liquid molar volume + activity-coefficient models.
    m.add_class::<crate::liquid_volume::VolumeModel>()?;
    m.add_function(wrap_pyfunction!(liquid_molar_volume, m)?)?;
    m.add_function(wrap_pyfunction!(activity_ln_gamma, m)?)?;
    m.add_function(wrap_pyfunction!(activity_excess_gibbs, m)?)?;
    m.add_function(wrap_pyfunction!(activity_excess_enthalpy, m)?)?;
    m.add_function(wrap_pyfunction!(activity_excess_entropy, m)?)?;

    // M8.3 mixture core — mixing rules, multicomponent fugacity, the exact
    // composition-derivative Jacobian, and 3-param / Chao-Seader paths.
    m.add_function(wrap_pyfunction!(mixture_z, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_ln_phi, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_d_ln_phi_d_n, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_chao_seader_ln_phi, m)?)?;

    // M8.4 mixture energy properties — departures + ideal + total assembly.
    m.add_function(wrap_pyfunction!(mixture_h_departure_rt, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_s_departure_r, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_ideal_enthalpy, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_ideal_entropy, m)?)?;
    m.add_function(wrap_pyfunction!(mixture_phase_enthalpy_entropy, m)?)?;

    // M9 flash / bubble / dew / stability.
    m.add_function(wrap_pyfunction!(rachford_rice, m)?)?;
    m.add_function(wrap_pyfunction!(flash_pt, m)?)?;
    m.add_function(wrap_pyfunction!(flash_k_values_py, m)?)?;
    m.add_function(wrap_pyfunction!(flash_stability, m)?)?;
    m.add_function(wrap_pyfunction!(bubble_pressure_py, m)?)?;
    m.add_function(wrap_pyfunction!(bubble_temperature_py, m)?)?;
    m.add_function(wrap_pyfunction!(dew_pressure_py, m)?)?;
    m.add_function(wrap_pyfunction!(dew_temperature_py, m)?)?;

    // M9 critical point, adiabatic flash, phase envelope, kij/Aij regression.
    m.add_function(wrap_pyfunction!(critical_point_py, m)?)?;
    m.add_function(wrap_pyfunction!(flash_adiabatic_py, m)?)?;
    m.add_function(wrap_pyfunction!(trace_envelope_py, m)?)?;
    m.add_function(wrap_pyfunction!(fit_kij_py, m)?)?;
    m.add_function(wrap_pyfunction!(fit_aij_py, m)?)?;

    // M10 — the persistent System handle + batch numpy API (Track D).
    m.add_class::<crate::py_system::System>()?;

    // M13 — IAPWS-IF97 steam tables (STEAM_TABLES_PLAN.md), surfaced as
    // `vle.steam`. Gated on the `steam` feature that `python` turns on.
    {
        use crate::py_steam::*;
        m.add_class::<SteamState>()?;
        m.add_class::<SatState>()?;
        m.add_function(wrap_pyfunction!(steam_tp, m)?)?;
        m.add_function(wrap_pyfunction!(steam_tx, m)?)?;
        m.add_function(wrap_pyfunction!(steam_px, m)?)?;
        m.add_function(wrap_pyfunction!(steam_ph, m)?)?;
        m.add_function(wrap_pyfunction!(steam_ps, m)?)?;
        m.add_function(wrap_pyfunction!(steam_sat_t, m)?)?;
        m.add_function(wrap_pyfunction!(steam_sat_p, m)?)?;
        m.add_function(wrap_pyfunction!(steam_psat, m)?)?;
        m.add_function(wrap_pyfunction!(steam_tsat, m)?)?;
        m.add_function(wrap_pyfunction!(steam_psat_derivative, m)?)?;
        m.add_function(wrap_pyfunction!(steam_latent_heat, m)?)?;
        m.add_function(wrap_pyfunction!(steam_viscosity, m)?)?;
        m.add_function(wrap_pyfunction!(steam_thermal_conductivity, m)?)?;
        m.add_function(wrap_pyfunction!(steam_surface_tension, m)?)?;
        m.add_function(wrap_pyfunction!(steam_tp_batch, m)?)?;
        m.add_function(wrap_pyfunction!(steam_ph_batch, m)?)?;
        m.add_function(wrap_pyfunction!(steam_sat_t_batch, m)?)?;
        m.add_function(wrap_pyfunction!(steam_transport_batch, m)?)?;
    }

    // M12.2 — Rust-side component database lookups (gap G3). Guarded on the
    // feature the `python` feature turns on, so a hypothetical wheel built
    // without it still compiles.
    #[cfg(feature = "component-db")]
    {
        m.add_function(wrap_pyfunction!(db_component, m)?)?;
        m.add_function(wrap_pyfunction!(db_available, m)?)?;
    }

    Ok(())
}