nereids_physics/

penetrability.rs

//! Hard-sphere penetrability, shift, and phase shift functions.
//!
//! These are the standard neutron penetrability factors for a hard-sphere
//! potential, needed for the R-matrix cross-section calculation.
//!
//! ## SAMMY Reference
//! - `rml/mrml07.f`: `Sinsix` (phase shifts), `Pgh` (penetrability + shift),
//!   `Pf` (penetrability only)
//! - SAMMY manual Section II.A, Table II A.1 (penetrability, level-shift, and
//!   phase-shift formulas)
//!
//! ## Definitions
//! For channel radius a and neutron wave number k:
//!   ρ = k·a
//!
//! The penetrability P_l(ρ), shift factor S_l(ρ), and hard-sphere phase
//! shift φ_l(ρ) depend on orbital angular momentum l.

use nereids_core::constants::NEAR_ZERO_FLOOR;

/// Penetrability factor P_l(ρ) for orbital angular momentum l.
///
/// Returns the ratio of the outgoing wave amplitude at the nuclear surface
/// to the asymptotic amplitude. For s-wave (l=0), P = ρ.
///
/// Reference: SAMMY `rml/mrml07.f` Pf function (lines 474-522)
pub fn penetrability(l: u32, rho: f64) -> f64 {
    let rho2 = rho * rho;
    match l {
        0 => rho,
        1 => {
            let denom = 1.0 + rho2;
            rho2 * rho / denom
        }
        2 => {
            let rho4 = rho2 * rho2;
            let denom = 9.0 + 3.0 * rho2 + rho4;
            rho4 * rho / denom
        }
        3 => {
            let rho4 = rho2 * rho2;
            let rho6 = rho4 * rho2;
            let denom = 225.0 + 45.0 * rho2 + 6.0 * rho4 + rho6;
            rho6 * rho / denom
        }
        4 => {
            let rho4 = rho2 * rho2;
            let rho6 = rho4 * rho2;
            let rho8 = rho4 * rho4;
            let denom = 11025.0 + 1575.0 * rho2 + 135.0 * rho4 + 10.0 * rho6 + rho8;
            rho8 * rho / denom
        }
        _ => penetrability_general(l, rho),
    }
}

/// Shift factor S_l(ρ).
///
/// The shift factor modifies the effective resonance energy.
/// For s-wave, S_0 = 0. For higher l, S_l is negative and energy-dependent.
///
/// Reference: SAMMY `rml/mrml07.f` Pgh function (lines 347-469)
pub fn shift_factor(l: u32, rho: f64) -> f64 {
    let rho2 = rho * rho;
    match l {
        0 => 0.0,
        1 => {
            let denom = 1.0 + rho2;
            -1.0 / denom
        }
        2 => {
            let rho4 = rho2 * rho2;
            let denom = 9.0 + 3.0 * rho2 + rho4;
            -(18.0 + 3.0 * rho2) / denom
        }
        3 => {
            let rho4 = rho2 * rho2;
            let rho6 = rho4 * rho2;
            let denom = 225.0 + 45.0 * rho2 + 6.0 * rho4 + rho6;
            -(675.0 + 90.0 * rho2 + 6.0 * rho4) / denom
        }
        4 => {
            let rho4 = rho2 * rho2;
            let rho6 = rho4 * rho2;
            let rho8 = rho4 * rho4;
            let denom = 11025.0 + 1575.0 * rho2 + 135.0 * rho4 + 10.0 * rho6 + rho8;
            -(44100.0 + 4725.0 * rho2 + 270.0 * rho4 + 10.0 * rho6) / denom
        }
        _ => shift_factor_general(l, rho),
    }
}

/// Hard-sphere phase shift φ_l(ρ).
///
/// Reference: SAMMY `rml/mrml07.f` Sinsix function (lines 254-342)
pub fn phase_shift(l: u32, rho: f64) -> f64 {
    let rho2 = rho * rho;
    match l {
        0 => rho,
        // For l = 1..=4 the phase is conceptually φ = ρ − atan(X/Y), but
        // forming `ρ − (X/Y).atan()` directly introduces a spurious −π jump
        // whenever the rational denominator Y crosses zero (the atan pole at
        // X/Y → ±∞).  SAMMY's `Sinsix` (rml/mrml07.f:254-342) never forms that
        // difference: it computes cos φ = (C·Y + S·X)/√G and sin φ =
        // (S·Y − C·X)/√G with G = X² + Y², C = cos ρ, S = sin ρ, which are
        // continuous through Y = 0.  Recovering φ = atan2(sin φ, cos φ) gives a
        // pole-free phase that is identical (mod π) to ρ − atan(X/Y) away from
        // the pole — the same continuous form the l > 4 path already uses.
        // l=1: φ = ρ − atan(ρ), i.e. X = ρ, Y = 1 (SAMMY uses G = 1 + X²,
        // Sinphi = (S − C·ρ)/√G, Cosphi = (C + S·ρ)/√G).  atan(ρ) has no pole,
        // so l=1 was already continuous; routing it through the same helper
        // keeps a single code path.
        1 => phase_shift_rational(rho, rho, 1.0),
        2 => phase_shift_rational(rho, 3.0 * rho, 3.0 - rho2),
        3 => phase_shift_rational(rho, rho * (15.0 - rho2), 15.0 - 6.0 * rho2),
        4 => {
            let rho4 = rho2 * rho2;
            phase_shift_rational(rho, rho * (105.0 - 10.0 * rho2), 105.0 - 45.0 * rho2 + rho4)
        }
        _ => phase_shift_general(l, rho),
    }
}

/// Continuous hard-sphere phase shift from the SAMMY `Sinsix` (X, Y) rational
/// form, avoiding the −π jump of `ρ − atan(X/Y)` at the Y = 0 pole.
///
/// Given the orbital rational numerator `x` and denominator `y` (the SAMMY `X`
/// and `Y` for this l), returns
///   φ = atan2(S·Y − C·X,  C·Y + S·X),   S = sin ρ,  C = cos ρ,
/// which equals ρ − atan(X/Y) for Y > 0 and differs by π for Y < 0 (so they
/// agree modulo π), and is continuous across Y = 0.
///
/// Because every consumer uses the phase only through π-invariant combinations
/// (sin²φ, sin 2φ, |U|²), removing the −π jump leaves all cross-sections
/// unchanged; it removes a discontinuity that would corrupt any future
/// phase-sensitive observable (e.g. angular distributions).
///
/// Reference: SAMMY `rml/mrml07.f` `Sinsix` (lines 254-342) —
/// `Cosphi = (C*Y+S*X)/√G`, `Sinphi = (S*Y-C*X)/√G`.
#[inline]
fn phase_shift_rational(rho: f64, x: f64, y: f64) -> f64 {
    let (s, c) = rho.sin_cos();
    (s * y - c * x).atan2(c * y + s * x)
}

/// Shift factor at imaginary channel argument S_l(iκ) for a closed channel.
///
/// For a channel with energy E_c < 0 (below threshold), the wave number is
/// imaginary: k_c = iκ with κ = sqrt(2μ|E_c|) / ħ. The penetrability is
/// zero by definition, but the shift factor S_l(iκ) is real and finite and
/// must be included in the level matrix diagonal `L_c = (S_c − B_c) + i·P_c`.
///
/// S_l(iκ) is obtained by analytic continuation of S_l(ρ) to imaginary
/// argument, i.e. substituting ρ² → −κ² in the Blatt-Weisskopf formula.
/// For l ≤ 4 the result is a closed-form rational function of κ².
/// For l > 4 a complex Bessel recursion is used.
///
/// Note: at κ = 1 for l = 1 the denominator vanishes (virtual-state pole).
/// This is a genuine physical singularity; the caller receives ±∞ from
/// floating-point arithmetic and should handle it via the normal `|L_c|`
/// guard in the Y-matrix construction.
///
/// Reference: Lane & Thomas, Rev. Mod. Phys. 30 (1958);
/// SAMMY `rml/mrml07.f` Pgh — PH = 1/(S−B+iP).
pub fn shift_factor_closed(l: u32, kappa: f64) -> f64 {
    // Substitute ρ² → −κ² in the explicit S_l(ρ) formulas.
    let k2 = kappa * kappa;
    match l {
        0 => 0.0,
        1 => -1.0 / (1.0 - k2),
        2 => {
            let k4 = k2 * k2;
            let denom = 9.0 - 3.0 * k2 + k4;
            -(18.0 - 3.0 * k2) / denom
        }
        3 => {
            let k4 = k2 * k2;
            let k6 = k4 * k2;
            let denom = 225.0 - 45.0 * k2 + 6.0 * k4 - k6;
            -(675.0 - 90.0 * k2 + 6.0 * k4) / denom
        }
        4 => {
            let k4 = k2 * k2;
            let k6 = k4 * k2;
            let k8 = k4 * k4;
            let denom = 11025.0 - 1575.0 * k2 + 135.0 * k4 - 10.0 * k6 + k8;
            -(44100.0 - 4725.0 * k2 + 270.0 * k4 - 10.0 * k6) / denom
        }
        _ => shift_factor_closed_general(l, kappa),
    }
}

/// General imaginary-argument shift factor for l > 4 via complex Bessel recursion.
///
/// Evaluates S_l(iκ) = Re[ρ(f∂f/∂ρ + g∂g/∂ρ)/(f²+g²)] at ρ = iκ
/// using the complex Bessel functions f_l(iκ) and g_l(iκ) computed via
/// upward recursion from the l=0 seed values:
///   f_0(iκ) = i·sinh(κ),  g_0(iκ) = cosh(κ)
fn shift_factor_closed_general(l: u32, kappa: f64) -> f64 {
    if kappa < NEAR_ZERO_FLOOR {
        return 0.0;
    }

    let (f, g) = bessel_fg_imaginary(l, kappa);
    let h = kappa * 1e-6 + 1e-12;
    let (fp, gp) = bessel_fg_imaginary(l, kappa + h);
    let (fm, gm) = bessel_fg_imaginary(l, kappa - h);

    // Numerical ∂/∂κ
    let df_re = (fp.0 - fm.0) / (2.0 * h);
    let df_im = (fp.1 - fm.1) / (2.0 * h);
    let dg_re = (gp.0 - gm.0) / (2.0 * h);
    let dg_im = (gp.1 - gm.1) / (2.0 * h);

    // ∂/∂ρ = (1/i)·∂/∂κ = −i·∂/∂κ
    // f·∂f/∂ρ = (f_re + i·f_im)·(df_im − i·df_re)
    //   Im part = f_im·df_im − f_re·df_re
    // ρ·(f·∂f/∂ρ + …): ρ = iκ, so Re[iκ·z] = −κ·Im[z]
    let num_im = (f.1 * df_im - f.0 * df_re) + (g.1 * dg_im - g.0 * dg_re);
    let numerator = -kappa * num_im;

    // f² + g² (complex square, not modulus squared); result is real
    let denom = f.0 * f.0 - f.1 * f.1 + g.0 * g.0 - g.1 * g.1;
    if denom.abs() < NEAR_ZERO_FLOOR {
        return 0.0;
    }
    numerator / denom
}

/// Bessel functions f_l(iκ) = iκ·j_l(iκ) and g_l(iκ) = −iκ·n_l(iκ)
/// evaluated at imaginary argument ρ = iκ (κ > 0), returned as
/// (f_re, f_im) and (g_re, g_im).
///
/// Seed values (l = 0):  f = (0, sinh κ),  g = (cosh κ, 0).
/// The upward factor (2n+1)/ρ = (2n+1)/(iκ) is purely imaginary.
fn bessel_fg_imaginary(l: u32, kappa: f64) -> ((f64, f64), (f64, f64)) {
    let sh = kappa.sinh();
    let ch = kappa.cosh();
    let inv_k = 1.0 / kappa;

    let mut f_prev = (0.0_f64, sh); // f_0 = i·sinh κ
    let mut g_prev = (ch, 0.0_f64); // g_0 = cosh κ

    if l == 0 {
        return (f_prev, g_prev);
    }

    // l = 1: f_1 = f_0/(iκ) − g_0 = sinh(κ)/κ − cosh(κ)  [real]
    //         g_1 = g_0/(iκ) + f_0 = i·(sinh(κ) − cosh(κ)/κ) [imaginary]
    let mut f_curr = (sh * inv_k - ch, 0.0);
    let mut g_curr = (0.0, sh - ch * inv_k);

    if l == 1 {
        return (f_curr, g_curr);
    }

    // Upward recursion: factor = (2n+1)/(iκ), purely imaginary with imaginary part a.
    // (a·i)·(r + i·s) = −a·s + i·a·r
    for n in 1..(l as usize) {
        let a = -((2 * n + 1) as f64) * inv_k; // imaginary part of (2n+1)/(iκ)
        let f_next = (-a * f_curr.1 - f_prev.0, a * f_curr.0 - f_prev.1);
        let g_next = (-a * g_curr.1 - g_prev.0, a * g_curr.0 - g_prev.1);
        f_prev = f_curr;
        g_prev = g_curr;
        f_curr = f_next;
        g_curr = g_next;
    }

    (f_curr, g_curr)
}

/// General penetrability via recursion for l > 4.
///
/// Uses the recursion relation:
///   P_l = ρ / |F_l + i·G_l|²
/// where F_l, G_l are regular/irregular Coulomb functions (neutral case).
///
/// For the hard-sphere case (no Coulomb), we use the recursion:
///   P_l = ρ·P_{l-1} / [(2l-1) - ρ²·P_{l-1}/P_l ... ]
/// implemented via backwards recursion of f_l, g_l.
fn penetrability_general(l: u32, rho: f64) -> f64 {
    // Compute via backwards recursion of spherical Bessel functions.
    // P_l = ρ / (f_l² + g_l²) where f_l, g_l are defined by:
    //   f_l(ρ) = ρ·j_l(ρ), g_l(ρ) = -ρ·n_l(ρ)
    // Use the standard recursion relations.
    let (fl, gl) = bessel_fg(l, rho);
    rho / (fl * fl + gl * gl)
}

/// General shift factor via recursion for l > 4.
fn shift_factor_general(l: u32, rho: f64) -> f64 {
    let (fl, gl) = bessel_fg(l, rho);
    let denom = fl * fl + gl * gl;
    // S_l = ρ·(f_l·f_l' + g_l·g_l') / (f_l² + g_l²)
    // Using the Wronskian relation.
    // For the hard-sphere: S_l = (l+1) - ρ²/(S_{l-1} + (l) - iP_{l-1}) ... complicated
    // Simpler: use recursion S_l = ρ² / (l as f64 - S_{l-1}) - l as f64
    // But this is unstable. Use numerical derivative instead.
    let h = rho * 1e-6 + 1e-12;
    let (fl_p, gl_p) = bessel_fg(l, rho + h);
    let (fl_m, gl_m) = bessel_fg(l, rho - h);
    let p_p = (rho + h) / (fl_p * fl_p + gl_p * gl_p);
    let p_m = (rho - h) / (fl_m * fl_m + gl_m * gl_m);
    let _dp_drho = (p_p - p_m) / (2.0 * h);
    // S_l = ρ·dP_l/dρ / P_l - P_l + l(l+1)/ρ ... no, use relation:
    // Easier: S_l(ρ) = ρ·(f_l·df_l/dρ + g_l·dg_l/dρ) / (f_l² + g_l²)
    let dfl = (fl_p - fl_m) / (2.0 * h);
    let dgl = (gl_p - gl_m) / (2.0 * h);
    rho * (fl * dfl + gl * dgl) / denom
}

/// General phase shift via recursion for l > 4.
fn phase_shift_general(l: u32, rho: f64) -> f64 {
    let (fl, gl) = bessel_fg(l, rho);
    fl.atan2(gl)
}

/// Compute f_l(ρ) = ρ·j_l(ρ) and g_l(ρ) = -ρ·n_l(ρ) via upward recursion.
///
/// j_l and n_l are spherical Bessel functions of the first and second kind.
fn bessel_fg(l: u32, rho: f64) -> (f64, f64) {
    // Near-zero guard.  The sign of g_0 is mathematically +1 (cos(0)),
    // but we return -1.0 because this path is only reachable for l == 0
    // (higher l returns NEG_INFINITY) and only P_l = rho / (f^2 + g^2)
    // uses the result, where g^2 = 1 regardless of sign.
    if rho.abs() < NEAR_ZERO_FLOOR {
        return if l == 0 {
            (0.0, -1.0)
        } else {
            (0.0, f64::NEG_INFINITY)
        };
    }

    // Start with l=0:
    // f_0 = ρ·j_0(ρ) = sin(ρ)
    // g_0 = -ρ·n_0(ρ) = cos(ρ)
    let mut f_prev = rho.sin();
    let mut g_prev = rho.cos();

    if l == 0 {
        return (f_prev, g_prev);
    }

    // l=1:
    // f_1 = ρ·j_1(ρ) = sin(ρ)/ρ - cos(ρ)
    // g_1 = -ρ·n_1(ρ) = cos(ρ)/ρ + sin(ρ)
    let mut f_curr = f_prev / rho - g_prev;
    let mut g_curr = g_prev / rho + f_prev;

    if l == 1 {
        return (f_curr, g_curr);
    }

    // Upward recursion for l >= 2:
    // f_{l+1} = ((2l+1)/ρ)·f_l - f_{l-1}
    // g_{l+1} = ((2l+1)/ρ)·g_l - g_{l-1}
    for n in 1..(l as i64) {
        let factor = (2 * n + 1) as f64 / rho;
        let f_next = factor * f_curr - f_prev;
        let g_next = factor * g_curr - g_prev;
        f_prev = f_curr;
        g_prev = g_curr;
        f_curr = f_next;
        g_curr = g_next;
    }

    (f_curr, g_curr)
}

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

    #[test]
    fn test_penetrability_l0() {
        // P_0(ρ) = ρ
        assert!((penetrability(0, 0.5) - 0.5).abs() < 1e-15);
        assert!((penetrability(0, 1.0) - 1.0).abs() < 1e-15);
        assert!((penetrability(0, 0.001) - 0.001).abs() < 1e-15);
    }

    #[test]
    fn test_penetrability_l1() {
        // P_1(ρ) = ρ³/(1 + ρ²)
        let rho = 0.5;
        let expected = 0.5_f64.powi(3) / (1.0 + 0.25);
        assert!((penetrability(1, rho) - expected).abs() < 1e-15);

        let rho = 1.0;
        let expected = 1.0 / 2.0;
        assert!((penetrability(1, rho) - expected).abs() < 1e-15);
    }

    #[test]
    fn test_penetrability_l2() {
        // P_2(ρ) = ρ⁵/(9 + 3ρ² + ρ⁴)
        let rho = 1.0;
        let expected = 1.0 / (9.0 + 3.0 + 1.0);
        assert!((penetrability(2, rho) - expected).abs() < 1e-15);
    }

    #[test]
    fn test_shift_factor_l0() {
        // S_0 = 0 always
        assert!((shift_factor(0, 0.5)).abs() < 1e-15);
        assert!((shift_factor(0, 5.0)).abs() < 1e-15);
    }

    #[test]
    fn test_shift_factor_l1() {
        // S_1(ρ) = -1/(1 + ρ²)
        let rho = 1.0;
        assert!((shift_factor(1, rho) - (-0.5)).abs() < 1e-15);
        let rho = 0.0;
        assert!((shift_factor(1, rho) - (-1.0)).abs() < 1e-15);
    }

    #[test]
    fn test_phase_shift_l0() {
        // φ_0(ρ) = ρ
        assert!((phase_shift(0, 0.5) - 0.5).abs() < 1e-15);
        assert!((phase_shift(0, 1.0) - 1.0).abs() < 1e-15);
    }

    #[test]
    fn test_phase_shift_l1() {
        // φ_1(ρ) = ρ - arctan(ρ)
        let rho = 1.0;
        let expected = 1.0 - 1.0_f64.atan();
        assert!((phase_shift(1, rho) - expected).abs() < 1e-14);
    }

    #[test]
    fn test_phase_shift_small_rho() {
        // For small ρ, φ_l ≈ ρ^(2l+1) / (1·3·5···(2l+1))²
        // φ_0(0.01) ≈ 0.01
        assert!((phase_shift(0, 0.01) - 0.01).abs() < 1e-10);
        // φ_1(0.01) ≈ 0.01 - arctan(0.01) ≈ 3.33e-7
        let expected = 0.01 - 0.01_f64.atan();
        assert!((phase_shift(1, 0.01) - expected).abs() < 1e-15);
    }

    #[test]
    fn test_general_matches_explicit_l2() {
        // Verify that the general recursion gives the same result as explicit l=2
        let rho = 2.0;
        let p_explicit = penetrability(2, rho);
        let p_general = penetrability_general(2, rho);
        assert!(
            (p_explicit - p_general).abs() < 1e-10,
            "explicit={}, general={}",
            p_explicit,
            p_general
        );
    }

    // ── shift_factor_closed tests ─────────────────────────────────────────────

    #[test]
    fn test_shift_factor_closed_l0() {
        // S_0 = 0 regardless of κ
        assert_eq!(shift_factor_closed(0, 0.0), 0.0);
        assert_eq!(shift_factor_closed(0, 2.0), 0.0);
    }

    #[test]
    fn test_shift_factor_closed_l1_continuity() {
        // At κ → 0: S_1(iκ) = −1/(1−κ²) → −1 = S_1(0) from the real-ρ formula.
        // Continuity: S_1(ρ=0) = −1/(1+0) = −1; S_1(i·0) = −1/(1−0) = −1.
        let s = shift_factor_closed(1, 1e-10);
        assert!((s - (-1.0)).abs() < 1e-8, "got {s}");
    }

    #[test]
    fn test_shift_factor_closed_l1_exact() {
        // S_1(i·0.5) = −1/(1−0.25) = −4/3
        let expected = -4.0 / 3.0;
        let got = shift_factor_closed(1, 0.5);
        assert!(
            (got - expected).abs() < 1e-12,
            "got {got}, expected {expected}"
        );
    }

    #[test]
    fn test_shift_factor_closed_l2_exact() {
        // S_2(i·0.5): k2=0.25, k4=0.0625
        // num = -(18−3·0.25) = -(18−0.75) = −17.25
        // den = 9−3·0.25+0.0625 = 9−0.75+0.0625 = 8.3125
        let expected = -17.25 / 8.3125;
        let got = shift_factor_closed(2, 0.5);
        assert!(
            (got - expected).abs() < 1e-12,
            "got {got}, expected {expected}"
        );
    }

    #[test]
    fn test_shift_factor_closed_general_matches_explicit_l4() {
        // Verify general recursion against closed-form for l=4.
        let kappa = 0.3;
        let explicit = shift_factor_closed(4, kappa);
        let general = shift_factor_closed_general(4, kappa);
        assert!(
            (explicit - general).abs() < 1e-8,
            "l=4 kappa={kappa}: explicit={explicit}, general={general}"
        );
    }

    #[test]
    fn test_general_matches_explicit_l3() {
        let rho = 1.5;
        let p_explicit = penetrability(3, rho);
        let p_general = penetrability_general(3, rho);
        assert!(
            (p_explicit - p_general).abs() < 1e-10,
            "explicit={}, general={}",
            p_explicit,
            p_general
        );
    }

    // ── Fix #4: phase_shift continuity across the rational-form pole ──────────

    /// The old one-argument `ρ − atan(X/Y)` form jumps by −π whenever the
    /// rational denominator Y crosses zero.  The Sinsix `atan2` form must be
    /// continuous there.  The poles are: l=2 at ρ=√3, l=3 at ρ=√2.5,
    /// l=4 at ρ²=(45−√1605)/2 (≈1.572) and ρ²=(45+√1605)/2 (≈6.52).
    #[test]
    fn test_phase_shift_continuous_across_pole() {
        // (l, ρ_pole) pairs.
        let poles = [
            (2u32, 3.0_f64.sqrt()),
            (3u32, 2.5_f64.sqrt()),
            (4u32, ((45.0 - 1605.0_f64.sqrt()) / 2.0).sqrt()),
            (4u32, ((45.0 + 1605.0_f64.sqrt()) / 2.0).sqrt()),
        ];
        for (l, rho_pole) in poles {
            let delta = 1e-6;
            let below = phase_shift(l, rho_pole - delta);
            let above = phase_shift(l, rho_pole + delta);
            let jump = (above - below).abs();
            // A continuous φ varies by O(delta) across the pole; the buggy
            // form jumps by ≈ π.  Require the jump to be far below π.
            assert!(
                jump < 1e-3,
                "phase_shift(l={l}) discontinuous at ρ={rho_pole}: \
                 below={below}, above={above}, jump={jump} (expected ≈ 0, not ≈ π)"
            );
        }
    }

    /// The continuous Sinsix form must agree with the original
    /// `ρ − atan(X/Y)` form modulo π (they coincide for Y > 0 and differ by π
    /// for Y < 0) everywhere away from the pole, so every
    /// cross-section observable (which consumes φ only via sin²φ / sin2φ /
    /// |U|²) is bit-for-bit unchanged.  We assert the π-invariant combinations
    /// match the old closed form to round-off.
    #[test]
    fn test_phase_shift_matches_old_form_on_invariants() {
        // Old (discontinuous) reference: ρ − atan(X/Y).
        fn old_phase(l: u32, rho: f64) -> f64 {
            let rho2 = rho * rho;
            match l {
                1 => rho - rho.atan(),
                2 => rho - (3.0 * rho / (3.0 - rho2)).atan(),
                3 => {
                    let num = rho * (15.0 - rho2);
                    let den = 15.0 - 6.0 * rho2;
                    rho - (num / den).atan()
                }
                4 => {
                    let rho4 = rho2 * rho2;
                    let num = rho * (105.0 - 10.0 * rho2);
                    let den = 105.0 - 45.0 * rho2 + rho4;
                    rho - (num / den).atan()
                }
                _ => unreachable!(),
            }
        }

        for l in 1u32..=4 {
            // Sample ρ across a range that brackets every pole, skipping a
            // tiny neighbourhood of each pole (where the OLD form is ill-defined
            // but the NEW form is fine).
            let poles = [3.0_f64.sqrt(), 2.5_f64.sqrt(), 1.572, 6.52];
            for i in 1..=800 {
                let rho = i as f64 * 0.01; // 0.01 .. 8.0
                if poles.iter().any(|&p| (rho - p).abs() < 5e-3) {
                    continue;
                }
                let new = phase_shift(l, rho);
                let old = old_phase(l, rho);
                // π-invariant observables consumed downstream:
                let d_sin2 = (new.sin().powi(2) - old.sin().powi(2)).abs();
                let d_sin2ph = ((2.0 * new).sin() - (2.0 * old).sin()).abs();
                let d_cos2ph = ((2.0 * new).cos() - (2.0 * old).cos()).abs();
                assert!(
                    d_sin2 < 1e-9 && d_sin2ph < 1e-9 && d_cos2ph < 1e-9,
                    "l={l} ρ={rho}: new φ disagrees with old on invariants \
                     (sin²Δ={d_sin2}, sin2φΔ={d_sin2ph}, cos2φΔ={d_cos2ph})"
                );
            }
        }
    }
}