Merge pull request #339 from lbl-anp/fix-stable-isoqty

Fix stable IsotopeQuantity
lbl-anp · May 16, 2022 · 1758ebb · 1758ebb
2 parents a54659b + fc51a70
commit 1758ebb
Show file tree

Hide file tree

Showing 3 changed files with 76 additions and 44 deletions.
diff --git a/becquerel/tools/isotope_qty.py b/becquerel/tools/isotope_qty.py
@@ -163,16 +163,24 @@ def _quantities_from_kwargs(self, **kwargs):
 
         # dictionary with functions that define how to calculate all quantities
         # in a circular manner
-        conversions = dict(
-            atoms=lambda: ref_quantities["g"] / self.isotope.A * N_AV,
-            bq=lambda: ref_quantities["atoms"] * self.decay_const,
-            uci=lambda: ref_quantities["bq"] / UCI_TO_BQ,
-            g=lambda: ref_quantities["uci"]
-            * UCI_TO_BQ
-            / self.decay_const
-            / N_AV
-            * self.isotope.A,
-        )
+        if self.is_stable:
+            conversions = dict(
+                atoms=lambda: ref_quantities["g"] / self.isotope.A * N_AV,
+                g=lambda: ref_quantities["atoms"] * self.isotope.A / N_AV,
+                bq=lambda: 0,
+                uci=lambda: 0,
+            )
+        else:
+            conversions = dict(
+                atoms=lambda: ref_quantities["g"] / self.isotope.A * N_AV,
+                bq=lambda: ref_quantities["atoms"] * self.decay_const,
+                uci=lambda: ref_quantities["bq"] / UCI_TO_BQ,
+                g=lambda: ref_quantities["uci"]
+                * UCI_TO_BQ
+                / self.decay_const
+                / N_AV
+                * self.isotope.A,
+            )
 
         # rotates the order of the list so that the provided kwarg is at [0]
         order = ["atoms", "bq", "uci", "g"]

diff --git a/examples/isotopes.ipynb b/examples/isotopes.ipynb
@@ -148,7 +148,7 @@
       "Abundance:   0.01+/-0.00%\n",
       "Stable?      False\n",
       "Half-life:   1.248e+09 years\n",
-      "Decay modes: (['B-', 'EC'], [89.28, 10.72])\n",
+      "Decay modes: (['EC', 'B-'], [10.72, 89.28])\n",
       "\n",
       "Isotope: K-41\n",
       "Spin-parity: 3/2+\n",
@@ -202,7 +202,7 @@
     {
      "data": {
       "text/plain": [
-       "8.474020543863288"
+       "8.00631769726367"
       ]
      },
      "execution_count": 6,
@@ -231,12 +231,12 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Ba-133: 4.707 uCi\n",
-      "Cd-109: 0.016 uCi\n",
+      "Ba-133: 4.004 uCi\n",
+      "Cd-109: 0.004 uCi\n",
       "Co-57: 0.000 uCi\n",
-      "Co-60: 1.978 uCi\n",
-      "Mn-54: 0.001 uCi\n",
-      "Na-22: 0.437 uCi\n"
+      "Co-60: 1.431 uCi\n",
+      "Mn-54: 0.000 uCi\n",
+      "Na-22: 0.226 uCi\n"
      ]
     }
    ],
@@ -268,9 +268,9 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "37000.00000000001 Bq of Co-60 (at 2017-03-16 14:00:00)\n",
-      "73926311.79156348 Bq of K-40 (at 2019-11-14 16:16:17.484865)\n",
-      "4104.821259246284 Bq of U-238 (at 2019-11-14 16:16:17.495673)\n"
+      "37000.0 Bq of Co-60 (at 2017-03-16 14:00:00)\n",
+      "73926311.7915635 Bq of K-40 (at 2022-05-02 14:53:32.265827)\n",
+      "4104.821259246285 Bq of U-238 (at 2022-05-02 14:53:32.272388)\n"
      ]
     }
    ],
@@ -297,7 +297,7 @@
     {
      "data": {
       "text/plain": [
-       "3.831751598190806"
+       "3.8317515981908064"
       ]
      },
      "execution_count": 9,
@@ -317,7 +317,7 @@
     {
      "data": {
       "text/plain": [
-       "576.9086036541194"
+       "149.24053998847538"
       ]
      },
      "execution_count": 10,
@@ -357,7 +357,7 @@
     {
      "data": {
       "text/plain": [
-       "229361876.6085581"
+       "23101885.23242553"
       ]
      },
      "execution_count": 12,
@@ -494,7 +494,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "11804946.77955213 Bq of Na-24 (at 2017-04-30 11:32:00)\n"
+      "11804946.779552134 Bq of Na-24 (at 2017-04-30 11:32:00)\n"
      ]
     }
    ],
@@ -524,7 +524,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "2.1776251393191668e-07 g of Na-23\n"
+      "2.177625139319166e-07 g of Na-23\n"
      ]
     }
    ],
@@ -536,30 +536,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 19,
    "metadata": {},
-   "outputs": [],
-   "source": []
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2.177625139319166e-07"
+      ]
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "na23_qty.g_at()"
+   ]
   }
  ],
  "metadata": {
-  "anaconda-cloud": {},
-  "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
-  },
   "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
-   "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "name": "python"
   }
  },
  "nbformat": 4,

diff --git a/tests/isotope_qty_test.py b/tests/isotope_qty_test.py
@@ -101,6 +101,12 @@ def test_isotopequantity_init_stable(stable_isotope, iq_date, iq_kwargs):
     assert iq.isotope is stable_isotope
     assert iq.half_life == stable_isotope.half_life
     assert iq.decay_const == stable_isotope.decay_const
+    assert iq.bq_at() == iq.bq_at(iq.ref_date) == 0.0
+    assert iq.uci_at() == iq.uci_at(iq.ref_date) == 0.0
+    if "g" in iq_kwargs:
+        assert np.isclose(iq.g_at(), iq_kwargs["g"])
+    if "atoms" in iq_kwargs:
+        assert np.isclose(iq.atoms_at(), iq_kwargs["atoms"])
 
 
 def test_isotopequantity_quantity_at(radioisotope, iq_kwargs):
@@ -545,6 +551,26 @@ def test_irradiation_str(start, stop, n_cm2, n_cm2_s):
     print(str(ni))
 
 
+def test_irradiation_nan_g():
+    """Regression test for https://github.com/lbl-anp/becquerel/issues/325.
+
+    From examples/isotopes.ipynb
+    """
+    iso = Isotope("Na-23")
+    iso2 = Isotope("Na-24")
+    barns = 2.3  # making this up for now
+
+    irradiation_start = "2017-04-30 10:32:00"
+    irradiation_stop = "2017-04-30 11:32:00"
+    flux = 3.1e11
+    ni = NeutronIrradiation(irradiation_start, irradiation_stop, n_cm2_s=flux)
+
+    activated_qty = IsotopeQuantity(iso2, date="2017-05-01", bq=103.2)
+    na23_qty = ni.activate(barns, initial=iso, activated=activated_qty)
+
+    assert not np.isnan(na23_qty.g_at())
+
+
 # ----------------------------------------------------
 #               decay_normalize
 # ----------------------------------------------------