Add dose coefficients from ICRP 74 (#3020)

Co-authored-by: matteo.zammataro <matteo.zammataro@newcleo.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

openmc/data/effective_dose/dose.py

@@ -2,40 +2,61 @@
 
 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
-
-    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.
+import openmc.checkvalue as cv
+
+_FILES = {
+    ('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
+    ('icrp74', 'photon'): Path('icrp74') / 'photons.txt',
+    ('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',
+}
+
+_DOSE_TABLES = {}
+
+
+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
+
+    """
+    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.
+
+    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>`.
+
+    For ICRP 74 data, the photon effective dose per fluence is determined by
+    multiplying the air kerma per fluence values (Table A.1) by the effective
+    dose per air kerma (Table A.17). The neutron effective dose per fluence is
+    found in Table A.41. For ICRP 116 data, the photon effective dose per
+    fluence is found in Table A.1 and the neutron effective dose per fluence is
+    found in Table A.5.
 
     Parameters
     ----------
@@ -44,6 +65,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 : {'icrp74', 'icrp116'}
+        The data source for the effective dose conversion coefficients.
 
     Returns
     -------
@@ -54,19 +77,24 @@ 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, {'icrp74', 'icrp116'})
+
+    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_ICRP116.get(particle)
-    if data is None:
-        raise ValueError(f"{particle} has no effective dose data")
+    data = _DOSE_TABLES[data_source, particle]
 
     # Determine index for selected geometry
     if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
-        index = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO').index(geometry)
+        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()

openmc/data/effective_dose/icrp74/generate_photon_effective_dose.py

@@ -0,0 +1,69 @@
+from prettytable import PrettyTable
+import numpy as np
+
+# Data from Table A.1 (air kerma per fluence)
+energy_a1 = np.array([
+    0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.15, 0.2,
+    0.3, 0.4, 0.5, 0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0
+])
+air_kerma = np.array([7.43, 3.12, 1.68, 0.721, 0.429, 0.323, 0.289, 0.307, 0.371, 0.599, 0.856, 1.38,
+               1.89, 2.38, 2.84, 3.69, 4.47, 6.14, 7.55, 9.96, 12.1, 14.1, 16.1, 20.1, 24.0])
+
+# Data from Table A.17 (effective dose per air kerma)
+energy_a17 = np.array([
+    0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1, 0.15, 0.2, 0.3,
+    0.4, 0.5, 0.6, 0.8, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0
+])
+dose_per_airkerma = {
+    'AP': np.array([
+        0.00653, 0.0402, 0.122, 0.416, 0.788, 1.106, 1.308, 1.407, 1.433, 1.394,
+        1.256, 1.173, 1.093, 1.056, 1.036, 1.024, 1.010, 1.003, 0.992, 0.993,
+        0.993, 0.991, 0.990
+    ]),
+    'PA': np.array([
+        0.00248, 0.00586, 0.0181, 0.128, 0.370, 0.640, 0.846, 0.966, 1.019,
+        1.030, 0.959, 0.915, 0.880, 0.871, 0.869, 0.870, 0.875, 0.880, 0.901,
+        0.918, 0.924, 0.927, 0.929
+    ]),
+    'RLAT': np.array([
+        0.00172, 0.00549, 0.0151, 0.0782, 0.205, 0.345, 0.455, 0.522, 0.554,
+        0.571, 0.551, 0.549, 0.557, 0.570, 0.585, 0.600, 0.628, 0.651, 0.728,
+        0.796, 0.827, 0.846, 0.860
+    ]),
+    'LLAT': np.array([
+        0.00172, 0.00549, 0.0155, 0.0904, 0.241, 0.405, 0.528, 0.598, 0.628,
+        0.641, 0.620, 0.615, 0.615, 0.623, 0.635, 0.648, 0.670, 0.691, 0.757,
+        0.813, 0.836, 0.850, 0.859
+    ]),
+    'ROT': np.array([
+        0.00326, 0.0153, 0.0462, 0.191, 0.426, 0.661, 0.828, 0.924, 0.961,
+        0.960, 0.892, 0.854, 0.824, 0.814, 0.812, 0.814, 0.821, 0.831, 0.871,
+        0.909, 0.925, 0.934, 0.941
+    ]),
+    'ISO': np.array([
+        0.00271, 0.0123, 0.0362, 0.143, 0.326, 0.511, 0.642, 0.720, 0.749,
+        0.748, 0.700, 0.679, 0.664, 0.667, 0.675, 0.684, 0.703, 0.719, 0.774,
+        0.824, 0.846, 0.859, 0.868
+    ])
+}
+
+# Interpolate air kerma onto energy grid for Table A.17
+air_kerma = np.interp(energy_a17, energy_a1, air_kerma)
+
+# Compute effective dose per fluence
+dose_per_fluence = {
+    geometry: air_kerma * dose_per_airkerma
+    for geometry, dose_per_airkerma in dose_per_airkerma.items()
+}
+
+# Create table
+table = PrettyTable()
+table.field_names = ['Energy (MeV)', 'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO']
+table.float_format = '.7'
+for i, energy in enumerate(energy_a17):
+    row = [energy]
+    for geometry in table.field_names[1:]:
+        row.append(dose_per_fluence[geometry][i])
+    table.add_row(row)
+print('Photons: Effective dose per fluence, in units of pSv cm², for monoenergetic particles incident in various geometries.\n')
+print(table.get_string(border=False))

openmc/data/effective_dose/icrp74/neutrons.txt

@@ -0,0 +1,50 @@
+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.00E-09        5.24    3.52    1.68        1.36        2.99    2.4
+1.00E-08        6.55    4.39    2.04        1.7         3.72    2.89
+2.50E-08        7.6     5.16    2.31        1.99        4.4     3.3
+1.00E-07        9.95    6.77    2.86        2.58        5.75    4.13
+2.00E-07        11.2    7.63    3.21        2.92        6.43    4.59
+5.00E-07        12.8    8.76    3.72        3.35        7.27    5.2
+1.00E-06        13.8    9.55    4.12        3.67        7.84    5.63
+2.00E-06        14.5    10.2    4.39        3.89        8.31    5.96
+5.00E-06        15      10.7    4.66        4.08        8.72    6.28
+1.00E-05        15.1    11      4.8         4.16        8.9     6.44
+2.00E-05        15.1    11.1    4.89        4.2         8.92    6.51
+5.00E-05        14.8    11.1    4.95        4.19        8.82    6.51
+1.00E-04        14.6    11      4.95        4.15        8.69    6.45
+2.00E-04        14.4    10.9    4.92        4.1         8.56    6.32
+5.00E-04        14.2    10.7    4.86        4.03        8.4     6.14
+1.00E-03        14.2    10.7    4.84        4           8.34    6.04
+2.00E-03        14.4    10.8    4.87        4           8.39    6.05
+5.00E-03        15.7    11.6    5.25        4.29        9.06    6.52
+1.00E-02        18.3    13.5    6.14        5.02        10.6    7.7
+2.00E-02        23.8    17.3    7.95        6.48        13.8    10.2
+3.00E-02        29      21      9.74        7.93        16.9    12.7
+5.00E-02        38.5    27.6    13.1        10.6        22.7    17.3
+7.00E-02        47.2    33.5    16.1        13.1        27.8    21.5
+1.00E-01        59.8    41.3    20.1        16.4        34.8    27.2
+1.50E-01        80.2    52.2    25.5        21.2        45.4    35.2
+2.00E-01        99      61.5    30.3        25.6        54.8    42.4
+3.00E-01        133     77.1    38.6        33.4        71.6    54.7
+5.00E-01        188     103     53.2        46.8        99.4    75
+7.00E-01        231     124     66.6        58.3        123     92.8
+9.00E-01        267     144     79.6        69.1        144     108
+1               282     154     86          74.5        154     116
+1.2             310     175     99.8        85.8        173     130
+2               383     247     153         129         234     178
+3               432     308     195         171         283     220
+4               458     345     224         198         315     250
+5               474     366     244         217         335     272
+6               483     380     261         232         348     282
+7               490     391     274         244         358     290
+8               494     399     285         253         366     297
+9               497     406     294         261         373     303
+1.00E+01        499     412     302         268         378     309
+1.20E+01        499     422     315         278         385     322
+1.40E+01        496     429     324         286         390     333
+1.50E+01        494     431     328         290         391     338
+1.60E+01        491     433     331         293         393     342
+1.80E+01        486     435     335         299         394     345
+2.00E+01        480     436     338         305         395     343

openmc/data/effective_dose/icrp74/photons.txt

@@ -0,0 +1,26 @@
+Photons: Effective dose per fluence, in units of pSv cm², for monoenergetic particles incident in various geometries.
+
+ Energy (MeV)      AP          PA         LLAT        RLAT        ROT         ISO
+  0.0100000    0.0485179   0.0184264   0.0127796   0.0127796   0.0242218   0.0201353
+  0.0150000    0.1254240   0.0182832   0.0171288   0.0171288   0.0477360   0.0383760
+  0.0200000    0.2049600   0.0304080   0.0260400   0.0253680   0.0776160   0.0608160
+  0.0300000    0.2999360   0.0922880   0.0651784   0.0563822   0.1377110   0.1031030
+  0.0400000    0.3380520   0.1587300   0.1033890   0.0879450   0.1827540   0.1398540
+  0.0500000    0.3572380   0.2067200   0.1308150   0.1114350   0.2135030   0.1650530
+  0.0600000    0.3780120   0.2444940   0.1525920   0.1314950   0.2392920   0.1855380
+  0.0700000    0.4192860   0.2878680   0.1782040   0.1555560   0.2753520   0.2145600
+  0.0800000    0.4399310   0.3128330   0.1927960   0.1700780   0.2950270   0.2299430
+  0.1000000    0.5171740   0.3821300   0.2378110   0.2118410   0.3561600   0.2775080
+  0.1500000    0.7523440   0.5744410   0.3713800   0.3300490   0.5343080   0.4193000
+  0.2000000    1.0040880   0.7832400   0.5264400   0.4699440   0.7310240   0.5812240
+  0.3000000    1.5083400   1.2144000   0.8487000   0.7686600   1.1371200   0.9163200
+  0.4000000    1.9958400   1.6461900   1.1774700   1.0773000   1.5384600   1.2606300
+  0.5000000    2.4656800   2.0682200   1.5113000   1.3923000   1.9325600   1.6065000
+  0.6000000    2.9081600   2.4708000   1.8403200   1.7040000   2.3117600   1.9425600
+  0.8000000    3.7269000   3.2287500   2.4723000   2.3173200   3.0294900   2.5940700
+  1.0000000    4.4834100   3.9336000   3.0887700   2.9099700   3.7145700   3.2139300
+  2.0000000    7.4896000   6.8025500   5.7153500   5.4964000   6.5760500   5.8437000
+  4.0000000    12.0153000  11.1078000  9.8373000   9.6316000   10.9989000  9.9704000
+  6.0000000    15.9873000  14.8764000  13.4596000  13.3147000  14.8925000  13.6206000
+  8.0000000    19.9191000  18.6327000  17.0850000  17.0046000  18.7734000  17.2659000
+  10.0000000   23.7600000  22.2960000  20.6160000  20.6400000  22.5840000  20.8320000

tests/unit_tests/test_data_dose.py

@@ -22,8 +22,22 @@ 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(7.43*0.00653)
+    assert energy[-1] == approx(10.0e6)
+    assert dose[-1] == approx(24.0*0.990)
+
+    energy, dose = dose_coefficients('neutron', 'LLAT', data_source='icrp74')
+    assert energy[0] == approx(1e-3)
+    assert dose[0] == approx(1.68)
+    assert energy[-1] == approx(20.0e6)
+    assert dose[-1] == approx(338.0)
+
     # 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')