diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index cb55bdf..a1e1dd6 100644 --- a/beta_spectrum/components/radiative.py +++ b/beta_spectrum/components/radiative.py @@ -1,40 +1,66 @@ +# components/radiative.py import numpy as np -from scipy.special import log_wright_bessel, spence +from scipy.special import spence from beta_spectrum.base import SpectrumComponent -from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA +from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA, HBAR_C_MEV_FM class RadiativeCorrection(SpectrumComponent): """ - Outer radiative correction R(W, W0) = 1 + delta(W, W0). + Outer radiative correction R(W, W0) = 1 + delta_1 + delta_2. - Implements Eq. (47)-(50) from Hayen et al. 2017, with endpoint - divergence handeled by the soft photon resummation of Eq. (53). + Implements: + - delta_1: O(alpha) universal Sirlin function g(W, W0) from + Sirlin (1967) Eq. (20b), with endpoint soft-photon resummation + per Sirlin (1987) Eq. (53). + - delta_2: O(Z*alpha^2) finite nuclear size correction from + Sirlin (1987) Eq. (1-10), using the modified Gaussian model. This is the INCLUSIVE correction suitable for: - - Magnetic spectrometers - - Detectors where bremmstrahlung is not detected - - Precision at 1% level or worse + - Magnetic spectrometers + - Detectors where bremsstrahlung is not detected + - Precision at 1% level or worse + + The O(Z^2*alpha^3) correction is neglected for Z < 30 (spec Section 5). + Inner radiative corrections are absorbed into effective coupling constants + and are not part of the spectrum shape calculation. """ - def __init__(self, W0: float, use_endpoint_resummation: bool = True): + def __init__( + self, + W0: float, + Z: int = 0, + A: int | None = None, + use_endpoint_resummation: bool = True, + delta_cut: float = 1e-3, + ): """ Args: W0: Endpoint energy in natural units m_e=c=h_bar=1 - use_endpoint_resummation: Apply eq. (53) to handle ln(W0 - W) divergence + Z: Nuclear charge of daughter nucleus + A: Mass number (optional, affects nuclear model correction) + use_endpoint_resummation: Apply Eq.(53) to handle ln(W0 - W) divergence + delta_cut: Only apply resummation when (W0 - W) < delta_cut (spec Section 3.4) """ self.W0 = W0 + self.Z = Z + self.A = A self.use_endpoint_resummation = use_endpoint_resummation - self.m_p = MP_MEV / ME_MEV # ~1836.15 + self.delta_cut = delta_cut + self.m_p = MP_MEV / ME_MEV # ~1836.15 (proton mass in m_e units) + self.gamma_E = 0.57721566490153286 # Euler-Mascheroni constant - # Precompute constant terms of Eq. (48) + # Precompute constant terms self._outer_normalization = ALPHA / (2.0 * np.pi) # ~0.00116 self._subtraction_term = 3.0 * np.log(self.m_p / (2.0 * self.W0)) + # Precompute O(Z*alpha^2) energy-independent nuclear part if A is given + self._delta_F = self._compute_nuclear_model_correction(Z, A) + def __call__(self, W: np.ndarray) -> np.ndarray: """ - Calculate R(W, W0) = 1 + delta_R(w, W0) + Calculate R(W, W0) = 1 + delta_1(W, W0) + delta_2(Z, W). Args: W: Electron total energy array (in natural m_e=c=h_bar=1 units) @@ -42,9 +68,9 @@ class RadiativeCorrection(SpectrumComponent): Returns: Radiative correction factor array """ - W = np.asarray(W) + W = np.asarray(W, dtype=float) - # Avoid numerical issues exactly at enpoint W0 + # Avoid numerical issues exactly at endpoint W0 eps = 1e-9 mask_finite = W < (self.W0 - eps) @@ -53,107 +79,196 @@ class RadiativeCorrection(SpectrumComponent): if np.any(mask_finite): W_finite = W[mask_finite] - if self.use_endpoint_resummation: - g_func = self._g_function_resummed(W_finite) - else: - g_func = self._g_function_standard(W_finite) + # O(alpha) correction (Sirlin 1967) + delta_1 = self._delta_1(W_finite) - delta_R_finite = self._outer_normalization * ( - g_func - self._subtraction_term - ) + # O(Z*alpha^2) correction (Sirlin 1987) + delta_2 = self._delta_2(W_finite) + + delta_R_finite = delta_1 + delta_2 delta_R[mask_finite] = delta_R_finite - # Apply full correction including inner and L factor return 1.0 + delta_R - def _g_function_standard(self, W: np.ndarray) -> np.ndarray: + def _delta_1(self, W: np.ndarray) -> np.ndarray: """ - Standard delta_R from Eq. (50), without endpoint resummation. - Has logarithmic divergence as W -> W0 - """ - beta = np.sqrt(1.0 - 1.0 / W**2) - beta = np.maximum(beta, 1e-12) # Avoid division by 0 + O(alpha) universal Sirlin function g(W, W0) per Sirlin (1967) Eq. (20b). - # Dilogarithm term + g(W, W0) = 3*ln(m_p/(2*W0)) - 3/4 + + 4*(tanh_inv_beta/beta - 1)*[(W0-W)/(3W) - 3/2 + ln(2*(W0-W))] + + (4/beta)*L(2*beta/(1+beta)) + + (tanh_inv_beta/beta)*[2*(1+beta^2) + (W0-W)^2/(6*W^2) - 4*tanh_inv_beta] + + where L(x) = -Li_2(x) is the Spence function (dilogarithm). + + The correction is delta_1 = (alpha/2pi) * (g(W, W0) - subtraction_term). + The subtraction_term cancels the constant part in the Hayen reformulation. + """ + beta = self._safe_beta(W) + tanh_inv_beta_beta = self._safe_tanh_inv_beta_over_beta(beta) + + # Dilogarithm term: L(2*beta/(1+beta)) = -Li_2(2*beta/(1+beta)) arg = 2.0 * beta / (1.0 + beta) arg = np.minimum(arg, 1.0 - 1e-12) # Avoid numerical problems at arg=1 dilog_term = (4.0 / beta) * spence(1.0 - arg) - # arctanh factor - atanh_beta = np.arctanh(beta) - atanh_factor = (atanh_beta / beta) - 1.0 + # arctanh factor: (tanh^{-1}(beta)/beta) - 1 + atanh_factor = tanh_inv_beta_beta - 1.0 # Energy difference delta_W = self.W0 - W delta_W = np.where(delta_W > 0, delta_W, 1e-12) - # Second line Eq. (50) - term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + np.log(2.0 * delta_W)) + if self.use_endpoint_resummation: + # Apply soft-photon resummation from Sirlin (1987) + # Replace C(beta)*ln(W0-W) with C(beta)*((W0-W)^{t(beta)}-1)/t(beta) + # Only near endpoint: (W0 - W) < delta_cut (spec Section 3.4) + t_beta = (2.0 * ALPHA / np.pi) * atanh_factor + small_t = np.abs(t_beta) < 1e-6 - # Third line Eq. (50) - term3 = (atanh_beta / beta) * ( - 2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta + resummed_log = np.empty_like(delta_W) + resummed_log[:] = np.log(delta_W) + + near_endpoint = delta_W < self.delta_cut + if np.any(near_endpoint): + resummed_log[near_endpoint] = ( + np.power(delta_W[near_endpoint], t_beta[near_endpoint]) - 1.0 + ) / t_beta[near_endpoint] + small_t_mask = near_endpoint & small_t + resummed_log[small_t_mask] = np.log(delta_W[small_t_mask]) + + log_term = np.log(2.0) + resummed_log + else: + log_term = np.log(2.0 * delta_W) + + # Energy-dependent term + term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term) + + # Beta-dependent term (no divergence) + term3 = tanh_inv_beta_beta * ( + 2.0 * (1.0 + beta**2) + + (delta_W**2) / (6.0 * W**2) + - 4.0 * tanh_inv_beta_beta ) - # Leading constants from Eq. (50) + # Leading constant const_term = -0.75 - return const_term + dilog_term + term2 + term3 + g_func = const_term + dilog_term + term2 + term3 - def _g_function_resummed(self, W: np.ndarray) -> np.ndarray: - """ - delta_R with soft photon resummation to handle endpoint divergence, - Implements the replacement from Eq. (53): + return self._outer_normalization * (g_func - self._subtraction_term) - t(beta) * ln(W0 - W) -> (W0 - W)^t(beta) - 1 + def _delta_2(self, W: np.ndarray) -> np.ndarray: """ - beta = np.sqrt(1.0 - 1.0 / W**2) + O(Z*alpha^2) finite nuclear size correction per Sirlin (1987). + + delta_2(Z, W) = Z * alpha^2 * (Delta_1_4(W) + Delta_F) + + where: + - Delta_1_4(W) = ln(m_p) - (5/12)*ln(2*W) + 43/18 + (energy-dependent relativistic part) + - Delta_F(Z, A) = nuclear-structure-dependent part from modified + Gaussian model (Eq. 10 of Sirlin 1987) + + For Z=0, this correction is zero. + """ + if self.Z == 0: + return np.zeros_like(W) + + # Energy-dependent part: Delta_1^0(W) + Delta_4(W) + # From Sirlin 1987 relativistic limit + Delta_1_4 = np.log(self.m_p) - (5.0 / 12.0) * np.log(2.0 * W) + 43.0 / 18.0 + + # Total O(Z*alpha^2) correction + delta_2 = self.Z * ALPHA**2 * (Delta_1_4 + self._delta_F) + + return delta_2 + + def _compute_nuclear_model_correction(self, Z: int, A: int | None) -> float: + """ + Compute the nuclear-structure-dependent part Delta_F of the O(Z*alpha^2) correction. + + Uses the modified Gaussian model from Sirlin (1987) Eq. (10): + + Delta_F = ln(Lambda/M) - kappa_1(Z) + + where: + - Lambda = sqrt(6)/a, with a = r_0 * A^{1/3} (rms charge radius) + - M = proton mass + - kappa_1(Z) = 0.5 * [gamma_E + ln(3/(2*k^2)) + 2*alpha/(2+3*alpha)] + - k^2 = (3/2) * (2+5*alpha)/(2+3*alpha) + - alpha = (Z-2)/3 + - r_0 = 1.2 fm (Wilkinson) + + If A is not specified, returns 0 (no nuclear model correction). + """ + if A is None: + return 0.0 + + # Nuclear radius from Wilkinson's r_0 + r_0 = 1.2 # fm + a_fm = r_0 * A ** (1.0 / 3.0) # rms charge radius in fm + a_me = a_fm / (HBAR_C_MEV_FM / ME_MEV) # convert to m_e^{-1} units + + # Lambda = sqrt(6)/a + Lambda = np.sqrt(6.0) / a_me + + # M/Lambda ratio for the log term + M_over_Lambda = self.m_p / Lambda + + # kappa_1(Z) from Sirlin 1987 + alpha_param = (Z - 2.0) / 3.0 + k_sq = 1.5 * (2.0 + 5.0 * alpha_param) / (2.0 + 3.0 * alpha_param) + kappa_1 = 0.5 * ( + self.gamma_E + + np.log(3.0 / (2.0 * k_sq)) + + 2.0 * ALPHA / (2.0 + 3.0 * ALPHA) + ) + + delta_F = np.log(M_over_Lambda) - kappa_1 + return delta_F + + def _safe_beta(self, W: np.ndarray) -> np.ndarray: + """ + Compute beta = p/W = sqrt(1 - 1/W^2) with safe handling near threshold. + + At W=1 (threshold), beta=0. We use a floor to avoid division by zero. + """ + beta = np.sqrt(np.maximum(1.0 - 1.0 / W**2, 0.0)) beta = np.maximum(beta, 1e-12) + return beta - # t(beta) from Eq. (54) - atanh_beta = np.arctanh(beta) - t_beta = (2.0 * ALPHA / np.pi) * (atanh_beta / beta - 1.0) + def _safe_tanh_inv_beta_over_beta(self, beta: np.ndarray) -> np.ndarray: + """ + Compute tanh^{-1}(beta)/beta with numerical stability for small beta. - # Dilogarithm term (same as _delta_R_standard) - arg = 2.0 * beta / (1.0 + beta) - arg = np.minimum(arg, 1.0 - 1e-12) - dilog_term = (4.0 / beta) * spence(1.0 - arg) + For beta < 1e-6, use Taylor expansion: + tanh^{-1}(beta)/beta = 1 + beta^2/3 + beta^4/5 + ... - # arctanh factor - atanh_factor = (atanh_beta / beta) - 1.0 + For beta >= 1e-6, use the direct formula: + tanh^{-1}(beta)/beta = 0.5 * ln((1+beta)/(1-beta)) / beta - # Energy difference - delta_W = self.W0 - W - delta_W = np.where(delta_W > 0, delta_W, 1e-12) + Per spec Section 3.3, this avoids division by zero and loss of precision. + """ + result = np.empty_like(beta) + small_beta = beta < 1e-6 - # Apply resummation from Eq. (53) - # Replace divergent ln(W0 - W) term with (W0 - W)^t(beta) - 1 - - # Stable handling of small t_beta - small_t = np.abs(t_beta) < 1e-6 - - resummed_log = np.empty_like(delta_W) - - # For normal values - resummed_log[~small_t] = ( - np.power(delta_W[~small_t], t_beta[~small_t]) - 1.0 - ) / t_beta[~small_t] - - # For very small t -> recover ln(delta_W) - resummed_log[small_t] = np.log(delta_W[small_t]) - - # Reconstruct ln(2*delta_W) properly - log_term = np.log(2.0) + resummed_log - - # Term2 with resummed log (instead of ln(2*(W0 - W))) - term2_resummed = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term) - - # Term3 (no divergence) - term3 = (atanh_beta / beta) * ( - 2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta + # Taylor series for small beta: 1 + beta^2/3 + beta^4/5 + beta^6/7 + beta_sq = beta**2 + result[small_beta] = ( + 1.0 + + beta_sq[small_beta] / 3.0 + + beta_sq[small_beta] ** 2 / 5.0 + + beta_sq[small_beta] ** 3 / 7.0 ) - # Leading constants - const_term = -0.75 + # Direct computation for larger beta + not_small = ~small_beta + if np.any(not_small): + result[not_small] = ( + 0.5 + * np.log((1.0 + beta[not_small]) / (1.0 - beta[not_small])) + / beta[not_small] + ) - return const_term + dilog_term + term2_resummed + term3 + return result diff --git a/docs/06-radiative-corrections.md b/docs/06-radiative-corrections.md index b04465a..92f6aee 100644 --- a/docs/06-radiative-corrections.md +++ b/docs/06-radiative-corrections.md @@ -100,12 +100,13 @@ Experimental branching ratios: $BR_{\beta\gamma} \sim 3\times 10^{-3}$ (neutron) ## Calculator Implementation > [!tip] Code mapping -> Module: `beta_spectrum/components/radiative.py` → `RadiativeCorrection(W0, use_endpoint_resummation=True)` class +> Module: `beta_spectrum/components/radiative.py` → `RadiativeCorrection(W0, Z, A, use_endpoint_resummation)` class > -> Implements outer radiative corrections per Hayen Eqs. (47)–(53): -> - Sirlin's $g(W_0,W)$ with `scipy.special.spence` for dilogarithm -> - Soft-photon resummation at endpoint via Eq. (53) -> - Order Zα² terms Δ₁ through Δ₄ +> Implements outer radiative corrections per Sirlin (1967) and Sirlin (1987): +> - **O(α)**: Universal Sirlin function $g(W,W_0)$ with `scipy.special.spence` for dilogarithm +> - **O(Zα²)**: Finite nuclear size correction $\delta_2(Z,W)$ using modified Gaussian model +> - Soft-photon resummation at endpoint via Sirlin (1987) prescription +> - Small-beta Taylor expansion for numerical stability at threshold > - Optional endpoint resummation flag ### API Usage @@ -113,10 +114,21 @@ Experimental branching ratios: $BR_{\beta\gamma} \sim 3\times 10^{-3}$ (neutron) ```python from beta_spectrum import RadiativeCorrection +# Light nucleus (tritium-like) rad = RadiativeCorrection( - W0=1.02, # endpoint in m_e units - use_endpoint_resummation=True # handles endpoint divergence + W0=1.02, + Z=1, + use_endpoint_resummation=True ) + +# Heavy nucleus (uranium) with nuclear model correction +rad = RadiativeCorrection( + W0=5.0, + Z=92, + A=238, # affects Delta_F nuclear model term + use_endpoint_resummation=True +) + values = rad(W_grid) ``` diff --git a/tests/test_radiative.py b/tests/test_radiative.py index 0b47414..9f233c5 100644 --- a/tests/test_radiative.py +++ b/tests/test_radiative.py @@ -26,7 +26,6 @@ Common practice: For functions with known singularities, always test both "handl """ import numpy as np -import pytest from beta_spectrum.components.radiative import RadiativeCorrection @@ -41,9 +40,9 @@ class TestRadiativeBasicProperties: W_test = np.array([2.0]) result = rc(W_test) - assert 0.5 < result[0] < 3.0, ( - f"Radiative correction should be near unity: got {result[0]}" - ) + assert ( + 0.5 < result[0] < 3.0 + ), f"Radiative correction should be near unity: got {result[0]}" def test_positive_values(self): """R(W,W₀) must always be positive — it's a multiplicative factor.""" @@ -70,7 +69,9 @@ class TestRadiativeEndpointHandling: result = rc(W_near_end) assert not np.isnan(result[0]), "Resummed correction must be finite at endpoint" - assert not np.isinf(result[0]), "Resummed correction must not diverge at endpoint" + assert not np.isinf( + result[0] + ), "Resummed correction must not diverge at endpoint" def test_standard_mode_diverges(self): """Without resummation, δᵣ → −∞ as W → W₀ (logarithmic divergence). @@ -81,16 +82,16 @@ class TestRadiativeEndpointHandling: """ rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) - W_far = np.array([4.9]) # ΔW = 0.1 + W_far = np.array([4.9]) # ΔW = 0.1 W_close = np.array([4.99]) # ΔW = 0.01 val_far = rc(W_far)[0] val_close = rc(W_close)[0] # The correction should decrease (become more negative) near endpoint - assert val_close < val_far, ( - "Standard mode must show divergent behavior near endpoint" - ) + assert ( + val_close < val_far + ), "Standard mode must show divergent behavior near endpoint" class TestRadiativeEnergyDependence: @@ -112,9 +113,9 @@ class TestRadiativeEnergyDependence: r_near = rc(near_end_W)[0] # δᵣ becomes more negative (or less positive) toward endpoint → R decreases - assert r_near <= r_mid, ( - f"Correction factor should decrease near endpoint: {r_mid} → {r_near}" - ) + assert ( + r_near <= r_mid + ), f"Correction factor should decrease near endpoint: {r_mid} → {r_near}" class TestRadiativeOutputShape: @@ -141,27 +142,238 @@ class TestRadiativeResummationSwitch: """Test that the resummation flag actually changes behavior.""" def test_resummed_differs_from_standard(self): - """Near endpoint, resummed and standard modes should give different results. + """Near endpoint (within delta_cut), resummed and standard modes should differ. - At high energy (far from endpoint) they converge — both reduce to the same δᵣ ≈ α/π × [constant]. - Near W₀ the difference grows because resummation replaces ln(ΔW) with a finite power law. + Far from endpoint, both modes are identical. Near W₀ the difference grows + because resummation replaces ln(ΔW) with a finite power law. """ rc_resummed = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) - # Far from endpoint: both should be similar (both finite) + # Far from endpoint: both should be identical (resummation not applied) W_far = np.array([3.0]) - diff_far = abs(rc_resummed(W_far)[0] - rc_standard(W_far)[0]) + assert np.isclose(rc_resummed(W_far)[0], rc_standard(W_far)[0], rtol=1e-10) - # Near endpoint: check that values are not identical - # With resummation, R stays near 1; without it, δᵣ diverges negatively - W_near = np.array([4.99]) + # Near endpoint (within delta_cut=1e-3): resummed stays finite, standard diverges + W_near = np.array([4.9995]) r_resummed = rc_resummed(W_near)[0] + r_standard = rc_standard(W_near)[0] assert not np.isnan(r_resummed), "Resummed mode must be finite" + assert not np.isnan( + r_standard + ), "Standard mode must be finite (mask protects endpoint)" + assert ( + r_resummed > r_standard + ), f"Near endpoint: resummed={r_resummed} should be > standard={r_standard}" - # The key test: resummed stays close to 1 while standard diverges - # We verify they produce different values (exact threshold depends on implementation) - diff_near = abs(rc_resummed(W_near)[0] - rc_standard(W_near)[0]) + def test_resummation_only_near_endpoint(self): + """Per spec Section 3.4, resummation should only apply when (W0 - W) < delta_cut. - assert diff_far > 0, "Even far from endpoint there should be some difference" + Far from endpoint (delta_W > delta_cut), resummed and standard modes must + produce identical results. Near endpoint (delta_W < delta_cut), resummed + mode stays finite while standard mode becomes more negative. + """ + rc_resummed = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) + + # delta_cut defaults to 1e-3, so W0 - W = 0.01 > delta_cut + W_far = np.array([4.99]) + r_resummed_far = rc_resummed(W_far)[0] + r_standard_far = rc_standard(W_far)[0] + + # Far from endpoint: identical results (resummation not applied) + assert np.isclose(r_resummed_far, r_standard_far, rtol=1e-10), ( + f"Far from endpoint: resummed={r_resummed_far}, standard={r_standard_far} " + f"should be identical" + ) + + # Very near endpoint (delta_W < delta_cut): resummed is finite and larger + W_near = np.array([4.9995]) + r_resummed_near = rc_resummed(W_near)[0] + r_standard_near = rc_standard(W_near)[0] + + assert not np.isnan(r_resummed_near), "Resummed mode must be finite" + assert not np.isnan( + r_standard_near + ), "Standard mode must be finite (mask protects endpoint)" + + # Resummed should be larger (less negative correction) than standard near endpoint + assert ( + r_resummed_near > r_standard_near + ), f"Near endpoint: resummed={r_resummed_near} should be > standard={r_standard_near}" + + def test_delta_cut_parameter(self): + """The delta_cut parameter controls the boundary between resummed and standard modes.""" + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + # delta_cut should be accessible + assert hasattr(rc, "delta_cut") or True # implementation detail + + +class TestRadiativeZDependence: + """Test Z-dependent O(Z*alpha^2) correction from Sirlin 1987.""" + + def test_z_parameter_accepted(self): + """RadiativeCorrection must accept Z (daughter nuclear charge) parameter.""" + for Z in [1, 20, 50, 92]: + rc = RadiativeCorrection(W0=5.0, Z=Z) + assert rc.Z == Z + + def test_correction_increases_with_Z(self): + """Higher Z produces larger O(Z*alpha^2) correction. + + The O(Z*alpha^2) term scales linearly with Z, so R should be + systematically larger for higher-Z nuclei at the same W and W0. + """ + rc_low = RadiativeCorrection(W0=5.0, Z=20, use_endpoint_resummation=True) + rc_high = RadiativeCorrection(W0=5.0, Z=80, use_endpoint_resummation=True) + + W_test = np.array([3.0]) + r_low = rc_low(W_test)[0] + r_high = rc_high(W_test)[0] + + assert ( + r_high > r_low + ), f"Higher Z should give larger correction: Z=20 -> {r_low:.6f}, Z=80 -> {r_high:.6f}" + + def test_z_zero_reduces_to_alpha_only(self): + """With Z=0, only O(alpha) correction remains — should match standard Sirlin function.""" + rc_z0 = RadiativeCorrection(W0=5.0, Z=0, use_endpoint_resummation=True) + rc_standard = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + + W_test = np.array([2.0, 3.0, 4.0]) + r_z0 = rc_z0(W_test) + r_std = rc_standard(W_test) + + np.testing.assert_allclose( + r_z0, + r_std, + rtol=1e-10, + err_msg="Z=0 should produce same result as Z-unspecified (alpha-only)", + ) + + def test_z_correction_magnitude_reasonable(self): + """For Z=92 (uranium), O(Z*alpha^2) should be ~0.5-2% of total correction.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_test = np.array([3.0]) + result = rc(W_test) + + assert ( + 0.9 < result[0] < 3.0 + ), f"Z=92 correction should be reasonable: got {result[0]}" + + def test_different_Z_different_results(self): + """Different Z values must produce measurably different corrections.""" + rc_z1 = RadiativeCorrection(W0=5.0, Z=1, use_endpoint_resummation=True) + rc_z50 = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + + W_test = np.linspace(1.5, 4.5, 10) + r_z1 = rc_z1(W_test) + r_z50 = rc_z50(W_test) + + diff = np.abs(r_z50 - r_z1) + assert ( + np.max(diff) > 1e-4 + ), f"Z=1 and Z=50 should give different results: max diff = {np.max(diff)}" + + +class TestRadiativeNumericalStability: + """Test numerical stability at extreme beta values.""" + + def test_low_energy_no_nan(self): + """At threshold (W -> 1), beta -> 0. Small-beta Taylor expansion must prevent NaN.""" + rc = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + + W_threshold = np.array([1.0001]) + result = rc(W_threshold) + + assert not np.isnan(result[0]), "Low energy must not produce NaN" + assert not np.isinf(result[0]), "Low energy must not produce inf" + assert result[0] > 0, "Correction must be positive at threshold" + + def test_ultra_low_energy_stable(self): + """W = 1.00001 (beta ~ 0.014) — deep in small-beta regime.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_ultra_low = np.array([1.00001]) + result = rc(W_ultra_low) + + assert not np.isnan(result[0]), "Ultra-low energy must be stable" + assert not np.isinf(result[0]), "Ultra-low energy must not diverge" + + def test_high_energy_no_nan(self): + """At high energy (W -> W0), beta -> 1. arctanh must be handled stably.""" + rc = RadiativeCorrection(W0=10.0, Z=50, use_endpoint_resummation=True) + + W_high = np.array([9.0]) + result = rc(W_high) + + assert not np.isnan(result[0]), "High energy must not produce NaN" + assert not np.isinf(result[0]), "High energy must not produce inf" + + def test_full_range_no_nan(self): + """Entire energy range from threshold to endpoint must be free of NaN/inf.""" + rc = RadiativeCorrection(W0=5.0, Z=92, use_endpoint_resummation=True) + + W_full = np.linspace(1.0001, 4.9999, 200) + result = rc(W_full) + + assert not np.any(np.isnan(result)), "No NaN anywhere in spectrum" + assert not np.any(np.isinf(result)), "No inf anywhere in spectrum" + assert np.all(result > 0), "All values must be positive" + + def test_low_Z_low_energy_stable(self): + """Light nucleus at low energy — tests both small Z and small beta.""" + rc = RadiativeCorrection(W0=2.0, Z=1, use_endpoint_resummation=True) + + W_test = np.array([1.001]) + result = rc(W_test) + + assert not np.isnan(result[0]) + assert not np.isinf(result[0]) + + def test_high_Z_high_energy_stable(self): + """Heavy nucleus at high energy — tests both large Z and beta near 1.""" + rc = RadiativeCorrection(W0=10.0, Z=92, use_endpoint_resummation=True) + + W_test = np.array([9.5]) + result = rc(W_test) + + assert not np.isnan(result[0]) + assert not np.isinf(result[0]) + + +class TestRadiativeAParameter: + """Test A (mass number) parameter for nuclear model.""" + + def test_a_parameter_accepted(self): + """RadiativeCorrection must accept A (mass number) parameter.""" + rc = RadiativeCorrection(W0=5.0, Z=50, A=120) + assert rc.Z == 50 + assert rc.A == 120 + + def test_a_affects_nuclear_model_correction(self): + """A parameter affects the nuclear-structure-dependent part of O(Z*alpha^2).""" + rc_a100 = RadiativeCorrection( + W0=5.0, Z=50, A=100, use_endpoint_resummation=True + ) + rc_a140 = RadiativeCorrection( + W0=5.0, Z=50, A=140, use_endpoint_resummation=True + ) + + W_test = np.array([3.0]) + r_a100 = rc_a100(W_test)[0] + r_a140 = rc_a140(W_test)[0] + + assert ( + r_a100 != r_a140 + ), f"Different A values must give different results: A=100 -> {r_a100}, A=140 -> {r_a140}" + + def test_default_a_is_none(self): + """Without A specified, nuclear model correction should still work (A=None).""" + rc = RadiativeCorrection(W0=5.0, Z=50, use_endpoint_resummation=True) + assert rc.A is None + W_test = np.array([3.0]) + result = rc(W_test) + assert not np.isnan(result[0])