Stabilize activation at coincident capture rates

periodictable/activation.py

@@ -227,6 +227,7 @@
 
 
 from math import exp, log, expm1
+from sys import float_info
 import os
 from collections.abc import Callable, Sequence
 from typing import cast
@@ -808,14 +809,17 @@ def activity(
             # Column S: resulting product [1/h]
             product_2n = lam if ai.reaction == '2n' else 0
             # Column T: activity if 2n mode
-            activity = root*lam*(parent_activity-parent_lam)*(
-                (exp(-lam_2n*exposure)
-                 / ((parent_activity-lam_2n)*(product_2n-lam_2n)))
-                + (exp(-parent_activity*exposure)
-                   / ((lam_2n-parent_activity)*(product_2n-parent_activity)))
-                + (exp(-product_2n*exposure)
-                   / ((lam_2n-product_2n)*(parent_activity-product_2n)))
-                )
+            #     activity = root*lam*(parent_activity-parent_lam)*(
+            #         (exp(-lam_2n*exposure)
+            #          / ((parent_activity-lam_2n)*(product_2n-lam_2n)))
+            #         + (exp(-parent_activity*exposure)
+            #            / ((lam_2n-parent_activity)*(product_2n-parent_activity)))
+            #         + (exp(-product_2n*exposure)
+            #            / ((lam_2n-product_2n)*(parent_activity-product_2n)))
+            #         )
+            # [CAOShurong 2026-10-15] Evaluate T as a three-rate convolution.
+            activity = root*lam*(parent_activity-parent_lam)*_three_rate_integral(
+                exposure, lam_2n, parent_activity, product_2n)
             #print("N", parent_lam, "P", effectiveXS, "Q", lam_2n,
             #      "R", parent_activity, "S", product_2n, "T", activity)
         else:
@@ -830,15 +834,24 @@ def activity(
             # Note: This cross-section always uses the total thermal flux
             effectiveXS = ai.thermalXS_parent + env.epithermal_reduction_factor*ai.resonance_parent
             # Column U: nv1s1t [neutrons]
-            U = (flux*3600)*(initialXS*1e-24)*exposure
+            #     U = (flux*3600)*(initialXS*1e-24)*exposure
             # Column V: nv2s2t+L2*t [neutrons]
-            V = ((env.fluence*3600)*(effectiveXS*1e-24)+lam)*exposure
+            #     V = ((env.fluence*3600)*(effectiveXS*1e-24)+lam)*exposure
             # Column W: L/(L-nvs1+nvs2) [unitless]
-            W = lam/(lam-(flux*3600)*(initialXS*1e-24)+(env.fluence*3600)*(effectiveXS*1e-24))
-            # Column X: W*(exp(-U)-exp(V)) if U,V > 1e-10 else W*(V-U+(V-U)*(V+U)/2)
-            # [PAK 2024-02-28] Rewrite the exponential difference using expm1()
-            #Xp = exp(-V)*expm1(V-U) if U>V else -exp(-U)*expm1(U-V)
-            X = W*exp(-V)*expm1(V-U) if U>V else -W*exp(-U)*expm1(U-V)
+            #     W = lam/(lam-(flux*3600)*(initialXS*1e-24)+(env.fluence*3600)*(effectiveXS*1e-24))
+            # Column X: W*(exp(-U)-exp(-V))
+            #     Original small-U,V form: W*(V-U+(V-U)*(V+U)/2)
+            # [PAK 2024-02-28] Rewrite the exponential difference using expm1():
+            #     X = W*exp(-V)*expm1(V-U) if U>V else -W*exp(-U)*expm1(U-V)
+            # [CAOShurong 2026-10-15] Rewrite U, V, W, X for higher accuracy.
+            # With a=production_rate, b=removal_rate, t=exposure:
+            #     U = a*t; V = b*t; W = lam/(b-a)
+            #     X = lam*(exp(-a*t)-exp(-b*t))/(b-a)
+            # The convolution also provides the finite equal-rate limit.
+            production_rate = (flux*3600)*(initialXS*1e-24)
+            removal_rate = (env.fluence*3600)*(effectiveXS*1e-24) + lam
+            # The exponential divided difference has a finite equal-rate limit.
+            X = lam*_two_rate_integral(exposure, production_rate, removal_rate)
             # Column Y: O if "b" else T if "2n" else L*X
             activity = root*X
             #print(f"{ai.isotope}=>{ai.daughter} {U=} {V=} {W=} {Xp=} {X=} {activity=} {lam=} {flux*initialXS*3600*1e-24+env.fluence*effectiveXS*3600*1e-24}")
@@ -877,6 +890,35 @@ def activity(
 
     return result
 
+def _two_rate_integral(t: float, a: float, b: float) -> float:
+    """Convolve exp(-a*t) and exp(-b*t), including coincident rates."""
+    if a > b:
+        a, b = b, a
+    x = (b-a)*t
+    # A zero or subnormal product loses precision before division by b-a.
+    if abs(x) < float_info.min:
+        return t*exp(-a*t)
+    # Divide first so a subnormal decay is not lost in the intermediate product.
+    return (-expm1(-x)/(b-a))*exp(-a*t)
+
+def _three_rate_integral(t: float, a: float, b: float, c: float) -> float:
+    """Convolve three exponentials without singular repeated-rate factors."""
+    a, b, c = sorted((a, b, c))
+    x, y = (b-a)*t, (c-a)*t
+    if y > 0.01:
+        # Use the largest separation in the denominator.
+        return (_two_rate_integral(t, a, b)-_two_rate_integral(t, b, c))/(c-a)
+
+    # Near the all-equal limit, expand the second divided difference.
+    # h_k(x,y) = sum(x**j*y**(k-j)); coefficient truncation error is < 3e-22.
+    coefficient, homogeneous, x_power, factorial = 0.5, 1., 1., 2
+    for order in range(1, 8):
+        x_power *= x
+        homogeneous = y*homogeneous + x_power
+        factorial *= order+2
+        coefficient += (-1)**order * homogeneous/factorial
+    return t*(t*exp(-a*t))*coefficient
+
 class ActivationResult:
     r"""
     *isotope* :

test/test_activation.py

@@ -69,7 +69,8 @@ def _get_Au_activity(fluence=1e5):
         ("Co-60m+", [5.08870610591228E+03,	9.57100873523175E+01,	1.95440969864725E-38]),
         ("Co-61", [7.29613649387052E-09,	4.79211894136011E-09,	3.03056972241711E-13]),
         ("Co-60", [1.54971558843742E-01,	1.54969234199432E-01,	1.54915777003567E-01]),
-        ("Co-61", [1.36558180150732E-09,	8.96917213990930E-10,	5.67216754320407E-14]),
+        # Stable two-capture result, checked against a 100-digit calculation.
+        ("Co-61", [1.3655818015418685E-09, 8.969172140136206E-10, 5.672167543347562E-14]),
     ]
     #print(list(sample.activity.keys()))
     #print(" ".join(k for k in dir(list(sample.activity.keys())[0]) if k[0] != '_'))

test/test_activation_capture.py

@@ -0,0 +1,109 @@
+from math import exp, log, nextafter, inf
+from decimal import Decimal, localcontext
+from itertools import permutations
+
+import numpy as np
+import pytest
+
+import periodictable as pt
+from periodictable.activation import ActivationEnvironment, Sample
+
+
+@pytest.mark.parametrize('offset', [-1, 0, 1])
+@pytest.mark.parametrize('exposure', [0.0, 1.0, 10000.0])
+def test_coincident_capture_rates(offset, exposure):
+    isotope = pt.Co[59]
+    product = next(row for row in isotope.neutron_activation
+                   if row.reaction == 'act' and row.daughter == 'Co-60')
+    decay_rate = log(2) / product.Thalf_hrs
+    flux = decay_rate / (3600e-24 * (product.thermalXS - product.thermalXS_parent))
+    if offset:
+        flux = nextafter(flux, 0.0 if offset < 0 else inf)
+
+    sample = Sample('Co', mass=1)
+    sample.calculate_activation(ActivationEnvironment(fluence=flux),
+                                exposure=exposure, rest_times=[0.0, 10.0])
+    activity = next(values for row, values in sample.activity.items()
+                    if row.reaction == 'act' and row.daughter == 'Co-60')
+
+    # Coincident depletion/removal rates have a finite two-state limit.
+    capture_rate = flux * product.thermalXS * 3600e-24
+    production = flux * product.thermalXS * 1e-24 / isotope.isotope * 1.6276e19
+    initial = production * decay_rate * exposure * exp(-capture_rate * exposure)
+    expected = [initial, initial * exp(-decay_rate * 10)]
+    np.testing.assert_allclose(activity, expected, rtol=1e-12, atol=0)
+
+
+@pytest.mark.parametrize('time', [0.0, 1.0, 1000.0])
+def test_repeated_integrals(time):
+    from periodictable.activation import _two_rate_integral, _three_rate_integral
+
+    rate = 0.2
+    np.testing.assert_allclose(_two_rate_integral(time, rate, rate),
+                               time * exp(-rate*time), rtol=1e-14, atol=0)
+    np.testing.assert_allclose(_three_rate_integral(time, rate, rate, rate),
+                               time**2 * exp(-rate*time) / 2, rtol=1e-14, atol=0)
+
+
+@pytest.mark.parametrize('time', [0., 1e-8, 1., 100.])
+@pytest.mark.parametrize('rates', [
+    (0., 2.),
+    (0.2, 1.2),
+    (1., nextafter(1., inf)),
+    (1e-6, 1e-6 + 1e-12),
+    (20., 20.),
+])
+def test_two_rate_integral(time, rates):
+    from periodictable.activation import _two_rate_integral
+
+    with localcontext() as context:
+        context.prec = 100
+        t = Decimal.from_float(time)
+        a, b = map(Decimal.from_float, rates)
+        expected = float(t*(-a*t).exp() if a == b else
+                         ((-a*t).exp() - (-b*t).exp()) / (b-a))
+    for ordering in permutations(rates):
+        np.testing.assert_allclose(_two_rate_integral(time, *ordering),
+                                   expected, rtol=2e-14, atol=0)
+
+
+@pytest.mark.parametrize('time, rate', [(1e-200, 1e-125), (1e-124, 3e-200)])
+def test_two_rate_integral_subnormal_product(time, rate):
+    from periodictable.activation import _two_rate_integral
+
+    # With one zero rate and a subnormal rate*time, the integral rounds to time.
+    for rates in [(0., rate), (rate, 0.)]:
+        assert _two_rate_integral(time, *rates) == time
+
+
+def test_two_rate_integral_subnormal_decay():
+    from periodictable.activation import _two_rate_integral
+
+    time = 740.
+    rates = (1., nextafter(1., inf))
+    # Preserve the nonzero repeated-rate limit when exp(-rate*time) is subnormal.
+    expected = time*exp(-time)
+    for ordering in permutations(rates):
+        assert _two_rate_integral(time, *ordering) == expected
+
+
+@pytest.mark.parametrize('rates', [
+    (1., 1.+1e-12, 1.+2e-12),
+    (0., 0.004, 0.009),
+    (0., 0.004, 0.010),
+    (0., 0.004, 0.011),
+    (0.2, 1.2, 3.),
+])
+def test_three_distinct_rates(rates):
+    from periodictable.activation import _three_rate_integral
+
+    # Evaluate the original divided difference at independent high precision.
+    with localcontext() as context:
+        context.prec = 100
+        a, b, c = map(Decimal.from_float, rates)
+        expected = float((-a).exp()/((b-a)*(c-a))
+                         + (-b).exp()/((a-b)*(c-b))
+                         + (-c).exp()/((a-c)*(b-c)))
+    for ordering in permutations(rates):
+        np.testing.assert_allclose(_three_rate_integral(1., *ordering),
+                                   expected, rtol=1e-12, atol=0)