From 05290e8e3a110c699f13b2d258038e54ca9a0570 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Fri, 5 Jun 2026 12:38:59 -0400 Subject: [PATCH 1/4] Update 150Nd to use the 3-stage Bateman equation for 151Sm levels --- periodictable/activation.py | 105 ++++++++++++++++++++++++++++++++---- 1 file changed, 96 insertions(+), 9 deletions(-) diff --git a/periodictable/activation.py b/periodictable/activation.py index 0a82583..b334839 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -336,6 +336,11 @@ def decay_time(self, target: float, tol: float=1e-10): """ After determining the activation, compute the number of hours required to achieve a total activation level after decay. + + To keep the code simple, the 150Nd activation is treated as 151Pm -> 151Sm + instead of 151Nd -> 151Pm -> 151Sm. The missing 151Nd activity is added + to 151Sm to compensate, leading to a slightly long decay time estimate + for short export times. """ if not self.rest_times or not self.activity: return 0 @@ -357,6 +362,7 @@ def decay_time(self, target: float, tol: float=1e-10): MIN_TIME = -1e100 # Allow the initial guess to be negative t_guess = MIN_TIME initial_activity = 0. # Accumulate the intensity at t=0. + activity_Nd151 = 0. # Linter doesn't know activity_Nd151 is defined before use for k, (a, Ia) in enumerate(self.activity.items()): intensity = Ia[0] @@ -372,6 +378,19 @@ def decay_time(self, target: float, tol: float=1e-10): Ip, Lp, _ = data[-1] # parent activity and decay constant #print(f"{intensity=} {Ip=} {Lp=} {Ia[0]=} {La=}") intensity -= Ip*(expm1((La-Lp)*t)/(1 - Lp/La)) if Ip > 0. else 0. + # Hack: Subtract scaled 151Nd activity from 151Sm so that we don't need + # to fit with the three stage Bateman equation. This works because there + # is a million-fold difference in half-lives. If there is enough activity + # in 151Nd to push 151Sm over the threshold then the 151Sm intensity will + # determine decay time, otherwise 151Pm will dominate. Because 151Sm activity + # is front-loaded, the integrated intensity for short times will be a little + # higher than expected, giving a slightly long decay time estimate when + # 151Sm activity is below the threshold. + if a.daughter.startswith("Nd-151"): # Nd-151t + activity_Nd151 = intensity + elif a.isotope == "Nd-150" and a.daughter.startswith("Sm-151"): # Sm-151 + # print(f"151Nd correction to 151Sm = {activity_Nd151 * a.Thalf_parent / a.Thalf_hrs}") + intensity += activity_Nd151 * a.Thalf_parent / a.Thalf_hrs initial_activity += intensity data.append((intensity, La, a.reaction)) @@ -674,7 +693,7 @@ def activity( rest period. Activations are listed in *isotope.neutron_activation*. Most of the - activations (n,g n,p n,a n,2n) define a single step process, where a neutron + activations (n,g n,p n,a) define a single step process, where a neutron is absorbed yielding the daughter and some prompt radiation. The daughter itself will decay during exposure, yielding a balance between production and emission. Any exposure beyound about eight halflives will not increase @@ -682,7 +701,16 @@ def activity( Activity for daughter products may undergo further neutron capture, reducing the activity of the daughter product but introducing a grand daughter with - its own activity. + its own activity (2n reactions). Active delayed beta products (b mode reactions) + accumulate during irradiation, eventually reaching equilibrium. Burn up, where + the activated product undergoes further capture before beta decay is not + accounted for. The delayed beta product increases in activity until the + directly activated product has decayed away following the 2-stage Bateman + equation. 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu has two transient + daughters. During irradiation, the activity for 151Sm is computed + directly from 151Nd, skipping 151Pm. After irradiation the 151Pm activity + is subtracted from 151Sm, and decay calculated using the 3-stage Bateman + equation. The data tables for activation are only precise to about three significant figures. Any changes to the calculations below this threshold, e.g., due to @@ -726,8 +754,8 @@ def activity( # Hack: delayed beta calculation needs activity and decay rate # of the parent. Since the "b" mode reaction line follows directly # after the parent line, save the last activity and the last decay rate - # at each step of the loop. - last_activity = last_lam = 0. + # at each step of the loop. 151Sm also needs the 151Nd activity + last_activity = last_lam = activity_151Nd = 0. for ai in isotope.neutron_activation: # Ignore fast neutron interactions if not using fast ratio if ai.fast and env.fast_ratio == 0: @@ -859,13 +887,38 @@ def activity( # - It does not include activity from 151Nd in 151Sm, but that should be insignificant. # TODO: check 151Eu -> 152m2Eu -> 152Eu is missing b-mode 152Gd # - It is present for 151Eu -> 152m1Eu isomer and 151Eu -> 152Eu ground state. - if ai.reaction == 'b': + + if ai.daughter.startswith("Nd-151"): # Nd-151t + # Track 151Nd activity at t=0 so we can correct 151Sm for short exposure + activity_151Nd = activity + + if ai.isotope == "Nd-150" and ai.daughter.startswith("Sm-151"): # Sm-151 + # Decay chain: 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu + # Because 151Sm uses 151Nd as its parent halflife, all the intensity from the + # transient 151Pm is accounted for in the 151Sm activity. We subtract the scaled 151Pm + # activity at t=0 from 151Sm. We use the three stage bateman equation to dribble activity + # from 151Nd and 151Pm into 151Sm. + # print(f"151Sm @ t=0 = {activity}") + # print(f" - (activity_151Pm={last_activity})*{lam}/{last_lam} = {last_activity * lam/last_lam}") + # print(f" = {activity - last_activity * lam/last_lam}") + activity -= last_activity * lam/last_lam + + # Show the approximate activity in 151Sm after 151Nd and 151Pm have decayed + # print(f"151Sm @ t=280 h = {activity}") + # print(f" + ({activity_151Nd=})*{lam}/{parent_lam} = {activity_151Nd * lam/parent_lam}") + # print(f" + (activity_151Pm={last_activity})*{lam}/{last_lam} = {last_activity * lam/last_lam}") + # print(f" = {activity + activity_151Nd * lam/parent_lam + last_activity * lam/last_lam}") + + # Note: table for 151Sm has 151Nd for parent half-life. + result[ai] = [ + bateman3(Ti, activity, lam, last_activity, last_lam, activity_151Nd, parent_lam) + for Ti in rest_times + ] + + elif ai.reaction == 'b': # Accumulate build up following the Bateman equation: - # A_d(t) = λ_d/(λ_d - λ_p) A_p(0) (exp(-λ_p t) - exp(-λ_d t)) - # Add this to decay of the granddaughter at the end of the exposure: - # + A_d(0) exp(-λ_d t) result[ai] = [ - activity*exp(-lam*Ti) + last_activity * (exp(-last_lam*Ti) - exp(-lam*Ti)) / (1 - last_lam/lam) + bateman2(Ti, activity, lam, last_activity, last_lam) for Ti in rest_times ] # print(f"{ai.daughter} {activity=} {last_activity=}") @@ -877,6 +930,40 @@ def activity( return result +def bateman2(t, act, lam, act_p, lam_p): + """ + Accumulate build up following the Bateman equation:: + + A1(t) = A(0) exp(-λ t) + A2(t) = λ/(λ - λ_p) A_p(0) (exp(-λ_p t) - exp(-λ t)) + + A(t) = A1(t) + A2(t) + + """ + # simple decay + term1 = act*exp(-lam*t) + # two-stage decay + term2 = lam / (lam - lam_p) * act_p * (exp(-lam_p*t) - exp(-lam*t)) + + return term1 + term2 + +def bateman3(t, act, lam, act_p, lam_p, act_gp, lam_gp): + """ + Accumulate build up following the three stage Bateman equation. + """ + # simple decay + term1 = act * exp(-lam*t) + # two-stage decay (written in same form as three stage) + term2 = act_p * lam * (exp(-lam*t)/(lam_p - lam) + exp(-lam_p*t)/(lam - lam_p)) + # three-stage decay + term3 = act_gp * lam * lam_p * ( + exp(-lam*t)/((lam_p - lam)*(lam_gp - lam)) + + exp(-lam_p*t)/((lam - lam_p)*(lam_gp - lam_p)) + + exp(-lam_gp*t)/((lam - lam_gp)*(lam_p - lam_gp)) + ) + + return term1 + term2 + term3 + class ActivationResult: r""" *isotope* : From 2e713e70c16e6cafdf2f68429d995d6077159903 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Fri, 5 Jun 2026 15:36:36 -0400 Subject: [PATCH 2/4] tweak docs and comments re: 150Nd activation --- periodictable/activation.py | 49 +++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/periodictable/activation.py b/periodictable/activation.py index b334839..e3e4989 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -337,10 +337,11 @@ def decay_time(self, target: float, tol: float=1e-10): After determining the activation, compute the number of hours required to achieve a total activation level after decay. - To keep the code simple, the 150Nd activation is treated as 151Pm -> 151Sm - instead of 151Nd -> 151Pm -> 151Sm. The missing 151Nd activity is added - to 151Sm to compensate, leading to a slightly long decay time estimate - for short export times. + To keep the code simple, the 151Sm decay from 150Nd activation is + treated as 151Pm -> 151Sm instead of 151Nd -> 151Pm -> 151Sm. The + missing 151Nd activity is added to 151Sm at t=0 to compensate, leading + to a slightly long decay time estimate for short exposure times, and + a slightly short decay time estimate for long exposure times. """ if not self.rest_times or not self.activity: return 0 @@ -451,7 +452,7 @@ def dfdt(t): #return 0. # Decay time failed to compute; silently fail #return 1e100*365*24 # Return 1e100 rather than raising an error msg = ( - f"Failed to compute decay time correctly ({percent_error:.1g} error)." + f"Failed to compute decay time correctly ({percent_error:.1g}% error)." f" Please report material and activation parameters.") raise RuntimeError(msg) @@ -693,24 +694,24 @@ def activity( rest period. Activations are listed in *isotope.neutron_activation*. Most of the - activations (n,g n,p n,a) define a single step process, where a neutron - is absorbed yielding the daughter and some prompt radiation. The daughter + activations (n,g n,p n,a) define a single step process, where a neutron is + absorbed yielding the daughter and some prompt radiation. The daughter itself will decay during exposure, yielding a balance between production and emission. Any exposure beyound about eight halflives will not increase - activity for that product. - - Activity for daughter products may undergo further neutron capture, reducing - the activity of the daughter product but introducing a grand daughter with - its own activity (2n reactions). Active delayed beta products (b mode reactions) - accumulate during irradiation, eventually reaching equilibrium. Burn up, where - the activated product undergoes further capture before beta decay is not - accounted for. The delayed beta product increases in activity until the - directly activated product has decayed away following the 2-stage Bateman - equation. 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu has two transient - daughters. During irradiation, the activity for 151Sm is computed - directly from 151Nd, skipping 151Pm. After irradiation the 151Pm activity - is subtracted from 151Sm, and decay calculated using the 3-stage Bateman - equation. + activity for that product. Activity for daughter products may undergo + further neutron capture, reducing the activity of the daughter product but + introducing a grand daughter with its own activity (n,2n). + + Active delayed beta products (b mode reactions) accumulate during + irradiation, eventually reaching equilibrium. Burn up, where the activated + product undergoes further capture before beta decay is not accounted for. + The delayed beta product increases in activity until the directly activated + product has decayed away following the 2-stage Bateman equation. + + 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu has two transient daughters. + During irradiation, the activity for 151Sm is computed directly from 151Nd, + skipping 151Pm. After irradiation the 151Pm activity is subtracted from + 151Sm, and decay calculated using the 3-stage Bateman equation. The data tables for activation are only precise to about three significant figures. Any changes to the calculations below this threshold, e.g., due to @@ -819,7 +820,6 @@ def activity( # = root * (lam*expm1(x2) - parent_lam*expm1(x1)) / (parent_lam - lam) # Checked for each b-mode production that small halflife results are # unchanged to four digits and Eu[151] => Gd[152] no longer fails. - # TODO: 150Nd => 151Nd => 151Pm => 151Sm => 151Eu activity = root/(parent_lam - lam) * ( lam*expm1(-parent_lam*exposure) - parent_lam*expm1(-lam*exposure)) #print("N", parent_lam, "O", activity) @@ -883,8 +883,6 @@ def activity( #print(" ".join("%.5e"%v for v in data)) # 2026-05-27 PAK: delayed beta decay such as 209Bi -> 210Bi -> 210Po -> 206Pb - # TODO: check 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu - # - It does not include activity from 151Nd in 151Sm, but that should be insignificant. # TODO: check 151Eu -> 152m2Eu -> 152Eu is missing b-mode 152Gd # - It is present for 151Eu -> 152m1Eu isomer and 151Eu -> 152Eu ground state. @@ -892,7 +890,7 @@ def activity( # Track 151Nd activity at t=0 so we can correct 151Sm for short exposure activity_151Nd = activity - if ai.isotope == "Nd-150" and ai.daughter.startswith("Sm-151"): # Sm-151 + elif ai.isotope == "Nd-150" and ai.daughter.startswith("Sm-151"): # Sm-151 # Decay chain: 150Nd -> 151Nd -> 151Pm -> 151Sm -> 151Eu # Because 151Sm uses 151Nd as its parent halflife, all the intensity from the # transient 151Pm is accounted for in the 151Sm activity. We subtract the scaled 151Pm @@ -909,7 +907,6 @@ def activity( # print(f" + (activity_151Pm={last_activity})*{lam}/{last_lam} = {last_activity * lam/last_lam}") # print(f" = {activity + activity_151Nd * lam/parent_lam + last_activity * lam/last_lam}") - # Note: table for 151Sm has 151Nd for parent half-life. result[ai] = [ bateman3(Ti, activity, lam, last_activity, last_lam, activity_151Nd, parent_lam) for Ti in rest_times From 1cb0317b3232476ae7bd312ef2efa1f0bf16661c Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Fri, 5 Jun 2026 18:50:44 -0400 Subject: [PATCH 3/4] clarify that only the decay time estimate uses the simplication --- periodictable/activation.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/periodictable/activation.py b/periodictable/activation.py index e3e4989..f0c478e 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -337,11 +337,13 @@ def decay_time(self, target: float, tol: float=1e-10): After determining the activation, compute the number of hours required to achieve a total activation level after decay. - To keep the code simple, the 151Sm decay from 150Nd activation is - treated as 151Pm -> 151Sm instead of 151Nd -> 151Pm -> 151Sm. The - missing 151Nd activity is added to 151Sm at t=0 to compensate, leading - to a slightly long decay time estimate for short exposure times, and - a slightly short decay time estimate for long exposure times. + To simplify the decay time estimate, the 151Sm decay from 150Nd activation is + treated as 151Pm -> 151Sm instead of 151Nd -> 151Pm -> 151Sm, with the + missing 151Nd activity added to 151Sm at t=0 to compensate. This is + different from the decay calculation, which uses the full 3-stage Bateman + equation for the 151Sm activity. As a result, the decay time estimate will + be slightly long for short exposures and slightly short for long exposures, + depending on whether 151Sm activity exceeds the target threshhold. """ if not self.rest_times or not self.activity: return 0 From 0f6a2f39feed17c6ae3f4335b808ce6ef1fef6ab Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Fri, 5 Jun 2026 18:52:29 -0400 Subject: [PATCH 4/4] clarify that only the decay time estimate uses the simplication --- periodictable/activation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/periodictable/activation.py b/periodictable/activation.py index f0c478e..eb3b0f9 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -343,7 +343,7 @@ def decay_time(self, target: float, tol: float=1e-10): different from the decay calculation, which uses the full 3-stage Bateman equation for the 151Sm activity. As a result, the decay time estimate will be slightly long for short exposures and slightly short for long exposures, - depending on whether 151Sm activity exceeds the target threshhold. + depending on whether 151Sm activity exceeds the target threshold. """ if not self.rest_times or not self.activity: return 0