00001 00005 #include <stdio.h> 00006 #include <bpm/bpm_messages.h> 00007 #include <bpm/bpm_nr.h> 00008 00019 double nr_rangauss( double mean, double std_dev ) { 00020 00021 static int iset = 0; 00022 static double gset; 00023 double fac, rsq, v1, v2; 00024 double val = -DBL_MAX; 00025 00026 if ( iset == 0 ) { 00027 do { 00028 v1 = 2.0*nr_ran1( &bpm_rseed ) - 1.0; 00029 v2 = 2.0*nr_ran1( &bpm_rseed ) - 1.0; 00030 rsq = v1 * v1 + v2 * v2; 00031 } while (rsq >= 1.0 || rsq == 0 ); 00032 00033 fac = sqrt(-2.0 * log(rsq) / rsq); 00034 gset = v1 * fac; 00035 iset = 1; 00036 val = (v2*fac * std_dev) + mean; 00037 } else { 00038 iset = 0; 00039 val = (gset * std_dev) + mean; 00040 } 00041 00042 return val; 00043 } 00044