Initial commit

Base modular structure (phase space, Fermi function)
Initial .gitignore

beta_spectrum/__init__.py

@@ -0,0 +1,17 @@
+# beta_spectrum/__init__.py
+
+from .constants import ALPHA
+from .units import T_to_W, W_to_T
+
+from .components.phase_space import PhaseSpace
+from .components.fermi import FermiFunction
+
+__all__ = [
+    "ALPHA",
+    "T_to_W",
+    "W_to_T",
+    "PhaseSpace",
+    "FermiFunction",
+]
+
+__version__ = "0.1.0"

beta_spectrum/base.py

@@ -0,0 +1,20 @@
+# base.py
+from __future__ import annotations
+from abc import ABC, abstractmethod
+import numpy as np
+
+
+ArrayLike = np.ndarray
+
+
+class SpectrumComponent(ABC):
+    """
+    Base class for multiplicative components of beta-spectrum
+    """
+
+    @abstractmethod
+    def __call__(self, W: ArrayLike) -> ArrayLike:
+        """
+        Evaluate component  at total energy W (in m_e units)
+        """
+        pass

beta_spectrum/components/fermi.py

@@ -0,0 +1,52 @@
+# components/fermi.py
+import numpy as np
+from scipy.special import loggamma
+
+from ..base import SpectrumComponent
+from ..constants import ALPHA
+from ..utils import nuclear_radius
+
+
+class FermiFunction(SpectrumComponent):
+    def __init__(self, Z: int, A: int):
+        self.Z = Z
+        self.A = A
+        self.alphaZ = ALPHA * Z
+        self.gamma = np.sqrt(1.0 - self.alphaZ**2)
+        self.R = nuclear_radius(A)
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        p = np.sqrt(W**2 - 1.0)
+
+        # Avoid division by zero
+        p = np.where(p == 0, 1e-30, p)
+
+        eta = self.alphaZ * W / p
+
+        # --- log components ---
+
+        # log prefactor
+        log_prefactor = np.log(2 * (1 + self.gamma))
+
+        # log power term
+        log_power = -2 * (1 - self.gamma) * np.log(2 * p * self.R)
+
+        # log gamma ratio
+        z = self.gamma + 1j * eta
+        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 exponential
+        log_exp = np.pi * eta
+
+        logF = (
+            log_prefactor
+            + log_power
+            + log_gamma_abs_sq
+            - log_gamma_denom
+            + log_exp
+        )
+
+        return np.exp(logF)

beta_spectrum/components/finite_size.py

@@ -0,0 +1,4 @@
+import numpy as np
+
+from beta_spectrum.constants import ALPHA
+from beta_spectrum.base import SpectrumComponent

beta_spectrum/components/phase_space.py

@@ -0,0 +1,19 @@
+# components/phase_space.py
+import numpy as np
+from ..base import SpectrumComponent
+
+
+class PhaseSpace(SpectrumComponent):
+    def __init__(self, W0: float, m_nu: float = 0.0):
+        self.W0 = W0
+        self.m_nu = m_nu
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        p = np.sqrt(W**2 - 1.0)
+        Wv = self.W0 - W
+
+        # neutrino mass generalization
+        if self.m_nu == 0:
+            return p * W * Wv**2
+        else:
+            return p * W * Wv * np.sqrt(Wv**2 - self.m_nu**2)

beta_spectrum/constants.py

@@ -0,0 +1,4 @@
+# constants.py
+ALPHA = 1 / 137.035999084  # fine-structure constant
+HBAR_C_MEV_FM = 197.3269804  # MeV * fm
+ME_MEV = 0.510998950  # MeV

beta_spectrum/spectrum.py

@@ -0,0 +1,17 @@
+# spectrum.py
+import numpy as np
+from typing import List
+from .base import SpectrumComponent
+
+
+class BetaSpectrum:
+    def __init__(self, components: List[SpectrumComponent]):
+        self.components = components
+
+    def __call__(self, W: np.ndarray) -> np.ndarray:
+        result = np.ones_like(W)
+
+        for comp in self.components:
+            result *= comp(W)
+
+        return result

beta_spectrum/units.py

@@ -0,0 +1,12 @@
+# 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

@@ -0,0 +1,9 @@
+# utils.py
+from .constants import HBAR_C_MEV_FM, ME_MEV
+
+
+def nuclear_radius(A: float) -> float:
+    r0_fm = 1.2  # fm
+    R_fm = r0_fm * A ** (1 / 3)
+
+    return R_fm / (HBAR_C_MEV_FM / ME_MEV)