diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 48797ef5d..9631f055f 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -60,6 +60,7 @@ Core Functions :template: myfunction.rst atomic_mass + dose_coefficients gnd_name linearize thin diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 277c14ca8..c2d35565a 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -33,3 +33,5 @@ from .resonance_covariance import * from .multipole import * from .grid import * from .function import * + +from .effective_dose.dose import dose_coefficients diff --git a/openmc/data/effective_dose/__init__.py b/openmc/data/effective_dose/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/openmc/data/effective_dose/dose.py b/openmc/data/effective_dose/dose.py new file mode 100644 index 000000000..9f8ce87ae --- /dev/null +++ b/openmc/data/effective_dose/dose.py @@ -0,0 +1,71 @@ +from pathlib import Path + +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) + data[:, 0] *= 1e-6 # 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 to dose conversion coefficients for effective + dose for various types of external exposures based on values in `ICRP + Publication 116 `_. + + Parameters + ---------- + particle : {'neutron', 'photon', 'electron', 'positron'} + Incident particle + geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'} + Irradiation geometry assumed. Refer to ICRP-116 for the meaning of the + options here. + + Returns + ------- + energy : numpy.ndarray + Energies at which dose conversion coefficients are given + dose : numpy.ndarray + Effective dose in [pSv cm^2] at provided energies + + """ + if not _DOSE_ICRP116: + _load_dose_icrp116() + + # Get all data for selected particle + data = _DOSE_ICRP116.get(particle) + if data is None: + raise ValueError("{} has no effective dose data".format(particle)) + + # Determine index for selected geometry + if particle in ('neutron', 'photon', 'proton'): + index = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO').index(geometry) + else: + index = ('AP', 'PA', 'ISO').index(geometry) + + # Pull out energy and dose from table + energy = data[:, 0].copy() + dose = data[:, index + 1].copy() + return energy, dose diff --git a/tests/unit_tests/test_data_dose.py b/tests/unit_tests/test_data_dose.py new file mode 100644 index 000000000..5b9aafc21 --- /dev/null +++ b/tests/unit_tests/test_data_dose.py @@ -0,0 +1,29 @@ +from openmc.data import dose_coefficients +from pytest import approx, raises + + +def test_dose_coefficients(): + # Spot checks on values from ICRP tables + energy, dose = dose_coefficients('photon', 'AP') + assert energy[0] == approx(0.01e-6) + assert dose[0] == approx(0.0685) + assert energy[-1] == approx(10e-3) + assert dose[-1] == approx(90.4) # updated in corrigendum + + energy, dose = dose_coefficients('neutron', 'LLAT') + assert energy[0] == approx(1e-15) + assert dose[0] == approx(1.04) + assert energy[-1] == approx(10e-3) + assert dose[-1] == approx(1.23e3) + + energy, dose = dose_coefficients('electron', 'ISO') + assert energy[0] == approx(0.01e-6) + assert dose[0] == approx(0.0188) + assert energy[-1] == approx(10e-3) + assert dose[-1] == approx(699.0) + + # Invalid particle/geometry should raise an exception + with raises(ValueError): + dose_coefficients('slime', 'LAT') + with raises(ValueError): + dose_coefficients('neutron', 'ZZ')