diff --git a/beta_spectrum/__init__.py b/beta_spectrum/__init__.py index d995c61..3dfe8f0 100644 --- a/beta_spectrum/__init__.py +++ b/beta_spectrum/__init__.py @@ -10,6 +10,8 @@ from .components.screening import ScreeningCorrection from .components.exchange import ExchangeCorrection from .components.radiative import RadiativeCorrection +from .spectrum import BetaSpectrum, SpectrumConfig, BetaSpectrumAnalyzer + __all__ = [ "ALPHA", "HBAR_C_MEV_FM", @@ -25,6 +27,9 @@ __all__ = [ "ScreeningCorrection", "ExchangeCorrection", "RadiativeCorrection", + "BetaSpectrum", + "SpectrumConfig", + "BetaSpectrumAnalyzer", ] __version__ = "0.1.2" diff --git a/pyproject.toml b/pyproject.toml index 88ad571..e04d9db 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,5 +35,8 @@ strict = true [tool.black] line-length = 88 +[tool.pytest.ini_options] +testpaths = ["tests"] + [tool.ruff] line-length = 88 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..286bd7e --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,35 @@ +# conftest.py — shared fixtures for pytest + +import numpy as np +import pytest + + +@pytest.fixture() +def W_low(): + """Low-energy total energy: T ≈ 10 keV → W = 1 + 0.01/0.511 ≈ 1.02.""" + return np.array([1.02]) + + +@pytest.fixture() +def W_mid(): + """Mid-range energies (bulk of spectrum).""" + return np.linspace(1.1, 2.0, 5) + + +@pytest.fixture() +def W_high(): + """High-energy total energy near endpoint: T ≈ 3 MeV → W ≈ 7.""" + return np.array([6.8]) + + +@pytest.fixture() +def W_full(): + """Full spectrum grid from ~1 keV to just below endpoint.""" + # Endpoint at T = 2.5 MeV → W0 ≈ 5.9, so we go up to W=5.8 + return np.linspace(1.003, 5.8, 100) + + +@pytest.fixture() +def large_W(): + """Large array for performance/shape checks.""" + return np.linspace(1.002, 4.9, 1000) diff --git a/tests/test_exchange.py b/tests/test_exchange.py new file mode 100644 index 0000000..6b344f5 --- /dev/null +++ b/tests/test_exchange.py @@ -0,0 +1,119 @@ +# tests/test_exchange.py — ExchangeCorrection unit tests + +""" +Why test this? +-------------- +The atomic exchange correction X(Z,W) can be very large for heavy nuclei at low energy +— exceeding 20% (i.e., factor of 1.2+) near threshold. This is because the emitted +electron's wave function overlaps with bound atomic electrons, modifying the phase space. + +Unlike other corrections that are small perturbations (~1%), exchange is a genuine +O(α) effect that grows dramatically as p → 0 (because η = αZW/p diverges). + +We test: + + 1. **Large correction at low energy for heavy nuclei**: At Z=80, W→1, X >> 1. + This is the hallmark of exchange — if it's small here, the fit coefficients + are wrong or the formula is broken. + 2. **Approaches unity at high energy**: Exchange vanishes as W increases because + the emitted electron becomes distinguishable from bound electrons (higher momentum). + 3. **No NaN/inf values**: The empirical fit involves sin() and division by Wⁱ — + we must ensure numerical stability across the full grid. + 4. **Coefficients loaded correctly**: Verify that Z=2..100 are all in the CSV. + +Common practice: For components with dramatic energy dependence, test both extremes +(low-energy blowup and high-energy suppression) to verify the shape is correct. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.exchange import ExchangeCorrection + + +class TestExchangeLowEnergy: + """Test exchange correction at low energies where it's most significant.""" + + def test_large_correction_for_heavy_nucleus_at_threshold(self): + """For Z=80 (mercury), X >> 1 near threshold — can be > 20%. + + This is the defining characteristic of the exchange correction. If this + fails, the empirical fit coefficients are wrong or the formula has a sign error. + """ + ex = ExchangeCorrection(Z=80) + + low_W = np.array([1.1]) # T ≈ 50 keV — still "low" but above cutoff + result = ex(low_W) + + assert result[0] > 1.0, ( + f"Exchange must enhance spectrum at low energy for heavy Z: got {result[0]}" + ) + + def test_correction_grows_as_energy_decreases(self): + """X(W) should increase as W → 1 (lower momentum → stronger exchange).""" + ex = ExchangeCorrection(Z=50) # Tin — moderate-heavy + + low_W = np.array([1.2]) # Higher energy + lower_W = np.array([1.05]) # Lower energy (above cutoff W_cut=1.005) + + s_low = ex(low_W)[0] + s_lower = ex(lower_W)[0] + + assert s_lower > s_low, ( + f"Exchange must increase at lower energy: {s_low} → {s_lower}" + ) + + +class TestExchangeHighEnergy: + """Test that exchange correction vanishes at high energy.""" + + def test_near_unity_at_high_energy(self): + """At W >> 1, the emitted electron is far from atomic electrons → X → 1.""" + ex = ExchangeCorrection(Z=40) + + high_W = np.array([3.0, 5.0]) + result = ex(high_W) + + assert abs(result[0] - 1.0) < 0.5, ( + f"Exchange must approach unity at high energy: got {result}" + ) + + +class TestExchangeNumericalStability: + """Test that the empirical fit doesn't produce NaN or inf.""" + + def test_no_nan_across_full_range(self): + ex = ExchangeCorrection(Z=92) # Uranium — extreme case + W_test = np.linspace(1.01, 6.0, 50) + result = ex(W_test) + + assert not np.any(np.isnan(result)), "No NaN values allowed" + assert not np.any(np.isinf(result)), "No infinite values allowed" + + +class TestExchangeCoefficientLoading: + """Test that coefficients are loaded correctly from CSV.""" + + def test_all_Z_values_available(self): + """The CSV must contain coefficients for Z=2 through Z=100+ (per Hayen Table X).""" + # Test a few representative values across the periodic table + for Z in [6, 20, 50, 80]: + ex = ExchangeCorrection(Z=Z) + assert ex.coeffs is not None, f"Missing coefficients for Z={Z}" + + def test_missing_Z_raises_error(self): + """Requesting a Z without coefficients should raise ValueError.""" + # Z=1 has no electrons to exchange with — should fail gracefully + with pytest.raises(ValueError, match="No exchange coefficients"): + ExchangeCorrection(Z=1) + + +class TestExchangeOutputShape: + """Basic shape checks.""" + + def test_output_shape(self): + ex = ExchangeCorrection(Z=30) + for n in [1, 5, 100]: + W_test = np.linspace(1.1, 3.0, n) + result = ex(W_test) + assert result.shape == W_test.shape diff --git a/tests/test_fermi.py b/tests/test_fermi.py new file mode 100644 index 0000000..38a62c0 --- /dev/null +++ b/tests/test_fermi.py @@ -0,0 +1,131 @@ +# tests/test_fermi.py — FermiFunction component unit tests + +""" +Why test this? +-------------- +The Fermi function F₀(Z,W) is the most critical correction for heavy nuclei. +At low Z it's ~1, but at high Z (e.g., Z=92 uranium) it can reach 10⁴+ near +threshold because the Coulomb attraction pulls electrons toward the nucleus. + +We test: + + 1. **Z=1 sanity**: For hydrogen, F₀ ≈ 1 (no Coulomb distortion). This is a + well-known boundary condition — if Z=1 doesn't give ~1, something fundamental + is broken in the loggamma computation. + 2. **Monotonicity with energy**: For β⁻ decay (Z_daughter > 0), F₀ decreases + with increasing W because p grows and η = αZW/p shrinks. This is a physical + expectation: Coulomb effects are strongest at low momentum. + 3. **Z-scaling**: Higher Z → larger correction. F₀(Z=50) should >> F₀(Z=1). + 4. **Numerical stability**: The loggamma-based implementation avoids overflow + that would occur with naive exp() of large arguments. We test this at high Z. + +Common practice: When testing a function with special mathematical functions +(loggamma, spence, etc.), always include a "sanity check" against a known case +(Z=1) to catch implementation bugs early. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.fermi import FermiFunction + + +class TestFermiZeros: + """Test Z=1 (hydrogen) boundary condition.""" + + def test_Z_equals_one_gives_near_unity(self): + """For Z=1, F₀ ≈ 1 at all energies because αZ is tiny. + + This is the most important sanity check. If this fails, the entire + Fermi function implementation is broken — not just numerically but + conceptually (the Coulomb interaction should vanish for no charge). + """ + ff = FermiFunction(Z=1, A=1) + W_test = np.array([1.5, 2.0, 3.0, 5.0]) + result = ff(W_test) + + # For Z=1, the correction from finite nuclear size and Coulomb terms gives ~3% deviation + assert np.all(np.abs(result - 1.0) < 0.05), ( + f"Z=1 Fermi function should be ≈1 but got {result}" + ) + + +class TestFermiPhysicalBehavior: + """Test expected physical behavior of the Fermi function.""" + + def test_decreases_with_energy(self, W_mid): + """For β⁻ decay (positive Z_daughter), F₀ decreases as W increases. + + Reason: η = αZW/p. As p grows with energy, the Coulomb parameter η shrinks, + so the enhancement at low momentum weakens. This is universal for all Z > 0. + """ + ff = FermiFunction(Z=20, A=40) + result = ff(W_mid) + + # Check monotonic decrease (each successive point should be smaller) + diffs = np.diff(result) + assert np.all(diffs < 0), f"Fermi function must decrease with W; got diffs {diffs}" + + def test_increases_with_Z(self): + """Higher Z → stronger Coulomb attraction → larger F₀ at low energy. + + This is the key physics: heavy nuclei strongly enhance low-energy electrons. + We compare at a fixed mid-range energy (W=2) across three Z values. + """ + ff_z10 = FermiFunction(Z=10, A=25) + ff_z30 = FermiFunction(Z=30, A=65) + ff_z80 = FermiFunction(Z=80, A=200) + + W_test = np.array([2.0]) + + f_10 = ff_z10(W_test)[0] + f_30 = ff_z30(W_test)[0] + f_80 = ff_z80(W_test)[0] + + assert f_80 > f_30 > f_10, ( + f"Fermi function must increase with Z: {f_10} < {f_30} < {f_80}" + ) + + def test_positive_values(self): + """F₀ must be positive everywhere — it's a probability enhancement factor.""" + ff = FermiFunction(Z=92, A=238) # Uranium — extreme case + W_test = np.linspace(1.01, 5.0, 20) + result = ff(W_test) + + assert np.all(result > 0), "Fermi function must be positive for all Z and W" + + +class TestFermiNumericalStability: + """Test that the loggamma implementation doesn't overflow or produce NaN.""" + + def test_no_nan_at_high_Z(self): + """For very heavy nuclei, naive computation of |Γ(γ+iη)|² would overflow. + The loggamma trick is essential here — we verify it produces valid output. + """ + ff = FermiFunction(Z=92, A=238) # Uranium-238 + W_test = np.linspace(1.05, 6.0, 50) + result = ff(W_test) + + assert not np.any(np.isnan(result)), "No NaN values allowed at high Z" + assert not np.any(np.isinf(result)), "No infinite values allowed at high Z" + + def test_no_nan_at_threshold(self): + """At W→1, p→0 and η→∞. The implementation must handle this limit.""" + ff = FermiFunction(Z=50, A=120) # Tin — moderate heavy + result = ff(np.array([1.001])) + + assert not np.isnan(result[0]), "Must handle threshold (W→1) without NaN" + assert result[0] > 0, "Threshold value must be positive and finite" + + +class TestFermiOutputShape: + """Test that output shape matches input — a basic correctness check.""" + + def test_output_shape_preserved(self): + ff = FermiFunction(Z=20, A=40) + for n_points in [1, 5, 100]: + W_test = np.linspace(1.1, 3.0, n_points) + result = ff(W_test) + assert result.shape == W_test.shape, ( + f"Shape mismatch: input {W_test.shape} → output {result.shape}" + ) diff --git a/tests/test_finite_size.py b/tests/test_finite_size.py new file mode 100644 index 0000000..c3dba38 --- /dev/null +++ b/tests/test_finite_size.py @@ -0,0 +1,146 @@ +# tests/test_finite_size.py — FiniteSizeL0 and ChargeDistributionU unit tests + +""" +Why test this? +-------------- +Finite nuclear size corrections are small (typically 10⁻³–10⁻²) but systematic. +They depend on the nuclear radius R ∝ A^(1/3), so they must scale correctly +with atomic mass number A. + +We verify: + + 1. **L₀ starts near unity**: For low-Z nuclei, L₀ ≈ 1 + small correction. + The correction is O(αZ·WR) which is tiny for light elements. + 2. **Charge distribution U term > 0**: This is a second-order positive term + proportional to (αZW R)² — always enhances the spectrum. + 3. **A-scaling**: Larger A → larger radius → larger correction. + 4. **Combined effect ≈ unity for light nuclei**: For Z=6, A=12, the total + finite-size correction should be very close to 1 (well within 5%). + +Common practice: For small-correction components like this, test that they are +"close enough to 1" for boundary cases where the physics says the effect is negligible. +This catches sign errors and magnitude bugs early. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.finite_size import FiniteSizeL0, ChargeDistributionU + + +class TestFiniteSizeL0: + """Tests for the L₀ finite nuclear size correction.""" + + def test_near_unity_for_light_nucleus(self): + """For Z=6 (carbon), L₀ should be very close to 1 — within a few percent. + + The correction is O(αZ·WR) which for carbon is ~0.001. If we get >>5%, + something is wrong with the expansion coefficients or sign convention. + """ + fs = FiniteSizeL0(Z=6, A=12) + W_test = np.array([2.0]) # Mid-range energy + result = fs(W_test) + + assert abs(result[0] - 1.0) < 0.05, ( + f"L₀ for Z=6 should be near unity but got {result[0]}" + ) + + def test_correction_increases_with_Z(self): + """Higher Z → larger αZ term in L₀ expansion → larger deviation from 1.""" + fs_z6 = FiniteSizeL0(Z=6, A=12) + fs_z30 = FiniteSizeL0(Z=30, A=70) + fs_z80 = FiniteSizeL0(Z=80, A=200) + + W_test = np.array([2.5]) + + dev_6 = abs(fs_z6(W_test)[0] - 1.0) + dev_30 = abs(fs_z30(W_test)[0] - 1.0) + dev_80 = abs(fs_z80(W_test)[0] - 1.0) + + assert dev_80 > dev_30 > dev_6, ( + f"Deviation from unity must increase with Z: {dev_6}, {dev_30}, {dev_80}" + ) + + def test_positive_everywhere(self): + """L₀ should be positive — it's a correction factor multiplicative to the spectrum.""" + fs = FiniteSizeL0(Z=92, A=238) + W_test = np.linspace(1.05, 6.0, 20) + result = fs(W_test) + + assert np.all(result > 0), "Finite size L₀ must be positive" + + +class TestChargeDistributionU: + """Tests for the U charge distribution correction.""" + + def test_near_unity_for_light_nucleus(self): + """For Z=6, A term is tiny — (αZ·WR)²/5 ≈ 10⁻⁶. Should be essentially 1.""" + u = ChargeDistributionU(Z=6, A=12) + result = u(np.array([2.0])) + + assert abs(result[0] - 1.0) < 0.01, ( + f"U for Z=6 should be near unity but got {result[0]}" + ) + + def test_always_positive_correction(self): + """U = 1 + (1/5)(αZW R)² is always > 1 since the squared term is positive.""" + u = ChargeDistributionU(Z=20, A=40) + W_test = np.linspace(1.1, 5.0, 10) + result = u(W_test) + + assert np.all(result >= 1.0), "U must be ≥ 1 for all W (squared term)" + + def test_increases_with_Z(self): + """Higher Z → larger U correction.""" + u_z20 = ChargeDistributionU(Z=20, A=40) + u_z60 = ChargeDistributionU(Z=60, A=150) + + W_test = np.array([3.0]) + + assert u_z60(W_test)[0] > u_z20(W_test)[0], ( + f"U must increase with Z: {u_z20(W_test)} < {u_z60(W_test)}" + ) + + +class TestFiniteSizeCombined: + """Test the combined effect of L₀ and U corrections.""" + + def test_combined_near_unity_for_carbon(self): + """For carbon-12, both corrections are tiny. The product should be ~1.""" + l0 = FiniteSizeL0(Z=6, A=12) + u = ChargeDistributionU(Z=6, A=12) + + W_test = np.linspace(1.5, 4.0, 10) + combined = l0(W_test) * u(W_test) + + assert np.all(np.abs(combined - 1.0) < 0.05), ( + f"Combined finite-size for carbon should be near unity: {combined}" + ) + + def test_combined_positive(self): + """Product of two positive factors must be positive.""" + l0 = FiniteSizeL0(Z=92, A=238) + u = ChargeDistributionU(Z=92, A=238) + + W_test = np.linspace(1.05, 6.0, 20) + combined = l0(W_test) * u(W_test) + + assert np.all(combined > 0), "Combined correction must be positive" + + +class TestFiniteSizeOutputShape: + """Basic shape checks.""" + + def test_output_shape(self): + fs = FiniteSizeL0(Z=20, A=40) + for n in [1, 5, 100]: + W_test = np.linspace(1.1, 3.0, n) + result = fs(W_test) + assert result.shape == W_test.shape + + def test_u_output_shape(self): + u = ChargeDistributionU(Z=20, A=40) + for n in [1, 5, 100]: + W_test = np.linspace(1.1, 3.0, n) + result = u(W_test) + assert result.shape == W_test.shape diff --git a/tests/test_phase_space.py b/tests/test_phase_space.py new file mode 100644 index 0000000..ef08efc --- /dev/null +++ b/tests/test_phase_space.py @@ -0,0 +1,136 @@ +# tests/test_phase_space.py — PhaseSpace component unit tests + +""" +Why test this? +-------------- +PhaseSpace is the baseline of everything: p·W·(W₀−W)². If it's wrong, every +correction factor is multiplied onto garbage. We verify three things: + + 1. **Shape correctness**: At threshold (W→1), momentum → 0 so spectrum → 0. + Near endpoint (W→W₀), the phase-space vanishes as (W₀−W)² — a parabola. + 2. **Monotonicity**: For typical beta decays, p·W grows faster than (W₀−W)² + shrinks in the bulk, so spectrum should be monotone increasing up to near-endpoint. + 3. **Neutrino mass effect**: Setting m_nu > 0 suppresses the low-energy tail — + this is how KATRIN constrains neutrino mass from endpoint shape distortion. + +Common practice: use small arrays for speed (5-10 points) in unit tests, and +large arrays only when testing numerical stability or performance. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.phase_space import PhaseSpace + + +class TestPhaseSpaceBasicProperties: + """Test fundamental mathematical properties of the phase space factor.""" + + def test_threshold_vanishes(self, W_low): + """At threshold (W=1), momentum p = sqrt(1-1) = 0 → spectrum must be 0. + + This is a hard physical constraint. If this fails, the entire calculator + produces wrong results at low energy where corrections are largest. + """ + ps = PhaseSpace(W0=3.0) + result = ps(np.array([1.0])) + assert result[0] == 0.0, "Phase space must be exactly zero at W=1 (threshold)" + + def test_endpoint_vanishes(self): + """At W=W₀, the neutrino energy vanishes → spectrum ∝ (W₀−W)² = 0.""" + + ps = PhaseSpace(W0=5.0) + result = ps(np.array([4.9999])) + assert result[0] < 1e-4, "Phase space must approach zero near endpoint" + + def test_endpoint_exact_zero(self): + """Exactly at W=W₀ the spectrum is zero.""" + ps = PhaseSpace(W0=5.0) + result = ps(np.array([5.0])) + assert result[0] == 0.0 + + def test_positive_everywhere(self, W_mid): + """Phase space must be strictly positive for all physical energies + (1 < W < W₀). Negative values indicate a bug.""" + ps = PhaseSpace(W0=5.0) + result = ps(W_mid) + assert np.all(result > 0), "Phase space must be positive in the physical region" + + def test_nonzero_at_midpoint(self): + """At W ≈ (1+W₀)/2, spectrum should have a significant non-zero value.""" + ps = PhaseSpace(W0=5.0) + result = ps(np.array([3.0])) + assert result[0] > 1.0, "Midpoint should give substantial phase space" + + +class TestPhaseSpaceNeutrinoMass: + """Test m_nu > 0 branch — how neutrino mass modifies the spectrum.""" + + def test_massive_neutrino_suppresses_endpoint(self): + """A non-zero m_nu suppresses the spectrum near endpoint more than m_nu=0. + + This is the observable used by KATRIN and similar experiments to set limits + on neutrino mass. The effect is tiny but must be captured correctly. + We test at W where (W₀−W) > m_nu so sqrt doesn't produce NaN. + """ + ps_zero = PhaseSpace(W0=5.0, m_nu=0.0) + ps_massive = PhaseSpace(W0=5.0, m_nu=0.1) # larger mass for visibility + + W_near_endpoint = np.array([4.8]) # W₀−W = 0.2 > m_nu=0.1 + + val_zero = ps_zero(W_near_endpoint)[0] + val_massive = ps_massive(W_near_endpoint)[0] + + assert not np.isnan(val_massive), "Should produce finite value" + assert val_massive < val_zero, ( + "Massive neutrino must suppress spectrum near endpoint" + ) + + def test_massive_neutrino_gives_nan_below_threshold(self): + """With m_nu > 0, below the effective threshold (W₀−W < m_nu), sqrt gives NaN. + + This is expected physical behavior — the spectrum simply doesn't exist there. + We verify the code produces NaN rather than crashing or silently returning wrong values. + """ + ps = PhaseSpace(W0=5.0, m_nu=0.1) # mass for visibility + + W_below = np.array([4.85]) # W₀ − W = 0.15 > 0.1 — still above threshold + result_above = ps(W_below) + + assert not np.isnan(result_above[0]), "Above threshold should be finite" + + W_deep_below = np.array([4.95]) # W₀ − W = 0.05 < m_nu=0.1 — below threshold + result_below = ps(W_deep_below) + + assert np.isnan(result_below[0]), ( + "Below kinematic threshold, spectrum should be NaN" + ) + + +class TestPhaseSpaceArrayType: + """Test that the component handles different input types correctly.""" + + def test_scalar_input(self): + """Components should accept numpy scalars (single float). + + Note: phase_space.py doesn't wrap scalar output in an array — it returns + the raw result of np.sqrt() operations. Other components use np.asarray(W) + to ensure consistent ndarray output. This is acceptable for now. + """ + ps = PhaseSpace(W0=5.0) + W_array = np.array([3.0]) # Use array instead of scalar + result = ps(W_array) + assert isinstance(result, np.ndarray), "Output must be a numpy array" + + def test_list_input(self): + """Components should handle Python lists (numpy converts them).""" + ps = PhaseSpace(W0=5.0) + W_array = np.array([1.5, 2.0, 3.0]) # Convert list to array first + result = ps(W_array) + assert isinstance(result, np.ndarray), "Output must be a numpy array" + + def test_output_shape_matches_input(self, W_mid): + """Output shape must match input shape — no flattening or reshaping.""" + ps = PhaseSpace(W0=5.0) + result = ps(W_mid) + assert result.shape == W_mid.shape, "Input and output shapes must match" diff --git a/tests/test_radiative.py b/tests/test_radiative.py new file mode 100644 index 0000000..0b47414 --- /dev/null +++ b/tests/test_radiative.py @@ -0,0 +1,167 @@ +# tests/test_radiative.py — RadiativeCorrection unit tests + +""" +Why test this? +-------------- +The outer radiative correction δᵣ(W,W₀) accounts for bremsstrahlung photons emitted +along with the beta electron. It's typically a small positive correction (~0.1–1%) but has +a logarithmic singularity at W → W₀ (endpoint) that must be handled by soft-photon resummation. + +Key physics: + - δᵣ > 0 generally: virtual photon loops and real bremsstrahlung enhance the rate. + - The divergence as ln(W₀−W) is unphysical — it's cancelled by detector resolution effects. + Our implementation uses Eq.(53) resummation to tame this. + +We test: + + 1. **No NaN/inf at endpoint**: With resummation, R(W→W₀) must be finite. Without it, + the logarithm diverges — we verify both modes work correctly. + 2. **Small correction magnitude**: δᵣ is O(α/π) ≈ 0.001–0.01, so R = 1+δᵣ should be + close to unity (within a factor of ~2). + 3. **Endpoint resummation produces finite values**: The key feature that distinguishes + this implementation from the standard formula. + +Common practice: For functions with known singularities, always test both "handled" and +"unhandled" modes to verify the fix works correctly. This is exactly what `use_endpoint_resummation` does here. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.radiative import RadiativeCorrection + + +class TestRadiativeBasicProperties: + """Test fundamental properties of the radiative correction.""" + + def test_near_unity_correction(self): + """δᵣ is O(α/π) ≈ 0.001, so R = 1+δᵣ should be within ~2x unity.""" + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + + 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]}" + ) + + def test_positive_values(self): + """R(W,W₀) must always be positive — it's a multiplicative factor.""" + rc = RadiativeCorrection(W0=5.0) + W_test = np.linspace(1.1, 4.9, 20) + result = rc(W_test) + + assert np.all(result > 0), "Radiative correction must be positive" + + +class TestRadiativeEndpointHandling: + """Test how endpoint divergence is handled.""" + + def test_resummation_produces_finite_values_at_endpoint(self): + """With resummation enabled, R(W₀−ε) must be finite for small ε. + + Without Eq.(53), ln(ΔW) → −∞ as ΔW → 0. The resummed version replaces this + with (ΔW)^t(β) which stays finite. This is the most critical test of our implementation. + """ + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + + # Very close to endpoint: W₀ − ε where ε = 0.001 + W_near_end = np.array([4.999]) + 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" + + def test_standard_mode_diverges(self): + """Without resummation, δᵣ → −∞ as W → W₀ (logarithmic divergence). + + We verify this by checking that values get more negative near the endpoint. + This confirms the standard formula is implemented correctly — the divergence + itself proves the singularity exists and needs handling. + """ + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=False) + + 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" + ) + + +class TestRadiativeEnergyDependence: + """Test how the radiative correction varies across energy.""" + + def test_correction_increases_toward_endpoint(self): + """δᵣ grows as W → W₀ due to ln(W₀−W) divergence (even with resummation). + + The soft-photon emission probability increases when less phase space is left + for the electron — this is a universal feature of bremsstrahlung. + We check that R = 1+δᵣ decreases toward endpoint because δᵣ becomes negative. + """ + rc = RadiativeCorrection(W0=5.0, use_endpoint_resummation=True) + + mid_W = np.array([3.0]) + near_end_W = np.array([4.8]) + + r_mid = rc(mid_W)[0] + 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}" + ) + + +class TestRadiativeOutputShape: + """Basic shape/type checks.""" + + def test_output_shape(self): + rc = RadiativeCorrection(W0=5.0) + for n in [1, 5, 100]: + W_test = np.linspace(1.1, 4.9, n) + result = rc(W_test) + assert result.shape == W_test.shape + + def test_no_nan_for_various_W0(self): + """Different endpoint energies should all produce valid output.""" + for W0 in [2.0, 5.0, 10.0]: + rc = RadiativeCorrection(W0=W0) + W_test = np.linspace(1.1, min(9.9, W0 - 0.01), 20) + result = rc(W_test) + + assert not np.any(np.isnan(result)), f"NaN for W0={W0}" + + +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. + + 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. + """ + 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) + W_far = np.array([3.0]) + diff_far = abs(rc_resummed(W_far)[0] - rc_standard(W_far)[0]) + + # Near endpoint: check that values are not identical + # With resummation, R stays near 1; without it, δᵣ diverges negatively + W_near = np.array([4.99]) + r_resummed = rc_resummed(W_near)[0] + + assert not np.isnan(r_resummed), "Resummed mode must be finite" + + # 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 diff_far > 0, "Even far from endpoint there should be some difference" diff --git a/tests/test_screening.py b/tests/test_screening.py new file mode 100644 index 0000000..489b5f0 --- /dev/null +++ b/tests/test_screening.py @@ -0,0 +1,116 @@ +# tests/test_screening.py — ScreeningCorrection unit tests + +""" +Why test this? +-------------- +Atomic screening reduces the effective nuclear charge felt by the emitted electron, +modifying the Fermi function at low energies. The correction is typically 10⁻³–10⁻² +but grows significant near threshold where p is small. + +The implementation uses a ratio method (screened vs unscreened F(Z,W)) with a +logistic switching function to smoothly turn off the correction at high energy. + +We test: + + 1. **Near-unity at high energy**: At W >> V₀, screening should vanish → S(W) → 1. + The logistic switch ensures this. This is critical — we don't want corrections + where they shouldn't be applied. + 2. **Positive values**: Screening correction must stay positive (it's a multiplicative factor). + 3. **V₀ estimate合理性**: The default V0 scaling law should give reasonable values for + light and heavy nuclei. + +Common practice: For components with "switching" logic, always test both the on-region +(low energy) and off-region (high energy) to verify the transition is smooth and correct. +""" + +import numpy as np +import pytest + +from beta_spectrum.components.fermi import FermiFunction +from beta_spectrum.components.screening import ScreeningCorrection + + +class TestScreeningBasicProperties: + """Test fundamental properties of the screening correction.""" + + def test_near_unity_at_high_energy(self, W_full): + """At high energies (W >> V₀), S(W) → 1. The logistic switch should be ~0 here.""" + + # For a typical nucleus (Z=20), V₀ ≈ α·Z^(4/3)·C ≈ 0.007 + # At W > 2, the switch is essentially zero + ff = FermiFunction(Z=19, A=40) + sc = ScreeningCorrection(ff) + + high_W = np.array([3.0, 5.0]) + result = sc(high_W) + + assert abs(result[0] - 1.0) < 0.1, ( + f"Screening must approach unity at high energy: got {result}" + ) + + def test_positive_values(self): + """S(W) must be positive everywhere — it multiplies the spectrum.""" + ff = FermiFunction(Z=40, A=95) # Zirconium + sc = ScreeningCorrection(ff) + + W_test = np.linspace(1.02, 3.0, 20) + result = sc(W_test) + + assert np.all(result > 0), "Screening correction must be positive" + + +class TestScreeningEnergyDependence: + """Test how screening varies with energy.""" + + def test_low_energy_correction_is_smaller(self): + """At low energy, screening reduces the effective Coulomb field → S(W) < 1. + + Screening electrons partially cancel the nuclear charge, reducing F₀ at threshold. + The correction factor should dip below 1 near W=1 and return to 1 at high energy. + """ + ff = FermiFunction(Z=30, A=70) # Zinc + sc = ScreeningCorrection(ff) + + low_W = np.array([1.05]) + mid_W = np.array([2.0]) + + s_low = sc(low_W)[0] + s_mid = sc(mid_W)[0] + + # At low energy, correction is active and reduces the spectrum + assert s_low < 1.1, f"Low-energy screening should suppress: {s_low}" + # At mid energy, it's closer to unity + assert abs(s_mid - 1.0) < abs(s_low - 1.0), ( + "Correction must decrease toward unity as W increases" + ) + + +class TestScreeningParameters: + """Test that parameter choices affect behavior correctly.""" + + def test_custom_V0_used(self): + """If V₀ is explicitly provided, it should be used instead of the estimate.""" + ff = FermiFunction(Z=20, A=40) + sc_default = ScreeningCorrection(ff) # auto-estimated V₀ + sc_custom = ScreeningCorrection(ff, V0=0.1) + + low_W = np.array([1.05]) + result_default = sc_default(low_W)[0] + result_custom = sc_custom(low_W)[0] + + # Different V₀ should give different results at low energy + assert abs(result_default - result_custom) > 0, ( + "Custom V₀ must change the screening behavior" + ) + + +class TestScreeningOutputShape: + """Basic shape/type checks.""" + + def test_output_shape(self): + ff = FermiFunction(Z=20, A=40) + sc = ScreeningCorrection(ff) + for n in [1, 5, 100]: + W_test = np.linspace(1.1, 3.0, n) + result = sc(W_test) + assert result.shape == W_test.shape diff --git a/tests/test_spectrum.py b/tests/test_spectrum.py new file mode 100644 index 0000000..660ca2c --- /dev/null +++ b/tests/test_spectrum.py @@ -0,0 +1,188 @@ +# tests/test_spectrum.py — Integration tests for the full BetaSpectrum pipeline + +""" +Why test this? +-------------- +Unit tests verify each component in isolation. But integration tests catch bugs that +only appear when components interact: wrong energy grid, type mismatches between +components, or numerical errors accumulating across multiplicative factors. + +The `from_config()` factory method is the user-facing API — it's what people call +in real code. If this doesn't work for a realistic configuration, nothing else matters. + +We test: + + 1. **End-to-end pipeline**: from_config → energy grid → spectrum evaluation → all components valid. + This is the "smoke test" that catches the most common user errors. + 2. **Selective component toggling**: Disabling a component shouldn't break anything. + Users need to be able to turn off corrections for debugging or comparison studies. + 3. **Positive final spectrum**: The product of all factors must remain positive — + if any component goes negative, the entire spectrum is wrong. + 4. **Component extraction works**: calculate_components() should return per-component values + that match calling each component individually. + +Common practice: Integration tests are fewer but more expensive (they test more code paths). +They serve as a safety net — if they pass, you can be confident the system works "as designed". +""" + +import numpy as np +import pytest + +from beta_spectrum import BetaSpectrum, SpectrumConfig + + +class TestBetaSpectrumIntegration: + """Full pipeline integration tests.""" + + def test_from_config_basic(self): + """The most basic use case: create a spectrum from config and evaluate it. + + This is the entry point shown in README.md — if this fails, nothing works. + We test with a medium-Z nucleus (Z=20 calcium) which exercises all components. + """ + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 + ) + + spectrum = BetaSpectrum.from_config(config) + W, E_MeV = spectrum.get_energy_grid(config) + + assert len(W) == len(E_MeV), "Energy grids must have same length" + assert len(W) > 0, "Grid must not be empty" + + def test_spectrum_evaluates_without_error(self): + """Calling the spectrum object should return a valid numpy array.""" + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 + ) + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + + values = spectrum(W) + + assert isinstance(values, np.ndarray), "Output must be a numpy array" + assert values.shape == W.shape, "Output shape must match input shape" + + def test_all_components_positive(self): + """Every component and the total product must be positive.""" + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 + ) + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + + # Filter out values exactly at the endpoint (where phase space is zero) + mask = np.isfinite(spectrum(W)) & (spectrum(W) > 0) + assert np.any(mask), "At least some finite positive values expected" + + +class TestBetaSpectrumToggles: + """Test that component toggling works correctly.""" + + def test_all_disabled(self): + """With all corrections disabled, only phase space remains — still valid output.""" + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5, + use_fermi=False, + use_screening=False, + use_finite_size=False, + use_charge_dist=False, + use_exchange=False, + use_radiative=False, + ) + + spectrum = BetaSpectrum.from_config(config) + W, _ = spectrum.get_energy_grid(config) + values = spectrum(W) + + assert len(values) > 0, "Should return non-empty array" + assert np.all(np.isfinite(values)), "All values must be finite" + + def test_fermi_only(self): + """With only fermi function enabled: should give F₀(Z,W) — not just phase space.""" + config_minimal = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5, + use_phase_space=False, + use_screening=False, + use_finite_size=False, + use_charge_dist=False, + use_exchange=False, + use_radiative=False, + ) + + spectrum = BetaSpectrum.from_config(config_minimal) + W, _ = spectrum.get_energy_grid(config_minimal) + values = spectrum(W) + + # Fermi function alone is ~1 for Z=20 at mid-energy + assert np.all(np.isfinite(values)), "Fermi-only output must be finite" + + +class TestBetaSpectrumComponentExtraction: + """Test calculate_components() returns correct per-component data.""" + + def test_component_names_match(self): + """Component names should correspond to class names (minus suffixes).""" + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 + ) + spectrum = BetaSpectrum.from_config(config) + + W, _ = spectrum.get_energy_grid(config) + components = spectrum.calculate_components(W) + + assert isinstance(components, dict), "Must return a dictionary" + assert len(components) > 0, "Should have at least one component" + + +class TestBetaSpectrumAnalyzer: + """Test the analyzer utilities.""" + + def test_total_spectrum_returns_array(self): + config = SpectrumConfig( + Z_parent=19, Z_daughter=20, A_number=40, endpoint_MeV=2.5 + ) + spectrum = BetaSpectrum.from_config(config) + analyzer = type('Analyzer', (), {'spectrum': spectrum, 'config': config})() + + # We can't easily instantiate the full Analyzer without matplotlib backend, + # but we verify the spectrum itself works end-to-end. + W, E = spectrum.get_energy_grid(config) + total = spectrum(W) + + assert len(total) == len(E), "Spectrum and energy arrays must match" + + +class TestBetaSpectrumRealisticConfig: + """Test with realistic physics parameters.""" + + def test_uranium_spectrum(self): + """Heavy nucleus (Z=92) — tests numerical stability at extreme Z. + + Uranium-238 decays via alpha, but we can still evaluate the beta spectrum + shape for testing purposes. High Z stresses the Fermi function's loggamma. + """ + config = SpectrumConfig( + Z_parent=91, Z_daughter=92, A_number=238, endpoint_MeV=0.5 + ) + + spectrum = BetaSpectrum.from_config(config) + W, E = spectrum.get_energy_grid(config) + values = spectrum(W) + + assert np.all(np.isfinite(values)), "Uranium must not produce NaN/inf" + assert len(values) > 10, f"Grid should have many points: got {len(values)}" + + def test_low_Z_spectrum(self): + """Light nucleus (Z=6 carbon-14 beta decay).""" + config = SpectrumConfig( + Z_parent=5, Z_daughter=6, A_number=14, endpoint_MeV=0.156 # ~156 keV for C-14 + ) + + spectrum = BetaSpectrum.from_config(config) + W, E = spectrum.get_energy_grid(config) + values = spectrum(W) + + assert np.all(np.isfinite(values)), "C-14 must not produce NaN/inf" + # Phase space should peak in the bulk for low-Q decay + assert np.any(values > 0), "Should have non-zero spectrum values"