Use more stable algorithm to estimate decay time.

ChangeLog.rst

@@ -38,6 +38,7 @@ Modified:
 * Updated bound coherent scattering length for H-1, H-2, He-4, C-12,
   O-16, O-17, O-18, Sn-119, Sm-154, Eu-153, Pb-207, Bi-209
 * Updated total cross section for He, Kr, Xe
+* Use log-log interpolation for X-ray f''
 
 Breaking changes:
 

periodictable/activation.py

@@ -165,46 +165,50 @@ def calculate_activation(self, environment, exposure=1,
                         A = activity(el[iso], iso_mass, environment, exposure, rest_times)
                         self._accumulate(A)
 
-    def decay_time(self, target):
+    def decay_time(self, target, tol=1e-10):
         """
         After determining the activation, compute the number of hours required to achieve
         a total activation level after decay.
         """
         if not self.rest_times or not self.activity:
             return 0
 
-        # Find the small rest time (probably 0 hr)
-        min_rest, To = min(enumerate(self.rest_times), key=lambda x: x[1])
+        # Find the smallest rest time (probably 0 hr)
+        k, t_k = min(enumerate(self.rest_times), key=lambda x: x[1])
         # Find the activity at that time, and the decay rate
         data = [
-            (Ia[min_rest], LN2/a.Thalf_hrs)
+            (Ia[k], LN2/a.Thalf_hrs)
             for a, Ia in self.activity.items()
             # TODO: not sure why Ia is zero, but it messes up the initial value guess if it is there
-            if Ia[min_rest] > 0.0
+            if Ia[k] > 0.0
             ]
-        # Build functions for total activity at time T - target and its derivative
-        # This will be zero when activity is at target
-        f = lambda t: sum(Ia*exp(-La*(t-To)) for Ia, La in data) - target
-        df = lambda t: sum(-La*Ia*exp(-La*(t-To)) for Ia, La in data)
-        # Return target time, or 0 if target time is negative
-        if f(0) < target:
-            return 0
         # Need an initial guess near the answer otherwise find_root gets confused.
         # Small but significant activation with an extremely long half-life will
         # dominate at long times, but at short times they will not affect the
         # derivative. Choosing a time that satisfies the longest half-life seems
         # to work well enough.
+        guess = max(-log(target/Ia)/La + t_k for Ia, La in data)
+        # With times far from zero the time resolution in the exponential is
+        # poor. Adjust the start time to the initial guess, rescaling intensities
+        # to the activity at that time.
+        adj = [(Ia*exp(-La*(guess-t_k)), La) for Ia, La in data]
+        #print(adj)
+        # Build f(t) = total activity at time T minus target activity and its
+        # derivative df/dt. f(t) will be zero when activity is at target
+        f = lambda t: sum(Ia*exp(-La*t) for Ia, La in adj) - target
+        df = lambda t: sum(-La*Ia*exp(-La*t) for Ia, La in adj)
         #print("data", data, [])
-        initial = max(-log(target/Ia)/La + To for Ia, La in data)
-        t, ft = find_root(initial, f, df)
+        t, ft = find_root(0, f, df, tol=tol)
         percent_error = 100*abs(ft)/target
         if percent_error > 0.1:
             #return 1e100*365*24 # Return 1e100 rather than raising an error
             msg = (
                 "Failed to compute decay time correctly (%.1g error). Please"
                 " report material, mass, flux and exposure.") % percent_error
             raise RuntimeError(msg)
-        return t
+        # Return time at least zero hours after removal from the beam. Correct
+        # for time adjustment we used to stablize the fit.
+        return max(t+guess, 0.0)
 
     def _accumulate(self, activity):
         for el, activity_el in activity.items():
@@ -287,8 +291,8 @@ def find_root(x, f, df, max=20, tol=1e-10):
     """
     fx = f(x)
     for _ in range(max):
-        #print(f"step {_}: {x=} {fx=} f'={df(x)}")
-        if abs(f(x)) < tol:
+        #print(f"step {_}: {x=} {fx=} df/dx={df(x)} dx={fx/df(x)}")
+        if abs(fx) < tol:
             break
         x -= fx / df(x)
         fx = f(x)
@@ -409,7 +413,7 @@ def activity(isotope, mass, env, exposure, rest_times):
     * AK1495 (Au-198 => Au-199 2n) target should be Au-197
     * AN1428 (Tm-169 => Tm-171 2n) t1/2 updated to Tm-171 rather than Tm-172
     * AN1420 (Er-162 => Ho-163 b) t1/2 updated to 4570 y from 10 y
-    * AT1508 (Pb-208 => Pb-209 act) Thermal (b) x1000 to convert from mbarns to barns
+    * AT1508 (Pb-208 => Pb-209 act) Thermal (b) x 1000 to convert from mbarns to barns
     """
     # TODO: is the table missing 1-H => 3-H ?
     # 0nly activations which produce radioactive daughter products are
@@ -650,8 +654,10 @@ def init(table, reload=False):
 class ActivationResult(object):
     def __init__(self, **kw):
         self.__dict__ = kw
+    def __repr__(self):
+        return f"ActivationResult({self.isotope},{self.reaction},{self.daughter})"
     def __str__(self):
-        return str(self.__dict__)
+        return f"{self.isotope}={self.reaction}=>{self.daughter}"
 
 
 def demo():  # pragma: nocover

test/test_activation.py

@@ -31,30 +31,17 @@ def _get_Au_activity(fluence=1e5):
     assert (act2 - 1000*act1) < 1e-8
 
     # Smoke test: does every element run to completion?
-    formula = "".join(str(el) for el in pt.elements)[1:]
+    formula = "".join(str(el) for el in pt.elements)
     # Use an enormous mass to force significant activation of rare isotopes
     mass, fluence, exposure = 1e15, 1e8, 10
     env = ActivationEnvironment(fluence=fluence, Cd_ratio=70, fast_ratio=50, location="BT-2")
     sample = Sample(formula, mass=mass)
     sample.calculate_activation(
-        env, exposure=exposure, rest_times=(0, 1, 24, 360),
+        env, exposure=exposure, rest_times=(0, 1, 24),
         abundance=IAEA1987_isotopic_abundance,
         #abundance=table_abundance,
         )
-
-    # TODO: Why aren't we matching the spreadsheet results?
-    results = {
-        # The following were labelled as spreadsheet results
-        #"Co-60m+":  [5186.888,  98.878,   2.751e-38],
-        #"Co-61":    [2.689e-8,  1.767e-8, 1.127e-12],
-        #"Co-60":    [1.550322,  1.550299, 1.549764],
-        #"Co-61":    [5.695e-8,  3.741e-8, 2.386e-12],
-        # Values as of R1.7.0
-        "Co-60m+": [5088.093018977685, 96.91335724994167, 2.6450954672880745e-38],
-        "Co-61": [7.309565456941425e-09, 4.802298351964856e-09, 3.0568451604601675e-13],
-        "Co-60": [0.15507562804846606, 0.15507330056693813, 0.15501977813208836],
-        "Co-61": [1.3649499351603886e-09, 8.967560195949139e-10, 5.708192406435261e-14],
-    }
+    #sample.show_table(cutoff=0)
 
     sample = Sample('Co', mass=10)
     env = ActivationEnvironment(fluence=1e8)
@@ -63,25 +50,66 @@ def _get_Au_activity(fluence=1e5):
         abundance=IAEA1987_isotopic_abundance,
         #abundance=table_abundance,
         )
-    print(sample.activity)
+    #sample.show_table(cutoff=0)
+    # ACT_2N_X.xls for 10g Co at Io=1e8 for t = [0, 1, 24]
+    # There are duplicate entries for Co-61 because there are different
+    # activation paths for 59Co + n -> 61Co.
+    # Results are limited to 3 digits because formulas use ln(2) = 0.693. There
+    # is also precision loss on the calculation of differences in exponentials,
+    # with exp(a) - exp(b) converted to exp(b)(exp(a-b) - 1) = exp(b)expm1(a-b)
+    # in activation.py.
+    # results = {
+    #     "C0-60m+":  [5.08800989419543E+03,	9.69933141298983E+01,	2.69898447909379E-38],
+    #     "Co-61":    [7.30937687202485E-09,	4.80260282859366E-09,	3.06331725449002E-13],
+    #     "Co-60":    [1.55042700053951E-01,	1.55040373560730E-01,	1.54986873850802E-01],
+    #     "Co-61":    [1.36469792582999E-09,	8.96670432174800E-10,	5.71936948464344E-14],
+    # }
+    # Results from R2.0.0-pre
+    results = {
+        "Co-60m+":	[5.08746786932552e+03,	9.69014499692868e+01,	2.64477047705988e-38],
+        "Co-61":	[7.30866736559643e-09,	4.80170831653643e-09,	3.05646958051663e-13],
+        "Co-60":	[1.55056574647797e-01,	1.55054247452235e-01,	1.55000731593431e-01],
+        "Co-61":	[1.36478223029061e-09,	8.96645839472098e-10,	5.70749106813738e-14],
+    }
+    #print(list(sample.activity.keys()))
+    #print(dir(list(sample.activity.keys())[0]))
     for product, activity in sample.activity.items():
-        print(f'        "{product.daughter}": {activity},')
+        #print(product)
+        #print(dir(product))
+        # Uncomment to show new table values
+        activity_str = ",\t".join(f"{Ia:.14E}" for Ia in activity)
+        #print(f'        "{product.daughter}":\t[{activity_str}],')
         assert product.daughter in results, f"Missing {product.daughter}"
-        print(results[product.daughter])
-        print(activity)
-        #assert np.allclose(results[product.daughter], activity, atol=0, rtol=1e-12)
-
+        # TODO: include duplicate decay paths in test, or identical daughters
+        # Test that results haven't changed since last update
+        if product.daughter != "Co-61":
+            assert np.allclose(results[product.daughter], activity, atol=0, rtol=1e-12)
 
     # 129-I has a long half-life (16 My) so any combination of exposure
     # fluence and mass that causes significant activation will yield a product
     # that is radioactive for a very long time.
     sample = Sample('Te', mass=1e13)
     env = ActivationEnvironment(fluence=1e8)
-    sample.calculate_activation(env, rest_times=[0])
+    sample.calculate_activation(env, rest_times=[1,10,100])
+    #sample.show_table(cutoff=0)
     target = 1e-5
-    decay = sample.decay_time(target)
-    sample.calculate_activation(env, rest_times=[0, decay])
+    t_decay = sample.decay_time(target)
+    #print(f"{t_decay=}")
+    sample.calculate_activation(env, rest_times=[t_decay])
     total = sum(v[-1] for v in sample.activity.values())
-    assert abs(total - target)/target < 1e-6, f"total activity {total} != {target}"
+    assert abs(total - target) < 1e-10, f"total activity {total} != {target}"
+
+    # Al and Si daughters have short half-lives
+    sample = Sample('AlSi', mass=1e13)
+    env = ActivationEnvironment(fluence=1e8)
+    sample.calculate_activation(env, rest_times=[100,200])
+    #sample.show_table(cutoff=0)
+    target = 1e-5
+    t_decay = sample.decay_time(target)
+    #print(f"{t_decay=}")
+    sample.calculate_activation(env, rest_times=[t_decay])
+    total = sum(v[-1] for v in sample.activity.values())
+    assert abs(total - target) < 1e-10, f"total activity {total} != {target}"
+
 
 test()