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

      v1 = 0.d0
      do while (v1 <= 0.d0)
        v1 = random()

      end do
      v2 = random()

      rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2)

    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):: v1

      rnd = 0.D0
      v1 = 0.D0
      do while (v1 == 0.D0)
        v1 = random()
      
      end do

      rnd = sqrt(-2.d0*log(v1))

    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
