Add air dose response to dose_coefficients

openmc/data/effective_dose/dose.py

@@ -1,8 +1,12 @@
 from pathlib import Path
 
 import numpy as np
+from scipy.interpolate import interp1d
+from scipy.constants import eV, gram, pico
 
 import openmc.checkvalue as cv
+from .nist126 import NIST126_AIR
+from openmc.data import linearize
 
 _FILES = {
     ('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
@@ -60,7 +64,7 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
 
     Parameters
     ----------
-    particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
+    particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron', 'air kerma'}
         Incident particle
     geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
         Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
@@ -73,10 +77,19 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
     energy : numpy.ndarray
         Energies at which dose conversion coefficients are given
     dose_coeffs : numpy.ndarray
-        Effective dose coefficients in [pSv cm^2] at provided energies. For
-        'photon kerma', the coefficients are given in [Sv/Gy].
+        Effective dose coefficients in [pSv cm^2] at provided energies. 
+        For 'photon kerma', the coefficients are given in [Sv/Gy].
+        For 'air kerma', the coefficients are given in [pGy cm^2].
 
     """
+
+    if particle == 'air kerma':
+      interp = interp1d(NIST126_AIR[:,0], NIST126_AIR[:,1])
+
+      def func(e):
+          return e*interp(e)*eV/gram/pico
+
+      return linearize(interp.x, func)
 
     cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
     cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})

openmc/data/effective_dose/nist126.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+from openmc.data import EV_PER_MEV
+
+# Embedded NIST-126 data
+# Air (Dry Near Sea Level) — NIST Standard Reference Database 126 Table 4 (doi: 10.18434/T4D01F)
+# Columns: Energy (MeV), μen/ρ (cm^2/g)
+NIST126_AIR = np.array(
+    [
+        [1.00000e-03, 3.599e03],
+        [1.50000e-03, 1.188e03],
+        [2.00000e-03, 5.262e02],
+        [3.00000e-03, 1.614e02],
+        [3.20290e-03, 1.330e02],
+        [3.20290e-03, 1.460e02],
+        [4.00000e-03, 7.636e01],
+        [5.00000e-03, 3.931e01],
+        [6.00000e-03, 2.270e01],
+        [8.00000e-03, 9.446e00],
+        [1.00000e-02, 4.742e00],
+        [1.50000e-02, 1.334e00],
+        [2.00000e-02, 5.389e-01],
+        [3.00000e-02, 1.537e-01],
+        [4.00000e-02, 6.833e-02],
+        [5.00000e-02, 4.098e-02],
+        [6.00000e-02, 3.041e-02],
+        [8.00000e-02, 2.407e-02],
+        [1.00000e-01, 2.325e-02],
+        [1.50000e-01, 2.496e-02],
+        [2.00000e-01, 2.672e-02],
+        [3.00000e-01, 2.872e-02],
+        [4.00000e-01, 2.949e-02],
+        [5.00000e-01, 2.966e-02],
+        [6.00000e-01, 2.953e-02],
+        [8.00000e-01, 2.882e-02],
+        [1.00000e00, 2.789e-02],
+        [1.25000e00, 2.666e-02],
+        [1.50000e00, 2.547e-02],
+        [2.00000e00, 2.345e-02],
+        [3.00000e00, 2.057e-02],
+        [4.00000e00, 1.870e-02],
+        [5.00000e00, 1.740e-02],
+        [6.00000e00, 1.647e-02],
+        [8.00000e00, 1.525e-02],
+        [1.00000e01, 1.450e-02],
+        [1.50000e01, 1.353e-02],
+        [2.00000e01, 1.311e-02],
+    ],
+    dtype=float,
+)
+NIST126_AIR[:, 0] *= EV_PER_MEV

tests/unit_tests/test_data_dose.py

@@ -27,6 +27,12 @@ def test_dose_coefficients():
     assert dose[0] == approx(7.43*0.00653)
     assert energy[-1] == approx(10.0e6)
     assert dose[-1] == approx(24.0*0.990)
+
+    energy, dose = dose_coefficients('air kerma')
+    assert energy[0] == approx(1e3)
+    assert dose[0] == approx(576.623)
+    assert energy[-1] == approx(20e6)
+    assert dose[-1] == approx(42.0091)
 
     energy, dose = dose_coefficients('neutron', 'LLAT', data_source='icrp74')
     assert energy[0] == approx(1e-3)