feat(radiative): implement O(Z*alpha^2) correction and improve numerical stability

- Add Z (nuclear charge) and A (mass number) parameters to RadiativeCorrection
- Implement O(Z*alpha^2) finite nuclear size correction per Sirlin 1987:
  - Energy-dependent part: ln(m_p) - (5/12)*ln(2W) + 43/18
  - Nuclear-structure-dependent part: Delta_F(Z,A) from modified Gaussian model
- Add small-beta Taylor expansion for numerical stability at threshold (beta < 1e-6)
- Extract _safe_beta() and _safe_tanh_inv_beta_over_beta() helper methods
- Remove unused log_wright_bessel import
- Add 14 new tests: Z-dependence, numerical stability, A parameter coverage
- Update docs/06-radiative-corrections.md with new API usage examples

All 70 tests pass. 14 new tests in 3 new test classes.

beta_spectrum/components/radiative.py

@@ -1,40 +1,63 @@
+# 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
+    Implements:
-    divergence handeled by the soft photon resummation of Eq. (53).
+      - 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
+      - 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,
+    ):
         """
         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
         """
         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.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 +65,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 +76,190 @@ class RadiativeCorrection(SpectrumComponent):
         if np.any(mask_finite):
             W_finite = W[mask_finite]
 
-            if self.use_endpoint_resummation:
+            # O(alpha) correction (Sirlin 1967)
-                g_func = self._g_function_resummed(W_finite)
+            delta_1 = self._delta_1(W_finite)
-            else:
-                g_func = self._g_function_standard(W_finite)
 
-            delta_R_finite = self._outer_normalization * (
+            # O(Z*alpha^2) correction (Sirlin 1987)
-                g_func - self._subtraction_term
+            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.
+        O(alpha) universal Sirlin function g(W, W0) per Sirlin (1967) Eq. (20b).
-        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
 
-        # 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
+        # arctanh factor: (tanh^{-1}(beta)/beta) - 1
-        atanh_beta = np.arctanh(beta)
+        atanh_factor = tanh_inv_beta_beta - 1.0
-        atanh_factor = (atanh_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)
+        if self.use_endpoint_resummation:
-        term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + np.log(2.0 * delta_W))
+            # Apply soft-photon resummation from Sirlin (1987)
-
+            # Replace C(beta)*ln(W0-W) with C(beta)*((W0-W)^{t(beta)}-1)/t(beta)
-        # Third line Eq. (50)
+            t_beta = (2.0 * ALPHA / np.pi) * atanh_factor
-        term3 = (atanh_beta / beta) * (
-            2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta
-        )
-
-        # Leading constants from Eq. (50)
-        const_term = -0.75
-
-        return 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):
-
-        t(beta) * ln(W0 - W) -> (W0 - W)^t(beta) - 1
-        """
-        beta = np.sqrt(1.0 - 1.0 / W**2)
-        beta = np.maximum(beta, 1e-12)
-
-        # t(beta) from Eq. (54)
-        atanh_beta = np.arctanh(beta)
-        t_beta = (2.0 * ALPHA / np.pi) * (atanh_beta / beta - 1.0)
-
-        # 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)
-
-        # arctanh factor
-        atanh_factor = (atanh_beta / beta) - 1.0
-
-        # Energy difference
-        delta_W = self.W0 - W
-        delta_W = np.where(delta_W > 0, delta_W, 1e-12)
-
-        # 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
+        else:
+            log_term = np.log(2.0 * delta_W)
 
-        # Term2 with resummed log (instead of ln(2*(W0 - W)))
+        # Energy-dependent term
-        term2_resummed = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term)
+        term2 = 4.0 * atanh_factor * (delta_W / (3.0 * W) - 1.5 + log_term)
 
-        # Term3 (no divergence)
+        # Beta-dependent term (no divergence)
-        term3 = (atanh_beta / beta) * (
+        term3 = tanh_inv_beta_beta * (
-            2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta
+            2.0 * (1.0 + beta**2)
+            + (delta_W**2) / (6.0 * W**2)
+            - 4.0 * tanh_inv_beta_beta
         )
 
-        # Leading constants
+        # Leading constant
         const_term = -0.75
 
-        return const_term + dilog_term + term2_resummed + term3
+        g_func = const_term + dilog_term + term2 + term3
+
+        return self._outer_normalization * (g_func - self._subtraction_term)
+
+    def _delta_2(self, W: np.ndarray) -> np.ndarray:
+        """
+        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
+
+    def _safe_tanh_inv_beta_over_beta(self, beta: np.ndarray) -> np.ndarray:
+        """
+        Compute tanh^{-1}(beta)/beta with numerical stability for small beta.
+
+        For beta < 1e-6, use Taylor expansion:
+          tanh^{-1}(beta)/beta = 1 + beta^2/3 + beta^4/5 + ...
+
+        For beta >= 1e-6, use the direct formula:
+          tanh^{-1}(beta)/beta = 0.5 * ln((1+beta)/(1-beta)) / beta
+
+        Per spec Section 3.3, this avoids division by zero and loss of precision.
+        """
+        result = np.empty_like(beta)
+        small_beta = beta < 1e-6
+
+        # 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
+        )
+
+        # 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 result

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):
+> Implements outer radiative corrections per Sirlin (1967) and Sirlin (1987):
-> - Sirlin's $g(W_0,W)$ with `scipy.special.spence` for dilogarithm
+> - **O(α)**: Universal Sirlin function $g(W,W_0)$ with `scipy.special.spence` for dilogarithm
-> - Soft-photon resummation at endpoint via Eq. (53)
+> - **O(Zα²)**: Finite nuclear size correction $\delta_2(Z,W)$ using modified Gaussian model
-> - Order Zα² terms Δ₁ through Δ₄
+> - 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
+    W0=1.02,
-    use_endpoint_resummation=True  # handles endpoint divergence
+    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)
 ```
 

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, (
+        assert (
-            f"Radiative correction should be near unity: got {result[0]}"
+            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).
@@ -88,9 +89,9 @@ class TestRadiativeEndpointHandling:
         val_close = rc(W_close)[0]
 
         # The correction should decrease (become more negative) near endpoint
-        assert val_close < val_far, (
+        assert (
-            "Standard mode must show divergent behavior near endpoint"
+            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, (
+        assert (
-            f"Correction factor should decrease near endpoint: {r_mid} → {r_near}"
+            r_near <= r_mid
-        )
+        ), f"Correction factor should decrease near endpoint: {r_mid} → {r_near}"
 
 
 class TestRadiativeOutputShape:
@@ -162,6 +163,175 @@ class TestRadiativeResummationSwitch:
 
         # 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])
+        assert not np.isnan(rc_resummed(W_near)[0]), "Resummed mode must be finite"
 
         assert diff_far > 0, "Even far from endpoint there should be some difference"
+
+
+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])