import itertools
import json
import os
import re
from pathlib import Path
from math import sqrt, log
from warnings import warn
from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL,
EV_PER_MEV, K_BOLTZMANN, gnds_name, zam)
gnds_name.__module__ = __name__
zam.__module__ = __name__
# Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
# of the elements 2013 (IUPAC Technical Report)", Pure. Appl. Chem. 88 (3),
# pp. 293-306 (2013). The "representative isotopic abundance" values from
# column 9 are used except where an interval is given, in which case the
# "best measurement" is used.
# Note that the abundances are given as atomic fractions!
NATURAL_ABUNDANCE = {
'H1': 0.99984426, 'H2': 0.00015574, 'He3': 0.000002,
'He4': 0.999998, 'Li6': 0.07589, 'Li7': 0.92411,
'Be9': 1.0, 'B10': 0.1982, 'B11': 0.8018,
'C12': 0.988922, 'C13': 0.011078, 'N14': 0.996337,
'N15': 0.003663, 'O16': 0.9976206, 'O17': 0.000379,
'O18': 0.0020004, 'F19': 1.0, 'Ne20': 0.9048,
'Ne21': 0.0027, 'Ne22': 0.0925, 'Na23': 1.0,
'Mg24': 0.78951, 'Mg25': 0.1002, 'Mg26': 0.11029,
'Al27': 1.0, 'Si28': 0.9222968, 'Si29': 0.0468316,
'Si30': 0.0308716, 'P31': 1.0, 'S32': 0.9504074,
'S33': 0.0074869, 'S34': 0.0419599, 'S36': 0.0001458,
'Cl35': 0.757647, 'Cl37': 0.242353, 'Ar36': 0.003336,
'Ar38': 0.000629, 'Ar40': 0.996035, 'K39': 0.932581,
'K40': 0.000117, 'K41': 0.067302, 'Ca40': 0.96941,
'Ca42': 0.00647, 'Ca43': 0.00135, 'Ca44': 0.02086,
'Ca46': 0.00004, 'Ca48': 0.00187, 'Sc45': 1.0,
'Ti46': 0.0825, 'Ti47': 0.0744, 'Ti48': 0.7372,
'Ti49': 0.0541, 'Ti50': 0.0518, 'V50': 0.0025,
'V51': 0.9975, 'Cr50': 0.04345, 'Cr52': 0.83789,
'Cr53': 0.09501, 'Cr54': 0.02365, 'Mn55': 1.0,
'Fe54': 0.05845, 'Fe56': 0.91754, 'Fe57': 0.02119,
'Fe58': 0.00282, 'Co59': 1.0, 'Ni58': 0.680769,
'Ni60': 0.262231, 'Ni61': 0.011399, 'Ni62': 0.036345,
'Ni64': 0.009256, 'Cu63': 0.6915, 'Cu65': 0.3085,
'Zn64': 0.4917, 'Zn66': 0.2773, 'Zn67': 0.0404,
'Zn68': 0.1845, 'Zn70': 0.0061, 'Ga69': 0.60108,
'Ga71': 0.39892, 'Ge70': 0.2052, 'Ge72': 0.2745,
'Ge73': 0.0776, 'Ge74': 0.3652, 'Ge76': 0.0775,
'As75': 1.0, 'Se74': 0.0086, 'Se76': 0.0923,
'Se77': 0.076, 'Se78': 0.2369, 'Se80': 0.498,
'Se82': 0.0882, 'Br79': 0.50686, 'Br81': 0.49314,
'Kr78': 0.00355, 'Kr80': 0.02286, 'Kr82': 0.11593,
'Kr83': 0.115, 'Kr84': 0.56987, 'Kr86': 0.17279,
'Rb85': 0.7217, 'Rb87': 0.2783, 'Sr84': 0.0056,
'Sr86': 0.0986, 'Sr87': 0.07, 'Sr88': 0.8258,
'Y89': 1.0, 'Zr90': 0.5145, 'Zr91': 0.1122,
'Zr92': 0.1715, 'Zr94': 0.1738, 'Zr96': 0.028,
'Nb93': 1.0, 'Mo92': 0.14649, 'Mo94': 0.09187,
'Mo95': 0.15873, 'Mo96': 0.16673, 'Mo97': 0.09582,
'Mo98': 0.24292, 'Mo100': 0.09744, 'Ru96': 0.0554,
'Ru98': 0.0187, 'Ru99': 0.1276, 'Ru100': 0.126,
'Ru101': 0.1706, 'Ru102': 0.3155, 'Ru104': 0.1862,
'Rh103': 1.0, 'Pd102': 0.0102, 'Pd104': 0.1114,
'Pd105': 0.2233, 'Pd106': 0.2733, 'Pd108': 0.2646,
'Pd110': 0.1172, 'Ag107': 0.51839, 'Ag109': 0.48161,
'Cd106': 0.01245, 'Cd108': 0.00888, 'Cd110': 0.1247,
'Cd111': 0.12795, 'Cd112': 0.24109, 'Cd113': 0.12227,
'Cd114': 0.28754, 'Cd116': 0.07512, 'In113': 0.04281,
'In115': 0.95719, 'Sn112': 0.0097, 'Sn114': 0.0066,
'Sn115': 0.0034, 'Sn116': 0.1454, 'Sn117': 0.0768,
'Sn118': 0.2422, 'Sn119': 0.0859, 'Sn120': 0.3258,
'Sn122': 0.0463, 'Sn124': 0.0579, 'Sb121': 0.5721,
'Sb123': 0.4279, 'Te120': 0.0009, 'Te122': 0.0255,
'Te123': 0.0089, 'Te124': 0.0474, 'Te125': 0.0707,
'Te126': 0.1884, 'Te128': 0.3174, 'Te130': 0.3408,
'I127': 1.0, 'Xe124': 0.00095, 'Xe126': 0.00089,
'Xe128': 0.0191, 'Xe129': 0.26401, 'Xe130': 0.04071,
'Xe131': 0.21232, 'Xe132': 0.26909, 'Xe134': 0.10436,
'Xe136': 0.08857, 'Cs133': 1.0, 'Ba130': 0.0011,
'Ba132': 0.001, 'Ba134': 0.0242, 'Ba135': 0.0659,
'Ba136': 0.0785, 'Ba137': 0.1123, 'Ba138': 0.717,
'La138': 0.0008881, 'La139': 0.9991119, 'Ce136': 0.00186,
'Ce138': 0.00251, 'Ce140': 0.88449, 'Ce142': 0.11114,
'Pr141': 1.0, 'Nd142': 0.27153, 'Nd143': 0.12173,
'Nd144': 0.23798, 'Nd145': 0.08293, 'Nd146': 0.17189,
'Nd148': 0.05756, 'Nd150': 0.05638, 'Sm144': 0.0308,
'Sm147': 0.15, 'Sm148': 0.1125, 'Sm149': 0.1382,
'Sm150': 0.0737, 'Sm152': 0.2674, 'Sm154': 0.2274,
'Eu151': 0.4781, 'Eu153': 0.5219, 'Gd152': 0.002,
'Gd154': 0.0218, 'Gd155': 0.148, 'Gd156': 0.2047,
'Gd157': 0.1565, 'Gd158': 0.2484, 'Gd160': 0.2186,
'Tb159': 1.0, 'Dy156': 0.00056, 'Dy158': 0.00095,
'Dy160': 0.02329, 'Dy161': 0.18889, 'Dy162': 0.25475,
'Dy163': 0.24896, 'Dy164': 0.2826, 'Ho165': 1.0,
'Er162': 0.00139, 'Er164': 0.01601, 'Er166': 0.33503,
'Er167': 0.22869, 'Er168': 0.26978, 'Er170': 0.1491,
'Tm169': 1.0, 'Yb168': 0.00123, 'Yb170': 0.02982,
'Yb171': 0.14086, 'Yb172': 0.21686, 'Yb173': 0.16103,
'Yb174': 0.32025, 'Yb176': 0.12995, 'Lu175': 0.97401,
'Lu176': 0.02599, 'Hf174': 0.0016, 'Hf176': 0.0526,
'Hf177': 0.186, 'Hf178': 0.2728, 'Hf179': 0.1362,
'Hf180': 0.3508, 'Ta180_m1': 0.0001201, 'Ta181': 0.9998799,
'W180': 0.0012, 'W182': 0.265, 'W183': 0.1431,
'W184': 0.3064, 'W186': 0.2843, 'Re185': 0.374,
'Re187': 0.626, 'Os184': 0.0002, 'Os186': 0.0159,
'Os187': 0.0196, 'Os188': 0.1324, 'Os189': 0.1615,
'Os190': 0.2626, 'Os192': 0.4078, 'Ir191': 0.373,
'Ir193': 0.627, 'Pt190': 0.00012, 'Pt192': 0.00782,
'Pt194': 0.32864, 'Pt195': 0.33775, 'Pt196': 0.25211,
'Pt198': 0.07356, 'Au197': 1.0, 'Hg196': 0.0015,
'Hg198': 0.1004, 'Hg199': 0.1694, 'Hg200': 0.2314,
'Hg201': 0.1317, 'Hg202': 0.2974, 'Hg204': 0.0682,
'Tl203': 0.29524, 'Tl205': 0.70476, 'Pb204': 0.014,
'Pb206': 0.241, 'Pb207': 0.221, 'Pb208': 0.524,
'Bi209': 1.0, 'Th230': 0.0002, 'Th232': 0.9998,
'Pa231': 1.0, 'U234': 0.000054, 'U235': 0.007204,
'U238': 0.992742
}
DADZ = {
'(n,2nd)': (-3, -1),
'(n,2n)': (-1, 0),
'(n,3n)': (-2, 0),
'(n,na)': (-4, -2),
'(n,n3a)': (-12, -6),
'(n,2na)': (-5, -2),
'(n,3na)': (-6, -2),
'(n,np)': (-1, -1),
'(n,n2a)': (-8, -4),
'(n,2n2a)': (-9, -4),
'(n,nd)': (-2, -1),
'(n,nt)': (-3, -1),
'(n,n3He)': (-3, -2),
'(n,nd2a)': (-10, -5),
'(n,nt2a)': (-11, -5),
'(n,4n)': (-3, 0),
'(n,2np)': (-2, -1),
'(n,3np)': (-3, -1),
'(n,n2p)': (-2, -2),
'(n,npa)': (-5, -3),
'(n,gamma)': (1, 0),
'(n,p)': (0, -1),
'(n,d)': (-1, -1),
'(n,t)': (-2, -1),
'(n,3He)': (-2, -2),
'(n,a)': (-3, -2),
'(n,2a)': (-7, -4),
'(n,3a)': (-11, -6),
'(n,2p)': (-1, -2),
'(n,pa)': (-4, -3),
'(n,t2a)': (-10, -5),
'(n,d2a)': (-9, -5),
'(n,pd)': (-2, -2),
'(n,pt)': (-3, -2),
'(n,da)': (-5, -3),
'(n,5n)': (-4, 0),
'(n,6n)': (-5, 0),
'(n,2nt)': (-4, -1),
'(n,ta)': (-6, -3),
'(n,4np)': (-4, -1),
'(n,3nd)': (-4, -1),
'(n,nda)': (-6, -3),
'(n,2npa)': (-6, -3),
'(n,7n)': (-6, 0),
'(n,8n)': (-7, 0),
'(n,5np)': (-5, -1),
'(n,6np)': (-6, -1),
'(n,7np)': (-7, -1),
'(n,4na)': (-7, -2),
'(n,5na)': (-8, -2),
'(n,6na)': (-9, -2),
'(n,7na)': (-10, -2),
'(n,4nd)': (-5, -1),
'(n,5nd)': (-6, -1),
'(n,6nd)': (-7, -1),
'(n,3nt)': (-5, -1),
'(n,4nt)': (-6, -1),
'(n,5nt)': (-7, -1),
'(n,6nt)': (-8, -1),
'(n,2n3He)': (-4, -2),
'(n,3n3He)': (-5, -2),
'(n,4n3He)': (-6, -2),
'(n,3n2p)': (-4, -2),
'(n,3n2a)': (-10, -4),
'(n,3npa)': (-7, -3),
'(n,dt)': (-4, -2),
'(n,npd)': (-3, -2),
'(n,npt)': (-4, -2),
'(n,ndt)': (-5, -2),
'(n,np3He)': (-4, -3),
'(n,nd3He)': (-5, -3),
'(n,nt3He)': (-6, -3),
'(n,nta)': (-7, -3),
'(n,2n2p)': (-3, -2),
'(n,p3He)': (-4, -3),
'(n,d3He)': (-5, -3),
'(n,3Hea)': (-6, -4),
'(n,4n2p)': (-5, -2),
'(n,4n2a)': (-11, -4),
'(n,4npa)': (-8, -3),
'(n,3p)': (-2, -3),
'(n,n3p)': (-3, -3),
'(n,3n2pa)': (-8, -4),
'(n,5n2p)': (-6, -2),
}
# Values here are from the Committee on Data for Science and Technology
# (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
# Unit conversions
JOULE_PER_EV = 1.602176634e-19
# Avogadro's constant
AVOGADRO = 6.02214076e23
# Neutron mass in units of amu
NEUTRON_MASS = 1.00866491595
# Used in atomic_mass function as a cache
_ATOMIC_MASS: dict[str, float] = {}
# Used in half_life function as a cache
_HALF_LIFE: dict[str, float] = {}
_LOG_TWO = log(2.0)
[docs]
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
Atomic mass data comes from the `Atomic Mass Evaluation 2020
<https://doi.org/10.1088/1674-1137/abddaf>`_.
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
Returns
-------
float
Atomic mass of isotope in [amu]
"""
if not _ATOMIC_MASS:
# Load data from AME2020 file
mass_file = os.path.join(os.path.dirname(__file__), 'mass_1.mas20.txt')
with open(mass_file, 'r') as ame:
# Read lines in file starting at line 37
for line in itertools.islice(ame, 36, None):
name = f'{line[20:22].strip()}{int(line[16:19])}'
mass = float(line[106:109]) + 1e-6*float(
line[110:116] + '.' + line[117:123])
_ATOMIC_MASS[name.lower()] = mass
# For isotopes found in some libraries that represent all natural
# isotopes of their element (e.g. C0), calculate the atomic mass as
# the sum of the atomic mass times the natural abundance of the isotopes
# that make up the element.
for element in ['C', 'Zn', 'Pt', 'Os', 'Tl', 'V']:
isotope_zero = element.lower() + '0'
_ATOMIC_MASS[isotope_zero] = 0.
for iso, abundance in isotopes(element):
_ATOMIC_MASS[isotope_zero] += abundance * _ATOMIC_MASS[iso.lower()]
# Get rid of metastable information
if '_' in isotope:
isotope = isotope[:isotope.find('_')]
return _ATOMIC_MASS[isotope.lower()]
[docs]
def atomic_weight(element):
"""Return atomic weight of an element in atomic mass units.
Computes an average of the atomic mass of each of element's naturally
occurring isotopes weighted by their relative abundance.
Parameters
----------
element : str
Element symbol (e.g., 'H') or name (e.g., 'helium')
Returns
-------
float
Atomic weight of element in [amu]
"""
weight = 0.
for nuclide, abundance in isotopes(element):
weight += atomic_mass(nuclide) * abundance
if weight > 0.:
return weight
else:
raise ValueError(f"No naturally-occurring isotopes for element '{element}'.")
[docs]
def half_life(isotope):
"""Return half-life of isotope in seconds or None if isotope is stable
Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
<https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
.. versionadded:: 0.13.1
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
Returns
-------
float
Half-life of isotope in [s]
"""
global _HALF_LIFE
if not _HALF_LIFE:
# Load ENDF/B-VIII.0 data from JSON file
half_life_path = Path(__file__).with_name('half_life.json')
_HALF_LIFE = json.loads(half_life_path.read_text())
return _HALF_LIFE.get(isotope.lower())
[docs]
def decay_constant(isotope):
"""Return decay constant of isotope in [s^-1]
Decay constants are based on half-life values from the
:func:`~openmc.data.half_life` function. When the isotope is stable, a decay
constant of zero is returned.
.. versionadded:: 0.13.1
Parameters
----------
isotope : str
Name of isotope, e.g., 'Pu239'
Returns
-------
float
Decay constant of isotope in [s^-1]
See also
--------
openmc.data.half_life
"""
t = half_life(isotope)
return _LOG_TWO / t if t else 0.0
[docs]
def water_density(temperature, pressure=0.1013):
"""Return the density of liquid water at a given temperature and pressure.
The density is calculated from a polynomial fit using equations and values
from the 2012 version of the IAPWS-IF97 formulation. Only the equations
for region 1 are implemented here. Region 1 is limited to liquid water
below 100 [MPa] with a temperature above 273.15 [K], below 623.15 [K], and
below saturation.
Reference: International Association for the Properties of Water and Steam,
"Revised Release on the IAPWS Industrial Formulation 1997 for the
Thermodynamic Properties of Water and Steam", IAPWS R7-97(2012).
Parameters
----------
temperature : float
Water temperature in units of [K]
pressure : float
Water pressure in units of [MPa]
Returns
-------
float
Water density in units of [g/cm^3]
"""
# Make sure the temperature and pressure are inside the min/max region 1
# bounds. (Relax the 273.15 bound to 273 in case a user wants 0 deg C data
# but they only use 3 digits for their conversion to K.)
if pressure > 100.0:
warn("Results are not valid for pressures above 100 MPa.")
elif pressure < 0.0:
raise ValueError("Pressure must be positive.")
if temperature < 273:
warn("Results are not valid for temperatures below 273.15 K.")
elif temperature > 623.15:
warn("Results are not valid for temperatures above 623.15 K.")
elif temperature <= 0.0:
raise ValueError('Temperature must be positive.')
# IAPWS region 4 parameters
n4 = [0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2,
0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2,
-0.48232657361591e4, 0.40511340542057e6, -0.23855557567849,
0.65017534844798e3]
# Compute the saturation temperature at the given pressure.
beta = pressure**(0.25)
E = beta**2 + n4[2] * beta + n4[5]
F = n4[0] * beta**2 + n4[3] * beta + n4[6]
G = n4[1] * beta**2 + n4[4] * beta + n4[7]
D = 2.0 * G / (-F - sqrt(F**2 - 4 * E * G))
T_sat = 0.5 * (n4[9] + D
- sqrt((n4[9] + D)**2 - 4.0 * (n4[8] + n4[9] * D)))
# Make sure we aren't above saturation. (Relax this bound by .2 degrees
# for deg C to K conversions.)
if temperature > T_sat + 0.2:
warn("Results are not valid for temperatures above saturation "
"(above the boiling point).")
# IAPWS region 1 parameters
R_GAS_CONSTANT = 0.461526 # kJ / kg / K
ref_p = 16.53 # MPa
ref_T = 1386 # K
n1f = [0.14632971213167, -0.84548187169114, -0.37563603672040e1,
0.33855169168385e1, -0.95791963387872, 0.15772038513228,
-0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3,
-0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1,
-0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3,
-0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5,
-0.72694996297594e-15, -0.31679644845054e-4, -0.28270797985312e-5,
-0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6,
-0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8,
-0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19,
0.26335781662795e-22, -0.11947622640071e-22, 0.18228094581404e-23,
-0.93537087292458e-25]
I1f = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4,
4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32]
J1f = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4,
0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41]
# Nondimensionalize the pressure and temperature.
pi = pressure / ref_p
tau = ref_T / temperature
# Compute the derivative of gamma (dimensionless Gibbs free energy) with
# respect to pi.
gamma1_pi = 0.0
for n, I, J in zip(n1f, I1f, J1f):
gamma1_pi -= n * I * (7.1 - pi)**(I - 1) * (tau - 1.222)**J
# Compute the leading coefficient. This sets the units at
# 1 [MPa] * [kg K / kJ] * [1 / K]
# = 1e6 [N / m^2] * 1e-3 [kg K / N / m] * [1 / K]
# = 1e3 [kg / m^3]
# = 1 [g / cm^3]
coeff = pressure / R_GAS_CONSTANT / temperature
# Compute and return the density.
return coeff / pi / gamma1_pi
def _get_element_symbol(element: str) -> str:
if len(element) > 2:
symbol = ELEMENT_SYMBOL.get(element.lower())
if symbol is None:
raise ValueError(f'Element name "{element}" not recognized')
return symbol
else:
return element
[docs]
def isotopes(element: str) -> list[tuple[str, float]]:
"""Return naturally occurring isotopes and their abundances
.. versionadded:: 0.12.1
Parameters
----------
element : str
Element symbol (e.g., 'H') or name (e.g., 'helium')
Returns
-------
list
A list of tuples of (isotope, abundance)
Raises
------
ValueError
If the element name is not recognized
"""
element = _get_element_symbol(element)
# Get the nuclides present in nature
result = []
for kv in NATURAL_ABUNDANCE.items():
if re.match(r'{}\d+'.format(element), kv[0]):
result.append(kv)
return result