Update 150Nd to use the 3-stage Bateman equation for 151Sm levels

[pkienzle]

During irradiation the 151Pm activity is counted twice: once in the calculation of 151Pm from 151Nd and again in the calculation of 151Sm which directly estimated activity from 151Nd ignoring the 151Pm transient.

This PR subtracts the 151Pm activity calculated during the irradiation from the 151Sm activity at t=0 so it is only
counted once. The post irradiation decay calculation now uses the 3-stage Bateman equation to compute 151Sm
activity from 151Nd, 151Pm and 151Sm activity at t=0.

For the decay time estimate we make a simplifying assumption that all 151Nd activity immediately transforms into 151Sm. This leads to some error in the decay time estimate. For example 1 kg 150Nd at 1e8 for 2.77 m estimates 250 days to decay rather than 253 days.

Tested by adjusting exposure time and mass so that 151Sm is just above/just below threshold and again for 151Pm is just above/just below threshold.

The change to the decay calculation happened over several PRs. The full diff is here:

diff --git a/periodictable/activation.py b/periodictable/activation.py
index 767eb88..e3e4989 100644
--- a/periodictable/activation.py
+++ b/periodictable/activation.py
@@ -336,46 +336,128 @@ class Sample:
         """
         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
 
-        # 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[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[k] > 0.0
-            ]
+        # TODO: save intensity at zero separate from self.rest_times
+
+        # Find I(0) for each isotope given I(t) at the first rest time.
+        # Use Bateman equation to compute backward for granddaughter products:
+        #     Id(t) = Id(0)exp(-λd t) + Ip(0)(exp(-λp t) - exp(-λd t))/(1 - λp/λd)
+        # Data is [(I(0), λ, reaction), ...] for each active isotope
+        t = self.rest_times[0]
+        data: list[tuple[float, float, str]] = []
+
         # 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
+        # derivative. Choosing a time that satisfies the longest decay time 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)
+        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]
+            La = LN2/a.Thalf_hrs
+            # print(f"  {k}: {a.daughter} {intensity=} {a.Thalf_hrs=} {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. and intensity > 0.:
+                # print(f"{a.isotope} {a.daughter} {a.Thalf_hrs=} {La=}; {intensity=} at {t=}")
+                intensity = intensity * exp(La*t)
+                if a.reaction == "b":
+                    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))
+
+            # Daughter intensity with delayed beta decay 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.
+            # Note: data[k-1] != data[-1] because we have already appended to data in the loop
+            # Note: "b" is never the first reaction, so data[k-1] is always valid
+            Ib = data[k-1][0]*La/data[k-1][1] if a.reaction == "b" else 0.
+            t_guess_k = -log(target/(intensity + Ib))/La if intensity + Ib > 0. else MIN_TIME
+            # print(f"  {t_guess_k=} {Ib=} {intensity+Ib=} {La=}")
+            t_guess = max(t_guess, t_guess_k)
+        # print("corrected at t=0", [Ia for Ia, La, reaction in data])
+        # print(f"{t_guess=}")
+
+        # No need to fit decay time if initial intensity is below target. Even with delayed
+        # beta decay, we are only looking at daughter products with longer half life than
+        # the parent, so beta decay of the parent decreases its activity 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
-        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, [])
-        t, ft = find_root(0, f, df, tol=tol)
+        # derivative df/dt. f(t) will be zero when activity is at target.
+        # 2026-05-27 PAK: include post exposure delayed beta activation
+        def f(t):
+            total = 0.
+            for k, (Ia, La, reaction) in enumerate(data):
+                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]
+                    intensity += Ip*(exp(-Lp*t) - exp(-La*t))/(1 - Lp/La)
+                # print(f" f{k}: {intensity}")
+                total += intensity
+            # print(f"f({t}) = {total} - {target}\n")
+            return total - target
+        def dfdt(t):
+            total = 0.
+            for k, (Ia, La, reaction) in enumerate(data):
+                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]
+                    delta += Ip*(La*exp(-La*t) - Lp*exp(-Lp*t))/(1 - Lp/La)
+                total += delta
+                # print(f" df{k}: {delta}")
+            # 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. Correct
-        # for time adjustment we used to stablize the fit.
-        return max(t+guess, 0.0)
+
+        # Return time at least zero hours after removal from the beam.
+        return max(t, 0.0)
 
     def _accumulate(self, activity: dict["ActivationResult", list[float]]):
         for el, activity_el in activity.items():
@@ -611,15 +694,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
@@ -660,6 +752,11 @@ def activity(
     if not hasattr(isotope, 'neutron_activation'):
         return result
 
+    # 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. 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:
@@ -707,7 +804,7 @@ def activity(
         if ai.reaction == 'b':
             # Column N: 0.69/t1/2 [1/h] lambda of parent nuclide
             parent_lam = LN2 / ai.Thalf_parent
-            # Column O: Activation if "b" mode production
+            # Column O: Activation if "b" mode production (i.e., delayed beta)
             # 2022-05-18 PAK: addressed the following
             #    Note: problems resulting from precision limitation not addressed
             #    in "b" mode production
@@ -785,12 +882,85 @@ def activity(
             #data = env.fluence, initialXS, flux, root, U, V, W, precision_correction
             #print(" ".join("%.5e"%v for v in data))
 
-        # TODO: chained activity (e.g., )
-        result[ai] = [activity*exp(-lam*Ti) for Ti in rest_times]
+        # 2026-05-27 PAK: delayed beta decay such as 209Bi -> 210Bi -> 210Po -> 206Pb
+        # TODO: check 151Eu -> 152m2Eu -> 152Eu is missing b-mode 152Gd
+        # - It is present for 151Eu -> 152m1Eu isomer and 151Eu -> 152Eu ground state.
+
+        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:
+            result[ai] = [
+                bateman2(Ti, activity, lam, last_activity, last_lam)
+                for Ti in rest_times
+            ]
+            # print(f"{ai.daughter} {activity=} {last_activity=}")
+        else:
+            result[ai] = [activity*exp(-lam*Ti) for Ti in rest_times]
+        last_lam = lam  # For delayed beta computations (with 151Sm, last_lam != parent_lam)
+        last_activity = activity  # For delayed beta computations
         #print([(Ti, Ai) for Ti, Ai in zip(rest_times, result[ai])])
 
     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* :

Fixes #127