diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index c6a8a6d18f..a3168fca1c 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -692,16 +692,19 @@ double maxwell_spectrum(double T, uint64_t* seed) {
 
 
 double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
-  // now correct method to sample a normal distubted number
-  double v1 = 2 * prn(seed) - 1.;
-  double v2 = 2 * prn(seed) - 1.;
-  double r = std::pow(v1, 2) + std::pow(v2, 2);
-  if ( r == 0 || r > 1 ) return normal_variate(mean,standard_deviation,seed);
-  double z = std::sqrt(-2.0 * std::log(r)/r);
-  return mean + standard_deviation*z*v1;
+  // Sample a normal variate using Marsaglia's polar method
+  double x, y, r2;
+  do {
+    x = 2 * prn(seed) - 1.;
+    y = 2 * prn(seed) - 1.;
+    r2 = x*x + y*y;
+  } while (r2 > 1 || r2 == 0);
+  double z = std::sqrt(-2.0 * std::log(r2)/r2);
+  return mean + standard_deviation*z*x;
 }
 
 double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) {
+  // https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
   double sigma = std::sqrt(4.*e0*kt/m_rat);
   return normal_variate(e0, sigma, seed);
 }
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 94bf51ce57..2e3f654009 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -1,10 +1,11 @@
 import numpy as np
+import pytest
 import scipy as sp
+from scipy.stats import shapiro
 
 import openmc
 import openmc.lib
 
-import pytest
 
 def test_t_percentile():
     # Permutations include 1 DoF, 2 DoF, and > 2 DoF
@@ -216,21 +217,16 @@ def test_normal_dist():
     # make sigma = 1.0
     b = 1.0
 
-    # import shapiro test 
-    from scipy.stats import shapiro
     samples = []
     num_samples = 10000
     # sample some numbers
     for i in range(num_samples):
         # sample the normal distribution from openmc
         samples.append(openmc.lib.math.normal_variate(a, b, prn_seed))
-        prn_seed = prn_seed + 1
+        prn_seed += 1
 
-    stat,p = shapiro(samples)
-    assert stat > 0.97
-
-    # cleanup
-    del samples
+    stat, p = shapiro(samples)
+    assert p > 0.05
 
 
 def test_broaden_wmp_polynomials():