diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 7d9cec6..156f5e4 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -47,22 +47,23 @@ MODULE moduleRandom END FUNCTION randomIntAB - !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + !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):: x, y + REAL(8):: v1, v2, Rsquare - rnd = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO - CALL RANDOM_NUMBER(y) + 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 - rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + end do + + rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) END FUNCTION randomMaxwellian diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 4fb855a..5e72835 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -289,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian() + v = sqrt(2.0)*self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -302,7 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomHalfMaxwellian() + v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian