Fix up integration of ICRP 74, add one test, remove ICRP 119

openmc/data/effective_dose/dose.py

@@ -18,40 +18,38 @@
     ('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt',
     ('icrp116', 'positron'): Path('icrp116') / 'positrons.txt',
     ('icrp116', 'proton'): Path('icrp116') / 'protons.txt',
-    ('icrp119', 'neutron'): Path('icrp119') / 'neutrons.txt',
 }
 
-_DOSE_TABLES = {key: None for key in _FILES.keys()}
+_DOSE_TABLES = {}
 
-def _load_dose_icrp(particle, data_source):
-    """Load effective dose tables from text files
+
+def _load_dose_icrp(data_source: str, particle: str):
+    """Load effective dose tables from text files.
 
     Parameters
     ----------
+    data_source : {'icrp74', 'icrp116'}
+        The dose conversion data source to use
     particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
         Incident particle
-    data_source : {'icrp74', 'icrp116', 'icrp119'}
-        The dose conversion data source to use
-        
+
     """
-    print(f'loading {particle} {data_source}')
     path = Path(__file__).parent / _FILES[data_source, particle]
     data = np.loadtxt(path, skiprows=3, encoding='utf-8')
     data[:, 0] *= 1e6   # Change energies to eV
     _DOSE_TABLES[data_source, particle] = data
 
 
-def dose_coefficients(particle, geometry='AP', data_source = 'icrp116'):
-
-    """Return effective dose conversion coefficients from ICRP
+def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
+    """Return effective dose conversion coefficients.
 
-    This function provides fluence (and air kerma) to effective or ambient dose (H*(10)) conversion
-    coefficients for various types of external exposures based on values in ICRP
-    publications. Corrected values found in a corrigendum are used rather than
-    the values in the original report. Available libraries include 
-    `ICRP Publication 74 <https://journals.sagepub.com/doi/pdf/10.1177/ANIB_26_3-4>`, 
-    `ICRP Publication 116 <https://doi.org/10.1016/j.icrp.2011.10.001>` and 
-    `ICRP Publication 119 <https://journals.sagepub.com/doi/pdf/10.1016/j.icrp.2013.05.003>`.
+    This function provides fluence (and air kerma) to effective or ambient dose
+    (H*(10)) conversion coefficients for various types of external exposures
+    based on values in ICRP publications. Corrected values found in a
+    corrigendum are used rather than the values in the original report.
+    Available libraries include `ICRP Publication 74
+    <https://doi.org/10.1016/S0146-6453(96)90010-X>` and `ICRP Publication 116
+    <https://doi.org/10.1016/j.icrp.2011.10.001>`.
 
 
     Parameters
@@ -61,8 +59,8 @@ def dose_coefficients(particle, geometry='AP', data_source = 'icrp116'):
     geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
         Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
         meaning of the options here.
-    data_source : {'icrp74', 'icrp116', 'icrp119'}
-        The dose conversion data source to use.
+    data_source : {'icrp74', 'icrp116'}
+        The data source for the effective dose conversion coefficients.
 
     Returns
     -------
@@ -75,39 +73,27 @@ def dose_coefficients(particle, geometry='AP', data_source = 'icrp116'):
     """
 
     cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
-    cv.check_value('data_source', data_source, {'icrp74', 'icrp116', 'icrp119'})
+    cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})
 
-    if _DOSE_TABLES[data_source, particle] is None:
-        _load_dose_icrp(data_source=data_source, particle=particle)
+    if (data_source, particle) not in _FILES:
+        raise ValueError(f"{particle} has no dose data in data source {data_source}.")
+    elif (data_source, particle) not in _DOSE_TABLES:
+        _load_dose_icrp(data_source, particle)
 
     # Get all data for selected particle
     data = _DOSE_TABLES[data_source, particle]
-    if data is None:
-        raise ValueError(f"{particle} has no dose data in data source {data_source}.")
-
-    # Determine index for selected geometry
-    if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
-        index = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO').index(geometry)
-    else:
-        index = ('AP', 'PA', 'ISO').index(geometry)
-
-    # Immediately return the ambient dose coefficients if the ICRP74 is selected
-    if data_source == "icrp74":
-        energy = data[:, 0].copy()
-        dose_coeffs = data[:,1].copy()
-        return energy, dose_coeffs        
 
     # Determine index for selected geometry
     if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
-        index = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO').index(geometry)
+        if data_source == 'icrp74':
+            columns = ('AP', 'PA', 'RLAT', 'LLAT', 'ROT', 'ISO')
+        elif data_source == 'icrp116':
+            columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
     else:
-        index = ('AP', 'PA', 'ISO').index(geometry)
+        columns = ('AP', 'PA', 'ISO')
+    index = columns.index(geometry)
 
     # Pull out energy and dose from table
     energy = data[:, 0].copy()
     dose_coeffs = data[:, index + 1].copy()
-    # icrp119 neutron does have NaN values in them
-    if data_source == 'icrp119' and particle == 'neutron' and geometry in ['ISO', 'RLAT']:
-        dose_coeffs = dose_coeffs[~np.isnan(dose_coeffs)]
-        energy = energy[:len(dose_coeffs)]
-    return energy, dose_coeffs
+    return energy, dose_coeffs

tests/unit_tests/test_data_dose.py

@@ -22,8 +22,16 @@ def test_dose_coefficients():
     assert energy[-1] == approx(10e9)
     assert dose[-1] == approx(699.0)
 
+    energy, dose = dose_coefficients('photon', data_source='icrp74')
+    assert energy[0] == approx(0.01e6)
+    assert dose[0] == approx(0.061)
+    assert energy[-1] == approx(10.0e6)
+    assert dose[-1] == approx(25.6)
+
     # Invalid particle/geometry should raise an exception
     with raises(ValueError):
         dose_coefficients('slime', 'LAT')
     with raises(ValueError):
         dose_coefficients('neutron', 'ZZ')
+    with raises(ValueError):
+        dose_coefficients('neutron', data_source='icrp7000')