source/unres/src_CSA/ran.f

   1 ccccccccccccccccccccccccccccccccccccccccccccccccc
   2       FUNCTION ran0(idum)
   3       INTEGER idum,IA,IM,IQ,IR,MASK
   4       REAL ran0,AM
   5       PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
   6      *MASK=123459876)
   7       INTEGER k
   8       idum=ieor(idum,MASK)
   9       k=idum/IQ
  10       idum=IA*(idum-k*IQ)-IR*k
  11       if (idum.lt.0) idum=idum+IM
  12       ran0=AM*idum
  13       idum=ieor(idum,MASK)
  14       return
  15       END
  16 C  (C) Copr. 1986-92 Numerical Recipes Software *11915
  17 ccccccccccccccccccccccccccccccccccccccccccccccccc
  18       FUNCTION ran1(idum)
  19       INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
  20       REAL ran1,AM,EPS,RNMX
  21       PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
  22      *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
  23       INTEGER j,k,iv(NTAB),iy
  24       SAVE iv,iy
  25       DATA iv /NTAB*0/, iy /0/
  26       if (idum.le.0.or.iy.eq.0) then
  27         idum=max(-idum,1)
  28         do 11 j=NTAB+8,1,-1
  29           k=idum/IQ
  30           idum=IA*(idum-k*IQ)-IR*k
  31           if (idum.lt.0) idum=idum+IM
  32           if (j.le.NTAB) iv(j)=idum
  33 11      continue
  34         iy=iv(1)
  35       endif
  36       k=idum/IQ
  37       idum=IA*(idum-k*IQ)-IR*k
  38       if (idum.lt.0) idum=idum+IM
  39       j=1+iy/NDIV
  40       iy=iv(j)
  41       iv(j)=idum
  42       ran1=min(AM*iy,RNMX)
  43       return
  44       END
  45 C  (C) Copr. 1986-92 Numerical Recipes Software *11915
  46 ccccccccccccccccccccccccccccccccccccccccccccccccc
  47       FUNCTION ran2(idum)
  48       INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
  49       REAL ran2,AM,EPS,RNMX
  50       PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
  51      *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
  52      *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
  53       INTEGER idum2,j,k,iv(NTAB),iy
  54       SAVE iv,iy,idum2
  55       DATA idum2/123456789/, iv/NTAB*0/, iy/0/
  56       if (idum.le.0) then
  57         idum=max(-idum,1)
  58         idum2=idum
  59         do 11 j=NTAB+8,1,-1
  60           k=idum/IQ1
  61           idum=IA1*(idum-k*IQ1)-k*IR1
  62           if (idum.lt.0) idum=idum+IM1
  63           if (j.le.NTAB) iv(j)=idum
  64 11      continue
  65         iy=iv(1)
  66       endif
  67       k=idum/IQ1
  68       idum=IA1*(idum-k*IQ1)-k*IR1
  69       if (idum.lt.0) idum=idum+IM1
  70       k=idum2/IQ2
  71       idum2=IA2*(idum2-k*IQ2)-k*IR2
  72       if (idum2.lt.0) idum2=idum2+IM2
  73       j=1+iy/NDIV
  74       iy=iv(j)-idum2
  75       iv(j)=idum
  76       if(iy.lt.1)iy=iy+IMM1
  77       ran2=min(AM*iy,RNMX)
  78       return
  79       END
  80 C  (C) Copr. 1986-92 Numerical Recipes Software *11915
  81 ccccccccccccccccccccccccccccccccccccccccccccccccc
  82       FUNCTION ran3(idum)
  83       INTEGER idum
  84       INTEGER MBIG,MSEED,MZ
  85 C     REAL MBIG,MSEED,MZ
  86       REAL ran3,FAC
  87       PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
  88 C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
  89       INTEGER i,iff,ii,inext,inextp,k
  90       INTEGER mj,mk,ma(55)
  91 C     REAL mj,mk,ma(55)
  92       SAVE iff,inext,inextp,ma
  93       DATA iff /0/
  94       if(idum.lt.0.or.iff.eq.0)then
  95         iff=1
  96         mj=MSEED-iabs(idum)
  97         mj=mod(mj,MBIG)
  98         ma(55)=mj
  99         mk=1
 100         do 11 i=1,54
 101           ii=mod(21*i,55)
 102           ma(ii)=mk
 103           mk=mj-mk
 104           if(mk.lt.MZ)mk=mk+MBIG
 105           mj=ma(ii)
 106 11      continue
 107         do 13 k=1,4
 108           do 12 i=1,55
 109             ma(i)=ma(i)-ma(1+mod(i+30,55))
 110             if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
 111 12        continue
 112 13      continue
 113         inext=0
 114         inextp=31
 115         idum=1
 116       endif
 117       inext=inext+1
 118       if(inext.eq.56)inext=1
 119       inextp=inextp+1
 120       if(inextp.eq.56)inextp=1
 121       mj=ma(inext)-ma(inextp)
 122       if(mj.lt.MZ)mj=mj+MBIG
 123       ma(inext)=mj
 124       ran3=mj*FAC
 125       return
 126       END
 127 C  (C) Copr. 1986-92 Numerical Recipes Software *11915
 128 ccccccccccccccccccccccccccccccccccccccccccccccccc