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
