分类:
2009-05-13 07:12:09
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