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)19 double 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