diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 3d22ccdac6..cd7ab40520 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -98,6 +98,30 @@ def calc_pn(n, x):
     return _dll.calc_pn(c_int(n), c_double(x))
 
 
+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_int(len(data)),
+                                  data_arr.ctypes.data_as(POINTER(c_double)),
+                                  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 (where -n <= m <= n).
@@ -151,30 +175,6 @@ def calc_zn(n, rho, phi):
     return zn
 
 
-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_int(len(data)),
-                                  data_arr.ctypes.data_as(POINTER(c_double)),
-                                  c_double(x))
-
-
 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.
diff --git a/src/math.F90 b/src/math.F90
index 40c2d09dc7..9e1f01de5a 100644
--- a/src/math.F90
+++ b/src/math.F90
@@ -148,7 +148,7 @@ contains
 ! the return value will be 1.0.
 !===============================================================================
 
-  pure function calc_pn(n,x) result(pnx) bind(C)
+  pure function calc_pn(n, x) result(pnx) bind(C)
 
     integer(C_INT), intent(in) :: n   ! Legendre order requested
     real(C_DOUBLE), intent(in) :: x   ! Independent variable the Legendre is to
@@ -157,39 +157,25 @@ contains
     real(C_DOUBLE)             :: 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
+    pnx = calc_pn_cc(n, x)
 
   end function calc_pn
 
+!===============================================================================
+! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
+! and the value of x
+!===============================================================================
+
+  pure function evaluate_legendre(n, data, x) result(val) bind(C)
+    integer(C_INT), intent(in) :: n
+    real(C_DOUBLE), intent(in) :: data(n)
+    real(C_DOUBLE), intent(in) :: x
+    real(C_DOUBLE)             :: val
+
+    val = evaluate_legendre_cc(size(data), data, x)
+
+  end function evaluate_legendre
+
 !===============================================================================
 ! 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)
@@ -202,387 +188,7 @@ contains
                                          ! assumed to be on unit sphere
     real(C_DOUBLE)             :: rn(2*n + 1) ! The resultant R_n(uvw)
 
-
-    real(C_DOUBLE) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw)
-    real(C_DOUBLE) :: 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
+    call calc_rn_cc(n, uvw, rn)
 
   end subroutine calc_rn
 
@@ -693,26 +299,6 @@ contains
     end do
   end subroutine calc_zn
 
-!===============================================================================
-! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
-! and the value of x
-!===============================================================================
-
-  pure function evaluate_legendre(n, data, x) result(val) bind(C)
-    integer(C_INT), intent(in) :: n
-    real(C_DOUBLE), intent(in) :: data(n)
-    real(C_DOUBLE), intent(in) :: x
-    real(C_DOUBLE)             :: val
-
-    integer(C_INT) :: l
-
-    val =  HALF * data(1)
-    do l = 1, n - 1
-      val = val + (real(l, 8) +  HALF) * data(l + 1) * calc_pn(l,x)
-    end do
-
-  end function evaluate_legendre
-
 !===============================================================================
 ! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
 ! mu and through an azimuthal angle sampled uniformly. Note that this is done
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 3b645c7724..6ddf944a76 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -161,6 +161,22 @@ double __attribute__ ((const)) calc_pn_c(int n, double x) {
   return pnx;
 }
 
+//==============================================================================
+// EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
+// and the value of x
+//==============================================================================
+
+double __attribute__ ((const)) evaluate_legendre_c(int n, double data[],
+                                                   double x) {
+  double val;
+
+  val = 0.5 * data[0];
+  for (int l = 1; l < n; l++) {
+    val += (static_cast<double>(l) + 0.5) * data[l] * calc_pn_c(l, x);
+  }
+  return val;
+}
+
 //==============================================================================
 // CALC_RN calculates the n-th order spherical harmonics for a given angle
 // (in terms of (u,v,w)).  All Rn,m values are provided (where -n<=m<=n)
@@ -561,21 +577,6 @@ void calc_rn_c(int n, double uvw[3], double rn[]){
   }
 }
 
-//==============================================================================
-// EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
-// and the value of x
-//==============================================================================
-
-double __attribute__ ((const)) evaluate_legendre_c(int n, double data[],
-                                                   double x) {
-  double val;
-
-  val = 0.5 * data[0];
-  for (int l = 1; l < n; l++) {
-    val += (static_cast<double>(l) + 0.5) * data[l] * calc_pn_c(l, x);
-  }
-
-}
 
 //==============================================================================
 // ROTATE_ANGLE rotates direction std::cosines through a polar angle whose
diff --git a/src/math_functions.h b/src/math_functions.h
index aa048e5f67..bea30cf916 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -32,13 +32,6 @@ extern "C" double t_percentile_c(double p, int df) __attribute__ ((const));
 
 extern "C" double calc_pn_c(int n, double x) __attribute__ ((const));
 
-//==============================================================================
-// CALC_RN calculates the n-th order spherical harmonics for a given angle
-// (in terms of (u,v,w)).  All Rn,m values are provided (where -n<=m<=n)
-//==============================================================================
-
-extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
-
 //==============================================================================
 // EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
 // and the value of x
@@ -47,6 +40,13 @@ extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
 extern "C" double evaluate_legendre_c(int n, double data[], double x)
      __attribute__ ((const));
 
+//==============================================================================
+// CALC_RN calculates the n-th order spherical harmonics for a given angle
+// (in terms of (u,v,w)).  All Rn,m values are provided (where -n<=m<=n)
+//==============================================================================
+
+extern "C" void calc_rn_c(int n, double uvw[3], double rn[]);
+
 //==============================================================================
 // ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
 // mu and through an azimuthal angle sampled uniformly. Note that this is done
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 98295aefbf..c8d33e6012 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -51,6 +51,25 @@ def test_calc_pn():
     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 includes the (2l+1)/2 term
+    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)])
@@ -99,23 +118,6 @@ def test_calc_zn():
     pass
 
 
-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
-    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_rotate_angle():
     uvw0 = np.array([1., 0., 0.])
     phi = 0.