002_Random_utilities.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of CosmicFish.
!
! Copyright (C) 2015-2016 by the CosmicFish authors
!
! The CosmicFish code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file LICENSE at
! the top level of the CosmicFish distribution.
!
!----------------------------------------------------------------------------------------



!----------------------------------------------------------------------------------------


module random_utilities

    use precision

    implicit none

    integer :: rand_inst = 0
    integer, parameter :: krand = kind(1.d0)
    integer :: rand_feedback = 1

contains

    subroutine initrandom(i, i2)
        implicit none
        integer, optional, intent(in) :: i
        integer, optional, intent(in) :: i2
        integer seed_in,kl,ij
        character(len=10) :: fred
        real(krand) :: klr

        if (present(i)) then
            seed_in = i
        else
            seed_in = -1
        end if
        if (seed_in /=-1) then
            if (present(i2)) then
                kl=i2
                if (i2 > 30081) stop 'initRandom:second seed too large'
            else
                kl = 9373
            end if
            ij = i
        else
            call system_clock(count=ij)
            ij = mod(ij + rand_inst*100, 31328)
            call date_and_time(time=fred)
            read (fred,'(e10.3)') klr
            kl = mod(int(klr*1000), 30081)
        end if

        if (rand_feedback > 0 ) write(*,'(" Random seeds:",1I6,",",1I6," rand_inst:",1I4)') ij,kl,rand_inst
        call rmarin(ij,kl)

    end subroutine initrandom

    ! ---------------------------------------------------------------------------------------------
    real(dl) function random_uniform(a,b)

        implicit none

        real(dl), intent(in) :: a
        real(dl), intent(in) :: b


        random_uniform = (b-a)*ranmar() + a

    end function random_uniform

    ! ---------------------------------------------------------------------------------------------
    real(dl) function gaussian_1()

        implicit none

        real(dl) :: R, V1, V2, FAC
        integer , save :: iset = 0
        real(dl), save :: gset

        if ( iset==0 ) then
            r = 2._dl
            do while (r >= 1._dl)
                v1 = 2.d0*ranmar()-1.d0
                v2 = 2.d0*ranmar()-1.d0
                r  = v1**2+v2**2
            end do
            fac        = sqrt(-2.d0*log(r)/r)
            gset       = v1*fac
            gaussian_1 = v2*fac
            iset = 1
        else
            gaussian_1 = gset
            iset = 0
        endif

    end function gaussian_1

    ! ---------------------------------------------------------------------------------------------
    real(dl)  function random_gaussian( mu, sigma )

        implicit none

        real(dl), intent(in) :: mu
        real(dl), intent(in) :: sigma

        random_gaussian = gaussian_1()*sigma + mu

    end function random_gaussian


    ! This random number generator originally appeared in ''Toward a Universal
    ! Random Number Generator'' by George Marsaglia and Arif Zaman.
    ! Florida State University Report: FSU-SCRI-87-50 (1987)
    !
    ! It was later modified by F. James and published in ''A Review of Pseudo-
    ! random Number Generators''
    !
    ! THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
    !    (However, a newly discovered technique can yield
    !        a period of 10^600. But that is still in the development stage.)
    !
    ! It passes ALL of the tests for random number generators and has a period
    !   of 2^144, is completely portable (gives bit identical results on all
    !   machines with at least 24-bit mantissas in the floating point
    !   representation).
    !
    ! The algorithm is a combination of a Fibonacci sequence (with lags of 97
    !   and 33, and operation "subtraction plus one, modulo one") and an
    !   "arithmetic sequence" (using subtraction).
    !
    ! On a Vax 11/780, this random number generator can produce a number in
    !    13 microseconds.
    !========================================================================
    !
    !      PROGRAM TstRAN
    !     INTEGER IJ, KL, I
    ! Thee are the seeds needed to produce the test case results
    !      IJ = 1802
    !      KL = 9373
    !
    !
    ! Do the initialization
    !      call rmarin(ij,kl)
    !
    ! Generate 20000 random numbers
    !      do 10 I = 1, 20000
    !         x = RANMAR()
    !10    continue
    !
    ! If the random number generator is working properly, the next six random
    !    numbers should be:
    !          6533892.0  14220222.0  7275067.0
    !    6172232.0  8354498.0   10633180.0
    !
    !
    !
    !      write(6,20) (4096.0*4096.0*RANMAR(), I=1,6)
    !20    format (3f12.1)
    !      end
    !
    subroutine rmarin(IJ,KL)
        ! This is the initialization routine for the random number generator RANMAR()
        ! NOTE: The seed variables can have values between:    0 <= IJ <= 31328
        !                                                      0 <= KL <= 30081
        !The random number sequences created by these two seeds are of sufficient
        ! length to complete an entire calculation with. For example, if sveral
        ! different groups are working on different parts of the same calculation,
        ! each group could be assigned its own IJ seed. This would leave each group
        ! with 30000 choices for the second seed. That is to say, this random
        ! number generator can create 900 million different subsequences -- with
        ! each subsequence having a length of approximately 10^30.
        !
        ! Use IJ = 1802 & KL = 9373 to test the random number generator. The
        ! subroutine RANMAR should be used to generate 20000 random numbers.
        ! Then display the next six random numbers generated multiplied by 4096*4096
        ! If the random number generator is working properly, the random numbers
        !    should be:
        !           6533892.0  14220222.0  7275067.0
        !           6172232.0  8354498.0   10633180.0
        double precision U(97), C, CD, CM, S, T
        integer I97, J97,i,j,k,l,m
        integer ij,kl
        integer ii,jj


        !      INTEGER IRM(103)

        common /raset1/ u, c, cd, cm, i97, j97
        if( ij < 0  .or.  ij > 31328  .or. &
            kl < 0  .or.  kl > 30081 ) then
            print '(A)', ' The first random number seed must have a value  between 0 and 31328'
            print '(A)',' The second seed must have a value between 0 and   30081'
            stop
        endif
        i = mod(ij/177, 177) + 2
        j = mod(ij    , 177) + 2
        k = mod(kl/169, 178) + 1
        l = mod(kl,     169)
        do ii = 1, 97
            s = 0.0
            t = 0.5
            do jj = 1, 24
                m = mod(mod(i*j, 179)*k, 179)
                i = j
                j = k
                k = m
                l = mod(53*l+1, 169)
                if (mod(l*m, 64) >= 32) then
                    s = s + t
                endif
                t = 0.5 * t
3           continue
            end do
            u(ii) = s
2       continue
        end do
        c = 362436.0 / 16777216.0
        cd = 7654321.0 / 16777216.0
        cm = 16777213.0 /16777216.0
        i97 = 97
        j97 = 33

    end subroutine rmarin

    ! ---------------------------------------------------------------------------------------------
    double precision function ranmar()

        double precision U(97), C, CD, CM
        integer I97, J97
        double precision uni

        common /raset1/ u, c, cd, cm, i97, j97
        !      INTEGER IVEC
        uni = u(i97) - u(j97)
        if( uni < 0.0 ) uni = uni + 1.0
        u(i97) = uni
        i97 = i97 - 1
        if(i97 == 0) i97 = 97
        j97 = j97 - 1
        if(j97 == 0) j97 = 97
        c = c - cd
        if( c < 0.d0 ) c = c + cm
        uni = uni - c
        if( uni < 0.d0 ) uni = uni + 1.0
        ranmar = uni

    end function ranmar

end module random_utilities