MVP Correction Implemented

Add implementation for the basic set of spectrum components:

Primary
-------
- Phase Space
- Fermi Function

Corrections
-----------

- Atomic:
  - Screening
  - Atomic exchange
- E/M and finite size:
  - Finite nucleus size and charge distribution
  - Radiative

beta_spectrum/__init__.py

@@ -1,17 +1,30 @@
 # beta_spectrum/__init__.py
 
-from .constants import ALPHA
+from .constants import ALPHA, HBAR_C_MEV_FM, ME_MEV, MP_MEV
-from .units import T_to_W, W_to_T
+from .utils import T_to_W, W_to_T, momentum
 
 from .components.phase_space import PhaseSpace
 from .components.fermi import FermiFunction
+from .components.finite_size import FiniteSizeL0, ChargeDistributionU
+from .components.screening import ScreeningCorrection
+from .components.exchange import ExchangeCorrection
+from .components.radiative import RadiativeCorrection
 
 __all__ = [
     "ALPHA",
+    "HBAR_C_MEV_FM",
+    "ME_MEV",
+    "MP_MEV",
     "T_to_W",
     "W_to_T",
+    "momentum",
     "PhaseSpace",
     "FermiFunction",
+    "FiniteSizeL0",
+    "ChargeDistributionU",
+    "ScreeningCorrection",
+    "ExchangeCorrection",
+    "RadiativeCorrection",
 ]
 
-__version__ = "0.1.0"
+__version__ = "0.1.1"

beta_spectrum/components/exchange.py

@@ -0,0 +1,66 @@
+from pathlib import Path
+import numpy as np
+from beta_spectrum.base import SpectrumComponent
+
+
+class ExchangeCorrection(SpectrumComponent):
+    """
+    Atomic xchange correction using Hayen 2018 Table X coefficients (App. G)
+    """
+
+    def __init__(self, Z: int, filename=None, eps: float = 1e-6):
+        self.Z = Z
+        self.eps = eps
+
+        if filename is None:
+            filename = self._default_coeff_path()
+
+        self.coeffs = self._load_coeffs(filename)
+
+    def _load_coeffs(self, filename):
+        data = np.genfromtxt(filename, delimiter=",", names=True)
+
+        for row in data:
+            if int(row["Z"]) == self.Z:
+                return {
+                    "a": row["a"],
+                    "b": row["b"],
+                    "c": row["c"],
+                    "d": row["d"],
+                    "e": row["e"],
+                    "f": row["f"],
+                    "g": row["g"],
+                    "h": row["h"],
+                    "i": row["i"],
+                }
+        raise ValueError(f"No exchange coefficients available for Z={self.Z}")
+
+    def _default_coeff_path(self):
+        return Path(__file__).resolve().parent.parent / "data" / "exchange_coeff.csv"
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        W = np.asarray(W)
+
+        # W' = W - 1 (threshold-sensitive)
+        Wp = np.maximum(W - 1.0, self.eps)
+
+        a = self.coeffs["a"]
+        b = self.coeffs["b"]
+        c = self.coeffs["c"]
+        d = self.coeffs["d"]
+        e = self.coeffs["e"]
+        f = self.coeffs["f"]
+        g = self.coeffs["g"]
+        h = self.coeffs["h"]
+        i = self.coeffs["i"]
+
+        # Terms
+        term1 = a / Wp
+        term2 = b / (Wp**2)
+        term3 = c * np.exp(-d * Wp)
+
+        # Numerical safe (W - f)^g
+        base = np.maximum(W - f, self.eps)
+        term4 = e * np.sin((base**g + h)) / (W**i)
+
+        return 1.0 + term1 + term2 + term3 + term4

beta_spectrum/components/fermi.py

@@ -19,7 +19,7 @@ class FermiFunction(SpectrumComponent):
         p = np.sqrt(W**2 - 1.0)
 
         # Avoid division by zero
-        p = np.where(p == 0, 1e-30, p)
+        p = np.clip(p, 1e-30, None)
 
         eta = self.alphaZ * W / p
 
@@ -36,17 +36,12 @@ class FermiFunction(SpectrumComponent):
         log_gamma_complex = loggamma(z)
         log_gamma_abs_sq = 2 * np.real(log_gamma_complex)
 
-        log_gamma_denom = 2 * loggamma(2 * self.gamma + 1).real
+        log_gamma_denom = 2 * np.real(loggamma(2 * self.gamma + 1))
 
         # log exponential
         log_exp = np.pi * eta
 
-        logF = (
+        logF = log_prefactor + log_power + log_gamma_abs_sq - log_gamma_denom + log_exp
-            log_prefactor
+
-            + log_power
+        return np.exp(logF)
-            + log_gamma_abs_sq
-            - log_gamma_denom
-            + log_exp
-        )
 
-        return np.exp(logF)

beta_spectrum/components/finite_size.py

@@ -2,3 +2,43 @@ import numpy as np
 
 from beta_spectrum.constants import ALPHA
 from beta_spectrum.base import SpectrumComponent
+from beta_spectrum.utils import nuclear_radius
+
+
+class FiniteSizeL0(SpectrumComponent):
+    """
+    Finite nuclear size correction L0(W)
+
+    Implements low-Z expansion (Hayen et al.)
+    Includes (1 + gamma)/2 prefactor
+    """
+
+    def __init__(self, Z: int, A: int):
+        self.Z = Z
+        self.A = A
+
+        self.R = nuclear_radius(A)
+
+        self.alphaZ = ALPHA * Z
+        self.gamma = np.sqrt(1 - self.alphaZ**2)
+
+        # beta-minus decay
+        self.sign = -1
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        # Terms from expansion
+        term1 = self.sign * self.alphaZ * W * self.R
+        term2 = (13 / 60) * self.alphaZ**2
+        term3 = -0.5 * (self.alphaZ * self.R) / W
+
+        return 1 + term1 + term2 + term3
+
+
+class ChargeDistributionU(SpectrumComponent):
+    def __init__(self, Z: int, A: int):
+        self.Z = Z
+        self.R = nuclear_radius(A)
+        self.alphaZ = ALPHA * Z
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        return 1.0 + (1 / 5) * (self.alphaZ * W * self.R) ** 2

beta_spectrum/components/radiative.py

@@ -0,0 +1,148 @@
+import numpy as np
+from scipy.special import spence
+
+from beta_spectrum.base import SpectrumComponent
+from beta_spectrum.constants import MP_MEV, ME_MEV, ALPHA
+
+
+class RadiativeCorrection(SpectrumComponent):
+    """
+    Outer radiative correction R(W, W0) = 1 + delta(W, W0).
+
+    Implements Eq. (47)-(50) from Hayen et al. 2017, with endpoint
+    divergence handeled by the soft photon resummation of Eq. (53).
+
+    This is the INCLUSIVE correction suitable for:
+    - Magnetic spectrometers
+    - Detectors where bremmstrahlung is not detected
+    - Precision at 1% level or worse
+    """
+
+    def __init__(self, W0: float, 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
+        """
+        self.W0 = W0
+        self.use_endpoint_resummation = use_endpoint_resummation
+        self.m_p = MP_MEV / ME_MEV  # ~1836.15
+
+        # Precompute constant terms of Eq. (48)
+        self._outer_normalization = ALPHA / (2.0 * np.pi)  # ~0.00116
+        self._subtraction_term = 3.0 * np.log(self.m_p / (2.0 * self.W0))
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        """
+        Calculate R(W, W0) = 1 + delta_R(w, W0)
+
+        Args:
+            W: Electron total energy array (in natural m_e=c=h_bar=1 units)
+
+        Returns:
+            Radiative correction factor array
+        """
+        W = np.asarray(W)
+
+        # Avoid numerical issues exactly at enpoint W0
+        eps = 1e-9
+        mask_finite = W < (self.W0 - eps)
+
+        delta_R = np.zeros_like(W, dtype=float)
+
+        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)
+
+            delta_R_finite = self._outer_normalization * (
+                g_func - self._subtraction_term
+            )
+            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:
+        """
+        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
+
+        # Dilogarithm term
+        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(arg)
+
+        # arctanh factor
+        atanh_beta = np.arctanh(beta)
+        atanh_factor = (atanh_beta / beta) - 1.0
+
+        # Energy difference
+        delta_W = self.W0 - W
+        delta_W = np.maximum(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))
+
+        # 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
+        )
+
+        # Leading constants from Eq. (50)
+        const_term = 3.0 * np.log(self.m_p) - 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(arg)
+
+        # arctanh factor
+        atanh_factor = (atanh_beta / beta) - 1.0
+
+        # Energy difference
+        delta_W = self.W0 - W
+        delta_W = np.maximum(delta_W, 1e-12)
+
+        # Apply resummation from Eq. (53)
+        # Replace divergent ln(W0 - W) term with (W0 - W)^t(beta) - 1
+        resummed_log_term = (delta_W**t_beta) - 1.0
+
+        # Term2 with resummed log (instead of ln(2*(W0 - W)))
+        term2_resummed = (
+            4.0
+            * atanh_factor
+            * (delta_W / (3.0 * W) - 1.5 + 0.5 * resummed_log_term + np.log(2.0))
+        )
+
+        # Term3 (no divergence)
+        term3 = (atanh_beta / beta) * (
+            2.0 * (1.0 + beta**2) + (delta_W**2) / (6.0 * W**2) - 4.0 * atanh_beta
+        )
+
+        # Leading constants
+        const_term = 3.0 * np.log(self.m_p) - 0.75
+
+        return const_term + dilog_term + term2_resummed + term3

beta_spectrum/components/screening.py

@@ -0,0 +1,97 @@
+import numpy as np
+
+from beta_spectrum.base import SpectrumComponent
+from beta_spectrum.utils import momentum
+from beta_spectrum.constants import ALPHA
+
+
+class ScreeningCorrection(SpectrumComponent):
+    """
+    Atomic screning correction S(Z, W).
+
+    Implements:
+        - Ratio method (non-screened to screened F(Z,W)*p*W) using shifted energy W -> W_tilde
+        - Smooth switching function to regularize low-energy behaviour
+
+    S(W) = 1 + f(W) * (S_raw(W) - 1
+    """
+
+    def __init__(
+        self,
+        fermi_function,
+        V0: float = None,
+        Ws: float = 1.01,
+        Delta: float = 0.005,
+        eps: float = 1e-6,
+        C=0.015,
+    ):
+        """
+        Parameters
+        ----------
+        Z : nuclear charge
+        fermi_function : callable
+            Instance of FermiFunction class
+        V0 : screening potential (in m_e units)
+            If None, uses simple scaling law
+        Ws : switching midpoint
+        Delta : switching width
+        eps : numerical floor for W_tilde
+        """
+        self.fermi = fermi_function
+        self.Z = fermi_function.Z
+        self.Ws = Ws
+        self.Delta = Delta
+        self.eps = eps
+        self._C = C
+        self.V0 = V0 if V0 is not None else self._estimate_V0()
+
+    # --------------------
+    # Core evaluation
+    # --------------------
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        W = np.asarray(W)
+
+        # Energy shift by screening potential V0
+        W_tilde = self._screened_energy(W)
+
+        # Momenta
+        p = momentum(W)
+        p_tilde = momentum(W_tilde)
+
+        # Fermi functions
+        F = self.fermi(W)
+        F_tilde = self.fermi(W_tilde)
+
+        # Raw ratio model (for W >> V0)
+        S_raw = (F_tilde / F) * (p_tilde * W_tilde) / (p * W)
+
+        # Smooth switching
+        f = self._switching_function(W)
+
+        return 1.0 + f * (S_raw - 1.0)
+
+    # --------------------
+    # Helpers
+    # --------------------
+    def _screened_energy(self, W):
+        """
+        Apply energy shift W -> W_tilde with safe numerical floor
+        """
+        W_tilde = W - self.V0
+        return np.maximum(W_tilde, 1.0 + self.eps)
+
+    def _switching_function(self, W):
+        """
+        Logistic switching function
+        """
+        return 1.0 / (1.0 + np.exp((self.Ws - W) / self.Delta))
+
+    def _estimate_V0(self):
+        """
+        Simple scaling estimate for screening potential.
+
+        V0 ~ alpha * Z^(4/3) * C
+
+        C ~ 0.01-0.02 (tunable constant)
+        """
+        return ALPHA * self.Z ** (4.0 / 3.0) * self._C

beta_spectrum/constants.py

@@ -1,4 +1,5 @@
 # constants.py
 ALPHA = 1 / 137.035999084  # fine-structure constant
 HBAR_C_MEV_FM = 197.3269804  # MeV * fm
-ME_MEV = 0.510998950  # MeV
+ME_MEV = 0.510998950  # MeV / c^2
+MP_MEV = 938.272088816  # MeV / c^2

beta_spectrum/units.py

@@ -1,12 +0,0 @@
-# units.py
-import numpy as np
-
-ME_KEV = 510.99895000  # keV electron mass
-
-
-def W_to_T(W: float | np.ndarray) -> np.ndarray:
-    return (W - 1.0) * ME_KEV
-
-
-def T_to_W(T_keV: float | np.ndarray) -> np.ndarray:
-    return T_keV / ME_KEV + 1.0

beta_spectrum/utils.py

@@ -1,9 +1,27 @@
 # utils.py
+import numpy as np
+
 from .constants import HBAR_C_MEV_FM, ME_MEV
 
 
-def nuclear_radius(A: float) -> float:
+def nuclear_radius(A: float, r0: float = 1.2) -> float:
     r0_fm = 1.2  # fm
     R_fm = r0_fm * A ** (1 / 3)
 
     return R_fm / (HBAR_C_MEV_FM / ME_MEV)
+
+
+def momentum(W):
+    return np.sqrt(W**2 - 1)
+
+
+def neutrino_energy(W, W0):
+    return W0 - W
+
+
+def W_to_T(W: float | np.ndarray) -> np.ndarray:
+    return (W - 1.0) * ME_MEV
+
+
+def T_to_W(T_keV: float | np.ndarray) -> np.ndarray:
+    return T_keV / ME_MEV + 1.0

pyproject.toml

@@ -4,10 +4,10 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "beta-spectrum"
-version = "0.1.0"
+version = "0.1.1"
 description = "High-precision beta spectrum calculation for nuclear physics"
 authors = [
-    { name = "Your Name" }
+    { name = "Evgeniy Unzhakov" }
 ]
 requires-python = ">=3.10"