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
|
C Taken from Press et al., Numerical Recipes
FUNCTION RAN1(IDUM)
INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
REAL RAN1,AM,EPS,RNMX
PARAMETER(IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
+NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2E-7,RNMX=1.-EPS)
INTEGER J,K,IV(NTAB),IY
SAVE IV,IY
DATA IV/NTAB*0/,IY/0/
IF (IDUM.LE.0.OR.IY.EQ.0) THEN
IDUM=MAX(-IDUM,1)
DO J=NTAB+8,1,-1
K=IDUM/IQ
IDUM=IA*(IDUM-K*IQ)-IR*K
IF (IDUM.LT.0) IDUM=IDUM+IM
IF (J.LE.NTAB) IV(J)=IDUM
ENDDO
IY=IV(1)
ENDIF
K=IDUM/IQ
IDUM=IA*(IDUM-K*IQ)-IR*K
IF (IDUM.LT.0) IDUM=IDUM+IM
J=1+IY/NDIV
IY=IV(J)
IV(J)=IDUM
RAN1=MIN(AM*IY,RNMX)
RETURN
END
FUNCTION GASDEV(IDUM)
INTEGER IDUM,ISET
REAL GASDEV
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.0 .OR. R.EQ.0.0) 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
|
