MODULE Random
!------------------------------------------------------------------------------
!  random.f90
!
!    Random Module  (Aug 2008)
!    Returns random numbers from a special distributions of popular algorithms
!
!------------------------------------------------------------------------------


  CONTAINS

!------------------------------------------------------------------------------
! Flat (uniform) distribution
! Returns a uniformly distributed random real number between [min,max]
!
   REAL FUNCTION Random_Flat(min, max)
   REAL, INTENT(IN), OPTIONAL :: min, max
   REAL :: x

     CALL RANDOM_NUMBER(x)

     IF(PRESENT(min) .AND. PRESENT(max)) THEN
        Random_Flat = min + ABS(max-min)*x
     ELSE
        Random_Flat = x
     END IF

   END FUNCTION Random_Flat

!------------------------------------------------------------------------------
! Integer (uniform) distribution
! Returns a uniformly distributed random integer number between [min,max]
!
   INTEGER FUNCTION Random_Integer(min, max)
   INTEGER, INTENT(IN), OPTIONAL :: min, max
   REAL :: x

     CALL RANDOM_NUMBER(x)

     IF(PRESENT(min) .AND. PRESENT(max)) THEN
        Random_Integer = min + INT(max*x)
     ELSE
        Random_Integer = INT(100*x)        
     END IF

   END FUNCTION Random_Integer

!------------------------------------------------------------------------------
! Gaussian Distribution
! Returns a normally distributed deviate with mean and sigma
! The routine uses the Box-Muller transformation of uniform deviates.
!
   REAL FUNCTION Random_Gauss(mean, sigma)
   INTEGER, INTENT(IN), OPTIONAL :: mean, sigma
   REAL :: x, y, z

     DO
        x = 2.0 * Random_Flat() - 1.0
        y = 2.0 * Random_Flat() - 1.0
        z = x*x + y*y
        if( z <= 1.0 ) exit
     END DO

     IF(PRESENT(mean) .AND. PRESENT(sigma)) THEN
        Random_Gauss = mean + sigma*x*sqrt(-2.0*log(z)/z)
     ELSE
        Random_Gauss = x*sqrt(-2.0*log(z)/z)
     END IF

   END FUNCTION Random_Gauss

!------------------------------------------------------------------------------
! Exponential (decay) distribution
! Returns a random number between times t1 and t2 
! according to f(t) = exp (-t/tau)
!
   REAL FUNCTION Random_Exponential(tau, tmin, tmax)
   REAL, INTENT(IN) :: tau
   REAL, INTENT(IN), OPTIONAL :: tmin, tmax
   REAL :: r1, r2

     IF(PRESENT(tmin) .AND. PRESENT(tmax)) THEN
        r1 =  exp(-tmin/tau)
        r2 =  exp(-tmax/tau)
     ELSE
        r1 = 1.0
        r2 = 0.0
     END IF

     Random_Exponential = -tau*log(r2 + Random_Flat() * (r1-r2) )
   
   END FUNCTION Random_Exponential

!------------------------------------------------------------------------------
! Breit-Wigner Distribution
! Returns a random number from a Breit-Wigner distribution 
! for center mean Full Width Half Maximum fwhm
!
   REAL FUNCTION Random_BreitWigner(mean, fwhm)
   REAL, INTENT(IN), OPTIONAL :: mean, fwhm
   REAL x, y, z
   
     DO
        x = 2.0 * Random_Flat() - 1.0
        y = 2.0 * Random_Flat() - 1.0
        z = x*x + y*y
        if( z <= 1.0 ) exit
     END DO

     IF(PRESENT(mean) .AND. PRESENT(fwhm)) THEN
        Random_BreitWigner = mean + 0.5*fwhm*x/y
     ELSE
        Random_BreitWigner = 0.5*x/y
     END IF

   END FUNCTION Random_BreitWigner

END MODULE Random