Skip to content

Debug failure of decay time estimate for trace N2 #129

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pkienzle merged 3 commits into from
Jun 3, 2026
Merged

Debug failure of decay time estimate for trace N2 #129

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b6e2a82
Debug failure of decay time estimate for trace N2
Jun 2, 2026
f7ea74b
raise an error in the calculator; we can catch it in the web callback
Jun 2, 2026
3ae56bf
fix potential bug in low activity decay time estimate
Jun 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
60 changes: 40 additions & 20 deletions periodictable/activation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -354,63 +354,82 @@ def decay_time(self, target: float, tol: float=1e-10):
# dominate at long times, but at short times they will not affect the
# derivative. Choosing a time that satisfies the longest decay time seems
# to work well enough.
t_guess = 0.
MIN_TIME = -1e100 # Allow the initial guess to be negative
t_guess = MIN_TIME
initial_activity = 0. # Accumulate the intensity at t=0.

for a, Ia in self.activity.items():
for k, (a, Ia) in enumerate(self.activity.items()):
intensity = Ia[0]
La = LN2/a.Thalf_hrs
# print(f"{k}: {a.daughter} {intensity=} {La=}")
# Estimate activity at t=0, with a check that it hasn't decayed to zero.
# The check is necessary because exp(λt) will exceed the floating point
# number range when I(t) = I(0)exp(-λt) = 0.
if t > 0.:
if t > 0. and intensity > 0.:
# print(f"{a.isotope} {a.daughter} {a.Thalf_hrs=} {La=}; {intensity=} at {t=}")
intensity = intensity * exp(La*t) if intensity > 0. else 0.
intensity = intensity * exp(La*t)
if a.reaction == "b":
Ip, Lp, _ = data[-1]
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.
initial_activity += intensity
data.append((intensity, La, a.reaction))

# Daughter intensity will be transformed to granddaughter intensity
# over time. Assume it is instantaneous at t=0 for the purpose of guessing
# the decay time of the granddaughter. The root finder will correct for
# the discrepency.
I_buildup = data[-1][0]*La/data[-1][1] if a.reaction == "b" else 0.
decay_estimate = -log(target/(intensity + I_buildup))/La if intensity > 0. else 0.
t_guess = max(t_guess, decay_estimate)
I_burnup = data[-1][0]*La/data[-1][1] if a.reaction == "b" else 0.
t_guess_k = -log(target/(intensity + I_burnup))/La if intensity + I_burnup > 0. else MIN_TIME
t_guess = max(t_guess, t_guess_k)
#print("corrected at t=0", [Ia for Ia, La, reaction in data])

# No need to fit decay time if initial intensity is below target. Even with two-stage
# decay chains, we are only looking at daughter products with longer half life than
# the parent, so activity of the parent during burnup will decrease more than the
# the activity of the daughter can increase.
if initial_activity <= target:
return 0.


# 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.
# 2026-05-27 PAK: include post exposure granddaughter buildup
# 2026-05-27 PAK: include post exposure daughter burnup
def f(t):
total = 0
total = 0.
for k, (Ia, La, reaction) in enumerate(data):
total += Ia*exp(-La*t)
intensity = Ia*exp(-La*t)
if reaction == "b":
# For delayed β- decay use Bateman equation with parent activity
# and half-life from the previous row in the table.
Ip, Lp, _ = data[k-1]
total += Ip*(exp(-Lp*t) - exp(-La*t))/(1 - Lp/La)
intensity += Ip*(exp(-Lp*t) - exp(-La*t))/(1 - Lp/La)
# print(f" f{k}: {activity}")
total += intensity
# print(f"f({t}) = {total} - {target}\n")
return total - target
def dfdt(t):
total = 0
total = 0.
for k, (Ia, La, reaction) in enumerate(data):
total += -La*Ia*exp(-La*t)
delta = -La*Ia*exp(-La*t)
if reaction == "b":
# For delayed β- decay use Bateman equation with parent activity
# and half-life from the previous row in the table.
Ip, Lp, _ = data[k-1]
total += Ip*(La*exp(-La*t) - Lp*exp(-Lp*t))/(1 - Lp/La)
delta += Ip*(La*exp(-La*t) - Lp*exp(-Lp*t))/(1 - Lp/La)
total += delta
# print(f" df{k}: {d_activity}")
# print(f"dfdt({t}) = {total}\n")
return total
t, ft = find_root(t_guess, f, dfdt, tol=tol)

percent_error = 100*abs(ft)/target
if percent_error > 0.1:
#return 0. # Decay time failed to compute; silently fail
#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
f"Failed to compute decay time correctly ({percent_error:.1g} error)."
f" Please report material and activation parameters.")
raise RuntimeError(msg)

# Return time at least zero hours after removal from the beam.
Expand Down Expand Up @@ -503,10 +522,11 @@ def find_root(
"""
fx = f(x)
for _ in range(maxiter):
# print(f"step {_}: {x=} {fx=} df/dx={df(x)} step={-fx/df(x)}")
if abs(fx) < tol:
break
x -= fx / df(x)
dfx = df(x)
# print(f"step {_}: {x=} {fx=} df/dx={dfx} step={-fx/dfx}")
x -= fx / dfx
fx = f(x)
return x, fx

Expand Down Expand Up @@ -845,7 +865,7 @@ def activity(
else:
result[ai] = [activity*exp(-lam*Ti) for Ti in rest_times]
# 2025-05-27 PAK: Hack to use 151Nd activation intensity rather than
# the 151Pm intensity when computing the buildup of the 151Sm granddaugter.
# the 151Pm intensity when computing the activity of the 151Sm granddaugter.
if ai.daughter != "Pm-151":
last_activity = activity
#print([(Ti, Ai) for Ti, Ai in zip(rest_times, result[ai])])
Expand Down
Loading