MODULE moduleCoulomb
  USE moduleSpecies
  IMPLICIT NONE

  INTEGER:: nCoulombPairs = 0

  !Type for Coulomb iteraction matrix
  TYPE:: interactionsCoulomb
    CLASS(speciesGeneric), POINTER:: sp_i
    CLASS(speciesGeneric), POINTER:: sp_j
    REAL(8):: one_plus_massRatio_ij
    REAL(8):: lnCoulomb !This can be done a function in the future
    REAL(8):: A_i
    REAL(8):: A_j
    REAL(8):: l2_j
    REAL(8):: l2_i
    CONTAINS
      PROCEDURE, PASS:: init => initInteractionCoulomb

  END TYPE interactionsCoulomb

  !Coulomb collision 'matrix'
  TYPE(interactionsCoulomb), ALLOCATABLE:: coulombMatrix(:)

  CONTAINS
    PURE REAL(8) FUNCTION G(x)
      USE moduleConstParam
      IMPLICIT NONE

      REAL(8), INTENT(in):: x

      G = 0.D0
      IF (x /= 0.D0) THEN
        G = (ERF(x) - x*2.D0/SQRT(PI)*EXP(-x**2))/(2.D0*x**2)

      END IF

    END FUNCTION G

    PURE REAL(8) FUNCTION H(x)
      IMPLICIT NONE

      REAL(8), INTENT(in):: x

      H = ERF(x) - G(x)

    END FUNCTION H

    SUBROUTINE initInteractionCoulomb(self, i, j)
      USE moduleSpecies
      USE moduleErrors
      USE moduleConstParam
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(interactionsCoulomb), INTENT(out):: self
      INTEGER, INTENT(in):: i, j
      REAL(8):: Z_i, Z_j
      REAL(8):: scaleFactor

      self%sp_i => species(i)%obj
      self%sp_j => species(j)%obj

      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m

      Z_i = 0.D0
      Z_j = 0.D0
      SELECT TYPE(sp => self%sp_i)
      TYPE IS (speciesCharged)
        Z_i = sp%q

      CLASS DEFAULT
        CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')

      END SELECT

      SELECT TYPE(sp => self%sp_j)
      TYPE IS (speciesCharged)
        Z_j = sp%q

      CLASS DEFAULT
        CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')

      END SELECT

      self%lnCoulomb = 10.0 !Make this function dependent

      scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)

      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent

      self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
      self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent

    END SUBROUTINE initInteractionCoulomb

END MODULE moduleCoulomb