From 9eeff1ede98895c4b35a6adfd1f9791dc81dd733 Mon Sep 17 00:00:00 2001
From: Olaf Schumann <o.j.schumann@gmail.com>
Date: Tue, 5 Oct 2021 20:04:21 +0000
Subject: [PATCH] Update to rational distribution   * store offset, span and
 ninv in class   * accessor functions calculate a, b, n

---
 include/openmc/distribution.h | 16 +++++++---------
 openmc/stats/univariate.py    |  2 +-
 src/distribution.cpp          | 23 +++++++++++------------
 3 files changed, 19 insertions(+), 22 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 6497078e47..48c167810c 100644
--- a/include/openmc/distribution.h
+++ b/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
 };
 
 //==============================================================================
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index cdc2f4c300..cfbd0535bc 100644
--- a/openmc/stats/univariate.py
+++ b/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
     ----------
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 8fedbf57a6..362c237fc6 100644
--- a/src/distribution.cpp
+++ b/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_);
 }
 
 //==============================================================================