! Functions that return a random velocity based on the distribution function
module velocityDistribution
  use moduleRandom
  !Generic type for velocity distribution function
  type, abstract:: velDistGeneric
    contains
      !Returns random velocity from distribution function
      procedure(randomVel_interface), deferred, pass:: randomVel

  end type velDistGeneric

  abstract interface
    function randomVel_interface(self) result(v)
      import velDistGeneric

      class(velDistGeneric), intent(in):: self
      real(8):: v

    end function randomVel_interface

  end interface

  !Container for velocity distributions
  type:: velDistCont
    class(velDistGeneric), allocatable:: obj

  end type velDistCont

  !Maxwellian distribution function
  type, extends(velDistGeneric):: velDistMaxwellian
    real(8):: vTh !Thermal Velocity
    contains
      procedure, pass:: randomVel => randomVelMaxwellian
  
  end type velDistMaxwellian

  type, extends(velDistGeneric):: velDistHalfMaxwellian
    real(8):: vTh !Thermal Velocity
    contains
      procedure, pass:: randomVel => randomVelHalfMaxwellian
  
  end type velDistHalfMaxwellian

  !Dirac's delta distribution function
  type, extends(velDistGeneric):: velDistDelta
    contains
      procedure, pass:: randomVel => randomVelDelta

  end type velDistDelta

  contains
    !Random velocity from Maxwellian distribution
    function randomVelMaxwellian(self) result (v)
      use moduleRandom
      implicit none

      class(velDistMaxwellian), intent(in):: self
      real(8):: v
      v = 0.D0

      v = self%vTh*randomMaxwellian()/sqrt(2.d0)

    end function randomVelMaxwellian

    !Random velocity from Half Maxwellian distribution
    function randomVelHalfMaxwellian(self) result (v)
      implicit none

      class(velDistHalfMaxwellian), intent(in):: self
      real(8):: v
      v = 0.D0

      v = self%vTh*randomHalfMaxwellian()

    end function randomVelHalfMaxwellian

    !Random velocity from Dirac's delta distribution
    PURE function randomVelDelta(self) result(v)
      implicit none

      class(velDistDelta), intent(in):: self
      real(8):: v

      v = 0.D0

    end function randomVelDelta

end module velocityDistribution
