nereids_physics/

channel.rs

//! Neutron channel calculations: wave number, ρ parameter, statistical weights.
//!
//! Computes channel-level quantities needed for R-matrix cross-section
//! evaluation: the neutron wave number k(E), channel parameter ρ = k·a,
//! and spin statistical weight factors g_J.
//!
//! ## SAMMY Reference
//! - `rml/mrml03.f` Fxradi: wave number and radius computations
//! - SAMMY manual Section II.A (scattering-theory equations and kinematics)
//!
//! ## Units
//! All quantities in natural units: energies in eV, lengths in fm.
//!
//! ## Key Relations
//! - k = √(2·μ·E_cm) / ℏc = [AWR/(1+AWR)] × √(2·m_n·E_lab) / ℏc
//! - ρ = k·a, where a = channel radius
//! - g_J = (2J+1) / ((2I+1)·(2s+1))

use nereids_core::constants;

/// Neutron wave number squared k² in fm⁻².
///
/// k² = 2·μ·E_cm / (ℏc)² = 2·m_n·[AWR/(1+AWR)]²·E_lab / (ℏc)²
///
/// where μ = m_n · AWR / (1 + AWR) is the reduced mass,
/// E_cm = E_lab · AWR / (1 + AWR) is the center-of-mass energy,
/// and the two factors of AWR/(1+AWR) come from μ and E_cm respectively.
///
/// ## SAMMY Reference
/// - `rml/mrml03.f` Fxradi: Zke = Twomhb × √(Redmas × Factor)
///   where Redmas = AWR/(1+AWR) and Factor = AWR/(1+AWR).
///
/// # Arguments
/// * `energy_ev` — Neutron energy in eV (lab frame, as stored in ENDF).
/// * `awr` — Ratio of target mass to neutron mass (from ENDF).
pub fn k_squared(energy_ev: f64, awr: f64) -> f64 {
    // Guard: bound-state resonances have negative E_r in ENDF, but k² is
    // evaluated at the lab energy which is always positive.  If a caller
    // passes E ≤ 0 (e.g. during extrapolation or misuse), return 0.0 to
    // prevent a negative k² from propagating NaN through sqrt downstream.
    if energy_ev <= 0.0 {
        return 0.0;
    }

    // μ/m_n = AWR / (1 + AWR)
    let mass_ratio = awr / (1.0 + awr);

    // k² = 2 · m_n · (AWR/(1+AWR))² · E_lab / (ℏc)²
    // One factor from reduced mass, one from lab→CM energy conversion.
    let mn_ev = constants::NEUTRON_MASS_MEV * 1e6; // MeV → eV
    let hbar_c = constants::HBAR_EV_S * constants::SPEED_OF_LIGHT * 1e15; // eV·fm

    2.0 * mn_ev * mass_ratio * mass_ratio * energy_ev / (hbar_c * hbar_c)
}

/// Neutron wave number k in fm⁻¹.
///
/// # Arguments
/// * `energy_ev` — Neutron energy in eV (must be positive).
/// * `awr` — Ratio of target mass to neutron mass.
pub fn wave_number(energy_ev: f64, awr: f64) -> f64 {
    k_squared(energy_ev, awr).sqrt()
}

/// Wave number k_c from reduced mass ratio and center-of-mass kinetic energy.
///
/// k_c = √(2·m_n·μ·E_cm) / ℏc
///
/// Used for LRF=7 (R-Matrix Limited) channels, whose particle pairs have
/// their own masses and Q-values.
///
/// For the elastic channel (MA=1 neutron, MB=AWR target):
///   μ = AWR/(1+AWR),  E_cm = E_lab·AWR/(1+AWR)
///   → k = `wave_number(E_lab, AWR)` (algebraically identical).
///
/// ## SAMMY Reference
/// - `rml/mrml03.f` Fxradi: `Zke = Twomhb × √(Redmas × Factor)`
///   where Redmas = MA·MB/(MA+MB) and Factor = E_cm for the channel.
///
/// # Arguments
/// * `e_cm` — Center-of-mass kinetic energy in this channel (eV). Must be ≥ 0.
/// * `reduced_mass_ratio` — μ = MA·MB/(MA+MB) in neutron mass units.
pub fn wave_number_from_cm(e_cm: f64, reduced_mass_ratio: f64) -> f64 {
    // Guard against non-positive CM energies, which can arise from numerical noise
    // near threshold energies. Zero wave number means no contribution from this channel.
    // Uses `<= 0.0` to match `k_squared` (line 41), which also catches zero.
    if e_cm <= 0.0 {
        return 0.0;
    }
    let mn_ev = constants::NEUTRON_MASS_MEV * 1e6; // MeV → eV
    let hbar_c = constants::HBAR_EV_S * constants::SPEED_OF_LIGHT * 1e15; // eV·fm
    (2.0 * mn_ev * reduced_mass_ratio * e_cm / (hbar_c * hbar_c)).sqrt()
}

/// Channel parameter ρ = k·a.
///
/// # Arguments
/// * `energy_ev` — Neutron energy in eV.
/// * `awr` — Mass ratio from ENDF.
/// * `channel_radius_fm` — Channel radius in fm.
pub fn rho(energy_ev: f64, awr: f64, channel_radius_fm: f64) -> f64 {
    wave_number(energy_ev, awr) * channel_radius_fm
}

/// Spin statistical weight factor g_J.
///
///   g_J = (2|J| + 1) / ((2I + 1) · (2s + 1))
///
/// where I = target spin, s = neutron spin (1/2), J = total angular momentum.
///
/// Uses `|J|` so that the SAMMY sign convention (negative J to distinguish
/// spin groups with the same |J|) produces the correct weight.  ENDF data
/// always has J ≥ 0, so the `abs()` is a no-op for standard ENDF inputs.
///
/// # Arguments
/// * `j_total` — Total angular momentum J (may be negative per SAMMY convention).
/// * `target_spin` — Target nucleus spin I.
pub fn statistical_weight(j_total: f64, target_spin: f64) -> f64 {
    let neutron_spin = 0.5;
    (2.0 * j_total.abs() + 1.0) / ((2.0 * target_spin + 1.0) * (2.0 * neutron_spin + 1.0))
}

/// ENDF-6 §2.2.1 default channel radius (fm).
///
/// a_c = (0.123 · A^(1/3) + 0.08) × 10  \[fm\]
///
/// The ENDF formula gives values in 10⁻¹² cm; the × 10 converts to fm
/// (matching the ENDF_RADIUS_TO_FM factor applied to AP during parsing).
pub fn endf_channel_radius_fm(awr: f64) -> f64 {
    (0.123 * awr.cbrt() + 0.08) * 10.0
}

/// Convert between lab and center-of-mass energy.
///
/// E_cm = E_lab · AWR / (1 + AWR)
pub fn lab_to_cm_energy(energy_lab: f64, awr: f64) -> f64 {
    energy_lab * awr / (1.0 + awr)
}

/// π/k² factor that appears in all cross-section formulas.
///
/// Returns π/k² in barns (1 barn = 100 fm²).
///
/// # Arguments
/// * `energy_ev` — Neutron energy in eV.
/// * `awr` — Mass ratio from ENDF.
pub fn pi_over_k_squared_barns(energy_ev: f64, awr: f64) -> f64 {
    // Catch non-positive energies in debug builds: k_squared(0.0) returns 0.0,
    // and pi_over_k_squared_barns would silently map 0.0 → 1e-20 via the
    // floor, producing a finite (wrong) result instead of the expected
    // infinity/zero. Callers should never pass non-positive energy here.
    //
    // Note on debug_assert vs assert: the check `energy_ev > 0.0` fires
    // for both zero and negative energies, but only in debug builds.
    // In release builds, a non-positive energy is silently clamped to 1e-20
    // by the floor below, yielding a large but finite π/k² (~10¹⁸ barns).
    // This is acceptable because:
    //   1. All call sites guarantee positive energies (ENDF grids, fitting
    //      bounds, transmission forward model).
    //   2. The 1e-20 floor is defense-in-depth: it prevents division by zero
    //      and NaN propagation without paying for a branch in hot loops.
    //   3. A full `assert!` would add a branch to the hottest inner loop of
    //      cross-section evaluation for a condition that cannot occur in
    //      correct usage.
    debug_assert!(
        energy_ev > 0.0,
        "pi_over_k_squared_barns called with non-positive energy: {energy_ev}"
    );
    // Apply a tiny energy floor (1e-20 eV) so that k² is never zero.
    // This preserves the correct 1/E functional form at very low energies
    // instead of returning an arbitrary sentinel value that could mislead
    // a fitting optimizer.
    let e_safe = energy_ev.max(1e-20);
    let k2 = k_squared(e_safe, awr);
    // π/k² in fm², convert to barns (1 barn = 100 fm²)
    std::f64::consts::PI / k2 / 100.0
}

/// Σ g_J over the channels of orbital angular momentum `l` that no J-group
/// of `present_j` holds.
///
/// A neutron (spin ½) on a target of spin I has channel spins
/// s = |I − ½| and I + ½, and each s allows J = |l − s| … l + s, so a J can
/// belong to two channels; over every channel Σ g_J = 2l + 1. A J-group holds
/// one channel of its |J|. A channel no J-group holds scatters as a hard
/// sphere, π/k²·g_J·4 sin²φ_l.
pub(crate) fn uncovered_channel_weight(
    l: u32,
    target_spin: f64,
    present_j: impl IntoIterator<Item = f64>,
) -> f64 {
    let twice_i = (2.0 * target_spin).round() as i64;
    let twice_l = 2 * i64::from(l);
    let mut open: Vec<i64> = [twice_i - 1, twice_i + 1]
        .into_iter()
        .filter(|&twice_s| twice_s >= 0)
        .flat_map(|twice_s| ((twice_l - twice_s).abs()..=twice_l + twice_s).step_by(2))
        .collect();
    for j in present_j {
        let twice_j = (2.0 * j.abs()).round() as i64;
        if let Some(at) = open.iter().position(|&o| o == twice_j) {
            open.swap_remove(at);
        }
    }
    open.into_iter()
        .map(|twice_j| statistical_weight(twice_j as f64 / 2.0, target_spin))
        .sum()
}

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

    #[test]
    fn every_channel_together_weighs_two_l_plus_one() {
        for twice_i in 0..=9 {
            let target_spin = f64::from(twice_i) / 2.0;
            for l in 0..=3 {
                let all = uncovered_channel_weight(l, target_spin, []);
                assert!(
                    (all - f64::from(2 * l + 1)).abs() < 1e-12,
                    "I = {target_spin}, l = {l}: {all}"
                );
            }
        }
    }

    #[test]
    fn a_j_group_holds_one_channel_of_its_j() {
        // I = 1/2, l = 1: s = 0 gives J = 1; s = 1 gives J = 0, 1, 2.
        let g = |j: f64| statistical_weight(j, 0.5);
        assert_eq!(uncovered_channel_weight(1, 0.5, [1.0]), 3.0 - g(1.0));
        assert_eq!(
            uncovered_channel_weight(1, 0.5, [1.0, -1.0]),
            g(0.0) + g(2.0)
        );
        assert_eq!(uncovered_channel_weight(1, 0.5, [0.0, 1.0, -1.0, 2.0]), 0.0);
        // A J the channels do not allow holds none of them.
        assert_eq!(uncovered_channel_weight(0, 3.5, [3.0, 4.0, 5.0]), 0.0);
        assert_eq!(uncovered_channel_weight(0, 3.5, [5.0]), 1.0);
    }

    #[test]
    fn test_wave_number_u238() {
        // U-238: AWR ≈ 236.006
        // k = sqrt(2 m_n / (ℏc)²) × AWR/(1+AWR) × sqrt(E)
        // At E = 1 eV: k_free = 2.197e-4 fm⁻¹, mass_ratio = 236.006/237.006 = 0.99578
        // k = 2.197e-4 × 0.99578 ≈ 2.188e-4 fm⁻¹
        let awr = 236.006;
        let k = wave_number(1.0, awr);
        assert!(
            (k - 2.188e-4).abs() < 1e-6,
            "k = {} fm⁻¹, expected ~2.188e-4",
            k
        );
    }

    #[test]
    fn test_rho_u238() {
        // U-238 at 6.674 eV, channel radius 9.4285 fm
        // ρ = k·a ≈ 2.188e-4 × √6.674 × 9.4285 ≈ 5.33e-3
        let awr = 236.006;
        let r = rho(6.674, awr, 9.4285);
        assert!((r - 5.33e-3).abs() < 1e-4, "ρ = {}, expected ~5.33e-3", r);
    }

    #[test]
    fn test_statistical_weight() {
        // U-238: I=0, s=1/2, J=1/2 → g = 2/2 = 1.0
        assert!((statistical_weight(0.5, 0.0) - 1.0).abs() < 1e-15);
        // J=3/2, I=0 → g = 4/2 = 2.0
        assert!((statistical_weight(1.5, 0.0) - 2.0).abs() < 1e-15);
        // Ta-181: I=7/2, J=3 → g = 7/(8×2) = 7/16
        assert!((statistical_weight(3.0, 3.5) - 7.0 / 16.0).abs() < 1e-15);
    }

    #[test]
    fn test_pi_over_k2_u238() {
        // At 6.674 eV for U-238:
        // k² = 2·m_n·(AWR/(1+AWR))²·E / (ℏc)²
        // With corrected formula: π/k² ≈ 9.86e4 barns
        // (factor of (1+AWR)/AWR = 1.0042 larger than uncorrected)
        let val = pi_over_k_squared_barns(6.674, 236.006);
        assert!(
            (val - 9.86e4).abs() < 1e3,
            "π/k² = {} barns, expected ~9.86e4",
            val
        );
    }

    /// k_squared must return 0.0 for negative energies (guard path).
    /// Negative E_r can appear for bound-state resonances in ENDF; the guard
    /// prevents a negative k² from propagating NaN through sqrt downstream.
    #[test]
    fn test_k_squared_negative_energy() {
        assert_eq!(k_squared(-1.0, 236.006), 0.0);
        assert_eq!(k_squared(-100.0, 10.0), 0.0);
        assert_eq!(k_squared(0.0, 236.006), 0.0);
    }

    /// wave_number_from_cm must return 0.0 for negative CM energies.
    /// Negative E_cm arises near threshold energies from numerical noise.
    #[test]
    fn test_wave_number_from_cm_negative_energy() {
        assert_eq!(wave_number_from_cm(-1.0, 0.5), 0.0);
        assert_eq!(wave_number_from_cm(-1e-10, 0.9412), 0.0);
        // Zero CM energy → zero wave number
        assert_eq!(wave_number_from_cm(0.0, 0.5), 0.0);
        // Positive CM energy → positive wave number
        assert!(wave_number_from_cm(1.0, 0.5) > 0.0);
    }
}