Updating my local OpenMC.

Merge branch 'develop' of https://github.com/mit-crpg/openmc into new-update
2026-07-27 21:55:41 -04:00 · 2018-05-23 19:18:36 +00:00 · 2018-05-23 19:18:36 +00:00 · 92f147d387
commit 92f147d387
parent 2452747a9b cacde822f9
37 changed files with 1678 additions and 916 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -435,6 +435,7 @@ set(LIBOPENMC_CXX_SRC
  src/initialize.cpp
  src/finalize.cpp
  src/hdf5_interface.cpp
+  src/math_functions.cpp
  src/message_passing.cpp
  src/plot.cpp
  src/random_lcg.cpp
--- a/include/openmc.h
+++ b/include/openmc.h
@ -17,6 +17,7 @@ extern "C" {
  };

  int openmc_calculate_volumes();
+  int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
  int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
  int openmc_cell_get_id(int32_t index, int32_t* id);
  int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
@ -31,7 +32,7 @@ extern "C" {
  int openmc_extend_sources(int32_t n, int32_t* index_start, int32_t* index_end);
  int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end);
  int openmc_filter_get_id(int32_t index, int32_t* id);
-  int openmc_filter_get_type(int32_t index, const char** type);
+  int openmc_filter_get_type(int32_t index, char* type);
  int openmc_filter_set_id(int32_t index, int32_t id);
  int openmc_filter_set_type(int32_t index, const char* type);
  int openmc_finalize();
--- a/openmc/capi/init.py
+++ b/openmc/capi/init.py
@ -47,3 +47,4 @@ from .mesh import *
 from .filter import *
 from .tally import *
 from .settings import settings
+from .math import *
--- a/openmc/capi/filter.py
+++ b/openmc/capi/filter.py
@ -17,11 +17,16 @@ from .mesh import Mesh
 __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
           'CellbornFilter', 'CellfromFilter', 'DistribcellFilter',
           'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter',
-           'EnergyFunctionFilter', 'MaterialFilter', 'MeshFilter',
-           'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SurfaceFilter',
-           'UniverseFilter', 'filters']
+           'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter',
+           'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
+           'SpatialLegendreFilter', 'SurfaceFilter',
+           'UniverseFilter', 'ZernikeFilter', 'filters']

 # Tally functions
+_dll.openmc_cell_filter_get_bins.argtypes = [
+    c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)]
+_dll.openmc_cell_filter_get_bins.restype = c_int
+_dll.openmc_cell_filter_get_bins.errcheck = _error_handler
 _dll.openmc_energy_filter_get_bins.argtypes = [
    c_int32, POINTER(POINTER(c_double)), POINTER(c_int32)]
 _dll.openmc_energy_filter_get_bins.restype = c_int
@ -47,6 +52,12 @@ _dll.openmc_filter_set_type.errcheck = _error_handler
 _dll.openmc_get_filter_index.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_get_filter_index.restype = c_int
 _dll.openmc_get_filter_index.errcheck = _error_handler
+_dll.openmc_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_legendre_filter_get_order.restype = c_int
+_dll.openmc_legendre_filter_get_order.errcheck = _error_handler
+_dll.openmc_legendre_filter_set_order.argtypes = [c_int32, c_int]
+_dll.openmc_legendre_filter_set_order.restype = c_int
+_dll.openmc_legendre_filter_set_order.errcheck = _error_handler
 _dll.openmc_material_filter_get_bins.argtypes = [
    c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)]
 _dll.openmc_material_filter_get_bins.restype = c_int
@ -66,7 +77,24 @@ _dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
 _dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_meshsurface_filter_set_mesh.restype = c_int
 _dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
-
+_dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_spatial_legendre_filter_get_order.restype = c_int
+_dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler
+_dll.openmc_spatial_legendre_filter_set_order.argtypes = [c_int32, c_int]
+_dll.openmc_spatial_legendre_filter_set_order.restype = c_int
+_dll.openmc_spatial_legendre_filter_set_order.errcheck = _error_handler
+_dll.openmc_sphharm_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_sphharm_filter_get_order.restype = c_int
+_dll.openmc_sphharm_filter_get_order.errcheck = _error_handler
+_dll.openmc_sphharm_filter_set_order.argtypes = [c_int32, c_int]
+_dll.openmc_sphharm_filter_set_order.restype = c_int
+_dll.openmc_sphharm_filter_set_order.errcheck = _error_handler
+_dll.openmc_zernike_filter_get_order.argtypes = [c_int32, POINTER(c_int)]
+_dll.openmc_zernike_filter_get_order.restype = c_int
+_dll.openmc_zernike_filter_get_order.errcheck = _error_handler
+_dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int]
+_dll.openmc_zernike_filter_set_order.restype = c_int
+_dll.openmc_zernike_filter_set_order.errcheck = _error_handler

 class Filter(_FortranObjectWithID):
    __instances = WeakValueDictionary()
@ -150,6 +178,13 @@ class AzimuthalFilter(Filter):
 class CellFilter(Filter):
    filter_type = 'cell'

+    @property
+    def bins(self):
+        cells = POINTER(c_int32)()
+        n = c_int32()
+        _dll.openmc_cell_filter_get_bins(self._index, cells, n)
+        return as_array(cells, (n.value,))
+

 class CellbornFilter(Filter):
    filter_type = 'cellborn'
@ -171,6 +206,25 @@ class EnergyFunctionFilter(Filter):
    filter_type = 'energyfunction'


+class LegendreFilter(Filter):
+    filter_type = 'legendre'
+
+    def __init__(self, order=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if order is not None:
+            self.order = order
+
+    @property
+    def order(self):
+        temp_order = c_int()
+        _dll.openmc_legendre_filter_get_order(self._index, temp_order)
+        return temp_order.value
+
+    @order.setter
+    def order(self, order):
+        _dll.openmc_legendre_filter_set_order(self._index, order)
+
+
 class MaterialFilter(Filter):
    filter_type = 'material'

@ -241,6 +295,44 @@ class PolarFilter(Filter):
    filter_type = 'polar'


+class SphericalHarmonicsFilter(Filter):
+    filter_type = 'sphericalharmonics'
+
+    def __init__(self, order=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if order is not None:
+            self.order = order
+
+    @property
+    def order(self):
+        temp_order = c_int()
+        _dll.openmc_sphharm_filter_get_order(self._index, temp_order)
+        return temp_order.value
+
+    @order.setter
+    def order(self, order):
+        _dll.openmc_sphharm_filter_set_order(self._index, order)
+
+
+class SpatialLegendreFilter(Filter):
+    filter_type = 'spatiallegendre'
+
+    def __init__(self, order=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if order is not None:
+            self.order = order
+
+    @property
+    def order(self):
+        temp_order = c_int()
+        _dll.openmc_spatial_legendre_filter_get_order(self._index, temp_order)
+        return temp_order.value
+
+    @order.setter
+    def order(self, order):
+        _dll.openmc_spatial_legendre_filter_set_order(self._index, order)
+
+
 class SurfaceFilter(Filter):
    filter_type = 'surface'

@ -249,6 +341,25 @@ class UniverseFilter(Filter):
    filter_type = 'universe'


+class ZernikeFilter(Filter):
+    filter_type = 'zernike'
+
+    def __init__(self, order=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if order is not None:
+            self.order = order
+
+    @property
+    def order(self):
+        temp_order = c_int()
+        _dll.openmc_zernike_filter_get_order(self._index, temp_order)
+        return temp_order.value
+
+    @order.setter
+    def order(self, order):
+        _dll.openmc_zernike_filter_set_order(self._index, order)
+
+
 _FILTER_TYPE_MAP = {
    'azimuthal': AzimuthalFilter,
    'cell': CellFilter,
@ -259,13 +370,17 @@ _FILTER_TYPE_MAP = {
    'energy': EnergyFilter,
    'energyout': EnergyoutFilter,
    'energyfunction': EnergyFunctionFilter,
+    'legendre': LegendreFilter,
    'material': MaterialFilter,
    'mesh': MeshFilter,
    'meshsurface': MeshSurfaceFilter,
    'mu': MuFilter,
    'polar': PolarFilter,
+    'sphericalharmonics': SphericalHarmonicsFilter,
+    'spatiallegendre': SpatialLegendreFilter,
    'surface': SurfaceFilter,
    'universe': UniverseFilter,
+    'zernike': ZernikeFilter
 }


--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@ -0,0 +1,250 @@
+from ctypes import (c_int, c_double, POINTER)
+
+import numpy as np
+from numpy.ctypeslib import ndpointer
+
+from . import _dll
+
+
+_dll.t_percentile_c.restype = c_double
+_dll.t_percentile_c.argtypes = [c_double, c_int]
+
+_dll.calc_pn_c.restype = None
+_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
+
+_dll.evaluate_legendre_c.restype = c_double
+_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
+
+_dll.calc_rn_c.restype = None
+_dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
+
+_dll.calc_zn_c.restype = None
+_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
+
+_dll.rotate_angle_c.restype = None
+_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
+                                POINTER(c_double)]
+_dll.maxwell_spectrum_c.restype = c_double
+_dll.maxwell_spectrum_c.argtypes = [c_double]
+
+_dll.watt_spectrum_c.restype = c_double
+_dll.watt_spectrum_c.argtypes = [c_double, c_double]
+
+_dll.broaden_wmp_polynomials_c.restype = None
+_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
+                                           ndpointer(c_double)]
+
+
+def t_percentile(p, df):
+    """ Calculate the percentile of the Student's t distribution with a
+    specified probability level and number of degrees of freedom
+
+    Parameters
+    ----------
+    p : float
+        Probability level
+    df : int
+        Degrees of freedom
+
+    Returns
+    -------
+    float
+        Corresponding t-value
+
+    """
+
+    return _dll.t_percentile_c(p, df)
+
+
+def calc_pn(n, x):
+    """ Calculate the n-th order Legendre polynomial at the value of x.
+
+    Parameters
+    ----------
+    n : int
+        Legendre order
+    x : float
+        Independent variable to evaluate the Legendre at
+
+    Returns
+    -------
+    float
+        Corresponding Legendre polynomial result
+
+    """
+
+    pnx = np.empty(n + 1, dtype=np.float64)
+    _dll.calc_pn_c(n, x, pnx)
+    return pnx
+
+
+def evaluate_legendre(data, x):
+    """ Finds the value of f(x) given a set of Legendre coefficients
+    and the value of x.
+
+    Parameters
+    ----------
+    data : iterable of float
+        Legendre coefficients
+    x : float
+        Independent variable to evaluate the Legendre at
+
+    Returns
+    -------
+    float
+        Corresponding Legendre expansion result
+
+    """
+
+    data_arr = np.array(data, dtype=np.float64)
+    return _dll.evaluate_legendre_c(len(data),
+                                    data_arr.ctypes.data_as(POINTER(c_double)),
+                                    x)
+
+
+def calc_rn(n, uvw):
+    """ Calculate the n-th order real Spherical Harmonics for a given angle;
+    all Rn,m values are provided for all n (where -n <= m <= n).
+
+    Parameters
+    ----------
+    n : int
+        Harmonics order
+    uvw : iterable of float
+        Independent variable to evaluate the Legendre at
+
+    Returns
+    -------
+    numpy.ndarray
+        Corresponding real harmonics value
+
+    """
+
+    num_nm = (n + 1) * (n + 1)
+    rn = np.empty(num_nm, dtype=np.float64)
+    uvw_arr = np.array(uvw, dtype=np.float64)
+    _dll.calc_rn_c(n, uvw_arr, rn)
+    return rn
+
+
+def calc_zn(n, rho, phi):
+    """ Calculate the n-th order modified Zernike polynomial moment for a
+    given angle (rho, theta) location in the unit disk. The normalization of
+    the polynomials is such that the integral of Z_pq*Z_pq over the unit disk
+    is exactly pi
+
+    Parameters
+    ----------
+    n : int
+        Maximum order
+    rho : float
+        Radial location in the unit disk
+    phi : float
+        Theta (radians) location in the unit disk
+
+    Returns
+    -------
+    numpy.ndarray
+        Corresponding resulting list of coefficients
+
+    """
+
+    num_bins = ((n + 1) * (n + 2)) // 2
+    zn = np.zeros(num_bins, dtype=np.float64)
+    _dll.calc_zn_c(n, rho, phi, zn)
+    return zn
+
+
+def rotate_angle(uvw0, mu, phi=None):
+    """ Rotates direction cosines through a polar angle whose cosine is
+    mu and through an azimuthal angle sampled uniformly.
+
+    Parameters
+    ----------
+    uvw0 : iterable of float
+        Original direction cosine
+    mu : float
+        Polar angle cosine to rotate
+    phi : float, optional
+        Azimuthal angle; if None, one will be sampled uniformly
+
+    Returns
+    -------
+    numpy.ndarray
+        Rotated direction cosine
+
+    """
+
+    uvw0_arr = np.array(uvw0, dtype=np.float64)
+
+    if phi is None:
+        _dll.rotate_angle_c(uvw0_arr, mu, None)
+    else:
+        _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi))
+    uvw = uvw0_arr
+
+    return uvw
+
+
+def maxwell_spectrum(T):
+    """ Samples an energy from the Maxwell fission distribution based
+    on a direct sampling scheme.
+
+    Parameters
+    ----------
+    T : float
+        Spectrum parameter
+
+    Returns
+    -------
+    float
+        Sampled outgoing energy
+
+    """
+
+    return _dll.maxwell_spectrum_c(T)
+
+
+def watt_spectrum(a, b):
+    """ Samples an energy from the Watt energy-dependent fission spectrum.
+
+    Parameters
+    ----------
+    a : float
+        Spectrum parameter a
+    b : float
+        Spectrum parameter b
+
+    Returns
+    -------
+    float
+        Sampled outgoing energy
+
+    """
+
+    return _dll.watt_spectrum_c(a, b)
+
+
+def broaden_wmp_polynomials(E, dopp, n):
+    """ Doppler broadens the windowed multipole curvefit.  The curvefit is a
+    polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ...
+
+    Parameters
+    ----------
+    E : float
+        Energy to evaluate at
+    dopp : float
+        sqrt(atomic weight ratio / kT), with kT given in eV
+    n : int
+        Number of components to the polynomial
+
+    Returns
+    -------
+    numpy.ndarray
+        Resultant leading coefficients
+
+    """
+
+    factors = np.zeros(n, dtype=np.float64)
+    _dll.broaden_wmp_polynomials_c(E, dopp, n, factors)
+    return factors
--- a/openmc/checkvalue.py
+++ b/openmc/checkvalue.py
@ -265,8 +265,9 @@ def check_filetype_version(obj, expected_type, expected_version):
        if this_version[0] != expected_version:
            raise IOError('{} file has a version of {} which is not '
                          'consistent with the version expected by OpenMC, {}'
-                          .format(this_filetype, '.'.join(this_version),
-                                  expected_version))
+                          .format(this_filetype, 
+                              '.'.join(str(v) for v in this_version),
+                              expected_version))
    except AttributeError:
        raise IOError('Could not read {} file. This most likely means the {} '
                      'file was produced by a different version of OpenMC than '
--- a/openmc/deplete/results_list.py
+++ b/openmc/deplete/results_list.py
@ -43,8 +43,8 @@ class ResultsList(list):
            Total number of atoms for specified nuclide

        """
-        time = np.empty_like(self)
-        concentration = np.empty_like(self)
+        time = np.empty_like(self, dtype=float)
+        concentration = np.empty_like(self, dtype=float)

        # Evaluate value in each region
        for i, result in enumerate(self):
@ -73,8 +73,8 @@ class ResultsList(list):
            Array of reaction rates

        """
-        time = np.empty_like(self)
-        rate = np.empty_like(self)
+        time = np.empty_like(self, dtype=float)
+        rate = np.empty_like(self, dtype=float)

        # Evaluate value in each region
        for i, result in enumerate(self):
@ -94,8 +94,8 @@ class ResultsList(list):
            k-eigenvalue at each time

        """
-        time = np.empty_like(self)
-        eigenvalue = np.empty_like(self)
+        time = np.empty_like(self, dtype=float)
+        eigenvalue = np.empty_like(self, dtype=float)

        # Get time/eigenvalue at each point
        for i, result in enumerate(self):
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@ -546,10 +546,15 @@ class RectLattice(Lattice):
        """
        if self.ndim == 2:
            nx, ny = self.shape
-            return np.broadcast(*np.ogrid[:nx, :ny])
+            for iy in range(ny):
+                for ix in range(nx):
+                    yield (ix, iy)
        else:
            nx, ny, nz = self.shape
-            return np.broadcast(*np.ogrid[:nx, :ny, :nz])
+            for iz in range(nz):
+                for iy in range(ny):
+                    for ix in range(nx):
+                        yield (ix, iy, iz)

    @property
    def lower_left(self):
--- a/src/api.F90
+++ b/src/api.F90
@ -10,6 +10,7 @@ module openmc_api
  use geometry_header
  use hdf5_interface
  use material_header
+  use math
  use mesh_header
  use message_passing
  use nuclide_header
@ -34,6 +35,7 @@ module openmc_api

  private
  public :: openmc_calculate_volumes
+  public :: openmc_cell_filter_get_bins
  public :: openmc_cell_get_id
  public :: openmc_cell_get_fill
  public :: openmc_cell_set_fill
--- a/src/distribution_multivariate.F90
+++ b/src/distribution_multivariate.F90
@ -119,7 +119,7 @@ contains
    else
      ! Sample azimuthal angle
      phi = this % phi % sample()
-      uvw(:) = rotate_angle(this % reference_uvw, mu, phi)
+      uvw = rotate_angle(this % reference_uvw, mu, phi)
    end if
  end function polar_azimuthal_sample

--- a/src/faddeeva/Faddeeva.h
+++ b/src/faddeeva/Faddeeva.h
@ -1,5 +1,5 @@
 /* Copyright (c) 2012 Massachusetts Institute of Technology
- * 
+ *
 * Permission is hereby granted, free of charge, to any person obtaining
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
@ -7,17 +7,17 @@
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
- * 
+ *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
- * 
+ *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
- * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

 /* Available at: http://ab-initio.mit.edu/Faddeeva
--- a/src/geometry.F90
+++ b/src/geometry.F90
@ -882,9 +882,9 @@ contains
        type is (RectLattice)

          ! Loop over lattice coordinates
-          do j = 1, lat % n_cells(1)
+          do m = 1, lat % n_cells(3)
            do k = 1, lat % n_cells(2)
-              do m = 1, lat % n_cells(3)
+              do j = 1, lat % n_cells(1)
                lat % offset(map, j, k, m) = offset
                next_univ => universes(lat % universes(j, k, m))
                offset = offset + &
@ -999,9 +999,9 @@ contains
        type is (RectLattice)

          ! Loop over lattice coordinates
-          do j = 1, lat % n_cells(1)
+          do m = 1, lat % n_cells(3)
            do k = 1, lat % n_cells(2)
-              do m = 1, lat % n_cells(3)
+              do j = 1, lat % n_cells(1)
                next_univ => universes(lat % universes(j, k, m))

                ! Found target - stop since target cannot contain itself
@ -1091,9 +1091,9 @@ contains
            type is (RectLattice)

              ! Loop over lattice coordinates
-              do j = 1, lat % n_cells(1)
+              do m = 1, lat % n_cells(3)
                do k = 1, lat % n_cells(2)
-                  do m = 1, lat % n_cells(3)
+                  do j = 1, lat % n_cells(1)
                    call count_instance(universes(lat % universes(j, k, m)))
                  end do
                end do
@ -1166,9 +1166,9 @@ contains
        type is (RectLattice)

          ! Loop over lattice coordinates
-          do j = 1, lat % n_cells(1)
+          do m = 1, lat % n_cells(3)
            do k = 1, lat % n_cells(2)
-              do m = 1, lat % n_cells(3)
+              do j = 1, lat % n_cells(1)
                next_univ => universes(lat % universes(j, k, m))
                levels_below = max(levels_below, maximum_levels(next_univ))
              end do
--- a/src/math.F90
+++ b/src/math.F90
@ -6,6 +6,18 @@ module math
  use random_lcg, only: prn

  implicit none
+  private
+  public :: t_percentile
+  public :: calc_pn
+  public :: calc_rn
+  public :: calc_zn
+  public :: evaluate_legendre
+  public :: rotate_angle
+  public :: maxwell_spectrum
+  public :: watt_spectrum
+  public :: faddeeva
+  public :: w_derivative
+  public :: broaden_wmp_polynomials

 !===============================================================================
 ! FADDEEVA_W evaluates the scaled complementary error function.  This
@ -13,6 +25,86 @@ module math
 !===============================================================================

  interface
+
+    pure function t_percentile(p, df) bind(C, name='t_percentile_c') &
+         result(t)
+      use ISO_C_BINDING
+      implicit none
+      real(C_DOUBLE), value, intent(in) :: p
+      integer(C_INT), value, intent(in) :: df
+      real(C_DOUBLE) :: t
+    end function t_percentile
+
+    pure subroutine calc_pn(n, x, pnx) bind(C, name='calc_pn_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), value, intent(in) :: x
+      real(C_DOUBLE), intent(out) :: pnx(n + 1)
+    end subroutine calc_pn
+
+    pure function evaluate_legendre_c_intfc(n, data, x) &
+         bind(C, name='evaluate_legendre_c') result(val)
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(in) :: data(n)
+      real(C_DOUBLE), value, intent(in) :: x
+      real(C_DOUBLE) :: val
+    end function evaluate_legendre_c_intfc
+
+    pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(in)  :: uvw(3)
+      real(C_DOUBLE), intent(out) :: rn(2 * n + 1)
+    end subroutine calc_rn
+
+    pure subroutine calc_zn(n, rho, phi, zn) bind(C, name='calc_zn_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), value, intent(in) :: rho
+      real(C_DOUBLE), value, intent(in) :: phi
+      real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
+    end subroutine calc_zn
+
+    subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
+      use ISO_C_BINDING
+      implicit none
+      real(C_DOUBLE), intent(inout) :: uvw(3)
+      real(C_DOUBLE), value, intent(in)    :: mu
+      real(C_DOUBLE), optional, intent(in) :: phi
+    end subroutine rotate_angle_c_intfc
+
+    function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') &
+         result(E_out)
+      use ISO_C_BINDING
+      implicit none
+      real(C_DOUBLE), value, intent(in) :: T
+      real(C_DOUBLE) :: E_out
+    end function maxwell_spectrum
+
+    function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') &
+         result(E_out)
+      use ISO_C_BINDING
+      implicit none
+      real(C_DOUBLE), value, intent(in) :: a
+      real(C_DOUBLE), value, intent(in) :: b
+      real(C_DOUBLE) :: E_out
+    end function watt_spectrum
+
+    subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
+           bind(C, name='broaden_wmp_polynomials_c')
+      use ISO_C_BINDING
+      implicit none
+      real(C_DOUBLE), value, intent(in) :: E
+      real(C_DOUBLE), value, intent(in) :: dopp
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), intent(inout) :: factors(n)
+    end subroutine broaden_wmp_polynomials
+
    function faddeeva_w(z, relerr) bind(C, name='Faddeeva_w') result(w)
      use ISO_C_BINDING
      implicit none
@ -24,696 +116,17 @@ module math

 contains

-!===============================================================================
-! NORMAL_PERCENTILE calculates the percentile of the standard normal
-! distribution with a specified probability level
-!===============================================================================
-
-  elemental function normal_percentile(p) result(z)
-
-    real(8), intent(in) :: p ! probability level
-    real(8)             :: z ! corresponding z-value
-
-    real(8)            :: q
-    real(8)            :: r
-    real(8), parameter :: p_low  = 0.02425_8
-    real(8), parameter :: a(6) = (/ &
-         -3.969683028665376e1_8, 2.209460984245205e2_8, -2.759285104469687e2_8, &
-         1.383577518672690e2_8, -3.066479806614716e1_8, 2.506628277459239e0_8 /)
-    real(8), parameter :: b(5) = (/ &
-         -5.447609879822406e1_8, 1.615858368580409e2_8, -1.556989798598866e2_8, &
-         6.680131188771972e1_8, -1.328068155288572e1_8 /)
-    real(8), parameter :: c(6) = (/ &
-         -7.784894002430293e-3_8, -3.223964580411365e-1_8, -2.400758277161838_8, &
-         -2.549732539343734_8, 4.374664141464968_8, 2.938163982698783_8 /)
-    real(8), parameter :: d(4) = (/ &
-         7.784695709041462e-3_8, 3.224671290700398e-1_8, &
-         2.445134137142996_8,    3.754408661907416_8 /)
-
-    ! The rational approximation used here is from an unpublished work at
-    ! http://home.online.no/~pjacklam/notes/invnorm/
-
-    if (p < p_low) then
-      ! Rational approximation for lower region.
-
-      q = sqrt(-TWO*log(p))
-      z = (((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
-           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE)
-
-    elseif (p <= ONE - p_low) then
-      ! Rational approximation for central region
-
-      q = p - HALF
-      r = q*q
-      z = (((((a(1)*r + a(2))*r + a(3))*r + a(4))*r + a(5))*r + a(6))*q / &
-           (((((b(1)*r + b(2))*r + b(3))*r + b(4))*r + b(5))*r + ONE)
-
-    else
-      ! Rational approximation for upper region
-
-      q = sqrt(-TWO*log(ONE - p))
-      z = -(((((c(1)*q + c(2))*q + c(3))*q + c(4))*q + c(5))*q + c(6)) / &
-           ((((d(1)*q + d(2))*q + d(3))*q + d(4))*q + ONE)
-    endif
-
-    ! Refinement based on Newton's method
-#ifndef NO_F2008
-    z = z - (HALF * erfc(-z/sqrt(TWO)) - p) * sqrt(TWO*PI) * exp(HALF*z*z)
-#endif
-
-  end function normal_percentile
-
-!===============================================================================
-! T_PERCENTILE calculates the percentile of the Student's t distribution with a
-! specified probability level and number of degrees of freedom
-!===============================================================================
-
-  elemental function t_percentile(p, df) result(t)
-
-    real(8), intent(in) :: p  ! probability level
-    integer, intent(in) :: df ! degrees of freedom
-    real(8)             :: t  ! corresponding t-value
-
-    real(8)            :: n  ! degrees of freedom as a real(8)
-    real(8)            :: k  ! n - 2
-    real(8)            :: z  ! percentile of normal distribution
-    real(8)            :: z2 ! z * z
-
-    if (df == 1) then
-      ! For one degree of freedom, the t-distribution becomes a Cauchy
-      ! distribution whose cdf we can invert directly
-
-      t = tan(PI*(p - HALF))
-
-    elseif (df == 2) then
-      ! For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
-      ! 2)). This can be directly inverted to yield the solution below
-
-      t = TWO*sqrt(TWO)*(p - HALF)/sqrt(ONE - FOUR*(p - HALF)**2)
-
-    else
-
-      ! This approximation is from E. Olusegun George and Meenakshi Sivaram, "A
-      ! modification of the Fisher-Cornish approximation for the student t
-      ! percentiles," Communication in Statistics - Simulation and Computation,
-      ! 16 (4), pp. 1123-1132 (1987).
-
-      n = real(df,8)
-      k = ONE/(n - TWO)
-      z = normal_percentile(p)
-      z2 = z * z
-      t = sqrt(n*k) * (z + (z2 - THREE)*z*k/FOUR + ((5._8*z2 - 56._8)*z2 + &
-           75._8)*z*k*k/96._8 + (((z2 - 27._8)*THREE*z2 + 417._8)*z2 - 315._8) &
-           *z*k*k*k/384._8)
-
-    end if
-
-  end function t_percentile
-
-!===============================================================================
-! CALC_PN calculates the n-th order Legendre polynomial at the value of x.
-! Since this function is called repeatedly during the neutron transport process,
-! neither n or x is checked to see if they are in the applicable range.
-! This is left to the client developer to use where applicable. x is to be in
-! the domain of [-1,1], and 0<=n<=5. If x is outside of the range, the return
-! value will be outside the expected range; if n is outside the stated range,
-! the return value will be 1.0.
-!===============================================================================
-
-  elemental function calc_pn(n,x) result(pnx)
-
-    integer, intent(in) :: n   ! Legendre order requested
-    real(8), intent(in) :: x   ! Independent variable the Legendre is to be
-                               ! evaluated at; x must be in the domain [-1,1]
-    real(8)             :: pnx ! The Legendre poly of order n evaluated at x
-
-    select case(n)
-    case(1)
-      pnx = x
-    case(2)
-      pnx = 1.5_8 * x * x - HALF
-    case(3)
-      pnx = 2.5_8 * x * x * x - 1.5_8 * x
-    case(4)
-      pnx = 4.375_8 * (x ** 4) - 3.75_8 * x * x + 0.375_8
-    case(5)
-      pnx = 7.875_8 * (x ** 5) - 8.75_8 * x * x * x + 1.875 * x
-    case(6)
-      pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + &
-           6.5625_8 * x * x - 0.3125_8
-    case(7)
-      pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + &
-           19.6875_8 * x * x * x - 2.1875_8 * x
-    case(8)
-      pnx = 50.2734375_8 * (x ** 8) - 93.84375_8 * (x ** 6) + &
-           54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8
-    case(9)
-      pnx = 94.9609375_8 * (x ** 9) - 201.09375_8 * (x ** 7) + &
-           140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x
-    case(10)
-      pnx = 180.42578125_8 * (x ** 10) - 427.32421875_8 * (x ** 8) + &
-           351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + &
-           13.53515625_8 * x * x - 0.24609375_8
-    case default
-      pnx = ONE ! correct for case(0), incorrect for the rest
-    end select
-
-  end function calc_pn
-
-!===============================================================================
-! CALC_RN calculates the n-th order real spherical harmonics for a given angle
-! (in terms of (u,v,w)).  All Rn,m values are provided (where -n<=m<=n)
-!===============================================================================
-
-  pure function calc_rn(n,uvw) result(rn)
-
-    integer, intent(in) :: n      ! Order requested
-    real(8), intent(in) :: uvw(3) ! Direction of travel, assumed to be on unit sphere
-    real(8)             :: rn(2*n + 1)     ! The resultant R_n(uvw)
-
-    real(8) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw)
-    real(8) :: w2m1   ! (w^2 - 1), frequently used in these
-
-    w = uvw(3) ! z = cos(polar)
-    if (uvw(1) == ZERO) then
-      phi = ZERO
-    else
-      phi = atan2(uvw(2), uvw(1))
-    end if
-
-    w2m1 = (ONE - w**2)
-    select case(n)
-    case (0)
-      ! l = 0, m = 0
-      rn(1) = ONE
-    case (1)
-      ! l = 1, m = -1
-      rn(1) = -(ONE*sqrt(w2m1) * sin(phi))
-      ! l = 1, m = 0
-      rn(2) = ONE * w
-      ! l = 1, m = 1
-      rn(3) = -(ONE*sqrt(w2m1) * cos(phi))
-    case (2)
-      ! l = 2, m = -2
-      rn(1) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * sin(TWO*phi)
-      ! l = 2, m = -1
-      rn(2) = -(1.73205080756888_8 * w*sqrt(w2m1) * sin(phi))
-      ! l = 2, m = 0
-      rn(3) = 1.5_8 * w**2 - HALF
-      ! l = 2, m = 1
-      rn(4) = -(1.73205080756888_8 * w*sqrt(w2m1) * cos(phi))
-      ! l = 2, m = 2
-      rn(5) = 0.288675134594813_8 * (-THREE * w**2 + THREE) * cos(TWO*phi)
-    case (3)
-      ! l = 3, m = -3
-      rn(1) = -(0.790569415042095_8 * (w2m1)**(THREE/TWO) * sin(THREE * phi))
-      ! l = 3, m = -2
-      rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi)
-      ! l = 3, m = -1
-      rn(3) = -(0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
-           sin(phi))
-      ! l = 3, m = 0
-      rn(4) = 2.5_8 * w**3 - 1.5_8 * w
-      ! l = 3, m = 1
-      rn(5) = -(0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
-           cos(phi))
-      ! l = 3, m = 2
-      rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi)
-      ! l = 3, m = 3
-      rn(7) = -(0.790569415042095_8 * (w2m1)**(THREE/TWO) * cos(THREE* phi))
-    case (4)
-      ! l = 4, m = -4
-      rn(1) = 0.739509972887452_8 * (w2m1)**2 * sin(4.0_8*phi)
-      ! l = 4, m = -3
-      rn(2) = -(2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi))
-      ! l = 4, m = -2
-      rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
-           sin(TWO*phi)
-      ! l = 4, m = -1
-      rn(4) = -(0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
-           * sin(phi))
-      ! l = 4, m = 0
-      rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8
-      ! l = 4, m = 1
-      rn(6) = -(0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
-           * cos(phi))
-      ! l = 4, m = 2
-      rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
-           cos(TWO*phi)
-      ! l = 4, m = 3
-      rn(8) = -(2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi))
-      ! l = 4, m = 4
-      rn(9) = 0.739509972887452_8 * (w2m1)**2 * cos(4.0_8*phi)
-    case (5)
-      ! l = 5, m = -5
-      rn(1) = -(0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * sin(5.0_8*phi))
-      ! l = 5, m = -4
-      rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi)
-      ! l = 5, m = -3
-      rn(3) = -(0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
-           ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi))
-      ! l = 5, m = -2
-      rn(4) = 0.0487950036474267_8 * (w2m1) &
-           * ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * sin(TWO*phi)
-      ! l = 5, m = -1
-      rn(5) = -(0.258198889747161_8*sqrt(w2m1)* &
-           ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
-           * sin(phi))
-      ! l = 5, m = 0
-      rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w
-      ! l = 5, m = 1
-      rn(7) = -(0.258198889747161_8*sqrt(w2m1)* &
-           ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
-           * cos(phi))
-      ! l = 5, m = 2
-      rn(8) = 0.0487950036474267_8 * (w2m1)* &
-           ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
-      ! l = 5, m = 3
-      rn(9) = -(0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
-           ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi))
-      ! l = 5, m = 4
-      rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi)
-      ! l = 5, m = 5
-      rn(11) = -(0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * cos(5.0_8* phi))
-    case (6)
-      ! l = 6, m = -6
-      rn(1) = 0.671693289381396_8 * (w2m1)**3 * sin(6.0_8*phi)
-      ! l = 6, m = -5
-      rn(2) = -(2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi))
-      ! l = 6, m = -4
-      rn(3) = 0.00104990131391452_8 * (w2m1)**2 * &
-           ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
-      ! l = 6, m = -3
-      rn(4) = -(0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
-           ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi))
-      ! l = 6, m = -2
-      rn(5) = 0.0345032779671177_8 * (w2m1) * &
-           ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
-           * sin(TWO*phi)
-      ! l = 6, m = -1
-      rn(6) = -(0.218217890235992_8*sqrt(w2m1) * &
-           ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
-           * sin(phi))
-      ! l = 6, m = 0
-      rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8
-      ! l = 6, m = 1
-      rn(8) = -(0.218217890235992_8*sqrt(w2m1) * &
-           ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
-           * cos(phi))
-      ! l = 6, m = 2
-      rn(9) = 0.0345032779671177_8 * (w2m1) * &
-           ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
-           * cos(TWO*phi)
-      ! l = 6, m = 3
-      rn(10) = -(0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
-           ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi))
-      ! l = 6, m = 4
-      rn(11) = 0.00104990131391452_8 * (w2m1)**2 * &
-           ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
-      ! l = 6, m = 5
-      rn(12) = -(2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi))
-      ! l = 6, m = 6
-      rn(13) = 0.671693289381396_8 * (w2m1)**3 * cos(6.0_8*phi)
-    case (7)
-      ! l = 7, m = -7
-      rn(1) = -(0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * sin(7.0_8*phi))
-      ! l = 7, m = -6
-      rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi)
-      ! l = 7, m = -5
-      rn(3) = -(9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi))
-      ! l = 7, m = -4
-      rn(4) = 0.000548293079133141_8 * (w2m1)**2* &
-           ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
-      ! l = 7, m = -3
-      rn(5) = -(0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
-                 ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
-                 sin(THREE*phi))
-      ! l = 7, m = -2
-      rn(6) = 0.025717224993682_8 * (w2m1)* &
-           ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
-           sin(TWO*phi)
-      ! l = 7, m = -1
-      rn(7) = -(0.188982236504614_8*sqrt(w2m1)* &
-                 ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
-                 (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi))
-      ! l = 7, m = 0
-      rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 &
-           * w
-      ! l = 7, m = 1
-      rn(9) = -(0.188982236504614_8*sqrt(w2m1)* &
-                 ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
-                  (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi))
-      ! l = 7, m = 2
-      rn(10) = 0.025717224993682_8 * (w2m1)* &
-           ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
-           cos(TWO*phi)
-      ! l = 7, m = 3
-      rn(11) = -(0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
-                 ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
-                 cos(THREE*phi))
-      ! l = 7, m = 4
-      rn(12) = 0.000548293079133141_8 * (w2m1)**2 * &
-           ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
-      ! l = 7, m = 5
-      rn(13) = -(9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi))
-      ! l = 7, m = 6
-      rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi)
-      ! l = 7, m = 7
-      rn(15) = -(0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * cos(7.0_8*phi))
-    case (8)
-      ! l = 8, m = -8
-      rn(1) = 0.626706654240044_8 * (w2m1)**4 * sin(8.0_8*phi)
-      ! l = 8, m = -7
-      rn(2) = -(2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi))
-      ! l = 8, m = -6
-      rn(3) = 6.77369783729086d-6*(w2m1)**3* &
-           ((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi)
-      ! l = 8, m = -5
-      rn(4) = -(4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi))
-      ! l = 8, m = -4
-      rn(5) = 0.000316557156832328_8 * (w2m1)**2* &
-           ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 &
-           + 10395.0_8/8.0_8) * sin(4.0_8*phi)
-      ! l = 8, m = -3
-      rn(6) = -(0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
-                 ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 &
-                 + (10395.0_8/8.0_8)*w) * sin(THREE*phi))
-      ! l = 8, m = -2
-      rn(7) = 0.0199204768222399_8 * (w2m1)* &
-           ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
-           (10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi)
-      ! l = 8, m = -1
-      rn(8) = -(0.166666666666667_8*sqrt(w2m1)* &
-                 ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
-                 (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi))
-      ! l = 8, m = 0
-      rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -&
-           9.84375_8 * w**2 + 0.2734375_8
-      ! l = 8, m = 1
-      rn(10) = -(0.166666666666667_8*sqrt(w2m1)* &
-                 ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
-                 (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi))
-      ! l = 8, m = 2
-      rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- &
-           45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
-           315.0_8/16.0_8) * cos(TWO*phi)
-      ! l = 8, m = 3
-      rn(12) = -(0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
-                 ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
-                 (10395.0_8/8.0_8)*w) * cos(THREE*phi))
-      ! l = 8, m = 4
-      rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - &
-           135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
-      ! l = 8, m = 5
-      rn(14) = -(4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 -&
-                 135135.0_8/TWO*w) * cos(5.0_8*phi))
-      ! l = 8, m = 6
-      rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - &
-           135135.0_8/TWO) * cos(6.0_8*phi)
-      ! l = 8, m = 7
-      rn(16) = -(2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi))
-      ! l = 8, m = 8
-      rn(17) = 0.626706654240044_8 * (w2m1)**4 * cos(8.0_8*phi)
-    case (9)
-      ! l = 9, m = -9
-      rn(1) = -(0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * sin(9.0_8*phi))
-      ! l = 9, m = -8
-      rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi)
-      ! l = 9, m = -7
-      rn(3) = -(4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
-                 ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi))
-      ! l = 9, m = -6
-      rn(4) = 3.02928976464514d-6*(w2m1)**3* &
-           ((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi)
-      ! l = 9, m = -5
-      rn(5) = -(2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + &
-                 135135.0_8/8.0_8) * sin(5.0_8*phi))
-      ! l = 9, m = -4
-      rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
-           675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
-      ! l = 9, m = -3
-      rn(7) = -(0.00173385495536766_8 * (w2m1)**(THREE/TWO)* &
-                 ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
-                 (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi))
-      ! l = 9, m = -2
-      rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- &
-           135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
-           - 3465.0_8/16.0_8 * w) * sin(TWO*phi)
-      ! l = 9, m = -1
-      rn(9) = -(0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
-                 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 &
-                 * w**2 + 315.0_8/128.0_8) * sin(phi))
-      ! l = 9, m = 0
-      rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- &
-           36.09375_8 * w**3 + 2.4609375_8 * w
-      ! l = 9, m = 1
-      rn(11) = -(0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
-                 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 &
-                 * w**2 + 315.0_8/128.0_8) * cos(phi))
-      ! l = 9, m = 2
-      rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - &
-           135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
-           - 3465.0_8/ 16.0_8 * w) * cos(TWO*phi)
-      ! l = 9, m = 3
-      rn(13) = -(0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)&
-                 *w**6 - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 &
-                 - 3465.0_8/16.0_8)* cos(THREE*phi))
-      ! l = 9, m = 4
-      rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
-           675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
-      ! l = 9, m = 5
-      rn(15) = -(2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* &
-                 w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi))
-      ! l = 9, m = 6
-      rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - &
-           2027025.0_8/TWO*w) * cos(6.0_8*phi)
-      ! l = 9, m = 7
-      rn(17) = -(4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
-                 ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi))
-      ! l = 9, m = 8
-      rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi)
-      ! l = 9, m = 9
-      rn(19) = -(0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * cos(9.0_8*phi))
-    case (10)
-      ! l = 10, m = -10
-      rn(1) = 0.593627917136573_8 * (w2m1)**5 * sin(10.0_8*phi)
-      ! l = 10, m = -9
-      rn(2) = -(2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi))
-      ! l = 10, m = -8
-      rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - &
-           34459425.0_8/TWO) * sin(8.0_8*phi)
-      ! l = 10, m = -7
-      rn(4) = -(1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
-                 ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi))
-      ! l = 10, m = -6
-      rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
-           34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi)
-      ! l = 10, m = -5
-      rn(6) = -(1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
-                 (2027025.0_8/8.0_8)*w) * sin(5.0_8*phi))
-      ! l = 10, m = -4
-      rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
-           11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
-           45045.0_8/16.0_8) * sin(4.0_8*phi)
-      ! l = 10, m = -3
-      rn(8) = -(0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
-                 ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
-                 (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi))
-      ! l = 10, m = -2
-      rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
-           765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
-           45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi)
-      ! l = 10, m = -1
-      rn(10) = -(0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
-                 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - &
-                 15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi))
-      ! l = 10, m = 0
-      rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 &
-           * w**6 - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8
-      ! l = 10, m = 1
-      rn(12) = -(0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
-                 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ &
-                 32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi))
-      ! l = 10, m = 2
-      rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
-           765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
-           45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
-      ! l = 10, m = 3
-      rn(14) = -(0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
-                 ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
-                 (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi))
-      ! l = 10, m = 4
-      rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 -&
-           11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
-           45045.0_8/16.0_8) * cos(4.0_8*phi)
-      ! l = 10, m = 5
-      rn(16) = -(1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
-                 ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
-                 (2027025.0_8/8.0_8)*w) * cos(5.0_8*phi))
-      ! l = 10, m = 6
-      rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
-           34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi)
-      ! l = 10, m = 7
-      rn(18) = -(1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
-                 ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi))
-      ! l = 10, m = 8
-      rn(19) = 2.49953651452314d-8*(w2m1)**4* &
-           ((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi)
-      ! l = 10, m = 9
-      rn(20) = -(2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi))
-      ! l = 10, m = 10
-      rn(21) = 0.593627917136573_8 * (w2m1)**5 * cos(10.0_8*phi)
-    case default
-      rn = ONE
-    end select
-
-  end function calc_rn
-
-!===============================================================================
-! CALC_ZN calculates the n-th order modified Zernike polynomial moment for a
-! given angle (rho, theta) location in the unit disk. The normlization of the
-! polynomials is such that the integral of Z_pq*Z_pq over the unit disk is
-! exactly pi
-!===============================================================================
-
-  subroutine calc_zn(n, rho, phi, zn)
-    ! This procedure uses the modified Kintner's method for calculating Zernike
-    ! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
-    ! R. (2003). A comparative analysis of algorithms for fast computation of
-    ! Zernike moments. Pattern Recognition, 36(3), 731-742.
-
-    integer, intent(in) :: n           ! Maximum order
-    real(8), intent(in) :: rho         ! Radial location in the unit disk
-    real(8), intent(in) :: phi         ! Theta (radians) location in the unit disk
-    real(8), intent(out) :: zn(:)      ! The resulting list of coefficients
-
-    real(8) :: sin_phi, cos_phi        ! Sine and Cosine of phi
-    real(8) :: sin_phi_vec(n+1)        ! Contains sin(n*phi)
-    real(8) :: cos_phi_vec(n+1)        ! Contains cos(n*phi)
-    real(8) :: zn_mat(n+1, n+1)        ! Matrix form of the coefficients which is
-                                       ! easier to work with
-    real(8) :: k1, k2, k3, k4          ! Variables for R_m_n calculation
-    integer :: i,p,q                   ! Loop counters
-
-    ! n == radial degree
-    ! m == azimuthal frequency
-
-    ! ==========================================================================
-    ! Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
-    ! following recurrence relations so that only a single sin/cos have to be
-    ! evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html)
-    !
-    ! sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
-    ! cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
-
-    sin_phi = sin(phi)
-    cos_phi = cos(phi)
-
-    sin_phi_vec(1) = 1.0_8
-    cos_phi_vec(1) = 1.0_8
-
-    sin_phi_vec(2) = 2.0_8 * cos_phi
-    cos_phi_vec(2) = cos_phi
-
-    do i = 3, n+1
-      sin_phi_vec(i) = 2.0_8 * cos_phi * sin_phi_vec(i-1) - sin_phi_vec(i-2)
-      cos_phi_vec(i) = 2.0_8 * cos_phi * cos_phi_vec(i-1) - cos_phi_vec(i-2)
-    end do
-
-    do i = 1, n+1
-      sin_phi_vec(i) = sin_phi_vec(i) * sin_phi
-    end do
-
-    ! ==========================================================================
-    ! Calculate R_pq(rho)
-
-    ! Fill the main diagonal first (Eq. 3.9 in Chong)
-    do p = 0, n
-      zn_mat(p+1, p+1) = rho**p
-    end do
-
-    ! Fill in the second diagonal (Eq. 3.10 in Chong)
-    do q = 0, n-2
-      zn_mat(q+2+1, q+1) = (q+2) * zn_mat(q+2+1, q+2+1) - (q+1) * zn_mat(q+1, q+1)
-    end do
-
-    ! Fill in the rest of the values using the original results (Eq. 3.8 in Chong)
-    do p = 4, n
-      k2 = 2 * p * (p - 1) * (p - 2)
-      do q = p-4, 0, -2
-        k1 = (p + q) * (p - q) * (p - 2) / 2
-        k3 = -q**2*(p - 1) - p * (p - 1) * (p - 2)
-        k4 = -p * (p + q - 2) * (p - q - 2) / 2
-        zn_mat(p+1, q+1) = ((k2 * rho**2 + k3) * zn_mat(p-2+1, q+1) + k4 * zn_mat(p-4+1, q+1)) / k1
-      end do
-    end do
-
-    ! Roll into a single vector for easier computation later
-    ! The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
-    ! (2, 2), ....   in (n,m) indices
-    ! Note that the cos and sin vectors are offset by one
-    ! sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
-    ! cos_phi_vec = [1.0, cos(x), cos(2x)... ]
-    i = 1
-    do p = 0, n
-      do q = -p, p, 2
-        if (q < 0) then
-          zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q))
-        else if (q == 0) then
-          zn(i) = zn_mat(p+1, q+1)
-        else
-          zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1)
-        end if
-        i = i + 1
-      end do
-    end do
-  end subroutine calc_zn
-
-!===============================================================================
-! EXPAND_HARMONIC expands a given series of real spherical harmonics
-!===============================================================================
-
-  pure function expand_harmonic(data, order, uvw) result(val)
-    real(8), intent(in) :: data(:)
-    integer, intent(in) :: order
-    real(8), intent(in) :: uvw(3)
-    real(8)             :: val
-
-    integer :: l, lm_lo, lm_hi
-
-    val = data(1)
-    lm_lo = 2
-    lm_hi = 4
-    do l = 1, order - 1
-      val = val + sqrt(TWO * real(l,8) + ONE) * &
-           dot_product(calc_rn(l,uvw), data(lm_lo:lm_hi))
-      lm_lo = lm_hi + 1
-      lm_hi = lm_lo + 2 * (l + 1)
-    end do
-
-  end function expand_harmonic
-
 !===============================================================================
 ! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
 ! and the value of x
 !===============================================================================

-  pure function evaluate_legendre(data, x) result(val)
-    real(8), intent(in) :: data(:)
-    real(8), intent(in) :: x
-    real(8)             :: val
+  pure function evaluate_legendre(data, x) result(val) bind(C)
+    real(C_DOUBLE), intent(in) :: data(:)
+    real(C_DOUBLE), intent(in) :: x
+    real(C_DOUBLE)             :: val

-    integer :: l
-
-    val =  HALF * data(1)
-    do l = 1, size(data) - 1
-      val = val + (real(l,8) +  HALF) * data(l + 1) * calc_pn(l,x)
-    end do
+    val = evaluate_legendre_c_intfc(size(data) - 1, data, x)

  end function evaluate_legendre

@ -724,108 +137,24 @@ contains
 !===============================================================================

  function rotate_angle(uvw0, mu, phi) result(uvw)
-    real(8), intent(in) :: uvw0(3) ! directional cosine
-    real(8), intent(in) :: mu      ! cosine of angle in lab or CM
-    real(8), optional   :: phi     ! azimuthal angle
-    real(8)             :: uvw(3)  ! rotated directional cosine
+    real(C_DOUBLE), intent(in) :: uvw0(3)       ! directional cosine
+    real(C_DOUBLE), intent(in) :: mu            ! cosine of angle in lab or CM
+    real(C_DOUBLE), intent(in), optional :: phi ! azimuthal angle

-    real(8) :: phi_   ! azimuthal angle
-    real(8) :: sinphi ! sine of azimuthal angle
-    real(8) :: cosphi ! cosine of azimuthal angle
-    real(8) :: a      ! sqrt(1 - mu^2)
-    real(8) :: b      ! sqrt(1 - w^2)
-    real(8) :: u0     ! original cosine in x direction
-    real(8) :: v0     ! original cosine in y direction
-    real(8) :: w0     ! original cosine in z direction
+    real(C_DOUBLE) :: uvw(3)  ! rotated directional cosine

-    ! Copy original directional cosines
-    u0 = uvw0(1)
-    v0 = uvw0(2)
-    w0 = uvw0(3)
-
-    ! Sample azimuthal angle in [0,2pi) if none provided
-    if (present(phi)) then
-      phi_ = phi
-    else
-      phi_ = TWO * PI * prn()
-    end if
-
-    ! Precompute factors to save flops
-    sinphi = sin(phi_)
-    cosphi = cos(phi_)
-    a = sqrt(max(ZERO, ONE - mu*mu))
-    b = sqrt(max(ZERO, ONE - w0*w0))
-
-    ! Need to treat special case where sqrt(1 - w**2) is close to zero by
-    ! expanding about the v component rather than the w component
-    if (b > 1e-10) then
-      uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b
-      uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b
-      uvw(3) = mu*w0 - a*b*cosphi
-    else
-      b = sqrt(ONE - v0*v0)
-      uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b
-      uvw(2) = mu*v0 - a*b*cosphi
-      uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b
-    end if
+    uvw = uvw0
+    call rotate_angle_c_intfc(uvw, mu, phi)

  end function rotate_angle

-!===============================================================================
-! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based
-! on a direct sampling scheme. The probability distribution function for a
-! Maxwellian is given as p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can
-! be sampled using rule C64 in the Monte Carlo Sampler LA-9721-MS.
-!===============================================================================
-
-  function maxwell_spectrum(T) result(E_out)
-
-    real(8), intent(in)  :: T     ! tabulated function of incoming E
-    real(8)              :: E_out ! sampled energy
-
-    real(8) :: r1, r2, r3  ! random numbers
-    real(8) :: c           ! cosine of pi/2*r3
-
-    r1 = prn()
-    r2 = prn()
-    r3 = prn()
-
-    ! determine cosine of pi/2*r
-    c = cos(PI/TWO*r3)
-
-    ! determine outgoing energy
-    E_out = -T*(log(r1) + log(r2)*c*c)
-
-  end function maxwell_spectrum
-
-!===============================================================================
-! WATT_SPECTRUM samples the outgoing energy from a Watt energy-dependent fission
-! spectrum. Although fitted parameters exist for many nuclides, generally the
-! continuous tabular distributions (LAW 4) should be used in lieu of the Watt
-! spectrum. This direct sampling scheme is an unpublished scheme based on the
-! original Watt spectrum derivation (See F. Brown's MC lectures).
-!===============================================================================
-
-  function watt_spectrum(a, b) result(E_out)
-
-    real(8), intent(in) :: a     ! Watt parameter a
-    real(8), intent(in) :: b     ! Watt parameter b
-    real(8)             :: E_out ! energy of emitted neutron
-
-    real(8) :: w ! sampled from Maxwellian
-
-    w     = maxwell_spectrum(a)
-    E_out = w + a*a*b/4. + (TWO*prn() - ONE)*sqrt(a*a*b*w)
-
-  end function watt_spectrum
-
 !===============================================================================
 ! FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
 !===============================================================================

-  function faddeeva(z) result(wv)
-    complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
-    complex(8)     :: wv     ! The resulting w(z) value
+  function faddeeva(z) result(wv) bind(C)
+    complex(C_DOUBLE_COMPLEX), intent(in) :: z  ! The point to evaluate Z at
+    complex(C_DOUBLE_COMPLEX)             :: wv ! The resulting w(z) value
    real(C_DOUBLE) :: relerr ! Target relative error in inner loop of MIT
                             !  Faddeeva

@ -853,10 +182,10 @@ contains

  end function faddeeva

-  recursive function w_derivative(z, order) result(wv)
+  recursive function w_derivative(z, order) result(wv) bind(C)
    complex(C_DOUBLE_COMPLEX), intent(in) :: z ! The point to evaluate Z at
-    integer,                   intent(in) :: order
-    complex(8)     :: wv     ! The resulting w(z) value
+    integer(C_INT),            intent(in) :: order
+    complex(C_DOUBLE_COMPLEX)     :: wv     ! The resulting w(z) value

    select case(order)
    case (0)
@ -869,62 +198,4 @@ contains
    end select
  end function w_derivative

-!===============================================================================
-! BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit.  The
-! curvefit is a polynomial of the form
-! a/E + b/sqrt(E) + c + d sqrt(E) ...
-!===============================================================================
-
-  subroutine broaden_wmp_polynomials(E, dopp, n, factors)
-    real(8), intent(in) :: E          ! Energy to evaluate at
-    real(8), intent(in) :: dopp       ! sqrt(atomic weight ratio / kT),
-                                      !  kT given in eV.
-    integer, intent(in) :: n          ! number of components to polynomial
-    real(8), intent(out):: factors(n) ! output leading coefficient
-
-    integer :: i
-
-    real(8) :: sqrtE               ! sqrt(energy)
-    real(8) :: beta                ! sqrt(atomic weight ratio * E / kT)
-    real(8) :: half_inv_dopp2      ! 0.5 / dopp**2
-    real(8) :: quarter_inv_dopp4   ! 0.25 / dopp**4
-    real(8) :: erf_beta            ! error function of beta
-    real(8) :: exp_m_beta2         ! exp(-beta**2)
-
-    sqrtE = sqrt(E)
-    beta = sqrtE * dopp
-    half_inv_dopp2 = HALF / dopp**2
-    quarter_inv_dopp4 = half_inv_dopp2**2
-
-    if (beta > 6.0_8) then
-      ! Save time, ERF(6) is 1 to machine precision.
-      ! beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
-      erf_beta = ONE
-      exp_m_beta2 = ZERO
-    else
-      erf_beta = erf(beta)
-      exp_m_beta2 = exp(-beta**2)
-    end if
-
-    ! Assume that, for sure, we'll use a second order (1/E, 1/V, const)
-    ! fit, and no less.
-
-    factors(1) = erf_beta / E
-    factors(2) = ONE / sqrtE
-    factors(3) = factors(1) * (half_inv_dopp2 + E) &
-         + exp_m_beta2 / (beta * SQRT_PI)
-
-    ! Perform recursive broadening of high order components
-    do i = 1, n-3
-      if (i /= 1) then
-        factors(i+3) = -factors(i-1) * (i - ONE) * i * quarter_inv_dopp4 &
-             + factors(i+1) * (E + (ONE + TWO * i) * half_inv_dopp2)
-      else
-        ! Although it's mathematically identical, factors(0) will contain
-        ! nothing, and we don't want to have to worry about memory.
-        factors(i+3) = factors(i+1)*(E + (ONE + TWO * i) * half_inv_dopp2)
-      end if
-    end do
-  end subroutine broaden_wmp_polynomials
-
 end module math
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@ -0,0 +1,693 @@
+#include "math_functions.h"
+
+namespace openmc {
+
+//==============================================================================
+// Mathematical methods
+//==============================================================================
+
+double normal_percentile_c(double p) {
+  constexpr double p_low = 0.02425;
+  constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
+                           -2.759285104469687e2, 1.383577518672690e2,
+                           -3.066479806614716e1, 2.506628277459239e0};
+  constexpr double b[5] = {-5.447609879822406e1, 1.615858368580409e2,
+                           -1.556989798598866e2, 6.680131188771972e1,
+                           -1.328068155288572e1};
+  constexpr double c[6] = {-7.784894002430293e-3, -3.223964580411365e-1,
+                           -2.400758277161838, -2.549732539343734,
+                           4.374664141464968, 2.938163982698783};
+  constexpr double d[4] = {7.784695709041462e-3, 3.224671290700398e-1,
+                           2.445134137142996, 3.754408661907416};
+
+  // The rational approximation used here is from an unpublished work at
+  // http://home.online.no/~pjacklam/notes/invnorm/
+
+  double z;
+  double q;
+
+  if (p < p_low) {
+    // Rational approximation for lower region.
+
+    q = std::sqrt(-2.0 * std::log(p));
+    z = (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) /
+          ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
+
+  } else if (p <= 1.0 - p_low) {
+    // Rational approximation for central region
+    q = p - 0.5;
+    double r = q * q;
+    z = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q /
+         (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0);
+
+  } else {
+    // Rational approximation for upper region
+
+    q = std::sqrt(-2.0*std::log(1.0 - p));
+    z = -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) /
+          ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
+  }
+
+  // Refinement based on Newton's method
+
+  z = z - (0.5 * std::erfc(-z / std::sqrt(2.0)) - p) * std::sqrt(2.0 * PI) *
+       std::exp(0.5 * z * z);
+
+  return z;
+
+}
+
+
+double t_percentile_c(double p, int df){
+  double t;
+
+  if (df == 1) {
+    // For one degree of freedom, the t-distribution becomes a Cauchy
+    // distribution whose cdf we can invert directly
+
+    t = std::tan(PI*(p - 0.5));
+  } else if (df == 2) {
+     // For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
+     // 2)). This can be directly inverted to yield the solution below
+
+    t = 2.0 * std::sqrt(2.0)*(p - 0.5) /
+         std::sqrt(1. - 4. * std::pow(p - 0.5, 2.));
+  } else {
+    // This approximation is from E. Olusegun George and Meenakshi Sivaram, "A
+    // modification of the Fisher-Cornish approximation for the student t
+    // percentiles," Communication in Statistics - Simulation and Computation,
+    // 16 (4), pp. 1123-1132 (1987).
+    double n = df;
+    double k = 1. / (n - 2.);
+    double z = normal_percentile_c(p);
+    double z2 = z * z;
+    t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 +
+         75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) *
+         z * k * k * k / 384.);
+  }
+
+  return t;
+}
+
+
+void calc_pn_c(int n, double x, double pnx[]) {
+  pnx[0] = 1.;
+  if (n >= 1) {
+    pnx[1] = x;
+  }
+
+  // Use recursion relation to build the higher orders
+  for (int l = 1; l < n; l ++) {
+    pnx[l + 1] = ((2 * l + 1) * x * pnx[l] - l * pnx[l - 1]) / (l + 1);
+  }
+}
+
+
+double evaluate_legendre_c(int n, const double data[], double x) {
+  double pnx[n + 1];
+  double val = 0.0;
+  calc_pn_c(n, x, pnx);
+  for (int l = 0; l <= n; l++) {
+    val += (l + 0.5) * data[l] * pnx[l];
+  }
+  return val;
+}
+
+
+void calc_rn_c(int n, const double uvw[3], double rn[]){
+  // rn[] is assumed to have already been allocated to the correct size
+
+  // Store the cosine of the polar angle and the azimuthal angle
+  double w = uvw[2];
+  double phi;
+  if (uvw[0] == 0.) {
+    phi = 0.;
+  } else {
+    phi = std::atan2(uvw[1], uvw[0]);
+  }
+
+  // Store the shorthand of 1-w * w
+  double w2m1 = 1. - w * w;
+
+  // Now evaluate the spherical harmonics function
+  rn[0] = 1.;
+  int i = 0;
+  for (int l = 1; l <= n; l++) {
+    // Set the index to the start of this order
+    i += 2 * (l - 1) + 1;
+
+    // Now evaluate each
+    switch (l) {
+      case 1:
+        // l = 1, m = -1
+        rn[i] = -(std::sqrt(w2m1) * std::sin(phi));
+        // l = 1, m = 0
+        rn[i + 1] = w;
+        // l = 1, m = 1
+        rn[i + 2] = -(std::sqrt(w2m1) * std::cos(phi));
+        break;
+      case 2:
+        // l = 2, m = -2
+        rn[i] = 0.288675134594813 * (-3. * w * w + 3.) * std::sin(2. * phi);
+        // l = 2, m = -1
+        rn[i + 1] = -(1.73205080756888 * w*std::sqrt(w2m1) * std::sin(phi));
+        // l = 2, m = 0
+        rn[i + 2] = 1.5 * w * w - 0.5;
+        // l = 2, m = 1
+        rn[i + 3] = -(1.73205080756888 * w*std::sqrt(w2m1) * std::cos(phi));
+        // l = 2, m = 2
+        rn[i + 4] = 0.288675134594813 * (-3. * w * w + 3.) * std::cos(2. * phi);
+        break;
+      case 3:
+        // l = 3, m = -3
+        rn[i] = -(0.790569415042095 * std::pow(w2m1, 1.5) * std::sin(3. * phi));
+        // l = 3, m = -2
+        rn[i + 1] = 1.93649167310371 * w*(w2m1) * std::sin(2.*phi);
+        // l = 3, m = -1
+        rn[i + 2] = -(0.408248290463863*std::sqrt(w2m1)*((7.5)*w * w - 3./2.) *
+             std::sin(phi));
+        // l = 3, m = 0
+        rn[i + 3] = 2.5 * std::pow(w, 3) - 1.5 * w;
+        // l = 3, m = 1
+        rn[i + 4] = -(0.408248290463863*std::sqrt(w2m1)*((7.5)*w * w - 3./2.) *
+             std::cos(phi));
+        // l = 3, m = 2
+        rn[i + 5] = 1.93649167310371 * w*(w2m1) * std::cos(2.*phi);
+        // l = 3, m = 3
+        rn[i + 6] = -(0.790569415042095 * std::pow(w2m1, 1.5) * std::cos(3.* phi));
+        break;
+      case 4:
+        // l = 4, m = -4
+        rn[i] = 0.739509972887452 * (w2m1 * w2m1) * std::sin(4.0*phi);
+        // l = 4, m = -3
+        rn[i + 1] = -(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::sin(3.* phi));
+        // l = 4, m = -2
+        rn[i + 2] = 0.074535599249993 * (w2m1)*(52.5 * w * w - 7.5) * std::sin(2. *phi);
+        // l = 4, m = -1
+        rn[i + 3] = -(0.316227766016838*std::sqrt(w2m1)*(17.5 * std::pow(w, 3) - 7.5 * w) *
+             std::sin(phi));
+        // l = 4, m = 0
+        rn[i + 4] = 4.375 * std::pow(w, 4) - 3.75 * w * w + 0.375;
+        // l = 4, m = 1
+        rn[i + 5] = -(0.316227766016838*std::sqrt(w2m1)*(17.5 * std::pow(w, 3) - 7.5*w) *
+             std::cos(phi));
+        // l = 4, m = 2
+        rn[i + 6] = 0.074535599249993 * (w2m1)*(52.5*w * w - 7.5) * std::cos(2.*phi);
+        // l = 4, m = 3
+        rn[i + 7] = -(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::cos(3.* phi));
+        // l = 4, m = 4
+        rn[i + 8] = 0.739509972887452 * w2m1 * w2m1 * std::cos(4.0*phi);
+        break;
+      case 5:
+        // l = 5, m = -5
+        rn[i] = -(0.701560760020114 * std::pow(w2m1, 2.5) * std::sin(5.0 * phi));
+        // l = 5, m = -4
+        rn[i + 1] = 2.21852991866236 * w * w2m1 * w2m1 * std::sin(4.0 * phi);
+        // l = 5, m = -3
+        rn[i + 2] = -(0.00996023841111995 * std::pow(w2m1, 1.5) *
+             ((945.0 /2.)* w * w - 52.5) * std::sin(3.*phi));
+        // l = 5, m = -2
+        rn[i + 3] = 0.0487950036474267 * (w2m1)
+             * ((315.0/2.)* std::pow(w, 3) - 52.5 * w) * std::sin(2.*phi);
+        // l = 5, m = -1
+        rn[i + 4] = -(0.258198889747161*std::sqrt(w2m1) *
+             (39.375 * std::pow(w, 4) - 105.0/4.0 * w * w + 15.0/8.0) * std::sin(phi));
+        // l = 5, m = 0
+        rn[i + 5] = 7.875 * std::pow(w, 5) - 8.75 * std::pow(w, 3) + 1.875 * w;
+        // l = 5, m = 1
+        rn[i + 6] = -(0.258198889747161 * std::sqrt(w2m1)*
+             (39.375 * std::pow(w, 4) - 105.0/4.0 * w * w + 15.0/8.0) * std::cos(phi));
+        // l = 5, m = 2
+        rn[i + 7] = 0.0487950036474267 * (w2m1) *
+             ((315.0 / 2.) * std::pow(w, 3) - 52.5*w) * std::cos(2.*phi);
+        // l = 5, m = 3
+        rn[i + 8] = -(0.00996023841111995 * std::pow(w2m1, 1.5) *
+             ((945.0 / 2.) * w * w - 52.5) * std::cos(3.*phi));
+        // l = 5, m = 4
+        rn[i + 9] = 2.21852991866236 * w * w2m1 * w2m1 * std::cos(4.0*phi);
+        // l = 5, m = 5
+        rn[i + 10] = -(0.701560760020114 * std::pow(w2m1, 2.5) * std::cos(5.0* phi));
+        break;
+      case 6:
+        // l = 6, m = -6
+        rn[i] = 0.671693289381396 * std::pow(w2m1, 3) * std::sin(6.0*phi);
+        // l = 6, m = -5
+        rn[i + 1] = -(2.32681380862329 * w*std::pow(w2m1, 2.5) * std::sin(5.0*phi));
+        // l = 6, m = -4
+        rn[i + 2] = 0.00104990131391452 * w2m1 * w2m1 *
+             ((10395.0/2.) * w * w - 945.0/2.) * std::sin(4.0 * phi);
+        // l = 6, m = -3
+        rn[i + 3] = -(0.00575054632785295 * std::pow(w2m1, 1.5) *
+             ((3465.0/2.) * std::pow(w, 3) - 945.0/2.*w) * std::sin(3.*phi));
+        // l = 6, m = -2
+        rn[i + 4] = 0.0345032779671177 * (w2m1) *
+             ((3465.0/8.0)* std::pow(w, 4) - 945.0/4.0 * w * w + 105.0/8.0) *
+             std::sin(2. * phi);
+        // l = 6, m = -1
+        rn[i + 5] = -(0.218217890235992*std::sqrt(w2m1) *
+             ((693.0/8.0)* std::pow(w, 5)- 315.0/4.0 * std::pow(w, 3) + (105.0/8.0)*w) *
+             std::sin(phi));
+        // l = 6, m = 0
+        rn[i + 6] = 14.4375 * std::pow(w, 6) - 19.6875 * std::pow(w, 4) + 6.5625 * w * w -
+             0.3125;
+        // l = 6, m = 1
+        rn[i + 7] = -(0.218217890235992*std::sqrt(w2m1) *
+             ((693.0/8.0)* std::pow(w, 5)- 315.0/4.0 * std::pow(w, 3) + (105.0/8.0)*w) *
+             std::cos(phi));
+        // l = 6, m = 2
+        rn[i + 8] = 0.0345032779671177 * w2m1 *
+             ((3465.0/8.0)* std::pow(w, 4) -945.0/4.0 * w * w + 105.0/8.0) *
+             std::cos(2.*phi);
+        // l = 6, m = 3
+        rn[i + 9] = -(0.00575054632785295 * std::pow(w2m1, 1.5) *
+             ((3465.0/2.) * std::pow(w, 3) - 945.0/2.*w) * std::cos(3.*phi));
+        // l = 6, m = 4
+        rn[i + 10] = 0.00104990131391452 * w2m1 * w2m1 *
+             ((10395.0/2.)*w * w - 945.0/2.) * std::cos(4.0*phi);
+        // l = 6, m = 5
+        rn[i + 11] = -(2.32681380862329 * w * std::pow(w2m1, 2.5) * std::cos(5.0*phi));
+        // l = 6, m = 6
+        rn[i + 12] = 0.671693289381396 * std::pow(w2m1, 3) * std::cos(6.0*phi);
+        break;
+      case 7:
+        // l = 7, m = -7
+        rn[i] = -(0.647259849287749 * std::pow(w2m1, 3.5) * std::sin(7.0*phi));
+        // l = 7, m = -6
+        rn[i + 1] = 2.42182459624969 * w*std::pow(w2m1, 3) * std::sin(6.0*phi);
+        // l = 7, m = -5
+        rn[i + 2] = -(9.13821798555235e-5*std::pow(w2m1, 2.5) *
+             ((135135.0/2.)*w * w - 10395.0/2.) * std::sin(5.0*phi));
+        // l = 7, m = -4
+        rn[i + 3] = 0.000548293079133141 * w2m1 * w2m1 *
+             ((45045.0/2.)*std::pow(w, 3) - 10395.0/2.*w) * std::sin(4.0*phi);
+        // l = 7, m = -3
+        rn[i + 4] = -(0.00363696483726654 * std::pow(w2m1, 1.5) *
+             ((45045.0/8.0)* std::pow(w, 4) - 10395.0/4.0 * w * w + 945.0/8.0) *
+             std::sin(3.*phi));
+        // l = 7, m = -2
+        rn[i + 5] = 0.025717224993682 * (w2m1) *
+             ((9009.0/8.0)* std::pow(w, 5) -3465.0/4.0 * std::pow(w, 3) + (945.0/8.0)*w) *
+             std::sin(2.*phi);
+        // l = 7, m = -1
+        rn[i + 6] = -(0.188982236504614*std::sqrt(w2m1) *
+             ((3003.0/16.0)* std::pow(w, 6) - 3465.0/16.0 * std::pow(w, 4) +
+             (945.0/16.0)*w * w - 35.0/16.0) * std::sin(phi));
+        // l = 7, m = 0
+        rn[i + 7] = 26.8125 * std::pow(w, 7) - 43.3125 * std::pow(w, 5) + 19.6875 * std::pow(w, 3) -
+             2.1875 * w;
+        // l = 7, m = 1
+        rn[i + 8] = -(0.188982236504614*std::sqrt(w2m1) * ((3003.0/16.0) * std::pow(w, 6) -
+             3465.0/16.0 * std::pow(w, 4) + (945.0/16.0)*w * w - 35.0/16.0) * std::cos(phi));
+        // l = 7, m = 2
+        rn[i + 9] = 0.025717224993682 * (w2m1) * ((9009.0/8.0)* std::pow(w, 5) -
+             3465.0/4.0 * std::pow(w, 3) + (945.0/8.0)*w) * std::cos(2.*phi);
+        // l = 7, m = 3
+        rn[i + 10] = -(0.00363696483726654 * std::pow(w2m1, 1.5) *
+             ((45045.0/8.0)* std::pow(w, 4) - 10395.0/4.0 * w * w + 945.0/8.0) *
+             std::cos(3.*phi));
+        // l = 7, m = 4
+        rn[i + 11] = 0.000548293079133141 * w2m1 * w2m1 *
+             ((45045.0/2.)*std::pow(w, 3) - 10395.0/2.*w) * std::cos(4.0*phi);
+        // l = 7, m = 5
+        rn[i + 12] = -(9.13821798555235e-5*std::pow(w2m1, 2.5) *
+             ((135135.0/2.)*w * w - 10395.0/2.) * std::cos(5.0*phi));
+        // l = 7, m = 6
+        rn[i + 13] = 2.42182459624969 * w*std::pow(w2m1, 3) * std::cos(6.0*phi);
+        // l = 7, m = 7
+        rn[i + 14] = -(0.647259849287749 * std::pow(w2m1, 3.5) * std::cos(7.0*phi));
+        break;
+      case 8:
+        // l = 8, m = -8
+        rn[i] = 0.626706654240044 * std::pow(w2m1, 4) * std::sin(8.0*phi);
+        // l = 8, m = -7
+        rn[i + 1] = -(2.50682661696018 * w*std::pow(w2m1, 3.5) * std::sin(7.0*phi));
+        // l = 8, m = -6
+        rn[i + 2] = 6.77369783729086e-6*std::pow(w2m1, 3)*
+             ((2027025.0/2.)*w * w - 135135.0/2.) * std::sin(6.0*phi);
+        // l = 8, m = -5
+        rn[i + 3] = -(4.38985792528482e-5*std::pow(w2m1, 2.5) *
+                  ((675675.0/2.)*std::pow(w, 3) - 135135.0/2.*w) * std::sin(5.0*phi));
+        // l = 8, m = -4
+        rn[i + 4] = 0.000316557156832328 * w2m1 * w2m1 *
+             ((675675.0/8.0)* std::pow(w, 4) - 135135.0/4.0 * w * w + 10395.0/8.0) *
+             std::sin(4.0*phi);
+        // l = 8, m = -3
+        rn[i + 5] = -(0.00245204119306875 * std::pow(w2m1, 1.5) * ((135135.0/8.0) *
+             std::pow(w, 5) - 45045.0/4.0 * std::pow(w, 3) + (10395.0/8.0)*w) * std::sin(3.*phi));
+        // l = 8, m = -2
+        rn[i + 6] = 0.0199204768222399 * (w2m1) *
+             ((45045.0/16.0)* std::pow(w, 6)- 45045.0/16.0 * std::pow(w, 4) +
+             (10395.0/16.0)*w * w - 315.0/16.0) * std::sin(2.*phi);
+        // l = 8, m = -1
+        rn[i + 7] = -(0.166666666666667*std::sqrt(w2m1) *
+             ((6435.0/16.0)* std::pow(w, 7) - 9009.0/16.0 * std::pow(w, 5) +
+             (3465.0/16.0)*std::pow(w, 3) - 315.0/16.0 * w) * std::sin(phi));
+        // l = 8, m = 0
+        rn[i + 8] = 50.2734375 * std::pow(w, 8) - 93.84375 * std::pow(w, 6) + 54.140625 *
+             std::pow(w, 4) - 9.84375 * w * w + 0.2734375;
+        // l = 8, m = 1
+        rn[i + 9] = -(0.166666666666667*std::sqrt(w2m1) *
+                   ((6435.0/16.0)* std::pow(w, 7) - 9009.0/16.0 * std::pow(w, 5) +
+                   (3465.0/16.0)*std::pow(w, 3) - 315.0/16.0 * w) * std::cos(phi));
+        // l = 8, m = 2
+        rn[i + 10] = 0.0199204768222399 * (w2m1)*((45045.0/16.0)* std::pow(w, 6)-
+             45045.0/16.0 * std::pow(w, 4) + (10395.0/16.0)*w * w -
+             315.0/16.0) * std::cos(2.*phi);
+        // l = 8, m = 3
+        rn[i + 11] = -(0.00245204119306875 * std::pow(w2m1, 1.5)*
+                   ((135135.0/8.0) * std::pow(w, 5) - 45045.0/4.0 * std::pow(w, 3) +
+                   (10395.0/8.0)*w) * std::cos(3.*phi));
+        // l = 8, m = 4
+        rn[i + 12] = 0.000316557156832328 * w2m1 * w2m1*((675675.0/8.0)* std::pow(w, 4) -
+             135135.0/4.0 * w * w + 10395.0/8.0) * std::cos(4.0*phi);
+        // l = 8, m = 5
+        rn[i + 13] = -(4.38985792528482e-5*std::pow(w2m1, 2.5)*((675675.0/2.)*std::pow(w, 3) -
+                   135135.0/2.*w) * std::cos(5.0*phi));
+        // l = 8, m = 6
+        rn[i + 14] = 6.77369783729086e-6*std::pow(w2m1, 3)*((2027025.0/2.)*w * w -
+             135135.0/2.) * std::cos(6.0*phi);
+        // l = 8, m = 7
+        rn[i + 15] = -(2.50682661696018 * w*std::pow(w2m1, 3.5) * std::cos(7.0*phi));
+        // l = 8, m = 8
+        rn[i + 16] = 0.626706654240044 * std::pow(w2m1, 4) * std::cos(8.0*phi);
+        break;
+      case 9:
+        // l = 9, m = -9
+        rn[i] = -(0.609049392175524 * std::pow(w2m1, 4.5) * std::sin(9.0 * phi));
+        // l = 9, m = -8
+        rn[i + 1] = 2.58397773170915 * w*std::pow(w2m1, 4) * std::sin(8.0 * phi);
+        // l = 9, m = -7
+        rn[i + 2] = -(4.37240315267812e-7*std::pow(w2m1, 3.5) *
+             ((34459425.0/2.)*w * w - 2027025.0/2.) * std::sin(7.0 * phi));
+        // l = 9, m = -6
+        rn[i + 3] = 3.02928976464514e-6*std::pow(w2m1, 3)*
+             ((11486475.0/2.)*std::pow(w, 3) - 2027025.0/2.*w) * std::sin(6.0 * phi);
+        // l = 9, m = -5
+        rn[i + 4] = -(2.34647776186144e-5*std::pow(w2m1, 2.5) *
+             ((11486475.0/8.0)* std::pow(w, 4) - 2027025.0 / 4.0 * w * w +
+              135135.0/8.0) * std::sin(5.0 * phi));
+        // l = 9, m = -4
+        rn[i + 5] = 0.000196320414650061 * w2m1 * w2m1*((2297295.0/8.0)* std::pow(w, 5) -
+             675675.0/4.0 * std::pow(w, 3) + (135135.0/8.0)*w) * std::sin(4.0*phi);
+        // l = 9, m = -3
+        rn[i + 6] = -(0.00173385495536766 * std::pow(w2m1, 1.5) *
+                  ((765765.0/16.0)* std::pow(w, 6) - 675675.0/16.0 * std::pow(w, 4) +
+                  (135135.0/16.0)*w * w - 3465.0/16.0) * std::sin(3. * phi));
+        // l = 9, m = -2
+        rn[i + 7] = 0.0158910431540932 * (w2m1)*((109395.0/16.0)* std::pow(w, 7)-
+             135135.0/16.0 * std::pow(w, 5) + (45045.0/16.0)*std::pow(w, 3) -
+             3465.0/16.0 * w) * std::sin(2. * phi);
+        // l = 9, m = -1
+        rn[i + 8] = -(0.149071198499986*std::sqrt(w2m1)*((109395.0/128.0)* std::pow(w, 8) -
+                  45045.0/32.0 * std::pow(w, 6) + (45045.0/64.0)* std::pow(w, 4) -
+                  3465.0/32.0 * w * w + 315.0/128.0) * std::sin(phi));
+        // l = 9, m = 0
+        rn[i + 9] = 94.9609375 * std::pow(w, 9) - 201.09375 * std::pow(w, 7) +
+             140.765625 * std::pow(w, 5)- 36.09375 * std::pow(w, 3) + 2.4609375 * w;
+        // l = 9, m = 1
+        rn[i + 10] = -(0.149071198499986*std::sqrt(w2m1)*((109395.0/128.0)* std::pow(w, 8) -
+                   45045.0/32.0 * std::pow(w, 6) + (45045.0/64.0)* std::pow(w, 4) -
+                   3465.0/32.0 * w * w + 315.0/128.0) * std::cos(phi));
+        // l = 9, m = 2
+        rn[i + 11] = 0.0158910431540932 * (w2m1)*((109395.0/16.0)* std::pow(w, 7) -
+             135135.0/16.0 * std::pow(w, 5) + (45045.0/16.0)*std::pow(w, 3) -
+             3465.0/ 16.0 * w) * std::cos(2. * phi);
+        // l = 9, m = 3
+        rn[i + 12] = -(0.00173385495536766 * std::pow(w2m1, 1.5)*((765765.0/16.0) *
+                   std::pow(w, 6) - 675675.0/16.0 * std::pow(w, 4) +
+                   (135135.0/16.0)* w * w - 3465.0/16.0)* std::cos(3. * phi));
+        // l = 9, m = 4
+        rn[i + 13] = 0.000196320414650061 * w2m1 * w2m1*((2297295.0/8.0) * std::pow(w, 5) -
+             675675.0/4.0 * std::pow(w, 3) + (135135.0/8.0)*w) * std::cos(4.0 * phi);
+        // l = 9, m = 5
+        rn[i + 14] = -(2.34647776186144e-5*std::pow(w2m1, 2.5)*((11486475.0/8.0) *
+                   std::pow(w, 4) - 2027025.0/4.0 * w * w + 135135.0/8.0) *
+                   std::cos(5.0 * phi));
+        // l = 9, m = 6
+        rn[i + 15] = 3.02928976464514e-6*std::pow(w2m1, 3)*((11486475.0/2.)*std::pow(w, 3) -
+             2027025.0/2. * w) * std::cos(6.0 * phi);
+        // l = 9, m = 7
+        rn[i + 16] = -(4.37240315267812e-7*std::pow(w2m1, 3.5)*
+                   ((34459425.0/2.) * w * w - 2027025.0/2.) * std::cos(7.0 * phi));
+        // l = 9, m = 8
+        rn[i + 17] = 2.58397773170915 * w*std::pow(w2m1, 4) * std::cos(8.0 * phi);
+        // l = 9, m = 9
+        rn[i + 18] = -(0.609049392175524 * std::pow(w2m1, 4.5) * std::cos(9.0 * phi));
+        break;
+      case 10:
+        // l = 10, m = -10
+        rn[i] = 0.593627917136573 * std::pow(w2m1, 5) * std::sin(10.0 * phi);
+        // l = 10, m = -9
+        rn[i + 1] = -(2.65478475211798 * w * std::pow(w2m1, 4.5) * std::sin(9.0 * phi));
+        // l = 10, m = -8
+        rn[i + 2] = 2.49953651452314e-8 * std::pow(w2m1, 4) *
+             ((654729075.0/2.) * w * w - 34459425.0/2.) * std::sin(8.0 * phi);
+        // l = 10, m = -7
+        rn[i + 3] = -(1.83677671621093e-7*std::pow(w2m1, 3.5)*
+                  ((218243025.0/2.)*std::pow(w, 3) - 34459425.0/2.*w) *
+                  std::sin(7.0 * phi));
+        // l = 10, m = -6
+        rn[i + 4] = 1.51464488232257e-6*std::pow(w2m1, 3)*((218243025.0/8.0)* std::pow(w, 4) -
+             34459425.0/4.0 * w * w + 2027025.0/8.0) * std::sin(6.0 * phi);
+        // l = 10, m = -5
+        rn[i + 5] = -(1.35473956745817e-5*std::pow(w2m1, 2.5)*
+                  ((43648605.0/8.0)* std::pow(w, 5) - 11486475.0/4.0 * std::pow(w, 3) +
+                  (2027025.0/8.0)*w) * std::sin(5.0 * phi));
+        // l = 10, m = -4
+        rn[i + 6] = 0.000128521880085575 * w2m1 * w2m1*((14549535.0/16.0)* std::pow(w, 6) -
+             11486475.0/16.0 * std::pow(w, 4) + (2027025.0/16.0)*w * w -
+             45045.0/16.0) * std::sin(4.0 * phi);
+        // l = 10, m = -3
+        rn[i + 7] = -(0.00127230170115096 * std::pow(w2m1, 1.5)*
+                  ((2078505.0/16.0)* std::pow(w, 7) - 2297295.0/16.0 * std::pow(w, 5) +
+                  (675675.0/16.0)*std::pow(w, 3) - 45045.0/16.0 * w) * std::sin(3. * phi));
+        // l = 10, m = -2
+        rn[i + 8] = 0.012974982402692 * (w2m1)*((2078505.0/128.0)* std::pow(w, 8) -
+             765765.0/32.0 * std::pow(w, 6) + (675675.0/64.0)* std::pow(w, 4) -
+             45045.0/32.0 * w * w + 3465.0/128.0) * std::sin(2. * phi);
+        // l = 10, m = -1
+        rn[i + 9] = -(0.134839972492648*std::sqrt(w2m1)*((230945.0/128.0)* std::pow(w, 9) -
+                   109395.0/32.0 * std::pow(w, 7) + (135135.0/64.0)* std::pow(w, 5) -
+                   15015.0/32.0 * std::pow(w, 3) + (3465.0/128.0)*w) * std::sin(phi));
+        // l = 10, m = 0
+        rn[i + 10] = 180.42578125 * std::pow(w, 10) - 427.32421875 * std::pow(w, 8) +351.9140625
+             * std::pow(w, 6) - 117.3046875 * std::pow(w, 4) + 13.53515625 * w * w -0.24609375;
+        // l = 10, m = 1
+        rn[i + 11] = -(0.134839972492648*std::sqrt(w2m1)*((230945.0/128.0)* std::pow(w, 9) -
+                   109395.0/32.0 * std::pow(w, 7) + (135135.0/64.0)* std::pow(w, 5) -15015.0/
+                   32.0 * std::pow(w, 3) + (3465.0/128.0)*w) * std::cos(phi));
+        // l = 10, m = 2
+        rn[i + 12] = 0.012974982402692 * (w2m1)*((2078505.0/128.0)* std::pow(w, 8) -
+             765765.0/32.0 * std::pow(w, 6) + (675675.0/64.0)* std::pow(w, 4) -
+             45045.0/32.0 * w * w + 3465.0/128.0) * std::cos(2. * phi);
+        // l = 10, m = 3
+        rn[i + 13] = -(0.00127230170115096 * std::pow(w2m1, 1.5)*
+                   ((2078505.0/16.0)* std::pow(w, 7) - 2297295.0/16.0 * std::pow(w, 5) +
+                   (675675.0/16.0)*std::pow(w, 3) - 45045.0/16.0 * w) * std::cos(3. * phi));
+        // l = 10, m = 4
+        rn[i + 14] = 0.000128521880085575 * w2m1 * w2m1*((14549535.0/16.0)* std::pow(w, 6) -
+             11486475.0/16.0 * std::pow(w, 4) + (2027025.0/16.0) * w * w -
+             45045.0/16.0) * std::cos(4.0 * phi);
+        // l = 10, m = 5
+        rn[i + 15] = -(1.35473956745817e-5*std::pow(w2m1, 2.5)*
+                   ((43648605.0/8.0)* std::pow(w, 5) - 11486475.0/4.0 * std::pow(w, 3) +
+                   (2027025.0/8.0)*w) * std::cos(5.0 * phi));
+        // l = 10, m = 6
+        rn[i + 16] = 1.51464488232257e-6*std::pow(w2m1, 3)*((218243025.0/8.0)* std::pow(w, 4) -
+             34459425.0/4.0 * w * w + 2027025.0/8.0) * std::cos(6.0 * phi);
+        // l = 10, m = 7
+        rn[i + 17] = -(1.83677671621093e-7*std::pow(w2m1, 3.5) *
+             ((218243025.0/2.)*std::pow(w, 3) - 34459425.0/2.*w) * std::cos(7.0 * phi));
+        // l = 10, m = 8
+        rn[i + 18] = 2.49953651452314e-8*std::pow(w2m1, 4)*
+             ((654729075.0/2.)*w * w - 34459425.0/2.) * std::cos(8.0 * phi);
+        // l = 10, m = 9
+        rn[i + 19] = -(2.65478475211798 * w*std::pow(w2m1, 4.5) * std::cos(9.0 * phi));
+        // l = 10, m = 10
+        rn[i + 20] = 0.593627917136573 * std::pow(w2m1, 5) * std::cos(10.0 * phi);
+    }
+  }
+}
+
+
+void calc_zn_c(int n, double rho, double phi, double zn[]) {
+  // ===========================================================================
+  // Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
+  // following recurrence relations so that only a single sin/cos have to be
+  // evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html)
+  //
+  // sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
+  // cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
+
+  double sin_phi = std::sin(phi);
+  double cos_phi = std::cos(phi);
+
+  double sin_phi_vec[n + 1]; // Sin[n * phi]
+  double cos_phi_vec[n + 1]; // Cos[n * phi]
+  sin_phi_vec[0] = 1.0;
+  cos_phi_vec[0] = 1.0;
+  sin_phi_vec[1] = 2.0 * cos_phi;
+  cos_phi_vec[1] = cos_phi;
+
+  for (int i = 2; i <= n; i++) {
+    sin_phi_vec[i] = 2. * cos_phi * sin_phi_vec[i - 1] - sin_phi_vec[i - 2];
+    cos_phi_vec[i] = 2. * cos_phi * cos_phi_vec[i - 1] - cos_phi_vec[i - 2];
+  }
+
+  for (int i = 0; i <= n; i++) {
+    sin_phi_vec[i] *=  sin_phi;
+  }
+
+  // ===========================================================================
+  // Calculate R_pq(rho)
+  double zn_mat[n + 1][n + 1]; // Matrix forms of the coefficients which are
+                               // easier to work with
+
+  // Fill the main diagonal first (Eq 3.9 in Chong)
+  for (int p = 0; p <= n; p++) {
+    zn_mat[p][p] = std::pow(rho, p);
+  }
+
+  // Fill the 2nd diagonal (Eq 3.10 in Chong)
+  for (int q = 0; q <= n - 2; q++) {
+    zn_mat[q][q+2] = (q + 2) * zn_mat[q+2][q+2] - (q + 1) * zn_mat[q][q];
+  }
+
+  // Fill in the rest of the values using the original results (Eq. 3.8 in Chong)
+  for (int p = 4; p <= n; p++) {
+    double k2 = 2 * p * (p - 1) * (p - 2);
+    for (int q = p - 4; q >= 0; q -= 2) {
+      double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;
+      double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);
+      double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;
+      zn_mat[q][p] =
+        ((k2 * rho * rho + k3) * zn_mat[q][p-2] + k4 * zn_mat[q][p-4]) / k1;
+    }
+  }
+
+  // Roll into a single vector for easier computation later
+  // The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
+  // (2, 2), ....   in (n,m) indices
+  // Note that the cos and sin vectors are offset by one
+  // sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
+  // cos_phi_vec = [1.0, cos(x), cos(2x)... ]
+  int i = 0;
+  for (int p = 0; p <= n; p++) {
+    for (int q = -p; q <= p; q += 2) {
+      if (q < 0) {
+        zn[i] = zn_mat[std::abs(q)][p] * sin_phi_vec[std::abs(q) - 1];
+      } else if (q == 0) {
+        zn[i] = zn_mat[q][p];
+      } else {
+        zn[i] = zn_mat[q][p] * cos_phi_vec[q];
+      }
+      i++;
+    }
+  }
+
+}
+
+
+void rotate_angle_c(double uvw[3], double mu, double* phi) {
+  // Copy original directional cosines
+  double u0 = uvw[0]; // original cosine in x direction
+  double v0 = uvw[1]; // original cosine in y direction
+  double w0 = uvw[2]; // original cosine in z direction
+
+  // Sample azimuthal angle in [0,2pi) if none provided
+  double phi_;
+  if (phi != nullptr) {
+    phi_ = (*phi);
+  } else {
+    phi_ = 2. * PI * prn();
+  }
+
+  // Precompute factors to save flops
+  double sinphi = std::sin(phi_);
+  double cosphi = std::cos(phi_);
+  double a = std::sqrt(std::fmax(0., 1. - mu * mu));
+  double b = std::sqrt(std::fmax(0., 1. - w0 * w0));
+
+  // Need to treat special case where sqrt(1 - w**2) is close to zero by
+  // expanding about the v component rather than the w component
+  if (b > 1e-10) {
+    uvw[0] = mu * u0 + a * (u0 * w0 * cosphi - v0 * sinphi) / b;
+    uvw[1] = mu * v0 + a * (v0 * w0 * cosphi + u0 * sinphi) / b;
+    uvw[2] = mu * w0 - a * b * cosphi;
+  } else {
+    b = std::sqrt(1. - v0 * v0);
+    uvw[0] = mu * u0 + a * (u0 * v0 * cosphi + w0 * sinphi) / b;
+    uvw[1] = mu * v0 - a * b * cosphi;
+    uvw[2] = mu * w0 + a * (v0 * w0 * cosphi - u0 * sinphi) / b;
+  }
+}
+
+
+double maxwell_spectrum_c(double T) {
+  // Set the random numbers
+  double r1 = prn();
+  double r2 = prn();
+  double r3 = prn();
+
+  // determine cosine of pi/2*r
+  double c = std::cos(PI / 2. * r3);
+
+  // Determine outgoing energy
+  double E_out = -T * (std::log(r1) + std::log(r2) * c * c);
+
+  return E_out;
+}
+
+
+double watt_spectrum_c(double a, double b) {
+  double w = maxwell_spectrum_c(a);
+  double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);
+
+  return E_out;
+}
+
+
+void broaden_wmp_polynomials_c(double E, double dopp, int n, double factors[]) {
+  // Factors is already pre-allocated
+  double sqrtE = std::sqrt(E);
+  double beta = sqrtE * dopp;
+  double half_inv_dopp2 = 0.5 / (dopp * dopp);
+  double quarter_inv_dopp4 = half_inv_dopp2 * half_inv_dopp2;
+
+  double erf_beta;    // error function of beta
+  double exp_m_beta2; // exp(-beta**2)
+  if (beta > 6.0) {
+    // Save time, ERF(6) is 1 to machine precision.
+    // beta/sqrtpi*exp(-beta**2) is also approximately 1 machine epsilon.
+    erf_beta = 1.;
+    exp_m_beta2 = 0.;
+  } else {
+    erf_beta = std::erf(beta);
+    exp_m_beta2 = std::exp(-beta * beta);
+  }
+
+  // Assume that, for sure, we'll use a second order (1/E, 1/V, const)
+  // fit, and no less.
+
+  factors[0] = erf_beta / E;
+  factors[1] = 1. / sqrtE;
+  factors[2] = factors[0] * (half_inv_dopp2 + E) + exp_m_beta2 /
+       (beta * SQRT_PI);
+
+  // Perform recursive broadening of high order components
+  for (int i = 0; i < n - 3; i++) {
+    double ip1_dbl = i + 1;
+    if (i != 0) {
+      factors[i + 3] = -factors[i - 1] * (ip1_dbl - 1.) * ip1_dbl *
+           quarter_inv_dopp4 + factors[i + 1] *
+           (E + (1. + 2. * ip1_dbl) * half_inv_dopp2);
+    } else {
+      // Although it's mathematically identical, factors[0] will contain
+      // nothing, and we don't want to have to worry about memory.
+      factors[i + 3] = factors[i + 1] *
+           (E + (1. + 2. * ip1_dbl) * half_inv_dopp2);
+    }
+  }
+}
+
+} // namespace openmc
--- a/src/math_functions.h
+++ b/src/math_functions.h
@ -0,0 +1,163 @@
+//! \file math_functions.h
+//! A collection of elementary math functions.
+
+#ifndef MATH_FUNCTIONS_H
+#define MATH_FUNCTIONS_H
+
+#include <cmath>
+#include <cstdlib>
+
+#include "random_lcg.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Module constants.
+//==============================================================================
+
+// TODO: cmath::M_PI has 3 more digits precision than the Fortran constant we
+// use so for now we will reuse the Fortran constant until we are OK with
+// modifying test results
+extern "C" constexpr double PI {3.1415926535898};
+
+extern "C" const double SQRT_PI {std::sqrt(PI)};
+
+//==============================================================================
+//! Calculate the percentile of the standard normal distribution with a
+//! specified probability level.
+//!
+//! @param p The probability level
+//! @return The requested percentile
+//==============================================================================
+
+extern "C" double normal_percentile_c(double p);
+
+//==============================================================================
+//! Calculate the percentile of the Student's t distribution with a specified
+//! probability level and number of degrees of freedom.
+//!
+//! @param p  The probability level
+//! @param df The degrees of freedom
+//! @return The requested percentile
+//==============================================================================
+
+extern "C" double t_percentile_c(double p, int df);
+
+//==============================================================================
+//! Calculate the n-th order Legendre polynomials at the value of x.
+//!
+//! @param n   The maximum order requested
+//! @param x   The value to evaluate at; x is expected to be within [-1,1]
+//! @param pnx The requested Legendre polynomials of order 0 to n (inclusive)
+//!   evaluated at x.
+//==============================================================================
+
+extern "C" void calc_pn_c(int n, double x, double pnx[]);
+
+//==============================================================================
+//! Find the value of f(x) given a set of Legendre coefficients and the value
+//! of x.
+//!
+//! @param n    The maximum order of the expansion
+//! @param data The polynomial expansion coefficient data; without the (2l+1)/2
+//!   factor.
+//! @param x    The value to evaluate at; x is expected to be within [-1,1]
+//! @return The requested Legendre polynomials of order 0 to n (inclusive)
+//!   evaluated at x
+//==============================================================================
+
+extern "C" double evaluate_legendre_c(int n, const double data[], double x);
+
+//==============================================================================
+//! Calculate the n-th order real spherical harmonics for a given angle (in
+//! terms of (u,v,w)) for all 0<=n and -m<=n<=n.
+//!
+//! @param n      The maximum order requested
+//! @param uvw[3] The direction the harmonics are requested at
+//! @param rn     The requested harmonics of order 0 to n (inclusive)
+//!   evaluated at uvw.
+//==============================================================================
+
+extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
+
+//==============================================================================
+//! Calculate the n-th order modified Zernike polynomial moment for a given
+//! angle (rho, theta) location on the unit disk.
+//!
+//! This procedure uses the modified Kintner's method for calculating Zernike
+//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
+//! R. (2003). A comparative analysis of algorithms for fast computation of
+//! Zernike moments. Pattern Recognition, 36(3), 731-742.
+//! The normalization of the polynomials is such that the integral of Z_pq^2
+//! over the unit disk is exactly pi.
+//!
+//! @param n   The maximum order requested
+//! @param rho The radial parameter to specify location on the unit disk
+//! @param phi The angle parameter to specify location on the unit disk
+//! @param zn  The requested moments of order 0 to n (inclusive)
+//!   evaluated at rho and phi.
+//==============================================================================
+
+extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
+
+//==============================================================================
+//! Rotate the direction cosines through a polar angle whose cosine is mu and
+//! through an azimuthal angle sampled uniformly.
+//!
+//! This is done with direct sampling rather than rejection sampling as is done
+//! in MCNP and Serpent.
+//!
+//! @param uvw[3] The initial, and final, direction vector
+//! @param mu     The cosine of angle in lab or CM
+//! @param phi    The azimuthal angle; will randomly chosen angle if a nullptr
+//!   is passed
+//==============================================================================
+
+extern "C" void rotate_angle_c(double uvw[3], double mu, double* phi);
+
+//==============================================================================
+//! Samples an energy from the Maxwell fission distribution based on a direct
+//! sampling scheme.
+//!
+//! The probability distribution function for a Maxwellian is given as
+//! p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can be sampled using
+//! rule C64 in the Monte Carlo Sampler LA-9721-MS.
+//!
+//! @param T The tabulated function of the incoming energy
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double maxwell_spectrum_c(double T);
+
+//==============================================================================
+//! Samples an energy from a Watt energy-dependent fission distribution.
+//!
+//! Although fitted parameters exist for many nuclides, generally the
+//! continuous tabular distributions (LAW 4) should be used in lieu of the Watt
+//! spectrum. This direct sampling scheme is an unpublished scheme based on the
+//! original Watt spectrum derivation (See F. Brown's MC lectures).
+//!
+//! @param a Watt parameter a
+//! @param b Watt parameter b
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double watt_spectrum_c(double a, double b);
+
+//==============================================================================
+//! Doppler broadens the windowed multipole curvefit.
+//!
+//! The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E)...
+//!
+//! @param E       The energy to evaluate the broadening at
+//! @param dopp    sqrt(atomic weight ratio / kT) with kT given in eV
+//! @param n       The number of components to the polynomial
+//! @param factors The output leading coefficient
+//==============================================================================
+
+extern "C" void broaden_wmp_polynomials_c(double E, double dopp, int n,
+                                          double factors[]);
+
+} // namespace openmc
+#endif // MATH_FUNCTIONS_H
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@ -7,7 +7,7 @@ module tally
  use dict_header,      only: EMPTY
  use error,            only: fatal_error
  use geometry_header
-  use math,             only: t_percentile, calc_pn, calc_rn
+  use math,             only: t_percentile
  use mesh_header,      only: RegularMesh, meshes
  use message_passing
  use mgxs_header
--- a/src/tallies/tally_filter_cell.F90
+++ b/src/tallies/tally_filter_cell.F90
@ -14,13 +14,14 @@ module tally_filter_cell

  implicit none
  private
+  public :: openmc_cell_filter_get_bins

 !===============================================================================
 ! CELLFILTER specifies which geometric cells tally events reside in.
 !===============================================================================

  type, public, extends(TallyFilter) :: CellFilter
-    integer, allocatable :: cells(:)
+    integer(C_INT32_T), allocatable :: cells(:)
    type(DictIntInt)     :: map
  contains
    procedure :: from_xml
@ -116,4 +117,28 @@ contains
    label = "Cell " // to_str(cells(this % cells(bin)) % id)
  end function text_label_cell

+!===============================================================================
+!                               C API FUNCTIONS
+!===============================================================================
+
+  function openmc_cell_filter_get_bins(index, cells, n) result(err) bind(C)
+    ! Return the cells associated with a cell filter
+    integer(C_INT32_T), value :: index
+    type(C_PTR), intent(out) :: cells
+    integer(C_INT32_T), intent(out) :: n
+    integer(C_INT) :: err
+
+    err = verify_filter(index)
+    if (err == 0) then
+       select type (f => filters(index) % obj)
+       type is (CellFilter)
+          cells = C_LOC(f % cells)
+          n = size(f % cells)
+       class default
+          err = E_INVALID_TYPE
+          call set_errmsg("Tried to get cells from a non-cell filter.")
+       end select
+    end if
+  end function openmc_cell_filter_get_bins
+
 end module tally_filter_cell
--- a/src/tallies/tally_filter_distribcell.F90
+++ b/src/tallies/tally_filter_distribcell.F90
@ -278,9 +278,9 @@ contains
            old_k = 1

            ! Loop over lattice coordinates
-            do k = 1, n_x
+            do m = 1, n_z
              do l = 1, n_y
-                do m = 1, n_z
+                do k = 1, n_x

                  if (target_offset >= lat % offset(map, k, l, m) + offset) then
                    if (k == n_x .and. l == n_y .and. m == n_z) then
--- a/src/tallies/tally_filter_legendre.F90
+++ b/src/tallies/tally_filter_legendre.F90
@ -51,13 +51,12 @@ contains
    type(TallyFilterMatch),   intent(inout) :: match

    integer :: i
-    real(8) :: wgt
+    real(C_DOUBLE) :: wgt(this % n_bins)

-    ! TODO: Use recursive formula to calculate higher orders
-    do i = 0, this % order
-      wgt = calc_pn(i, p % mu)
-      call match % bins % push_back(i + 1)
-      call match % weights % push_back(wgt)
+    call calc_pn(this % order, p % mu, wgt)
+    do i = 1, this % n_bins
+      call match % bins % push_back(i)
+      call match % weights % push_back(wgt(i))
    end do
  end subroutine get_all_bins

--- a/src/tallies/tally_filter_sph_harm.F90
+++ b/src/tallies/tally_filter_sph_harm.F90
@ -74,30 +74,32 @@ contains

    integer :: i, j, n
    integer :: num_nm
-    real(8) :: wgt
-    real(8) :: rn(2*this % order + 1)
+    real(C_DOUBLE) :: wgt(this % order + 1)
+    real(C_DOUBLE) :: rn(this % n_bins)
+
+    ! Determine cosine term for scatter expansion if necessary
+    if (this % cosine == COSINE_SCATTER) then
+      call calc_pn(this % order, p % mu, wgt)
+    else
+      wgt = ONE
+    end if
+
+    ! Find the Rn,m values
+    call calc_rn(this % order, p % last_uvw, rn)

-    ! TODO: Use recursive formula to calculate higher orders
    j = 0
    do n = 0, this % order
-      ! Determine cosine term for scatter expansion if necessary
-      if (this % cosine == COSINE_SCATTER) then
-        wgt = calc_pn(n, p % mu)
-      else
-        wgt = ONE
-      end if
-
      ! Calculate n-th order spherical harmonics for (u,v,w)
      num_nm = 2*n + 1
-      rn(1:num_nm) = calc_rn(n, p % last_uvw)

      ! Append matching (bin,weight) for each moment
      do i = 1, num_nm
        j = j + 1
        call match % bins % push_back(j)
-        call match % weights % push_back(wgt * rn(i))
+        call match % weights % push_back(wgt(n + 1) * rn(j))
      end do
    end do
+
  end subroutine get_all_bins

  subroutine to_statepoint(this, filter_group)
--- a/src/tallies/tally_filter_sptl_legendre.F90
+++ b/src/tallies/tally_filter_sptl_legendre.F90
@ -75,20 +75,19 @@ contains
    type(TallyFilterMatch),   intent(inout) :: match

    integer :: i
-    real(8) :: wgt
-    real(8) :: x       ! Position on specified axis
-    real(8) :: x_norm  ! Normalized position
+    real(C_DOUBLE) :: wgt(this % n_bins)
+    real(C_DOUBLE) :: x       ! Position on specified axis
+    real(C_DOUBLE) :: x_norm  ! Normalized position

    x = p % coord(1) % xyz(this % axis)
    if (this % min <= x .and. x <= this % max) then
      ! Calculate normalized position between min and max
      x_norm = TWO*(x - this % min)/(this % max - this % min) - ONE

-      ! TODO: Use recursive formula to calculate higher orders
-      do i = 0, this % order
-        wgt = calc_pn(i, x_norm)
-        call match % bins % push_back(i + 1)
-        call match % weights % push_back(wgt)
+      call calc_pn(this % order, x_norm, wgt)
+      do i = 1, this % n_bins
+        call match % bins % push_back(i)
+        call match % weights % push_back(wgt(i))
      end do
    end if
  end subroutine get_all_bins
--- a/src/tallies/tally_filter_zernike.F90
+++ b/src/tallies/tally_filter_zernike.F90
@ -61,7 +61,7 @@ contains

    integer :: i
    real(8) :: x, y, r, theta
-    real(8) :: zn(this % n_bins)
+    real(C_DOUBLE) :: zn(this % n_bins)

    ! Determine normalized (r,theta) positions
    x = p % coord(1) % xyz(1) - this % x
--- a/tests/regression_tests/asymmetric_lattice/results_true.dat
+++ b/tests/regression_tests/asymmetric_lattice/results_true.dat
@ -1 +1 @@
-3b76b468b9c0e7df6508189b75748a1b7b4f2b37486396749e1a15e536526336f70b60bb207095c39ecbd46822e8c8705ea81184a3c8546e7da09bb905d74637
+bd3fd10177bc0c7b8542f15228decf8608de013872d201e6ae885cf9b5b6d7516ddb32044ee2f7fd5b207e7d37d5ce3f03e347a2e850953bba4fc26b150c89d2
--- a/tests/regression_tests/distribmat/inputs_true.dat
+++ b/tests/regression_tests/distribmat/inputs_true.dat
@ -1,7 +1,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <geometry>
  <cell id="1" material="1" universe="1" />
-  <cell id="11" material="2 3 void 2" region="-9" universe="11" />
+  <cell id="11" material="2 void 3 2" region="-9" universe="11" />
  <cell id="12" material="1" region="9" universe="11" />
  <cell fill="101" id="101" region="10 -11 12 -13" universe="0" />
  <lattice id="101">
--- a/tests/regression_tests/distribmat/results_true.dat
+++ b/tests/regression_tests/distribmat/results_true.dat
@ -3,7 +3,7 @@ k-combined:
 Cell
 	ID             =	11
 	Name           =	
-	Fill           =	[2, 3, None, 2]
+	Fill           =	[2, None, 3, 2]
 	Region         =	-9
 	Rotation       =	None
 	Translation    =	None
--- a/tests/regression_tests/distribmat/test.py
+++ b/tests/regression_tests/distribmat/test.py
@ -34,7 +34,7 @@ class DistribmatTestHarness(PyAPITestHarness):

        r0 = openmc.ZCylinder(R=0.3)
        c11 = openmc.Cell(cell_id=11, region=-r0)
-        c11.fill = [dense_fuel, light_fuel, None, dense_fuel]
+        c11.fill = [dense_fuel, None, light_fuel, dense_fuel]
        c12 = openmc.Cell(cell_id=12, region=+r0, fill=moderator)
        fuel_univ = openmc.Universe(universe_id=11, cells=[c11, c12])

--- a/tests/regression_tests/filter_distribcell/case-1/results_true.dat
+++ b/tests/regression_tests/filter_distribcell/case-1/results_true.dat
@ -3,10 +3,10 @@ k-combined:
 tally 1:
 1.548980E-02
 2.399339E-04
-1.278781E-02
-1.635280E-04
 1.426319E-02
 2.034385E-04
+1.278781E-02
+1.635280E-04
 1.018927E-02
 1.038213E-04
 tally 2:
--- a/tests/regression_tests/filter_distribcell/case-3/results_true.dat
+++ b/tests/regression_tests/filter_distribcell/case-3/results_true.dat
@ -1 +1 @@
-386147796cf64c908002a81bc47af6763a14811ba6c457ca790653ceb8ef1c6b0376b94783f9746baf7389f18321b788c72758ec56ffd05fa3146139e6e53308
+a1f6ccf0bef1075ffc12e7fa058a44c83360801ba6e2c76f00e1bd7e9543fe1832777b122321dcb1cfb5fc1165ce84f711a65ca83903dd57cd09652e003daf62
--- a/tests/regression_tests/multipole/inputs_true.dat
+++ b/tests/regression_tests/multipole/inputs_true.dat
@ -1,7 +1,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <geometry>
  <cell id="1" material="1" universe="1" />
-  <cell id="11" material="2" region="-1" temperature="500 0 700 800" universe="11" />
+  <cell id="11" material="2" region="-1" temperature="500 700 0 800" universe="11" />
  <cell id="12" material="1" region="1" universe="11" />
  <cell fill="101" id="101" region="2 -3 4 -5" universe="0" />
  <lattice id="101">
--- a/tests/regression_tests/multipole/results_true.dat
+++ b/tests/regression_tests/multipole/results_true.dat
@ -37,5 +37,5 @@ Cell
 	Fill           =	Material 2
 	Region         =	-1
 	Rotation       =	None
-	Temperature    =	[ 500.    0.  700.  800.]
+	Temperature    =	[ 500.  700.    0.  800.]
 	Translation    =	None
--- a/tests/regression_tests/multipole/test.py
+++ b/tests/regression_tests/multipole/test.py
@ -29,7 +29,7 @@ def make_model():

    r0 = openmc.ZCylinder(R=0.3)
    c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0)
-    c11.temperature = [500, 0, 700, 800]
+    c11.temperature = [500, 700, 0, 800]
    c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0)
    fuel_univ = openmc.Universe(universe_id=11, cells=(c11, c12))

--- a/tests/regression_tests/tally_aggregation/results_true.dat
+++ b/tests/regression_tests/tally_aggregation/results_true.dat
@ -1 +1 @@
-eb2002dd2f3016154e7014501629227eb2e5b2a9655b5c0cb3b9d4d681f918fae4197bf4756fe9865c8d34f392b981b287a15e9b2fab5a46b2a5b8a33a8ae770
+8294c7481b3433a4beb25541ae72c1ea915def0c0d0f393f7585371877187f4a2a25e70bb602232e2bc1e877e78db0e9384591e5c71575659f95abb10fae0af3
--- a/tests/regression_tests/test_deplete_full.py
+++ b/tests/regression_tests/test_deplete_full.py
@ -22,6 +22,7 @@ def test_full(run_in_tmpdir):

    This test runs a complete OpenMC simulation and tests the outputs.
    It will take a while.
+
    """

    n_rings = 2
--- a/tests/regression_tests/test_reference.h5
+++ b/tests/regression_tests/test_reference.h5
--- a/tests/unit_tests/test_capi.py
+++ b/tests/unit_tests/test_capi.py
@ -26,6 +26,15 @@ def pincell_model():
    mat_tally.scores = ['total', 'elastic', '(n,gamma)']
    pincell.tallies.append(mat_tally)

+    # Add an expansion tally
+    zernike_tally = openmc.Tally()
+    filter3 = openmc.ZernikeFilter(5, r=.63)
+    cells = pincell.geometry.root_universe.cells
+    filter4 = openmc.CellFilter(list(cells.values()))
+    zernike_tally.filters = [filter3, filter4]
+    zernike_tally.scores = ['fission']
+    pincell.tallies.append(zernike_tally)
+
    # Write XML files in tmpdir
    with cdtemp():
        pincell.export_to_xml()
@ -132,7 +141,7 @@ def test_settings(capi_init):
 def test_tally_mapping(capi_init):
    tallies = openmc.capi.tallies
    assert isinstance(tallies, Mapping)
-    assert len(tallies) == 1
+    assert len(tallies) == 2
    for tally_id, tally in tallies.items():
        assert isinstance(tally, openmc.capi.Tally)
        assert tally_id == tally.id
@ -168,6 +177,14 @@ def test_tally(capi_init):
    t.active = True
    assert t.active

+    t2 = openmc.capi.tallies[2]
+    t2.id = 2
+    assert len(t2.filters) == 2
+    assert isinstance(t2.filters[0], openmc.capi.ZernikeFilter)
+    assert isinstance(t2.filters[1], openmc.capi.CellFilter)
+    assert len(t2.filters[1].bins) == 3
+    assert t2.filters[0].order == 5
+

 def test_new_tally(capi_init):
    with pytest.raises(exc.AllocationError):
@ -176,7 +193,7 @@ def test_new_tally(capi_init):
    new_tally.scores = ['flux']
    new_tally_with_id = openmc.capi.Tally(10)
    new_tally_with_id.scores = ['flux']
-    assert len(openmc.capi.tallies) == 3
+    assert len(openmc.capi.tallies) == 4


 def test_tally_results(capi_run):
@ -187,6 +204,10 @@ def test_tally_results(capi_run):
    assert np.all(t.std_dev[nonzero] >= 0)
    assert np.all(t.ci_width()[nonzero] >= 1.95*t.std_dev[nonzero])

+    t2 = openmc.capi.tallies[2]
+    n = 5
+    assert t2.mean.size == (n + 1) * (n + 2) // 2 * 3 # Number of Zernike coeffs * 3 cells
+

 def test_global_tallies(capi_run):
    assert openmc.capi.num_realizations() == 5
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@ -22,8 +22,8 @@ def test_get_atoms(res):
    n_ref = [6.6747328233649218e+08, 3.5421791038348462e+14,
             3.6208592242443462e+14, 3.3799758969347038e+14]

-    np.testing.assert_array_equal(t, t_ref)
-    np.testing.assert_array_equal(n, n_ref)
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(n, n_ref)

 def test_get_reaction_rate(res):
    """Tests evaluating reaction rate."""
@ -35,8 +35,8 @@ def test_get_reaction_rate(res):
    xs_ref = np.array([4.0594392323131994e-05, 3.9249546927524987e-05,
                       3.8394587728581798e-05, 4.1521845978371697e-05])

-    np.testing.assert_array_equal(t, t_ref)
-    np.testing.assert_array_equal(r, n_ref * xs_ref)
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(r, n_ref * xs_ref)


 def test_get_eigenvalue(res):
@ -47,5 +47,5 @@ def test_get_eigenvalue(res):
    k_ref = [1.181281798790367, 1.1798750921988739, 1.1965943696058159,
             1.2207119847790813]

-    np.testing.assert_array_equal(t, t_ref)
-    np.testing.assert_array_equal(k, k_ref)
+    np.testing.assert_allclose(t, t_ref)
+    np.testing.assert_allclose(k, k_ref)
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@ -0,0 +1,212 @@
+import numpy as np
+import scipy as sp
+
+import openmc
+import openmc.capi
+
+
+def test_t_percentile():
+    # Permutations include 1 DoF, 2 DoF, and > 2 DoF
+    # We will test 5 p-values at 3-DoF values
+    test_ps = [0.02, 0.4, 0.5, 0.6, 0.98]
+    test_dfs = [1, 2, 5]
+
+    # The reference solutions come from Scipy
+    ref_ts = [[sp.stats.t.ppf(p, df) for p in test_ps] for df in test_dfs]
+
+    test_ts = [[openmc.capi.math.t_percentile(p, df) for p in test_ps]
+               for df in test_dfs]
+
+    # The 5 DoF approximation in openmc.capi.math.t_percentile is off by up to
+    # 8e-3 from the scipy solution, so test that one separately with looser
+    # tolerance
+    assert np.allclose(ref_ts[:-1], test_ts[:-1])
+    assert np.allclose(ref_ts[-1], test_ts[-1], atol=1e-2)
+
+
+def test_calc_pn():
+    max_order = 10
+    test_xs = np.linspace(-1., 1., num=5, endpoint=True)
+
+    # Reference solutions from scipy
+    ref_vals = np.array([sp.special.eval_legendre(n, test_xs)
+                         for n in range(0, max_order + 1)])
+
+    test_vals = []
+    for x in test_xs:
+        test_vals.append(openmc.capi.math.calc_pn(max_order, x).tolist())
+
+    test_vals = np.swapaxes(np.array(test_vals), 0, 1)
+
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_evaluate_legendre():
+    max_order = 10
+    # Coefficients are set to 1, but will incorporate the (2l+1)/2 norm factor
+    # for the reference solution
+    test_coeffs = [0.5 * (2. * l + 1.) for l in range(max_order + 1)]
+    test_xs = np.linspace(-1., 1., num=5, endpoint=True)
+
+    ref_vals = np.polynomial.legendre.legval(test_xs, test_coeffs)
+
+    # Set the coefficients back to 1s for the test values since
+    # evaluate legendre incorporates the (2l+1)/2 term on its own
+    test_coeffs = [1. for l in range(max_order + 1)]
+
+    test_vals = np.array([openmc.capi.math.evaluate_legendre(test_coeffs, x)
+                          for x in test_xs])
+
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_calc_rn():
+    max_order = 10
+    test_ns = np.array([i for i in range(0, max_order + 1)])
+    azi = 0.1  # Longitude
+    pol = 0.2  # Latitude
+    test_uvw = np.array([np.sin(pol) * np.cos(azi),
+                         np.sin(pol) * np.sin(azi),
+                         np.cos(pol)])
+
+    # Reference solutions from the equations
+    ref_vals = []
+
+    def coeff(n, m):
+        return np.sqrt((2. * n + 1) * sp.special.factorial(n - m) /
+                       (sp.special.factorial(n + m)))
+
+    def pnm_bar(n, m, mu):
+        val = coeff(n, m)
+        if m != 0:
+            val *= np.sqrt(2.)
+        val *= sp.special.lpmv([m], [n], [mu])
+        return val[0]
+
+    ref_vals = []
+    for n in test_ns:
+        for m in range(-n, n + 1):
+            if m < 0:
+                ylm = pnm_bar(n, np.abs(m), np.cos(pol)) * \
+                    np.sin(np.abs(m) * azi)
+            else:
+                ylm = pnm_bar(n, m, np.cos(pol)) * np.cos(m * azi)
+
+            # Un-normalize for comparison
+            ylm /= np.sqrt(2. * n + 1.)
+            ref_vals.append(ylm)
+
+    test_vals = []
+    test_vals = openmc.capi.math.calc_rn(max_order, test_uvw)
+
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_calc_zn():
+    n = 10
+    rho = 0.5
+    phi = 0.5
+
+    # Reference solution from running the Fortran implementation
+    ref_vals = np.array([
+        1.00000000e+00, 2.39712769e-01, 4.38791281e-01,
+        2.10367746e-01, -5.00000000e-01, 1.35075576e-01,
+        1.24686873e-01, -2.99640962e-01, -5.48489101e-01,
+        8.84215021e-03, 5.68310892e-02, -4.20735492e-01,
+        -1.25000000e-01, -2.70151153e-01, -2.60091773e-02,
+        1.87022545e-02, -3.42888902e-01, 1.49820481e-01,
+        2.74244551e-01, -2.43159131e-02, -2.50357380e-02,
+        2.20500013e-03, -1.98908812e-01, 4.07587508e-01,
+        4.37500000e-01, 2.61708929e-01, 9.10321205e-02,
+        -1.54686328e-02, -2.74049397e-03, -7.94845816e-02,
+        4.75368705e-01, 7.11647284e-02, 1.30266162e-01,
+        3.37106977e-02, 1.06401886e-01, -7.31606787e-03,
+        -2.95625975e-03, -1.10250006e-02, 3.55194307e-01,
+        -1.44627826e-01, -2.89062500e-01, -9.28644588e-02,
+        -1.62557358e-01, 7.73431638e-02, -2.55329539e-03,
+        -1.90923851e-03, 1.57578403e-02, 1.72995854e-01,
+        -3.66267690e-01, -1.81657333e-01, -3.32521518e-01,
+        -2.59738162e-02, -2.31580576e-01, 4.20673902e-02,
+        -4.11710546e-04, -9.36449487e-04, 1.92156884e-02,
+        2.82515641e-02, -3.90713738e-01, -1.69280296e-01,
+        -8.98437500e-02, -1.08693628e-01, 1.78813094e-01,
+        -1.98191857e-01, 1.65964201e-02, 2.77013853e-04])
+
+    test_vals = openmc.capi.math.calc_zn(n, rho, phi)
+
+    assert np.allclose(ref_vals, test_vals)
+
+
+def test_rotate_angle():
+    uvw0 = np.array([1., 0., 0.])
+    phi = 0.
+    mu = 0.
+
+    # reference: mu of 0 pulls the vector the bottom, so:
+    ref_uvw = np.array([0., 0., -1.])
+
+    test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi)
+
+    assert np.array_equal(ref_uvw, test_uvw)
+
+    # Repeat for mu = 1 (no change)
+    mu = 1.
+    ref_uvw = np.array([1., 0., 0.])
+
+    test_uvw = openmc.capi.math.rotate_angle(uvw0, mu, phi)
+
+    assert np.array_equal(ref_uvw, test_uvw)
+
+    # Now to test phi is None
+    mu = 0.9
+    settings = openmc.capi.settings
+    settings.seed = 1
+
+    # When seed = 1, phi will be sampled as 1.9116495709698769
+    # The resultant reference is from hand-calculations given the above
+    ref_uvw = [0.9, 0.410813051297112, 0.1457142302040]
+    test_uvw = openmc.capi.math.rotate_angle(uvw0, mu)
+
+    assert np.allclose(ref_uvw, test_uvw)
+
+
+def test_maxwell_spectrum():
+    settings = openmc.capi.settings
+    settings.seed = 1
+    T = 0.5
+    ref_val = 0.6129982175261098
+    test_val = openmc.capi.math.maxwell_spectrum(T)
+
+    assert ref_val == test_val
+
+
+def test_watt_spectrum():
+    settings = openmc.capi.settings
+    settings.seed = 1
+    a = 0.5
+    b = 0.75
+    ref_val = 0.6247242713640233
+    test_val = openmc.capi.math.watt_spectrum(a, b)
+
+    assert ref_val == test_val
+
+
+def test_broaden_wmp_polynomials():
+    # Two branches of the code to worry about, beta > 6 and otherwise
+    # beta = sqrtE * dopp
+    # First lets do beta > 6
+    test_E = 0.5
+    test_dopp = 100.  # approximately U235 at room temperature
+    n = 6
+
+    ref_val = [2., 1.41421356, 1.0001, 0.70731891, 0.50030001, 0.353907]
+    test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n)
+
+    assert np.allclose(ref_val, test_val)
+
+    # now beta < 6
+    test_dopp = 5.
+    ref_val = [1.99999885, 1.41421356, 1.04, 0.79195959, 0.6224, 0.50346003]
+    test_val = openmc.capi.math.broaden_wmp_polynomials(test_E, test_dopp, n)
+
+    assert np.allclose(ref_val, test_val)