source: trunk/MagicSoft/Simulation/Detector/Starfield/rand_un_gen.cxx

Last change on this file was 341, checked in by petry, 25 years ago
First version of the Starfield Generator in this repository. Fully functional version using the SKY2000 star catalog.
File size: 2.3 KB
Line 
1/////////////////////////////////////////////////////////////////////////////////////////////////
2// This function is from 'Numerical Recipes in C: The Art of Scientific Computing
3// William H. Press, Brian P. Flannery, Saul A. Teulosky, William T. Vetteling; New York:
4// Cambridge University Press, 1990.'
5//
6// It returns uniform random numbers between 0 and 1.
7//
8//////////////////////////////////////////////////////////////////////////////////////////////////
9
10#define MBIG 1000000000
11#define MSEED 161803398
12#define MZ 0
13#define FAC (1.0/MBIG)
14
15
16// According to Knuth, any large MBIG, and any smaller ( but still large ) MSEED can be substituted for the above values.
17
18float rand_un_gen(long *idum_local)
19
20 //Returns a uniform random deviate between 0.0 and 1.0. Set idum_local to any negative value to initialize or reanitialize the sequence.
21
22{
23
24
25 static int inext, inextp;
26 static long ma[56]; //The value 56 is special and should not be modified
27 static int iff=0;
28 long mj, mk;
29 int i, ii, k;
30
31
32 if (*idum_local < 0 || iff == 0) { //Initialization
33 iff=1;
34 mj=MSEED-(*idum_local > 0 ? -*idum_local : *idum_local); // Initialize ma[55] using the seed idum_local and the large number MSEED.
35 mj %= MBIG;
36 ma[55]=mj;
37 mk=1;
38
39 //Now initialize the rest of the table, in a slightly random order, with numbers that are not specially random.
40
41
42
43 for (i=1; i<=54; i++) {
44 ii=(21*i) % 55;
45 ma[ii]=mk;
46 mk=mj-mk;
47 if ( mk < MZ ) mk += MBIG;
48 mj=ma[ii];
49 }
50
51 for (k=1;k<=4;k++)
52 for(i=1;i<=55;i++){
53 ma[i] -= ma[1+(i+30) % 5];
54 if(ma[i] <MZ) ma[i] +=MBIG;
55 }
56
57
58
59
60 // Prepare indices for our first generated number. The constant 31 is special.
61
62 inext=0;
63 inextp=31;
64 *idum_local=1;
65 }
66 // Here is where we start, excep on initialization.
67 // Increment inext and inextp, wrapping around 56 to 1.
68
69 if (++inext == 56) inext=1;
70 if (++inextp == 56) inextp=1;
71
72 // Generate a new random number subtractively.
73
74 mj=ma[inext]-ma[inextp];
75
76 // Be sure that it is in range.
77
78 if(mj < MZ) mj += MBIG;
79
80 //Store it and output the derived uniform deviate.
81
82 ma[inext]=mj;
83
84
85 return mj*FAC;
86}
87
88
89
90
91
92
Note: See TracBrowser for help on using the repository browser.