Allowing other dose coefficients and adding ICRP119 neutron

openmc/data/effective_dose/dose.py

@@ -2,40 +2,55 @@
 
 import numpy as np
 
-_FILES = (
-    ('electron', 'electrons.txt'),
-    ('helium', 'helium_ions.txt'),
-    ('mu-', 'negative_muons.txt'),
-    ('pi-', 'negative_pions.txt'),
-    ('neutron', 'neutrons.txt'),
-    ('photon', 'photons.txt'),
-    ('photon kerma', 'photons_kerma.txt'),
-    ('mu+', 'positive_muons.txt'),
-    ('pi+', 'positive_pions.txt'),
-    ('positron', 'positrons.txt'),
-    ('proton', 'protons.txt')
-)
-
-_DOSE_ICRP116 = {}
-
-
-def _load_dose_icrp116():
-    """Load effective dose tables from text files"""
-    for particle, filename in _FILES:
-        path = Path(__file__).parent / filename
-        data = np.loadtxt(path, skiprows=3, encoding='utf-8')
-        data[:, 0] *= 1e6   # Change energies to eV
-        _DOSE_ICRP116[particle] = data
-
-
-def dose_coefficients(particle, geometry='AP'):
-    """Return effective dose conversion coefficients from ICRP-116
+import openmc.checkvalue as cv
+
+
+_FILES = {
+    ('icrp116', 'electron'): Path('icrp116') / 'electrons.txt',
+    ('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt',
+    ('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt',
+    ('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt',
+    ('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt',
+    ('icrp116', 'photon'): Path('icrp116') / 'photons.txt',
+    ('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt',
+    ('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt',
+    ('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()}
+
+
+def _load_dose(particle, data_source='icrp116'):
+    """Load effective dose tables from text files
+
+    Parameters
+    ----------
+    particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
+        Incident particle
+    data_source : {'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
 
     This function provides fluence (and air kerma) to effective dose conversion
-    coefficients for various types of external exposures based on values in
-    `ICRP Publication 116 <https://doi.org/10.1016/j.icrp.2011.10.001>`_.
-    Corrected values found in a correigendum are used rather than the values in
-    theoriginal report.
+    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 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>`_
+
 
     Parameters
     ----------
@@ -44,6 +59,8 @@ def dose_coefficients(particle, geometry='AP'):
     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 : {'icrp116', 'icrp119'}
+        The dose conversion data source to use.
 
     Returns
     -------
@@ -54,13 +71,17 @@ def dose_coefficients(particle, geometry='AP'):
         'photon kerma', the coefficients are given in [Sv/Gy].
 
     """
-    if not _DOSE_ICRP116:
-        _load_dose_icrp116()
+
+    cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
+    cv.check_value('data_source', data_source, {'icrp116', 'icrp119'})
+
+    if _DOSE_TABLES[data_source, particle] is None:
+        _load_dose(data_source=data_source, particle=particle)
 
     # Get all data for selected particle
-    data = _DOSE_ICRP116.get(particle)
+    data = _DOSE_TABLES[data_source, particle]
     if data is None:
-        raise ValueError(f"{particle} has no effective dose data")
+        raise ValueError(f"{particle} has no effective dose data in data source {data_source}.")
 
     # Determine index for selected geometry
     if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
@@ -71,4 +92,8 @@ def dose_coefficients(particle, geometry='AP'):
     # 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

openmc/data/effective_dose/icrp119/neutrons.txt

@@ -0,0 +1,57 @@
+Neutrons: Effective dose per fluence, in units of pSv cm², for monoenergetic particles incident in various geometries.
+
+Energy (MeV) AP PA LLAT RLAT ROT ISO
+1.0E-9 5.24 3.52 1.36 1.68 2.99 2.40
+1.0E-8 6.55 4.39 1.7 2.04 3.72 2.89
+2.5E-8 7.60 5.16 1.99 2.31 4.40 3.30
+1.0E-7 9.95 6.77 2.58 2.86 5.75 4.13
+2.0E-7 11.2 7.63 2.92 3.21 6.43 4.59
+5.0E-7 12.8 8.76 3.35 3.72 7.27 5.20
+1.0E-6 13.8 9.55 3.67 4.12 7.84 5.63
+2.0E-6 14.5 10.2 3.89 4.39 8.31 5.96
+5.0E-6 15.0 10.7 4.08 4.66 8.72 6.28
+1.0E-5 15.1 11.0 4.16 4.80 8.90 6.44
+2.0E-5 15.1 11.1 4.20 4.89 8.92 6.51
+5.0E-5 14.8 11.1 4.19 4.95 8.82 6.51
+1.0E-4 14.6 11.0 4.15 4.95 8.69 6.45
+2.0E-4 14.4 10.9 4.10 4.92 8.56 6.32
+5.0E-4 14.2 10.7 4.03 4.86 8.40 6.14
+0.001 14.2 10.7 4.00 4.84 8.34 6.04
+0.002 14.4 10.8 4.00 4.87 8.39 6.05
+0.005 15.7 11.6 4.29 5.25 9.06 6.52
+0.01 18.3 13.5 5.02 6.14 10.6 7.70
+0.02 23.8 17.3 6.48 7.95 13.8 10.2
+0.03 29.0 21.0 7.93 9.74 16.9 12.7
+0.05 38.5 27.6 10.6 13.1 22.7 17.3
+0.07 47.2 33.5 13.1 16.1 27.8 21.5
+0.1 59.8 41.3 16.4 20.1 34.8 27.2
+0.15 80.2 52.2 21.2 25.5 45.4 35.2
+0.2 99.0 61.5 25.6 30.3 54.8 42.4
+0.3 133 77.1 33.4 38.6 71.6 54.7
+0.5 188 103 46.8 53.2 99.4 75.0
+0.7 231 124 58.3 66.6 123 92.8
+0.9 267 144 69.1 79.6 144 108
+1.0 282 154 74.5 86.0 154 116
+1.2 310 175 85.8 99.8 173 130
+2.0 383 247 129 153 234 178
+3.0 432 308 171 195 283 220
+4.0 458 345 198 224 315 250
+5.0 474 366 217 244 335 272
+6.0 483 380 232 261 348 282
+7.0 490 391 244 274 358 290
+8.0 494 399 253 285 366 297
+9.0 497 406 261 294 373 303
+10.0 499 412 268 302 378 309
+12.0 499 422 278 315 385 322
+14.0 496 429 286 324 390 333
+15.0 494 431 290 328 391 338
+16.0 491 433 293 331 393 342
+18.0 486 435 299 335 394 345
+20.0 480 436 305 338 395 343
+30.0 458 437 324 nan 395 nan
+50.0 437 444 358 nan 404 nan
+75.0 429 459 397 nan 422 nan
+100.0 429 477 433 nan 443 nan
+130.0 432 495 467 nan 465 nan
+150.0 438 514 501 nan 489 nan
+180.0 445 535 542 nan 517 nan

setup.py

@@ -29,7 +29,8 @@
     'package_data': {
         'openmc.lib': ['libopenmc.{}'.format(suffix)],
         'openmc.data': ['mass_1.mas20.txt', 'BREMX.DAT', 'half_life.json', '*.h5'],
-        'openmc.data.effective_dose': ['*.txt']
+        'openmc.data.effective_dose.icrp116': ['*.txt'],
+        'openmc.data.effective_dose.icrp119': ['*.txt']
     },
 
     # Metadata

tests/unit_tests/test_data_dose.py

@@ -4,19 +4,31 @@
 
 def test_dose_coefficients():
     # Spot checks on values from ICRP tables
-    energy, dose = dose_coefficients('photon', 'AP')
+    energy, dose = dose_coefficients('photon', 'AP', 'icrp116')
     assert energy[0] == approx(0.01e6)
     assert dose[0] == approx(0.0685)
     assert energy[-1] == approx(10e9)
     assert dose[-1] == approx(90.4)  # updated in corrigendum
 
-    energy, dose = dose_coefficients('neutron', 'LLAT')
+    energy, dose = dose_coefficients('neutron', 'LLAT', 'icrp116')
     assert energy[0] == approx(1e-3)
     assert dose[0] == approx(1.04)
     assert energy[-1] == approx(10e9)
     assert dose[-1] == approx(1.23e3)
 
-    energy, dose = dose_coefficients('electron', 'ISO')
+    energy, dose = dose_coefficients('neutron', 'LLAT', 'icrp119')
+    assert energy[0] == approx(1e-3)
+    assert dose[0] == approx(1.36)
+    assert energy[-1] == approx(180e6)
+    assert dose[-1] == approx(542.0)
+
+    # the ISO column in icrp119 has NaN values, this test checks they are removed
+    energy, dose = dose_coefficients('neutron', 'ISO', 'icrp119')
+    assert energy[-1] == approx(20e6)
+    assert dose[-1] == approx(343)
+    assert len(energy) == len(dose)
+
+    energy, dose = dose_coefficients('electron', 'ISO', 'icrp116')
     assert energy[0] == approx(0.01e6)
     assert dose[0] == approx(0.0188)
     assert energy[-1] == approx(10e9)
@@ -27,3 +39,5 @@ def test_dose_coefficients():
         dose_coefficients('slime', 'LAT')
     with raises(ValueError):
         dose_coefficients('neutron', 'ZZ')
+    with raises(ValueError):
+        dose_coefficients('neutron', 'ISO', 'foo')