You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
134 lines
2.8 KiB
Fortran
134 lines
2.8 KiB
Fortran
c FILE : random.f
|
|
c ----------------------------------------------------------------------
|
|
c ------ Uniform Distribution Random Number Generator ------------------
|
|
c ----------------------------------------------------------------------
|
|
c From Numerical Recipies.
|
|
|
|
FUNCTION ran3(idum)
|
|
INTEGER*2 i,ii,k,inext,inextp,iff
|
|
INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
|
|
REAL*4 fac, ran3
|
|
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
|
|
save inext,inextp,ma,iff
|
|
data iff /0/
|
|
C
|
|
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 i=1,54
|
|
ii=mod(21*i,55)
|
|
ma(ii)=mk
|
|
mk=mj-mk
|
|
if(mk.lt.mz) mk=mk+mbig
|
|
mj=ma(ii)
|
|
end do
|
|
do k=1,4
|
|
do i=1,55
|
|
ma(i)=ma(i)-ma(1+mod(i+30,55))
|
|
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
|
|
end do
|
|
end do
|
|
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
|
|
ran3=real(mj)*fac
|
|
return
|
|
end
|
|
|
|
|
|
FUNCTION ran4(idum)
|
|
INTEGER*2 i,ii,k,inext,inextp,iff
|
|
INTEGER*4 idum,mbig,mseed,mz,mj,ma(55),mk
|
|
REAL*4 fac, ran4
|
|
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
|
|
save inext,inextp,ma,iff
|
|
data iff /0/
|
|
C
|
|
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 i=1,54
|
|
ii=mod(21*i,55)
|
|
ma(ii)=mk
|
|
mk=mj-mk
|
|
if(mk.lt.mz) mk=mk+mbig
|
|
mj=ma(ii)
|
|
end do
|
|
do k=1,4
|
|
do i=1,55
|
|
ma(i)=ma(i)-ma(1+mod(i+30,55))
|
|
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
|
|
end do
|
|
end do
|
|
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
|
|
ran4=real(mj)*fac
|
|
return
|
|
end
|
|
|
|
|
|
c ---------------------------------------------------------------------
|
|
c -------- Normally Distributed Random number generator ---------------
|
|
c ---------------------------------------------------------------------
|
|
c From numerical recipies
|
|
c Mean=0, Stdev=1
|
|
|
|
FUNCTION normdis(idum)
|
|
INTEGER iset, idum
|
|
REAL v1, v2, r, fac, normdis, ran4, gset
|
|
SAVE gset, iset
|
|
DATA iset/0/
|
|
c idum=0
|
|
IF (iset.eq.0) THEN
|
|
1 v1=2.*RAN4(idum)-1.
|
|
v2=2.*RAN4(idum)-1.
|
|
r=v1**2+v2**2
|
|
IF (r.ge.1.) GOTO 1
|
|
fac=SQRT(-2.*LOG(r)/r)
|
|
gset=v1*fac
|
|
normdis=v2*fac
|
|
iset=1
|
|
ELSE
|
|
normdis=gset
|
|
iset=0
|
|
END IF
|
|
RETURN
|
|
END
|
|
|
|
c ---------------------------------------------------------------
|
|
|
|
FUNCTION rbit()
|
|
INTEGER rbit, isd
|
|
REAL r, ran3
|
|
r=ran3(isd)
|
|
IF (r.le.0.5) THEN
|
|
rbit=1
|
|
ELSE
|
|
rbit=-1
|
|
ENDIF
|
|
RETURN
|
|
END
|