1 /** 2 @file 3 @ingroup nr 4 */ 5 #include <stdio.h> 6 #include <bpm/bpm_messages.h> 7 #include <bpm/bpm_nr.h> 8 9 /** 10 Sample a given Gaussian distribution using ran1 as the source of the 11 uniform deviate between 0 and 1. Nicked from numerical recipes, C7.2, p289 12 13 @param mean the mean of the gaussian 14 @param std_dev the standard deviation of the gaussian 15 16 @return a gaussian deviate, returns -DBL_MAX if the random seed is not 17 set properly before 18 */ nr_rangauss(double mean,double std_dev)19double nr_rangauss( double mean, double std_dev ) { 20 21 static int iset = 0; 22 static double gset; 23 double fac, rsq, v1, v2; 24 double val = -DBL_MAX; 25 26 if ( iset == 0 ) { 27 do { 28 v1 = 2.0*nr_ran1( &bpm_rseed ) - 1.0; 29 v2 = 2.0*nr_ran1( &bpm_rseed ) - 1.0; 30 rsq = v1 * v1 + v2 * v2; 31 } while (rsq >= 1.0 || rsq == 0 ); 32 33 fac = sqrt(-2.0 * log(rsq) / rsq); 34 gset = v1 * fac; 35 iset = 1; 36 val = (v2*fac * std_dev) + mean; 37 } else { 38 iset = 0; 39 val = (gset * std_dev) + mean; 40 } 41 42 return val; 43 } 44 45