1 /* ************************************************************************* *
2 bugsx - (C) Copyright 1990-1997 Joshua R. Smith (jrs@media.mit.edu)
3 http://physics.www.media.mit.edu/~jrs
4
5 (C) Copyright 1995-1997 Robert Gasch (Robert_Gasch@peoplesoft.com)
6 http://www.peoplesoft.com/peoplepages/g/robert_gasch/index.htm
7
8 Permission to use, copy, modify and distribute this software for any
9 purpose and without fee is hereby granted, provided that this copyright
10 notice appear in all copies as well as supporting documentation. All
11 work developed as a consequence of the use of this program should duly
12 acknowledge such use.
13
14 No representations are made about the suitability of this software for
15 any purpose. This software is provided "as is" without express or implied
16 warranty.
17
18 See the GNU General Public Licence for more information.
19 * ************************************************************************* */
20
21
22
23 /* ************************************************************************
24 *
25 * Gaussian Random Number Generator
26 *
27 * Donald H. House July 1, 1982
28 *
29 * This function takes as parameters real valued mean and standard-deviation,
30 * and an integer valued seed. It returns a real number which may be
31 * interpreted as a sample of a normally distributed (Gaussian) random variable
32 * with the specified mean and standard deviation. As a side effect the value
33 * of the integer seed is modified.
34 *
35 * The computational technique used is to pass a uniformly distributed random
36 * number through the inverse of the Normal Distribution function.
37 *
38 * Translated into C by Robert Allen January 25, 1989
39 * ************************************************************************ */
40
41 #include "bugsx.h"
42
43 #define ITBLMAX 20
44 #define didu (20.0 / 0.5)
45 /* Where it says 20.0, put the value of ITBLMAX */
46
47
48 /* ************************************************************************ */
49 /* **************************** noise generator ************************** */
50 /* ************************************************************************ */
noise(mean,std)51 double noise (mean, std)
52 double mean, std;
53 {
54 double u, du;
55 int index, sign;
56 double temp;
57 static double tbl[ITBLMAX+1]
58 = {0.00000E+00, 6.27500E-02, 1.25641E-01, 1.89000E-01,
59 2.53333E-01, 3.18684E-01, 3.85405E-01, 4.53889E-01,
60 5.24412E-01, 5.97647E-01, 6.74375E-01, 7.55333E-01,
61 8.41482E-01, 9.34615E-01, 1.03652E+00, 1.15048E+00,
62 1.28167E+00, 1.43933E+00, 1.64500E+00, 1.96000E+00,
63 3.87000E+00 };
64
65 u = drand48();
66 if (u >= 0.5)
67 {
68 sign = 1;
69 u = u - 0.5;
70 }
71 else
72 sign = -1;
73
74 index = (int)(didu * u);
75 #ifdef DEBUG
76 /* *** leftover as a result of foretting #include <stdlib.h> *** */
77 if (index >= ITBLMAX+1)
78 {
79 printf ("Index Sanity Check in noise() failed!!!\n");
80 printf ("Index = (int)(didu * u)\n");
81 printf ("didu = %f\n", didu);
82 printf ("u = %f\n", u);
83 printf ("index = %d\n", (int)(didu * u));
84 printf ("Exiting ...\n");
85 exit (1);
86 }
87 #endif
88 du = u - index / didu;
89 temp = mean + sign * (tbl [index] + (tbl [index + 1] - tbl [index]) *
90 du * didu) * std;
91 return (temp);
92 }
93
94
95 /* ************************************************************************ */
96 /* ************************ test the noise generator ********************* */
97 /* ************************************************************************ */
98 /*
99 void test_noise ()
100 {
101 int i;
102 double mean;
103
104 for (mean=0.0; mean < 5.0; mean += 1.0)
105 {
106 printf("mean: %g ",mean);
107 for (i=0; i<5; i++)
108 printf("%g, ", noise(mean,0.1));
109 printf("\n");
110 }
111 }
112 */
113