Move distributions into random_dist.h/cpp

src/distribution.cpp

@@ -9,6 +9,7 @@
 
 #include "openmc/error.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
 

src/distribution_energy.cpp

@@ -9,6 +9,7 @@
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 

src/math_functions.cpp

@@ -2,6 +2,9 @@
 
 #include "Faddeeva.hh"
 
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
 namespace openmc {
 
 //==============================================================================
@@ -669,49 +672,6 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* seed
 }
 
 
-double maxwell_spectrum(double T, uint64_t* seed) {
-  // Set the random numbers
-  double r1 = prn(seed);
-  double r2 = prn(seed);
-  double r3 = prn(seed);
-
-  // determine cosine of pi/2*r
-  double c = std::cos(PI / 2. * r3);
-
-  // Determine outgoing energy
-  double E_out = -T * (std::log(r1) + std::log(r2) * c * c);
-
-  return E_out;
-}
-
-
-double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
-  // Sample a normal variate using Marsaglia's polar method
-  double x, y, r2;
-  do {
-    x = 2 * prn(seed) - 1.;
-    y = 2 * prn(seed) - 1.;
-    r2 = x*x + y*y;
-  } while (r2 > 1 || r2 == 0);
-  double z = std::sqrt(-2.0 * std::log(r2)/r2);
-  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);
-}
-
-
-double watt_spectrum(double a, double b, uint64_t* seed) {
-  double w = maxwell_spectrum(a, seed);
-  double E_out = w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w);
-
-  return E_out;
-}
-
-
 void spline(int n, const double x[], const double y[], double z[])
 {
   std::vector<double> c_new(n-1);

src/random_dist.cpp

@@ -0,0 +1,48 @@
+#include "openmc/random_dist.h"
+
+#include <cmath>
+
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+double maxwell_spectrum(double T, uint64_t* seed) {
+  // Set the random numbers
+  double r1 = prn(seed);
+  double r2 = prn(seed);
+  double r3 = prn(seed);
+
+  // determine cosine of pi/2*r
+  double c = std::cos(PI / 2. * r3);
+
+  // Determine outgoing energy
+  return -T * (std::log(r1) + std::log(r2) * c * c);
+}
+
+double watt_spectrum(double a, double b, uint64_t* seed) {
+  double w = maxwell_spectrum(a, seed);
+  return w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w);
+}
+
+
+double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
+  // Sample a normal variate using Marsaglia's polar method
+  double x, y, r2;
+  do {
+    x = 2 * prn(seed) - 1.;
+    y = 2 * prn(seed) - 1.;
+    r2 = x*x + y*y;
+  } while (r2 > 1 || r2 == 0);
+  double z = std::sqrt(-2.0 * std::log(r2)/r2);
+  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

src/random_lcg.cpp

@@ -39,9 +39,7 @@ double prn(uint64_t* seed)
   uint64_t result = (word >> 43u) ^ word;
 
   // Convert output from unsigned integer to double
-  double result_d = ldexp(result, -64);
-
-  return result_d;
+  return ldexp(result, -64);
 }
 
 //==============================================================================
@@ -51,7 +49,7 @@ double prn(uint64_t* seed)
 double future_prn(int64_t n, uint64_t seed)
 {
   uint64_t fseed = future_seed(static_cast<uint64_t>(n), seed);
-  return prn( &fseed );
+  return prn(&fseed);
 }
 
 //==============================================================================

src/secondary_nbody.cpp

@@ -5,6 +5,7 @@
 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 
 namespace openmc {