diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst
index b4f115854..10be9454f 100644
--- a/docs/source/pythonapi/stats.rst
+++ b/docs/source/pythonapi/stats.rst
@@ -22,7 +22,7 @@ Univariate Probability Distributions
    openmc.stats.Legendre
    openmc.stats.Mixture
    openmc.stats.Normal
-   openmc.stats.Muir
+   openmc.stats.muir
 
 Angular Distributions
 ---------------------
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 8403f2c67..d3fe958ca 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -31,7 +31,6 @@ using UPtrDist = unique_ptr<Distribution>;
 //! \return Unique pointer to distribution
 UPtrDist distribution_from_xml(pugi::xml_node node);
 
-
 //==============================================================================
 //! A discrete distribution (probability mass function)
 //==============================================================================
@@ -99,11 +98,12 @@ public:
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
   double n() const { return 1 / ninv_ - 1; }
+
 private:
   //! Store processed values in object to allow for faster sampling
   double offset_; //!< a^(n+1)
-  double span_; //!< b^(n+1) - a^(n+1)
-  double ninv_; //!< 1/(n+1)
+  double span_;   //!< b^(n+1) - a^(n+1)
+  double ninv_;   //!< 1/(n+1)
 };
 
 //==============================================================================
@@ -172,35 +172,6 @@ private:
   double std_dev_;    //!< standard deviation [eV]
 };
 
-//==============================================================================
-//! Muir (fusion) spectrum derived from Normal 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
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const;
-
-  double e0() const { return e0_; }
-  double m_rat() const { return m_rat_; }
-  double kt() const { return kt_; }
-
-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
 //==============================================================================
@@ -277,7 +248,6 @@ private:
     distribution_; //!< sub-distributions + cummulative probabilities
 };
 
-
 } // namespace openmc
 
 #endif // OPENMC_DISTRIBUTION_H
diff --git a/include/openmc/random_dist.h b/include/openmc/random_dist.h
index 9bf55aa93..4c9ab3c67 100644
--- a/include/openmc/random_dist.h
+++ b/include/openmc/random_dist.h
@@ -64,23 +64,6 @@ extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
 
 extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed);
 
-//==============================================================================
-//! 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
-//!
-//! \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 seed A pointer to the pseudorandom seed
-//! \result The sampled outgoing energy
-//==============================================================================
-
-extern "C" double muir_spectrum(
-  double e0, double m_rat, double kt, uint64_t* seed);
-
 } // namespace openmc
 
 #endif // OPENMC_RANDOM_DIST_H
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index 5298499ce..28f60de6d 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -2,7 +2,9 @@ from abc import ABC, abstractmethod
 from collections import defaultdict
 from collections.abc import Iterable
 from copy import deepcopy
+import math
 from numbers import Real
+from warnings import warn
 from xml.etree import ElementTree as ET
 
 import numpy as np
@@ -53,7 +55,9 @@ class Univariate(EqualityMixin, ABC):
         elif distribution == 'normal':
             return Normal.from_xml_element(elem)
         elif distribution == 'muir':
-            return Muir.from_xml_element(elem)
+            # Support older files where Muir had its own class
+            params = [float(x) for x in get_text(elem, 'parameters').split()]
+            return muir(*params)
         elif distribution == 'tabular':
             return Tabular.from_xml_element(elem)
         elif distribution == 'legendre':
@@ -717,119 +721,45 @@ class Normal(Univariate):
         return cls(*map(float, params))
 
 
-class Muir(Univariate):
-    """Muir energy spectrum.
+def muir(e0, m_rat, kt):
+    """Generate a 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.
+    The Muir energy spectrum is a normal distribution, but for convenience
+    reasons allows the user to specify three parameters to define the
+    distribution: the mean energy of particles ``e0``, the mass of reactants
+    ``m_rat``, and the ion temperature ``kt``.
+
+    .. versionadded:: 0.14.0
 
     Parameters
     ----------
     e0 : float
-        Mean of the Muir distribution in units of eV
+        Mean of the Muir distribution in [eV]
     m_rat : float
-        Ratio of the sum of the masses of the reaction inputs to an
-        AMU
+        Ratio of the sum of the masses of the reaction inputs to 1 amu
     kt : float
-         Ion temperature for the Muir distribution in units of eV
+         Ion temperature for the Muir distribution in [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
+    Returns
+    -------
+    openmc.stats.Normal
+        Corresponding normal distribution
 
     """
+    # https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+    std_dev = math.sqrt(4 * e0 * kt / m_rat)
+    return Normal(e0, std_dev)
 
-    def __init__(self, e0=14.08e6, m_rat = 5., kt = 20000.):
-        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
-
-    @property
-    def std_dev(self):
-        return np.sqrt(4.*self.e0*self.kt/self.m_rat)
-
-    def sample(self, n_samples=1, seed=None):
-        # Based on LANL report LA-05411-MS
-        np.random.seed(seed)
-        return np.random.normal(self.e0, self.std_dev, n_samples)
-
-    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", "muir")
-        element.set("parameters", '{} {} {}'.format(self._e0, self._m_rat, self._kt))
-        return element
-
-    @classmethod
-    def from_xml_element(cls, elem):
-        """Generate Muir distribution from an XML element
-
-        Parameters
-        ----------
-        elem : xml.etree.ElementTree.Element
-            XML element
-
-        Returns
-        -------
-        openmc.stats.Muir
-            Muir distribution generated from XML element
-
-        """
-        params = get_text(elem, 'parameters').split()
-        return cls(*map(float, params))
+# Retain deprecated name for the time being
+def Muir(*args, **kwargs):
+    # warn of name change
+    warn(
+        "The Muir(...) class has been replaced by the muir(...) function and "
+        "will be removed in a future version of OpenMC. Use muir(...) instead.",
+        FutureWarning
+    )
+    return muir(*args, **kwargs)
 
 
 class Tabular(Univariate):
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 901852e0f..c3e1cc1b3 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -165,27 +165,6 @@ double Normal::sample(uint64_t* seed) const
   return normal_variate(mean_value_, std_dev_, seed);
 }
 
-//==============================================================================
-// 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(uint64_t* seed) const
-{
-  return muir_spectrum(e0_, m_rat_, kt_, seed);
-}
-
 //==============================================================================
 // Tabular implementation
 //==============================================================================
@@ -324,8 +303,10 @@ Mixture::Mixture(pugi::xml_node node)
   double cumsum = 0.0;
   for (pugi::xml_node pair : node.children("pair")) {
     // Check that required data exists
-    if (!pair.attribute("probability")) fatal_error("Mixture pair element does not have probability.");
-    if (!pair.child("dist")) fatal_error("Mixture pair element does not have a distribution.");
+    if (!pair.attribute("probability"))
+      fatal_error("Mixture pair element does not have probability.");
+    if (!pair.child("dist"))
+      fatal_error("Mixture pair element does not have a distribution.");
 
     // cummulative sum of probybilities
     cumsum += std::stod(pair.attribute("probability").value());
@@ -381,8 +362,6 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
     dist = UPtrDist {new Watt(node)};
   } else if (type == "normal") {
     dist = UPtrDist {new Normal(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/random_dist.cpp b/src/random_dist.cpp
index c8d2380c3..1aa35a689 100644
--- a/src/random_dist.cpp
+++ b/src/random_dist.cpp
@@ -46,11 +46,4 @@ double normal_variate(double mean, double standard_deviation, uint64_t* seed)
   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);
-}
-
 } // namespace openmc
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index dbf06d0b5..e752c004f 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -376,16 +376,14 @@ def test_muir():
     mean = 10.0
     mass = 5.0
     temp = 20000.
-    d = openmc.stats.Muir(mean,mass,temp)
+    d = openmc.stats.muir(mean, mass, temp)
+    assert isinstance(d, openmc.stats.Normal)
 
     elem = d.to_xml_element('energy')
-    assert elem.attrib['type'] == 'muir'
+    assert elem.attrib['type'] == 'normal'
 
-    d = openmc.stats.Muir.from_xml_element(elem)
-    assert d.e0 == pytest.approx(mean)
-    assert d.m_rat == pytest.approx(mass)
-    assert d.kt == pytest.approx(temp)
-    assert len(d) == 3
+    d = openmc.stats.Univariate.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.Normal)
 
     # sample muir distribution
     n_samples = 10000