Update to rational distribution

* store offset, span and ninv in class
  * accessor functions calculate a, b, n

include/openmc/distribution.h

@@ -87,22 +87,20 @@ private:
 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)} {};
+  Rational(double a, double b, double n) : offset_{std::pow(a, n+1)}, span_{std::pow(b, n+1) - offset_}, ninv_{1/(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_; }
+  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:
-  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
+  double offset_; //!< Lower bound of distribution
+  double span_; //!< Upper bound of distribution
+  double ninv_; //!< inverse Exponent of distribution
 };
 
 //==============================================================================

openmc/stats/univariate.py

@@ -254,7 +254,7 @@ class Rational(Univariate):
     b : float, optional
         Upper bound of the sampling interval. Defaults to unity.
     n : float, optional
-        power law exponent. Defaults to zero -> Uniform distribution
+        power law exponent. Defaults to zero, which is eqivalent to an uniform distribution
 
     Attributes
     ----------

src/distribution.cpp

@@ -72,8 +72,8 @@ Uniform::Uniform(pugi::xml_node node)
 {
   auto params = get_node_array<double>(node, "parameters");
   if (params.size() != 2) {
-    openmc::fatal_error("Uniform distribution must have two "
-                        "parameters specified.");
+    fatal_error("Uniform distribution must have two "
+                "parameters specified.");
   }
 
   a_ = params.at(0);
@@ -93,22 +93,21 @@ 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.");
+    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);
+  const double a = params.at(0);
+  const double b = params.at(1);
+  const double n = params.at(2);
+  offset_ = std::pow(a, n+1);
+  span_ = std::pow(a, n+1) - offset_;
+  ninv_ = 1/(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;
+  return std::pow(offset_ + prn(seed) * span_, ninv_);
 }
 
 //==============================================================================