nereids_physics/

rmatrix_limited.rs

//! R-Matrix Limited (LRF=7) cross-section calculation.
//!
//! Computes energy-dependent cross-sections (total, elastic, capture, fission)
//! from ENDF LRF=7 resonance parameters using the full multi-channel R-matrix
//! formalism.
//!
//! ## Relationship to Reich-Moore (LRF=3)
//!
//! Reich-Moore is a special case of R-matrix theory that *eliminates* the
//! capture channel via the Wigner-Eisenbud approximation, collapsing the
//! level matrix to a scalar per spin group. LRF=7 retains all channels
//! (elastic, capture, fission) explicitly, requiring an NCH×NCH complex
//! matrix inversion per spin group per energy point.
//!
//! ## Key Formulas
//!
//! (ENDF-6 Formats Manual Appendix D; SAMMY manual §3.1; SAMMY `rml/mrml07.f` Pgh + `mrml11.f` Setxqx)
//!
//! ```text
//! R-matrix:
//!   R_cc'(E) = Σ_n γ_nc · γ_nc' / (E_n - E)   [complex for KRM=3, NCH×NCH, symmetric]
//!   γ_nc = reduced width amplitude for resonance n in channel c (eV^{1/2})
//!
//! Level denominator per channel:
//!   L_c(E) = S_c(E) - B_c + i·P_c(E)
//!   P_c = penetrability,  S_c = shift factor,  B_c = boundary condition
//!
//! Reduced level matrix (SAMMY "Ymat"):
//!   Ỹ_cc'(E) = δ_cc' / L_c(E) - R_cc'(E)   [complex, NCH×NCH]
//!   → Ỹinv = Ỹ⁻¹  (invert Ỹ)
//!
//! Intermediate matrix (SAMMY "XXXX"):
//!   Ξ_cc'(E) = (√P_c / L_c) · (Ỹinv · R)_cc' · √P_c'
//!
//! Collision matrix (SAMMY manual eq. III.D.4):
//!   W_cc' = δ_cc' + 2i·Ξ_cc'
//!   U_cc' = Ω_c · W_cc' · Ω_c'    where Ω_c = exp(-iφ_c)
//!
//!   TRUTH SOURCE: SAMMY rml/mrml11.f lines 14-18 (W = I + 2i·XXXX)
//!   and lines 84-88 (elastic formula consistent with e^{-2iφ}).
//!   Unitarity: |U| ≤ 1 always; hard sphere (R=0) → U = exp(2iφ)·I  ✓
//!
//! Cross sections per spin group (J,π), summed over entrance neutron channels c0:
//!   σ_total   = Σ_{c0} 2·(π/k²)·g_J·(1 - Re(U_{c0,c0}))
//!   σ_elastic = Σ_{c0} (π/k²)·g_J·|1 - U_{c0,c0}|²
//!   σ_fission = Σ_{c0} Σ_{c'∈fission} (π/k²)·g_J·|U_{c0,c'}|²
//!   σ_capture = σ_total - σ_elastic - σ_fission
//! ```
//!
//! ## Photon Channels
//!
//! For channels involving zero-rest-mass particles (photons, MA=0), the
//! penetrability is set to 1.0 and phase shifts to 0.0, following the
//! standard R-matrix convention (ENDF-6 Formats Manual §2.2.1.6 Note 4).
//!
//! ## SAMMY Reference
//! - `rml/mrml01.f` — LRF=7 reader (Scan_File_2, particle pair loop)
//! - `rml/mrml09.f` — Level matrix inversion (Yinvrs, Xspfa, Xspsl)
//! - `rml/mrml11.f` — Cross-section calculation (Sectio, Setxqx)
//! - SAMMY manual §3.1 (multi-channel R-matrix)

use num_complex::Complex64;

use nereids_core::constants::{LOG_FLOOR, NEAR_ZERO_FLOOR, PIVOT_FLOOR, QUANTUM_NUMBER_EPS};
use nereids_endf::resonance::{ParticlePair, RmlData, SpinGroup};

use crate::{channel, coulomb, penetrability};

/// Pre-allocated workspace for RML cross-section evaluation.
///
/// Eliminates per-energy-point allocation of NCH×NCH complex matrices
/// (`r_cplx`, `y_tilde`, `y_inv`, `xq`, `xxxx`, `u`) and per-channel
/// vectors (`p_c`, `s_c`, `phi_c`, flags, etc.).
///
/// Flat `Vec<Complex64>` buffers store matrices in row-major order.
/// For typical NCH=3-6, this avoids ~12 small heap allocations per
/// energy point per spin group.
struct RmlWorkspace {
    // ── NCH×NCH complex matrix buffers (flat, row-major) ──────────────
    r_cplx: Vec<Complex64>,
    y_tilde: Vec<Complex64>,
    y_inv: Vec<Complex64>,
    xq: Vec<Complex64>,
    xxxx: Vec<Complex64>,
    u: Vec<Complex64>,
    // ── Augmented matrix for Gauss-Jordan inversion (NCH × 2·NCH) ─────
    aug: Vec<Complex64>,
    // ── Temp row for elimination ───────────────────────────────────────
    aug_tmp: Vec<Complex64>,
    // ── Per-channel vectors ───────────────────────────────────────────
    p_c: Vec<f64>,
    s_c: Vec<f64>,
    phi_c: Vec<f64>,
    l_c: Vec<Complex64>,
    sqrt_p: Vec<f64>,
    omega: Vec<Complex64>,
    is_entrance: Vec<bool>,
    is_fission: Vec<bool>,
    is_capture: Vec<bool>,
    is_inelastic: Vec<bool>,
    is_closed: Vec<bool>,
    // ── Per-resonance scratch ─────────────────────────────────────────
    gamma_vals: Vec<f64>,
}

impl RmlWorkspace {
    fn new() -> Self {
        Self {
            r_cplx: Vec::new(),
            y_tilde: Vec::new(),
            y_inv: Vec::new(),
            xq: Vec::new(),
            xxxx: Vec::new(),
            u: Vec::new(),
            aug: Vec::new(),
            aug_tmp: Vec::new(),
            p_c: Vec::new(),
            s_c: Vec::new(),
            phi_c: Vec::new(),
            l_c: Vec::new(),
            sqrt_p: Vec::new(),
            omega: Vec::new(),
            is_entrance: Vec::new(),
            is_fission: Vec::new(),
            is_capture: Vec::new(),
            is_inelastic: Vec::new(),
            is_closed: Vec::new(),
            gamma_vals: Vec::new(),
        }
    }

    /// Resize all buffers for a spin group with `nch` channels and
    /// `max_widths` resonance width entries, then zero them out.
    fn resize_and_clear(&mut self, nch: usize, max_widths: usize) {
        let nn = nch * nch;
        let aug_len = nch * 2 * nch;

        // Complex matrix buffers — resize then fill with zero.
        resize_and_zero(&mut self.r_cplx, nn);
        resize_and_zero(&mut self.y_tilde, nn);
        resize_and_zero(&mut self.y_inv, nn);
        resize_and_zero(&mut self.xq, nn);
        resize_and_zero(&mut self.xxxx, nn);
        resize_and_zero(&mut self.u, nn);
        resize_and_zero(&mut self.aug, aug_len);
        resize_and_zero(&mut self.aug_tmp, 2 * nch);

        // Per-channel real/bool vectors.
        resize_and_zero_f64(&mut self.p_c, nch);
        resize_and_zero_f64(&mut self.s_c, nch);
        resize_and_zero_f64(&mut self.phi_c, nch);
        resize_and_zero(&mut self.l_c, nch);
        resize_and_zero_f64(&mut self.sqrt_p, nch);
        resize_and_zero(&mut self.omega, nch);
        resize_and_zero_bool(&mut self.is_entrance, nch);
        resize_and_zero_bool(&mut self.is_fission, nch);
        resize_and_zero_bool(&mut self.is_capture, nch);
        resize_and_zero_bool(&mut self.is_inelastic, nch);
        resize_and_zero_bool(&mut self.is_closed, nch);

        // Per-resonance scratch.
        resize_and_zero_f64(&mut self.gamma_vals, max_widths);
    }
}

fn resize_and_zero(buf: &mut Vec<Complex64>, len: usize) {
    buf.clear();
    buf.resize(len, Complex64::ZERO);
}

fn resize_and_zero_f64(buf: &mut Vec<f64>, len: usize) {
    buf.clear();
    buf.resize(len, 0.0);
}

fn resize_and_zero_bool(buf: &mut Vec<bool>, len: usize) {
    buf.clear();
    buf.resize(len, false);
}

/// Compute cross-section contributions from an LRF=7 energy range.
///
/// Returns `(total, elastic, capture, fission)` in barns.
///
/// Iterates over all spin groups (J,π), sums their contributions.
/// A single `RmlWorkspace` is allocated once and reused across spin groups.
pub fn cross_sections_for_rml_range(rml: &RmlData, energy_ev: f64) -> (f64, f64, f64, f64) {
    let mut total = 0.0;
    let mut elastic = 0.0;
    let mut capture = 0.0;
    let mut fission = 0.0;

    let mut ws = RmlWorkspace::new();

    for sg in &rml.spin_groups {
        let (t, e, cap, fis) = spin_group_cross_sections(
            sg,
            &rml.particle_pairs,
            energy_ev,
            rml.awr,
            rml.target_spin,
            rml.krm,
            &mut ws,
        );
        total += t;
        elastic += e;
        capture += cap;
        fission += fis;
    }

    (total, elastic, capture, fission)
}

/// Cross-section contribution from a single spin group (J,π).
///
/// Returns (total, elastic, capture, fission) in barns.
///
/// The caller-owned `ws` workspace is resized and reused across spin groups
/// for a given energy evaluation to avoid per-call heap allocations.
fn spin_group_cross_sections(
    sg: &SpinGroup,
    particle_pairs: &[ParticlePair],
    energy_ev: f64,
    awr: f64,
    target_spin: f64,
    krm: u32,
    ws: &mut RmlWorkspace,
) -> (f64, f64, f64, f64) {
    let nch = sg.channels.len();
    if nch == 0 {
        return (0.0, 0.0, 0.0, 0.0);
    }

    // KRM guard: the parser rejects KRM values other than 2 and 3 at load time,
    // so reaching here with an unsupported KRM indicates a programming error.
    // Panic rather than silently returning zero physics, which would look valid
    // to callers.  Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml01.f KRM field.
    assert!(
        krm == 2 || krm == 3,
        "spin_group_cross_sections called with unsupported KRM={krm}; \
         should have been rejected at parse time"
    );

    // P2b: NRS=0 is valid — the R-matrix is zero but the hard-sphere phase shift
    // still produces nonzero potential-scattering cross sections (σ_el = 4π|Ω-1|²/k²
    // per spin group).  Do NOT return early here; the resonance loop simply executes
    // zero iterations, leaving R = 0, which is exactly the hard-sphere limit.

    // Resize workspace buffers for this spin group's channel count.
    let max_widths = sg
        .resonances
        .iter()
        .map(|r| r.widths.len())
        .max()
        .unwrap_or(nch)
        .max(nch);
    ws.resize_and_clear(nch, max_widths);

    let g_j = channel::statistical_weight(sg.j, target_spin);
    let pok2 = channel::pi_over_k_squared_barns(energy_ev, awr);

    // Entrance-channel CM energy: E_cm = E_lab × AWR/(1+AWR).
    // Each exit channel adds its Q-value to get its own available energy.
    // Reference: SAMMY rml/mrml03.f Fxradi — channel thresholds via Q.
    let e_cm = channel::lab_to_cm_energy(energy_ev, awr);

    for (c, ch) in sg.channels.iter().enumerate() {
        // P3: particle_pair_idx must be a valid index. The old `.min(len-1)` clamped
        // silently, misclassifying any channel with an OOB index as the last pair.
        // An OOB value indicates corrupted ENDF data; let Rust's bounds check panic.
        let pp = &particle_pairs[ch.particle_pair_idx];
        ws.is_entrance[c] = pp.mt == 2;
        ws.is_fission[c] = pp.mt == 18;
        ws.is_capture[c] = pp.mt == 102;
        // Inelastic neutron channels (MT=51+): massive particle, not elastic/fission/capture.
        // Their flux appears in σ_total (optical theorem) but must not be assigned to capture.
        // Reference: ENDF MT number conventions, §3.4; SAMMY rml/mrml11.f Sectio.
        ws.is_inelastic[c] =
            pp.ma >= 0.5 && !ws.is_entrance[c] && !ws.is_fission[c] && !ws.is_capture[c];

        if pp.ma < 0.5 {
            // Photon channel (MA = 0): P=1, S=0, φ=0.
            // Convention per ENDF-6 Formats Manual §2.2.1.6 Note 4.
            // SAMMY: rml/mrml07.f sets penetrability = 1 for massless particles.
            ws.p_c[c] = 1.0;
            ws.s_c[c] = 0.0;
            ws.phi_c[c] = 0.0;
        } else {
            // Massive particle channel: channel-specific kinematics (P1).
            // E_c = E_cm + Q (CM kinetic energy in this exit channel).
            // Reference: SAMMY rml/mrml03.f Fxradi — Zke = Twomhb*sqrt(Redmas*Factor)
            let e_c = e_cm + pp.q;
            if e_c <= 0.0 {
                // Closed channel (below threshold): P_c = 0, φ_c = 0.
                // S_c depends on SHF:
                //   SHF=0: convention is S_c = B_c; L_c = 0 when B_c = 0 (common).
                //   SHF=1: S_c is the analytic shift factor at imaginary argument
                //     ρ = iκ, which is real and finite.  L_c = (S_c − B_c) is generally
                //     non-zero and its dispersive contribution must be preserved.
                // Reference: SAMMY rml/mrml07.f Pgh — PH = 1/(S−B+iP).
                ws.p_c[c] = 0.0;
                ws.phi_c[c] = 0.0;
                ws.is_closed[c] = true;
                // SHF=0: S_c = B_c so (S_c − B_c) = 0 in the level matrix.
                // SHF=1: S_c is the analytic shift at imaginary argument ρ = iκ.
                //   For non-Coulomb channels we use the Blatt-Weisskopf formula.
                //   For Coulomb channels the imaginary-argument Coulomb shift is
                //   not yet implemented; fall back to S_c = B_c.  This matches
                //   SAMMY's convention: mrml07.f ELSE branch (Su ≤ Echan) sets
                //   Elinvr=1/Elinvi=0 (i.e. L_c = 1) for all closed channels,
                //   Coulomb and non-Coulomb alike, without calling Pghcou.
                let is_coulomb = pp.za.abs() > 0.5 && pp.zb.abs() > 0.5;
                ws.s_c[c] = if pp.shf == 1 && !is_coulomb {
                    let redmas = pp.ma * pp.mb / (pp.ma + pp.mb);
                    let kappa = channel::wave_number_from_cm(e_c.abs(), redmas);
                    penetrability::shift_factor_closed(ch.l, kappa * ch.effective_radius)
                } else {
                    ch.boundary
                };
            } else {
                // Channel wave number from reduced mass μ = MA·MB/(MA+MB).
                // For elastic (MA=1, MB=AWR): k_c = wave_number(E_lab, AWR) [identical].
                let redmas = pp.ma * pp.mb / (pp.ma + pp.mb);
                let k_c = channel::wave_number_from_cm(e_c, redmas);
                // APE (effective radius) for P_c, S_c; APT (true radius) for φ_c.
                // Reference: SAMMY rml/mrml07.f (Pgh, Sinsix, Pf subroutines)
                let rho_eff = k_c * ch.effective_radius;
                let rho_true = k_c * ch.true_radius;
                // ── Coulomb vs hard-sphere routing ───────────────────────────
                // SAMMY rml/mrml07.f Pgh — `if (Zeta(I).NE.Zero)` branch.
                // Both particles charged → Coulomb wave functions F_L / G_L.
                // One neutral (za=0 or zb=0) → hard-sphere Blatt-Weisskopf.
                if pp.za.abs() > 0.5 && pp.zb.abs() > 0.5 {
                    // Coulomb channel (e.g. n+α→p+X, (n,p), fission fragments).
                    // Two CF1+CF2 solves: rho_eff for P_c/S_c, rho_true for φ_c
                    // (different radii — cannot reuse the same (F,G) pair).
                    // Reference: SAMMY rml/mrml07.f Pgh — Pghcou, then Sinsix/Pf.
                    let eta = coulomb::sommerfeld_eta(pp.za, pp.zb, pp.ma, pp.mb, e_c);
                    // P_c/S_c depend only on rho_eff; φ_c depends only on rho_true.
                    // The two solves are independent: a rho_true failure must not
                    // close a channel that rho_eff confirmed is open.
                    // Reference: SAMMY rml/mrml07.f Pgh — Pghcou (rho_eff),
                    //   then Sinsix/Pf (rho_true).
                    match coulomb::coulomb_wave_functions(ch.l, eta, rho_eff) {
                        Some((f, g, fp, gp)) => {
                            // rho_eff succeeded: channel is genuinely open.
                            let fg_sq = f * f + g * g;
                            ws.p_c[c] = rho_eff / fg_sq;
                            // SHF=1: Coulomb shift ρ(F·F'+G·G')/(F²+G²).
                            // SHF=0: S_c = B_c so (S_c − B_c) = 0 in level matrix.
                            // Note: parser rejects Coulomb + SHF=1, so this arm is
                            // only reachable if that validation is later relaxed.
                            ws.s_c[c] = if pp.shf == 1 {
                                rho_eff * (f * fp + g * gp) / fg_sq
                            } else {
                                ch.boundary
                            };
                            // φ_c from rho_true; if rho_true ≤ acch, default to 0
                            // (hard-sphere limit φ → 0 as ρ → 0) without closing
                            // the channel.
                            ws.phi_c[c] = coulomb::coulomb_wave_functions(ch.l, eta, rho_true)
                                .map_or(0.0, |(fl_t, gl_t, _, _)| fl_t.atan2(gl_t));
                        }
                        None => {
                            // rho_eff ≤ acch (≈ 1e-8, SAMMY Coulfg threshold):
                            // penetrability → 0 at threshold; treat as closed.
                            // Reference: SAMMY coulomb/mrml08.f90 Coulfg — acch.
                            ws.p_c[c] = 0.0;
                            ws.s_c[c] = ch.boundary;
                            ws.phi_c[c] = 0.0;
                            ws.is_closed[c] = true;
                        }
                    }
                } else {
                    // Hard-sphere (Blatt-Weisskopf) channel.
                    ws.p_c[c] = penetrability::penetrability(ch.l, rho_eff);
                    // SHF=0: shift factor not calculated; S_c = B_c so (S_c - B_c) = 0
                    // in the level matrix diagonal.
                    // SHF=1: calculate S_c analytically (Blatt-Weisskopf).
                    // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml07.f Pgh (Ishift check)
                    ws.s_c[c] = if pp.shf == 1 {
                        penetrability::shift_factor(ch.l, rho_eff)
                    } else {
                        ch.boundary
                    };
                    ws.phi_c[c] = penetrability::phase_shift(ch.l, rho_true);
                }
            }
        }
    }

    // ── R-matrix (complex for KRM=3, real for KRM=2) ─────────────────────────
    // KRM=2 (standard R-matrix):
    //   R_cc'(E) = Σ_n γ_nc · γ_nc' / (E_n - E)   [real, reduced amplitude widths]
    //
    // KRM=3 (Reich-Moore approximation):
    //   R_cc'(E) = Σ_n γ_nc · γ_nc' / (Ẽ_n - E)   [complex, Ẽ_n = E_n - i·Γ_γn/2]
    //   where γ_nc = √(Γ_nc / (2·P_c(E_n))) (partial width → reduced amplitude).
    //   The imaginary shift makes capture implicit — |U| < 1, with missing flux
    //   going to capture.
    //
    // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml07.f Setr subroutine
    // ws.r_cplx is already zeroed by resize_and_clear.
    for res in &sg.resonances {
        let e_tilde = if krm == 3 {
            // KRM=3: convert formal partial widths to reduced amplitudes.
            // γ_nc = √(|Γ_nc| / (2·P_c(E_n))).  Sign preserved from Γ_nc.
            // For closed channels or P=0 (e.g. bound states at E_n<0): use
            // γ_nc = √(|Γ_nc|) directly (SAMMY convention for bound states).
            // Complex energy Ẽ_n = E_n - i·Γ_γ/2.
            // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml01.f
            let e_tilde = Complex64::new(res.energy, -res.gamma_gamma / 2.0);
            // P_c must be evaluated at the resonance energy E_n, not at the
            // incident energy E.  γ_nc = √(Γ_nc / (2·P_c(E_n))) is a property
            // of the resonance and must be energy-independent.  Using P_c(E)
            // would make γ_nc depend on the evaluation point, distorting the
            // resonance shape away from the peak.
            // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml01.f (reads Γ_nc then
            // converts via P at resonance energy in Pgh subroutine).
            for c in 0..nch {
                let gamma_formal = res.widths[c];
                let ch = &sg.channels[c];
                let pp_c = &particle_pairs[ch.particle_pair_idx];
                // P_c at the resonance energy E_n.  Returns None when the channel
                // is closed at E_n (bound-state resonance), Some(P) otherwise.
                //
                // The fallback must be gated on whether e_cm_n ≤ 0, NOT on
                // whether P is numerically small.  For genuinely open channels
                // near threshold (high-l, small ρ), P is positive but tiny;
                // a magnitude guard would replace √(|Γ|/(2P)) ≫ 1 with √|Γ| ≪ 1,
                // underestimating the reduced amplitude by orders of magnitude.
                // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml01.f Pgh.
                let p_at_en: Option<f64> = if pp_c.ma < 0.5 {
                    Some(1.0) // photon: P = 1 by convention
                } else {
                    // P1: res.energy is lab-frame; convert to CM before adding Q.
                    // Reference: SAMMY rml/mrml03.f Fxradi; ENDF-6 §2.2.1.6.
                    let e_cm_n = channel::lab_to_cm_energy(res.energy, awr) + pp_c.q;
                    if e_cm_n <= 0.0 {
                        None // channel closed at resonance energy (bound state)
                    } else {
                        let redmas = pp_c.ma * pp_c.mb / (pp_c.ma + pp_c.mb);
                        let k_cn = channel::wave_number_from_cm(e_cm_n, redmas);
                        let rho_eff_n = k_cn * ch.effective_radius;
                        // Must use the same penetrability type as the open-channel
                        // block: Coulomb P_c(E_n) for charged pairs, hard-sphere
                        // otherwise.  Mixing them produces inconsistent γ_nc
                        // normalisation: γ_nc = √(Γ_nc / (2·P_c(E_n))).
                        // SAMMY rml/mrml07.f Pgh — same Zeta check applies here.
                        let p = if pp_c.za.abs() > 0.5 && pp_c.zb.abs() > 0.5 {
                            // If rho_eff_n ≤ acch (SAMMY Coulfg threshold,
                            // ≈ 1e-8) wave functions return None and P = 0.0
                            // is used.  The filter(p > 0.0) below then maps
                            // P = 0 → None, applying the closed-channel
                            // √(|Γ|) normalisation — correct physical limit
                            // for a resonance barely above its Coulomb
                            // channel threshold.
                            // Reference: SAMMY coulomb/mrml08.f90 Coulfg.
                            let eta =
                                coulomb::sommerfeld_eta(pp_c.za, pp_c.zb, pp_c.ma, pp_c.mb, e_cm_n);
                            coulomb::coulomb_wave_functions(ch.l, eta, rho_eff_n)
                                .map_or(0.0, |(fl, gl, _, _)| rho_eff_n / (fl * fl + gl * gl))
                        } else {
                            penetrability::penetrability(ch.l, rho_eff_n)
                        };
                        Some(p)
                    }
                };
                // Guard: prevent division by zero if P_c = 0.  Covers both
                // non-Coulomb channels at rho → 0 and Coulomb channels below
                // SAMMY's acch threshold (coulomb_wave_functions returns None
                // → P = 0.0).  Collapsing to None uses the closed-channel
                // √(|Γ|) normalisation, which is the correct limit when the
                // channel has effectively zero penetrability.
                let p_at_en = p_at_en.filter(|&p| p > 0.0);
                ws.gamma_vals[c] = match p_at_en {
                    None => {
                        // Closed or non-computable channel at E_n: formal width used
                        // directly as reduced amplitude (SAMMY bound-state convention).
                        gamma_formal.abs().sqrt().copysign(gamma_formal)
                    }
                    Some(p) => {
                        // Open channel: γ = √(|Γ| / (2·P_c(E_n))) with sign of Γ.
                        let magnitude = (gamma_formal.abs() / (2.0 * p)).sqrt();
                        magnitude.copysign(gamma_formal)
                    }
                };
            }
            e_tilde
        } else {
            // KRM=2: widths are already reduced amplitudes; real denominator.
            // P2: Guard only against exact IEEE 754 zero; complex infrastructure
            // handles the Lorentzian width naturally via i·P_c in level matrix.
            ws.gamma_vals[..nch].copy_from_slice(&res.widths[..nch]);
            Complex64::new(res.energy, 0.0)
        };

        let denom = e_tilde - energy_ev;
        // Near-pole regularization: add a tiny imaginary offset ε to the denominator
        // so that evaluating exactly at E = E_n (where denom → 0 for real KRM=2 poles)
        // gives a finite, physically meaningful result via the Cauchy principal value.
        // For KRM=3, e_tilde already carries −iΓ_γ/2 so the denominator is never zero;
        // the correction is negligible (ε << Γ_γ/2).
        // Reference: Cauchy PV regularisation; SAMMY avoids the exact pole by perturbing
        // the resonance energy during input processing.
        let inv_denom = if denom.norm() < QUANTUM_NUMBER_EPS {
            (denom + Complex64::new(0.0, QUANTUM_NUMBER_EPS)).inv()
        } else {
            denom.inv()
        };
        for c in 0..nch {
            let gc = ws.gamma_vals[c];
            for cp in 0..nch {
                ws.r_cplx[c * nch + cp] += gc * ws.gamma_vals[cp] * inv_denom;
            }
        }
    }

    // ── L_c = (S_c - B_c) + i·P_c (per-channel level denominator) ───────────
    // Reference: SAMMY rml/mrml07.f Pgh subroutine, "PH = 1/(S-B+IP)"
    for c in 0..nch {
        ws.l_c[c] = Complex64::new(ws.s_c[c] - sg.channels[c].boundary, ws.p_c[c]);
    }

    // ── Reduced level matrix Ỹ = L⁻¹ - R (SAMMY "Ymat") ─────────────────────
    // Ỹ_cc'(E) = (1/L_c)·δ_cc' - R_cc'
    // Reference: SAMMY rml/mrml07.f — "Ymat = (1/(S-B+IP) - Rmat)"
    //
    // NOTE: This is NOT (L - R). The SAMMY formulation inverts L⁻¹ - R, not
    // L - R. Using L - R gives |U| = 3 for the hard sphere (R=0, A=iP,
    // A⁻¹·P = -i/P·P = -i, W = 1+2i²·(−1)=3) — catastrophically wrong.
    // Using L⁻¹ - R gives |U| = 1 for R=0 (Ỹ = 1/L, Ỹinv = L, XQ = L·0 = 0,
    // XXXX = 0, W = 1, U = exp(2iφ)) — correct hard-sphere limit.
    for c in 0..nch {
        // L_c = (S_c − B_c) + i·P_c.
        // For SHF=0 closed channels: S_c = B_c and P_c = 0 ⇒ L_c = 0.
        // Correct limit: 1/L_c → ∞ ⇒ Ỹ[c,c] >> R[c,c] ⇒ Ỹ⁻¹[c,c] ≈ 0
        // ⇒ channel decouples from U.  Setting 1/L_c = 0 (old bug) removes
        // the diagonal and lets R dominate — wrong coupling / Ỹ singular.
        //
        // For SHF=1 or non-matching B_c, L_c is generally finite even when
        // P_c = 0; the dispersive (real) shift must be preserved.  Do NOT
        // force the sentinel just because the channel is sub-threshold; check
        // whether |L_c| is actually near zero.
        //
        // Reference: SAMMY rml/mrml07.f — PH = 1/(S−B+iP).
        let inv_l = if ws.l_c[c].norm_sqr() < NEAR_ZERO_FLOOR {
            // |L_c|² < NEAR_ZERO_FLOOR: use finite-but-large sentinel so the diagonal
            // dominates and the channel decouples without overflow in inversion.
            Complex64::new(1e30, 0.0)
        } else {
            Complex64::new(1.0, 0.0) / ws.l_c[c]
        };
        for cp in 0..nch {
            let diag = if c == cp { inv_l } else { Complex64::ZERO };
            ws.y_tilde[c * nch + cp] = diag - ws.r_cplx[c * nch + cp];
        }
    }

    // ── Invert Ỹ to get Ỹinv (SAMMY "Yinv") ─────────────────────────────────
    // Reference: SAMMY rml/mrml09.f Yinvrs subroutine
    if !invert_complex_matrix_flat(
        &ws.y_tilde,
        nch,
        &mut ws.y_inv,
        &mut ws.aug,
        &mut ws.aug_tmp,
    ) {
        // Singular Ỹ matrix — regularize by adding a small real epsilon to
        // the diagonal and retry.  This can happen when channels are
        // near-degenerate (e.g. two channels at the same threshold).
        // The epsilon is real-only to preserve Hermitian symmetry of the
        // level matrix; an imaginary perturbation would break unitarity.
        //
        // Use a *relative* epsilon: ε = |diag| × QUANTUM_NUMBER_EPS, with a floor of
        // NEAR_ZERO_FLOOR for zero diagonals.  A fixed absolute epsilon
        // could be comparable to or larger than the diagonal value itself
        // for high-L channels with very small penetrabilities (where
        // 1/L_c ~ 1/P_c can be enormous, but R_cc' is also large, making
        // the net diagonal small).  The relative approach perturbs the
        // matrix by a fraction of its natural scale.
        //
        // Per-diagonal (not matrix-norm) regularization is intentional:
        // each channel's diagonal element lives on its own physical scale
        // (set by 1/L_c − R_cc), which can differ by orders of magnitude
        // across channels (e.g. an s-wave elastic channel vs. a high-L
        // fission channel).  A single matrix-norm epsilon would be
        // dominated by the largest channel and could either over-perturb
        // small channels or under-perturb large ones.  Per-diagonal
        // epsilon ensures each channel is nudged proportionally to its
        // own scale.

        // Copy y_tilde into y_inv as a temp buffer for the regularized matrix.
        ws.y_inv.copy_from_slice(&ws.y_tilde);
        for i in 0..nch {
            let diag_norm = ws.y_inv[i * nch + i].norm();
            // Relative regularization (QUANTUM_NUMBER_EPS × diagonal) with an
            // absolute floor of NEAR_ZERO_FLOOR for near-zero diagonals.
            let eps = (diag_norm * QUANTUM_NUMBER_EPS).max(NEAR_ZERO_FLOOR);
            ws.y_inv[i * nch + i] += Complex64::new(eps, 0.0);
        }
        // y_inv now holds the regularized y_tilde; copy it to y_tilde so the
        // inversion reads from the regularized version.
        ws.y_tilde.copy_from_slice(&ws.y_inv);
        if !invert_complex_matrix_flat(
            &ws.y_tilde,
            nch,
            &mut ws.y_inv,
            &mut ws.aug,
            &mut ws.aug_tmp,
        ) {
            return (0.0, 0.0, 0.0, 0.0); // truly degenerate
        }
    }

    // ── XQ = Ỹinv · R (matrix product, SAMMY "Xqr/Xqi") ─────────────────────
    // Reference: SAMMY rml/mrml11.f Setxqx — "Xqr(k,i) = (L**-1-R)**-1 * R"
    for c in 0..nch {
        for cp in 0..nch {
            let mut sum = Complex64::ZERO;
            for k in 0..nch {
                sum += ws.y_inv[c * nch + k] * ws.r_cplx[k * nch + cp];
            }
            ws.xq[c * nch + cp] = sum;
        }
    }

    // ── XXXX = (√P_c / L_c) · XQ · √P_c' ────────────────────────────────────
    // Reference: SAMMY rml/mrml11.f Setxqx — "Xxxx = sqrt(P)/L * xq * sqrt(P)"
    for c in 0..nch {
        ws.sqrt_p[c] = ws.p_c[c].sqrt();
    }
    for c in 0..nch {
        // For a closed channel: sqrt_p[c] = 0 and L_c = 0 (0/0 indeterminate).
        // The full XXXX[c,cp] = (√P_c / L_c) · XQ[c,cp] · √P_c'.
        // Since √P_c = 0 for any closed channel c, the entire row is zero
        // regardless of the value of L_c. Setting sqrt_p_over_l = 0 is correct.
        // (The Ỹ sentinel handles Ỹ inversion correctly; this row zeroing is
        //  consistent: a closed channel contributes nothing to XXXX/U.)
        // Guard: at exact channel threshold, both √P_c and L_c can be
        // zero simultaneously, producing 0/0 = NaN.  The is_closed flag
        // catches most cases, but a channel right at threshold might not
        // be flagged closed yet have |L_c| ≈ 0.  The extra norm check
        // prevents NaN propagation.
        let sqrt_p_over_l = if ws.is_closed[c] || ws.l_c[c].norm() < PIVOT_FLOOR {
            Complex64::ZERO
        } else {
            ws.sqrt_p[c] / ws.l_c[c]
        };
        for cp in 0..nch {
            ws.xxxx[c * nch + cp] = sqrt_p_over_l * ws.xq[c * nch + cp] * ws.sqrt_p[cp];
        }
    }

    // ── Collision matrix U = Ω · W · Ω, W = I + 2i·Ξ ────────────────────────
    //
    // Phase convention: Ω_c = exp(-iφ_c), NOT exp(+iφ_c).
    //
    // TRUTH SOURCE: SAMMY rml/mrml11.f lines 14-18:
    //   "W(c,c') = delta(c,c') + 2i XXXX(c,c')" (eq. III.D.4)
    // And mrml11.f lines 84-88, the elastic formula:
    //   "sin²(φ)·(1-2Xi) - sin(2φ)·Xr + Xr²+Xi²"
    // is consistent ONLY with U = e^{-2iφ}·(I+2iX), not e^{+2iφ}.
    //
    // For hard sphere (W=I): |1-e^{-2iφ}|² = |1-e^{+2iφ}|² = 4sin²φ,
    // so the sign error is invisible in unitarity tests. It ONLY manifests
    // when resonances are present (W ≠ I), producing wrong interference
    // patterns in elastic and incorrect total from optical theorem.
    //
    // History: same class of bug as the MLBW e^{+2iφ} error fixed in
    // slbw.rs (commit f0eadc1). The negative exponent is the ENDF/SAMMY
    // convention; the positive exponent is a common error.
    for c in 0..nch {
        ws.omega[c] = Complex64::from_polar(1.0, -ws.phi_c[c]);
    }

    for c in 0..nch {
        for cp in 0..nch {
            let delta = if c == cp {
                Complex64::new(1.0, 0.0)
            } else {
                Complex64::ZERO
            };
            let w_cc = delta + Complex64::new(0.0, 2.0) * ws.xxxx[c * nch + cp];
            ws.u[c * nch + cp] = ws.omega[c] * w_cc * ws.omega[cp];
        }
    }

    // ── Cross-sections (sum over entrance channels) ───────────────────────────
    // Optical theorem gives σ_total = 2·(π/k²)·g_J·(1 - Re(U_cc)) per channel.
    // Reference: SAMMY rml/mrml11.f Sectio subroutine; SAMMY manual §3.1 Eq. 3.4
    let mut tot = 0.0;
    let mut elas = 0.0;
    let mut cap = 0.0;
    let mut fis = 0.0;
    let mut inel = 0.0; // inelastic neutron channels (MT=51+): tracked separately

    // Whether this spin group has explicit capture (photon) channels in the
    // level matrix.  KRM=2 with photon channels: yes.  KRM=3: no (capture is
    // implicit via complex poles; no MT=102 channel appears in NCH).
    let has_explicit_capture = ws.is_capture[..nch].iter().any(|&x| x);

    for c0 in 0..nch {
        if !ws.is_entrance[c0] {
            continue;
        }
        let u_diag = ws.u[c0 * nch + c0];
        // σ_total (optical theorem, per entrance channel)
        tot += 2.0 * pok2 * g_j * (1.0 - u_diag.re);
        // σ_elastic: |1 - U_{c0,c0}|²
        elas += pok2 * g_j * (Complex64::new(1.0, 0.0) - u_diag).norm_sqr();

        for cp in 0..nch {
            if ws.is_fission[cp] {
                // σ_fission: |U_{c0,c'}|² for fission channels c'
                fis += pok2 * g_j * ws.u[c0 * nch + cp].norm_sqr();
            }
            if has_explicit_capture && ws.is_capture[cp] {
                // σ_capture (explicit): |U_{c0,c'}|² for photon channels c'.
                // Avoids lumping inelastic neutron channels (MT=51+) into capture.
                // Reference: SAMMY rml/mrml11.f Sectio — explicit sum over γ channels.
                cap += pok2 * g_j * ws.u[c0 * nch + cp].norm_sqr();
            }
            if ws.is_inelastic[cp] {
                // σ_inelastic: |U_{c0,c'}|² for inelastic neutron channels (MT=51+).
                // Tracked separately so KRM=3 capture residual excludes this flux.
                // Reference: ENDF MT conventions §3.4; SAMMY rml/mrml11.f Sectio.
                inel += pok2 * g_j * ws.u[c0 * nch + cp].norm_sqr();
            }
        }
    }

    if krm == 3 && !has_explicit_capture {
        // KRM=3 (Reich-Moore approximation): capture is implicit via complex poles
        // (Ẽ_n = E_n - i·Γγ/2).  Flux not going to elastic, fission, or inelastic
        // channels is capture.  Inelastic flux must be excluded; folding it into
        // capture would mislabel σ_capture when MT=51+ channels are present.
        // Clamp to ≥0 for floating-point safety near pole energies.
        // Reference: ENDF-6 §2.2.1.6; SAMMY rml/mrml11.f Sectio.
        cap = (tot - elas - fis - inel).max(0.0);
    }
    // KRM=2: capture was accumulated explicitly above from MT=102 channels.
    // Do NOT add residual flux — it may include inelastic (MT=51+) contributions
    // and would mislabel them as capture, biasing channel-resolved fits.
    // Reference: SAMMY rml/mrml11.f Sectio (explicit γ-channel sum for KRM=2).

    (tot, elas, cap, fis)
}

// ── Complex Gauss-Jordan Elimination (flat-buffer version) ──────────────────
//
// Inverts an n×n complex matrix using Gauss-Jordan elimination with partial
// pivoting. Returns false if the matrix is singular (pivot magnitude < LOG_FLOOR).
//
// For LRF=7 isotopes relevant to VENUS imaging, NCH ≤ 6, so O(n³) is fast.
// SAMMY uses a specialized complex symmetric factorization (Xspfa/Xspsl in
// rml/mrml10.f), but Gauss-Jordan is correct and sufficient for our purposes.
//
// All buffers are caller-provided to avoid per-call allocation:
// - `a`: input matrix (flat row-major, n×n)
// - `out`: output inverse (flat row-major, n×n)
// - `aug`: augmented matrix workspace (flat row-major, n×2n)
// - `tmp`: temporary row buffer (length 2n)

fn invert_complex_matrix_flat(
    a: &[Complex64],
    n: usize,
    out: &mut [Complex64],
    aug: &mut [Complex64],
    tmp: &mut [Complex64],
) -> bool {
    let w = 2 * n; // augmented row width

    // Build augmented matrix [A | I] of size n × 2n
    for r in 0..n {
        for c in 0..n {
            aug[r * w + c] = a[r * n + c];
        }
        for c in 0..n {
            aug[r * w + n + c] = if c == r {
                Complex64::new(1.0, 0.0)
            } else {
                Complex64::ZERO
            };
        }
    }

    for col in 0..n {
        // Partial pivoting: find row with largest magnitude in this column
        let mut best = col;
        let mut best_norm = aug[col * w + col].norm();
        for r in (col + 1)..n {
            let norm = aug[r * w + col].norm();
            if norm > best_norm {
                best = r;
                best_norm = norm;
            }
        }
        // Swap rows col and best in aug
        if best != col {
            for j in 0..w {
                aug.swap(col * w + j, best * w + j);
            }
        }

        let pivot = aug[col * w + col];
        if pivot.norm() < LOG_FLOOR {
            return false; // singular
        }

        // Scale pivot row so leading entry becomes 1
        let inv_pivot = pivot.inv();
        for j in 0..w {
            aug[col * w + j] *= inv_pivot;
        }

        // Eliminate this column from all other rows.
        // Copy pivot row to tmp to avoid aliasing issues.
        tmp[..w].copy_from_slice(&aug[col * w..col * w + w]);
        for row in 0..n {
            if row == col {
                continue;
            }
            let factor = aug[row * w + col];
            if factor.norm() < LOG_FLOOR {
                continue;
            }
            for j in 0..w {
                aug[row * w + j] -= factor * tmp[j];
            }
        }
    }

    // Extract the right half (the inverse)
    for r in 0..n {
        for c in 0..n {
            out[r * n + c] = aug[r * w + n + c];
        }
    }
    true
}

// ── Legacy wrapper for tests (allocating version) ───────────────────────────
#[cfg(test)]
fn invert_complex_matrix(a: &[Vec<Complex64>], n: usize) -> Option<Vec<Vec<Complex64>>> {
    let flat_a: Vec<Complex64> = a.iter().flat_map(|row| row.iter().copied()).collect();
    let mut flat_out = vec![Complex64::ZERO; n * n];
    let mut aug = vec![Complex64::ZERO; n * 2 * n];
    let mut tmp = vec![Complex64::ZERO; 2 * n];
    if invert_complex_matrix_flat(&flat_a, n, &mut flat_out, &mut aug, &mut tmp) {
        Some(
            (0..n)
                .map(|r| flat_out[r * n..(r + 1) * n].to_vec())
                .collect(),
        )
    } else {
        None
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_identity_inversion() {
        let n = 3;
        let a: Vec<Vec<Complex64>> = (0..n)
            .map(|r| {
                (0..n)
                    .map(|c| {
                        if r == c {
                            Complex64::new(1.0, 0.0)
                        } else {
                            Complex64::ZERO
                        }
                    })
                    .collect()
            })
            .collect();
        let inv = invert_complex_matrix(&a, n).unwrap();
        for (r, row) in inv.iter().enumerate() {
            for (c, val) in row.iter().enumerate() {
                let expected = if r == c { 1.0 } else { 0.0 };
                assert!(
                    (val.re - expected).abs() < 1e-12,
                    "inv[{r}][{c}].re = {}, expected {expected}",
                    val.re
                );
                assert!(
                    val.im.abs() < 1e-12,
                    "inv[{r}][{c}].im = {} should be 0",
                    val.im
                );
            }
        }
    }

    #[test]
    fn test_2x2_complex_inversion() {
        // A = [[2+i, 1], [0, 3-2i]]  → A⁻¹ = [[1/(2+i), -1/((2+i)(3-2i))], [0, 1/(3-2i)]]
        let a00 = Complex64::new(2.0, 1.0);
        let a01 = Complex64::new(1.0, 0.0);
        let a11 = Complex64::new(3.0, -2.0);
        let a = vec![vec![a00, a01], vec![Complex64::ZERO, a11]];
        let inv = invert_complex_matrix(&a, 2).unwrap();

        // Verify A · A⁻¹ ≈ I
        let i00 = a00 * inv[0][0] + a01 * inv[1][0];
        let i01 = a00 * inv[0][1] + a01 * inv[1][1];
        let i11 = a11 * inv[1][1];
        assert!((i00.re - 1.0).abs() < 1e-12, "i00.re = {}", i00.re);
        assert!(i00.im.abs() < 1e-12, "i00.im = {}", i00.im);
        assert!(i01.norm() < 1e-12, "i01 = {}", i01);
        assert!((i11.re - 1.0).abs() < 1e-12, "i11.re = {}", i11.re);
    }

    #[test]
    fn test_singular_returns_none() {
        let a = vec![
            vec![Complex64::new(1.0, 0.0), Complex64::new(2.0, 0.0)],
            vec![Complex64::new(2.0, 0.0), Complex64::new(4.0, 0.0)],
        ];
        assert!(invert_complex_matrix(&a, 2).is_none());
    }

    /// Hard-sphere unitarity check: with R = 0, U must equal exp(2iφ)·I
    /// and σ_total = 2·(π/k²)·g_J·(1 − cos 2φ) ≥ 0.
    ///
    /// This test is purely local (no network) and guards against the
    /// classic sign error where U = 3·exp(2iφ)·I (|U| = 3) that arises
    /// when using A = L − R instead of SAMMY's Ỹ = L⁻¹ − R.
    #[test]
    fn test_hard_sphere_unitarity() {
        use nereids_endf::resonance::{ParticlePair, RmlChannel, RmlData, SpinGroup};

        // Minimal synthetic LRF=7 / KRM=2 with a single elastic channel
        // and NO resonances.  Result must be pure hard-sphere scattering.
        let pp = ParticlePair {
            ma: 1.0,
            mb: 184.0,
            za: 0.0,  // neutron charge Z=0
            zb: 74.0, // W-184 charge Z=74 (ENDF LRF=7 stores charge directly)
            ia: 0.5,
            ib: 0.0,
            q: 0.0,
            pnt: 1,
            shf: 0, // SHF=0 → S_c = B_c → L_c = iP_c
            mt: 2,
            pa: 1.0,
            pb: 1.0,
        };
        let channel = RmlChannel {
            particle_pair_idx: 0,
            l: 0,
            channel_spin: 0.5,
            boundary: 0.0,
            effective_radius: 8.3,
            true_radius: 8.3,
        };
        let sg = SpinGroup {
            j: 0.5,
            parity: 1.0,
            channels: vec![channel],
            resonances: vec![], // no resonances: pure hard sphere
            has_background_correction: false,
        };
        let rml = RmlData {
            target_spin: 0.0,
            awr: 183.0,
            scattering_radius: 8.3,
            krm: 2,
            particle_pairs: vec![pp],
            spin_groups: vec![sg],
        };

        // Evaluate at several energies.  σ_total must be non-negative,
        // and σ_capture/σ_fission must be zero (no absorption channels).
        //
        // Note: cap is computed as (tot - elas - fis); since tot and elas are
        // calculated via different floating-point paths they may differ by
        // ~1e-15 × pok2 (where pok2 can be ~1e5 b at low energy), giving a
        // residual |cap| ~ 1e-10 b.  Use a relative tolerance of 1e-9.
        for &e_ev in &[10.0, 50.0, 100.0, 500.0, 1000.0] {
            let (tot, elas, cap, fis) = cross_sections_for_rml_range(&rml, e_ev);
            assert!(
                tot >= 0.0,
                "hard sphere σ_total < 0 at {e_ev} eV: {tot:.6} b"
            );
            let tol = 1e-9 * tot.abs().max(1.0);
            assert!(
                cap.abs() < tol,
                "hard sphere σ_capture ≠ 0 at {e_ev} eV: {cap:.3e} b (tol={tol:.3e})"
            );
            assert!(
                fis.abs() < tol,
                "hard sphere σ_fission ≠ 0 at {e_ev} eV: {fis:.3e} b (tol={tol:.3e})"
            );
            // σ_elastic ≈ σ_total (capture=fission=0)
            assert!(
                (tot - elas).abs() < tol,
                "σ_total ≠ σ_elastic at {e_ev} eV: tot={tot:.6}, elas={elas:.6}"
            );
        }
    }

    /// Coulomb exit channels route through coulomb::coulomb_penetrability,
    /// not the hard-sphere Blatt-Weisskopf functions.
    ///
    /// Constructs a 2-channel spin group:
    ///   ch0: neutron (za=0)  + target — hard-sphere entrance channel
    ///   ch1: α (za=2)        + O-16 (zb=8) — Coulomb exit, Q=+50 eV (always open)
    ///
    /// ENDF LRF=7 stores charge Z directly in ZA/ZB (neutron=0, alpha=2, O-16=8).
    ///
    /// Verifies σ_total ≥ 0 (physics sanity) and no panic at both an open
    /// and a closed Coulomb channel (Q very negative).
    ///
    /// SAMMY ref: rml/mrml07.f Pgh — `if (Zeta(I).NE.Zero)` branch.
    #[test]
    fn test_coulomb_channel_open_and_closed_no_panic() {
        use nereids_endf::resonance::{ParticlePair, RmlChannel, RmlData, SpinGroup};

        // Entrance channel: neutron (Z=0) + W-184 target (Z=74).
        // ENDF LRF=7 stores charge directly: neutron za=0, W-184 zb=74.
        let pp_entrance = ParticlePair {
            ma: 1.0,
            mb: 184.0,
            za: 0.0,  // neutron charge Z=0
            zb: 74.0, // W-184 charge Z=74
            ia: 0.5,
            ib: 0.0,
            q: 0.0,
            pnt: 1,
            shf: 0,
            mt: 2,
            pa: 1.0,
            pb: 1.0,
        };

        // Coulomb exit channel: α(Z=2) + O-16(Z=8), Q=+50 eV → always open.
        // sommerfeld_eta(2, 8, ...) gives η > 0, confirming Coulomb branch.
        let pp_coulomb_open = ParticlePair {
            ma: 4.0,
            mb: 16.0,
            za: 2.0, // alpha charge Z=2
            zb: 8.0, // O-16 charge Z=8
            ia: 0.0,
            ib: 0.0,
            q: 50.0, // Q > 0 → e_c = e_cm + 50 > 0 for all positive energies
            pnt: 1,
            shf: 0,
            mt: 22, // (n,α)
            pa: 1.0,
            pb: 1.0,
        };

        // Coulomb exit channel with Q very negative → closed at all reasonable energies.
        let pp_coulomb_closed = ParticlePair {
            ma: 4.0,
            mb: 16.0,
            za: 2.0, // alpha charge Z=2
            zb: 8.0, // O-16 charge Z=8
            ia: 0.0,
            ib: 0.0,
            q: -1e6, // far below threshold
            pnt: 1,
            shf: 0,
            mt: 22,
            pa: 1.0,
            pb: 1.0,
        };

        // Build and evaluate the open-channel case.
        for (desc, pp_exit, expect_positive_total) in [
            ("open Coulomb exit", &pp_coulomb_open, true),
            ("closed Coulomb exit", &pp_coulomb_closed, false),
        ] {
            let ch0 = RmlChannel {
                particle_pair_idx: 0,
                l: 0,
                channel_spin: 0.5,
                boundary: 0.0,
                effective_radius: 8.3,
                true_radius: 8.3,
            };
            let ch1 = RmlChannel {
                particle_pair_idx: 1,
                l: 0,
                channel_spin: 0.5,
                boundary: 0.0,
                effective_radius: 5.0,
                true_radius: 5.0,
            };
            let sg = SpinGroup {
                j: 0.5,
                parity: 1.0,
                channels: vec![ch0, ch1],
                resonances: vec![],
                has_background_correction: false,
            };
            let rml = RmlData {
                target_spin: 0.0,
                awr: 183.0,
                scattering_radius: 8.3,
                krm: 2,
                particle_pairs: vec![pp_entrance.clone(), pp_exit.clone()],
                spin_groups: vec![sg],
            };

            let (tot, _elas, _cap, _fis) = cross_sections_for_rml_range(&rml, 100.0);
            assert!(tot >= 0.0, "{desc}: σ_total = {tot:.6} b must be ≥ 0");
            if expect_positive_total {
                // Hard-sphere entrance channel alone gives positive σ_total
                // (the Coulomb channel merely adds a second channel but no resonances).
                assert!(
                    tot > 0.0,
                    "{desc}: σ_total = {tot} b should be > 0 (hard-sphere entrance channel)"
                );
            }
        }
    }

    /// Singular Y-matrix regularization: construct a KRM=2 spin group where
    /// the R-matrix nearly cancels L⁻¹ at the resonance energy, making
    /// the Y-matrix nearly singular.  The relative-epsilon regularization
    /// must produce finite, non-negative cross-sections (no NaN/Inf).
    #[test]
    fn test_rml_singular_y_matrix_regularization() {
        use nereids_endf::resonance::{ParticlePair, RmlChannel, RmlData, RmlResonance, SpinGroup};

        let pp = ParticlePair {
            ma: 1.0,
            mb: 184.0,
            za: 0.0,
            zb: 74.0,
            ia: 0.5,
            ib: 0.0,
            q: 0.0,
            pnt: 1,
            shf: 0,
            mt: 2,
            pa: 1.0,
            pb: 1.0,
        };
        let ch = RmlChannel {
            particle_pair_idx: 0,
            l: 0,
            channel_spin: 0.5,
            boundary: 0.0,
            effective_radius: 8.3,
            true_radius: 8.3,
        };
        // KRM=2: widths are reduced amplitudes γ_nc.
        // A very large reduced amplitude γ at E_r makes R ≈ γ²/(E_r - E)
        // huge near E_r, potentially cancelling L⁻¹ on the diagonal.
        let sg = SpinGroup {
            j: 0.5,
            parity: 1.0,
            channels: vec![ch],
            resonances: vec![RmlResonance {
                energy: 100.0,
                gamma_gamma: 0.0,  // KRM=2 has no implicit capture
                widths: vec![1e5], // very large reduced amplitude
            }],
            has_background_correction: false,
        };
        let rml = RmlData {
            target_spin: 0.0,
            awr: 183.0,
            scattering_radius: 8.3,
            krm: 2,
            particle_pairs: vec![pp],
            spin_groups: vec![sg],
        };

        // Evaluate exactly at the resonance energy (worst case for singularity).
        let (tot, elas, cap, _fis) = cross_sections_for_rml_range(&rml, 100.0);
        assert!(
            tot.is_finite() && tot >= 0.0,
            "σ_total must be finite and ≥ 0 near singular Y-matrix, got {tot}"
        );
        assert!(
            elas.is_finite() && elas >= 0.0,
            "σ_elastic must be finite and ≥ 0, got {elas}"
        );
        assert!(cap.is_finite(), "σ_capture must be finite, got {cap}");
    }

    /// Verify W-184 cross-sections show resonance structure.
    ///
    /// Downloads W-184 ENDF/B-VIII.0, parses LRF=7 parameters,
    /// then checks that σ_total at the first resonance (~101.9 eV) is
    /// significantly larger than the background off-resonance.
    ///
    /// Run with: cargo test -p nereids-physics -- --ignored test_w184_cross_section_peak
    #[test]
    #[ignore = "requires network: downloads W-184 ENDF from IAEA (~50 kB)"]
    fn test_w184_cross_section_peak() {
        use crate::reich_moore::cross_sections_at_energy;
        use nereids_core::types::Isotope;
        use nereids_endf::parser::parse_endf_file2;
        use nereids_endf::retrieval::{EndfLibrary, EndfRetriever};

        let retriever = EndfRetriever::new();
        let isotope = Isotope::new(74, 184).unwrap();
        let (_, text) = retriever
            .get_endf_file(&isotope, EndfLibrary::EndfB8_0, 7437)
            .expect("Failed to download W-184 ENDF/B-VIII.0");

        let data = parse_endf_file2(&text).expect("Failed to parse W-184 ENDF");

        // W-184 ENDF/B-VIII.0 (KRM=3, 3 spin groups J=1/2+, J=1/2-, J=3/2-):
        //   First positive resonance in J=1/2+ spin group: ~101.9 eV
        //   Background between the -386 eV bound state and 101.9 eV resonance: ~0.07 b
        // Test that σ_total at the 101.9 eV resonance peak >> background at 50 eV.
        let xs_on_res = cross_sections_at_energy(&data, 101.9);
        let xs_off_res = cross_sections_at_energy(&data, 50.0);

        // Background must be non-negative (guards against the A=L-R sign error
        // that gave σ(50 eV) = −0.228 b).
        assert!(
            xs_off_res.total >= 0.0,
            "σ_total at 50 eV must be non-negative (hard-sphere background), got {:.4} b",
            xs_off_res.total
        );
        assert!(
            xs_on_res.total > 0.0,
            "σ_total at 101.9 eV should be positive, got {}",
            xs_on_res.total
        );
        assert!(
            xs_on_res.capture >= 0.0,
            "σ_capture at 101.9 eV should be non-negative, got {}",
            xs_on_res.capture
        );
        assert!(
            xs_on_res.total > xs_off_res.total * 5.0,
            "Resonance peak at 101.9 eV should be >5× the 50 eV background: \
             σ(101.9)={:.3} vs σ(50.0)={:.3} barns",
            xs_on_res.total,
            xs_off_res.total
        );

        println!(
            "W-184 σ_total: {:.3} b at 101.9 eV (resonance), {:.3} b at 50.0 eV (background)",
            xs_on_res.total, xs_off_res.total
        );
        println!(
            "W-184 σ_capture: {:.4} b, σ_elastic: {:.4} b at 101.9 eV",
            xs_on_res.capture, xs_on_res.elastic
        );
    }
}