1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
|
#if HAVE_CONFIG_H
# include "config.fh"
#endif
c
c This is really ran3
c
function ran1(idum)
c
c Random number generator from Numerical Recipes
c
c implicit real*4(m)
c parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=2.5e-7)
implicit real*8 (a-h,o-z)
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
dimension ma(55)
data iff /0/
save iff, mj, ma, mk, inext, inextp
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do 11 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
11 continue
do 13 k=1,4
do 12 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12 continue
13 continue
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran1=mj*fac
return
end
c
function gasdev(idum)
implicit real*8 (a-h,o-z)
data iset/0/
save iset, gset
if (iset.eq.0) then
1 v1=2.*ran1(idum)-1.
v2=2.*ran1(idum)-1.
r=v1**2+v2**2
if(r.ge.1.)go to 1
fac=sqrt(-2.*log(r)/r)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
endif
return
end
|
