diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 58fd7bedff..fecddc3461 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -99,24 +99,24 @@ private:
 };
 
 //==============================================================================
-//! Gaussian (fusion) spectrum with form 1/std_dev*sqrt(pi) exp (-(e-E0)/std_dev)^2
+//! Normal distributions with form 1/2*std_dev*sqrt(pi) exp (-(e-E0)/2*std_dev)^2
 //==============================================================================
 
-class Gaussian : public Distribution {
+class Normal : public Distribution {
 public:
-  explicit Gaussian(pugi::xml_node node);
-  Gaussian(double mean, double std_dev) : mean_{mean}, std_dev_{std_dev} { };
+  explicit Normal(pugi::xml_node node);
+  Normal(double mean_value, double std_dev) : mean_value_{mean_value}, std_dev_{std_dev} { };
 
   //! Sample a value from the distribution
   //! \return Sampled value
   double sample() const;
 private:
-  double mean_;    //!< middle of distribution [eV]
+  double mean_value_;    //!< middle of distribution [eV]
   double std_dev_; //!< standard deviation [eV]
 };
 
 //==============================================================================
-//! Muit (fusion) spectrum derived from Gaussian with extra params e0 is mean
+//! Muir (fusion) spectrum derived from Normal with extra params e0 is mean
 //! std dev is sqrt(4*e0*kt/m)
 //==============================================================================
 
@@ -134,7 +134,7 @@ private:
   // mean neutron energy 14.08e6 eV
   double e0_;    //!< mean neutron energy [eV]
   double m_rat_; //!< ratio of reactant masses relative to atomic mass unit
-  double kt_;    //!< ion temperature [ev]
+  double kt_;    //!< ion temperature [eV]
 };
 
 //==============================================================================
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index b1faa4f923..fc801813b3 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -164,8 +164,8 @@ extern "C" double watt_spectrum(double a, double b);
 //==============================================================================
 //! Samples an energy from the Gaussian energy-dependent fission distribution.
 //!
-//! Samples from a Gaussian spectrum with a given mean and standard deviation
-//! The PDF is defined as s(x) = (1/sigma*sqrt(2) * e-((mu-x)/sigma)^2
+//! Samples from a Normal distribution with a given mean and standard deviation
+//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
 //! Its sampled according to 
 //! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
 //! section 33.4.4
@@ -175,12 +175,14 @@ extern "C" double watt_spectrum(double a, double b);
 //! @result The sampled outgoing energy
 //==============================================================================
 
-extern "C" double gaussian_spectrum(double mean, double std_dev);
+extern "C" double normal_variate(double mean, double std_dev);
 
 //==============================================================================
 //! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
 //!
 //! This is another form of the Gaussian distribution but with more easily 
+//! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+//!  
 //! modifyable parameters
 //!
 //! @param e0 peak neutron energy [ev]
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 0ff9c4fa7b..f1c73adec8 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -239,7 +239,7 @@ class Maxwell(Univariate):
 
 
 class Watt(Univariate):
-    """Watt fission energy spectrum.
+    r"""Watt fission energy spectrum.
 
     The Watt fission energy spectrum is characterized by two parameters
     :math:`a` and :math:`b` and has density function :math:`p(E) dE = c e^{-E/a}
@@ -308,58 +308,58 @@ class Watt(Univariate):
         element.set("parameters", '{} {}'.format(self.a, self.b))
         return element
 
-class Gaussian(Univariate):
-    """Gaussian energy spectrum.
+class Normal(Univariate):
+    r"""Normally distributed sampling.
 
-    The Gaussian energy spectrum is characterized by two parameters
+    The Normal Distribution is characterized by two parameters
     :math:`\mu` and :math:`\sigma` and has density function :math:
-    `p(E) dE = 1/\sigma\sqrt{\pi} * e^{(E-\mu/\sigma}`
+    `p(X) dX = 1/\2*sigma\sqrt{\pi} * e^{(X-\mu/\2*sigma}`
 
     Parameters
     ----------
-    mean : float
-        Mean of the Gaussian distribution in units of eV
+    mean_value : float
+        Mean value of the  distribution [dimensionless]
     std_dev : float
-        Standard deviation of the Gaussian distribution in units of eV
+        Standard deviation of the Normal distribution [dimensionsless]
 
     Attributes
     ----------
-    mean : flo1at
-        Mean of the Gaussian distribution in units of eV
+    mean_value : float
+        Mean of the Normal distribution [dimensionless]
     std_dev : float
-        Standard deviation of the Gaussian distribution in units of eV
+        Standard deviation of the Normal distribution [dimensionless]
     """
 
-    def __init__(self, mean=14.08e6, std_dev=4.74636e5):
+    def __init__(self, mean_value, std_dev):
         super().__init__()
-        self.mean = mean
+        self.mean_value = mean_value
         self.std_dev = std_dev
 
     def __len__(self):
         return 2
 
     @property
-    def mean(self):
-        return self._mean
+    def mean_value(self):
+        return self._mean_value
 
     @property
     def std_dev(self):
         return self._std_dev
 
-    @mean.setter
-    def mean(self, mean):
-        cv.check_type('Gaussian mean', mean, Real)
-        cv.check_greater_than('Gaussian mean', mean, 0.0)
-        self._mean = mean
+    @mean_value.setter
+    def mean_value(self, mean_value):
+        cv.check_type('Normal mean_value', mean_value, Real)
+        cv.check_greater_than('Normal mean_value', mean_value, 0.0)
+        self._mean_value = mean_value
 
     @std_dev.setter
     def std_dev(self, std_dev):
-        cv.check_type('Gaussian std_dev', std_dev, Real)
-        cv.check_greater_than('Gaussian std_dev', std_dev, 0.0)
+        cv.check_type('Normal std_dev', std_dev, Real)
+        cv.check_greater_than('Normal std_dev', std_dev, 0.0)
         self._std_dev = std_dev
 
     def to_xml_element(self, element_name):
-        """Return XML representation of the Watt distribution
+        """Return XML representation of the Normal distribution
 
         Parameters
         ----------
@@ -373,8 +373,8 @@ class Gaussian(Univariate):
 
         """
         element = ET.Element(element_name)
-        element.set("type", "gaussian")
-        element.set("parameters", '{} {}'.format(self.mean, self.std_dev))
+        element.set("type", "normal")
+        element.set("parameters", '{} {}'.format(self.mean_value, self.std_dev))
         return element
 
 class Muir(Univariate):
@@ -461,8 +461,8 @@ class Muir(Univariate):
 
         """
         element = ET.Element(element_name)
-        element.set("type", "gaussian")
-        element.set("parameters", '{} {} {}'.format(self.mean, self.std_dev, self.kt))
+        element.set("type", "muir")
+        element.set("parameters", '{} {} {}'.format(self._e0, self._m_rat, self._kt))
         return element
 
 
diff --git a/src/distribution.cpp b/src/distribution.cpp
index f5a621df71..c74d286808 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -114,22 +114,22 @@ double Watt::sample() const
 }
 
 //==============================================================================
-// Gaussian implementation
+// Normal implementation
 //==============================================================================
-Gaussian::Gaussian(pugi::xml_node node) 
+Normal::Normal(pugi::xml_node node) 
 {
   auto params = get_node_array<double>(node,"parameters");
   if (params.size() != 2)
-    openmc::fatal_error("Gaussian energy distribution must have two "
+    openmc::fatal_error("Normal energy distribution must have two "
                         "parameters specified.");
 
-  mean_ = params.at(0);
+  mean_value_ = params.at(0);
   std_dev_ = params.at(1);
 }
 
-double Gaussian::sample() const
+double Normal::sample() const
 {
-  return gaussian_spectrum(mean_, std_dev_);
+  return normal_variate(mean_value_, std_dev_);
 }
 
 //==============================================================================
@@ -296,7 +296,7 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
   } else if (type == "watt") {
     dist = UPtrDist{new Watt(node)};
   } else if (type == "gaussian") {
-    dist = UPtrDist{new Gaussian(node)};
+    dist = UPtrDist{new Normal(node)};
   } else if (type == "muir") {
     dist = UPtrDist{new Muir(node)};
   } else if (type == "discrete") {
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index edbbdbe327..4c675927b6 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -674,32 +674,30 @@ double maxwell_spectrum(double T) {
 }
 
 
-double gaussian_spectrum(double mean, double standard_deviation) {
-  double r2 = 1.01;
-  double v1 = 0.;
-  double v2 = 0.;
-  
+double normal_variate(double mean, double standard_deviation) {
   // perhaps there should be a limit to the number of resamples
-  while ( r2 > 1 ) {
+  while ( true ) {
+    double v1 = 0.;
+    double v2 = 0.;
+
     v1 = 2 * prn() - 1.;
     v2 = 2 * prn() - 1.;
-    double r = std::pow(v1,2) + std::pow(v2,2);
-    r2 = std::pow(r,2);
-    if ( r2 < 1 ) {
-      double z1 = v1 * std::sqrt(-2.0 * std::log(r2)/r2);
-      double z2 = v2 * std::sqrt(-2.0 * std::log(r2)/r2); 
-      // sample each side of the bell curve
-      double z = 0;
-      if ( prn() <= 0.5 ) z = z1;
-      else z = z2;
+    
+    double r = std::pow(v1, 2) + std::pow(v2, 2);
+    double r2 = std::pow(r, 2);
+    if (r2 < 1) {
+      double z = std::sqrt(-2.0 * std::log(r2)/r2);
+      z *= (prn() <= 0.5) ? v1 : v2;
       return mean + standard_deviation*z;
     }
   }
 }
 
 double muir_spectrum(double e0, double m_rat, double kt) {
-  double sigma = std::sqrt(4.*e0*kt/m_rat);
-  return gaussian_spectrum(e0,sigma);
+  // note sigma here is a factor of 2 shy of equation
+  // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+  double sigma = std::sqrt(2.*e0*kt/m_rat);
+  return normal_variate(e0, sigma);
 }