JGonzalez
4b040c35c3
After reading some works and reviewing what I had, I've done some corrections to how the randomb velicities in Maxwellian distributions are calculated. These should be correct now.
112 lines
2.4 KiB
Fortran
112 lines
2.4 KiB
Fortran
MODULE moduleRandom
|
|
!Interface for random number generator
|
|
INTERFACE random
|
|
MODULE PROCEDURE randomReal, randomRealAB, randomIntAB
|
|
|
|
END INTERFACE random
|
|
|
|
CONTAINS
|
|
!Returns a Real random number between 0 and 1
|
|
FUNCTION randomReal() RESULT(rnd)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8):: rnd
|
|
|
|
rnd = 0.D0
|
|
CALL RANDOM_NUMBER(rnd)
|
|
|
|
END FUNCTION randomReal
|
|
|
|
!Returns a Real random number between a and b
|
|
FUNCTION randomRealAB(a, b) RESULT(rnd)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: a, b
|
|
REAL(8):: rnd
|
|
REAL(8):: rnd01 !random real between 0 and 1
|
|
|
|
rnd = 0.D0
|
|
CALL RANDOM_NUMBER(rnd01)
|
|
|
|
rnd = (b - a) * rnd01 + a
|
|
|
|
END FUNCTION randomRealAB
|
|
|
|
!Returns an Integer random numnber between a and b
|
|
FUNCTION randomIntAB(a, b) RESULT(rnd)
|
|
IMPLICIT NONE
|
|
|
|
INTEGER, INTENT(in):: a, b
|
|
INTEGER:: rnd
|
|
REAL(8):: rnd01
|
|
|
|
rnd = 0
|
|
CALL RANDOM_NUMBER(rnd01)
|
|
|
|
rnd = a + FLOOR((b+1-a)*rnd01)
|
|
|
|
END FUNCTION randomIntAB
|
|
|
|
!Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
|
|
FUNCTION randomMaxwellian() RESULT(rnd)
|
|
USE moduleConstParam, ONLY: PI
|
|
IMPLICIT NONE
|
|
|
|
REAL(8):: rnd
|
|
REAL(8):: v1, v2, Rsquare
|
|
|
|
Rsquare = 1.D0
|
|
do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0)
|
|
v1 = 2.D0 * random() - 1.D0
|
|
v2 = 2.D0 * random() - 1.D0
|
|
Rsquare = v1**2 + v2**2
|
|
|
|
end do
|
|
|
|
rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare)
|
|
|
|
END FUNCTION randomMaxwellian
|
|
|
|
!Returns a random number in a Maxwellian distribution of mean 0 and width 1
|
|
FUNCTION randomHalfMaxwellian() RESULT(rnd)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8):: rnd
|
|
REAL(8):: x
|
|
|
|
rnd = 0.D0
|
|
x = 0.D0
|
|
DO WHILE (x == 0.D0)
|
|
CALL RANDOM_NUMBER(x)
|
|
END DO
|
|
|
|
rnd = DSQRT(-DLOG(x))
|
|
|
|
END FUNCTION randomHalfMaxwellian
|
|
|
|
!Returns a random number weighted with the cumWeight array
|
|
FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: cumWeight(1:)
|
|
REAL(8), INTENT(in):: sumWeight
|
|
REAL(8):: rnd0b
|
|
INTEGER:: rnd, i
|
|
|
|
rnd0b = random()
|
|
i = 1
|
|
DO
|
|
IF (rnd0b <= cumWeight(i)/sumWeight) THEN
|
|
rnd = i
|
|
EXIT
|
|
|
|
ELSE
|
|
i = i +1
|
|
|
|
END IF
|
|
END DO
|
|
! rnd = MINLOC(DABS(rnd0b - cumWeight), 1)
|
|
|
|
END FUNCTION randomWeighted
|
|
|
|
END MODULE moduleRandom
|