MODULE moduleCollisions
  USE moduleSpecies
  USE moduleTable

  !Integer for when collisions are computed
  INTEGER:: everyColl

  !Abstract type for collision between two particles
  TYPE, ABSTRACT:: collisionBinary
    REAL(8):: rMass !Reduced mass
    REAL(8):: sMassInv !Summed mass
    TYPE(table1D):: crossSec !cross section of collision
    CONTAINS
      PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide

  END TYPE collisionBinary

  ABSTRACT INTERFACE
    SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j)
      USE moduleSpecies
      IMPORT:: collisionBinary

      CLASS(collisionBinary), INTENT(in):: self
      REAL(8), INTENT(in):: sigmaVrelMax
      REAL(8), INTENT(inout):: sigmaVrelMaxNew
      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j

    END SUBROUTINE

  END INTERFACE

  !Container for binary collisions
  TYPE:: collisionCont
    CLASS(collisionBinary), ALLOCATABLE:: obj

  END TYPE collisionCont

  !Binary elastic interaction
  TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic
    CONTAINS
      PROCEDURE, PASS:: collide => collideBinaryElastic

  END TYPE collisionBinaryElastic

  !Ionization binary interaction
  TYPE, EXTENDS(collisionBinary):: collisionBinaryIonization
    REAL(8):: eThreshold !Minimum energy (non-dimensional units) required for ionization
    REAL(8):: deltaV !Change in velocity due to exchange of eThreshold
    CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
    CONTAINS
      PROCEDURE, PASS:: collide => collideBinaryIonization

  END TYPE collisionBinaryIonization

  TYPE, EXTENDS(collisionBinary):: collisionBinaryRecombination
    REAL(8):: eBinding !binding energy of free electron in recombining ion
    REAL(8):: deltaV !Change in velocity due to energy exchange
    CLASS(speciesCharged), POINTER:: electron !Pointer to species considerer as electrons
    CONTAINS
      PROCEDURE, PASS:: collide => collideBinaryRecombination

  END TYPE collisionBinaryRecombination

  !Resonant charge-exchange
  TYPE, EXTENDS(collisionBinary):: collisionBinaryChargeExchange
    CONTAINS
      PROCEDURE, PASS:: collide => collideBinaryChargeExchange

  END TYPE collisionBinaryChargeExchange

  !Type for interaction matrix
  TYPE:: interactionsBinary
    INTEGER:: amount
    TYPE(collisionCont), ALLOCATABLE:: collisions(:)
    CONTAINS
      PROCEDURE, PASS:: init => initInteractionBinary

  END TYPE interactionsBinary

  !Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
  TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
  !Folder for collision cross section tables
  CHARACTER(:), ALLOCATABLE:: pathCollisions

  CONTAINS
    !Velocity of center of mass of two particles
    PURE FUNCTION velocityCM(m_i, v_i, m_j, v_j) RESULT(vCM)
      IMPLICIT NONE

      REAL(8), INTENT(in):: m_i, m_j
      REAL(8), INTENT(in), DIMENSION(1:3):: v_i, v_j
      REAL(8):: vCM(1:3)

      vCM = (m_i*v_i + m_j*v_j)/(m_i + m_j)

    END FUNCTION velocityCM

    !Random direction for hard sphere collisions
    FUNCTION randomDirectionVHS() RESULT(n)
      USE moduleConstParam
      USE moduleRandom
      IMPLICIT NONE

      REAL(8):: n(1:3)
      REAL(8):: cosXii, sinXii, eps

      cosXii = random(-1.D0, 1.D0)
      sinXii = DSQRT(1.D0 - cosXii**2)
      eps    = random(0.D0, PI2)

      n = (/ cosXii, sinXii*DCOS(eps), sinXii*DSIN(eps) /)

    END FUNCTION randomDirectionVHS

    !Inits the interaction matrix
    SUBROUTINE initInteractionMatrix(interactionMatrix)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:)
      INTEGER:: nInteractions

      nInteractions = (nSpecies*(nSpecies+1))/2
      ALLOCATE(interactionMatrix(1:nInteractions))

    END SUBROUTINE initInteractionMatrix

    !Gets the interaction index from the collision matrix from index i,j
    FUNCTION interactionIndex(i,j) RESULT(k)

      INTEGER:: i, j
      INTEGER:: p
      INTEGER:: k

      k = i + j
      p = (k + ABS(i - j))/2
      k = k + (p*(p-3))/2

    END FUNCTION interactionIndex

    !Inits the binary interaction
    SUBROUTINE initInteractionBinary(self, amount)
      IMPLICIT NONE

      CLASS(interactionsBinary), INTENT(inout):: self
      INTEGER, INTENT(in):: amount

      self%amount = amount

      ALLOCATE(self%collisions(1:self%amount))

    END SUBROUTINE initInteractionBinary

    !ELASTIC COLLISIONS
    !Inits binary elastic collision
    SUBROUTINE initBinaryElastic(collision, crossSectionFilename, mass_i, mass_j)
      USE moduleTable
      USE moduleRefParam
      USE moduleConstParam
      IMPLICIT NONE

      CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
      REAL(8), INTENT(in):: mass_i, mass_j

      ALLOCATE(collisionBinaryElastic:: collision)

      !Reads data from file
      CALL collision%crossSec%init(crossSectionFilename)

      !Convert to no-dimensional units
      CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)

      !Calculates reduced mass
      collision%sMassInv =  1.D0/(mass_i+mass_j)
      collision%rMass    = (mass_i*mass_j)*collision%sMassInv

    END SUBROUTINE initBinaryElastic

    !Binary elastic process
    SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                    part_i, part_j)
      USE moduleSpecies
      USE moduleConstParam
      USE moduleRandom
      USE moduleMath
      IMPLICIT NONE

      CLASS(collisionBinaryElastic), INTENT(in):: self
      REAL(8), INTENT(in):: sigmaVrelMax
      REAL(8), INTENT(inout):: sigmaVrelMaxNew
      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
      REAL(8):: sigmaVrel
      REAL(8):: vRel !relative velocity
      REAL(8):: eRel !relative energy
      REAL(8):: m_i, m_j
      REAL(8), DIMENSION(1:3):: vCM
      REAL(8):: vp(1:3)

      vRel = NORM2(part_i%v-part_j%v)
      eRel = self%rMass*vRel**2
      sigmaVrel = self%crossSec%get(eRel)*vRel
      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
        m_i = part_i%species%m
        m_j = part_j%species%m
        !Applies the collision
        vCM = velocityCM(m_i, part_i%v, m_j, part_j%v)
        vp  = vRel*randomDirectionVHS()

        !Assign velocities to particles
        part_i%v = vCM + m_j*vp*self%sMassInv
        part_j%v = vCM - m_i*vp*self%sMassInv

      END IF

    END SUBROUTINE collideBinaryElastic

    !ELECTRON IMPACT IONIZATION
    !Inits electron impact ionization
    SUBROUTINE initBinaryIonization(collision, crossSectionFilename, energyThreshold, mass_i, mass_j, electron)
      USE moduleTable
      USE moduleRefParam
      USE moduleConstParam
      USE moduleSpecies
      USE moduleErrors
      IMPLICIT NONE

      CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
      REAL(8), INTENT(in):: energyThreshold
      REAL(8), INTENT(in):: mass_i, mass_j
      CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
      INTEGER:: electronIndex

      ALLOCATE(collisionBinaryIonization:: collision)

      !Reads data from file
      CALL collision%crossSec%init(crossSectionFilename)

      !Convert to no-dimensional units
      CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)

      !Calculates reduced mass
      collision%sMassInv =  1.D0/(mass_i+mass_j)
      collision%rMass    = (mass_i*mass_j)*collision%sMassInv

      !Specific parameters for ionization collision
      SELECT TYPE(collision)
      TYPE IS(collisionBinaryIonization)
        !Assign the energy threshold
        !Input energy is in eV. Convert to J with ev2J and then to
        !non-dimensional units.
        collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2)
        electronIndex = speciesName2Index(electron)
        SELECT TYPE(sp => species(electronIndex)%obj)
        TYPE IS(speciesCharged)
          collision%electron => sp

        CLASS DEFAULT
          CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')

        END SELECT

      END SELECT

    END SUBROUTINE initBinaryIonization

    !Binary electron impact ionization process
    SUBROUTINE collideBinaryIonization(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                    part_i, part_j)
      USE moduleSpecies
      USE moduleErrors
      USE moduleList
      USE moduleRandom
      USE moduleMath
      USE OMP_LIB
      IMPLICIT NONE

      CLASS(collisionBinaryIonization), INTENT(in):: self
      REAL(8), INTENT(in):: sigmaVrelMax
      REAL(8), INTENT(inout):: sigmaVrelMaxNew
      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
      TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
      TYPE(particle), POINTER:: newElectron
      REAL(8):: vRel, eRel
      REAL(8):: sigmaVrel
      REAL(8), DIMENSION(1:3):: vp_e, vp_n

      !eRel (in units of [m][L]^2[t]^-2
      vRel = NORM2(part_i%v-part_j%v)
      eRel = self%rMass*vRel**2
      !Relative energy must be higher than threshold
      IF (eRel > self%eThreshold) THEN
        sigmaVrel = self%crossSec%get(eRel)*vRel
        sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
        IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
          !Find which particle is the ionizing electron
          IF    (ASSOCIATED(part_i%species, self%electron)) THEN
            electron => part_i
            neutral  => part_j

          ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
            electron => part_j
            neutral  => part_i

          ELSE
            CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')

          END IF

          !Exchange energy between 
          vp_e = electron%v*(1.D0 - self%deltaV/NORM2(electron%v))
          vp_n = neutral%v* (1.D0 + self%deltaV/NORM2(neutral%v) )
          
          !Changes velocity of impacting electron
          electron%v = vp_e

          !Creates a new electron from ionization
          ALLOCATE(newElectron)
          newElectron%species => electron%species
          newElectron%v       =  vp_n
          newElectron%r       =  neutral%r
          newElectron%xi      = neutral%xi
          newElectron%n_in    = .TRUE.
          newElectron%vol     = neutral%vol
          newElectron%volColl = neutral%volColl
          newElectron%weight  = neutral%weight
          newElectron%qm      = electron%qm

          !Ionize neutral particle
          SELECT TYPE(sp => neutral%species)
          TYPE IS(speciesNeutral)
            CALL sp%ionize(neutral)

          CLASS DEFAULT
            CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')

          END SELECT

          !Adds new electron to list of new particles from collisions
          CALL OMP_SET_LOCK(lockCollisions)
          CALL partCollisions%add(newElectron)
          CALL OMP_UNSET_LOCK(lockCollisions)

        END IF

      END IF

    END SUBROUTINE collideBinaryIonization

    !ELECTRON ION RESONANT RECOMBINATION
    !Inits electron ion recombination
    SUBROUTINE initBinaryRecombination(collision, crossSectionFilename, energyBinding, mass_i, mass_j, electron)
      USE moduleTable
      USE moduleRefParam
      USE moduleConstParam
      USE moduleSpecies
      USE moduleErrors
      IMPLICIT NONE

      CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
      REAL(8), INTENT(in):: energyBinding
      REAL(8), INTENT(in):: mass_i, mass_j
      CHARACTER(:), ALLOCATABLE, INTENT(in):: electron
      INTEGER:: electronIndex

      ALLOCATE(collisionBinaryRecombination:: collision)

      !Reads data from file
      CALL collision%crossSec%init(crossSectionFilename)

      !Convert to no-dimensional units
      CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)

      !Calculates reduced mass
      collision%sMassInv =  1.D0/(mass_i+mass_j)
      collision%rMass    = (mass_i*mass_j)*collision%sMassInv

      !Specific parameters for ionization collision
      SELECT TYPE(collision)
      TYPE IS(collisionBinaryRecombination)
        !Assign the energy threshold
        !Input energy is in eV. Convert to J with ev2J and then to
        !non-dimensional units.
        collision%eBinding = energyBinding*eV2J/(m_ref*v_ref**2)
        electronIndex = speciesName2Index(electron)
        SELECT TYPE(sp => species(electronIndex)%obj)
        TYPE IS(speciesCharged)
          collision%electron => sp

        CLASS DEFAULT
          CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization')

        END SELECT

      END SELECT

    END SUBROUTINE initBinaryRecombination

    !Binary electron impact ionization process
    SUBROUTINE collideBinaryRecombination(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                          part_i, part_j)
      USE moduleSpecies
      USE moduleErrors
      USE moduleList
      USE moduleRandom
      USE moduleMath
      IMPLICIT NONE

      CLASS(collisionBinaryRecombination), INTENT(in):: self
      REAL(8), INTENT(in):: sigmaVrelMax
      REAL(8), INTENT(inout):: sigmaVrelMaxNew
      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
      TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
      REAL(8):: vRel, eRel
      REAL(8):: sigmaVrel
      REAL(8), DIMENSION(1:3):: vp_i

      !eRel (in units of [m][L]^2[t]^-2
      vRel = NORM2(part_i%v-part_j%v)
      eRel = self%rMass*vRel**2
      !Relative energy must be higher than threshold
      sigmaVrel = self%crossSec%get(eRel)*vRel
      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
        !Find which particle is the ionizing electron
        IF    (ASSOCIATED(part_i%species, self%electron)) THEN
          electron => part_i
          ion      => part_j

        ELSEIF(ASSOCIATED(part_j%species, self%electron)) THEN
          electron => part_j
          ion      => part_i

        ELSE
          CALL criticalError("No matching between input particles and ionizing species", 'collideBinaryIonization')

        END IF

        !Excess energy
        !TODO: This energy should be transformed into photons
        vp_i = ion%v*     (1.D0 -  (vRel + self%deltaV)/NORM2(ion%v))

        !Remove electron from simulation
        electron%n_in = .FALSE.

        !Neutralize ion particle
        SELECT TYPE(sp => ion%species)
        TYPE IS(speciesCharged)
          CALL sp%neutralize(ion)

        CLASS DEFAULT
          CALL criticalError(sp%name // " is not a charge", 'collideBinaryRecombination')

        END SELECT

      END IF

    END SUBROUTINE collideBinaryRecombination

    !RESONANT CHARGE EXCHANGE
    !Inits resonant charge exchange
    SUBROUTINE initBinaryChargeExchange(collision, crossSectionFilename, mass_i, mass_j)
      USE moduleTable
      USE moduleRefParam
      USE moduleConstParam
      IMPLICIT NONE

      CLASS(collisionBinary), INTENT(out), ALLOCATABLE:: collision
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
      REAL(8), INTENT(in):: mass_i, mass_j

      ALLOCATE(collisionBinaryChargeExchange:: collision)

      !Reads data from file
      CALL collision%crossSec%init(crossSectionFilename)

      !Convert to no-dimensional units
      CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)

      !Calculates reduced mass
      collision%sMassInv =  1.D0/(mass_i+mass_j)
      collision%rMass    = (mass_i*mass_j)/collision%sMassInv

    END SUBROUTINE initBinaryChargeExchange

    SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, &
                                           part_i, part_j)
      USE moduleSpecies
      USE moduleRandom
      USE moduleMath
      IMPLICIT NONE

      CLASS(collisionBinaryChargeExchange), INTENT(in):: self
      REAL(8), INTENT(in):: sigmaVrelMax
      REAL(8), INTENT(inout):: sigmaVrelMaxNew
      TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
      REAL(8):: sigmaVrel
      REAL(8):: vRel !relative velocity
      REAL(8):: eRel !relative energy

      !eRel (in units of [m][L]^2[t]^-2
      vRel = NORM2(part_i%v-part_j%v)
      eRel = self%rMass*vRel**2
      sigmaVrel = self%crossSec%get(eRel)*vRel
      sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
      IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN
        SELECT TYPE(sp => part_i%species)
        TYPE IS (speciesNeutral)
          !Species i is neutral, ionize particle i
          CALL sp%ionize(part_i)

        TYPE IS (speciesCharged)
          !Species i is charged, neutralize particle i
          CALL sp%neutralize(part_i)

        END SELECT

        SELECT TYPE(sp => part_j%species)
        TYPE IS (speciesNeutral)
          !Species j is neutral, ionize particle j
          CALL sp%ionize(part_j)

        TYPE IS (speciesCharged)
          !Species j is charged, neutralize particle j
          CALL sp%neutralize(part_j)

        END SELECT

      END IF

    END SUBROUTINE collideBinaryChargeExchange

END MODULE moduleCollisions