From fa80bfe613e6f555d68694408ac3169936ea2117 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Mon, 28 Sep 2020 08:54:36 +0100 Subject: [PATCH] Fixed erroneous normal distro sampling --- src/math_functions.cpp | 24 ++++++++---------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 860ff2e8c..c6a8a6d18 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -692,25 +692,17 @@ double maxwell_spectrum(double T, uint64_t* seed) { double normal_variate(double mean, double standard_deviation, uint64_t* seed) { - // perhaps there should be a limit to the number of resamples - while ( true ) { - double v1 = 2 * prn(seed) - 1.; - double v2 = 2 * prn(seed) - 1.; - - double r = std::pow(v1, 2) + std::pow(v2, 2); - double r2 = std::pow(r, 2); - if (r2 < 1) { - double z = std::sqrt(-2.0 * std::log(r2)/r2); - z *= (prn(seed) <= 0.5) ? v1 : v2; - return mean + standard_deviation*z; - } - } + // now correct method to sample a normal distubted number + double v1 = 2 * prn(seed) - 1.; + double v2 = 2 * prn(seed) - 1.; + double r = std::pow(v1, 2) + std::pow(v2, 2); + if ( r == 0 || r > 1 ) return normal_variate(mean,standard_deviation,seed); + double z = std::sqrt(-2.0 * std::log(r)/r); + return mean + standard_deviation*z*v1; } double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) { - // note sigma here is a factor of 2 shy of equation - // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS - double sigma = std::sqrt(2.*e0*kt/m_rat); + double sigma = std::sqrt(4.*e0*kt/m_rat); return normal_variate(e0, sigma, seed); }