From 2bc1ffc2d48ef75f2aced1a7dc197300f568010a Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Thu, 24 Oct 2024 17:30:11 -0400 Subject: [PATCH] Use more stable algorithm to estimate decay time. --- ChangeLog.rst | 1 + periodictable/activation.py | 44 +++++++++++--------- test/test_activation.py | 80 +++++++++++++++++++++++++------------ 3 files changed, 80 insertions(+), 45 deletions(-) diff --git a/ChangeLog.rst b/ChangeLog.rst index 1473e3d..57d34e3 100644 --- a/ChangeLog.rst +++ b/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: diff --git a/periodictable/activation.py b/periodictable/activation.py index a5a48b4..1ba1acd 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -165,7 +165,7 @@ 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. @@ -173,30 +173,32 @@ def decay_time(self, target): 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 @@ -204,7 +206,9 @@ def decay_time(self, target): "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 diff --git a/test/test_activation.py b/test/test_activation.py index bf0816a..77b8a6e 100644 --- a/test/test_activation.py +++ b/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()