diff --git a/periodictable/activation.py b/periodictable/activation.py index 0a82583..eb3b0f9 100644 --- a/periodictable/activation.py +++ b/periodictable/activation.py @@ -336,6 +336,14 @@ 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 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 threshold. """ if not self.rest_times or not self.activity: return 0 @@ -357,6 +365,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 +381,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)) @@ -432,7 +454,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) @@ -674,15 +696,24 @@ 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 - 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 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. - 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. + 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 +757,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: @@ -791,7 +822,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) @@ -855,17 +885,39 @@ 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. - 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 + + 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 + # 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}") + + 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 +929,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* :