crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fixed erroneous normal distro sampling

Browse source
This commit is contained in:
Andrew Davis 2020-09-28 08:54:36 +01:00 • committed by Paul Romano
parent f97ed95559
commit fa80bfe613
1 changed files with 8 additions and 16 deletions
Download patch file Download diff file Expand all files Collapse all files

24
src/math_functions.cpp
View file

@ -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);
}