randoms.f90

Go to the documentation of this file.

! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16  Time: 11:09:33



SUBROUTINE gbrshi(n,a)
    IMPLICIT NONE
    INTEGER :: i
    INTEGER :: ian
    INTEGER :: iboost
    INTEGER :: ic
    INTEGER :: idum
    INTEGER :: irotor
    INTEGER :: iseed
    INTEGER :: it
    INTEGER :: iwarm
    INTEGER :: j
    INTEGER :: jseed
    INTEGER :: jwarm
    INTEGER :: k
    INTEGER :: m
    INTEGER :: mbuff

    INTEGER, INTENT(IN)                      :: n
    REAL, INTENT(OUT)                        :: a(*)
    INTEGER, PARAMETER :: nb=511
    INTEGER, PARAMETER :: ia=16807
    INTEGER, PARAMETER :: im=2147483647
    INTEGER, PARAMETER :: iq=127773
    INTEGER, PARAMETER :: ir=2836
    REAL, PARAMETER :: aeps=1.0e-10
    REAL, PARAMETER :: scalin=4.6566125e-10
    common/ranbuf/mbuff(0:nb),ian,ic,iboost

    INTEGER :: istart

    irotor(m,n)=ieor(ior(ishft(m,17),ishft(m,-15)),n)
    DATA istart/0/,iwarm/10/,iseed/4711/
    IF(istart /= 0) go to 20
    WRITE(*,*) ' Automatic GBRSHI initialization using:'
    !     initialize buffer
10  idum=iseed+9876543          ! prevent damage, if iseed=0
    WRITE(*,*) '           ISEED=',iseed,'   IWARM=',iwarm
    DO j=0,nb+1                 ! fill buffer
        k=idum/iq                  ! minimal standard generator
        idum=ia*(idum-k*iq)-ir*k   !    with Schrages method
        IF(idum < 0) idum=idum+im !
        mbuff(j)=ishft(idum,1)     ! fill in leading bit
    END DO
    ian=iand(ian,nb)            ! mask angle
    ic=1                        ! set pointer
    iboost=0
    DO j=1,iwarm*nb             ! warm up a few times
        it=mbuff(ian)              ! hit ball angle
        mbuff(ian)=irotor(it,ic)   ! new spin
        ic=it                      ! replace red spin
        ian=iand(it+iboost,nb)     ! boost and mask angle
        iboost=iboost+1            ! increment boost
    END DO
    IF(istart < 0) return      ! return for RBNVIN
    istart=1                    ! set done-flag
    !     generate array of r.n.
    20   DO i=1,n
        it=mbuff(ian)              ! hit ball angle
        mbuff(ian)=irotor(it,ic)   ! new spin
        ic=it                      ! replace red spin
        ian=iand(it+iboost,nb)     ! boost and mask angle
        a(i)=float(ishft(it,-1))*scalin+aeps ! avoid zero output
        iboost=iboost+1            ! increment boost
    END DO
    iboost=iand(iboost,nb)
    return

    entry gbrvin(jseed,jwarm)   ! initialize, but only once
    IF(istart == 0) THEN
        WRITE(*,*) ' Gbrshi initialization by GBRVIN-call using:'
        iseed=jseed              ! copy seed and
        iwarm=jwarm              ! warm-up parameter
        istart=-1                ! start flag
        go to 10
    END IF
END SUBROUTINE gbrshi

SUBROUTINE gbrtim
    IMPLICIT NONE
    INTEGER :: jseed
    REAL :: time

    LOGICAL :: done
    DATA    done/.false./
    IF(done) return
    jseed=time()
    WRITE(*,*) ' Gbrshi initialialization using Time()'
    CALL gbrvin(jseed,10)
    done=.true.
END SUBROUTINE gbrtim


REAL FUNCTION uran()     ! U(0,1)
    IMPLICIT NONE
    INTEGER :: indx
    INTEGER :: ndim

    parameter(ndim=100)
    REAL :: buffer(ndim)
    DATA indx/ndim/
    SAVE indx,buffer
    indx=mod(indx,ndim)+1
    IF(indx == 1) CALL gbrshi(ndim,buffer)
    uran=buffer(indx)
END FUNCTION uran


REAL FUNCTION gran()     ! N(0,1)
    IMPLICIT NONE
    REAL :: al
    REAL :: cs
    INTEGER :: indx
    INTEGER :: kn
    INTEGER :: ndim
    REAL :: radsq
    REAL :: rn1
    REAL :: rn2
    REAL :: sn

    parameter(ndim=100)
    REAL:: buffer(ndim)
    DATA indx/ndim/,kn/1/
    SAVE indx,buffer,kn,cs,al
    !     ...
    IF(kn <= 1) THEN
        !        two U(-1,+1) random numbers
10      indx=mod(indx,ndim)+2
        IF(indx == 2) CALL gbrshi(ndim,buffer)
        rn1=buffer(indx-1)-1.0+buffer(indx-1)
        rn2=buffer(indx  )-1.0+buffer(indx)
        radsq=rn1*rn1+rn2*rn2
        IF(radsq > 1.0) go to 10 ! test point inside circle?
        !        sine and cosine for random phi
        sn=rn1/sqrt(radsq)
        cs=rn2/sqrt(radsq)
        !        transform to gaussians
        al=sqrt(-2.0*alog(radsq))
        kn =2
        gran=sn*al
    ELSE
        kn =1
        gran=cs*al
    END IF
END FUNCTION gran