Merge pull request #1023 from amandalund/photon-new

Photon Transport

.gitignore

@@ -59,14 +59,12 @@ src/install_manifest.txt
 scripts/nndc
 scripts/nndc_hdf5
 scripts/wmp
-scripts/multipole_lib.tar.gz
-scripts/ENDF-B-VII.1-*.tar.gz
-scripts/JEFF32-ACE-*.tar.gz
 scripts/JEFF32-ACE-*.zip
 scripts/TSLs.tar.gz
 scripts/jeff-3.2
 scripts/jeff-3.2-hdf5
-scripts/*.tar.xz
+scripts/*.tar.*
+scripts/G4EMLOW*/
 
 # Images
 *.ppm

CMakeLists.txt

@@ -319,6 +319,8 @@ add_library(libopenmc SHARED
   src/output.F90
   src/particle_header.F90
   src/particle_restart.F90
+  src/photon_header.F90
+  src/photon_physics.F90
   src/physics_common.F90
   src/physics.F90
   src/physics_mg.F90
@@ -371,6 +373,7 @@ add_library(libopenmc SHARED
   src/tallies/tally_filter_mesh.F90
   src/tallies/tally_filter_meshsurface.F90
   src/tallies/tally_filter_mu.F90
+  src/tallies/tally_filter_particle.F90
   src/tallies/tally_filter_polar.F90
   src/tallies/tally_filter_sph_harm.F90
   src/tallies/tally_filter_sptl_legendre.F90

docs/source/io_formats/nuclear_data.rst

@@ -1,8 +1,8 @@
 .. _io_nuclear_data:
 
-========================
-Nuclear Data File Format
-========================
+=========================
+Nuclear Data File Formats
+=========================
 
 ---------------------
 Incident Neutron Data
@@ -10,7 +10,7 @@ Incident Neutron Data
 
 **/**
 
-:Attributes:
+:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
              - **version** (*int[2]*) -- Major and minor version of the data
 
 **/<nuclide name>/**
@@ -22,7 +22,9 @@ Incident Neutron Data
              - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses
              - **n_reaction** (*int*) -- Number of reactions
 
-:Datasets: - **energy** (*double[]*) -- Energy points at which cross sections are tabulated
+:Datasets:
+           - **energy** (*double[]*) -- Energies in [eV] at which cross sections
+             are tabulated
 
 **/<nuclide name>/kTs/**
 
@@ -31,7 +33,7 @@ temperature-dependent data set.  For example, the data set corresponding to
 300 Kelvin would be located at `300K`.
 
 :Datasets:
-           - **<TTT>K** (*double*) -- kT values (in eV) for each temperature
+           - **<TTT>K** (*double*) -- kT values in [eV] for each temperature
              TTT (in Kelvin)
 
 **/<nuclide name>/reactions/reaction_<mt>/**
@@ -113,6 +115,92 @@ temperature-dependent data set.  For example, the data set corresponding to
              :ref:`tabulated <1d_tabulated>`) -- The recoverable fission Q-value
              (Q_prompt + delayed neutrons + delayed photons + betas)
 
+--------------------
+Incident Photon Data
+--------------------
+
+**/**
+
+:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
+             - **version** (*int[2]*) -- Major and minor version of the data
+
+**/<element>/**
+
+:Attributes: - **Z** (*int*) -- Atomic number
+
+:Datasets:
+           - **energy** (*double[]*) -- Energies in [eV] at which cross sections
+             are tabulated
+
+**/<element>/bremsstrahlung/**
+
+:Datasets: - **electron_energy** (*double[]*) -- Incident electron energy in [eV]
+           - **photon_energy** (*double[]*) -- Outgoing photon energy as
+             fraction of incident electron energy
+           - **dcs** (*double[][]*) -- Bremsstrahlung differential cross section
+             at each incident energy in [mb/eV]
+
+**/<element>/coherent/**
+
+:Datasets: - **xs** (*double[]*) -- Coherent scattering cross section in [b]
+           - **integrated_scattering_factor** (:ref:`tabulated <1d_tabulated>`)
+             -- Integrated coherent scattering form factor
+           - **anomalous_real** (:ref:`tabulated <1d_tabulated>`) -- Real part
+             of the anomalous scattering factor
+           - **anomalous_imag** (:ref:`tabulated <1d_tabulated>`) -- Imaginary
+             part of the anomalous scattering factor
+
+**/<element>/compton_profiles/**
+
+:Datasets: - **binding_energy** (*double[]*) -- Binding energy for each subshell in [eV]
+           - **num_electrons** (*double[]*) -- Number of electrons in each subshell
+           - **pz** (*double[]*) -- Projection of the electron momentum on the
+             scattering vector in units of :math:`me^2 / \hbar` where :math:`m`
+             is the electron rest mass and :math:`e` is the electron charge
+           - **J** (*double[][]*) -- Compton profile for each subshell in units
+             of :math:`\hbar / (me^2)`
+
+**/<element>/incoherent/**
+
+:Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b]
+           - **scattering_factor** (:ref:`tabulated <1d_tabulated>`) --
+
+**/<element>/pair_production_electron/**
+
+:Datasets: - **xs** (*double[]*) -- Pair production (electron field) cross section in [b]
+
+**/<element>/pair_production_nuclear/**
+
+:Datasets: - **xs** (*double[]*) -- Pair production (nuclear field) cross section in [b]
+
+**/<element>/photoelectric/**
+
+:Datasets: - **xs** (*double[]*) -- Total photoionization cross section in [b]
+
+**/<element>/stopping_powers/**
+
+:Datasets: - **I** (*double*) -- Mean excitation energy in [eV]
+           - **energy** (*double[]*) -- Energies in [eV]
+           - **s_collision** (*double[]*) -- Collision stopping power in [eV-cm\ :sup:`2`\ /g]
+           - **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g]
+
+**/<element>/subshells/**
+
+:Attributes: - **designators** (*char[][]*) -- Designator for each shell, e.g. 'M2'
+
+**/<element>/subshells/<designator>/**
+
+:Attributes: - **binding_energy** (*double*) -- Binding energy of the subshell in [eV]
+             - **num_electrons** (*double*) -- Number of electrons in the subshell
+
+:Datasets: - **transitions** (*double[][]*) -- Atomic relaxation data
+           - **xs** (*double[]*) -- Photoionization cross section for subshell
+             in [b] tabulated against the main energy grid
+
+             :Attributes:
+                          - **threshold_idx** (*int*) -- Index on the energy
+                            grid of the reaction threshold
+
 -------------------------------
 Thermal Neutron Scattering Data
 -------------------------------

docs/source/io_formats/settings.rst

@@ -32,6 +32,19 @@ standard deviation.
 
   *Default*: false
 
+-------------------------------------
+``<create_fission_neutrons>`` Element
+-------------------------------------
+
+The ``<create_fission_neutrons>`` element indicates whether fission neutrons
+should be created or not.  If this element is set to "true", fission neutrons
+will be created; otherwise the fission is treated as capture and no fission
+neutron will be created. Note that this option is only applied to fixed source
+calculation. For eigenvalue calculation, fission will always be treated as real
+fission.
+
+  *Default*: true
+
 --------------------
 ``<cutoff>`` Element
 --------------------
@@ -55,31 +68,35 @@ you care. This element has the following attributes/sub-elements:
 
     *Default*: 1.0
 
-  :energy:
-    The energy under which particles will be killed.
+  :energy_neutron:
+    The energy under which neutrons will be killed.
 
     *Default*: 0.0
 
--------------------------
-``<energy_grid>`` Element
--------------------------
+  :energy_photon:
+    The energy under which photons will be killed.
 
-The ``<energy_grid>`` element determines the treatment of the energy grid during
-a simulation. The valid options are "nuclide", "logarithm", and
-"material-union". Setting this element to "nuclide" will cause OpenMC to use a
-nuclide's energy grid when determining what points to interpolate between for
-determining cross sections (i.e. non-unionized energy grid). Setting this
-element to "logarithm" causes OpenMC to use a logarithmic mapping technique
-described in LA-UR-14-24530_. Setting this element to "material-union" will
-cause OpenMC to create energy grids that are unionized material-by-material and
-use these grids when determining the energy-cross section pairs to interpolate
-cross section values between.
+    *Default*: 1000.0
 
-  *Default*: logarithm
+  :energy_electron:
+    The energy under which electrons will be killed.
 
-  .. note:: This element is not used in the multi-group :ref:`energy_mode`.
+    *Default*: 0.0
 
-.. _LA-UR-14-24530: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
+  :energy_positron:
+    The energy under which positrons will be killed.
+
+    *Default*: 0.0
+
+--------------------------------
+``<electron_treatment>`` Element
+--------------------------------
+
+When photon transport is enabled, the ``<electron_treatment>`` element tells
+OpenMC whether to deposit all energy from electrons locally (``led``) or create
+secondary bremsstrahlung photons (``ttb``).
+
+  *Default*: ttb
 
 .. _energy_mode:
 
@@ -153,8 +170,7 @@ the estimated eigenvalue. It has the following attributes/sub-elements:
 
     *Default*: None
 
-.. note:: See section on the :ref:`trigger` for more information.
-
+  .. note:: See section on the :ref:`trigger` for more information.
 
 ---------------------------
 ``<log_grid_bins>`` Element
@@ -169,6 +185,8 @@ based on the recommended value in LA-UR-14-24530_.
 
   .. note:: This element is not used in the multi-group :ref:`energy_mode`.
 
+.. _LA-UR-14-24530: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
+
 ---------------------------
 ``<max_order>`` Element
 ---------------------------
@@ -259,11 +277,21 @@ out the file and "false" will not.
 -----------------------
 
 This element indicates the number of neutrons to simulate per fission source
-iteration when a k-eigenvalue calculation is performed or the number of neutrons
-per batch for a fixed source simulation.
+iteration when a k-eigenvalue calculation is performed or the number of
+particles per batch for a fixed source simulation.
 
   *Default*: None
 
+------------------------------
+``<photon_transport>`` Element
+------------------------------
+
+The ``<photon_transport>`` element determines whether photon transport is
+enabled. This element has no attributes or sub-elements and can be set to
+either "false" or "true".
+
+  *Default*: false
+
 ---------------------
 ``<ptables>`` Element
 ---------------------
@@ -379,6 +407,11 @@ attributes/sub-elements:
 
     *Default*: 1.0
 
+  :particle:
+    The source particle type, either ``neutron`` or ``photon``.
+
+    *Default*: neutron
+
   :file:
     If this attribute is given, it indicates that the source is to be read from
     a binary source file whose path is given by the value of this element. Note,
@@ -812,20 +845,6 @@ and 10. The verbosity levels are defined as follows:
 
   *Default*: 7
 
--------------------------------------
-``<create_fission_neutrons>`` Element
--------------------------------------
-
-The ``<create_fission_neutrons>`` element indicates whether fission neutrons
-should be created or not.  If this element is set to "true", fission neutrons
-will be created; otherwise the fission is treated as capture and no fission
-neutron will be created. Note that this option is only applied to fixed source
-calculation. For eigenvalue calculation, fission will always be treated as real
-fission.
-
-  *Default*: true
-
-
 -------------------------
 ``<volume_calc>`` Element
 -------------------------

docs/source/io_formats/tallies.rst

@@ -125,8 +125,8 @@ attributes/sub-elements:
   :type:
     The type of the filter. Accepted options are "cell", "cellfrom", 
     "cellborn", "surface", "material", "universe", "energy", "energyout", "mu", 
-    "polar", "azimuthal", "mesh", "distribcell", "delayedgroup", and
-    "energyfunction".
+    "polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
+    "energyfunction", and "particle".
 
   :bins:
      A description of the bins for each type of filter can be found in
@@ -299,6 +299,9 @@ should be set to:
   ``energyfunction`` filters do not use the ``bins`` entry.  Instead
   they use ``energy`` and ``y``.
 
+:particle:
+  A list of integers indicating the type of particles to tally ('neutron' = 1,
+  'photon' = 2, 'electron' = 3, 'positron' = 4).
 
 ------------------
 ``<mesh>`` Element

docs/source/methods/index.rst

@@ -12,7 +12,8 @@ Theory and Methodology
     geometry
     cross_sections
     random_numbers
-    physics
+    neutron_physics
+    photon_physics
     tallies
     eigenvalue
     parallelization

docs/source/methods/neutron_physics.rst

@@ -1,8 +1,8 @@
-.. _methods_physics:
+.. _methods_neutron_physics:
 
-=======
-Physics
-=======
+===============
+Neutron Physics
+===============
 
 There are limited differences between physics treatments used in the
 continuous-energy and multi-group modes.  If distinctions are necessary, each

docs/source/methods/photon_physics.rst

@@ -0,0 +1,275 @@
+.. _methods_photon_physics:
+
+==============
+Photon Physics
+==============
+
+Photons, being neutral particles, behave much in the same manner as neutrons,
+traveling in straight lines and experiencing occasional collisions which change
+their energy and direction. Photons undergo four basic interactions as they pass
+through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering,
+photoelectric effect, and pair/triplet production. Photons with energy in the
+MeV range may also undergo photonuclear reactions with an atomic nucleus. In
+addition to these primary interaction mechanisms, all processes other than
+coherent scattering can result in the excitation/ionization of atoms. The
+de-excitation of these atoms can result in the emission of electrons and
+photons. Electrons themselves also can produce photons by means of
+bremsstrahlung radiation.
+
+-------------------
+Photon Interactions
+-------------------
+
+Coherent (Rayleigh) Scattering
+------------------------------
+
+The elastic scattering of a photon off a free charged particle is known as
+Thomson scattering. The differential cross section is independent of the energy
+of the incident photon. For scattering off a free electron, the differential
+cross section is
+
+.. math::
+    :label: thomson
+
+    \frac{d\sigma}{d\mu} = \pi r_e^2 ( 1 + \mu^2 )
+
+where :math:`\mu` is the cosine of the scattering angle and :math:`r_e` is the
+classical electron radius. Thomson scattering can generally occur when the
+photon energy is much less than the rest mass energy of the particle.
+
+In practice, most elastic scattering of photons off electrons happens not with
+free electrons but those bound in atoms. This process is known as Rayleigh
+scattering. The radiation scattered off of individual bound electrons combines
+coherently, and thus Rayleigh scattering is also known as coherent
+scattering. Even though conceptually we think of the photon interacting with a
+single electron, because the wave functions combine constructively it is really
+as though the photon is interacting with the entire atom.
+
+The differential cross section for Rayleigh scattering is given by
+
+.. math::
+    :label: coherent-xs
+
+    \frac{d\sigma(E,E',\mu)}{d\mu} &= \pi r_e^2 ( 1 + \mu^2 )~\left| F(x,Z)
+                                      + F' + iF'' \right|^2 \\
+                                   &= \pi r_e^2 ( 1 + \mu^2 ) \left [ ( F(x,Z)
+                                      + F'(E) )^2 + F''(E)^2 \right ]
+
+where :math:`F(x,Z)` is a form factor as a function of the momentum transfer
+:math:`x` and the atomic number :math:`Z` and the term :math:`F' + iF''`
+accounts for `anomalous scattering`_ which can occur near absorption edges. In
+a Monte Carlo simulation, when coherent scattering occurs, we only need to
+sample the scattering angle using the differential cross section in
+:eq:`coherent-xs` since the energy of the photon does not change. In OpenMC,
+anomalous scattering is ignored such that the differential cross section
+becomes
+
+.. math::
+    :label: coherent-xs-openmc
+
+    \frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_e^2 ( 1 + \mu^2 ) F(x, Z)^2
+
+To construct a proper probability density, we need to normalize the
+differential cross section in :eq:`coherent-xs-openmc` by the integrated
+coherent scattering cross section:
+
+.. math::
+    :label: coherent-pdf-1
+
+    p(\mu) d\mu = \frac{\pi r_e^2}{\sigma(E)} ( 1 + \mu^2 ) F(x, Z)^2 d\mu.
+
+Since the form factor is given in terms of the momentum transfer, it is more
+convenient to change variables of the probability density to :math:`x^2`. The
+momentum transfer is traditionally expressed as
+
+.. math::
+    :label: momentum-transfer
+
+    x = \kappa \alpha \sqrt{1 - \mu}
+
+where :math:`\alpha` is the ratio of the photon energy to the electron rest
+mass, and the coefficient :math:`\kappa` can be shown to be
+
+.. math::
+    :label: kappa
+
+    \kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329,
+
+where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light
+in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`,
+we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 =
+-1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is
+
+.. math::
+    :label: coherent-pdf-x2
+
+    p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi
+    r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left (
+    \frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2
+
+where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for
+:math:`\mu=-1`,
+
+.. math::
+    :label: xmax
+
+    \bar{x} = \kappa \alpha \sqrt{2} = \frac{m_e c^2}{hc} \alpha,
+
+and :math:`A(x^2, Z)` is the integral of the square of the form factor:
+
+.. math::
+    :label: coherent-int-ff
+
+    A(x^2, Z) = \int_0^{x^2} F(x,Z)^2 dx^2.
+
+As you see, we have multiplied and divided the probability density by the
+integral of the squared form factor so that the density in :eq:`coherent-pdf-x2`
+is expressed as the product of two separate densities in parentheses. In OpenMC,
+a table of :math:`A(x^2, Z)` versus :math:`x^2` is pre-generated and used at
+run-time to do a table search on the cumulative distribution function:
+
+.. math::
+    :label: coherent-form-factor-cdf
+
+    \frac{\int_0^{x^2} F(x,Z)^2 dx^2}{\int_0^{\bar{x}^2} F(x,Z)^2 dx^2}
+
+Once a trial :math:`x^2` value has been selected, we can calculate :math:`\mu`
+and perform rejection sampling using the Thomson scattering differential cross
+section. The complete algorithm is as follows:
+
+1. Determine :math:`\bar{x}^2` using :eq:`xmax`.
+
+2. Determine :math:`A_{max} = A(\bar{x}^2, Z)` using the pre-generated
+   tabulated data.
+
+3. Sample the cumulative density by calculating :math:`A' = \xi_1 A_{max}` where
+   :math:`\xi_1` is a uniformly distributed random number.
+
+4. Perform a binary search to determine the value of :math:`x^2` which satisfies
+   :math:`A(x^2, Z) = A'`.
+
+5. By combining :eq:`momentum-transfer` and :eq:`xmax`, calculate :math:`\mu =
+   1 - 2x^2/\bar{x}^2`.
+
+6. If :math:`\xi_2 < (1 + \mu^2)/2`, accept :math:`\mu`. Otherwise, repeat the
+   sampling at step 3.
+
+Incoherent (Compton) Scattering
+-------------------------------
+
+Before we noted that the Thomson cross section gives the behavior for photons
+scattering off of free electrons valid at low energies. The formula for photon
+scattering off of free electrons that is valid for all energies can be found
+using quantum electrodynamics and is known as the Klein-Nishina_ formula after
+the two authors who discovered it:
+
+.. math::
+    :label: klein-nishina
+
+    \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right
+    )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1
+    \right ]
+
+where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and
+exiting photon energies to the electron rest mass energy equivalent (0.511 MeV),
+respectively. Although it appears that the outgoing energy and angle are
+separate, there is actually a one-to-one relationship between them such that
+only one needs to be sampled:
+
+.. math::
+    :label: compton-energy-angle
+
+    \alpha' = \frac{\alpha}{1 + \alpha(1 - \mu)}.
+
+Note that when :math:`\alpha'/\alpha` goes to one, i.e., scattering is elastic,
+the Klein-Nishina cross section becomes identical to the Thomson cross
+section. In general though, the scattering is inelastic and is known as Compton
+scattering. When a photon interacts with a bound electron in an atom, the
+Klein-Nishina formula must be modified to account for the binding effects. As in
+the case of coherent scattering, this is done by means of a form factor. The
+differential cross section for incoherent scattering is given by
+
+.. math::
+    :label: incoherent-xs
+
+    \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left (
+    \frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} +
+    \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z)
+
+where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first
+sample the Klein-Nishina cross section and then perform rejection sampling on
+the form factor. As in other codes, `Kahn's rejection method`_ is used for
+:math:`\alpha < 3` and a direct method by Koblinger_ is used for :math:`\alpha
+\ge 3`. The complete algorithm is as follows:
+
+1. If :math:`\alpha < 3`, sample :math:`\mu` from the Klein-Nishina cross
+   section using Kahn's rejection method. Otherwise, use Koblinger's direct
+   method.
+
+2. Calculate :math:`x` and :math:`\bar{x}` using :eq:`momentum-transfer` and
+   :eq:`xmax`, respectively.
+
+3. If :math:`\xi < S(x, Z)/S(\bar{x}, Z)`, accept :math:`\mu`. Otherwise repeat
+   from step 1.
+
+Doppler Energy Broadening
++++++++++++++++++++++++++
+
+LA-UR-04-0487_ and LA-UR-04-0488_
+
+Compton Electrons
++++++++++++++++++
+
+
+Photoelectric Effect
+--------------------
+
+
+Pair Production
+---------------
+
+
+-------------------
+Secondary Processes
+-------------------
+
+New photons may be produced in secondary processes related to the main photon
+interactions discussed above. A Compton-scattered photon transfers a portion of
+its energy to the kinetic energy of the recoil electron, which in turn may lose
+the energy as bremsstrahlung radiation. The vacancy left in the shell by the
+ejected electron is filled through atomic relaxation, creating a shower of
+electrons and fluorescence photons. Similarly, the vacancy left by the electron
+emitted in the photoelectric effect is filled through atomic relaxation. Pair
+production generates an electron and a positron, both of which can emit
+bremsstrahlung radiation before the positron eventually collides with an
+electron, resulting in annihilation of the pair and the creation of two
+additional photons.
+
+Atomic Relaxation
+-----------------
+
+
+Electron-Positron Annihilation
+------------------------------
+
+
+Bremsstrahlung
+--------------
+
+.. _ttb:
+
+Thick-Target Bremsstrahlung Approximation
++++++++++++++++++++++++++++++++++++++++++
+
+
+.. _Koblinger: https://doi.org/10.13182/NSE75-A26663
+
+.. _anomalous scattering: http://pd.chem.ucl.ac.uk/pdnn/diff1/anomscat.htm
+
+.. _Kahn's rejection method: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/aecu-3259_kahn.pdf
+
+.. _Klein-Nishina: https://en.wikipedia.org/wiki/Klein%E2%80%93Nishina_formula
+
+.. _LA-UR-04-0487: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0487.pdf
+
+.. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf

docs/source/pythonapi/base.rst

@@ -122,6 +122,7 @@ Constructing Tallies
    openmc.SpatialLegendreFilter
    openmc.SphericalHarmonicsFilter
    openmc.ZernikeFilter
+   openmc.ParticleFilter
    openmc.Mesh
    openmc.Trigger
    openmc.TallyDerivative

docs/source/pythonapi/data.rst

@@ -19,6 +19,9 @@ Core Classes
     openmc.data.CoherentElastic
     openmc.data.FissionEnergyRelease
     openmc.data.DataLibrary
+    openmc.data.IncidentPhoton
+    openmc.data.PhotonReaction
+    openmc.data.AtomicRelaxation
     openmc.data.Decay
     openmc.data.FissionProductYields
     openmc.data.WindowedMultipole

docs/source/usersguide/beginners.rst

@@ -8,19 +8,19 @@ A Beginner's Guide to OpenMC
 What does OpenMC do?
 --------------------
 
-In a nutshell, OpenMC simulates neutral particles (presently only neutrons)
-moving stochastically through an arbitrarily defined model that represents an
-real-world experimental setup. The experiment could be as simple as a sphere of
-metal or as complicated as a full-scale `nuclear reactor`_. This is what's known
-as `Monte Carlo`_ simulation. In the case of a nuclear reactor model, neutrons
-are especially important because they are the particles that induce `fission`_
-in isotopes of uranium and other elements. Knowing the behavior of neutrons
-allows one to determine how often and where fission occurs. The amount of energy
-released is then directly proportional to the fission reaction rate since most
-heat is produced by fission. By simulating many neutrons (millions or billions),
-it is possible to determine the average behavior of these neutrons (or the
-behavior of the energy produced, or any other quantity one is interested in)
-very accurately.
+In a nutshell, OpenMC simulates neutral particles (presently neutrons and
+photons) moving stochastically through an arbitrarily defined model that
+represents an real-world experimental setup. The experiment could be as simple
+as a sphere of metal or as complicated as a full-scale `nuclear reactor`_. This
+is what's known as `Monte Carlo`_ simulation. In the case of a nuclear reactor
+model, neutrons are especially important because they are the particles that
+induce `fission`_ in isotopes of uranium and other elements. Knowing the
+behavior of neutrons allows one to determine how often and where fission
+occurs. The amount of energy released is then directly proportional to the
+fission reaction rate since most heat is produced by fission. By simulating
+many neutrons (millions or billions), it is possible to determine the average
+behavior of these neutrons (or the behavior of the energy produced, or any
+other quantity one is interested in) very accurately.
 
 Using Monte Carlo methods to determine the average behavior of various physical
 quantities in a system is quite different from other means of solving the same

docs/source/usersguide/cross_sections.rst

@@ -258,6 +258,49 @@ method using :attr:`Settings.resonance_scattering`.
           running the :meth:`IncidentNeutron.add_elastic_0K_from_endf` method
           may take several minutes to complete.
 
+Photon Cross Sections
+---------------------
+
+Photon interaction data is needed to run OpenMC with photon transport enabled.
+Some of this data, namely bremsstrahlung cross sections from `Seltzer and
+Berger`_, stopping powers from the `NIST ESTAR database`_, and Compton profiles
+calculated by `Biggs et al.`_ and available in the Geant4 G4EMLOW data file, is
+distributed with OpenMC. The rest is available from the NNDC, which provides
+ENDF data from the photo-atomic and atomic relaxation sublibraries of the
+ENDF/B-VII.1 library. By default, the :ref:`scripts_nndc` script will download
+the ENDF data in addition to the neutron and thermal scattering data, extract
+it, combine it with the data from other sources, and convert it to an HDF5
+library. Alternatively, the :ref:`scripts_photon` script can be used to
+download the photon data on its own and create the HDF5 library:
+
+.. code-block:: sh
+
+    openmc-get-photon-data
+
+As with neutrons and thermal scattering, it is possible to use the Python API
+directly to convert photon interaction data from an ENDF or ACE file to an HDF5
+file. The :class:`openmc.data.IncidentPhoton` class contains an
+:meth:`IncidentPhoton.from_ace` method that will generate photon data from an
+ACE table and an :meth:`IncidentPhoton.export_to_hdf5` method that writes the
+data to an HDF5 file:
+
+::
+
+  u = openmc.data.IncidentPhoton.from_ace('92000.12p')
+  u.export_to_hdf5('U.h5')
+
+Similarly, the :meth:`IncidentPhoton.from_endf` method can be used to read
+photon data from an ENDF file. In the case, both the photo-atomic and atomic
+relaxation sublibrary files are required:
+
+::
+
+  u = openmc.data.IncidentPhoton.from_endf('photoat-092_U_000.endf',
+                                           'atom-092_U_000.endf')
+
+Once the HDF5 files have been generated, a library can be created using the
+:class:`DataLibrary` class as described in :ref:`create_xs_library`.
+
 -----------------------
 Windowed Multipole Data
 -----------------------
@@ -291,3 +334,6 @@ For an example of how to create a multi-group library, see
 .. _MCNP: http://mcnp.lanl.gov
 .. _Serpent: http://montecarlo.vtt.fi
 .. _TENDL: https://tendl.web.psi.ch/tendl_2015/tendl2015.html
+.. _Seltzer and Berger: https://www.sciencedirect.com/science/article/pii/0092640X86900148?via%3Dihub
+.. _NIST ESTAR database: https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
+.. _Biggs et al.: https://www.sciencedirect.com/science/article/pii/0092640X75900303

docs/source/usersguide/scripts.rst

@@ -166,12 +166,57 @@ ENDF/B-VII.1. It has the following optional arguments:
 ``openmc-get-nndc-data``
 ------------------------
 
-This script downloads `ENDF/B-VII.1 ACE data
-<http://www.nndc.bnl.gov/endf/b7.1/acefiles.html>`_ from NNDC and converts it to
-an HDF5 library for use with OpenMC. This script has the following optional
-arguments:
+This script downloads `ENDF/B-VII.1
+<http://www.nndc.bnl.gov/endf/b7.1/acefiles.html>`_ incident neutron ACE data
+and incident photon ENDF data from NNDC and converts it to an HDF5 library for
+use with OpenMC. This script has the following optional arguments:
 
--b, --batch     Suppress standard in
+-b, --batch
+    Suppress standard in
+
+-n, --neutron_only
+    Whether to exclude photon interaction/atomic data
+
+.. _scripts_photon:
+
+--------------------------
+``openmc-get-photon-data``
+--------------------------
+
+This script downloads `ENDF data <http://www.nndc.bnl.gov/endf/b7.1/zips/>`_
+from NNDC for photo-atomic and atomic relaxation sublibraries and converts it
+to an HDF5 library for use with photon transport in OpenMC. This script has the
+following optional arguments:
+
+-b, --batch
+    Suppress standard in
+
+-c, --cross-sections
+    cross_sections.xml file to append libraries to
+
+.. _scripts_compton:
+
+-----------------------
+``openmc-make-compton``
+-----------------------
+
+This script generates an HDF5 file called ``compton_profiles.h5`` that contains
+Compton profile data using an existing data library from `Geant4
+<http://geant4.cern.ch/>`_. Note that OpenMC includes this data file by default
+so it should not be necessary in practice to generate it yourself.
+
+.. _scripts_stopping:
+
+-------------------------------
+``openmc-make-stopping-powers``
+-------------------------------
+
+This script generates an HDF5 file called ``stopping_power.h5`` that contains
+radiative and collision stopping powers and mean excitation energy pulled from
+the `NIST ESTAR database
+<https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html>`_. Note that OpenMC
+includes this data file by default so it should not be necessary in practice to
+generate it yourself.
 
 .. _scripts_plot:
 

docs/source/usersguide/settings.rst

@@ -151,10 +151,19 @@ time and another that should be sampled 30% of the time::
 
   settings.source = [src1, src2]
 
+Finally, the :attr:`Source.particle` attribute can be used to indicate the
+source should be composed of particles other than neutrons. For example, the
+following would generate a photon source::
+
+  source = openmc.Source()
+  source.particle = 'photon'
+  ...
+
+  settings.source = source
+
 For a full list of all classes related to statistical distributions, see
 :ref:`pythonapi_stats`.
 
-
 ---------------
 Shannon Entropy
 ---------------
@@ -190,6 +199,52 @@ property::
 
   settings.entropy_mesh = m
 
+----------------
+Photon Transport
+----------------
+
+In addition to neutrons, OpenMC is also capable of simulating the passage of
+photons through matter. This allows the modeling of photon production from
+neutrons as well as pure photon calculations. The
+:attr:`Settings.photon_transport` attribute can be used to enable photon
+transport::
+
+  settings.photon_transport = True
+
+The way in which OpenMC handles secondary charged particles can be specified
+with the :attr:`Settings.electron_treatment` attribute. By default, the
+:ref:`thick-target bremsstrahlung <ttb>` (TTB) approximation is used to generate
+bremsstrahlung radiation emitted by electrons and positrons created in photon
+interactions. To neglect secondary bremsstrahlung photons and instead deposit
+all energy from electrons locally, the local energy deposition option can be
+selected::
+
+  settings.electron_treatment = 'led'
+
+.. warning::
+   Currently, collision stopping powers used in the TTB approximation come from
+   the `NIST ESTAR database`_, which provides data for each element calculated
+   using by default the material density at standard temperature and pressure.
+   In OpenMC, stopping powers for compounds are calculated from this elemental
+   data using Bragg's additivity rule. However, this is not a good
+   approximation --- the collision stopping power is a function of certain
+   quantities, such as the mean excitation energy and particularly the density
+   effect correction, that depend on material properties. Data for constituent
+   elements in a compound cannot simply be summed together, but rather these
+   quantities should be calculated for the material. This treatment will be
+   especially poor when the density of a material is different from the
+   densities used in the NIST data.
+
+.. note::
+   Some features related to photon transport are not currently implemented,
+   including:
+
+     * Tallying photon energy deposition.
+     * Properly accounting for energy deposition in coupled n-p calculations.
+     * Generating a photon source from a neutron calculation that can be used
+       for a later fixed source photon calculation.
+     * Photoneutron reactions.
+
 --------------------------
 Generation of Output Files
 --------------------------
@@ -224,3 +279,5 @@ As an example, to write a statepoint file every five batches::
 
   settings.batches = n
   settings.statepoint = {'batches': range(5, n + 5, 5)}
+
+.. _NIST ESTAR database: https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html

openmc/data/__init__.py

@@ -9,12 +9,13 @@ WMP_VERSION = 'v0.2'
 
 from .data import *
 from .neutron import *
+from .photon import *
 from .decay import *
 from .reaction import *
-from .ace import *
+from . import ace
 from .angle_distribution import *
 from .function import *
-from .endf import *
+from . import endf
 from .energy_distribution import *
 from .product import *
 from .angle_energy import *

openmc/data/ace.py

@@ -16,13 +16,82 @@ generates ACE-format cross sections.
 """
 
 from os import SEEK_CUR
+from pathlib import PurePath
 import struct
 import sys
 
 import numpy as np
 
 from openmc.mixin import EqualityMixin
-from openmc.data.endf import _ENDF_FLOAT_RE
+import openmc.checkvalue as cv
+from .data import ATOMIC_SYMBOL, gnd_name
+from .endf import ENDF_FLOAT_RE
+
+
+def get_metadata(zaid, metastable_scheme='nndc'):
+    """Return basic identifying data for a nuclide with a given ZAID.
+
+    Parameters
+    ----------
+    zaid : int
+        ZAID (1000*Z + A) obtained from a library
+    metastable_scheme : {'nndc', 'mcnp'}
+        Determine how ZAID identifiers are to be interpreted in the case of
+        a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
+        encode metastable information, different conventions are used among
+        different libraries. In MCNP libraries, the convention is to add 400
+        for a metastable nuclide except for Am242m, for which 95242 is
+        metastable and 95642 (or 1095242 in newer libraries) is the ground
+        state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
+
+    Returns
+    -------
+    name : str
+        Name of the table
+    element : str
+        The atomic symbol of the isotope in the table; e.g., Zr.
+    Z : int
+        Number of protons in the nucleus
+    mass_number : int
+        Number of nucleons in the nucleus
+    metastable : int
+        Metastable state of the nucleus. A value of zero indicates ground state.
+
+    """
+
+    cv.check_type('zaid', zaid, int)
+    cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp'])
+
+    Z = zaid // 1000
+    mass_number = zaid % 1000
+
+    if metastable_scheme == 'mcnp':
+        if zaid > 1000000:
+            # New SZA format
+            Z = Z % 1000
+            if zaid == 1095242:
+                metastable = 0
+            else:
+                metastable = zaid // 1000000
+        else:
+            if zaid == 95242:
+                metastable = 1
+            elif zaid == 95642:
+                metastable = 0
+            else:
+                metastable = 1 if mass_number > 300 else 0
+    elif metastable_scheme == 'nndc':
+        metastable = 1 if mass_number > 300 else 0
+
+    while mass_number > 3 * Z:
+        mass_number -= 100
+
+    # Determine name
+    element = ATOMIC_SYMBOL[Z]
+    name = gnd_name(Z, mass_number, metastable)
+
+    return (name, element, Z, mass_number, metastable)
+
 
 def ascii_to_binary(ascii_file, binary_file):
     """Convert an ACE file in ASCII format (type 1) to binary format (type 2).
@@ -160,7 +229,7 @@ class Library(EqualityMixin):
 
         # Determine whether file is ASCII or binary
         try:
-            fh = open(filename, 'rb')
+            fh = open(str(filename), 'rb')
             # Grab 10 lines of the library
             sb = b''.join([fh.readline() for i in range(10)])
 
@@ -349,7 +418,7 @@ class Library(EqualityMixin):
             # after it). If it's too short, then we apply the ENDF float regular
             # expression. We don't do this by default because it's expensive!
             if xss.size != nxs[1] + 1:
-                datastr = _ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
+                datastr = ENDF_FLOAT_RE.sub(r'\1e\2', datastr)
                 xss = np.fromstring(datastr, sep=' ')
                 assert xss.size == nxs[1] + 1
 

openmc/data/endf.py

@@ -22,10 +22,31 @@ from .function import Tabulated1D, INTERPOLATION_SCHEME
 from openmc.stats.univariate import Uniform, Tabular, Legendre
 
 
-LIBRARIES = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF',
-             4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL',
-             31: 'INDL/V', 32: 'INDL/A', 33: 'FENDL', 34: 'IRDF',
-             35: 'BROND', 36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'}
+_LIBRARY = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF',
+            4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL',
+            17: 'TENDL', 18: 'ROSFOND', 21: 'SG-21', 31: 'INDL/V',
+            32: 'INDL/A', 33: 'FENDL', 34: 'IRDF', 35: 'BROND',
+            36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'}
+
+_SUBLIBRARY = {
+    0: 'Photo-nuclear data',
+    1: 'Photo-induced fission product yields',
+    3: 'Photo-atomic data',
+    4: 'Radioactive decay data',
+    5: 'Spontaneous fission product yields',
+    6: 'Atomic relaxation data',
+    10: 'Incident-neutron data',
+    11: 'Neutron-induced fission product yields',
+    12: 'Thermal neutron scattering data',
+    19: 'Neutron standards',
+    113: 'Electro-atomic data',
+    10010: 'Incident-proton data',
+    10011: 'Proton-induced fission product yields',
+    10020: 'Incident-deuteron data',
+    10030: 'Incident-triton data',
+    20030: 'Incident-helion (3He) data',
+    20040: 'Incident-alpha data'
+}
 
 SUM_RULES = {1: [2, 3],
              3: [4, 5, 11, 16, 17, 22, 23, 24, 25, 27, 28, 29, 30, 32, 33, 34, 35,
@@ -45,7 +66,7 @@ SUM_RULES = {1: [2, 3],
              106: list(range(750, 800)),
              107: list(range(800, 850))}
 
-_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
+ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
 
 
 def float_endf(s):
@@ -68,7 +89,7 @@ def float_endf(s):
         The number
 
     """
-    return float(_ENDF_FLOAT_RE.sub(r'\1e\2', s))
+    return float(ENDF_FLOAT_RE.sub(r'\1e\2', s))
 
 
 def _int_endf(s):
@@ -410,6 +431,14 @@ class Evaluation(object):
 
         self._read_header()
 
+    def __repr__(self):
+        if 'zsymam' in self.target:
+            name = self.target['zsymam'].replace(' ', '')
+        else:
+            name = 'Unknown'
+        return '<{} for {} {}>'.format(self.info['sublibrary'], name,
+                                       self.info['library'])
+
     def _read_header(self):
         file_obj = io.StringIO(self.section[1, 451])
 
@@ -422,8 +451,7 @@ class Evaluation(object):
         self._LRP = items[2]
         self.target['fissionable'] = (items[3] == 1)
         try:
-            global LIBRARIES
-            library = LIBRARIES[items[4]]
+            library = _LIBRARY[items[4]]
         except KeyError:
             library = 'Unknown'
         self.info['modification'] = items[5]
@@ -446,7 +474,7 @@ class Evaluation(object):
         self.projectile['mass'] = items[0]
         self.info['energy_max'] = items[1]
         library_release = items[2]
-        self.info['sublibrary'] = items[4]
+        self.info['sublibrary'] = _SUBLIBRARY[items[4]]
         library_version = items[5]
         self.info['library'] = (library, library_version, library_release)
 

openmc/data/library.py

@@ -52,14 +52,17 @@ class DataLibrary(EqualityMixin):
             Path to the file to be registered.
 
         """
-        h5file = h5py.File(filename, 'r')
+        with h5py.File(filename, 'r') as h5file:
 
-        materials = []
-        filetype = 'neutron'
-        for name in h5file:
-            if name.startswith('c_'):
-                filetype = 'thermal'
-            materials.append(name)
+            materials = []
+            if 'filetype' in h5file.attrs:
+                filetype = h5file.attrs['filetype'].decode().lstrip('data_')
+            else:
+                filetype = 'neutron'
+            for name in h5file:
+                if name.startswith('c_'):
+                    filetype = 'thermal'
+                materials.append(name)
 
         library = {'path': filename, 'type': filetype, 'materials': materials}
         self.libraries.append(library)
@@ -81,10 +84,9 @@ class DataLibrary(EqualityMixin):
         if common_dir == '':
             common_dir = '.'
 
-        directory = os.path.relpath(common_dir, os.path.dirname(path))
-        if directory != '.':
+        if os.path.relpath(common_dir, os.path.dirname(path)) != '.':
             dir_element = ET.SubElement(root, "directory")
-            dir_element.text = directory
+            dir_element.text = os.path.realpath(common_dir)
 
         for library in self.libraries:
             lib_element = ET.SubElement(root, "library")

openmc/data/neutron.py

@@ -14,8 +14,8 @@ import numpy as np
 import h5py
 
 from . import HDF5_VERSION, HDF5_VERSION_MAJOR
-from .ace import Library, Table, get_table
-from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnd_name
+from .ace import Library, Table, get_table, get_metadata
+from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
 from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record
 from .fission_energy import FissionEnergyRelease
 from .function import Tabulated1D, Sum, ResonancesWithBackground
@@ -34,88 +34,26 @@ from openmc.mixin import EqualityMixin
 _RESONANCE_ENERGY_GRID = np.logspace(-3, 3, 61)
 
 
-def _get_metadata(zaid, metastable_scheme='nndc'):
-    """Return basic identifying data for a nuclide with a given ZAID.
-
-    Parameters
-    ----------
-    zaid : int
-        ZAID (1000*Z + A) obtained from a library
-    metastable_scheme : {'nndc', 'mcnp'}
-        Determine how ZAID identifiers are to be interpreted in the case of
-        a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
-        encode metastable information, different conventions are used among
-        different libraries. In MCNP libraries, the convention is to add 400
-        for a metastable nuclide except for Am242m, for which 95242 is
-        metastable and 95642 (or 1095242 in newer libraries) is the ground
-        state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
-
-    Returns
-    -------
-    name : str
-        Name of the table
-    element : str
-        The atomic symbol of the isotope in the table; e.g., Zr.
-    Z : int
-        Number of protons in the nucleus
-    mass_number : int
-        Number of nucleons in the nucleus
-    metastable : int
-        Metastable state of the nucleus. A value of zero indicates ground state.
-
-    """
-
-    cv.check_type('zaid', zaid, int)
-    cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp'])
-
-    Z = zaid // 1000
-    mass_number = zaid % 1000
-
-    if metastable_scheme == 'mcnp':
-        if zaid > 1000000:
-            # New SZA format
-            Z = Z % 1000
-            if zaid == 1095242:
-                metastable = 0
-            else:
-                metastable = zaid // 1000000
-        else:
-            if zaid == 95242:
-                metastable = 1
-            elif zaid == 95642:
-                metastable = 0
-            else:
-                metastable = 1 if mass_number > 300 else 0
-    elif metastable_scheme == 'nndc':
-        metastable = 1 if mass_number > 300 else 0
-
-    while mass_number > 3 * Z:
-        mass_number -= 100
-
-    # Determine name
-    element = ATOMIC_SYMBOL[Z]
-    name = gnd_name(Z, mass_number, metastable)
-
-    return (name, element, Z, mass_number, metastable)
-
-
 class IncidentNeutron(EqualityMixin):
     """Continuous-energy neutron interaction data.
 
-    Instances of this class are not normally instantiated by the user but rather
-    created using the factory methods :meth:`IncidentNeutron.from_hdf5` and
-    :meth:`IncidentNeutron.from_ace`.
+    This class stores data derived from an ENDF-6 format neutron interaction
+    sublibrary. Instances of this class are not normally instantiated by the
+    user but rather created using the factory methods
+    :meth:`IncidentNeutron.from_hdf5`, :meth:`IncidentNeutron.from_ace`, and
+    :meth:`IncidentNeutron.from_endf`.
 
     Parameters
     ----------
     name : str
         Name of the nuclide using the GND naming convention
     atomic_number : int
-        Number of protons in the nucleus
+        Number of protons in the target nucleus
     mass_number : int
-        Number of nucleons in the nucleus
+        Number of nucleons in the target nucleus
     metastable : int
-        Metastable state of the nucleus. A value of zero indicates ground state.
+        Metastable state of the target nucleus. A value of zero indicates ground
+        state.
     atomic_weight_ratio : float
         Atomic mass ratio of the target nuclide.
     kTs : Iterable of float
@@ -125,22 +63,19 @@ class IncidentNeutron(EqualityMixin):
     Attributes
     ----------
     atomic_number : int
-        Number of protons in the nucleus
+        Number of protons in the target nucleus
     atomic_symbol : str
         Atomic symbol of the nuclide, e.g., 'Zr'
     atomic_weight_ratio : float
         Atomic weight ratio of the target nuclide.
-    energy : dict of numpy.ndarray
-        The energy values (eV) at which reaction cross-sections are tabulated.
-        They keys of the dict are the temperature string ('294K') for each
-        set of energies
     fission_energy : None or openmc.data.FissionEnergyRelease
         The energy released by fission, tabulated by component (e.g. prompt
         neutrons or beta particles) and dependent on incident neutron energy
     mass_number : int
-        Number of nucleons in the nucleus
+        Number of nucleons in the target nucleus
     metastable : int
-        Metastable state of the nucleus. A value of zero indicates ground state.
+        Metastable state of the target nucleus. A value of zero indicates ground
+        state.
     name : str
         Name of the nuclide using the GND naming convention
     reactions : collections.OrderedDict
@@ -512,6 +447,7 @@ class IncidentNeutron(EqualityMixin):
 
         # Open file and write version
         f = h5py.File(path, mode, libver=libver)
+        f.attrs['filetype'] = np.string_('data_neutron')
         f.attrs['version'] = np.array(HDF5_VERSION)
 
         # Write basic data
@@ -686,7 +622,7 @@ class IncidentNeutron(EqualityMixin):
         # If mass number hasn't been specified, make an educated guess
         zaid, xs = ace.name.split('.')
         name, element, Z, mass_number, metastable = \
-            _get_metadata(int(zaid), metastable_scheme)
+            get_metadata(int(zaid), metastable_scheme)
 
         # Assign temperature to the running list
         kTs = [ace.temperature*EV_PER_MEV]

openmc/data/reaction.py

@@ -1093,7 +1093,7 @@ class Reaction(EqualityMixin):
         ev : openmc.data.endf.Evaluation
             ENDF evaluation
         mt : int
-            The MT value of the reaction to get angular distributions for
+            The MT value of the reaction to get data for
 
         Returns
         -------

openmc/data/thermal.py

@@ -279,6 +279,7 @@ class ThermalScattering(EqualityMixin):
         """
         # Open file and write version
         f = h5py.File(path, mode, libver=libver)
+        f.attrs['filetype'] = np.string_('data_thermal')
         f.attrs['version'] = np.array(HDF5_VERSION)
 
         # Write basic data

openmc/filter.py

@@ -22,7 +22,7 @@ _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike'
+    'sphericalharmonics', 'zernike', 'particle'
 )
 
 _CURRENT_NAMES = (
@@ -31,6 +31,7 @@ _CURRENT_NAMES = (
     'z-min out', 'z-min in', 'z-max out', 'z-max in'
 )
 
+_PARTICLE_IDS = {'neutron': 1, 'photon': 2, 'electron': 3, 'positron': 4}
 
 class FilterMeta(ABCMeta):
     def __new__(cls, name, bases, namespace, **kwargs):
@@ -540,6 +541,46 @@ class SurfaceFilter(WithIDFilter):
     expected_type = Surface
 
 
+class ParticleFilter(Filter):
+    """Bins tally events based on the Particle type.
+
+    Parameters
+    ----------
+    bins : str, int, or iterable of Integral
+        The Particles to tally. Either str with particle type or their
+        ID numbers can be used ('neutron' = 1, 'photon' = 2, 'electron' = 3,
+        'positron' = 4).
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    bins : Iterable of Integral
+        The Particles to tally
+    id : int
+        Unique identifier for the filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+    @property
+    def bins(self):
+        return self._bins
+
+    @bins.setter
+    def bins(self, bins):
+        bins = np.atleast_1d(bins)
+        cv.check_iterable_type('filter bins', bins, (Integral, str))
+        for edge in bins:
+            if isinstance(edge, Integral):
+                cv.check_value('filter bin', edge, _PARTICLE_IDS.values())
+            else:
+                cv.check_value('filter bin', edge, _PARTICLE_IDS.keys())
+        bins = np.atleast_1d([b if isinstance(b, Integral) else _PARTICLE_IDS[b]
+                              for b in bins])
+        self._bins = bins
+
+
 class MeshFilter(Filter):
     """Bins tally event locations onto a regular, rectangular mesh.
 

openmc/particle_restart.py

@@ -27,6 +27,8 @@ class Particle(object):
         Type of simulation (criticality or fixed source)
     id : long
         Identifier of the particle
+    type : int
+        Particle type (1 = neutron, 2 = photon, 3 = electron, 4 = positron)
     weight : float
         Weight of the particle
     energy : float
@@ -65,6 +67,10 @@ class Particle(object):
     def id(self):
         return self._f['id'].value
 
+    @property
+    def type(self):
+        return self._f['type'].value
+
     @property
     def n_particles(self):
         return self._f['n_particles'].value

openmc/settings.py

@@ -30,12 +30,16 @@ class Settings(object):
         Indicate whether fission neutrons should be created or not.
     cutoff : dict
         Dictionary defining weight cutoff and energy cutoff. The dictionary may
-        have three keys, 'weight', 'weight_avg' and 'energy'. Value for 'weight'
+        have six keys, 'weight', 'weight_avg', 'energy_neutron', 'energy_photon',
+        'energy_electron', and 'energy_positron'. Value for 'weight'
         should be a float indicating weight cutoff below which particle undergo
         Russian roulette. Value for 'weight_avg' should be a float indicating
         weight assigned to particles that are not killed after Russian
         roulette. Value of energy should be a float indicating energy in eV
-        below which particle will be killed.
+        below which particle type will be killed.
+    electron_treatment : {'led', 'ttb'}
+        Whether to deposit all energy from electrons locally ('led') or create
+        secondary bremsstrahlung photons ('ttb').
     energy_mode : {'continuous-energy', 'multi-group'}
         Set whether the calculation should be continuous-energy or multi-group.
     entropy_mesh : openmc.Mesh
@@ -68,6 +72,8 @@ class Settings(object):
         :tallies: Whether the 'tallies.out' file should be written (bool)
     particles : int
         Number of particles per generation
+    photon_transport : bool
+        Whether to use photon transport.
     ptables : bool
         Determine whether probability tables are used.
     resonance_scattering : dict
@@ -171,7 +177,9 @@ class Settings(object):
 
         self._confidence_intervals = None
         self._cross_sections = None
+        self._electron_treatment = None
         self._multipole_library = None
+        self._photon_transport = None
         self._ptables = None
         self._run_cmfd = None
         self._seed = None
@@ -256,10 +264,18 @@ class Settings(object):
     def confidence_intervals(self):
         return self._confidence_intervals
 
+    @property
+    def electron_treatment(self):
+        return self._electron_treatment
+
     @property
     def ptables(self):
         return self._ptables
 
+    @property
+    def photon_transport(self):
+        return self._photon_transport
+
     @property
     def run_cmfd(self):
         return self._run_cmfd
@@ -485,6 +501,16 @@ class Settings(object):
         cv.check_type('confidence interval', confidence_intervals, bool)
         self._confidence_intervals = confidence_intervals
 
+    @electron_treatment.setter
+    def electron_treatment(self, electron_treatment):
+        cv.check_value('electron treatment', electron_treatment, ['led', 'ttb'])
+        self._electron_treatment = electron_treatment
+
+    @photon_transport.setter
+    def photon_transport(self, photon_transport):
+        cv.check_type('photon transport', photon_transport, bool)
+        self._photon_transport = photon_transport
+
     @ptables.setter
     def ptables(self, ptables):
         cv.check_type('probability tables', ptables, bool)
@@ -514,16 +540,16 @@ class Settings(object):
             raise ValueError(msg)
         for key in cutoff:
             if key == 'weight':
-                cv.check_type('weight cutoff', cutoff['weight'], Real)
-                cv.check_greater_than('weight cutoff', cutoff['weight'], 0.0)
+                cv.check_type('weight cutoff', cutoff[key], Real)
+                cv.check_greater_than('weight cutoff', cutoff[key], 0.0)
             elif key == 'weight_avg':
-                cv.check_type('average survival weight', cutoff['weight_avg'],
-                              Real)
+                cv.check_type('average survival weight', cutoff[key], Real)
                 cv.check_greater_than('average survival weight',
-                                      cutoff['weight_avg'], 0.0)
-            elif key == 'energy':
-                cv.check_type('energy cutoff', cutoff['energy'], Real)
-                cv.check_greater_than('energy cutoff', cutoff['energy'], 0.0)
+                                      cutoff[key], 0.0)
+            elif key in ['energy_neutron', 'energy_photon', 'energy_electron',
+                         'energy_positron']:
+                cv.check_type('energy cutoff', cutoff[key], Real)
+                cv.check_greater_than('energy cutoff', cutoff[key], 0.0)
             else:
                 msg = 'Unable to set cutoff to "{0}" which is unsupported by '\
                       'OpenMC'.format(key)
@@ -777,6 +803,16 @@ class Settings(object):
             element = ET.SubElement(root, "confidence_intervals")
             element.text = str(self._confidence_intervals).lower()
 
+    def _create_electron_treatment_subelement(self, root):
+        if self._electron_treatment is not None:
+            element = ET.SubElement(root, "electron_treatment")
+            element.text = str(self._electron_treatment)
+
+    def _create_photon_transport_subelement(self, root):
+        if self._photon_transport is not None:
+            element = ET.SubElement(root, "photon_transport")
+            element.text = str(self._photon_transport).lower()
+
     def _create_ptables_subelement(self, root):
         if self._ptables is not None:
             element = ET.SubElement(root, "ptables")
@@ -800,17 +836,9 @@ class Settings(object):
     def _create_cutoff_subelement(self, root):
         if self._cutoff is not None:
             element = ET.SubElement(root, "cutoff")
-            if 'weight' in self._cutoff:
-                subelement = ET.SubElement(element, "weight")
-                subelement.text = str(self._cutoff['weight'])
-
-            if 'weight_avg' in self._cutoff:
-                subelement = ET.SubElement(element, "weight_avg")
-                subelement.text = str(self._cutoff['weight_avg'])
-
-            if 'energy' in self._cutoff:
-                subelement = ET.SubElement(element, "energy")
-                subelement.text = str(self._cutoff['energy'])
+            for key, value in self._cutoff.items():
+                subelement = ET.SubElement(element, key)
+                subelement.text = str(value)
 
     def _create_entropy_mesh_subelement(self, root):
         if self.entropy_mesh is not None:
@@ -941,8 +969,10 @@ class Settings(object):
         self._create_statepoint_subelement(root_element)
         self._create_sourcepoint_subelement(root_element)
         self._create_confidence_intervals(root_element)
+        self._create_electron_treatment_subelement(root_element)
         self._create_energy_mode_subelement(root_element)
         self._create_max_order_subelement(root_element)
+        self._create_photon_transport_subelement(root_element)
         self._create_ptables_subelement(root_element)
         self._create_run_cmfd_subelement(root_element)
         self._create_seed_subelement(root_element)

openmc/source.py

@@ -22,6 +22,8 @@ class Source(object):
         Source file from which sites should be sampled
     strength : Real
         Strength of the source
+    particle : {'neutron', 'photon'}
+        Source particle type
 
     Attributes
     ----------
@@ -35,10 +37,13 @@ class Source(object):
         Source file from which sites should be sampled
     strength : Real
         Strength of the source
+    particle : {'neutron', 'photon'}
+        Source particle type
 
     """
 
-    def __init__(self, space=None, angle=None, energy=None, filename=None, strength=1.0):
+    def __init__(self, space=None, angle=None, energy=None, filename=None,
+                 strength=1.0, particle='neutron'):
         self._space = None
         self._angle = None
         self._energy = None
@@ -53,6 +58,7 @@ class Source(object):
         if filename is not None:
             self.file = filename
         self.strength = strength
+        self.particle = particle
 
     @property
     def file(self):
@@ -74,6 +80,10 @@ class Source(object):
     def strength(self):
         return self._strength
 
+    @property
+    def particle(self):
+        return self._particle
+
     @file.setter
     def file(self, filename):
         cv.check_type('source file', filename, str)
@@ -100,6 +110,11 @@ class Source(object):
         cv.check_greater_than('source strength', strength, 0.0, True)
         self._strength = strength
 
+    @particle.setter
+    def particle(self, particle):
+        cv.check_value('source particle', particle, ['neutron', 'photon'])
+        self._particle = particle
+
     def to_xml_element(self):
         """Return XML representation of the source
 
@@ -111,6 +126,8 @@ class Source(object):
         """
         element = ET.Element("source")
         element.set("strength", str(self.strength))
+        if self.particle != 'neutron':
+            element.set("particle", self.particle)
         if self.file is not None:
             element.set("file", self.file)
         if self.space is not None:

openmc/statepoint.py

@@ -83,6 +83,8 @@ class StatePoint(object):
         Number of tally realizations
     path : str
         Working directory for simulation
+    photon_transport : bool
+        Indicate whether photon transport is active
     run_mode : str
         Simulation run mode, e.g. 'eigenvalue'
     runtime : dict
@@ -322,6 +324,10 @@ class StatePoint(object):
     def path(self):
         return self._f.attrs['path'].decode()
 
+    @property
+    def photon_transport(self):
+        return self._f.attrs['photon_transport'] > 0
+
     @property
     def run_mode(self):
         return self._f['run_mode'].value.decode()

scripts/openmc-get-nndc-data

@@ -1,5 +1,11 @@
 #!/usr/bin/env python
 
+"""
+Download ENDF/B-VII.1 incident neutron ACE data and incident photon ENDF data
+from NNDC and convert it to an HDF5 library for use with OpenMC. This data is
+used for OpenMC's regression test suite.
+"""
+
 import os
 import shutil
 import subprocess
@@ -13,31 +19,25 @@ from urllib.request import urlopen
 import openmc.data
 
 
-description = """
-Download ENDF/B-VII.1 ACE data from NNDC and convert it to an HDF5 library for
-use with OpenMC. This data is used for OpenMC's regression test suite.
-
-"""
-
-
 class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
                       argparse.RawDescriptionHelpFormatter):
     pass
 
 parser = argparse.ArgumentParser(
-    description=description,
+    description=__doc__,
     formatter_class=CustomFormatter
 )
 parser.add_argument('-b', '--batch', action='store_true',
                     help='supresses standard in')
+parser.add_argument('-n', '--neutron-only', action='store_true',
+                    help='Whether to exclude photon interaction/atomic data')
 parser.add_argument('--libver', choices=['earliest', 'latest'],
                     default='earliest', help="Output HDF5 versioning. Use "
                     "'earliest' for backwards compatibility or 'latest' for "
                     "performance")
 args = parser.parse_args()
 
-
-baseUrl = 'http://www.nndc.bnl.gov/endf/b7.1/aceFiles/'
+base_url = 'http://www.nndc.bnl.gov/endf/b7.1/aceFiles/'
 files = ['ENDF-B-VII.1-neutron-293.6K.tar.gz',
          'ENDF-B-VII.1-tsl.tar.gz']
 checksums = ['9729a17eb62b75f285d8a7628ace1449',
@@ -47,32 +47,29 @@ block_size = 16384
 # ==============================================================================
 # DOWNLOAD FILES FROM NNDC SITE
 
-filesComplete = []
+files_complete = []
 for f in files:
     # Establish connection to URL
-    url = baseUrl + f
+    url = base_url + f
     req = urlopen(url)
 
     # Get file size from header
-    if sys.version_info[0] < 3:
-        file_size = int(req.info().getheaders('Content-Length')[0])
-    else:
-        file_size = req.length
+    file_size = req.length
     downloaded = 0
 
     # Check if file already downloaded
     if os.path.exists(f):
         if os.path.getsize(f) == file_size:
             print('Skipping ' + f)
-            filesComplete.append(f)
+            files_complete.append(f)
             continue
         else:
-            overwrite = input('Overwrite {0}? ([y]/n) '.format(f))
+            overwrite = input('Overwrite {}? ([y]/n) '.format(f))
             if overwrite.lower().startswith('n'):
                 continue
 
     # Copy file to disk
-    print('Downloading {0}... '.format(f), end='')
+    print('Downloading {}... '.format(f), end='')
     with open(f, 'wb') as fh:
         while True:
             chunk = req.read(block_size)
@@ -83,7 +80,7 @@ for f in files:
                 downloaded, downloaded * 100. / file_size)
             print(status + chr(8)*len(status), end='')
         print('')
-        filesComplete.append(f)
+        files_complete.append(f)
 
 # ==============================================================================
 # VERIFY MD5 CHECKSUMS
@@ -101,18 +98,18 @@ for f, checksum in zip(files, checksums):
 # EXTRACT FILES FROM TGZ
 
 for f in files:
-    if f not in filesComplete:
+    if f not in files_complete:
         continue
 
     # Extract files
     suffix = f[f.rindex('-') + 1:].rstrip('.tar.gz')
     with tarfile.open(f, 'r') as tgz:
-        print('Extracting {0}...'.format(f))
+        print('Extracting {}...'.format(f))
         tgz.extractall(path='nndc/' + suffix)
 
 # Move ACE files down one level
 for filename in glob.glob('nndc/293.6K/ENDF-B-VII.1-neutron-293.6K/*'):
-    shutil.move(filename, 'nndc/293.6K/')
+    shutil.move(filename, 'nndc/293.6K/' + os.path.basename(filename))
 
 # ==============================================================================
 # FIX ZAID ASSIGNMENTS FOR VARIOUS S(A,B) TABLES
@@ -142,7 +139,7 @@ else:
 if not response or response.lower().startswith('y'):
     for f in files:
         if os.path.exists(f):
-            print('Removing {0}...'.format(f))
+            print('Removing {}...'.format(f))
             os.remove(f)
 
 # ==============================================================================
@@ -162,3 +159,10 @@ subprocess.call([ace2hdf5,
                  '-d', 'nndc_hdf5',
                  '--fission_energy_release', fer_file,
                  '--libver', args.libver] + ace_files)
+
+# Generate photo interaction library files
+if not args.neutron_only:
+    pwd = os.path.dirname(os.path.realpath(__file__))
+    photo_endf = os.path.join(pwd, 'openmc-get-photon-data')
+    subprocess.call([photo_endf, '-c', 'cross_sections.xml'],
+                    cwd='nndc_hdf5')

scripts/openmc-get-photon-data

@@ -0,0 +1,106 @@
+#!/usr/bin/env python
+
+"""
+Download ENDF/B-VII.1 ENDF data from NNDC for photo-atomic and atomic
+relaxation data and convert it to an HDF5 library for use with OpenMC.
+This data is used for photon transport in OpenMC.
+"""
+
+import os
+import sys
+import shutil
+import zipfile
+import argparse
+from io import BytesIO
+from urllib.request import urlopen
+
+import openmc.data
+
+
+class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
+                      argparse.RawDescriptionHelpFormatter):
+    pass
+
+parser = argparse.ArgumentParser(
+    description=__doc__,
+    formatter_class=CustomFormatter
+)
+parser.add_argument('-c', '--cross-sections',
+                    help='cross_sections.xml file to append libraries to')
+args = parser.parse_args()
+
+base_url = 'http://www.nndc.bnl.gov/endf/b7.1/zips/'
+files = ['ENDF-B-VII.1-photoat.zip', 'ENDF-B-VII.1-atomic_relax.zip']
+block_size = 16384
+
+# ==============================================================================
+# DOWNLOAD FILES FROM NNDC SITE
+
+if not os.path.exists('photon_hdf5'):
+    os.mkdir('photon_hdf5')
+
+for f in files:
+    # Establish connection to URL
+    url = base_url + f
+    req = urlopen(url)
+
+    # Get file size from header
+    file_size = req.length
+    downloaded = 0
+
+    # Check if file already downloaded
+    if os.path.exists(f):
+        if os.path.getsize(f) == file_size:
+            print('Skipping ' + f)
+            continue
+        else:
+            overwrite = input('Overwrite {}? ([y]/n) '.format(f))
+            if overwrite.lower().startswith('n'):
+                continue
+
+    # Copy file to disk
+    print('Downloading {}... '.format(f), end='')
+    with open(f, 'wb') as fh:
+        while True:
+            chunk = req.read(block_size)
+            if not chunk: break
+            fh.write(chunk)
+            downloaded += len(chunk)
+            status = '{0:10}  [{1:3.2f}%]'.format(
+                downloaded, downloaded * 100. / file_size)
+            print(status + chr(8)*len(status), end='')
+        print('')
+
+# ==============================================================================
+# EXTRACT FILES
+
+for f in files:
+    print('Extracting {0}...'.format(f))
+    zipfile.ZipFile(f).extractall()
+
+# ==============================================================================
+# GENERATE HDF5 DATA LIBRARY
+
+# If previous cross_sections.xml was specified, load it in
+if args.cross_sections is not None:
+    lib_path = args.cross_sections
+    library = openmc.data.DataLibrary.from_xml(lib_path)
+else:
+    lib_path = os.path.join('photon_hdf5', 'cross_sections.xml')
+    library = openmc.data.DataLibrary()
+
+for z in range(1, 101):
+    element = openmc.data.ATOMIC_SYMBOL[z]
+    print('Generating HDF5 file for Z={} ({})...'.format(z, element))
+
+    # Generate instance of IncidentPhoton
+    photo_file = os.path.join('photoat', 'photoat-{:03}_{}_000.endf'.format(z, element))
+    atom_file = os.path.join('atomic_relax', 'atom-{:03}_{}_000.endf'.format(z, element))
+    f = openmc.data.IncidentPhoton.from_endf(photo_file, atom_file)
+
+    # Write HDF5 file and register it
+    hdf5_file = os.path.join('photon_hdf5', element + '.h5')
+    f.export_to_hdf5(hdf5_file, 'w')
+    library.register_file(hdf5_file)
+
+library.export_to_xml(lib_path)

scripts/openmc-make-compton

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import tarfile
+from urllib.request import urlopen
+
+import numpy as np
+import h5py
+
+
+base_url = 'http://geant4.cern.ch/support/source/'
+filename = 'G4EMLOW.6.48.tar.gz'
+block_size = 16384
+
+# ==============================================================================
+# DOWNLOAD FILES FROM GEANT4 SITE
+
+# Establish connection to URL
+req = urlopen(base_url + filename)
+
+# Get file size from header
+file_size = req.length
+downloaded = 0
+
+# Check if file already downloaded
+download = True
+if os.path.exists(filename):
+    if os.path.getsize(filename) == file_size:
+        print('Already downloaded ' + filename)
+        download = False
+    else:
+        overwrite = input('Overwrite {}? ([y]/n) '.format(filename))
+        if overwrite.lower().startswith('n'):
+            download = False
+
+if download:
+    # Copy file to disk
+    print('Downloading {}... '.format(filename), end='')
+    with open(filename, 'wb') as fh:
+        while True:
+            chunk = req.read(block_size)
+            if not chunk: break
+            fh.write(chunk)
+            downloaded += len(chunk)
+            status = '{0:10}  [{1:3.2f}%]'.format(
+                downloaded, downloaded * 100. / file_size)
+            print(status + chr(8)*len(status), end='')
+        print('')
+
+# ==============================================================================
+# EXTRACT FILES FROM TGZ
+
+if not os.path.isdir('G4EMLOW6.48'):
+    with tarfile.open(filename, 'r') as tgz:
+        print('Extracting {0}...'.format(filename))
+        tgz.extractall()
+
+# ==============================================================================
+# GENERATE COMPTON PROFILE HDF5 FILE
+
+print('Generating compton_profiles.h5...')
+
+shell_file = os.path.join('G4EMLOW6.48', 'doppler', 'shell-doppler.dat')
+
+with open(shell_file, 'r') as shell:
+    with h5py.File('compton_profiles.h5', 'w') as f:
+        # Read/write electron momentum values
+        pz = np.loadtxt(os.path.join('G4EMLOW6.48', 'doppler', 'p-biggs.dat'))
+        f.create_dataset('pz', data=pz)
+
+        for Z in range(1, 101):
+            # Create group for this element
+            group = f.create_group('{:03}'.format(Z))
+
+            # Read data into one long array
+            path = os.path.join('G4EMLOW6.48', 'doppler', 'profile-{}.dat'.format(Z))
+            J = np.fromstring(open(path, 'r').read(), sep=' ')
+
+            # Determine number of electron shells and reshape
+            n_shells = J.size // 31
+            J.shape = (n_shells, 31)
+
+            # Write Compton profile for this Z
+            group.create_dataset('J', data=J)
+
+            # Determine binding energies and number of electrons for each shell
+            num_electrons = []
+            binding_energy = []
+            while True:
+                words = shell.readline().split()
+                if words[0] == '-1':
+                    break
+                num_electrons.append(float(words[0]))
+                binding_energy.append(float(words[1]))
+
+            # Write binding energies and number of electrons
+            group.create_dataset('num_electrons', data=num_electrons)
+            group.create_dataset('binding_energy', data=binding_energy)

scripts/openmc-make-stopping-powers

@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+
+from urllib.parse import urlencode
+from urllib.request import urlopen
+from lxml import html
+
+import numpy as np
+import h5py
+from openmc.data import ATOMIC_SYMBOL
+
+
+base_url = 'https://physics.nist.gov/cgi-bin/Star/e_table-t.pl'
+energies = np.logspace(-3, 3, 200)
+data = {'matno': '', 'Energies': '\n'.join(str(x) for x in energies)}
+columns = {1: 's_collision', 2: 's_radiative'}
+
+# ==============================================================================
+# SCRAPE DATA FROM ESTAR SITE AND GENERATE STOPPING POWER HDF5 FILE
+
+print('Generating stopping_powers.h5...')
+
+with h5py.File('stopping_powers.h5', 'w') as f:
+
+    # Write energies
+    f.create_dataset('energy', data=energies)
+
+    for Z in range(1, 99):
+        print('Processing {} data...'.format(ATOMIC_SYMBOL[Z]))
+
+        # Update form-encoded data to send in POST request for this element
+        data['matno'] = '{:03}'.format(Z)
+        payload = urlencode(data).encode("utf-8")
+
+        # Retrieve data from ESTAR site
+        r = urlopen(url=base_url, data=payload).read()
+
+        # Remove text and reformat data
+        r = html.fromstring(r).xpath('//pre//text()')
+        values = np.fromstring(' '.join(r[12:-5]), sep=' ').reshape((-1, 5)).T
+
+        # Create group for this element
+        group = f.create_group('{:03}'.format(Z))
+
+        # Write the mean excitation energy
+        attributes = np.fromstring(r[3], sep=' ')
+        group.attrs['I'] = attributes[2]
+
+        # Write collision and radiative stopping powers
+        for i in columns:
+            group.create_dataset(columns[i], data=values[i])

scripts/openmc-update-inputs

@@ -248,8 +248,7 @@ def update_materials(root):
                 # If a nuclide name is in the ZAID notation (e.g., a number),
                 # convert it to the proper nuclide name.
                 if nucname.strip().isnumeric():
-                    nucname = \
-                      openmc.data.neutron._get_metadata(int(nucname))[0]
+                    nucname = openmc.data.ace.get_metadata(int(nucname))[0]
                 nucname = nucname.replace('Nat', '0')
                 if nucname.endswith('m'):
                     nucname = nucname[:-1] + '_m1'

setup.py

@@ -32,7 +32,7 @@ kwargs = {
     # Data files and librarries
     'package_data': {
         'openmc.capi': ['libopenmc.{}'.format(suffix)],
-        'openmc.data': ['mass.mas12', 'fission_Q_data_endfb71.h5']
+        'openmc.data': ['mass.mas12', '*.h5']
     },
 
     # Metadata

src/api.F90

@@ -120,9 +120,10 @@ contains
     check_overlaps = .false.
     confidence_intervals = .false.
     create_fission_neutrons = .true.
-    energy_cutoff = ZERO
-    energy_max_neutron = INFINITY
-    energy_min_neutron = ZERO
+    electron_treatment = ELECTRON_LED
+    energy_cutoff(:) = [ZERO, 1000.0_8, ZERO, ZERO]
+    energy_max(:) = [INFINITY, INFINITY]
+    energy_min(:) = [ZERO, ZERO]
     entropy_on = .false.
     gen_per_batch = 1
     index_entropy_mesh = -1
@@ -139,6 +140,7 @@ contains
     output_summary = .true.
     output_tallies = .true.
     particle_restart_run = .false.
+    photon_transport = .false.
     pred_batches = .false.
     reduce_tallies = .true.
     res_scat_on = .false.
@@ -310,6 +312,7 @@ contains
   subroutine free_memory()
 
     use cmfd_header
+    use photon_header
     use plot_header
     use sab_header
     use settings
@@ -326,6 +329,7 @@ contains
     call free_memory_volume()
     call free_memory_simulation()
     call free_memory_nuclide()
+    call free_memory_photon()
     call free_memory_settings()
     call free_memory_sab()
     call free_memory_source()

src/bank_header.F90

@@ -18,6 +18,7 @@ module bank_header
     real(C_DOUBLE) :: uvw(3)        ! diretional cosines
     real(C_DOUBLE) :: E             ! energy / energy group if in MG mode.
     integer(C_INT) :: delayed_group ! delayed group
+    integer(C_INT) :: particle      ! particle type (neutron, photon, etc.)
   end type Bank
 
   ! Source and fission bank

src/cmfd_input.F90

@@ -116,7 +116,7 @@ contains
       end if
     else
       if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
-      cmfd % egrid = [ ZERO, energy_max_neutron ]
+      cmfd % egrid = [ ZERO, energy_max(NEUTRON) ]
       cmfd % indices(4) = 1 ! one energy group
     end if
 

src/constants.F90

@@ -78,6 +78,9 @@ module constants
        MASS_NEUTRON     = 1.00866491588_8,   & ! mass of a neutron in amu
        MASS_NEUTRON_EV  = 939.5654133e6_8,   & ! mass of a neutron in eV/c^2
        MASS_PROTON      = 1.007276466879_8,  & ! mass of a proton in amu
+       MASS_ELECTRON_EV = 0.5109989461e6_8,  & ! electron mass energy equivalent in eV/c^2
+       FINE_STRUCTURE   = 137.035999139_8,   & ! inverse fine structure constant
+       PLANCK_C         = 1.2398419739062977e4_8,& ! Planck's constant times c in eV-Angstroms
        AMU              = 1.660539040e-27_8, & ! 1 amu in kg
        C_LIGHT          = 2.99792458e8_8,    & ! speed of light in m/s
        N_AVOGADRO       = 0.6022140857_8,    & ! Avogadro's number in 10^24/mol
@@ -91,6 +94,14 @@ module constants
        FOUR             = 4.0_8
   complex(8), parameter :: ONEI = (ZERO, ONE)
 
+  ! Electron subshell labels
+  character(3), parameter :: SUBSHELLS(39) = [ &
+       'K  ', 'L1 ', 'L2 ', 'L3 ', 'M1 ', 'M2 ', 'M3 ', 'M4 ', 'M5 ', &
+       'N1 ', 'N2 ', 'N3 ', 'N4 ', 'N5 ', 'N6 ', 'N7 ', 'O1 ', 'O2 ', &
+       'O3 ', 'O4 ', 'O5 ', 'O6 ', 'O7 ', 'O8 ', 'O9 ', 'P1 ', 'P2 ', &
+       'P3 ', 'P4 ', 'P5 ', 'P6 ', 'P7 ', 'P8 ', 'P9 ', 'P10', 'P11', &
+       'Q1 ', 'Q2 ', 'Q3 ']
+
   ! ============================================================================
   ! GEOMETRY-RELATED CONSTANTS
 
@@ -172,7 +183,8 @@ module constants
   integer, parameter :: &
        NEUTRON  = 1, &
        PHOTON   = 2, &
-       ELECTRON = 3
+       ELECTRON = 3, &
+       POSITRON = 4
 
   ! Angular distribution type
   integer, parameter :: &
@@ -219,7 +231,9 @@ module constants
        N_3HEA  = 193, N_4N2P  = 194, N_4N2A = 195, N_4NPA  = 196, N_3P    = 197, &
        N_N3P   = 198, N_3N2PA = 199, N_5N2P = 200, N_P0    = 600, N_PC    = 649, &
        N_D0    = 650, N_DC    = 699, N_T0   = 700, N_TC    = 749, N_3HE0  = 750, &
-       N_3HEC  = 799, N_A0    = 800, N_AC   = 849, N_2N0   = 875, N_2NC   = 891
+       N_3HEC  = 799, N_A0    = 800, N_AC   = 849, N_2N0   = 875, N_2NC   = 891, &
+       COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, &
+       PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517
 
   ! Depletion reactions
   integer, parameter :: DEPLETION_RX(6) = [N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N]
@@ -343,7 +357,7 @@ module constants
   integer, parameter :: NO_BIN_FOUND = -1
 
   ! Tally filter and map types
-  integer, parameter :: N_FILTER_TYPES = 20
+  integer, parameter :: N_FILTER_TYPES = 21
   integer, parameter :: &
        FILTER_UNIVERSE       = 1,  &
        FILTER_MATERIAL       = 2,  &
@@ -364,7 +378,8 @@ module constants
        FILTER_LEGENDRE       = 17, &
        FILTER_SPH_HARMONICS  = 18, &
        FILTER_SPTL_LEGENDRE  = 19, &
-       FILTER_ZERNIKE        = 20
+       FILTER_ZERNIKE        = 20, &
+       FILTER_PARTICLE       = 21
 
   ! Mesh types
   integer, parameter :: &
@@ -433,6 +448,7 @@ module constants
   integer(C_INT), bind(C, name='STREAM_SOURCE') :: STREAM_SOURCE
   integer(C_INT), bind(C, name='STREAM_URR_PTABLE') :: STREAM_URR_PTABLE
   integer(C_INT), bind(C, name='STREAM_VOLUME') :: STREAM_VOLUME
+  integer(C_INT), bind(C, name='STREAM_PHOTON') :: STREAM_PHOTON
   integer(C_INT64_T), parameter :: DEFAULT_SEED = 1_8
 
   ! ============================================================================
@@ -455,6 +471,11 @@ module constants
        MODE_PARTICLE    = 4, & ! Particle restart mode
        MODE_VOLUME      = 5    ! Volume calculation mode
 
+  ! Electron treatments
+  integer, parameter :: &
+       ELECTRON_LED     = 1, & ! Local Energy Deposition
+       ELECTRON_TTB     = 2    ! Thick Target Bremsstrahlung
+
   !=============================================================================
   ! CMFD CONSTANTS
 

src/endf.F90

@@ -160,6 +160,20 @@ contains
       string = '(n,Xa)'
     case (444)
       string = '(damage)'
+    case (COHERENT)
+      string = 'coherent scatter'
+    case (INCOHERENT)
+      string = 'incoherent scatter'
+    case (PAIR_PROD_ELEC)
+      string = 'pair production, electron'
+    case (PAIR_PROD)
+      string = 'pair production'
+    case (PAIR_PROD_NUC)
+      string = 'pair production, nuclear'
+    case (PHOTOELECTRIC)
+      string = 'photoelectric'
+    case (534 : 572)
+      string = 'photoelectric, ' // trim(SUBSHELLS(MT - 533)) // ' subshell'
     case (600 : 648)
       string = '(n,p' // trim(to_str(MT-600)) // ')'
     case (649)

src/energy_distribution.F90

@@ -331,8 +331,8 @@ contains
       c_k = c_k1
     end do
 
-    ! Check to make sure k is <= NP - 1
-    k = min(k, n_energy_out - 1)
+    ! Check to make sure 1 <= k <= NP - 1
+    k = max(1, min(k, n_energy_out - 1))
 
     E_l_k = this%distribution(l)%e_out(k)
     p_l_k = this%distribution(l)%p(k)
@@ -359,7 +359,7 @@ contains
     end if
 
     ! Now interpolate between incident energy bins i and i + 1
-    if (.not. histogram_interp) then
+    if (.not. histogram_interp .and. n_energy_out > 1) then
       if (l == i) then
         E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
       else

src/geometry.F90

@@ -117,8 +117,14 @@ contains
     end do
     j = p % n_coord
 
-    ! set size of list to search
+    ! Determine universe (if not yet set, use root universe)
     i_universe = p % coord(j) % universe
+    if (i_universe == NONE) then
+      p % coord(j) % universe = root_universe
+      i_universe = root_universe
+    end if
+
+    ! set size of list to search
     if (present(search_cells)) then
       use_search_cells = .true.
       n = size(search_cells)

src/input_xml.F90

@@ -22,6 +22,7 @@ module input_xml
   use mgxs_interface
   use nuclide_header
   use output,           only: title, header, print_plot
+  use photon_header
   use plot_header
   use random_lcg,       only: prn, openmc_set_seed
   use surface_header
@@ -389,6 +390,30 @@ contains
       call openmc_set_seed(seed)
     end if
 
+    ! Check for electron treatment
+    if (check_for_node(root, "electron_treatment")) then
+      call get_node_value(root, "electron_treatment", temp_str)
+      select case (to_lower(temp_str))
+      case ("led")
+        electron_treatment = ELECTRON_LED
+      case ("ttb")
+        electron_treatment = ELECTRON_TTB
+      case default
+        call fatal_error("Unrecognized electron treatment: " // &
+             trim(temp_str) // ".")
+      end select
+    end if
+
+    ! Check for photon transport
+    if (check_for_node(root, "photon_transport")) then
+      call get_node_value(root, "photon_transport", photon_transport)
+
+      if (.not. run_CE .and. photon_transport) then
+        call fatal_error("Photon transport is not currently supported &
+             &in Multi-group mode")
+      end if
+    end if
+
     ! Number of bins for logarithmic grid
     if (check_for_node(root, "log_grid_bins")) then
       call get_node_value(root, "log_grid_bins", n_log_bins)
@@ -425,6 +450,7 @@ contains
     if (n == 0) then
       ! Default source is isotropic point source at origin with Watt spectrum
       allocate(external_source(1))
+      external_source % particle = NEUTRON
       external_source % strength = ONE
 
       allocate(SpatialPoint :: external_source(1) % space)
@@ -476,8 +502,21 @@ contains
       if (check_for_node(node_cutoff, "weight_avg")) then
         call get_node_value(node_cutoff, "weight_avg", weight_survive)
       end if
-      if (check_for_node(node_cutoff, "energy")) then
-        call get_node_value(node_cutoff, "energy", energy_cutoff)
+      if (check_for_node(node_cutoff, "energy_neutron")) then
+        call get_node_value(node_cutoff, "energy_neutron", energy_cutoff(1))
+      elseif (check_for_node(node_cutoff, "energy")) then
+        call warning("The use of an <energy> cutoff is deprecated and should &
+             &be replaced by <energy_neutron>.")
+        call get_node_value(node_cutoff, "energy", energy_cutoff(1))
+      end if
+      if (check_for_node(node_cutoff, "energy_photon")) then
+        call get_node_value(node_cutoff, "energy_photon", energy_cutoff(2))
+      end if
+      if (check_for_node(node_cutoff, "energy_electron")) then
+        call get_node_value(node_cutoff, "energy_electron", energy_cutoff(3))
+      end if
+      if (check_for_node(node_cutoff, "energy_positron")) then
+        call get_node_value(node_cutoff, "energy_positron", energy_cutoff(4))
       end if
     end if
 
@@ -1450,9 +1489,11 @@ contains
     integer :: n_sab          ! number of sab tables for a material
     integer :: i_library      ! index in libraries array
     integer :: index_nuclide  ! index in nuclides
+    integer :: index_element  ! index in elements
     integer :: index_sab      ! index in sab_tables
     logical :: file_exists    ! does materials.xml exist?
-    character(20)           :: name         ! name of nuclide, e.g. 92235.03c
+    character(20)           :: name         ! name of nuclide, e.g. U235
+    character(3)            :: element      ! name of element, e.g. Zr
     character(MAX_WORD_LEN) :: units        ! units on density
     character(MAX_LINE_LEN) :: filename     ! absolute path to materials.xml
     character(MAX_LINE_LEN) :: temp_str     ! temporary string when reading
@@ -1501,6 +1542,7 @@ contains
 
     ! Initialize count for number of nuclides/S(a,b) tables
     index_nuclide = 0
+    index_element = 0
     index_sab = 0
 
     do i = 1, n_materials
@@ -1723,6 +1765,7 @@ contains
       mat % n_nuclides = n
       allocate(mat % names(n))
       allocate(mat % nuclide(n))
+      allocate(mat % element(n))
       allocate(mat % atom_density(n))
 
       ALL_NUCLIDES: do j = 1, mat % n_nuclides
@@ -1753,6 +1796,27 @@ contains
           mat % nuclide(j) = nuclide_dict % get(to_lower(name))
         end if
 
+        ! If the corresponding element hasn't been encountered yet and photon
+        ! transport will be used, we need to add its symbol to the element_dict
+        if (photon_transport) then
+          element = name(1:scan(name, '0123456789') - 1)
+
+          ! Make sure photon cross section data is available
+          if (.not. library_dict % has(to_lower(element))) then
+            call fatal_error("Could not find element " // trim(element) &
+                 // " in cross_sections data file!")
+          end if
+
+          if (.not. element_dict % has(element)) then
+            index_element = index_element + 1
+            mat % element(j) = index_element
+
+            call element_dict % set(element, index_element)
+          else
+            mat % element(j) = element_dict % get(element)
+          end if
+        end if
+
         ! Copy name and atom/weight percent
         mat % names(j) = name
         mat % atom_density(j) = densities % data(j)
@@ -1868,6 +1932,7 @@ contains
 
     ! Set total number of nuclides and S(a,b) tables
     n_nuclides = index_nuclide
+    n_elements = index_element
     n_sab_tables = index_sab
 
     ! Close materials XML file
@@ -1885,6 +1950,7 @@ contains
     integer :: i             ! loop over user-specified tallies
     integer :: j             ! loop over words
     integer :: k             ! another loop index
+    integer :: l             ! loop over bins
     integer :: filter_id     ! user-specified identifier for filter
     integer :: i_filt        ! index in filters array
     integer :: i_elem        ! index of entry in dictionary
@@ -2510,6 +2576,47 @@ contains
             end if
           end do
         end do
+
+        ! Check if tally is compatible with particle type
+        if (photon_transport) then
+          if (t % find_filter(FILTER_PARTICLE) == 0) then
+            do j = 1, n_scores
+              select case (t % score_bins(j))
+              case (SCORE_INVERSE_VELOCITY)
+                call fatal_error("Particle filter must be used with photon &
+                     &transport on and inverse velocity score")
+              case (SCORE_FLUX, SCORE_TOTAL, SCORE_SCATTER, SCORE_NU_SCATTER, &
+                   SCORE_ABSORPTION, SCORE_FISSION, SCORE_NU_FISSION, &
+                   SCORE_CURRENT, SCORE_EVENTS, SCORE_DELAYED_NU_FISSION, &
+                   SCORE_PROMPT_NU_FISSION, SCORE_DECAY_RATE)
+                call warning("Particle filter is not used with photon transport&
+                     & on and " // trim(to_str(t % score_bins(j))) // " score")
+              end select
+            end do
+          else
+            select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
+            type is (ParticleFilter)
+              do l = 1, filt % n_bins
+                if (filt % particles(l) == ELECTRON .or. filt % particles(l) == POSITRON) then
+                  t % estimator = ESTIMATOR_ANALOG
+                end if
+              end do
+            end select
+          end if
+        else
+          if (t % find_filter(FILTER_PARTICLE) > 0) then
+            select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
+            type is (ParticleFilter)
+              do l = 1, filt % n_bins
+                if (filt % particles(l) /= NEUTRON) then
+                  call warning("Particle filter other than NEUTRON used with &
+                       &photon transport turned off. All tallies for particle &
+                       &type " // trim(to_str(filt % particles(l))) // " will have no scores")
+                end if
+              end do
+            end select
+          end if
+        end if
       else
         call fatal_error("No <scores> specified on tally " &
              // trim(to_str(t % id)) // ".")
@@ -3319,6 +3426,8 @@ contains
           libraries(i) % type = LIBRARY_NEUTRON
         case ('thermal')
           libraries(i) % type = LIBRARY_THERMAL
+        case ('photon')
+          libraries(i) % type = LIBRARY_PHOTON
         end select
       else
         call fatal_error("Missing library type")
@@ -3408,8 +3517,8 @@ contains
     end do
 
     ! Get the minimum and maximum energies
-    energy_min_neutron = energy_bins(num_energy_groups + 1)
-    energy_max_neutron = energy_bins(1)
+    energy_min(NEUTRON) = energy_bins(num_energy_groups + 1)
+    energy_max(NEUTRON) = energy_bins(1)
 
     ! Get the datasets present in the library
     call get_groups(file_id, names)
@@ -3519,15 +3628,22 @@ contains
     integer :: i, j
     integer :: i_library
     integer :: i_nuclide
+    integer :: i_element
     integer :: i_sab
     integer(HID_T) :: file_id
     integer(HID_T) :: group_id
     logical :: mp_found     ! if windowed multipole libraries were found
     character(MAX_WORD_LEN) :: name
+    character(3) :: element
     type(SetChar) :: already_read
+    type(SetChar) :: element_already_read
 
     allocate(nuclides(n_nuclides))
+    allocate(elements(n_elements))
     allocate(sab_tables(n_sab_tables))
+    if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
+      allocate(ttb(n_materials))
+    end if
 
     ! Read cross sections
     do i = 1, size(materials)
@@ -3559,15 +3675,50 @@ contains
           ! Determine if minimum/maximum energy for this nuclide is greater/less
           ! than the previous
           if (size(nuclides(i_nuclide) % grid) >= 1) then
-            energy_min_neutron = max(energy_min_neutron, &
+            energy_min(NEUTRON) = max(energy_min(NEUTRON), &
                  nuclides(i_nuclide) % grid(1) % energy(1))
-            energy_max_neutron = min(energy_max_neutron, nuclides(i_nuclide) % &
+            energy_max(NEUTRON) = min(energy_max(NEUTRON), nuclides(i_nuclide) % &
                  grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy)))
           end if
 
           ! Add name and alias to dictionary
           call already_read % add(name)
 
+          ! Check if elemental data has been read, if needed
+          element = name(1:scan(name, '0123456789') - 1)
+          if (photon_transport) then
+            if (.not. element_already_read % contains(element)) then
+              ! Read photon interaction data from HDF5 photon library
+              i_library = library_dict % get(to_lower(element))
+              i_element = element_dict % get(element)
+              call write_message('Reading ' // trim(element) // ' from ' // &
+                   trim(libraries(i_library) % path), 6)
+
+              ! Open file and make sure version is sufficient
+              file_id = file_open(libraries(i_library) % path, 'r')
+              call check_data_version(file_id)
+
+              ! Read element data from HDF5
+              group_id = open_group(file_id, element)
+              call elements(i_element) % from_hdf5(group_id)
+              call close_group(group_id)
+              call file_close(file_id)
+
+              ! Determine if minimum/maximum energy for this element is
+              ! greater/less than the previous
+              if (size(elements(i_element) % energy) >= 1) then
+                energy_min(PHOTON) = max(energy_min(PHOTON), &
+                     exp(elements(i_element) % energy(1)))
+                energy_max(PHOTON) = min(energy_max(PHOTON), &
+                     exp(elements(i_element) % energy(size(elements(i_element) &
+                     % energy))))
+              end if
+
+              ! Add element to set
+              call element_already_read % add(element)
+            end if
+          end if
+
           ! Read multipole file into the appropriate entry on the nuclides array
           if (temperature_multipole) call read_multipole_data(i_nuclide)
         end if
@@ -3577,14 +3728,43 @@ contains
           materials(i) % fissionable = .true.
         end if
       end do
+
+      ! Generate material bremsstrahlung data for electrons and positrons
+      if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
+        call bremsstrahlung_init(ttb(i) % electron, i, ELECTRON)
+        call bremsstrahlung_init(ttb(i) % positron, i, POSITRON)
+      end if
     end do
 
+    if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
+      ! Deallocate element bremsstrahlung DCS and stopping power data since
+      ! only the material bremsstrahlung data is needed
+      do i = 1, size(elements)
+        if (allocated(elements(i) % stopping_power_collision)) &
+             deallocate(elements(i) % stopping_power_collision)
+        if (allocated(elements(i) % stopping_power_radiative)) &
+             deallocate(elements(i) % stopping_power_radiative)
+        if (allocated(elements(i) % dcs)) deallocate(elements(i) % dcs)
+        if (allocated(ttb_k_grid)) deallocate(ttb_k_grid)
+      end do
+
+      ! Determine if minimum/maximum energy for bremsstrahlung is greater/less
+      ! than the current minimum/maximum
+      if (size(ttb_e_grid) >= 1) then
+        energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1))
+        energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid)))
+      end if
+
+      ! Take logarithm of energies since they are log-log interpolated
+      ttb_e_grid = log(ttb_e_grid)
+    end if
+
     ! Set up logarithmic grid for nuclides
     do i = 1, size(nuclides)
-      call nuclides(i) % init_grid(energy_min_neutron, &
-           energy_max_neutron, n_log_bins)
+      call nuclides(i) % init_grid(energy_min(NEUTRON), &
+           energy_max(NEUTRON), n_log_bins)
     end do
-    log_spacing = log(energy_max_neutron/energy_min_neutron) / n_log_bins
+    log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins
 
     do i = 1, size(materials)
       ! Skip materials with no S(a,b) tables
@@ -3627,9 +3807,9 @@ contains
       ! grid has not been allocated
       if (size(nuclides(i) % grid) > 0) then
         if (nuclides(i) % grid(1) % energy(size(nuclides(i) % grid(1) % energy)) &
-             == energy_max_neutron) then
+             == energy_max(NEUTRON)) then
           call write_message("Maximum neutron transport energy: " // &
-               trim(to_str(energy_max_neutron)) // " eV for " // &
+               trim(to_str(energy_max(NEUTRON))) // " eV for " // &
                trim(adjustl(nuclides(i) % name)), 7)
           exit
         end if

src/material_header.F90

@@ -5,7 +5,10 @@ module material_header
   use constants
   use dict_header, only: DictIntInt
   use error
+  use math, only: spline, spline_integrate
   use nuclide_header
+  use particle_header, only: Particle
+  use photon_header
   use sab_header
   use simulation_header, only: log_spacing
   use stl_vector, only: VectorReal, VectorInt
@@ -14,6 +17,7 @@ module material_header
   implicit none
 
   private
+  public :: bremsstrahlung_init
   public :: free_memory_material
   public :: openmc_extend_materials
   public :: openmc_get_material_index
@@ -33,6 +37,7 @@ module material_header
     character(len=104)   :: name = ""       ! User-defined name
     integer              :: n_nuclides = 0  ! number of nuclides
     integer, allocatable :: nuclide(:)      ! index in nuclides array
+    integer, allocatable :: element(:)      ! index in elements array
     real(8)              :: density         ! total atom density in atom/b-cm
     real(C_DOUBLE), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
     real(8)              :: density_gpcc    ! total density in g/cm^3
@@ -66,6 +71,8 @@ module material_header
     procedure :: init_nuclide_index => material_init_nuclide_index
     procedure :: assign_sab_tables => material_assign_sab_tables
     procedure :: calculate_xs => material_calculate_xs
+    procedure, private :: calculate_neutron_xs
+    procedure, private :: calculate_photon_xs
   end type Material
 
   integer(C_INT32_T), public, bind(C) :: n_materials ! # of materials
@@ -251,18 +258,36 @@ contains
   end subroutine material_assign_sab_tables
 
 !===============================================================================
-! MATERIAL_CALCULATE_XS determines the macroscopic cross sections for the material the
-! particle is currently traveling through.
+! MATERIAL_CALCULATE_XS determines the macroscopic cross sections for the
+! material the particle is currently traveling through.
 !===============================================================================
 
-  subroutine material_calculate_xs(this, E, sqrtkT, micro_xs, nuclides, &
-                                   material_xs)
+  subroutine material_calculate_xs(this, p)
     class(Material), intent(in) :: this
-    real(8),         intent(in) :: E       ! Particle energy
-    real(8),         intent(in) :: sqrtkT  ! Last temperature sampled
-    type(Nuclide), allocatable, intent(in)           :: nuclides(:)
-    type(NuclideMicroXS), allocatable, intent(inout) :: micro_xs(:)  ! Cache for each nuclide
-    type(MaterialMacroXS), intent(inout)             :: material_xs  ! Cache for current material
+    type(Particle),  intent(in) :: p
+
+    ! Set all material macroscopic cross sections to zero
+    material_xs % total           = ZERO
+    material_xs % absorption      = ZERO
+    material_xs % fission         = ZERO
+    material_xs % nu_fission      = ZERO
+
+    if (p % type == NEUTRON) then
+      call this % calculate_neutron_xs(p)
+    elseif (p % type == PHOTON) then
+      call this % calculate_photon_xs(p)
+    end if
+
+  end subroutine material_calculate_xs
+
+!===============================================================================
+! CALCULATE_NEUTRON_XS determines the neutron cross section for the material the
+! particle is traveling through
+!===============================================================================
+
+  subroutine calculate_neutron_xs(this, p)
+    class(Material), intent(in) :: this
+    type(Particle),  intent(in) :: p
 
     integer :: i             ! loop index over nuclides
     integer :: i_nuclide     ! index into nuclides array
@@ -274,14 +299,8 @@ contains
     real(8) :: sab_frac      ! fraction of atoms affected by S(a,b)
     logical :: check_sab     ! should we check for S(a,b) table?
 
-    ! Set all material macroscopic cross sections to zero
-    material_xs % total          = ZERO
-    material_xs % absorption     = ZERO
-    material_xs % fission        = ZERO
-    material_xs % nu_fission     = ZERO
-
     ! Find energy index on energy grid
-    i_grid = int(log(E/energy_min_neutron)/log_spacing)
+    i_grid = int(log(p % E/energy_min(NEUTRON))/log_spacing)
 
     ! Determine if this material has S(a,b) tables
     check_sab = (this % n_sab > 0)
@@ -307,7 +326,7 @@ contains
 
           ! If particle energy is greater than the highest energy for the
           ! S(a,b) table, then don't use the S(a,b) table
-          if (E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
+          if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
             i_sab = 0
           end if
 
@@ -326,12 +345,12 @@ contains
       i_nuclide = this % nuclide(i)
 
       ! Calculate microscopic cross section for this nuclide
-      if (E /= micro_xs(i_nuclide) % last_E &
-           .or. sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
+      if (p % E /= micro_xs(i_nuclide) % last_E &
+           .or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
            .or. i_sab /= micro_xs(i_nuclide) % index_sab &
            .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
-        call nuclides(i_nuclide) % calculate_xs(i_sab, E, i_grid, &
-             sqrtkT, sab_frac, micro_xs(i_nuclide))
+        call nuclides(i_nuclide) % calculate_xs(i_sab, p % E, i_grid, &
+             p % sqrtkT, sab_frac, micro_xs(i_nuclide))
       end if
 
       ! ======================================================================
@@ -357,7 +376,68 @@ contains
            atom_density * micro_xs(i_nuclide) % nu_fission
     end do
 
-  end subroutine material_calculate_xs
+  end subroutine calculate_neutron_xs
+
+!===============================================================================
+! CALCULATE_PHOTON_XS determines the macroscopic photon cross sections for the
+! material the particle is currently traveling through.
+!===============================================================================
+
+  subroutine calculate_photon_xs(this, p)
+    class(Material), intent(in) :: this
+    type(Particle),  intent(in) :: p
+
+    integer :: i             ! loop index over nuclides
+    integer :: i_element     ! index into elements array
+    real(8) :: atom_density  ! atom density of a nuclide
+
+    material_xs % coherent        = ZERO
+    material_xs % incoherent      = ZERO
+    material_xs % photoelectric   = ZERO
+    material_xs % pair_production = ZERO
+
+    ! Add contribution from each nuclide in material
+    do i = 1, this % n_nuclides
+      ! ========================================================================
+      ! CALCULATE MICROSCOPIC CROSS SECTION
+
+      ! Determine microscopic cross sections for this nuclide
+      i_element = this % element(i)
+
+      ! Calculate microscopic cross section for this nuclide
+      if (p % E /= micro_photon_xs(i_element) % last_E) then
+        call elements(i_element) % calculate_xs(&
+             p % E, micro_photon_xs(i_element))
+      end if
+
+      ! ========================================================================
+      ! ADD TO MACROSCOPIC CROSS SECTION
+
+      ! Copy atom density of nuclide in material
+      atom_density = this % atom_density(i)
+
+      ! Add contributions to material macroscopic total cross section
+      material_xs % total = material_xs % total + &
+           atom_density * micro_photon_xs(i_element) % total
+
+      ! Add contributions to material macroscopic coherent cross section
+      material_xs % coherent = material_xs % coherent + &
+           atom_density * micro_photon_xs(i_element) % coherent
+
+      ! Add contributions to material macroscopic incoherent cross section
+      material_xs % incoherent = material_xs % incoherent + &
+           atom_density * micro_photon_xs(i_element) % incoherent
+
+      ! Add contributions to material macroscopic photoelectric cross section
+      material_xs % photoelectric = material_xs % photoelectric + &
+           atom_density * micro_photon_xs(i_element) % photoelectric
+
+      ! Add contributions to material macroscopic pair production cross section
+      material_xs % pair_production = material_xs % pair_production + &
+           atom_density * micro_photon_xs(i_element) % pair_production
+    end do
+
+  end subroutine calculate_photon_xs
 
 !===============================================================================
 ! FREE_MEMORY_MATERIAL deallocates global arrays defined in this module
@@ -620,4 +700,215 @@ contains
 
   end function openmc_material_set_densities
 
+  subroutine bremsstrahlung_init(this, i_material, particle)
+    class(BremsstrahlungData), intent(inout) :: this
+    integer, intent(in) :: i_material
+    integer, intent(in) :: particle
+
+    integer                 :: i, j
+    integer                 :: i_k
+    integer                 :: n, n_e, n_k
+    real(8)                 :: c
+    real(8)                 :: k, k_l, k_r
+    real(8)                 :: e, e_l, e_r
+    real(8)                 :: w, w_l, w_r
+    real(8)                 :: x, x_l, x_r
+    real(8)                 :: t
+    real(8)                 :: r
+    real(8)                 :: awr
+    real(8)                 :: beta
+    real(8)                 :: Z_eq_sq
+    real(8)                 :: atom_density
+    real(8)                 :: mass_density
+    real(8)                 :: sum_density
+    real(8), allocatable    :: stopping_power_collision(:)
+    real(8), allocatable    :: stopping_power_radiative(:)
+    real(8), allocatable    :: stopping_power(:)
+    real(8), allocatable    :: dcs(:,:)
+    real(8), allocatable    :: f(:)
+    real(8), allocatable    :: z(:)
+    logical                 :: positron_
+    type(Material), pointer :: mat
+    type(PhotonInteraction), pointer :: elm
+
+    ! Get pointer to this material
+    mat => materials(i_material)
+
+    ! Determine whether we are generating electron or positron data
+    positron_ = (particle == POSITRON)
+
+    ! Get the size of the energy grids
+    n_k = size(ttb_k_grid)
+    n_e = size(ttb_e_grid)
+
+    ! Allocate arrays for TTB data
+    allocate(this % pdf(n_e, n_e), source=ZERO)
+    allocate(this % cdf(n_e, n_e), source=ZERO)
+    allocate(this % yield(n_e))
+
+    ! Allocate temporary arrays
+    allocate(stopping_power_collision(n_e), source=ZERO)
+    allocate(stopping_power_radiative(n_e), source=ZERO)
+    allocate(stopping_power(n_e))
+    allocate(dcs(n_k, n_e), source=ZERO)
+    allocate(f(n_e))
+    allocate(z(n_e))
+
+    Z_eq_sq = ZERO
+    sum_density = ZERO
+
+    ! Calculate the molecular DCS and the molecular total stopping power using
+    ! Bragg's additivity rule.
+    ! TODO: The collision stopping power cannot be accurately calculated using
+    ! Bragg's additivity rule since the mean excitation energies and the
+    ! density effect corrections cannot simply be summed together. Bragg's
+    ! additivity rule fails especially when a higher-density compound is
+    ! composed of elements that are in lower-density form at normal temperature
+    ! and pressure (at which the NIST stopping powers are given). It will be
+    ! used to approximate the collision stopping powers for now, but should be
+    ! fixed in the future.
+    do i = 1, mat % n_nuclides
+      ! Get pointer to current element
+      elm => elements(mat % element(i))
+
+      awr = nuclides(mat % nuclide(i)) % awr
+
+      ! Get atomic density and mass density of nuclide given atom percent
+      if (mat % atom_density(1) > ZERO) then
+        atom_density = mat % atom_density(i)
+        mass_density = mat % atom_density(i) * awr
+      ! Given weight percent
+      else
+        atom_density = -mat % atom_density(i) / awr
+        mass_density = -mat % atom_density(i)
+      end if
+
+      ! Calculate the "equivalent" atomic number Zeq of the material
+      Z_eq_sq = Z_eq_sq + atom_density * elm % Z**2
+      sum_density = sum_density + atom_density
+
+      ! Accumulate material DCS
+      dcs = dcs + atom_density * elm % Z**2 * elm % dcs
+
+      ! Accumulate material collision stopping power
+      stopping_power_collision = stopping_power_collision + mass_density &
+           * MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_collision
+
+      ! Accumulate material radiative stopping power
+      stopping_power_radiative = stopping_power_radiative + mass_density &
+           * MASS_NEUTRON / N_AVOGADRO * elm % stopping_power_radiative
+    end do
+    Z_eq_sq = Z_eq_sq / sum_density
+
+    ! Calculate the positron DCS and radiative stopping power. These are
+    ! obtained by multiplying the electron DCS and radiative stopping powers by
+    ! a factor r, which is a numerical approximation of the ratio of the
+    ! radiative stopping powers for positrons and electrons. Source: F. Salvat,
+    ! J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
+    ! Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
+    ! Issy-les-Moulineaux, France (2011).
+    if (positron_) then
+      do i = 1, n_e
+        t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV))
+        r = ONE - exp(-1.2359e-1_8*t + 6.1274e-2_8*t**2 - 3.1516e-2_8*t**3 + &
+             7.7446e-3_8*t**4 - 1.0595e-3_8*t**5 + 7.0568e-5_8*t**6 - &
+             1.808e-6_8*t**7)
+        stopping_power_radiative(i) = r*stopping_power_radiative(i)
+        dcs(:,i) = r*dcs(:,i)
+      end do
+    end if
+
+    ! Total material stopping power
+    stopping_power = stopping_power_collision + stopping_power_radiative
+
+    ! Loop over photon energies
+    do i = 1, n_e - 1
+      w = ttb_e_grid(i)
+
+      ! Loop over incident particle energies
+      do j = i, n_e
+        e = ttb_e_grid(j)
+
+        ! Reduced photon energy
+        k = w / e
+
+        ! Find the lower bounding index of the reduced photon energy
+        i_k = binary_search(ttb_k_grid, n_k, k)
+
+        ! Get the interpolation bounds
+        k_l = ttb_k_grid(i_k)
+        k_r = ttb_k_grid(i_k+1)
+        x_l = dcs(i_k, j)
+        x_r = dcs(i_k+1, j)
+
+        ! Find the value of the DCS using linear interpolation in reduced
+        ! photon energy k
+        x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l)
+
+        ! Ratio of the velocity of the charged particle to the speed of light
+        beta = sqrt(e*(e + TWO*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV)
+
+        ! Compute the integrand of the PDF
+        f(j) = x / (beta**2 * stopping_power(j) * w)
+      end do
+
+      ! Number of points to integrate
+      n = n_e - i + 1
+
+      ! Integrate the PDF using cubic spline integration over the incident
+      ! particle energy
+      if (n > 2) then
+        call spline(n, ttb_e_grid(i:), f(i:), z(i:))
+
+        c = ZERO
+        do j = i, n_e - 1
+          c = c + spline_integrate(n, ttb_e_grid(i:), f(i:), z(i:), &
+               ttb_e_grid(j), ttb_e_grid(j+1))
+          this % pdf(i,j+1) = c
+        end do
+
+      ! Integrate the last two points using trapezoidal rule in log-log space
+      else
+        e_l = log(ttb_e_grid(i))
+        e_r = log(ttb_e_grid(i+1))
+        x_l = log(f(i))
+        x_r = log(f(i+1))
+
+        this % pdf(i,i+1) = HALF * (e_r - e_l) * (exp(e_l + x_l) + exp(e_r + x_r))
+      end if
+    end do
+
+    ! Loop over incident particle energies
+    do j = 2, n_e
+      ! Set last element of PDF to small non-zero value to enable log-log
+      ! interpolation
+      this % pdf(j,j) = exp(-500.0_8)
+
+      ! Loop over photon energies
+      c = ZERO
+      do i = 1, j - 1
+        ! Integrate the CDF from the PDF using the trapezoidal rule in log-log
+        ! space
+        w_l = log(ttb_e_grid(i))
+        w_r = log(ttb_e_grid(i+1))
+        x_l = log(this % pdf(i,j))
+        x_r = log(this % pdf(i+1,j))
+
+        c = c + HALF * (w_r - w_l) * (exp(w_l + x_l) + exp(w_r + x_r))
+        this % cdf(i+1,j) = c
+      end do
+
+      ! Set photon number yield
+      this % yield(j) = c
+    end do
+
+    ! Use logarithm of number yield since it is log-log interpolated
+    where (this % yield > ZERO)
+      this % yield = log(this % yield)
+    elsewhere
+      this % yield = -500.0_8
+    end where
+
+  end subroutine bremsstrahlung_init
+
 end module material_header

src/math.F90

@@ -18,11 +18,9 @@ module math
   public :: faddeeva
   public :: w_derivative
   public :: broaden_wmp_polynomials
-
-!===============================================================================
-! FADDEEVA_W evaluates the scaled complementary error function.  This
-! interfaces with the MIT C library
-!===============================================================================
+  public :: spline
+  public :: spline_interpolate
+  public :: spline_integrate
 
   interface
 
@@ -105,6 +103,40 @@ module math
       real(C_DOUBLE), intent(inout) :: factors(n)
     end subroutine broaden_wmp_polynomials
 
+    subroutine spline(n, x, y, z) bind(C, name='spline_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(in)        :: x(n)
+      real(C_DOUBLE), intent(in)        :: y(n)
+      real(C_DOUBLE), intent(in)        :: z(n)
+    end subroutine spline
+
+    function spline_interpolate(n, x, y, z, xint) &
+         bind(C, name='spline_interpolate_c') result(yint)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(in)        :: x(n)
+      real(C_DOUBLE), intent(in)        :: y(n)
+      real(C_DOUBLE), intent(in)        :: z(n)
+      real(C_DOUBLE), value, intent(in) :: xint
+      real(C_DOUBLE)                    :: yint
+    end function spline_interpolate
+
+    function spline_integrate(n, x, y, z, xa, xb) &
+         bind(C, name='spline_integrate_c') result(s)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(in)        :: x(n)
+      real(C_DOUBLE), intent(in)        :: y(n)
+      real(C_DOUBLE), intent(in)        :: z(n)
+      real(C_DOUBLE), value, intent(in) :: xa
+      real(C_DOUBLE), value, intent(in) :: xb
+      real(C_DOUBLE)                    :: s
+    end function spline_integrate
+
     function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w)
       use ISO_C_BINDING
       implicit none

src/math_functions.cpp

@@ -690,4 +690,96 @@ void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
   }
 }
 
+
+void spline_c(int n, const double x[], const double y[], double z[])
+{
+  double c_new[n-1];
+
+  // Set natural boundary conditions
+  c_new[0] = 0.0;
+  z[0] = 0.0;
+  z[n-1] = 0.0;
+
+  // Solve using tridiagonal matrix algorithm; first do forward sweep
+  for (int i = 1; i < n - 1; i++) {
+    double a = x[i] - x[i-1];
+    double c = x[i+1] - x[i];
+    double b = 2.0*(a + c);
+    double d = 6.0*((y[i+1] - y[i])/c - (y[i] - y[i-1])/a);
+
+    c_new[i] = c/(b - a*c_new[i-1]);
+    z[i] = (d - a*z[i-1])/(b - a*c_new[i-1]);
+  }
+
+  // Back substitution
+  for (int i = n - 2; i >= 0; i--) {
+    z[i] = z[i] - c_new[i]*z[i+1];
+  }
+}
+
+
+double spline_interpolate_c(int n, const double x[], const double y[],
+                            const double z[], double xint)
+{
+  // Find the lower bounding index in x of xint
+  int i = n - 1;
+  while (--i) {
+    if (xint >= x[i]) break;
+  }
+
+  double h = x[i+1] - x[i];
+  double r = xint - x[i];
+
+  // Compute the coefficients
+  double b = (y[i+1] - y[i])/h - (h/6.0)*(z[i+1] + 2.0*z[i]);
+  double c = z[i]/2.0;
+  double d = (z[i+1] - z[i])/(h*6.0);
+
+  return y[i] + b*r + c*r*r + d*r*r*r;
+}
+
+
+double spline_integrate_c(int n, const double x[], const double y[],
+                          const double z[], double xa, double xb)
+{
+  // Find the lower bounding index in x of the lower limit of integration.
+  int ia = n - 1;
+  while (--ia) {
+    if (xa >= x[ia]) break;
+  }
+
+  // Find the lower bounding index in x of the upper limit of integration.
+  int ib = n - 1;
+  while (--ib) {
+    if (xb >= x[ib]) break;
+  }
+
+  // Evaluate the integral
+  double s = 0.0;
+  for (int i = ia; i <= ib; i++) {
+    double h = x[i+1] - x[i];
+
+    // Compute the coefficients
+    double b = (y[i+1] - y[i])/h - (h/6.0)*(z[i+1] + 2.0*z[i]);
+    double c = z[i]/2.0;
+    double d = (z[i+1] - z[i])/(h*6.0);
+
+    // Subtract the integral from x[ia] to xa
+    if (i == ia) {
+      double r = xa - x[ia];
+      s = s - (y[i]*r + b/2.0*r*r + c/3.0*r*r*r + d/4.0*r*r*r*r);
+    }
+
+    // Integrate from x[ib] to xb in final interval
+    if (i == ib) {
+      h = xb - x[ib];
+    }
+
+    // Accumulate the integral
+    s = s + y[i]*h + b/2.0*h*h + c/3.0*h*h*h + d/4.0*h*h*h*h;
+  }
+
+  return s;
+}
+
 } // namespace openmc

src/math_functions.h

@@ -149,5 +149,57 @@ extern "C" double watt_spectrum_c(double a, double b);
 extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
                                           double factors[]);
 
+//==============================================================================
+//! Constructs a natural cubic spline.
+//!
+//! Given a tabulated function y_i = f(x_i), this computes the second
+//! derivative of the interpolating function at each x_i, which can then be
+//! used in any subsequent calls to spline_interpolate or spline_integrate for
+//! the same set of x and y values.
+//!
+//! @param n       Number of points
+//! @param x       Values of the independent variable, which must be strictly
+//!   increasing.
+//! @param y       Values of the dependent variable.
+//! @param[out] z  The second derivative of the interpolating function at each
+//!   value of x.
+//==============================================================================
+
+extern "C" void spline_c(int n, const double x[], const double y[], double z[]);
+
+//==============================================================================
+//! Determine the cubic spline interpolated y-value for a given x-value.
+//!
+//! @param n     Number of points
+//! @param x     Values of the independent variable, which must be strictly
+//!   increasing.
+//! @param y     Values of the dependent variable.
+//! @param z     The second derivative of the interpolating function at each
+//!   value of x.
+//! @param xint  Point at which to evaluate the cubic spline polynomial
+//! @result      Interpolated value
+//==============================================================================
+
+extern "C" double spline_interpolate_c(int n, const double x[], const double y[],
+                                       const double z[], double xint);
+
+//==============================================================================
+//! Evaluate the definite integral of the interpolating cubic spline between
+//! the given endpoints.
+//!
+//! @param n   Number of points
+//! @param x   Values of the independent variable, which must be strictly
+//!   increasing.
+//! @param y   Values of the dependent variable.
+//! @param z   The second derivative of the interpolating function at each
+//!   value of x.
+//! @param xa  Lower limit of integration
+//! @param xb  Upper limit of integration
+//! @result    Integral
+//==============================================================================
+
+extern "C" double spline_integrate_c(int n, const double x[], const double y[],
+                                     const double z[], double xa, double xb);
+
 } // namespace openmc
-#endif // MATH_FUNCTIONS_H
+#endif // MATH_FUNCTIONS_H

src/nuclide_header.F90

@@ -42,14 +42,15 @@ module nuclide_header
 
   ! Positions for first dimension of Nuclide % xs
   integer, parameter :: &
-       XS_TOTAL      = 1, &
-       XS_ABSORPTION = 2, &
-       XS_FISSION    = 3, &
-       XS_NU_FISSION = 4
+       XS_TOTAL       = 1, &
+       XS_ABSORPTION  = 2, &
+       XS_FISSION     = 3, &
+       XS_NU_FISSION  = 4, &
+       XS_PHOTON_PROD = 5
 
-  ! The array within SumXS is of shape (4, n_energy) where the first dimension
+  ! The array within SumXS is of shape (5, n_energy) where the first dimension
   ! corresponds to the following values: 1) total, 2) absorption (MT > 100), 3)
-  ! fission, 4) neutron production
+  ! fission, 4) neutron production, 5) photon production
   type SumXS
     real(8), allocatable :: value(:,:)
   end type SumXS
@@ -104,8 +105,8 @@ module nuclide_header
     integer :: reaction_index(891)
 
     ! Fission energy release
-    class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas
-    class(Function1D), allocatable :: fission_q_recov  ! neutrons, gammas, betas
+    class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas
+    class(Function1D), allocatable :: fission_q_recov  ! fragments, neutrons, gammas, betas
 
   contains
     procedure :: assign_0K_elastic_scattering
@@ -138,6 +139,7 @@ module nuclide_header
                                 !   averaged over bound and non-bound nuclei
     real(8) :: thermal          ! Bound thermal elastic & inelastic scattering
     real(8) :: thermal_elastic  ! Bound thermal elastic scattering
+    real(8) :: photon_prod      ! microscopic photon production xs
 
     ! Cross sections for depletion reactions (note that these are not stored in
     ! macroscopic cache)
@@ -169,6 +171,13 @@ module nuclide_header
     real(C_DOUBLE) :: absorption    ! macroscopic absorption xs
     real(C_DOUBLE) :: fission       ! macroscopic fission xs
     real(C_DOUBLE) :: nu_fission    ! macroscopic production xs
+    real(C_DOUBLE) :: photon_prod   ! macroscopic photon production xs
+
+    ! Photon cross sections
+    real(C_DOUBLE) :: coherent        ! macroscopic coherent xs
+    real(C_DOUBLE) :: incoherent      ! macroscopic incoherent xs
+    real(C_DOUBLE) :: photoelectric   ! macroscopic photoelectric xs
+    real(C_DOUBLE) :: pair_production ! macroscopic pair production xs
   end type MaterialMacroXS
 
 !===============================================================================
@@ -196,8 +205,8 @@ module nuclide_header
 !$omp threadprivate(micro_xs, material_xs)
 
   ! Minimum/maximum energies
-  real(8) :: energy_min_neutron = ZERO
-  real(8) :: energy_max_neutron = INFINITY
+  real(8) :: energy_min(2) = [ZERO, ZERO]
+  real(8) :: energy_max(2) = [INFINITY, INFINITY]
 
 contains
 
@@ -549,34 +558,36 @@ contains
     if (object_exists(group_id, 'fission_energy_release')) then
       fer_group = open_group(group_id, 'fission_energy_release')
 
-      ! Check to see if this is polynomial or tabulated data
+      ! Q-PROMPT
       fer_dset = open_dataset(fer_group, 'q_prompt')
       call read_attribute(temp_str, fer_dset, 'type')
       if (temp_str == 'Polynomial') then
-        ! Read the prompt Q-value
         allocate(Polynomial :: this % fission_q_prompt)
         call this % fission_q_prompt % from_hdf5(fer_dset)
         call close_dataset(fer_dset)
-
-        ! Read the recoverable energy Q-value
-        allocate(Polynomial :: this % fission_q_recov)
-        fer_dset = open_dataset(fer_group, 'q_recoverable')
-        call this % fission_q_recov % from_hdf5(fer_dset)
-        call close_dataset(fer_dset)
       else if (temp_str == 'Tabulated1D') then
-        ! Read the prompt Q-value
         allocate(Tabulated1D :: this % fission_q_prompt)
         call this % fission_q_prompt % from_hdf5(fer_dset)
         call close_dataset(fer_dset)
+      else
+        call fatal_error('Unrecognized fission prompt energy release format.')
+      end if
 
-        ! Read the recoverable energy Q-value
+      ! Q-RECOV
+      fer_dset = open_dataset(fer_group, 'q_recoverable')
+      call read_attribute(temp_str, fer_dset, 'type')
+      if (temp_str == 'Polynomial') then
+        allocate(Polynomial :: this % fission_q_recov)
+        call this % fission_q_recov % from_hdf5(fer_dset)
+        call close_dataset(fer_dset)
+      else if (temp_str == 'Tabulated1D') then
         allocate(Tabulated1D :: this % fission_q_recov)
-        fer_dset = open_dataset(fer_group, 'q_recoverable')
         call this % fission_q_recov % from_hdf5(fer_dset)
         call close_dataset(fer_dset)
       else
-        call fatal_error('Unrecognized fission energy release format.')
+        call fatal_error('Unrecognized fission recoverable energy release format.')
       end if
+
       call close_group(fer_group)
     end if
 
@@ -591,7 +602,7 @@ contains
   subroutine nuclide_create_derived(this)
     class(Nuclide), intent(inout) :: this
 
-    integer :: i, j, k
+    integer :: i, j, k, l
     integer :: t
     integer :: m
     integer :: n
@@ -606,7 +617,7 @@ contains
     do i = 1, n_temperature
       ! Allocate and initialize derived cross sections
       n_grid = size(this % grid(i) % energy)
-      allocate(this % xs(i) % value(4,n_grid))
+      allocate(this % xs(i) % value(5,n_grid))
       this % xs(i) % value(:,:) = ZERO
     end do
 
@@ -638,6 +649,18 @@ contains
           this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % &
                value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value
 
+          ! Calculate photon production cross section
+          do k = 1, size(rx % products)
+            if (rx % products(k) % particle == PHOTON) then
+              do l = 1, n
+                this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = &
+                     this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + &
+                     rx % xs(t) % value(l) * rx % products(k) % &
+                     yield % evaluate(this % grid(t) % energy(l+j-1))
+              end do
+            end if
+          end do
+
           ! Add contribution to absorption cross section
           if (is_disappearance(rx % MT)) then
             this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
@@ -664,10 +687,6 @@ contains
             this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
                  value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
 
-            ! If total fission reaction is present, there's no need to store the
-            ! reaction cross-section since it was copied to this % fission
-            if (rx % MT == N_FISSION) deallocate(rx % xs(t) % value)
-
             ! Keep track of this reaction for easy searching later
             if (t == 1) then
               i_fission = i_fission + 1
@@ -975,6 +994,10 @@ contains
           micro_xs % fission         = ZERO
           micro_xs % nu_fission      = ZERO
         end if
+
+        ! Calculate microscopic nuclide photon production cross section
+        micro_xs % photon_prod = (ONE - f) * xs % value(XS_PHOTON_PROD,i_grid) &
+             + f * xs % value(XS_PHOTON_PROD,i_grid + 1)
       end associate
 
       ! Depletion-related reactions
@@ -1695,8 +1718,8 @@ contains
         if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering()
 
         ! Initialize nuclide grid
-        call nuclides(n) % init_grid(energy_min_neutron, &
-             energy_max_neutron, n_log_bins)
+        call nuclides(n) % init_grid(energy_min(NEUTRON), &
+             energy_max(NEUTRON), n_log_bins)
       else
         err = E_DATA
         call set_errmsg("Nuclide '" // trim(name_) // "' is not present &

src/output.F90

@@ -547,7 +547,7 @@ contains
 
     ! format for write statements
 100 format (1X,A,T36,"= ",ES11.4," seconds")
-101 format (1X,A,T36,"=  ",A," neutrons/second")
+101 format (1X,A,T36,"=  ",A," particles/second")
 
   end subroutine print_runtime
 

src/particle_header.F90

@@ -3,7 +3,6 @@ module particle_header
   use bank_header,     only: Bank, source_bank
   use constants
   use error,           only: fatal_error, warning
-  use geometry_header, only: root_universe
   use hdf5_interface
   use settings
   use simulation_header
@@ -159,9 +158,10 @@ contains
 ! the secondary bank and increments the number of sites in the secondary bank.
 !===============================================================================
 
-  subroutine create_secondary(this, uvw, type, run_CE)
+  subroutine create_secondary(this, uvw, E, type, run_CE)
     class(Particle), intent(inout) :: this
     real(8),         intent(in)    :: uvw(3)
+    real(8),         intent(in)    :: E
     integer,         intent(in)    :: type
     logical,         intent(in)    :: run_CE
 
@@ -174,14 +174,15 @@ contains
     end if
 
     n = this % n_secondary + 1
-    this % secondary_bank(n) % wgt    = this % wgt
+    this % secondary_bank(n) % particle = type
+    this % secondary_bank(n) % wgt = this % wgt
     this % secondary_bank(n) % xyz(:) = this % coord(1) % xyz
     this % secondary_bank(n) % uvw(:) = uvw
-    this % n_secondary = n
-    this % secondary_bank(this % n_secondary) % E = this % E
+    this % secondary_bank(n) % E = E
     if (.not. run_CE) then
-      this % secondary_bank(this % n_secondary) % E = real(this % g, 8)
+      this % secondary_bank(n) % E = real(this % g, 8)
     end if
+    this % n_secondary = n
 
   end subroutine create_secondary
 
@@ -219,7 +220,7 @@ contains
     this % g = NONE
 
     ! Set up base level coordinates
-    this % coord(1) % universe = root_universe
+    this % coord(1) % universe = NONE
     this % n_coord = 1
     this % last_n_coord = 1
 
@@ -241,6 +242,7 @@ contains
     call this % initialize()
 
     ! copy attributes from source bank site
+    this % type             = src % particle
     this % wgt              = src % wgt
     this % last_wgt         = src % wgt
     this % coord(1) % xyz   = src % xyz
@@ -335,6 +337,7 @@ contains
         call write_dataset(file_id, 'run_mode', 'particle restart')
       end select
       call write_dataset(file_id, 'id', this % id)
+      call write_dataset(file_id, 'type', this % type)
       call write_dataset(file_id, 'weight', src % wgt)
       call write_dataset(file_id, 'energy', src % E)
       call write_dataset(file_id, 'xyz', src % xyz)

src/particle_restart.F90

@@ -100,6 +100,7 @@ contains
       previous_run_mode = MODE_FIXEDSOURCE
     end select
     call read_dataset(p % id, file_id, 'id')
+    call read_dataset(p % type, file_id, 'type')
     call read_dataset(p % wgt, file_id, 'weight')
     call read_dataset(p % E, file_id, 'energy')
     call read_dataset(p % coord(1) % xyz, file_id, 'xyz')

src/photon_header.F90

@@ -0,0 +1,506 @@
+module photon_header
+
+  use algorithm,        only: binary_search
+  use constants
+  use dict_header,      only: DictIntInt, DictCharInt
+  use endf_header,      only: Tabulated1D
+  use hdf5_interface
+  use nuclide_header,   only: nuclides
+  use settings
+
+  real(8), allocatable :: compton_profile_pz(:)
+  real(8), allocatable :: ttb_e_grid(:) ! energy T of incident electron
+  real(8), allocatable :: ttb_k_grid(:) ! reduced energy W/T of emitted photon
+
+  type ElectronSubshell
+    integer :: index_subshell  ! index in SUBSHELLS
+    integer :: threshold
+    real(8) :: n_electrons
+    real(8) :: binding_energy
+    real(8), allocatable :: cross_section(:)
+
+    ! Transition data
+    integer :: n_transitions
+    integer, allocatable :: transition_subshells(:,:)
+    real(8), allocatable :: transition_energy(:)
+    real(8), allocatable :: transition_probability(:)
+  end type ElectronSubshell
+
+  type PhotonInteraction
+    character(3) :: name  ! atomic symbol, e.g. 'Zr'
+    integer      :: Z     ! atomic number
+
+    ! Microscopic cross sections
+    real(8), allocatable :: energy(:)
+    real(8), allocatable :: coherent(:)
+    real(8), allocatable :: incoherent(:)
+    real(8), allocatable :: photoelectric_total(:)
+    real(8), allocatable :: pair_production_total(:)
+    real(8), allocatable :: pair_production_electron(:)
+    real(8), allocatable :: pair_production_nuclear(:)
+
+    ! Form factors
+    type(Tabulated1D) :: incoherent_form_factor
+    type(Tabulated1D) :: coherent_int_form_factor
+    type(Tabulated1D) :: coherent_anomalous_real
+    type(Tabulated1D) :: coherent_anomalous_imag
+
+    ! Photoionization and atomic relaxation data
+    type(DictIntInt) :: shell_dict ! Given a shell designator, e.g. 3, this
+                                   ! dictionary gives an index in shells(:)
+    type(ElectronSubshell), allocatable :: shells(:)
+
+    ! Compton profile data
+    real(8), allocatable :: profile_pdf(:,:)
+    real(8), allocatable :: profile_cdf(:,:)
+    real(8), allocatable :: binding_energy(:)
+    real(8), allocatable :: electron_pdf(:)
+
+    ! Stopping power data
+    real(8) :: I ! mean excitation energy
+    real(8), allocatable :: stopping_power_collision(:)
+    real(8), allocatable :: stopping_power_radiative(:)
+
+    ! Bremsstrahlung scaled DCS
+    real(8), allocatable :: dcs(:,:)
+
+  contains
+    procedure :: from_hdf5 => photon_from_hdf5
+    procedure :: calculate_xs => photon_calculate_xs
+  end type PhotonInteraction
+
+  type BremsstrahlungData
+    real(8), allocatable :: pdf(:,:) ! Bremsstrahlung energy PDF
+    real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
+    real(8), allocatable :: yield(:) ! Photon number yield
+  end type BremsstrahlungData
+
+  type Bremsstrahlung
+    type(BremsstrahlungData) :: electron
+    type(BremsstrahlungData) :: positron
+  end type Bremsstrahlung
+
+  type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections
+  integer :: n_elements       ! Number of photon cross section tables
+
+  type(DictCharInt) :: element_dict
+
+  type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung data
+
+!===============================================================================
+! ELEMENTMICROXS contains cached microscopic photon cross sections for a
+! particular element at the current energy
+!===============================================================================
+
+  type ElementMicroXS
+    integer :: index_grid      ! index on element energy grid
+    real(8) :: last_E = ZERO   ! last evaluated energy
+    real(8) :: interp_factor   ! interpolation factor on energy grid
+    real(8) :: total           ! microscropic total photon xs
+    real(8) :: coherent        ! microscopic coherent xs
+    real(8) :: incoherent      ! microscopic incoherent xs
+    real(8) :: photoelectric   ! microscopic photoelectric xs
+    real(8) :: pair_production ! microscopic pair production xs
+  end type ElementMicroXS
+
+  type(ElementMicroXS), allocatable :: micro_photon_xs(:) ! Cache for each element
+!$omp threadprivate(micro_photon_xs)
+
+contains
+
+  subroutine photon_from_hdf5(this, group_id)
+    class(PhotonInteraction), intent(inout) :: this
+    integer(HID_T), intent(in) :: group_id
+
+    integer          :: i, j
+    integer(HID_T)   :: rgroup, tgroup
+    integer(HID_T)   :: dset_id
+    integer(HSIZE_T) :: dims(1), dims2(2)
+    integer          :: n_energy
+    integer          :: n_shell
+    integer          :: n_profile
+    integer          :: n_transition
+    integer          :: n_k
+    integer          :: n_e
+    character(3), allocatable :: designators(:)
+    real(8)          :: c
+    real(8)          :: f
+    real(8)          :: y
+    real(8), allocatable :: electron_energy(:)
+    real(8), allocatable :: matrix(:,:)
+    real(8), allocatable :: dcs(:,:)
+
+    ! Get name of nuclide from group
+    this % name = get_name(group_id)
+
+    ! Get rid of leading '/'
+    this % name = trim(this % name(2:))
+
+    ! Get atomic number
+    call read_attribute(this % Z, group_id, 'Z')
+
+    ! Determine number of energies and read energy grid
+    dset_id = open_dataset(group_id, 'energy')
+    call get_shape(dset_id, dims)
+    n_energy = int(dims(1), 4)
+    allocate(this % energy(dims(1)))
+    call read_dataset(this % energy, dset_id)
+    call close_dataset(dset_id)
+
+    ! Allocate arrays
+    allocate(this % coherent(n_energy))
+    allocate(this % incoherent(n_energy))
+    allocate(this % pair_production_total(n_energy))
+    allocate(this % pair_production_nuclear(n_energy))
+    allocate(this % pair_production_electron(n_energy))
+    allocate(this % photoelectric_total(n_energy))
+
+    ! Read coherent scattering
+    rgroup = open_group(group_id, 'coherent')
+    call read_dataset(this % coherent, rgroup, 'xs')
+
+    dset_id = open_dataset(rgroup, 'integrated_scattering_factor')
+    call this % coherent_int_form_factor % from_hdf5(dset_id)
+    call close_dataset(dset_id)
+
+    dset_id = open_dataset(rgroup, 'anomalous_real')
+    call this % coherent_anomalous_real % from_hdf5(dset_id)
+    call close_dataset(dset_id)
+
+    dset_id = open_dataset(rgroup, 'anomalous_imag')
+    call this % coherent_anomalous_imag % from_hdf5(dset_id)
+    call close_dataset(dset_id)
+    call close_group(rgroup)
+
+    ! Read incoherent scattering
+    rgroup = open_group(group_id, 'incoherent')
+    call read_dataset(this % incoherent, rgroup, 'xs')
+    dset_id = open_dataset(rgroup, 'scattering_factor')
+    call this % incoherent_form_factor % from_hdf5(dset_id)
+    call close_dataset(dset_id)
+    call close_group(rgroup)
+
+    ! Read pair production
+    rgroup = open_group(group_id, 'pair_production_electron')
+    call read_dataset(this % pair_production_electron, rgroup, 'xs')
+    call close_group(rgroup)
+
+    ! Read pair production
+    if (object_exists(group_id, 'pair_production_nuclear')) then
+      rgroup = open_group(group_id, 'pair_production_nuclear')
+      call read_dataset(this % pair_production_nuclear, rgroup, 'xs')
+      call close_group(rgroup)
+    else
+      this % pair_production_nuclear(:) = ZERO
+    end if
+
+    ! Read photoelectric
+    rgroup = open_group(group_id, 'photoelectric')
+    call read_dataset(this % photoelectric_total, rgroup, 'xs')
+    call close_group(rgroup)
+
+    ! Read subshell photoionization cross section and atomic relaxation data
+    rgroup = open_group(group_id, 'subshells')
+    call read_attribute(designators, rgroup, 'designators')
+    n_shell = size(designators)
+    allocate(this % shells(n_shell))
+    do i = 1, n_shell
+      ! Create mapping from designator to index
+      do j = 1, size(SUBSHELLS)
+        if (designators(i) == SUBSHELLS(j)) then
+          call this % shell_dict % set(j, i)
+          this % shells(i) % index_subshell = j
+          exit
+        end if
+      end do
+
+      ! Read binding energy and number of electrons
+      tgroup = open_group(rgroup, trim(designators(i)))
+      call read_attribute(this % shells(i) % binding_energy, tgroup, &
+           'binding_energy')
+      call read_attribute(this % shells(i) % n_electrons, tgroup, &
+           'num_electrons')
+
+      ! Read subshell cross section
+      dset_id = open_dataset(tgroup, 'xs')
+      call read_attribute(j, dset_id, 'threshold_idx')
+      this % shells(i) % threshold = j
+      allocate(this % shells(i) % cross_section(n_energy - j))
+      call read_dataset(this % shells(i) % cross_section, dset_id)
+      call close_dataset(dset_id)
+      where (this % shells(i) % cross_section > ZERO)
+        this % shells(i) % cross_section = log(this % shells(i) % cross_section)
+      elsewhere
+        this % shells(i) % cross_section = -500.0_8
+      end where
+
+      if (object_exists(tgroup, 'transitions')) then
+        dset_id = open_dataset(tgroup, 'transitions')
+        call get_shape(dset_id, dims2)
+        n_transition = int(dims2(2), 4)
+        this % shells(i) % n_transitions = n_transition
+        if (n_transition > 0) then
+          allocate(this % shells(i) % transition_subshells(2, n_transition))
+          allocate(this % shells(i) % transition_energy(n_transition))
+          allocate(this % shells(i) % transition_probability(n_transition))
+
+          allocate(matrix(dims2(1), dims2(2)))
+          call read_dataset(matrix, dset_id)
+
+          this % shells(i) % transition_subshells(:,:) = int(matrix(1:2, :), 4)
+          this % shells(i) % transition_energy(:) = matrix(3, :)
+          this % shells(i) % transition_probability(:) = matrix(4, :) &
+               / sum(matrix(4, :))
+          deallocate(matrix)
+        end if
+        call close_dataset(dset_id)
+      else
+        this % shells(i) % n_transitions = 0
+      end if
+      call close_group(tgroup)
+    end do
+    call close_group(rgroup)
+    deallocate(designators)
+
+    ! Determine number of electron shells
+    rgroup = open_group(group_id, 'compton_profiles')
+
+    ! Determine number of shells
+    dset_id = open_dataset(rgroup, 'num_electrons')
+    call get_shape(dset_id, dims)
+    n_shell = int(dims(1), 4)
+
+    ! Read electron shell PDF and binding energies
+    allocate(this % electron_pdf(n_shell), this % binding_energy(n_shell))
+    call read_dataset(this % electron_pdf, dset_id)
+    call close_dataset(dset_id)
+    call read_dataset(this % binding_energy, rgroup, 'binding_energy')
+    this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf)
+
+    ! Read Compton profiles
+    dset_id = open_dataset(rgroup, 'J')
+    call get_shape(dset_id, dims2)
+    n_profile = int(dims2(1), 4)
+    allocate(this % profile_pdf(n_profile, n_shell))
+    call read_dataset(this % profile_pdf, dset_id)
+    call close_dataset(dset_id)
+
+    ! Get Compton profile momentum grid
+    if (.not. allocated(compton_profile_pz)) then
+      allocate(compton_profile_pz(n_profile))
+      call read_dataset(compton_profile_pz, rgroup, 'pz')
+    end if
+    call close_group(rgroup)
+
+    ! Create Compton profile CDF
+    allocate(this % profile_cdf(n_profile, n_shell))
+    do i = 1, n_shell
+      c = ZERO
+      this % profile_cdf(1,i) = ZERO
+      do j = 1, n_profile - 1
+        c = c + HALF*(compton_profile_pz(j+1) - compton_profile_pz(j)) * &
+             (this%profile_pdf(j,i) + this%profile_pdf(j+1,i))
+        this % profile_cdf(j+1,i) = c
+      end do
+    end do
+
+    ! Calculate total pair production
+    this % pair_production_total(:) = this % pair_production_nuclear + &
+         this % pair_production_electron
+
+    if (electron_treatment == ELECTRON_TTB) then
+      ! Read bremsstrahlung scaled DCS
+      rgroup = open_group(group_id, 'bremsstrahlung')
+      dset_id = open_dataset(rgroup, 'dcs')
+      call get_shape(dset_id, dims2)
+      n_k = int(dims2(1), 4)
+      n_e = int(dims2(2), 4)
+      allocate(this % dcs(n_k, n_e))
+      call read_dataset(this % dcs, dset_id)
+      call close_dataset(dset_id)
+
+      ! Get energy grids used for bremsstrahlung DCS and for stopping powers
+      allocate(electron_energy(n_e))
+      call read_dataset(electron_energy, rgroup, 'electron_energy')
+      if (.not. allocated(ttb_k_grid)) then
+        allocate(ttb_k_grid(n_k))
+        call read_dataset(ttb_k_grid, rgroup, 'photon_energy')
+      end if
+      call close_group(rgroup)
+
+      ! Read stopping power data
+      if (this % Z < 99) then
+        rgroup = open_group(group_id, 'stopping_powers')
+        allocate(this % stopping_power_collision(n_e))
+        allocate(this % stopping_power_radiative(n_e))
+        call read_dataset(this % stopping_power_collision, rgroup, 's_collision')
+        call read_dataset(this % stopping_power_radiative, rgroup, 's_radiative')
+        call read_attribute(this % I, rgroup, 'I')
+        call close_group(rgroup)
+      end if
+
+      ! Truncate the bremsstrahlung data at the cutoff energy
+      if (energy_cutoff(PHOTON) > electron_energy(1)) then
+        i_grid = binary_search(electron_energy, n_e, energy_cutoff(PHOTON))
+
+        ! calculate interpolation factor
+        f = (log(energy_cutoff(PHOTON)) - log(electron_energy(i_grid))) / &
+             (log(electron_energy(i_grid+1)) - log(electron_energy(i_grid)))
+
+        ! Interpolate collision stopping power at the cutoff energy and
+        ! truncate
+        y = exp(log(this % stopping_power_collision(i_grid)) + &
+             f*(log(this % stopping_power_collision(i_grid+1)) - &
+             log(this % stopping_power_collision(i_grid))))
+        this % stopping_power_collision = &
+             [y, this % stopping_power_collision(i_grid+1:n_e)]
+
+        ! Interpolate radiative stopping power at the cutoff energy and
+        ! truncate
+        y = exp(log(this % stopping_power_radiative(i_grid)) + &
+             f*(log(this % stopping_power_radiative(i_grid+1)) - &
+             log(this % stopping_power_radiative(i_grid))))
+        this % stopping_power_radiative = &
+             [y, this % stopping_power_radiative(i_grid+1:n_e)]
+
+        ! Interpolate bremsstrahlung DCS at the cutoff energy and truncate
+        allocate(dcs(n_k, n_e-i_grid+1))
+        do i = 1, n_k
+          y = exp(log(this % dcs(i,i_grid)) + &
+               f*(log(this % dcs(i,i_grid+1)) - log(this % dcs(i,i_grid))))
+          dcs(i,:) = [y, this % dcs(i,i_grid+1:n_e)]
+        end do
+        call move_alloc(dcs, this % dcs)
+
+        electron_energy = [energy_cutoff(PHOTON), electron_energy(i_grid+1:n_e)]
+      end if
+
+      ! Set incident particle energy grid
+      if (.not. allocated(ttb_e_grid)) then
+        call move_alloc(electron_energy, ttb_e_grid)
+      end if
+    end if
+
+    ! Take logarithm of energies and cross sections since they are log-log
+    ! interpolated
+    this % energy = log(this % energy)
+
+    where (this % coherent > ZERO)
+      this % coherent = log(this % coherent)
+    elsewhere
+      this % coherent = -500.0_8
+    end where
+
+    where (this % incoherent > ZERO)
+      this % incoherent = log(this % incoherent)
+    elsewhere
+      this % incoherent = -500.0_8
+    end where
+
+    where (this % photoelectric_total > ZERO)
+      this % photoelectric_total = log(this % photoelectric_total)
+    elsewhere
+      this % photoelectric_total = -500.0_8
+    end where
+
+    where (this % pair_production_total > ZERO)
+      this % pair_production_total = log(this % pair_production_total)
+    elsewhere
+      this % pair_production_total = -500.0_8
+    end where
+
+  end subroutine photon_from_hdf5
+
+!===============================================================================
+! CALCULATE_ELEMENT_XS determines microscopic photon cross sections for an
+! element of a given index in the elements array at the energy of the given
+! particle
+!===============================================================================
+
+  subroutine photon_calculate_xs(this, E, xs)
+    class(PhotonInteraction), intent(in)    :: this ! index into elements array
+    real(8),                  intent(in)    :: E ! energy
+    type(ElementMicroXS),     intent(inout) :: xs
+
+    integer :: i_grid  ! index on element energy grid
+    integer :: i_shell ! index in subshells
+    integer :: i_start ! threshold index
+    integer :: n_grid  ! number of grid points
+    real(8) :: f       ! interp factor on element energy grid
+    real(8) :: log_E   ! logarithm of the energy
+
+    ! Perform binary search on the element energy grid in order to determine
+    ! which points to interpolate between
+    n_grid = size(this % energy)
+    log_E = log(E)
+    if (log_E <= this % energy(1)) then
+      i_grid = 1
+    elseif (log_E > this % energy(n_grid)) then
+      i_grid = n_grid - 1
+    else
+      i_grid = binary_search(this % energy, n_grid, log_E)
+    end if
+
+    ! check for case where two energy points are the same
+    if (this % energy(i_grid) == this % energy(i_grid+1)) i_grid = i_grid + 1
+
+    ! calculate interpolation factor
+    f = (log_E - this % energy(i_grid)) / &
+         (this % energy(i_grid+1) - this % energy(i_grid))
+
+    xs % index_grid    = i_grid
+    xs % interp_factor = f
+
+    ! Calculate microscopic coherent cross section
+    xs % coherent = exp(this % coherent(i_grid) + f * &
+         (this % coherent(i_grid+1) - this % coherent(i_grid)))
+
+    ! Calculate microscopic incoherent cross section
+    xs % incoherent = exp(this % incoherent(i_grid) + &
+         f*(this % incoherent(i_grid+1) - this % incoherent(i_grid)))
+
+    ! Calculate microscopic photoelectric cross section
+    xs % photoelectric = ZERO
+    do i_shell = 1, size(this % shells)
+      ! Check threshold of reaction
+      i_start = this % shells(i_shell) % threshold
+      if (i_grid <= i_start) cycle
+
+      ! Evaluation subshell photoionization cross section
+      xs % photoelectric = xs % photoelectric + &
+           exp(this % shells(i_shell) % cross_section(i_grid-i_start) + &
+           f*(this % shells(i_shell) % cross_section(i_grid+1-i_start) - &
+           this % shells(i_shell) % cross_section(i_grid-i_start)))
+    end do
+
+    ! Calculate microscopic pair production cross section
+    xs % pair_production = exp(&
+         this % pair_production_total(i_grid) + f*(&
+         this % pair_production_total(i_grid+1) - &
+         this % pair_production_total(i_grid)))
+
+    ! Calculate microscopic total cross section
+    xs % total = xs % coherent + xs % incoherent + xs % photoelectric + &
+         xs % pair_production
+
+    xs % last_E = E
+
+  end subroutine photon_calculate_xs
+
+!===============================================================================
+! FREE_MEMORY_PHOTON deallocates/resets global variables in this module
+!===============================================================================
+
+  subroutine free_memory_photon()
+    ! Deallocate photon cross section data
+    if (allocated(elements)) deallocate(elements)
+    if (allocated(compton_profile_pz)) deallocate(compton_profile_pz)
+    n_elements = 0
+    call element_dict % clear()
+
+    ! Clear TTB-related arrays
+    if (allocated(ttb_e_grid)) deallocate(ttb_e_grid)
+    if (allocated(ttb)) deallocate(ttb)
+  end subroutine free_memory_photon
+
+end module photon_header

src/photon_physics.F90

@@ -0,0 +1,640 @@
+module photon_physics
+
+  use algorithm,       only: binary_search
+  use constants
+  use particle_header, only: Particle
+  use photon_header,   only: PhotonInteraction, BremsstrahlungData, &
+                             compton_profile_pz, ttb_e_grid, ttb
+  use random_lcg,      only: prn
+  use settings
+
+contains
+
+!===============================================================================
+! KLEIN_NISHINA
+!===============================================================================
+
+  subroutine klein_nishina(alpha, alpha_out, mu)
+    real(8), intent(in)  :: alpha
+    real(8), intent(out) :: alpha_out
+    real(8), intent(out) :: mu
+
+    real(8) :: beta ! 1 + 2a
+    real(8) :: t    ! (1 + 2a)/(9 + 2a)
+    real(8) :: r, s, x
+    real(8) :: gamma
+
+    beta = ONE + TWO*alpha
+    if (alpha < THREE) then
+      ! Kahn's rejection method
+      t = beta/(beta + 8.0_8)
+      do
+        if (prn() < t) then
+          ! Left branch of flow chart
+          r = TWO*prn()
+          x = ONE + alpha*r
+          if (prn() < FOUR/x*(ONE - ONE/x)) then
+            mu = 1 - r
+            exit
+          end if
+        else
+          ! Right branch of flow chart
+          x = beta/(ONE + TWO*alpha*prn())
+          mu = ONE + (ONE - x)/alpha
+          if (prn() < HALF*(mu**2 + ONE/x)) exit
+        end if
+      end do
+      alpha_out = alpha/x
+
+    else
+      ! Koblinger's direct method
+      gamma = ONE - beta**(-2)
+      s = prn()*(FOUR/alpha + HALF*gamma + &
+           (ONE - (ONE + beta)/alpha**2)*log(beta))
+      if (s <= 2./alpha) then
+        ! For first term, x = 1 + 2ar
+        ! Therefore, a' = a/(1 + 2ar)
+        alpha_out = alpha/(ONE + TWO*alpha*prn())
+      elseif (s <= FOUR/alpha) then
+        ! For third term, x = beta/(1 + 2ar)
+        ! Therefore, a' = a(1 + 2ar)/beta
+        alpha_out = alpha*(ONE + TWO*alpha*prn())/beta
+      elseif (s <= FOUR/alpha + HALF*gamma) then
+        ! For fourth term, x = 1/sqrt(1 - gamma*r)
+        ! Therefore, a' = a*sqrt(1 - gamma*r)
+        alpha_out = alpha*sqrt(ONE - gamma*prn())
+      else
+        ! For third term, x = beta^r
+        ! Therefore, a' = a/beta^r
+        alpha_out = alpha/beta**prn()
+      end if
+
+      ! Calculate cosine of scattering angle based on basic relation
+      mu = ONE + ONE/alpha - ONE/alpha_out
+    end if
+
+  end subroutine klein_nishina
+
+!===============================================================================
+! COMPTON_SCATTER
+!===============================================================================
+
+  subroutine compton_scatter(el, alpha, alpha_out, mu, i_shell, use_doppler)
+    type(PhotonInteraction), intent(in) :: el
+    real(8), intent(in)  :: alpha
+    real(8), intent(out) :: alpha_out
+    real(8), intent(out) :: mu
+    integer, intent(out) :: i_shell
+    logical, intent(in), optional :: use_doppler
+
+    real(8) :: x
+    real(8) :: form_factor_xmax
+    real(8) :: form_factor_x
+    real(8) :: e_out
+    logical :: use_doppler_
+
+    if (present(use_doppler)) then
+      use_doppler_ = use_doppler
+    else
+      use_doppler_ = .false.
+    end if
+
+    form_factor_xmax = ZERO
+    do
+      ! Sample Klein-Nishina distribution for trial energy and angle
+      call klein_nishina(alpha, alpha_out, mu)
+
+      ! Note that the parameter used here does not correspond exactly to the
+      ! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
+      ! parameter as defined by Hubbell, where the actual data comes from
+      x = MASS_ELECTRON_EV/PLANCK_C*alpha*sqrt(HALF*(ONE - mu))
+
+      ! Calculate S(x, Z) and S(x_max, Z)
+      form_factor_x = el % incoherent_form_factor % evaluate(x)
+      if (form_factor_xmax == ZERO) then
+        form_factor_xmax = el % incoherent_form_factor % evaluate(&
+             MASS_ELECTRON_EV/PLANCK_C*alpha)
+      end if
+
+      ! Perform rejection on form factor
+      if (prn() < form_factor_x / form_factor_xmax) then
+        if (use_doppler_) then
+          call compton_doppler(el, alpha, mu, e_out, i_shell)
+          alpha_out = e_out/MASS_ELECTRON_EV
+        else
+          i_shell = 0
+        end if
+        exit
+      end if
+    end do
+
+  end subroutine compton_scatter
+
+!===============================================================================
+! COMPTON_DOPPLER
+!===============================================================================
+
+  subroutine compton_doppler(el, alpha, mu, e_out, i_shell)
+    type(PhotonInteraction), intent(in) :: el
+    real(8), intent(in)  :: alpha
+    real(8), intent(in)  :: mu
+    real(8), intent(out) :: e_out
+    integer, intent(out) :: i_shell
+
+    integer :: i
+    integer :: n
+    real(8) :: rn, m
+    real(8) :: c, c_l, c_max
+    real(8) :: pz_l, pz_r, pz, pz_max
+    real(8) :: p_l, p_r
+    real(8) :: e, e_b
+    real(8) :: e_out1, e_out2
+    real(8) :: a, b, quad
+    real(8) :: f
+    real(8) :: momentum_sq
+
+    n = size(compton_profile_pz)
+
+    do
+      ! Sample electron shell
+      rn = prn()
+      c = ZERO
+      do i_shell = 1, size(el % electron_pdf)
+        c = c + el % electron_pdf(i_shell)
+        if (rn < c) exit
+      end do
+
+      ! Determine binding energy of shell
+      e_b = el % binding_energy(i_shell)
+
+      ! Determine p_z,max
+      e = alpha*MASS_ELECTRON_EV
+      if (e < e_b) then
+        e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
+        exit
+      end if
+
+      pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / &
+           sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2)
+      if (pz_max < ZERO) then
+        e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
+        exit
+      end if
+
+      ! Determine profile cdf value corresponding to p_z,max
+      if (pz_max > compton_profile_pz(n)) then
+        c_max = el % profile_cdf(n, i_shell)
+      else
+        i = binary_search(compton_profile_pz, n, pz_max)
+        pz_l = compton_profile_pz(i)
+        pz_r = compton_profile_pz(i + 1)
+        p_l = el % profile_pdf(i, i_shell)
+        p_r = el % profile_pdf(i + 1, i_shell)
+        c_l = el % profile_cdf(i, i_shell)
+        if (pz_l == pz_r) then
+          c_max = c_l
+        elseif (p_l == p_r) then
+          c_max = c_l + (pz_max - pz_l)*p_l
+        else
+          m = (p_l - p_r)/(pz_l - pz_r)
+          c_max = c_l + ((m*(pz_max - pz_l) + p_l)**2 - p_l**2)/(TWO*m)
+        end if
+      end if
+
+      ! Sample value on bounded cdf
+      c = prn()*c_max
+
+      ! Determine pz corresponding to sampled cdf value
+      i = binary_search(el % profile_cdf(:, i_shell), n, c)
+      pz_l = compton_profile_pz(i)
+      pz_r = compton_profile_pz(i + 1)
+      p_l = el % profile_pdf(i, i_shell)
+      p_r = el % profile_pdf(i + 1, i_shell)
+      c_l = el % profile_cdf(i, i_shell)
+      if (pz_l == pz_r) then
+        pz = pz_l
+      elseif (p_l == p_r) then
+        pz = pz_l + (c - c_l)/p_l
+      else
+        m = (p_l - p_r)/(pz_l - pz_r)
+        pz = pz_l + (sqrt(p_l**2 + TWO*m*(c - c_l)) - p_l)/m
+      end if
+
+      ! Determine outgoing photon energy corresponding to electron momentum
+      ! (solve Eq. 39 in LA-UR-04-0487 for E')
+      momentum_sq = (pz/FINE_STRUCTURE)**2
+      f = ONE + alpha*(ONE - mu)
+      a = momentum_sq - f*f
+      b = TWO*e*(f - momentum_sq*mu)
+      c = e**2*(momentum_sq - ONE)
+
+      quad = b**2 - FOUR*a*c
+      if (quad < 0) then
+        e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
+        exit
+      end if
+      quad = sqrt(quad)
+      e_out1 = -(b + quad)/(TWO*a)
+      e_out2 = -(b - quad)/(TWO*a)
+
+      ! Determine solution to quadratic equation that is positive
+      if (e_out1 > ZERO) then
+        if (e_out2 > ZERO) then
+          ! If both are positive, pick one at random
+          if (prn() < HALF) then
+            e_out = e_out1
+          else
+            e_out = e_out2
+          end if
+        else
+          e_out = e_out1
+        end if
+      else
+        if (e_out2 > ZERO) then
+          e_out = e_out2
+        else
+          ! No positive solution -- resample
+          cycle
+        end if
+      end if
+      if (e_out < e - e_b) exit
+    end do
+
+  end subroutine compton_doppler
+
+!===============================================================================
+! RAYLEIGH_SCATTER
+!===============================================================================
+
+  subroutine rayleigh_scatter(el, alpha, mu)
+    type(PhotonInteraction), intent(in) :: el
+    real(8), intent(in)  :: alpha
+    real(8), intent(out) :: mu
+
+    integer :: i
+    real(8) :: F
+    real(8) :: F_max
+    real(8) :: x2
+    real(8) :: x2_max
+    real(8) :: r
+
+    do
+      ! Determine maximum value of x^2
+      x2_max = (MASS_ELECTRON_EV/PLANCK_C*alpha)**2
+
+      ! Determine F(x^2_max, Z)
+      F_max = el % coherent_int_form_factor % evaluate(x2_max)
+
+      ! Sample cumulative distribution
+      F = prn()*F_max
+
+      ! Determine x^2 corresponding to F
+      i = binary_search(el%coherent_int_form_factor%y, &
+           size(el%coherent_int_form_factor%y), F)
+      r = (F - el%coherent_int_form_factor%y(i)) / &
+           (el%coherent_int_form_factor%y(i+1) - el%coherent_int_form_factor%y(i))
+      x2 = el%coherent_int_form_factor%x(i) + r*(el%coherent_int_form_factor%x(i+1) - &
+           el%coherent_int_form_factor%x(i))
+
+      ! Calculate mu
+      mu = ONE - TWO*x2/x2_max
+
+      if (prn() < HALF*(ONE + mu**2)) exit
+    end do
+
+  end subroutine rayleigh_scatter
+
+!===============================================================================
+! ATOMIC_RELAXATION
+!===============================================================================
+
+  recursive subroutine atomic_relaxation(p, elm, i_shell)
+    type(Particle), intent(inout) :: p
+    type(PhotonInteraction), intent(in) :: elm
+    integer, intent(in) :: i_shell
+
+    integer :: i_hole
+    integer :: i_transition
+    integer :: primary
+    integer :: secondary
+    real(8) :: c
+    real(8) :: rn
+    real(8) :: E
+    real(8) :: mu
+    real(8) :: phi
+    real(8) :: uvw(3)
+
+    ! If no transitions, assume fluorescent photon from captured free electron
+    if (elm % shells(i_shell) % n_transitions == 0) then
+      mu = TWO*prn() - ONE
+      phi = TWO*PI*prn()
+      uvw(1) = mu
+      uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
+      uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
+      E = elm % shells(i_shell) % binding_energy
+      call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
+      return
+    end if
+
+    ! Sample transition
+    rn = prn()
+    c = ZERO
+    do i_transition = 1, elm % shells(i_shell) % n_transitions
+      c = c + elm % shells(i_shell) % &
+           transition_probability(i_transition)
+      if (rn < c) exit
+    end do
+
+    ! Get primary and secondary subshell designators
+    primary = elm % shells(i_shell) % transition_subshells(1, i_transition)
+    secondary = elm % shells(i_shell) % transition_subshells(2, i_transition)
+
+    ! Sample angle isotropically
+    mu = TWO*prn() - ONE
+    phi = TWO*PI*prn()
+    uvw(1) = mu
+    uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
+    uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
+
+    ! Get the transition energy
+    E = elm % shells(i_shell) % transition_energy(i_transition)
+
+    if (secondary /= 0) then
+      ! Non-radiative transition -- Auger/Coster-Kronig effect
+
+      ! Create auger electron
+      call p % create_secondary(uvw, E, ELECTRON, run_ce=.true.)
+
+      ! Fill hole left by emitted auger electron
+      i_hole = elm % shell_dict % get(secondary)
+      call atomic_relaxation(p, elm, i_hole)
+    else
+      ! Radiative transition -- get X-ray energy
+
+      ! Create fluorescent photon
+      call p % create_secondary(uvw, E, PHOTON, run_ce=.true.)
+
+    end if
+
+    ! Fill hole created by electron transitioning to the photoelectron hole
+    i_hole = elm % shell_dict % get(primary)
+    call atomic_relaxation(p, elm, i_hole)
+
+  end subroutine atomic_relaxation
+
+!===============================================================================
+! PAIR_PRODUCTION samples the kinetic energy and direction of the electron and
+! positron created when a photon is absorbed near an atomic nucleus. The
+! simulation procedure follows the semiempirical model outlined in F. Salvat, J.
+! M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for Monte
+! Carlo Simulation of Electron and Photon Transport," OECD-NEA,
+! Issy-les-Moulineaux, France (2011).
+!===============================================================================
+
+  subroutine pair_production(elm, alpha, E_electron, E_positron, mu_electron, &
+       mu_positron)
+    type(PhotonInteraction), intent(in) :: elm
+    real(8), intent(in)  :: alpha
+    real(8), intent(out) :: E_electron
+    real(8), intent(out) :: E_positron
+    real(8), intent(out) :: mu_electron
+    real(8), intent(out) :: mu_positron
+
+    integer :: i
+    real(8) :: f
+    real(8) :: c
+    real(8) :: a
+    real(8) :: b
+    real(8) :: q
+    real(8) :: rn
+    real(8) :: beta
+    real(8) :: e, e_min, e_max
+    real(8) :: t1, t2, t3, t4
+    real(8) :: u1, u2
+    real(8) :: phi1, phi2
+    real(8) :: phi1_max, phi2_max
+    real(8), parameter :: r(99) = (/ &
+         122.81_8, 73.167_8, 69.228_8, 67.301_8, 64.696_8, 61.228_8, &
+         57.524_8, 54.033_8, 50.787_8, 47.851_8, 46.373_8, 45.401_8, &
+         44.503_8, 43.815_8, 43.074_8, 42.321_8, 41.586_8, 40.953_8, &
+         40.524_8, 40.256_8, 39.756_8, 39.144_8, 38.462_8, 37.778_8, &
+         37.174_8, 36.663_8, 35.986_8, 35.317_8, 34.688_8, 34.197_8, &
+         33.786_8, 33.422_8, 33.068_8, 32.740_8, 32.438_8, 32.143_8, &
+         31.884_8, 31.622_8, 31.438_8, 31.142_8, 30.950_8, 30.758_8, &
+         30.561_8, 30.285_8, 30.097_8, 29.832_8, 29.581_8, 29.411_8, &
+         29.247_8, 29.085_8, 28.930_8, 28.721_8, 28.580_8, 28.442_8, &
+         28.312_8, 28.139_8, 27.973_8, 27.819_8, 27.675_8, 27.496_8, &
+         27.285_8, 27.093_8, 26.911_8, 26.705_8, 26.516_8, 26.304_8, &
+         26.108_8, 25.929_8, 25.730_8, 25.577_8, 25.403_8, 25.245_8, &
+         25.100_8, 24.941_8, 24.790_8, 24.655_8, 24.506_8, 24.391_8, &
+         24.262_8, 24.145_8, 24.039_8, 23.922_8, 23.813_8, 23.712_8, &
+         23.621_8, 23.523_8, 23.430_8, 23.331_8, 23.238_8, 23.139_8, &
+         23.048_8, 22.967_8, 22.833_8, 22.694_8, 22.624_8, 22.545_8, &
+         22.446_8, 22.358_8, 22.264_8 /)
+
+    ! The reduced screening radius r is the ratio of the screening radius to
+    ! the Compton wavelength of the electron, where the screening radius is
+    ! obtained under the assumption that the Coulomb field of the nucleus is
+    ! exponentially screened by atomic electrons. This allows us to use a
+    ! simplified atomic form factor and analytical approximations of the
+    ! screening functions in the pair production DCS instead of computing the
+    ! screening functions numerically. The reduced screening radii above for
+    ! Z = 1-99 come from F. Salvat, J. M. Fernández-Varea, and J. Sempau,
+    ! "PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and
+    ! Photon Transport," OECD-NEA, Issy-les-Moulineaux, France (2011).
+
+    ! Compute the minimum and maximum values of the electron reduced energy,
+    ! i.e. the fraction of the photon energy that is given to the electron
+    e_min = ONE/alpha
+    e_max = ONE - ONE/alpha
+
+    ! Compute the high-energy Coulomb correction
+    a = elm % Z / FINE_STRUCTURE
+    c = a**2*(ONE/(ONE + a**2) + 0.202059_8 - 0.03693_8*a**2 + 0.00835_8*a**4 &
+         - 0.00201_8*a**6 + 0.00049_8*a**8 - 0.00012_8*a**10 + 0.00003_8*a**12)
+
+    ! The analytical approximation of the DCS underestimates the cross section
+    ! at low energies. The correction factor f compensates for this.
+    q = sqrt(TWO/alpha)
+    f = q*(-0.1774_8 - 12.10_8*a + 11.18_8*a**2) &
+         + q**2*(8.523_8 + 73.26_8*a - 44.41_8*a**2) &
+         + q**3*(-13.52_8 - 121.1_8*a + 96.41_8*a**2) &
+         + q**4*(8.946_8 + 62.05_8*a - 63.41_8*a**2)
+
+    ! Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced
+    ! energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and
+    ! phi_2 are non-negative and maximum at e = 1/2.
+    b = TWO*r(elm % Z)/alpha
+    t1 = TWO*log(ONE + b**2)
+    t2 = b*atan(ONE/b)
+    t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2))
+    t4 = FOUR*log(r(elm % Z)) - FOUR*c + f
+    phi1_max = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4
+    phi2_max = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4
+
+    ! To aid sampling, the unnormalized PDF can be expressed as
+    ! p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are
+    ! normalized PDFs on the interval (e_min, e_max) from which values of e can
+    ! be sampled using the inverse transform method, and
+    ! U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid
+    ! rejection functions. The reduced energy can now be sampled using a
+    ! combination of the composition and rejection methods.
+    u1 = TWO/THREE*(HALF - ONE/alpha)**2*phi1_max
+    u2 = phi2_max
+    do
+      rn = prn()
+
+      ! Sample the index i in (1, 2) using the point probabilities
+      ! p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2)
+      if (prn() < u1/(u1 + u2)) then
+        i = 1
+
+        ! Sample e from pi_1 using the inverse transform method
+        if (rn >= HALF) then
+          e = HALF + (HALF - ONE/alpha)*(TWO*rn - ONE)**(ONE/THREE)
+        else
+          e = HALF - (HALF - ONE/alpha)*(ONE - TWO*rn)**(ONE/THREE)
+        end if
+      else
+        i = 2
+
+        ! Sample e from pi_2 using the inverse transform method
+        e = ONE/alpha + (HALF - ONE/alpha)*TWO*rn
+      end if
+
+      ! Calculate phi_i(e) and deliver e if rn <= U_i(e)
+      b = r(elm % Z)/(TWO*alpha*e*(ONE - e))
+      t1 = TWO*log(ONE + b**2)
+      t2 = b*atan(ONE/b)
+      t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2))
+      if (i == 1) then
+        phi1 = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4
+        if (prn() <= phi1/phi1_max) exit
+      else
+        phi2 = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4
+        if (prn() <= phi2/phi2_max) exit
+      end if
+    end do
+
+    ! Compute the kinetic energy of the electron and the positron
+    E_electron = (alpha*e - ONE)*MASS_ELECTRON_EV
+    E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON_EV
+
+    ! Sample the scattering angle of the electron. The cosine of the polar
+    ! angle of the direction relative to the incident photon is sampled from
+    ! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
+    beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON_EV)) &
+         / (E_electron + MASS_ELECTRON_EV)
+    rn = TWO*prn() - ONE
+    mu_electron = (rn + beta)/(rn*beta + ONE)
+
+    ! Sample the scattering angle of the positron
+    beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON_EV)) &
+         / (E_positron + MASS_ELECTRON_EV)
+    rn = TWO*prn() - ONE
+    mu_positron = (rn + beta)/(rn*beta + ONE)
+
+  end subroutine pair_production
+
+!===============================================================================
+! THICK_TARGET_BREMSSTRAHLUNG
+!===============================================================================
+
+  subroutine thick_target_bremsstrahlung(p, E_lost)
+    type(Particle), intent(inout) :: p
+    real(8),        intent(inout) :: E_lost
+
+    integer :: i, j
+    integer :: i_e, i_w
+    integer :: n
+    integer :: n_e
+    real(8) :: a
+    real(8) :: f
+    real(8) :: e, e_l, e_r
+    real(8) :: y, y_l, y_r
+    real(8) :: w, w_l, w_r
+    real(8) :: p_l, p_r
+    real(8) :: c, c_l, c_max
+    type(BremsstrahlungData), pointer :: mat
+
+    if (p % E < energy_cutoff(PHOTON)) return
+
+    ! Get bremsstrahlung data for this material and particle type
+    if (p % type == POSITRON) then
+      mat => ttb(p % material) % positron
+    else
+      mat => ttb(p % material) % electron
+    end if
+
+    e = log(p % E)
+    n_e = size(ttb_e_grid)
+
+    ! Find the lower bounding index of the incident electron energy
+    j = binary_search(ttb_e_grid, n_e, e)
+    if (j == n_e) j = j - 1
+
+    ! Get the interpolation bounds
+    e_l = ttb_e_grid(j)
+    e_r = ttb_e_grid(j+1)
+    y_l = mat % yield(j)
+    y_r = mat % yield(j+1)
+
+    ! Calculate the interpolation weight w_j+1 of the bremsstrahlung energy PDF
+    ! interpolated in log energy, which can be interpreted as the probability
+    ! of index j+1
+    f = (e - e_l)/(e_r - e_l)
+
+    ! Get the photon number yield for the given energy using linear
+    ! interpolation on a log-log scale
+    y = exp(y_l + (y_r - y_l)*f)
+
+    ! Sample number of secondary bremsstrahlung photons
+    n = int(y + prn())
+
+    E_lost = ZERO
+    if (n == 0) return
+
+    ! Sample index of the tabulated PDF in the energy grid, j or j+1
+    if (prn() <= f .or. j == 1) then
+      i_e = j + 1
+
+      ! Interpolate the maximum value of the CDF at the incoming particle
+      ! energy on a log-log scale
+      p_l = mat % pdf(i_e-1, i_e)
+      p_r = mat % pdf(i_e, i_e)
+      c_l = mat % cdf(i_e-1, i_e)
+      a = log(p_r/p_l)/(e_r - e_l) + ONE
+      c_max = c_l + exp(e_l)*p_l/a*(exp(a*(e - e_l)) - ONE)
+    else
+      i_e = j
+
+      ! Maximum value of the CDF
+      c_max = mat % cdf(i_e, i_e)
+    end if
+
+    ! Sample the energies of the emitted photons
+    do i = 1, n
+      ! Generate a random number r and determine the index i for which
+      ! cdf(i) <= r*cdf,max <= cdf(i+1)
+      c = prn()*c_max
+      i_w = binary_search(mat % cdf(:i_e,i_e), i_e, c)
+
+      ! Sample the photon energy
+      w_l = ttb_e_grid(i_w)
+      w_r = ttb_e_grid(i_w+1)
+      p_l = mat % pdf(i_w, i_e)
+      p_r = mat % pdf(i_w+1, i_e)
+      c_l = mat % cdf(i_w, i_e)
+      a = log(p_r/p_l)/(w_r - w_l) + ONE
+      w = exp(w_l)*(a*(c - c_l)/(exp(w_l)*p_l) + ONE)**(ONE/a)
+
+      if (w > energy_cutoff(PHOTON)) then
+        ! Create secondary photon
+        call p % create_secondary(p % coord(1) % uvw, w, PHOTON, run_ce=.true.)
+        E_lost = E_lost + w
+      end if
+    end do
+
+  end subroutine thick_target_bremsstrahlung
+
+end module photon_physics

src/physics.F90

@@ -10,6 +10,10 @@ module physics
   use message_passing
   use nuclide_header
   use particle_header,        only: Particle
+  use photon_header
+  use photon_physics,         only: rayleigh_scatter, compton_scatter, &
+                                    atomic_relaxation, pair_production, &
+                                    thick_target_bremsstrahlung
   use physics_common
   use random_lcg,             only: prn, advance_prn_seed, prn_set_stream
   use reaction_header,        only: Reaction
@@ -36,40 +40,48 @@ contains
     ! Add to collision counter for particle
     p % n_collision = p % n_collision + 1
 
-    ! Sample nuclide/reaction for the material the particle is in
-    call sample_reaction(p)
+    ! Sample reaction for the material the particle is in
+    if (p % type == NEUTRON) then
+      call sample_neutron_reaction(p)
+    else if (p % type == PHOTON) then
+      call sample_photon_reaction(p)
+    else if (p % type == ELECTRON) then
+      call sample_electron_reaction(p)
+    else if (p % type == POSITRON) then
+      call sample_positron_reaction(p)
+    end if
+
+    ! Kill particle if energy falls below cutoff
+    if (p % E < energy_cutoff(p % type)) then
+      p % alive = .false.
+      p % wgt = ZERO
+      p % last_wgt = ZERO
+    end if
 
     ! Display information about collision
     if (verbosity >= 10 .or. trace) then
-      call write_message("    " // trim(reaction_name(p % event_MT)) &
-           &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) &
-           &// ". Energy = " // trim(to_str(p % E)) // " eV.")
+      if (p % type == NEUTRON) then
+        call write_message("    " // trim(reaction_name(p % event_MT)) &
+             &// " with " // trim(adjustl(nuclides(p % event_nuclide) % name)) &
+             &// ". Energy = " // trim(to_str(p % E)) // " eV.")
+      else
+        call write_message("    " // trim(reaction_name(p % event_MT)) &
+             &// " with " // trim(adjustl(elements(p % event_nuclide) % name)) &
+             &// ". Energy = " // trim(to_str(p % E)) // " eV.")
+      end if
     end if
 
-    ! check for very low energy
-    if (p % E < 1.0e-100_8) then
-      p % alive = .false.
-      if (master) call warning("Killing neutron with extremely low energy")
-    end if
-
-    ! Advance URR seed stream 'N' times after energy changes
-    if (p % E /= p % last_E) then
-      call prn_set_stream(STREAM_URR_PTABLE)
-      call advance_prn_seed(size(nuclides, kind=8))
-      call prn_set_stream(STREAM_TRACKING)
-    endif
-
   end subroutine collision
 
 !===============================================================================
-! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
-! each nuclide within a material and then samples a reaction for that nuclide
-! and calls the appropriate routine to process the physics. Note that there is
-! special logic when suvival biasing is turned on since fission and
-! disappearance are treated implicitly.
+! SAMPLE_NEUTRON_REACTION samples a nuclide based on the macroscopic cross
+! sections for each nuclide within a material and then samples a reaction for
+! that nuclide and calls the appropriate routine to process the physics. Note
+! that there is special logic when suvival biasing is turned on since fission
+! and disappearance are treated implicitly.
 !===============================================================================
 
-  subroutine sample_reaction(p)
+  subroutine sample_neutron_reaction(p)
 
     type(Particle), intent(inout) :: p
 
@@ -102,6 +114,13 @@ contains
       end if
     end if
 
+    ! Create secondary photons
+    if (photon_transport) then
+      call prn_set_stream(STREAM_PHOTON)
+      call sample_secondary_photons(p, i_nuclide)
+      call prn_set_stream(STREAM_TRACKING)
+    end if
+
     ! If survival biasing is being used, the following subroutine adjusts the
     ! weight of the particle. Otherwise, it checks to see if absorption occurs
 
@@ -116,21 +135,258 @@ contains
     ! exiting neutron
     call scatter(p, i_nuclide, i_nuc_mat)
 
-    ! Play russian roulette if survival biasing is turned on
+    ! Advance URR seed stream 'N' times after energy changes
+    if (p % E /= p % last_E) then
+      call prn_set_stream(STREAM_URR_PTABLE)
+      call advance_prn_seed(size(nuclides, kind=8))
+      call prn_set_stream(STREAM_TRACKING)
+    end if
 
+    ! Play russian roulette if survival biasing is turned on
     if (survival_biasing) then
       call russian_roulette(p)
       if (.not. p % alive) return
     end if
 
-    ! Kill neutron under certain energy
-    if (p % E < energy_cutoff) then
+  end subroutine sample_neutron_reaction
+
+!===============================================================================
+! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross
+! sections for each nuclide within a material and then samples a reaction for
+! that element and calls the appropriate routine to process the physics.
+!===============================================================================
+
+  subroutine sample_photon_reaction(p)
+    type(Particle), intent(inout) :: p
+
+    integer :: i_shell      ! index in subshells
+    integer :: i_grid       ! index on energy grid
+    integer :: i_element    ! index in nuclides array
+    integer :: i_start      ! threshold index
+    real(8) :: prob         ! cumulative probability
+    real(8) :: cutoff       ! sampled total cross section
+    real(8) :: f            ! interpolation factor
+    real(8) :: xs           ! photoionization cross section
+    real(8) :: r            ! random number
+    real(8) :: prob_after
+    real(8) :: alpha        ! photon energy divided by electron rest mass
+    real(8) :: alpha_out    ! outgoing photon energy over electron rest mass
+    real(8) :: mu           ! scattering cosine
+    real(8) :: mu_electron  ! electron scattering cosine
+    real(8) :: mu_positron  ! positron scattering cosine
+    real(8) :: phi          ! azimuthal angle
+    real(8) :: uvw(3)       ! new direction
+    real(8) :: rel_vel      ! relative velocity of electron
+    real(8) :: e_b          ! binding energy of electron
+    real(8) :: E_electron   ! electron energy
+    real(8) :: E_positron   ! positron energy
+
+    ! Kill photon if below energy cutoff -- an extra check is made here because
+    ! photons with energy below the cutoff may have been produced by neutrons
+    ! reactions or atomic relaxation
+    if (p % E < energy_cutoff(PHOTON)) then
+      p % E = ZERO
       p % alive = .false.
-      p % wgt = ZERO
-      p % last_wgt = ZERO
+      return
     end if
 
-  end subroutine sample_reaction
+    ! Sample element within material
+    i_element = sample_element(p)
+    p % event_nuclide = i_element
+
+    ! Calculate photon energy over electron rest mass equivalent
+    alpha = p % E/MASS_ELECTRON_EV
+
+    ! For tallying purposes, this routine might be called directly. In that
+    ! case, we need to sample a reaction via the cutoff variable
+    prob = ZERO
+    cutoff = prn() * micro_photon_xs(i_element) % total
+
+    associate (elm => elements(i_element))
+      ! Coherent (Rayleigh) scattering
+      prob = prob + micro_photon_xs(i_element) % coherent
+      if (prob > cutoff) then
+        call rayleigh_scatter(elm, alpha, mu)
+        p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
+        p % event_MT = COHERENT
+        return
+      end if
+
+      ! Incoherent (Compton) scattering
+      prob = prob + micro_photon_xs(i_element) % incoherent
+      if (prob > cutoff) then
+        call compton_scatter(elm, alpha, alpha_out, mu, i_shell, .true.)
+
+        ! Determine binding energy of shell. The binding energy is zero if
+        ! doppler broadening is not used.
+        if (i_shell == 0) then
+          e_b = ZERO
+        else
+          e_b = elm % binding_energy(i_shell)
+        end if
+
+        ! Create Compton electron
+        E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b
+        mu_electron = (alpha - alpha_out*mu) &
+             / sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu)
+        phi = TWO*PI*prn()
+        uvw = rotate_angle(p % coord(1) % uvw, mu_electron, phi)
+        call p % create_secondary(uvw, E_electron, ELECTRON, .true.)
+
+        ! TODO: Compton subshell data does not match atomic relaxation data
+        ! Allow electrons to fill orbital and produce auger electrons
+        ! and fluorescent photons
+        if (i_shell > 0) then
+          call atomic_relaxation(p, elm, i_shell)
+        end if
+
+        phi = phi + PI
+        p % E = alpha_out*MASS_ELECTRON_EV
+        p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu, phi)
+        p % event_MT = INCOHERENT
+        return
+      end if
+
+      ! Photoelectric effect
+      prob_after = prob + micro_photon_xs(i_element) % photoelectric
+      if (prob_after > cutoff) then
+        do i_shell = 1, size(elm % shells)
+          ! Get grid index and interpolation factor
+          i_grid = micro_photon_xs(i_element) % index_grid
+          f      = micro_photon_xs(i_element) % interp_factor
+
+          ! Check threshold of reaction
+          i_start = elm % shells(i_shell) % threshold
+          if (i_grid <= i_start) cycle
+
+          ! Evaluation subshell photoionization cross section
+          xs = exp(elm % shells(i_shell) % cross_section(i_grid - i_start) + &
+               f*(elm % shells(i_shell) % cross_section(i_grid + 1 - i_start) - &
+               elm % shells(i_shell) % cross_section(i_grid - i_start)))
+
+          prob = prob + xs
+          if (prob > cutoff) then
+            E_electron = p % E - elm % shells(i_shell) % binding_energy
+
+            ! Sample mu using non-relativistic Sauter distribution.
+            ! See Eqns 3.19 and 3.20 in "Implementing a photon physics
+            ! model in Serpent 2" by Toni Kaltiaisenaho
+            SAMPLE_MU: do
+              r = prn()
+              if (FOUR * (ONE - r) * r >= prn()) then
+                rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON_EV))&
+                     / (E_electron + MASS_ELECTRON_EV)
+                mu = (TWO * r + rel_vel - ONE) / &
+                     (TWO * rel_vel * r - rel_vel + ONE)
+                exit SAMPLE_MU
+              end if
+            end do SAMPLE_MU
+
+            phi = TWO*PI*prn()
+            uvw(1) = mu
+            uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
+            uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
+
+            ! Create secondary electron
+            call p % create_secondary(uvw, E_electron, ELECTRON, run_CE=.true.)
+
+            ! Allow electrons to fill orbital and produce auger electrons
+            ! and fluorescent photons
+            call atomic_relaxation(p, elm, i_shell)
+            p % event_MT = 533 + elm % shells(i_shell) % index_subshell
+            p % alive = .false.
+            p % E = ZERO
+
+            return
+          end if
+        end do
+      end if
+      prob = prob_after
+
+      ! Pair production
+      prob = prob + micro_photon_xs(i_element) % pair_production
+      if (prob > cutoff) then
+        call pair_production(elm, alpha, E_electron, E_positron, mu_electron, &
+             mu_positron)
+
+        ! Create secondary electron
+        uvw = rotate_angle(p % coord(1) % uvw, mu_electron)
+        call p % create_secondary(uvw, E_electron, ELECTRON, .true.)
+
+        ! Create secondary positron
+        uvw = rotate_angle(p % coord(1) % uvw, mu_positron)
+        call p % create_secondary(uvw, E_positron, POSITRON, .true.)
+
+        p % event_MT = PAIR_PROD
+        p % alive = .false.
+        p % E = ZERO
+      end if
+
+    end associate
+
+  end subroutine sample_photon_reaction
+
+!===============================================================================
+! SAMPLE_ELECTRON_REACTION terminates the particle and either deposits all
+! energy locally (electron_treatment = ELECTRON_LED) or creates secondary
+! bremsstrahlung photons from electron deflections with charged particles
+! (electron_treatment = ELECTRON_TTB).
+!===============================================================================
+
+  subroutine sample_electron_reaction(p)
+    type(Particle), intent(inout) :: p
+
+    real(8) :: E_lost  ! energy lost to bremsstrahlung photons
+
+    ! TODO: create reaction types
+
+    if (electron_treatment == ELECTRON_TTB) then
+      call thick_target_bremsstrahlung(p, E_lost)
+    end if
+
+    p % E = ZERO
+    p % alive = .false.
+
+  end subroutine sample_electron_reaction
+
+!===============================================================================
+! SAMPLE_POSITRON_REACTION terminates the particle and either deposits all
+! energy locally (electron_treatment = ELECTRON_LED) or creates secondary
+! bremsstrahlung photons from electron deflections with charged particles
+! (electron_treatment = ELECTRON_TTB). Two annihilation photons of energy
+! MASS_ELECTRON_EV (0.511 MeV) are created and travel in opposite directions.
+!===============================================================================
+
+  subroutine sample_positron_reaction(p)
+    type(Particle), intent(inout) :: p
+
+    real(8) :: mu     ! scattering cosine
+    real(8) :: phi    ! azimuthal angle
+    real(8) :: uvw(3) ! new direction
+
+    real(8) :: E_lost  ! energy lost to bremsstrahlung photons
+
+    ! TODO: create reaction types
+
+    if (electron_treatment == ELECTRON_TTB) then
+      call thick_target_bremsstrahlung(p, E_lost)
+    end if
+
+    ! Sample angle isotropically
+    mu = TWO*prn() - ONE
+    phi = TWO*PI*prn()
+    uvw(1) = mu
+    uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
+    uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
+
+    ! Create annihilation photon pair traveling in opposite directions
+    call p % create_secondary( uvw, MASS_ELECTRON_EV, PHOTON, .true.)
+    call p % create_secondary(-uvw, MASS_ELECTRON_EV, PHOTON, .true.)
+
+    p % E = ZERO
+    p % alive = .false.
+
+  end subroutine sample_positron_reaction
 
 !===============================================================================
 ! SAMPLE_NUCLIDE
@@ -194,6 +450,49 @@ contains
 
   end subroutine sample_nuclide
 
+!===============================================================================
+! SAMPLE_ELEMENT
+!===============================================================================
+
+  function sample_element(p) result(i_element)
+    type(Particle), intent(in) :: p
+    integer                    :: i_element
+
+    integer :: i
+    real(8) :: prob
+    real(8) :: cutoff
+    real(8) :: atom_density ! atom density of nuclide in atom/b-cm
+    real(8) :: sigma        ! microscopic total xs for nuclide
+
+    associate (mat => materials(p % material))
+      ! Sample cumulative distribution function
+      cutoff = prn() * material_xs % total
+
+      i = 0
+      prob = ZERO
+      do while (prob < cutoff)
+        i = i + 1
+
+        ! Check to make sure that a nuclide was sampled
+        if (i > mat % n_nuclides) then
+          call p % write_restart()
+          call fatal_error("Did not sample any element during collision.")
+        end if
+
+        ! Find atom density
+        i_element    = mat % element(i)
+        atom_density = mat % atom_density(i)
+
+        ! Determine microscopic cross section
+        sigma = atom_density * micro_photon_xs(i_element) % total
+
+        ! Increment probability to compare to cutoff
+        prob = prob + sigma
+      end do
+    end associate
+
+  end function sample_element
+
 !===============================================================================
 ! SAMPLE_FISSION
 !===============================================================================
@@ -260,6 +559,65 @@ contains
 
   end subroutine sample_fission
 
+!===============================================================================
+! SAMPLE_PHOTON_PRODUCT
+!===============================================================================
+
+  subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product)
+    integer, intent(in)  :: i_nuclide  ! index in nuclides array
+    real(8), intent(in)  :: E          ! energy of neutron
+    integer, intent(out) :: i_reaction ! index in nuc % reactions array
+    integer, intent(out) :: i_product  ! index in reaction % products array
+
+    integer :: i_grid
+    integer :: i_temp
+    integer :: threshold
+    integer :: last_valid_reaction
+    integer :: last_valid_product
+    real(8) :: f
+    real(8) :: prob
+    real(8) :: cutoff
+    real(8) :: yield
+
+    ! Get pointer to nuclide
+    associate (nuc => nuclides(i_nuclide))
+
+      ! Get grid index and interpolation factor and sample photon production cdf
+      i_temp = micro_xs(i_nuclide) % index_temp
+      i_grid = micro_xs(i_nuclide) % index_grid
+      f      = micro_xs(i_nuclide) % interp_factor
+      cutoff = prn() * micro_xs(i_nuclide) % photon_prod
+      prob   = ZERO
+
+      ! Loop through each reaction type
+      REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
+        associate (rx => nuc % reactions(i_reaction))
+          threshold = rx % xs(i_temp) % threshold
+
+          ! if energy is below threshold for this reaction, skip it
+          if (i_grid < threshold) cycle
+
+          do i_product = 1, size(rx % products)
+            if (rx % products(i_product) % particle == PHOTON) then
+              ! add to cumulative probability
+              yield = rx % products(i_product) % yield % evaluate(E)
+              prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) &
+                   + f*(rx % xs(i_temp) % value(i_grid - threshold + 2))) * yield
+
+              if (prob > cutoff) return
+              last_valid_reaction = i_reaction
+              last_valid_product = i_product
+            end if
+          end do
+        end associate
+      end do REACTION_LOOP
+    end associate
+
+    i_reaction = last_valid_reaction
+    i_product = last_valid_product
+
+  end subroutine sample_photon_product
+
 !===============================================================================
 ! ABSORPTION
 !===============================================================================
@@ -1145,6 +1503,9 @@ contains
       ! Bank source neutrons by copying particle data
       bank_array(i) % xyz = p % coord(1) % xyz
 
+      ! Set particle as neutron
+      bank_array(i) % particle = NEUTRON
+
       ! Set weight of fission bank site
       bank_array(i) % wgt = ONE/weight
 
@@ -1241,7 +1602,7 @@ contains
         call rxn % products(1 + group) % sample(E_in, site % E, mu)
 
         ! resample if energy is greater than maximum neutron energy
-        if (site % E < energy_max_neutron) exit
+        if (site % E < energy_max(NEUTRON)) exit
 
         ! check for large number of resamples
         n_sample = n_sample + 1
@@ -1265,7 +1626,7 @@ contains
         call rxn % products(1) % sample(E_in, site % E, mu)
 
         ! resample if energy is greater than maximum neutron energy
-        if (site % E < energy_max_neutron) exit
+        if (site % E < energy_max(NEUTRON)) exit
 
         ! check for large number of resamples
         n_sample = n_sample + 1
@@ -1334,7 +1695,8 @@ contains
     if (mod(yield, ONE) == ZERO) then
       ! If yield is integral, create exactly that many secondary particles
       do i = 1, nint(yield) - 1
-        call p % create_secondary(p % coord(1) % uvw, NEUTRON, run_CE=.true.)
+        call p % create_secondary(p % coord(1) % uvw, p % E, &
+             NEUTRON, run_CE=.true.)
       end do
     else
       ! Otherwise, change weight of particle based on yield
@@ -1343,4 +1705,50 @@ contains
 
   end subroutine inelastic_scatter
 
+!===============================================================================
+! SAMPLE_SECONDARY_PHOTONS
+!===============================================================================
+
+  subroutine sample_secondary_photons(p, i_nuclide)
+    type(Particle), intent(inout) :: p
+    integer, intent(in)           :: i_nuclide
+
+    integer :: i_reaction   ! index in nuc % reactions array
+    integer :: i_product    ! index in nuc % reactions % products array
+
+    real(8) :: nu_t
+    real(8) :: mu
+    real(8) :: E
+    real(8) :: uvw(3)
+    integer :: nu
+    integer :: i
+
+    ! Sample the number of photons produced
+    nu_t =  p % wgt * micro_xs(i_nuclide) % photon_prod / &
+         micro_xs(i_nuclide) % total
+    if (prn() > nu_t - int(nu_t)) then
+      nu = int(nu_t)
+    else
+      nu = int(nu_t) + 1
+    end if
+
+    ! Sample each secondary photon
+    do i = 1, nu
+
+      ! Sample the reaction and product
+      call sample_photon_product(i_nuclide, p % E, i_reaction, i_product)
+
+      ! Sample the outgoing energy and angle
+      call nuclides(i_nuclide) % reactions(i_reaction) % products(i_product) &
+           % sample(p % E, E, mu)
+
+      ! Sample the new direction
+      uvw = rotate_angle(p % coord(1) % uvw, mu)
+
+      ! Create the secondary photon
+      call p % create_secondary(uvw, E, PHOTON, run_CE=.true.)
+    end do
+
+  end subroutine sample_secondary_photons
+
 end module physics

src/physics_mg.F90

@@ -237,6 +237,9 @@ contains
       ! Bank source neutrons by copying particle data
       bank_array(i) % xyz = p % coord(1) % xyz
 
+      ! Set particle as neutron
+      bank_array(i) % particle = NEUTRON
+
       ! Set weight of fission bank site
       bank_array(i) % wgt = ONE/weight
 

src/random_lcg.cpp

@@ -6,12 +6,13 @@ namespace openmc {
 
 
 // Constants
-extern "C" const int N_STREAMS         {5};
+extern "C" const int N_STREAMS         {6};
 extern "C" const int STREAM_TRACKING   {0};
 extern "C" const int STREAM_TALLIES    {1};
 extern "C" const int STREAM_SOURCE     {2};
 extern "C" const int STREAM_URR_PTABLE {3};
 extern "C" const int STREAM_VOLUME     {4};
+extern "C" const int STREAM_PHOTON     {5};
 
 // Starting seed
 int64_t seed {1};

src/random_lcg.h

@@ -16,6 +16,7 @@ extern "C" const int STREAM_TALLIES;
 extern "C" const int STREAM_SOURCE;
 extern "C" const int STREAM_URR_PTABLE;
 extern "C" const int STREAM_VOLUME;
+extern "C" const int STREAM_PHOTON;
 
 //==============================================================================
 //! Generate a pseudo-random number using a linear congruential generator.

src/relaxng/settings.rnc

@@ -66,6 +66,7 @@ element settings {
   element source {
     grammar {
       start =
+        (element particle { xsd:string } | attribute particle { xsd:string })? &
         (element strength { xsd:double } | attribute strength { xsd:double })? &
         (element file { xsd:string } | attribute file { xsd:string })? &
         element space {

src/relaxng/settings.rng

@@ -274,6 +274,16 @@
         <grammar>
           <start>
             <interleave>
+              <optional>
+                <choice>
+                  <element name="particle">
+                    <data type="string"/>
+                  </element>
+                  <attribute name="particle">
+                    <data type="string"/>
+                  </attribute>
+                </choice>
+              </optional>
               <optional>
                 <choice>
                   <element name="strength">

src/settings.F90

@@ -26,6 +26,9 @@ module settings
 
   integer :: n_log_bins  ! number of bins for logarithmic grid
 
+  logical :: photon_transport = .false.
+  integer :: electron_treatment = ELECTRON_TTB
+
   ! ============================================================================
   ! MULTI-GROUP CROSS SECTION RELATED VARIABLES
 
@@ -71,7 +74,7 @@ module settings
   ! Variance reduction settins
   logical :: survival_biasing = .false.
   real(8) :: weight_cutoff = 0.25_8
-  real(8) :: energy_cutoff = ZERO
+  real(8) :: energy_cutoff(4) = [ZERO, 1000.0_8, ZERO, ZERO]
   real(8) :: weight_survive = ONE
 
   ! Mode to run in (fixed source, eigenvalue, plotting, etc)

src/simulation.F90

@@ -26,6 +26,7 @@ module simulation
                              print_batch_keff, print_generation, print_runtime, &
                              print_results, print_overlap_check, write_tallies
   use particle_header, only: Particle
+  use photon_header,   only: micro_photon_xs, n_elements
   use random_lcg,      only: set_particle_seed
   use settings
   use simulation_header
@@ -428,6 +429,7 @@ contains
 !$omp parallel
     ! Allocate array for microscopic cross section cache
     allocate(micro_xs(n_nuclides))
+    allocate(micro_photon_xs(n_elements))
 
     ! Allocate array for matching filter bins
     allocate(filter_matches(n_filters))
@@ -499,7 +501,7 @@ contains
       deallocate(materials(i) % mat_nuclide_index)
     end do
 !$omp parallel
-    deallocate(micro_xs, filter_matches)
+    deallocate(micro_xs, micro_photon_xs, filter_matches)
 !$omp end parallel
 
     ! Increment total number of generations

src/source_header.F90

@@ -9,8 +9,9 @@ module source_header
   use error
   use geometry, only: find_cell
   use material_header, only: materials
-  use nuclide_header, only: energy_min_neutron, energy_max_neutron
+  use nuclide_header, only: energy_min, energy_max
   use particle_header, only: Particle
+  use settings, only: photon_transport
   use string, only: to_lower
   use xml_interface
 
@@ -29,6 +30,7 @@ module source_header
 !===============================================================================
 
   type, public :: SourceDistribution
+    integer :: particle ! particle type
     real(8) :: strength = ONE ! source strength
     class(SpatialDistribution),    allocatable :: space  ! spatial distribution
     class(UnitSphereDistribution), allocatable :: angle  ! angle distribution
@@ -53,11 +55,27 @@ contains
 
     integer :: n
     logical :: file_exists
-    character(MAX_WORD_LEN) :: type
+    character(MAX_WORD_LEN) :: type, temp_str
     type(XMLNode) :: node_space
     type(XMLNode) :: node_angle
     type(XMLNode) :: node_dist
 
+    ! Check for particle type
+    if (check_for_node(node, "particle")) then
+      call get_node_value(node, "particle", temp_str)
+      select case (to_lower(temp_str))
+      case ('neutron')
+        this % particle = NEUTRON
+      case ('photon')
+        this % particle = PHOTON
+        photon_transport = .true.
+      case default
+        call fatal_error('Unknown source particle type: ' // trim(temp_str))
+      end select
+    else
+      this % particle = NEUTRON
+    end if
+
     ! Check for source strength
     if (check_for_node(node, "strength")) then
       call get_node_value(node, "strength", this % strength)
@@ -231,6 +249,9 @@ contains
       ! Set particle defaults
       call p % initialize()
 
+      ! Set particle type
+      site % particle = this % particle
+
       ! Sample spatial distribution
       site % xyz(:) = this % space % sample()
 
@@ -274,13 +295,13 @@ contains
     ! Sample angle
     site % uvw(:) = this % angle % sample()
 
-    ! Check for monoenergetic source above maximum neutron energy
+    ! Check for monoenergetic source above maximum particle energy
     select type (energy => this % energy)
     type is (Discrete)
-      if (any(energy % x > energy_max_neutron)) then
+      if (any(energy % x > energy_max(this % particle))) then
         call fatal_error("Source energy above range of energies of at least &
              &one cross section table")
-      else if (any(energy % x < energy_min_neutron)) then
+      else if (any(energy % x < energy_min(this % particle))) then
         call fatal_error("Source energy below range of energies of at least &
              &one cross section table")
       end if
@@ -290,8 +311,9 @@ contains
       ! Sample energy spectrum
       site % E = this % energy % sample()
 
-      ! Resample if energy falls outside minimum or maximum neutron energy
-      if (site % E < energy_max_neutron .and. site % E > energy_min_neutron) exit
+      ! Resample if energy falls outside minimum or maximum particle energy
+      if (site % E < energy_max(this % particle) .and. &
+           site % E > energy_min(this % particle)) exit
     end do
 
     ! Set delayed group

src/state_point.F90

@@ -126,6 +126,11 @@ contains
       case (MODE_EIGENVALUE)
         call write_dataset(file_id, "run_mode", "eigenvalue")
       end select
+      if (photon_transport) then
+        call write_attribute(file_id, "photon_transport", 1)
+      else
+        call write_attribute(file_id, "photon_transport", 0)
+      end if
       call write_dataset(file_id, "n_particles", n_particles)
       call write_dataset(file_id, "n_batches", n_batches)
 
@@ -678,6 +683,12 @@ contains
     case ('eigenvalue')
       run_mode = MODE_EIGENVALUE
     end select
+    call read_attribute(int_array(1), file_id, "photon_transport")
+    if (int_array(1) == 1) then
+      photon_transport = .true.
+    else
+      photon_transport = .false.
+    end if
     call read_dataset(n_particles, file_id, "n_particles")
     call read_dataset(int_array(1), file_id, "n_batches")
 

src/tallies/tally.F90

@@ -99,6 +99,7 @@ contains
     real(8) :: f                    ! interpolation factor
     real(8) :: score                ! analog tally score
     real(8) :: E                    ! particle energy
+    real(8) :: xs                   ! cross section
 
     ! Pre-collision energy of particle
     E = p % last_E
@@ -129,7 +130,12 @@ contains
           else
             score = p % last_wgt
           end if
-          score = score / material_xs % total * flux
+
+          if (p % type == NEUTRON .or. p % type == PHOTON) then
+            score = score / material_xs % total * flux
+          else
+            score = ZERO
+          end if
 
         else
           ! For flux, we need no cross section
@@ -250,6 +256,9 @@ contains
 
 
       case (SCORE_FISSION)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         if (t % estimator == ESTIMATOR_ANALOG) then
           if (survival_biasing) then
             ! No fission events occur if survival biasing is on -- need to
@@ -281,6 +290,9 @@ contains
 
 
       case (SCORE_NU_FISSION)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         if (t % estimator == ESTIMATOR_ANALOG) then
           if (survival_biasing .or. p % fission) then
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@@ -324,6 +336,9 @@ contains
 
 
       case (SCORE_PROMPT_NU_FISSION)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         if (t % estimator == ESTIMATOR_ANALOG) then
           if (survival_biasing .or. p % fission) then
             if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@@ -386,6 +401,9 @@ contains
         end if
 
       case (SCORE_DELAYED_NU_FISSION)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         ! Set the delayedgroup filter index and the number of delayed group bins
         dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
 
@@ -577,6 +595,9 @@ contains
         end if
 
       case (SCORE_DECAY_RATE)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         ! Set the delayedgroup filter index
         dg_filter = t % find_filter(FILTER_DELAYEDGROUP)
 
@@ -856,6 +877,9 @@ contains
         end if
 
       case (SCORE_KAPPA_FISSION)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         ! Determine kappa-fission cross section on the fly. The ENDF standard
         ! (ENDF-102) states that MT 18 stores the fission energy as the Q_value
         ! (fission(1))
@@ -957,7 +981,10 @@ contains
           end if
         end if
 
-      case (SCORE_FISS_Q_PROMPT)
+      case (SCORE_FISS_Q_PROMPT, SCORE_FISS_Q_RECOV)
+
+        if (material_xs % absorption == ZERO) cycle SCORE_LOOP
+
         score = ZERO
 
         if (t % estimator == ESTIMATOR_ANALOG) then
@@ -968,10 +995,15 @@ contains
             associate (nuc => nuclides(p % event_nuclide))
               if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
                    allocated(nuc % fission_q_prompt)) then
-                score = p % absorb_wgt &
-                     * nuc % fission_q_prompt % evaluate(E) &
+                if (score_bin == SCORE_FISS_Q_PROMPT) then
+                  xs = nuc % fission_q_prompt % evaluate(p % last_E)
+                else if (score_bin == SCORE_FISS_Q_RECOV) then
+                  xs = nuc % fission_q_recov % evaluate(p % last_E)
+                end if
+
+                score = p % absorb_wgt * xs * flux &
                      * micro_xs(p % event_nuclide) % fission &
-                     / micro_xs(p % event_nuclide) % absorption * flux
+                     / micro_xs(p % event_nuclide) % absorption
               end if
             end associate
           else
@@ -980,86 +1012,56 @@ contains
             ! All fission events will contribute, so again we can use
             ! particle's weight entering the collision as the estimate for
             ! the fission energy production rate
-            associate (nuc => nuclides(p % event_nuclide))
-              if (allocated(nuc % fission_q_prompt)) then
-                score = p % last_wgt &
-                     * nuc % fission_q_prompt % evaluate(E) &
-                     * micro_xs(p % event_nuclide) % fission &
-                     / micro_xs(p % event_nuclide) % absorption * flux
-              end if
-            end associate
-          end if
-
-        else
-          if (i_nuclide > 0) then
-            if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then
-              score = micro_xs(i_nuclide) % fission * atom_density * flux &
-                      * nuclides(i_nuclide) % fission_q_prompt % evaluate(E)
-            end if
-          else
-            if (p % material /= MATERIAL_VOID) then
-              do l = 1, materials(p % material) % n_nuclides
-                atom_density_ = materials(p % material) % atom_density(l)
-                i_nuc = materials(p % material) % nuclide(l)
-                if (allocated(nuclides(i_nuc) % fission_q_prompt)) then
-                  score = score + micro_xs(i_nuc) % fission * atom_density_ &
-                          * flux &
-                          * nuclides(i_nuc) % fission_q_prompt % evaluate(E)
-                end if
-              end do
-            end if
-          end if
-        end if
-
-      case (SCORE_FISS_Q_RECOV)
-        score = ZERO
-
-        if (t % estimator == ESTIMATOR_ANALOG) then
-          if (survival_biasing) then
-            ! No fission events occur if survival biasing is on -- need to
-            ! calculate fraction of absorptions that would have resulted in
-            ! fission scaled by Q-value
             associate (nuc => nuclides(p % event_nuclide))
               if (micro_xs(p % event_nuclide) % absorption > ZERO .and. &
-                   allocated(nuc % fission_q_recov)) then
-                score = p % absorb_wgt &
-                     * nuc % fission_q_recov % evaluate(E) &
+                   allocated(nuc % fission_q_prompt)) then
+
+                if (score_bin == SCORE_FISS_Q_PROMPT) then
+                  xs = nuc % fission_q_prompt % evaluate(p % last_E)
+                else if (score_bin == SCORE_FISS_Q_RECOV) then
+                  xs = nuc % fission_q_recov % evaluate(p % last_E)
+                end if
+
+                score = p % last_wgt * xs * flux &
                      * micro_xs(p % event_nuclide) % fission &
-                     / micro_xs(p % event_nuclide) % absorption * flux
-              end if
-            end associate
-          else
-            ! Skip any non-absorption events
-            if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
-            ! All fission events will contribute, so again we can use
-            ! particle's weight entering the collision as the estimate for
-            ! the fission energy production rate
-            associate (nuc => nuclides(p % event_nuclide))
-              if (allocated(nuc % fission_q_recov)) then
-                score = p % last_wgt &
-                     * nuc % fission_q_recov % evaluate(E) &
-                     * micro_xs(p % event_nuclide) % fission &
-                     / micro_xs(p % event_nuclide) % absorption * flux
+                     / micro_xs(p % event_nuclide) % absorption
               end if
             end associate
           end if
 
         else
           if (i_nuclide > 0) then
-            if (allocated(nuclides(i_nuclide) % fission_q_recov)) then
-              score = micro_xs(i_nuclide) % fission * atom_density * flux &
-                      * nuclides(i_nuclide) % fission_q_recov % evaluate(E)
-            end if
+            associate (nuc => nuclides(i_nuclide))
+              if (allocated(nuc % fission_q_prompt)) then
+
+                if (score_bin == SCORE_FISS_Q_PROMPT) then
+                  xs = nuc % fission_q_prompt % evaluate(E)
+                else if (score_bin == SCORE_FISS_Q_RECOV) then
+                  xs = nuc % fission_q_recov % evaluate(E)
+                end if
+
+                score = micro_xs(i_nuclide) % fission * atom_density * flux * xs
+              end if
+            end associate
           else
             if (p % material /= MATERIAL_VOID) then
               do l = 1, materials(p % material) % n_nuclides
                 atom_density_ = materials(p % material) % atom_density(l)
                 i_nuc = materials(p % material) % nuclide(l)
-                if (allocated(nuclides(i_nuc) % fission_q_recov)) then
-                  score = score + micro_xs(i_nuc) % fission * atom_density_ &
-                          * flux &
-                          * nuclides(i_nuc) % fission_q_recov % evaluate(E)
-                end if
+
+                associate (nuc => nuclides(i_nuc))
+                  if (allocated(nuc % fission_q_prompt)) then
+
+                    if (score_bin == SCORE_FISS_Q_PROMPT) then
+                      xs = nuc % fission_q_prompt % evaluate(E)
+                    else if (score_bin == SCORE_FISS_Q_RECOV) then
+                      xs = nuc % fission_q_recov % evaluate(E)
+                    end if
+
+                    score = score + micro_xs(i_nuc) % fission * atom_density_ &
+                         * flux * xs
+                  end if
+                end associate
               end do
             end if
           end if

src/tallies/tally_filter.F90

@@ -20,6 +20,7 @@ module tally_filter
   use tally_filter_mesh
   use tally_filter_meshsurface
   use tally_filter_mu
+  use tally_filter_particle
   use tally_filter_polar
   use tally_filter_sph_harm
   use tally_filter_sptl_legendre
@@ -88,6 +89,8 @@ contains
         type_ = 'universe'
       type is (ZernikeFilter)
         type_ = 'zernike'
+      type is (ParticleFilter)
+        type_ = 'particle'
       end select
 
       ! Convert Fortran string to null-terminated C string. We assume the
@@ -147,6 +150,8 @@ contains
           allocate(MeshSurfaceFilter :: filters(index) % obj)
         case ('mu')
           allocate(MuFilter :: filters(index) % obj)
+        case ('particle')
+          allocate(ParticleFilter :: filters(index) % obj)
         case ('polar')
           allocate(PolarFilter :: filters(index) % obj)
         case ('sphericalharmonics')

src/tallies/tally_filter_particle.F90

@@ -0,0 +1,78 @@
+module tally_filter_particle
+
+  use, intrinsic :: ISO_C_BINDING
+
+  use constants,          only: ONE, MAX_LINE_LEN
+  use hdf5_interface
+  use particle_header,    only: Particle
+  use string,             only: to_str
+  use tally_filter_header
+  use xml_interface
+
+  implicit none
+  private
+
+!===============================================================================
+! PARTICLEFILTER specifies what particles can score to a tally
+!===============================================================================
+
+  type, public, extends(TallyFilter) :: ParticleFilter
+    integer, allocatable :: particles(:)
+  contains
+    procedure :: from_xml
+    procedure :: get_all_bins
+    procedure :: to_statepoint
+    procedure :: text_label
+  end type ParticleFilter
+
+contains
+
+  subroutine from_xml(this, node)
+    class(ParticleFilter), intent(inout) :: this
+    type(XMLNode), intent(in) :: node
+
+    integer :: n
+
+    ! Determine how many bins were given
+    n = node_word_count(node, "bins")
+
+    ! Allocate and store bins
+    this % n_bins = n
+    allocate(this % particles(n))
+    call get_node_array(node, "bins", this % particles)
+  end subroutine from_xml
+
+  subroutine get_all_bins(this, p, estimator, match)
+    class(ParticleFilter), intent(in)  :: this
+    type(Particle),    intent(in)  :: p
+    integer,           intent(in)  :: estimator
+    type(TallyFilterMatch), intent(inout) :: match
+
+    integer :: i
+
+    do i = 1, this % n_bins
+      if (this % particles(i) == p % type) then
+        call match % bins % push_back(i)
+        call match % weights % push_back(ONE)
+      end if
+    end do
+  end subroutine get_all_bins
+
+  subroutine to_statepoint(this, filter_group)
+    class(ParticleFilter), intent(in) :: this
+    integer(HID_T),    intent(in) :: filter_group
+
+    call write_dataset(filter_group, "type", "particle")
+    call write_dataset(filter_group, "n_bins", this % n_bins)
+    call write_dataset(filter_group, "bins", this % particles)
+  end subroutine to_statepoint
+
+  function text_label(this, bin) result(label)
+    class(ParticleFilter), intent(in) :: this
+    integer,           intent(in) :: bin
+    character(MAX_LINE_LEN)       :: label
+
+    label = "Particle " // to_str(this % particles(bin))
+  end function text_label
+
+end module tally_filter_particle

src/tallies/tally_header.F90

@@ -339,6 +339,8 @@ contains
       type is (ZernikeFilter)
         j = FILTER_ZERNIKE
         this % estimator = ESTIMATOR_COLLISION
+      type is (ParticleFilter)
+        j = FILTER_PARTICLE
       end select
       this % find_filter(j) = i
     end do

src/tracking.F90

@@ -14,7 +14,7 @@ module tracking
   use particle_header,    only: LocalCoord, Particle
   use physics,            only: collision
   use physics_mg,         only: collision_mg
-  use random_lcg,         only: prn
+  use random_lcg,         only: prn, prn_set_stream
   use settings
   use simulation_header
   use string,             only: to_str
@@ -75,6 +75,13 @@ contains
     if (active_tallies % size() > 0) call zero_flux_derivs()
 
     EVENT_LOOP: do
+      ! Set the random number stream
+      if (p % type == NEUTRON) then
+        call prn_set_stream(STREAM_TRACKING)
+      else
+        call prn_set_stream(STREAM_PHOTON)
+      end if
+
       ! Store pre-collision particle properties
       p % last_wgt = p % wgt
       p % last_E   = p % E
@@ -109,8 +116,7 @@ contains
             ! If the material is the same as the last material and the
             ! temperature hasn't changed, we don't need to lookup cross
             ! sections again.
-            call materials(p % material) % calculate_xs(p % E, p % sqrtkT, &
-                 micro_xs, nuclides, material_xs)
+            call materials(p % material) % calculate_xs(p)
           end if
         else
           ! Get the MG data
@@ -135,7 +141,9 @@ contains
            lattice_translation, next_level)
 
       ! Sample a distance to collision
-      if (material_xs % total == ZERO) then
+      if (p % type == ELECTRON .or. p % type == POSITRON) then
+        d_collision = ZERO
+      else if (material_xs % total == ZERO) then
         d_collision = INFINITY
       else
         d_collision = -log(prn()) / material_xs % total
@@ -155,7 +163,7 @@ contains
       end if
 
       ! Score track-length estimate of k-eff
-      if (run_mode == MODE_EIGENVALUE) then
+      if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
         global_tally_tracklength = global_tally_tracklength + p % wgt * &
              distance * material_xs % nu_fission
       end if
@@ -196,7 +204,7 @@ contains
         ! PARTICLE HAS COLLISION
 
         ! Score collision estimate of keff
-        if (run_mode == MODE_EIGENVALUE) then
+        if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
           global_tally_collision = global_tally_collision + p % wgt * &
                material_xs % nu_fission / material_xs % total
         end if

tests/regression_tests/energy_cutoff/inputs_true.dat

@@ -27,7 +27,7 @@
     <energy parameters="988000.0 2.249e-06" type="watt" />
   </source>
   <cutoff>
-    <energy>4.0</energy>
+    <energy_neutron>4.0</energy_neutron>
   </cutoff>
 </settings>
 <?xml version='1.0' encoding='utf-8'?>

tests/regression_tests/energy_cutoff/test.py

@@ -37,7 +37,7 @@ class EnergyCutoffTestHarness(PyAPITestHarness):
         settings_file.run_mode = 'fixed source'
         settings_file.batches = 10
         settings_file.particles = 100
-        settings_file.cutoff = {'energy': energy_cutoff}
+        settings_file.cutoff = {'energy_neutron': energy_cutoff}
         bounds = [-1, -1, -1, 1, 1, 1]
         uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
         watt_dist = openmc.stats.Watt()

tests/regression_tests/photon_production/inputs_true.dat

@@ -0,0 +1,50 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="13" material="void" region="-9 10 -11" universe="9" />
+  <cell id="14" material="13" region="-9 11 -12" universe="9" />
+  <cell id="15" material="void" region="~(-9 10 -12)" universe="9" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 1e-06" id="9" type="x-cylinder" />
+  <surface boundary="vacuum" coeffs="-1.0" id="10" type="x-plane" />
+  <surface coeffs="1.0" id="11" type="x-plane" />
+  <surface boundary="vacuum" coeffs="11.0" id="12" type="x-plane" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="13">
+    <density units="g/cm3" value="2.6989" />
+    <nuclide ao="1.0" name="Al27" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>10000</particles>
+  <batches>1</batches>
+  <source strength="1.0">
+    <space type="point">
+      <parameters>0 0 0</parameters>
+    </space>
+    <angle reference_uvw="1.0 0.0 0.0" type="monodirectional" />
+    <energy type="discrete">
+      <parameters>14.0 1.0</parameters>
+    </energy>
+  </source>
+  <electron_treatment>ttb</electron_treatment>
+  <photon_transport>true</photon_transport>
+  <cutoff>
+    <energy_photon>1000.0</energy_photon>
+  </cutoff>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>14</bins>
+  </filter>
+  <filter id="2" type="particle">
+    <bins>2</bins>
+  </filter>
+  <tally id="1">
+    <filters>1 2</filters>
+    <scores>flux</scores>
+  </tally>
+</tallies>

tests/regression_tests/photon_production/results_true.dat

@@ -0,0 +1,3 @@
+tally 1:
+sum = 1.371553E-08
+sum_sq = 1.881158E-16

tests/regression_tests/photon_production/test.py

@@ -0,0 +1,70 @@
+from math import pi
+
+import numpy as np
+import openmc
+
+from tests.testing_harness import PyAPITestHarness
+
+
+class SourceTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        mat = openmc.Material()
+        mat.set_density('g/cm3', 2.6989)
+        mat.add_nuclide('Al27', 1.0)
+        materials = openmc.Materials([mat])
+        materials.export_to_xml()
+
+        cyl = openmc.XCylinder(boundary_type='vacuum', R=1.0e-6)
+        x_plane_left = openmc.XPlane(boundary_type='vacuum', x0=-1.0)
+        x_plane_center = openmc.XPlane(boundary_type='transmission', x0=1.0)
+        x_plane_right = openmc.XPlane(boundary_type='vacuum', x0=11.0)
+
+        inner_cyl_left = openmc.Cell()
+        inner_cyl_right = openmc.Cell()
+        outer_cyl = openmc.Cell()
+
+        inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center
+        inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
+        outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
+        inner_cyl_right.fill = mat
+        geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
+        geometry.export_to_xml()
+
+        source = openmc.Source()
+        source.space = openmc.stats.Point((0,0,0))
+        source.angle = openmc.stats.Monodirectional()
+        source.energy = openmc.stats.Discrete([14.0], [1.0])
+        source.particle = 'neutron'
+
+        settings = openmc.Settings()
+        settings.particles = 10000
+        settings.run_mode = 'fixed source'
+        settings.batches = 1
+        settings.photon_transport = True
+        settings.electron_treatment = 'ttb'
+        settings.cutoff = {'energy_photon' : 1000.0}
+        settings.source = source
+        settings.export_to_xml()
+
+        cell_filter = openmc.CellFilter(inner_cyl_right)
+        particle_filter = openmc.ParticleFilter('photon')
+        tally = openmc.Tally()
+        tally.filters = [cell_filter, particle_filter]
+        tally.scores = ['flux']
+        tallies = openmc.Tallies([tally])
+        tallies.export_to_xml()
+
+    def _get_results(self):
+        with openmc.StatePoint(self._sp_name) as sp:
+            outstr = ''
+            t = sp.get_tally()
+            outstr += 'tally {}:\n'.format(t.id)
+            outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0])
+            outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0])
+
+            return outstr
+
+
+def test_source():
+    harness = SourceTestHarness('statepoint.1.h5')
+    harness.main()

tests/regression_tests/photon_source/inputs_true.dat

@@ -0,0 +1,45 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="13" material="13" region="-9" universe="9" />
+  <surface boundary="reflective" coeffs="0.0 0.0 0.0 1000000000.0" id="9" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="13">
+    <density units="g/cm3" value="0.998207" />
+    <nuclide ao="0.11187657362844" name="H1" />
+    <nuclide ao="1.7426371559999997e-05" name="H2" />
+    <nuclide ao="0.8877694078259999" name="O16" />
+    <nuclide ao="0.000336592174" name="O17" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>10000</particles>
+  <batches>1</batches>
+  <source particle="photon" strength="1.0">
+    <space type="point">
+      <parameters>0 0 0</parameters>
+    </space>
+    <angle type="isotropic" />
+    <energy type="discrete">
+      <parameters>10000000.0 1.0</parameters>
+    </energy>
+  </source>
+  <electron_treatment>ttb</electron_treatment>
+  <photon_transport>true</photon_transport>
+  <cutoff>
+    <energy_photon>1000.0</energy_photon>
+  </cutoff>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="particle">
+    <bins>2</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <scores>flux</scores>
+  </tally>
+</tallies>

tests/regression_tests/photon_source/results_true.dat

@@ -0,0 +1,3 @@
+tally 1:
+sum = 2.254985E+02
+sum_sq = 5.084955E+04

tests/regression_tests/photon_source/test.py

@@ -0,0 +1,61 @@
+from math import pi
+
+import numpy as np
+import openmc
+
+from tests.testing_harness import PyAPITestHarness
+
+
+class SourceTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+        mat = openmc.Material()
+        mat.set_density('g/cm3', 0.998207)
+        mat.add_element('H', 0.111894)
+        mat.add_element('O', 0.888106)
+        materials = openmc.Materials([mat])
+        materials.export_to_xml()
+
+        sphere = openmc.Sphere(R=1.0e9, boundary_type='reflective')
+        inside_sphere = openmc.Cell()
+        inside_sphere.region = -sphere
+        inside_sphere.fill = mat
+        geometry = openmc.Geometry([inside_sphere])
+        geometry.export_to_xml()
+ 
+        source = openmc.Source()
+        source.space = openmc.stats.Point((0, 0, 0))
+        source.angle = openmc.stats.Isotropic()
+        source.energy = openmc.stats.Discrete([10.0e6], [1.0])
+        source.particle = 'photon'
+ 
+        settings = openmc.Settings()
+        settings.particles = 10000
+        settings.batches = 1
+        settings.photon_transport = True
+        settings.electron_treatment = 'ttb'
+        settings.cutoff = {'energy_photon' : 1000.0}
+        settings.run_mode = 'fixed source'
+        settings.source = source
+        settings.export_to_xml()
+ 
+        particle_filter = openmc.ParticleFilter('photon')
+        tally = openmc.Tally()
+        tally.filters = [particle_filter]
+        tally.scores = ['flux']
+        tallies = openmc.Tallies([tally])
+        tallies.export_to_xml()
+
+    def _get_results(self):
+        with openmc.StatePoint(self._sp_name) as sp:
+            outstr = ''
+            t = sp.get_tally()
+            outstr += 'tally {}:\n'.format(t.id)
+            outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0])
+            outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0])
+
+            return outstr
+
+
+def test_source():
+    harness = SourceTestHarness('statepoint.1.h5')
+    harness.main()

tests/unit_tests/test_data_photon.py

@@ -0,0 +1,138 @@
+#!/usr/bin/env python
+
+from collections import Mapping, Callable
+import os
+
+import numpy as np
+import pandas as pd
+import pytest
+import openmc.data
+
+
+_ENDF_DATA = os.environ['OPENMC_ENDF_DATA']
+
+
+@pytest.fixture(scope='module')
+def elements_endf():
+    """Dictionary of element ENDF data indexed by atomic symbol."""
+    elements = {'H': 1, 'O': 8, 'Al': 13, 'Cu': 29, 'Ag': 47, 'U': 92, 'Pu': 94}
+    data = {}
+    for symbol, Z in elements.items():
+        p_file = 'photoat-{:03}_{}_000.endf'.format(Z, symbol)
+        p_path = os.path.join(_ENDF_DATA, 'photoat', p_file)
+        a_file = 'atom-{:03}_{}_000.endf'.format(Z, symbol)
+        a_path = os.path.join(_ENDF_DATA, 'atomic_relax', a_file)
+        data[symbol] = openmc.data.IncidentPhoton.from_endf(p_path, a_path)
+    return data
+
+
+@pytest.fixture()
+def element(request, elements_endf):
+    """Element ENDF data"""
+    return elements_endf[request.param]
+
+
+@pytest.mark.parametrize(
+    'element, atomic_number', [
+        ('Al', 13),
+        ('Cu', 29),
+        ('Pu', 94)
+    ],
+    indirect=['element']
+)
+def test_attributes(element, atomic_number):
+    assert element.atomic_number == atomic_number
+
+
+@pytest.mark.parametrize(
+    'element, subshell, binding_energy, num_electrons', [
+        ('H', 'K', 13.61, 1.0),
+        ('O', 'L3', 14.15, 2.67),
+        ('U', 'P2', 34.09, 2.0)
+    ],
+    indirect=['element']
+)
+def test_atomic_relaxation(element, subshell, binding_energy, num_electrons):
+    atom_relax = element.atomic_relaxation
+    assert isinstance(atom_relax, openmc.data.photon.AtomicRelaxation)
+    assert subshell in atom_relax.subshells
+    assert atom_relax.binding_energy[subshell] == binding_energy
+    assert atom_relax.num_electrons[subshell] == num_electrons
+
+
+@pytest.mark.parametrize('element', ['Al', 'Cu', 'Pu'], indirect=True)
+def test_transitions(element):
+    transitions = element.atomic_relaxation.transitions
+    assert transitions
+    assert isinstance(transitions, Mapping)
+    for matrix in transitions.values():
+        assert isinstance(matrix, pd.core.frame.DataFrame)
+        assert len(matrix.columns) == 4
+        assert sum(matrix['probability']) == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize('element', ['H', 'Al', 'Ag'], indirect=True)
+def test_bremsstrahlung(element):
+    brems = element.bremsstrahlung
+    assert isinstance(brems, Mapping)
+    assert np.all(np.diff(brems['electron_energy']) > 0.0)
+    assert np.all(np.diff(brems['photon_energy']) > 0.0)
+    assert brems['photon_energy'][0] == 0.0
+    assert brems['photon_energy'][-1] == 1.0
+    assert brems['dcs'].shape == (200, 30)
+
+
+@pytest.mark.parametrize(
+    'element, n_shell', [
+        ('H', 1),
+        ('O', 3),
+        ('Al', 5)
+    ],
+    indirect=['element']
+)
+def test_compton_profiles(element, n_shell):
+    profile = element.compton_profiles
+    assert profile
+    assert isinstance(profile, Mapping)
+    assert all(isinstance(x, Callable) for x in profile['J'])
+    assert all(len(x) == n_shell for x in profile.values())
+
+
+@pytest.mark.parametrize(
+    'element, reaction', [
+        ('Cu', 541),
+        ('Ag', 502),
+        ('Pu', 504)
+    ],
+    indirect=['element']
+)
+def test_reactions(element, reaction):
+    reactions = element.reactions
+    assert all(isinstance(x, openmc.data.PhotonReaction) for x in reactions.values())
+    assert reaction in reactions
+    with pytest.raises(KeyError):
+        reactions[18]
+
+
+@pytest.mark.parametrize(
+    'element, I', [
+        ('H', 19.2),
+        ('O', 95.0),
+        ('U', 890.0)
+    ],
+    indirect=['element']
+)
+def test_stopping_powers(element, I):
+    stopping_powers = element.stopping_powers
+    assert isinstance(stopping_powers, Mapping)
+    assert stopping_powers['I'] == I
+    assert np.all(np.diff(stopping_powers['energy']) > 0.0)
+    assert len(stopping_powers['s_collision']) == 200
+    assert len(stopping_powers['s_radiative']) == 200
+
+
+@pytest.mark.parametrize('element', ['Pu'], indirect=True)
+def test_export_to_hdf5(tmpdir, element):
+    filename = str(tmpdir.join('tmp.h5'))
+    element.export_to_hdf5(filename)
+    assert os.path.exists(filename)

tools/ci/travis-before-script.sh

@@ -7,11 +7,12 @@ sh -e /etc/init.d/xvfb start
 
 # Download NNDC HDF5 data
 if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
-    wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ
+    wget https://anl.box.com/shared/static/na85do11dfh0lb9utye2il5o6yaxx8hi.xz -O - | tar -C $HOME -xvJ
 fi
 
 # Download ENDF/B-VII.1 distribution
-if [[ ! -d $HOME/endf-b-vii.1/neutrons ]]; then
+ENDF=$HOME/endf-b-vii.1/
+if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
     wget https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz -O - | tar -C $HOME -xvJ
 fi