Add uniform_distribution function for clarity

include/openmc/random_dist.h

@@ -5,6 +5,17 @@
 
 namespace openmc {
 
+//==============================================================================
+//! Sample a uniform distribution [a, b)
+//
+//! \param a Lower bound of uniform distribution
+//! \param b Upper bound of uniform distribtion
+//! \param seed A pointer to the pseudorandom seed
+//! \return Sampled variate
+//==============================================================================
+
+double uniform_distribution(double a, double b, uint64_t* seed);
+
 //==============================================================================
 //! Samples an energy from the Maxwell fission distribution based on a direct
 //! sampling scheme.

src/distribution_multi.cpp

@@ -6,6 +6,7 @@
 #include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
 
@@ -71,8 +72,8 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
 
 Direction isotropic_direction(uint64_t* seed)
 {
-  double phi = 2.0*PI*prn(seed);
-  double mu = 2.0*prn(seed) - 1.0;
+  double phi = uniform_distribution(0., 2.0*PI, seed);
+  double mu = uniform_distribution(-1., 1., seed);
   return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi),
       std::sqrt(1.0 - mu*mu) * std::sin(phi)};
 }

src/photon.cpp

@@ -6,6 +6,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
@@ -649,13 +650,13 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
   // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
   double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV))
     / (*E_electron + MASS_ELECTRON_EV)  ;
-  double rn = 2.0*prn(seed) - 1.0;
+  double rn = uniform_distribution(-1., 1., seed);
   *mu_electron = (rn + beta)/(rn*beta + 1.0);
 
   // Sample the scattering angle of the positron
   beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV))
     / (*E_positron + MASS_ELECTRON_EV);
-  rn = 2.0*prn(seed) - 1.0;
+  rn = uniform_distribution(-1., 1., seed);
   *mu_positron = (rn + beta)/(rn*beta + 1.0);
 }
 
@@ -726,7 +727,7 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t* seed)
     while (true) {
       if (prn(seed) < t) {
         // Left branch of flow chart
-        double r = 2.0*prn(seed);
+        double r = uniform_distribution(0.0, 2.0, seed);
         x = 1.0 + alpha*r;
         if (prn(seed) < 4.0/x*(1.0 - 1.0/x)) {
           mu = 1 - r;

src/physics.cpp

@@ -13,6 +13,7 @@
 #include "openmc/nuclide.h"
 #include "openmc/photon.h"
 #include "openmc/physics_common.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/reaction.h"
 #include "openmc/secondary_uncorrelated.h"
@@ -294,7 +295,7 @@ void sample_photon_reaction(Particle& p)
     }
 
     // Create Compton electron
-    double phi = 2.0*PI*prn(p.current_seed());
+    double phi = uniform_distribution(0., 2.0*PI, p.current_seed());
     double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
     int electron = static_cast<int>(Particle::Type::electron);
     if (E_electron >= settings::energy_cutoff[electron]) {
@@ -355,7 +356,7 @@ void sample_photon_reaction(Particle& p)
           }
         }
 
-        double phi = 2.0*PI*prn(p.current_seed());
+        double phi = uniform_distribution(0., 2.0*PI, p.current_seed());
         Direction u;
         u.x = mu;
         u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@@ -748,7 +749,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   if (!d_->angle().empty()) {
     mu_cm = d_->angle().sample(p.E_, p.current_seed());
   } else {
-    mu_cm = 2.0*prn(p.current_seed()) - 1.0;
+    mu_cm = uniform_distribution(-1., 1., p.current_seed());
   }
 
   // Determine direction cosines in CM
@@ -977,7 +978,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
     double beta_vt = std::sqrt(beta_vt_sq);
 
     // Sample cosine of angle between neutron and target velocity
-    mu = 2.0*prn(seed) - 1.0;
+    mu = uniform_distribution(-1., 1., seed);
 
     // Determine rejection probability
     double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq -

src/random_dist.cpp

@@ -7,6 +7,12 @@
 
 namespace openmc {
 
+double uniform_distribution(double a, double b, uint64_t* seed)
+{
+    return a + (b - a)*prn(seed);
+}
+
+
 double maxwell_spectrum(double T, uint64_t* seed) {
   // Set the random numbers
   double r1 = prn(seed);

src/secondary_kalbach.cpp

@@ -10,6 +10,7 @@
 #include "xtensor/xview.hpp"
 
 #include "openmc/hdf5_interface.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
@@ -222,7 +223,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* seed)
 
   // Sampled correlated angle from Kalbach-Mann parameters
   if (prn(seed) > km_r) {
-    double T = (2.0*prn(seed) - 1.0) * std::sinh(km_a);
+    double T = uniform_distribution(-1., 1., seed) * std::sinh(km_a);
     mu = std::log(T + std::sqrt(T*T + 1.0))/km_a;
   } else {
     double r1 = prn(seed);

src/secondary_nbody.cpp

@@ -27,7 +27,7 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu,
 {
   // By definition, the distribution of the angle is isotropic for an N-body
   // phase space distribution
-  mu = 2.0*prn(seed) - 1.0;
+  mu = uniform_distribution(-1., 1., seed);
 
   // Determine E_max parameter
   double Ap = mass_ratio_;

src/secondary_uncorrelated.cpp

@@ -6,7 +6,7 @@
 
 #include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
-#include "openmc/random_lcg.h"
+#include "openmc/random_dist.h"
 
 namespace openmc {
 
@@ -59,7 +59,7 @@ UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
     mu = angle_.sample(E_in, seed);
   } else {
     // no angle distribution given => assume isotropic for all energies
-    mu = 2.0*prn(seed) - 1.0;
+    mu = uniform_distribution(-1., 1., seed);
   }
 
   // Sample outgoing energy