TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
random.c
Go to the documentation of this file.
1 #include "turtle.h"
2 #include "t_random.h"
3 
4 /* For the Numerical Recipes Random Generator ran1 */
5 #define PI 3.141592654
6 #define IA 16807
7 #define IM 2147483647
8 #define AM (1.0/IM)
9 #define IQ 127773
10 #define IR 2836
11 #define M1 259200
12 #define IA1 7141
13 #define IC1 54773
14 #define RM1 (1.0/M1)
15 #define M2 134456
16 #define IA2 8121
17 #define IC2 28411
18 #define RM2 (1.0/M2)
19 #define M3 243000
20 #define IA3 4561
21 #define IC3 51349
22 #define NR_END 1
23 /* For the Numerical Recipes Random Generator ran2 */
24 
25 #define IM1 2147483563
26 #define IM2 2147483399
27 #define AM1 (1.0/IM1)
28 #define IMM1 (IM1-1)
29 #define IAA1 40014
30 #define IAA2 40692
31 #define IQ1 53668
32 #define IQ2 52774
33 #define IR1 12211
34 #define IR2 3791
35 #define NTAB 32
36 #define NDIV (1+IMM1/NTAB)
37 #define EPS 1.2e-7
38 #define RNMX (1.0 - EPS)
39 
40 
41 /*--------------------------------------------------------------------------*/
42 
43 /* Quick random number genrator modified from Numerical Recipes */
44 
45 long urand(long *idum,long range)
46 {
47 static long y, maxran, v[98];
48 long dum;
49 static long iff=0;
50 long j=0;
51 long i;
52 
53 if(*idum <0 || iff==0){
54 iff=1; i=2;
55 
56 maxran=RAND_MAX+1.0;
57 srand(*idum);
58 *idum=1;
59 for(i=1;j<=97;j++) dum=rand();
60 for(i=1;i<=97;i++) v[i]=range*(long)rand()/maxran;
61 
62 y=1+97*(long)rand()/maxran;
63 
64 }
65 j=y;
66 y=v[j];
67 v[j]=range*(long)rand()/maxran;
68 
69 return y;
70 }
71 /*--------------------------------------------------------------------------*/
72 /* Another random number generator from Numerical Recipes with some
73 modifications */
74 
75 double ran1(long *idum)
76 {
77 static long ix1,ix2,ix3;
78 static double r[98];
79 double temp;
80 static long iff;
81 long j;
82 
83 if(*idum<0 || iff==0){
84  iff=1;
85  ix1=(IC1-(*idum)) % M1;
86  ix1=(IA1*ix1+IC1) % M1;
87  ix2=ix1 % M2;
88  ix1=(IA1*ix1+IC1) % M1;
89  ix3=ix1 % M3;
90  for(j=0;j<=96;j++){
91  ix1=(IA1*ix1+IC1) % M1;
92  ix2=(IA2*ix2+IC2) % M2;
93  r[j]=(ix1+ix2*RM2)*RM1;
94  }
95  *idum=1;
96 }
97 
98 ix1=(IA1*ix1+IC1) % M1;
99 ix2=(IA2*ix2+IC2) % M2;
100 ix3=(IA3*ix3+IC3) % M3;
101 j=1+((96*ix3)/M3);
102 if(j > 96 || j<0) t_error("This cannot happen.");
103 temp=r[j];
104 r[j]=(ix1+ix2*RM2)*RM1;
105 return temp;
106 
107 }
108 
109 
110 /* As above a random generator from NR */
111 
112 double ran2(long *idum)
113 
114 {
115 
116 
117 long j;
118 long k;
119 static long idum2=123456789;
120 static long iy=0;
121 static long iv[NTAB];
122 double temp;
123 
124 
125 if(*idum <=0){
126  if(-(*idum) <1) *idum=1;
127  else *idum=-(*idum);
128  idum2=(*idum);
129  for(j=NTAB+7;j>=0;j--){
130  k=(*idum)/IQ1;
131  *idum=IAA1*(*idum-k*IQ1)-k*IR1;
132  if (*idum <0) *idum +=IM1;
133  if ( j < NTAB) iv[j] =* idum;
134 
135  }
136  iy=iv[0];
137 }
138 
139 k=(*idum)/IQ1;
140 *idum=IAA1*(*idum-k*IQ1)-k*IR1;
141 if(*idum <0) *idum += IM1;
142 k=idum2/IQ2;
143 idum2=IAA2*(idum2-k*IQ2)-k*IR2;
144 if(idum2 < 0) idum += IM2;
145 j=iy/NDIV;
146 iy=iv[j]-idum2;
147 iv[j] = *idum;
148 if( iy <1 ) iy +=IMM1;
149 
150 /* printf("rnd-> %f %f\n",RNMX,temp);*/
151 if ((temp=AM1*iy) > RNMX) return RNMX;
152 else return temp;
153 
154 }
155 
156 double ran3( long *idum)
157 /* "Minimal" random number generator of Park and Miller...
158 see Press et al."Numerical Recipes in C", 1992. */
159 { int j;
160 long k;
161 static long iy=0;
162 static long iv[NTAB];
163 float temp;
164 
165 if(*idum<=0 || !iy){
166 if(-(*idum)<1) *idum=1;
167 else *idum= -(*idum);
168 for (j=NTAB+7;j>=0; j--){
169 k=(*idum)/IQ;
170 *idum=IA*(*idum-k*IQ)-IR*k;
171 if(*idum<0) *idum +=IM;
172 if(j<NTAB) iv[j]=*idum;
173 }iy=iv[0];
174 }
175 k=(*idum)/IQ;
176 *idum=IA*(*idum-k*IQ)-IR*k;
177 if (*idum<0) *idum+=IM;
178 j=iy/NDIV;
179 iy=iv[j];
180 iv[j]=*idum;
181 if((temp=AM*iy)>RNMX)return RNMX;
182 else return temp;
183 }
184 
185 /*--------------------------------------------------------------------------*/
186 
187 /* commentata per warning: integer overflow in expression
188 void prand(long *pto, long row, long col,long toggle )
189 
190 {
191 
192 long i;
193 static long vx[98],vy[99];
194 long dnr;
195 if(toggle==0){
196  srand(time(NULL));
197  for(i=0;i<=97;i++) rand();
198  for(i=0;i<=97;i++) vx[i]=(row-2)*(long)rand()/(RAND_MAX+1);
199  for(i=0;i<=97;i++) vy[i]=(col-2)*(long)rand()/(RAND_MAX+1);
200  toggle=1;
201  }
202 
203 dnr=98*(long)rand()/(RAND_MAX+1);
204 if(dnr>97) {printf("dnr was %d",dnr); t_error("Something is wrong here");}
205 pto[0]=vx[dnr]+2;
206 vx[dnr]=(row-2)*(long)rand()/(RAND_MAX+1);
207 dnr=98*(long)rand()/(RAND_MAX+1);
208 if(dnr>97) {printf("dnr was %d",dnr); t_error("Something is wrong here");}
209 pto[1]=vy[dnr]+2;
210 vy[dnr]=(col-2)*(long)rand()/(RAND_MAX+1);
211 
212 }
213 */
214 
215 /*--------------------------------------------------------------------------*/
216 
217 /* commentata per warning: integer overflow in expression
218 long rrand(long row )
219 
220 {
221 long i;
222 static double vx[98];
223 static long toggle;
224 long pto,dnr; */
225 /* char ch; */
226 /* Initializing the randon number generator */
227 /* commentata per warning: integer overflow in expression
228 if(toggle==0){
229  srand(time(NULL));
230  for(i=0;i<=97;i++) rand();
231  for(i=0;i<=97;i++) vx[i]=(double)rand()/(RAND_MAX+1);
232  toggle=1;
233  }
234 
235 dnr=98*(long)rand()/(RAND_MAX+1);
236 if(dnr>97) {printf("dnr was %d",dnr); t_error("merda");}
237 pto=abs(row*vx[dnr]);
238 vx[dnr]=(double)rand()/(RAND_MAX+1);*/
239 /* printf("pto = %d",pto); scanf("%c",&ch);*/
240 /* commentata per warning: integer overflow in expression
241 return pto;
242 
243 }
244 
245 */
246 /*--------------------------------------------------------------------------*/
247 /* generator of random variables with lognormal distribution,
248 average vmed, standard deviation vvar */
249 
250 
251 double uuvel(float vmed,float vvar)
252 {
253 static long idum=3;//ibumb=-2;
254 float umed, uvar,u,v;
255 double y;
256 uvar=log(1+pow((vvar/vmed),2));
257 umed=0.5*log(pow(vmed,2)/(1+pow((vvar/vmed),2)));
258 y=gasdev(&idum);
259 u=umed+y*sqrt(uvar);
260 v=exp(u);
261 return v;
262 }
263 
264 
265 
266 /*--------------------------------------------------------------------------*/
267 double gasdev(long *idum)
268 
269 /* returns a normally distributed random deviate with zero mean mean
270 and unit variance, using rann1(idum) as source of uniform deviates */
271 {
272 
273 /* double ran1(long *idum); */
274 static int iset=0;
275 static double gset;
276 double fac,rsq,v1,v2;
277 
278 if(iset==0){
279 
280 /* we don't have an extra deviate handy, so pick two uniform numbers in
281 the square extending from -1 to +1 in each direction, see if they are in
282 the unit circle and if they are not, try again */
283 
284 do{
285 v1=2.0*ran1(idum)-1.0;
286 v2=2.0*ran1(idum)-1.0;
287 /*printf("v1=%f,v2=%f",(v1+1.0)/2.0,(v2+1.0)/2.0);*/
288 rsq=v1*v1+v2*v2;
289 } while (rsq>=1.0 || rsq==0.0);
290 fac=sqrt(-2.0*log(rsq)/rsq);
291 
292 /* now make the Box-Muller transformation to get two normal deviates.
293 return one and save the other for the next time. */
294 
295 gset=v1*fac;
296 
297 /* set flag */
298 
299 iset=1;
300 return v2*fac;
301 }else{
302 
303 /* we have an extra deviate handy, so unset the flag and return it */
304 
305 iset=0;
306 return gset;
307 }
308 }
309 
310 
311