Rational univariate distribution, that sampled a random variable with a pdf of p = x^n

Usefull for samplig homogenous cylindrical or spherical sources
Simmilar to MCNPs "sp -21 n" functionality

include/openmc/distribution.h

@@ -80,6 +80,31 @@ private:
   double b_; //!< Upper bound of distribution
 };
 
+//==============================================================================
+//! Rational distribution over the interval [a,b] with exponent n
+//==============================================================================
+
+class Rational : public Distribution {
+public:
+  explicit Rational(pugi::xml_node node);
+  Rational(double a, double b, double n) : a_{a}, b_{b}, n_{n}, an_{std::pow(a, n+1)}, bn_{std::pow(b, n+1)} {};
+
+  //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample(uint64_t* seed) const;
+
+  double a() const { return a_; }
+  double b() const { return b_; }
+  double n() const { return n_; }
+private:
+  double a_; //!< Lower bound of distribution
+  double b_; //!< Upper bound of distribution
+  double n_; //!< Exponent of distribution
+  double an_; //!< Lower bound of distribution
+  double bn_; //!< Upper bound of distribution
+};
+
 //==============================================================================
 //! Maxwellian distribution of form c*sqrt(E)*exp(-E/theta)
 //==============================================================================

openmc/stats/univariate.py

@@ -42,6 +42,8 @@ class Univariate(EqualityMixin, ABC):
             return Discrete.from_xml_element(elem)
         elif distribution == 'uniform':
             return Uniform.from_xml_element(elem)
+        elif distribution == 'rational':
+            return Rational.from_xml_element(elem)
         elif distribution == 'maxwell':
             return Maxwell.from_xml_element(elem)
         elif distribution == 'watt':
@@ -242,6 +244,101 @@ class Uniform(Univariate):
         params = get_text(elem, 'parameters').split()
         return cls(*map(float, params))
 
+class Rational(Univariate):
+    """Distribution with power law probability over a finite interval [a,b]
+
+    Parameters
+    ----------
+    a : float, optional
+        Lower bound of the sampling interval. Defaults to zero.
+    b : float, optional
+        Upper bound of the sampling interval. Defaults to unity.
+    n : float, optional
+        power law exponent. Defaults to zero -> Uniform distribution
+
+    Attributes
+    ----------
+    a : float
+        Lower bound of the sampling interval
+    b : float
+        Upper bound of the sampling interval
+    n : float
+        Power law exponent
+
+    """
+
+    def __init__(self, a=0.0, b=1.0, n=0):
+        self.a = a
+        self.b = b
+        self.n = n
+
+    def __len__(self):
+        return 3
+
+    @property
+    def a(self):
+        return self._a
+
+    @property
+    def b(self):
+        return self._b
+
+    @property
+    def n(self):
+        return self._n
+
+    @a.setter
+    def a(self, a):
+        cv.check_type('Uniform a', a, Real)
+        self._a = a
+
+    @b.setter
+    def b(self, b):
+        cv.check_type('Uniform b', b, Real)
+        self._b = b
+
+    @n.setter
+    def n(self, n):
+        cv.check_type('Uniform n', n, Real)
+        self._n = n
+
+    def to_xml_element(self, element_name):
+        """Return XML representation of the rational distribution
+
+        Parameters
+        ----------
+        element_name : str
+            XML element name
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing uniform distribution data
+
+        """
+        element = ET.Element(element_name)
+        element.set("type", "rational")
+        element.set("parameters", '{} {} {}'.format(self.a, self.b, self.n))
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate rational distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Rational
+            Uniform distribution generated from XML element
+
+        """
+        params = get_text(elem, 'parameters').split()
+        return cls(*map(float, params))
+
 
 class Maxwell(Univariate):
     r"""Maxwellian distribution in energy.

src/distribution.cpp

@@ -85,6 +85,32 @@ double Uniform::sample(uint64_t* seed) const
   return a_ + prn(seed) * (b_ - a_);
 }
 
+//==============================================================================
+// Rational implementation
+//==============================================================================
+
+Rational::Rational(pugi::xml_node node)
+{
+  auto params = get_node_array<double>(node, "parameters");
+  if (params.size() != 3) {
+    openmc::fatal_error("Rational distribution must have three "
+                        "parameters specified.");
+  }
+
+  a_ = params.at(0);
+  b_ = params.at(1);
+  n_ = params.at(2);
+  an_ = std::pow(a_, n_+1);
+  bn_ = std::pow(b_, n_+1);
+}
+
+double Rational::sample(uint64_t* seed) const
+{
+  const double u = prn(seed);
+  const double r = std::pow(an_+u*(bn_-an_), 1./(n_+1));
+  return r;
+}
+
 //==============================================================================
 // Maxwell implementation
 //==============================================================================
@@ -347,6 +373,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
   UPtrDist dist;
   if (type == "uniform") {
     dist = UPtrDist {new Uniform(node)};
+  } else if (type == "rational") {
+    dist = UPtrDist{new Rational(node)};
   } else if (type == "maxwell") {
     dist = UPtrDist {new Maxwell(node)};
   } else if (type == "watt") {