!****** V E L O C I T Y   D I S T R I B U T I O N *************************************************
!
!      F. Nattino, 05/2011
!
!     Select incidence energy according to flux weighted velocity distribution as 
!      in Michelsen and al. J Chem. Phys. 94, 7502 (1991) 
!    
!     G(E)dE = ( 1 / N ) * E * exp[ - 4Es * ( SQRT(E) - SQRT(Es) ) ** 2 / DelEs**2 ] dE
!
!     DelEs / Es = 2 DelVs / Vs = 2 / S
!
!     N = (DelEs / 2 S)**2 * { ( 1 + S**2 ) * exp( -S**2) + 
!         SQRT(PI) * S * ( 1.5 + S**2 )[ 1 + ERF(S) ] }
!
!     Units: Mass = kg, Vs and DelVs in m/sec, Etransmax in eV
!
!     USE: Constants (J2eV and ev2kjmol) and Random ( function RAN0)

      MODULE VelocityDistribution

       USE Constants
       USE Random

      IMPLICIT NONE

      REAL, SAVE              :: Gmax, Es, SqEs, DelEs, FacNor 
      INTEGER, PARAMETER      :: StepIntegration = 100000
      CHARACTER(*), PARAMETER :: VelocityDistributionFile = 'VelocityDistribution.dat'

        CONTAINS  

!     ---------------------------------------------------------------------------------------------

      SUBROUTINE SetupVelocityDistribution( Mass, Vs, DelVs, Etransmax)

      IMPLICIT NONE

      REAL, INTENT(IN) :: Vs, DelVs, Mass, Etransmax
      REAL             :: srma, EGmax, Expon, Expfac
      REAL             :: Eaver, Dnorm, E, DelE, GE
      INTEGER          :: i

      PRINT'(1/,28("*"),A,28("*")1/)','  V E L O C I T Y  D I S T R I B U T I O N  '

      Es = 0.5 * Mass * Vs**2. * J2eV
      DelEs = 2. * Es * DelVs / Vs
      SqEs = SQRT( Es )

      PRINT'(A,f7.4)', ' Energy corresponding to Stream Velocity:   E0 (eV) = ',Es
      PRINT'(A,f7.4)', ' DeltaE (eV) = ', DelEs
      PRINT'(A,f7.4)', ' Mass (amu) = ', Mass / amu2kg 

      srma = Vs / DelVs
      ! Now normalizing factor 
      FacNor = ( DelEs /( 2. * srma) )**2. * (( 1. + srma**2. )*EXP( -1. * srma**2. ) + &
                SQRT( PI ) * srma * ( 1.5 + srma**2. ) * ( 1. + ERF( srma ) ) )
      EGmax = (( DelEs**2. ) + 2. * ( Es**2. ) + 2. * Es * SQRT( ( Es**2. ) +  &
             ( DelEs ** 2.)))/(4. * Es)
      Expon = 4.* Es * ( sqrt(EGmax) - SqEs )**2. / DelEs**2.
      Expfac = EXP( -1.* Expon)

      PRINT '(A,f7.4)', ' Energy at Maximum Distribution (eV) = ',EGmax 
      Gmax = EGmax / FacNor * Expfac

      ! Integrate the distribution to find average collision energy 
      OPEN(1,FILE=VelocityDistributionFile)
      WRITE(1,*) '# Normalized Translational Energy Distribution. E in eV.'
 
      DelE = Etransmax / REAL( StepIntegration )

      DO i = 1, StepIntegration
        E = REAL(i) * DelE
        Expon = 4. * Es * (SQRT( E ) - SqEs)**2. / DelEs **2.
        Expfac = EXP ( - Expon )
        GE = Expfac * E
        Eaver = Eaver + GE * E * DelE
        Dnorm = Dnorm + GE * DelE
        WRITE(1,*) E, GE / FacNor
      ENDDO

      CLOSE(1)

      PRINT '(1/,3A,1/)',' Velocity Distribution printed in ', VelocityDistributionFile,' !'

      Eaver = Eaver / FacNor
      Dnorm = Dnorm / FacNor

      PRINT'(A,f8.4)',' The average collision energy is (eV) = '    ,  Eaver
      PRINT'(A,f8.4)',' The average collision energy is (kJ/mol) = ' ,  Eaver * ev2kjmol
      PRINT'(A,f9.7)',' Verify normalization of the distribution: ',  DNorm
      
      IF( ABS( DNorm - 1. ) > 1.E-4 ) PRINT*,' WARNING! Increase Etransmax!!'

      RETURN
      END SUBROUTINE

!     ---------------------------------------------------------------------------------------------
      SUBROUTINE SelectEnergyFromDistribution( Etransmax, Etrans0) 

      IMPLICIT NONE

      REAL, INTENT(IN)    :: Etransmax 
      REAL, INTENT(OUT)   :: Etrans0
      REAL                :: Expon, Expfac, Dice, G, Gran

2002  CONTINUE

      Dice = RAN0()

      Etrans0 = Dice * Etransmax
      Expon = 4.* Es * ( SQRT( Etrans0 ) - SqEs )**2. / DelEs**2.
      Expfac = EXP( -1.* Expon)
      G = Etrans0 / FacNor * Expfac
      Dice = RAN0()
      Gran = Dice * Gmax

      IF ( Gran > G ) then
       GOTO 2002
      ENDIF

      RETURN
      END SUBROUTINE

      END MODULE