Merge remote-tracking branch 'upstream/develop' into c_cell

CMakeLists.txt

@@ -438,6 +438,7 @@ set(LIBOPENMC_CXX_SRC
   src/geometry_aux.cpp
   src/hdf5_interface.cpp
   src/lattice.cpp
+  src/math_functions.cpp
   src/message_passing.cpp
   src/plot.cpp
   src/random_lcg.cpp

openmc/capi/__init__.py

@@ -47,3 +47,4 @@ from .mesh import *
 from .filter import *
 from .tally import *
 from .settings import settings
+from .math import *

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

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 '

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):

openmc/filter_expansion.py

@@ -325,15 +325,15 @@ class ZernikeFilter(ExpansionFilter):
 
     This filter allows scores to be multiplied by Zernike polynomials of the
     particle's position normalized to a given unit circle, up to a
-    user-specified order. The Zernike polynomials follow the definition by `Noll
-    <https://doi.org/10.1364/JOSA.66.000207>`_ and are defined as
+    user-specified order. The standard Zernike polynomials follow the definition by
+    Born and Wolf, *Principles of Optics* and are defined as
 
     .. math::
-        Z_n^m(\rho, \theta) = \sqrt{2n + 2} R_n^m(\rho) \cos (m\theta), \quad m > 0
+        Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0
 
-        Z_n^{m}(\rho, \theta) = \sqrt{2n + 2} R_n^{m}(\rho) \sin (m\theta), \quad m < 0
+        Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0
 
-        Z_n^{m}(\rho, \theta) = \sqrt{n + 1} R_n^{m}(\rho), \quad m = 0
+        Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0
 
     where the radial polynomials are
 
@@ -342,7 +342,8 @@ class ZernikeFilter(ExpansionFilter):
         \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
 
     With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
-    is exactly :math:`\pi` for each polynomial.
+    is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is
+    2 if :math:`m` equals 0 and 1 otherwise.
 
     Specifying a filter with order N tallies moments for all :math:`n` from 0 to
     N and each value of :math:`m`. The ordering of the Zernike polynomial

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):

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

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
 

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

src/geometry.F90

@@ -581,9 +581,9 @@ contains
         select type (lat)
 
         type is (RectLattice)
-          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
                 offset = offset + count_universe_instances(&
                      lat % get([j-1, k-1, m-1]), univ_id)

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

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

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

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

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
 

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)

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

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

tests/regression_tests/asymmetric_lattice/results_true.dat

@@ -1 +1 @@
-3b76b468b9c0e7df6508189b75748a1b7b4f2b37486396749e1a15e536526336f70b60bb207095c39ecbd46822e8c8705ea81184a3c8546e7da09bb905d74637
+bd3fd10177bc0c7b8542f15228decf8608de013872d201e6ae885cf9b5b6d7516ddb32044ee2f7fd5b207e7d37d5ce3f03e347a2e850953bba4fc26b150c89d2

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">

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

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])
 

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:

tests/regression_tests/filter_distribcell/case-3/results_true.dat

@@ -1 +1 @@
-386147796cf64c908002a81bc47af6763a14811ba6c457ca790653ceb8ef1c6b0376b94783f9746baf7389f18321b788c72758ec56ffd05fa3146139e6e53308
+a1f6ccf0bef1075ffc12e7fa058a44c83360801ba6e2c76f00e1bd7e9543fe1832777b122321dcb1cfb5fc1165ce84f711a65ca83903dd57cd09652e003daf62

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">

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

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))
 

tests/regression_tests/tally_aggregation/results_true.dat

@@ -1 +1 @@
-eb2002dd2f3016154e7014501629227eb2e5b2a9655b5c0cb3b9d4d681f918fae4197bf4756fe9865c8d34f392b981b287a15e9b2fab5a46b2a5b8a33a8ae770
+8294c7481b3433a4beb25541ae72c1ea915def0c0d0f393f7585371877187f4a2a25e70bb602232e2bc1e877e78db0e9384591e5c71575659f95abb10fae0af3

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

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)

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)