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
|
function xran1(idum)
dimension r(97)
parameter (m1=259200,ia1=7141,ic1=54773,rm1=3.8580247e-6)
parameter (m2=134456,ia2=8121,ic2=28411,rm2=7.4373773e-6)
parameter (m3=243000,ia3=4561,ic3=51349)
data iff /0/
save r, ix1,ix2,ix3
if (idum.lt.0.or.iff.eq.0) then
iff=1
ix1=mod(ic1-idum,m1)
ix1=mod(ia1*ix1+ic1,m1)
ix2=mod(ix1,m2)
ix1=mod(ia1*ix1+ic1,m1)
ix3=mod(ix1,m3)
do 11 j=1,97
ix1=mod(ia1*ix1+ic1,m1)
ix2=mod(ia2*ix2+ic2,m2)
r(j)=(float(ix1)+float(ix2)*rm2)*rm1
11 continue
idum=1
endif
ix1=mod(ia1*ix1+ic1,m1)
ix2=mod(ia2*ix2+ic2,m2)
ix3=mod(ia3*ix3+ic3,m3)
j=1+(97*ix3)/m3
if(j.gt.97.or.j.lt.1)then
write(*,*) 'j is bad in ran1.f',j, 97d0*ix3/m3
STOP
endif
xran1=r(j)
r(j)=(float(ix1)+float(ix2)*rm2)*rm1
return
end
|