From 9946db547c782521a0866d60075df39e39713791 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Fri, 30 Nov 2018 09:03:54 +0000
Subject: [PATCH 1/8] First commit of Muir/Gaussian sampling

---
 include/openmc/distribution.h   |  39 ++++++++
 include/openmc/math_functions.h |  32 +++++++
 openmc/stats/univariate.py      | 159 +++++++++++++++++++++++++++++++-
 src/distribution.cpp            |  43 +++++++++
 src/math_functions.cpp          |  29 ++++++
 5 files changed, 301 insertions(+), 1 deletion(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 6db3c1618..58fd7bedf 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -98,6 +98,45 @@ private:
   double b_; //!< Factor in square root [1/eV]
 };
 
+//==============================================================================
+//! Gaussian (fusion) spectrum with form 1/std_dev*sqrt(pi) exp (-(e-E0)/std_dev)^2
+//==============================================================================
+
+class Gaussian : public Distribution {
+public:
+  explicit Gaussian(pugi::xml_node node);
+  Gaussian(double mean, double std_dev) : mean_{mean}, std_dev_{std_dev} { };
+
+  //! Sample a value from the distribution
+  //! \return Sampled value
+  double sample() const;
+private:
+  double mean_;    //!< middle of distribution [eV]
+  double std_dev_; //!< standard deviation [eV]
+};
+
+//==============================================================================
+//! Muit (fusion) spectrum derived from Gaussian with extra params e0 is mean
+//! std dev is sqrt(4*e0*kt/m)
+//==============================================================================
+
+class Muir : public Distribution {
+public:
+  explicit Muir(pugi::xml_node node);
+  Muir(double e0, double m_rat, double kt) : e0_{e0}, m_rat_{m_rat}, kt_{kt} { };
+
+  //! Sample a value from the distribution
+  //! \return Sampled value
+  double sample() const;
+private:
+  // example DT fusion m_rat = 5 (D = 2 + T = 3)
+  // ion temp = 20000 eV
+  // 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]
+};
+
 //==============================================================================
 //! Histogram or linear-linear interpolated tabular distribution
 //==============================================================================
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index ca0646961..b1faa4f92 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -161,6 +161,38 @@ extern "C" double maxwell_spectrum(double T);
 
 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
+//! Its sampled according to 
+//! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
+//! section 33.4.4
+//!
+//! @param mean mean of the Gaussian distribution
+//! @param std_dev standard deviation of the Gaussian distribution
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double gaussian_spectrum(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 
+//! modifyable parameters
+//!
+//! @param e0 peak neutron energy [ev]
+//! @param m_rat ratio of the fusion reactants to AMU
+//! @param kt the ion temperature of the reactants [ev]
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double muir_spectrum(double e0, double m_rat, double kt);
+
+
+
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.
 //!
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 16a9dbb9e..0ff9c4fa7 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -239,7 +239,7 @@ class Maxwell(Univariate):
 
 
 class Watt(Univariate):
-    r"""Watt fission energy spectrum.
+    """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,6 +308,163 @@ class Watt(Univariate):
         element.set("parameters", '{} {}'.format(self.a, self.b))
         return element
 
+class Gaussian(Univariate):
+    """Gaussian energy spectrum.
+
+    The Gaussian energy spectrum 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}`
+
+    Parameters
+    ----------
+    mean : float
+        Mean of the Gaussian distribution in units of eV
+    std_dev : float
+        Standard deviation of the Gaussian distribution in units of eV
+
+    Attributes
+    ----------
+    mean : flo1at
+        Mean of the Gaussian distribution in units of eV
+    std_dev : float
+        Standard deviation of the Gaussian distribution in units of eV
+    """
+
+    def __init__(self, mean=14.08e6, std_dev=4.74636e5):
+        super().__init__()
+        self.mean = mean
+        self.std_dev = std_dev
+
+    def __len__(self):
+        return 2
+
+    @property
+    def mean(self):
+        return self._mean
+
+    @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
+
+    @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)
+        self._std_dev = std_dev
+
+    def to_xml_element(self, element_name):
+        """Return XML representation of the Watt distribution
+
+        Parameters
+        ----------
+        element_name : str
+            XML element name
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing Watt distribution data
+
+        """
+        element = ET.Element(element_name)
+        element.set("type", "gaussian")
+        element.set("parameters", '{} {}'.format(self.mean, self.std_dev))
+        return element
+
+class Muir(Univariate):
+    """Muir energy spectrum.
+
+    The Muir energy spectrum is a Gaussian spectrum, but for 
+    convenience reasons allows the user 3 parameters to define
+    the distribution, e0 the mean energy of particles, the mass
+    of reactants m_rat, and the ion temperature kt. 
+
+    Parameters
+    ----------
+    e0 : float
+        Mean of the Muir distribution in units of eV
+    m_rat : float
+        Ratio of the sum of the masses of the reaction inputs to an 
+        AMU
+    kt : float
+         Ion temperature for the Muir distribution in units of eV
+
+    Attributes
+    ----------
+    e0 : float
+        Mean of the Muir distribution in units of eV
+    m_rat : float
+        Ratio of the sum of the masses of the reaction inputs to an 
+        AMU
+    kt : float
+         Ion temperature for the Muir distribution in units of eV
+
+    """
+
+    def __init__(self, e0=14.08e6, m_rat = 5., kt = 20000.):
+        super().__init__()
+        self.e0 = e0
+        self.m_rat = m_rat
+        self.kt = kt
+
+    def __len__(self):
+        return 3
+
+    @property
+    def e0(self):
+        return self._e0
+
+    @property
+    def m_rat(self):
+        return self._m_rat
+
+    @property
+    def kt(self):
+        return self._kt
+
+    @e0.setter
+    def e0(self, e0):
+        cv.check_type('Muir e0', e0, Real)
+        cv.check_greater_than('Muir e0', e0, 0.0)
+        self._e0 = e0
+
+    @m_rat.setter
+    def m_rat(self, m_rat):
+        cv.check_type('Muir m_rat', m_rat, Real)
+        cv.check_greater_than('Muir m_rat', m_rat, 0.0)
+        self._m_rat = m_rat
+
+    @kt.setter
+    def kt(self, kt):
+        cv.check_type('Muir kt', kt, Real)
+        cv.check_greater_than('Muir kt', kt, 0.0)
+        self._kt = kt
+
+    def to_xml_element(self, element_name):
+        """Return XML representation of the Watt distribution
+
+        Parameters
+        ----------
+        element_name : str
+            XML element name
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing Watt distribution data
+
+        """
+        element = ET.Element(element_name)
+        element.set("type", "gaussian")
+        element.set("parameters", '{} {} {}'.format(self.mean, self.std_dev, self.kt))
+        return element
+
 
 class Tabular(Univariate):
     """Piecewise continuous probability distribution.
diff --git a/src/distribution.cpp b/src/distribution.cpp
index e546b11dd..f5a621df7 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -113,6 +113,45 @@ double Watt::sample() const
   return watt_spectrum(a_, b_);
 }
 
+//==============================================================================
+// Gaussian implementation
+//==============================================================================
+Gaussian::Gaussian(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 "
+                        "parameters specified.");
+
+  mean_ = params.at(0);
+  std_dev_ = params.at(1);
+}
+
+double Gaussian::sample() const
+{
+  return gaussian_spectrum(mean_, std_dev_);
+}
+
+//==============================================================================
+// Muir implementation
+//==============================================================================
+Muir::Muir(pugi::xml_node node) 
+{
+  auto params = get_node_array<double>(node,"parameters");
+  if (params.size() != 3)
+    openmc::fatal_error("Muir energy distribution must have three "
+                        "parameters specified.");
+
+  e0_ = params.at(0);
+  m_rat_ = params.at(1);
+  kt_ = params.at(2);
+}
+
+double Muir::sample() const
+{
+  return muir_spectrum(e0_, m_rat_, kt_);
+}
+
 //==============================================================================
 // Tabular implementation
 //==============================================================================
@@ -256,6 +295,10 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
     dist = UPtrDist{new Maxwell(node)};
   } else if (type == "watt") {
     dist = UPtrDist{new Watt(node)};
+  } else if (type == "gaussian") {
+    dist = UPtrDist{new Gaussian(node)};
+  } else if (type == "muir") {
+    dist = UPtrDist{new Muir(node)};
   } else if (type == "discrete") {
     dist = UPtrDist{new Discrete(node)};
   } else if (type == "tabular") {
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index a8646f09d..edbbdbe32 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -674,6 +674,35 @@ double maxwell_spectrum(double T) {
 }
 
 
+double gaussian_spectrum(double mean, double standard_deviation) {
+  double r2 = 1.01;
+  double v1 = 0.;
+  double v2 = 0.;
+  
+  // perhaps there should be a limit to the number of resamples
+  while ( r2 > 1 ) {
+    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;
+      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);
+}  
+
+
 double watt_spectrum(double a, double b) {
   double w = maxwell_spectrum(a);
   double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);

From 7a04968845d33801e211a14c219eab74c0550bed Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Thu, 13 Dec 2018 15:10:00 +0000
Subject: [PATCH 2/8] Addressed P. Romano review comments

---
 include/openmc/distribution.h   | 14 ++++-----
 include/openmc/math_functions.h |  8 +++--
 openmc/stats/univariate.py      | 54 ++++++++++++++++-----------------
 src/distribution.cpp            | 14 ++++-----
 src/math_functions.cpp          | 32 +++++++++----------
 5 files changed, 61 insertions(+), 61 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 58fd7bedf..fecddc346 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 b1faa4f92..fc801813b 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 0ff9c4fa7..f1c73adec 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 f5a621df7..c74d28680 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 edbbdbe32..4c675927b 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);
 }  
 
 

From a13d80e2bdd27e6d183c16a612fba7dac3b217f6 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Thu, 13 Dec 2018 15:37:47 +0000
Subject: [PATCH 3/8] Missing one review comment, and followed style for
 generic normal distribution from uniform distribuion

---
 docs/source/pythonapi/stats.rst | 2 ++
 openmc/stats/univariate.py      | 8 ++++----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst
index f5551a681..6ad0e6b03 100644
--- a/docs/source/pythonapi/stats.rst
+++ b/docs/source/pythonapi/stats.rst
@@ -20,6 +20,8 @@ Univariate Probability Distributions
    openmc.stats.Tabular
    openmc.stats.Legendre
    openmc.stats.Mixture
+   openmc.stats.Normal
+   openmc.stats.Muir
 
 Angular Distributions
 ---------------------
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index f1c73adec..f1565dbbf 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -318,16 +318,16 @@ class Normal(Univariate):
     Parameters
     ----------
     mean_value : float
-        Mean value of the  distribution [dimensionless]
+        Mean value of the  distribution
     std_dev : float
-        Standard deviation of the Normal distribution [dimensionsless]
+        Standard deviation of the Normal distribution
 
     Attributes
     ----------
     mean_value : float
-        Mean of the Normal distribution [dimensionless]
+        Mean of the Normal distribution 
     std_dev : float
-        Standard deviation of the Normal distribution [dimensionless]
+        Standard deviation of the Normal distribution 
     """
 
     def __init__(self, mean_value, std_dev):

From e1cdd5a5c20145f1b1401bf7c59f02fde5e266a3 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Thu, 13 Dec 2018 16:38:20 +0000
Subject: [PATCH 4/8] Added unit tests for the Normal & Muir distributions in
 test_stats.py

---
 tests/unit_tests/test_stats.py | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index 388408bb2..553f2410e 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -177,3 +177,25 @@ def test_point():
     assert elem.tag == 'space'
     assert elem.attrib['type'] == 'point'
     assert elem.find('parameters') is not None
+
+def test_normal():
+    mean = 10.0
+    std_dev = 2.0
+    d = openmc.stats.Normal(mean,std_dev)
+    assert d.mean_value == pytest.approx(mean)
+    assert d.std_dev == pytest.approx(std_dev)
+    assert len(d) == 2
+    elem = d.to_xml_element('distribution')
+    assert elem.attrib['type'] == 'normal'
+
+def test_muir():
+    mean = 10.0
+    mass = 5.0
+    temp = 20000.
+    d = openmc.stats.Muir(mean,mass,temp)
+    assert d.e0 == pytest.approx(mean)
+    assert d.m_rat == pytest.approx(mass)
+    assert d.kt == pytest.approx(temp)
+    assert len(d) == 3
+    elem = d.to_xml_element('energy')
+    assert elem.attrib['type'] == 'muir'

From 32674887e5ae71044e388d143936ec2190f66913 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Thu, 13 Dec 2018 16:51:47 +0000
Subject: [PATCH 5/8] Included missing CAPI wappers and unit tests for
 test_math.py

---
 openmc/capi/math.py           | 22 ++++++++++++++++++++++
 tests/unit_tests/test_math.py | 19 +++++++++++++++++++
 2 files changed, 41 insertions(+)

diff --git a/openmc/capi/math.py b/openmc/capi/math.py
index 34beb585f..126638c9a 100644
--- a/openmc/capi/math.py
+++ b/openmc/capi/math.py
@@ -37,6 +37,8 @@ _dll.broaden_wmp_polynomials_c.restype = None
 _dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
                                            ndpointer(c_double)]
 
+_dll.normal_variate.restype = c_double
+_dll.normal_variate.argtypes = [c_double, c_double]
 
 def t_percentile(p, df):
     """ Calculate the percentile of the Student's t distribution with a
@@ -253,6 +255,26 @@ def watt_spectrum(a, b):
     return _dll.watt_spectrum(a, b)
 
 
+def normal_variate(mean_value, std_dev):
+    """ Samples an energy from the Normal distribution.
+
+    Parameters
+    ----------
+    mean_value : float
+        Mean of the Normal distribution 
+    std_dev : float
+        Standard deviation of the normal distribution
+
+    Returns
+    -------
+    float
+        Sampled outgoing normally distributed value
+
+    """
+
+    return _dll.normal_variate(mean_value, std_dev)
+
+
 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) ...
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 7ef749dec..c891d97e3 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -205,6 +205,25 @@ def test_watt_spectrum():
     assert ref_val == test_val
 
 
+def test_normal_dist():
+    settings = openmc.capi.settings
+    settings.seed = 1
+    a = 14.08
+    b = 0.0
+    ref_val = 14.08
+    test_val = openmc.capi.math.normal_variate(a, b)
+
+    assert ref_val == test_val
+
+    settings.seed = 1
+    a = 14.08
+    b = 1.0
+    ref_val = 16.436645416691427
+    test_val = openmc.capi.math.normal_variate(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

From e254cd0aea00f827338e5be965dcb3b605699b1e Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <Jonathan.Shimwell@ukaea.uk>
Date: Tue, 18 Dec 2018 12:37:44 +0000
Subject: [PATCH 6/8] fixed some of @smharper comments

---
 include/openmc/math_functions.h | 4 ++--
 src/distribution.cpp            | 2 +-
 src/math_functions.cpp          | 7 ++-----
 tests/unit_tests/test_math.py   | 4 ++--
 4 files changed, 7 insertions(+), 10 deletions(-)

diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index fc801813b..a6654dd05 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -185,9 +185,9 @@ extern "C" double normal_variate(double mean, double std_dev);
 //!  
 //! modifyable parameters
 //!
-//! @param e0 peak neutron energy [ev]
+//! @param e0 peak neutron energy [eV]
 //! @param m_rat ratio of the fusion reactants to AMU
-//! @param kt the ion temperature of the reactants [ev]
+//! @param kt the ion temperature of the reactants [eV]
 //! @result The sampled outgoing energy
 //==============================================================================
 
diff --git a/src/distribution.cpp b/src/distribution.cpp
index c74d28680..8b77d86b8 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -295,7 +295,7 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
     dist = UPtrDist{new Maxwell(node)};
   } else if (type == "watt") {
     dist = UPtrDist{new Watt(node)};
-  } else if (type == "gaussian") {
+  } else if (type == "normal") {
     dist = UPtrDist{new Normal(node)};
   } else if (type == "muir") {
     dist = UPtrDist{new Muir(node)};
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 4c675927b..849fa0677 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -677,11 +677,8 @@ double maxwell_spectrum(double T) {
 double normal_variate(double mean, double standard_deviation) {
   // perhaps there should be a limit to the number of resamples
   while ( true ) {
-    double v1 = 0.;
-    double v2 = 0.;
-
-    v1 = 2 * prn() - 1.;
-    v2 = 2 * prn() - 1.;
+    double v1 = 2 * prn() - 1.;
+    double v2 = 2 * prn() - 1.;
     
     double r = std::pow(v1, 2) + std::pow(v2, 2);
     double r2 = std::pow(r, 2);
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index c891d97e3..5b8af979a 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -213,7 +213,7 @@ def test_normal_dist():
     ref_val = 14.08
     test_val = openmc.capi.math.normal_variate(a, b)
 
-    assert ref_val == test_val
+    assert ref_val == pytest.approx(test_val)
 
     settings.seed = 1
     a = 14.08
@@ -221,7 +221,7 @@ def test_normal_dist():
     ref_val = 16.436645416691427
     test_val = openmc.capi.math.normal_variate(a, b)
 
-    assert ref_val == test_val
+    assert ref_val == pytest.approx(test_val)
 
 
 def test_broaden_wmp_polynomials():

From 3ac6c6be636db4c858d8ac61c691581913249cb0 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Tue, 18 Dec 2018 16:33:51 +0000
Subject: [PATCH 7/8] Further review comments and errors fixed

---
 include/openmc/math_functions.h | 3 +--
 openmc/stats/univariate.py      | 4 ++--
 tests/unit_tests/test_math.py   | 1 +
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index a6654dd05..6f94ee4ca 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -181,9 +181,8 @@ 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 
+//! modifiable parameters
 //! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
-//!  
-//! modifyable parameters
 //!
 //! @param e0 peak neutron energy [eV]
 //! @param m_rat ratio of the fusion reactants to AMU
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index f1565dbbf..f0ce6f83f 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -312,8 +312,8 @@ class Normal(Univariate):
     r"""Normally distributed sampling.
 
     The Normal Distribution is characterized by two parameters
-    :math:`\mu` and :math:`\sigma` and has density function :math:
-    `p(X) dX = 1/\2*sigma\sqrt{\pi} * e^{(X-\mu/\2*sigma}`
+    :math:`\mu` and :math:`\sigma` and has density function 
+    :math:`p(X) dX = 1/\2*sigma\sqrt{\pi} * e^{(X-\mu/\2*sigma}`
 
     Parameters
     ----------
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 5b8af979a..1397d87eb 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -4,6 +4,7 @@ import scipy as sp
 import openmc
 import openmc.capi
 
+import pytest
 
 def test_t_percentile():
     # Permutations include 1 DoF, 2 DoF, and > 2 DoF

From 9d35b2beb06d6b4f6e171116e350b61dd69dd533 Mon Sep 17 00:00:00 2001
From: Andrew Davis <andrew.davis@ukaea.uk>
Date: Mon, 7 Jan 2019 07:48:31 +0000
Subject: [PATCH 8/8] Last review comment

---
 openmc/stats/univariate.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index f0ce6f83f..c172c7c64 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -313,7 +313,7 @@ class Normal(Univariate):
 
     The Normal Distribution is characterized by two parameters
     :math:`\mu` and :math:`\sigma` and has density function 
-    :math:`p(X) dX = 1/\2*sigma\sqrt{\pi} * e^{(X-\mu/\2*sigma}`
+    :math:`p(X) dX = 1/(\sqrt{2\pi}\sigma) e^{-(X-\mu)^2/(2\sigma^2)}`
 
     Parameters
     ----------