diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 95e30eff3..e5ba5fdd6 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -1,39 +1,38 @@
-from ctypes import (c_int, c_double, POINTER, CFUNCTYPE)
+from ctypes import (c_int, c_double, POINTER)
 
 import numpy as np
 from numpy.ctypeslib import ndpointer
 
 from . import _dll
 
-# _dll.normal_percentile.restype = c_double
-# _dll.normal_percentile.argtypes = [POINTER(c_double)]
 _dll.t_percentile.restype = c_double
 _dll.t_percentile.argtypes = [POINTER(c_double), POINTER(c_int)]
+
 _dll.calc_pn.restype = None
 _dll.calc_pn.argtypes = [POINTER(c_int), POINTER(c_double),
                          ndpointer(c_double)]
+
 _dll.evaluate_legendre.restype = c_double
 _dll.evaluate_legendre.argtypes = [POINTER(c_int), POINTER(c_double),
                                    POINTER(c_double)]
+
 _dll.calc_rn.restype = None
 _dll.calc_rn.argtypes = [POINTER(c_int), ndpointer(c_double),
                          ndpointer(c_double)]
+
 _dll.calc_zn.restype = None
 _dll.calc_zn.argtypes = [POINTER(c_int), POINTER(c_double), POINTER(c_double),
                          ndpointer(c_double)]
+
 _dll.rotate_angle.restype = None
 _dll.rotate_angle.argtypes = [ndpointer(c_double), POINTER(c_double),
                               ndpointer(c_double), POINTER(c_double)]
 _dll.maxwell_spectrum.restype = c_double
 _dll.maxwell_spectrum.argtypes = [POINTER(c_double)]
+
 _dll.watt_spectrum.restype = c_double
 _dll.watt_spectrum.argtypes = [POINTER(c_double), POINTER(c_double)]
-# COMPLEX DATA TYPES?
-# _dll.faddeeva.restype = c_double
-# _dll.faddeeva.argtypes = [POINTER(c_double)]
 
-# _dll.w_derivative.restype = c_double
-# _dll.w_derivative.argtypes = [POINTER(c_double)]
 _dll.broaden_wmp_polynomials.restype = None
 _dll.broaden_wmp_polynomials.argtypes = [POINTER(c_double), POINTER(c_double),
                                          POINTER(c_int), ndpointer(c_double)]
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 78ff08b05..9cd52096d 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -2,6 +2,16 @@
 
 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" const double PI {3.1415926535898};
+
+extern "C" const double SQRT_PI {std::sqrt(PI)};
 
 //==============================================================================
 // NORMAL_PERCENTILE calculates the percentile of the standard normal
@@ -748,59 +758,6 @@ double watt_spectrum_c(const double a, const double b) {
   return E_out;
 }
 
-//==============================================================================
-// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
-//==============================================================================
-
-// double complex __attribute__ ((const)) faddeeva_c(double complex z) {
-//   double complex wv; // The resultant w(z) value
-//   double relerr;      // Target relative error in the inner loop of MIT Faddeeva
-
-//   // Technically, the value we want is given by the equation:
-//   // w(z) = I/Pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
-//   // as shown in Equation 63 from Hwang, R. N. "A rigorous pole
-//   // representation of multilevel cross sections and its practical
-//   // applications." Nuclear Science and Engineering 96.3 (1987): 192-209.
-//   //
-//   // The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
-//   // two forms of the Faddeeva function are related by a transformation.
-//   //
-//   // If we call the integral form w_int, and the function form w_fun:
-//   // For imag(z) > 0, w_int(z) = w_fun(z)
-//   // For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
-
-//   // Note that faddeeva_w will interpret zero as machine epsilon
-
-//   relerr = 0.;
-//   if (cimag(z) > 0.) {
-//       wv = Faddeeva_w(z, relerr);
-//   } else {
-//     wv = -conj(Faddeeva_w(conj(z), relerr));
-//   }
-
-//   return wv;
-// }
-
-// double complex w_derivative_c(double complex z, int order){
-//   double complex wv; // The resultant w(z) value
-
-//   const double complex twoi_sqrtpi = 2.0 / SQRT_PI * I;
-
-//   switch(order) {
-//     case 0:
-//       wv = faddeeva_c(z);
-//       break;
-//     case 1:
-//       wv = -2. * z * faddeeva_c(z) + twoi_sqrtpi;
-//       break;
-//     default:
-//       wv = -2. * z * w_derivative_c(z, order - 1) - 2. * (order - 1) *
-//            w_derivative_c(z, order - 2);
-//   }
-
-//   return wv;
-// }
-
 //==============================================================================
 // BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit.
 // The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ...
diff --git a/src/math_functions.h b/src/math_functions.h
index 5417281ea..431e9ddee 100644
--- a/src/math_functions.h
+++ b/src/math_functions.h
@@ -2,21 +2,22 @@
 #define MATH_FUNCTIONS_H
 
 #include <cmath>
-#include <complex.h>
-#include <iostream>
 
 #include "random_lcg.h"
-// #include "faddeeva/Faddeeva.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
-const double PI = 3.1415926535898;
+extern "C" const double PI;
 
-const double SQRT_PI = std::sqrt(PI);
+extern "C" const double SQRT_PI;
 
 //==============================================================================
 // NORMAL_PERCENTILE calculates the percentile of the standard normal
@@ -94,14 +95,6 @@ extern "C" double maxwell_spectrum_c(const double T);
 
 extern "C" double watt_spectrum_c(const double a, const double b);
 
-//==============================================================================
-// FADDEEVA the Faddeeva function, using Stephen Johnson's implementation
-//==============================================================================
-
-// extern "C" double complex faddeeva_c(double complex z) __attribute__ ((const));
-
-// extern "C" double complex w_derivative_c(double complex z, int order);
-
 //==============================================================================
 // BROADEN_WMP_POLYNOMIALS Doppler broadens the windowed multipole curvefit.
 // The curvefit is a polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ...
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 487d02b10..6825f6b93 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -5,19 +5,6 @@ import openmc
 import openmc.capi
 
 
-# def test_normal_percentile():
-#     # normal_percentile has three branches to consider:
-#     # p < 0.02425; 0.02425 <= p <= 0.97575; and p > 0.97575
-#     test_ps = [0.02, 0.4, 0.5, 0.6, 0.98]
-
-#     # The reference solutions come from Scipy
-#     ref_zs = [sp.stats.norm.ppf(p) for p in test_ps]
-
-#     test_zs = [openmc.capi.math.normal_percentile(p) for p in test_ps]
-
-#     assert np.allclose(ref_zs, test_zs)
-
-
 def test_t_percentile():
     # Permutations include 1 DoF, 2 DoF, and > 2 DoF
     # We will test 5 p-values at 3-DoF values