!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Program file name: random_exp.f90                                      !
!                                                                         !
!  © Tao Pang 2006                                                        !
!                                                                         !
!  Last modified: June 28, 2006                                           !
!                                                                         !
!  (1) This F90 program is created for the book, "An Introduction to      !
!      Computational Physics, 2nd Edition," written by Tao Pang and       !
!      published by Cambridge University Press on January 19, 2006.       !
!                                                                         !
!  (2) No warranties, express or implied, are made for this program.      !
!                                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
MODULE CSEED
  INTEGER :: SEED
END MODULE CSEED
!
PROGRAM RANDOM_EXP
!
! An example of generating an array of exponential
! random numbers.
!
  USE CSEED
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=100
  INTEGER :: time, STIME, I, IC
  INTEGER, DIMENSION (9) :: T
  REAL :: RANE
!
! Initiate the seed from the current time
!
  STIME = time(%REF(0))
  CALL gmtime(STIME,T)
  SEED = T(6)+70*(T(5)+12*(T(4)+31*(T(3)+23*(T(2)+59*T(1)))))
  IF (MOD(SEED,2).EQ.0) SEED = SEED-1
!
  IC = 0
  DO I = 1, N
    WRITE(6, *) "The ", I, "th exponential number is: ", RANE()
  END DO
END PROGRAM RANDOM_EXP
!
FUNCTION RANE() RESULT(E)
!
! Exponential random number generator from a uniform random
! number generator.
!
  REAL E, RANF
!
  E = -ALOG(1.0-RANF())
END FUNCTION RANE
!
FUNCTION RANF() RESULT (CR)
!
! Function to generate a uniform random number in [0,1]
! following x(i+1)=a*x(i) mod c with a=7** 5 and
! c=2** 31-1.  Here the seed is a global variable.
!
  USE CSEED
  IMPLICIT NONE
  INTEGER :: H, L, T, A, C, Q, R
  DATA A/16807/, C/2147483647/, Q/127773/, R/2836/
  REAL :: CR
!
  H = SEED/Q
  L = MOD(SEED, Q)
  T = A*L - R*H
  IF (T .GT. 0) THEN
    SEED = T
  ELSE
    SEED = C + T
  END IF
  CR = SEED/FLOAT(C)
END FUNCTION RANF