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
    FUNCTION randomMaxwellian() RESULT(rnd)
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE

      REAL(8):: rnd
      REAL(8):: x, y

      rnd = 0.D0
      x = 0.D0
      DO WHILE (x == 0.D0)
        CALL RANDOM_NUMBER(x)
      END DO
      CALL RANDOM_NUMBER(y)

      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)

    END FUNCTION randomMaxwellian

    FUNCTION randomHalfMaxwellian() RESULT(rnd)
      USE moduleConstParam, ONLY: PI
      IMPLICIT NONE

      REAL(8):: rnd
      REAL(8):: x, y

      rnd = 0.D0
      x = 0.D0
      DO WHILE (x == 0.D0)
        CALL RANDOM_NUMBER(x)
      END DO
      y = random(-PI, PI)

      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y)

    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

    END FUNCTION randomWeighted

END MODULE moduleRandom