submodule(moduleMesh) elements
  CONTAINS
    !Reset the output of node
    PURE SUBROUTINE resetOutput(self)
      USE moduleSpecies
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshNode), INTENT(inout):: self
      INTEGER:: k

      DO k = 1, nSpecies
        self%output(k)%den     = 0.D0
        self%output(k)%mom     = 0.D0
        self%output(k)%tensorS = 0.D0

      END DO

    END SUBROUTINE resetOutput

    function meshNodePointer_equal_type_type(self, other) result(isEqual)
      implicit none

      class(meshNodePointer), intent(in):: self, other
      logical:: isEqual

      isEqual = self%obj%n == other%obj%n

    end function meshNodePointer_equal_type_type

    function meshNodePointer_equal_type_int(self, other) result(isEqual)
      implicit none

      class(meshNodePointer), intent(in):: self
      integer, intent(in):: other
      logical:: isEqual

      isEqual = self%obj%n == other

    end function meshNodePointer_equal_type_int

    function meshEdgePointer_equal_type_type(self, other) result(isEqual)
      implicit none

      class(meshEdgePointer), intent(in):: self, other
      logical:: isEqual

      isEqual = self%obj%n == other%obj%n

    end function meshEdgePointer_equal_type_type

    function meshEdgePointer_equal_type_int(self, other) result(isEqual)
      implicit none

      class(meshEdgePointer), intent(in):: self
      integer, intent(in):: other
      logical:: isEqual

      isEqual = self%obj%n == other

    end function meshEdgePointer_equal_type_int

    !Constructs the global K matrix
    SUBROUTINE constructGlobalK(self)
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(inout):: self
      INTEGER:: c
      INTEGER, ALLOCATABLE:: nodes(:)
      REAL(8), ALLOCATABLE:: localK(:,:)
      INTEGER:: i, j
      integer:: n, b, ni
      INTEGER:: info
      EXTERNAL:: dgetrf

      DO c = 1, self%numCells
        associate(nNodes => self%cells(c)%obj%nNodes)
          ALLOCATE(nodes(1:nNodes))
          ALLOCATE(localK(1:nNodes, 1:nNodes))
          nodes  = self%cells(c)%obj%getNodes(nNodes)
          localK = self%cells(c)%obj%elemK(nNodes)

          DO i = 1, nNodes
            DO j = 1, nNodes
              self%K(nodes(i), nodes(j)) = self%K(nodes(i), nodes(j)) + localK(i, j)

            END DO

          END DO

          DEALLOCATE(nodes, localK)

        end associate

      END DO

      ! Modify K matrix due to EM boundary conditions
      DO b = 1, nBoundariesEM
        SELECT TYPE(boundary => boundariesEM(b)%obj)
        TYPE IS(boundaryEMDirichlet)
          DO n = 1, boundary%nNodes
            ni = boundary%nodes(n)%obj%n
            self%K(ni, :)  = 0.D0
            self%K(ni, ni) = 1.D0

          END DO

        TYPE IS(boundaryEMDirichletTime)
          DO n = 1, boundary%nNodes
            ni = boundary%nodes(n)%obj%n
            self%K(ni, :)  = 0.D0
            self%K(ni, ni) = 1.D0

          END DO

        END SELECT

      END DO

      !Compute the PLU factorization of K once boundary conditions have been read
      CALL dgetrf(self%numNodes, self%numNodes, self%K, self%numNodes, self%IPIV, info)
      IF (info /= 0) THEN
        CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM')

      END IF

    END SUBROUTINE constructGlobalK

    ! Gather the value of valNodes at position Xi of an edge
    pure function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f)
      implicit none

      class(meshEdge), intent(in):: self
      real(8), intent(in):: Xi(1:3)
      integer, intent(in):: nNodes
      real(8), intent(in):: valNodes(1:nNodes)
      real(8):: f
      real(8):: fPsi(1:nNodes)

      fPsi = self%fPsi(Xi, nNodes)
      f = dot_product(fPsi, valNodes)

    end function gatherF_edge_scalar

    !Gather the value of valNodes (scalar) at position Xi
    PURE FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f)
      IMPLICIT NONE

      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: valNodes(1:nNodes)
      REAL(8):: f
      REAL(8):: fPsi(1:nNodes)

      fPsi = self%fPsi(Xi, nNodes)
      f = DOT_PRODUCT(fPsi, valNodes)

    END FUNCTION gatherF_cell_scalar

    !Gather the value of valNodes (array) at position Xi
    PURE FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f)
      IMPLICIT NONE

      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3)
      REAL(8):: f(1:3)
      REAL(8):: fPsi(1:nNodes)

      fPsi = self%fPsi(Xi, nNodes)
      f = MATMUL(fPsi, valNodes)

    END FUNCTION gatherF_cell_array

    !Gather the spatial derivative of valNodes (scalar) at position Xi
    PURE FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)
      IMPLICIT NONE

      CLASS(meshCell), INTENT(in):: self
      REAL(8), INTENT(in):: Xi(1:3)
      INTEGER, INTENT(in):: nNodes
      REAL(8), INTENT(in):: valNodes(1:nNodes)
      REAL(8):: df(1:3)
      REAL(8):: dPsi(1:3, 1:nNodes)
      REAL(8):: pDer(1:3,1:3)
      REAL(8):: dPsiR(1:3, 1:nNodes)
      REAL(8):: invJ(1:3, 1:3), detJ

      dPsi  = self%dPsi(Xi, nNodes)
      pDer  = self%partialDer(nNodes, dPsi)
      detJ  = self%detJac(pDer)
      invJ  = self%invJac(pDer)
      dPsiR =  MATMUL(invJ, dPsi)/detJ
      df    = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), &
                 DOT_PRODUCT(dPsiR(2,:), valNodes), &
                 DOT_PRODUCT(dPsiR(3,:), valNodes) /)

    END FUNCTION gatherDF_cell_scalar

    !Scatters particle properties into cell nodes
    SUBROUTINE scatter(self, nNodes, part)
      USE moduleMath
      USE moduleSpecies
      USE OMP_LIB
      IMPLICIT NONE

      CLASS(meshCell), INTENT(inout):: self
      INTEGER, INTENT(in):: nNodes
      CLASS(particle), INTENT(in):: part
      REAL(8):: fPsi(1:nNodes)
      INTEGER:: cellNodes(1:nNodes)
      REAL(8):: tensorS(1:3, 1:3)
      INTEGER:: sp
      INTEGER:: i
      CLASS(meshNode), POINTER:: node
      REAL(8):: pFraction !Particle fraction

      cellNodes = self%getNodes(nNodes)
      fPsi = self%fPsi(part%Xi, nNodes)

      tensorS = outerProduct(part%v, part%v)

      sp = part%species%n

      DO i = 1, nNodes
        node => mesh%nodes(cellNodes(i))%obj
        pFraction = fPsi(i)*part%weight
        CALL OMP_SET_LOCK(node%lock)
        node%output(sp)%den          = node%output(sp)%den          + pFraction
        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + pFraction*part%v(:)
        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
        CALL OMP_UNSET_LOCK(node%lock)

      END DO

    END SUBROUTINE scatter

    !Find next cell for particle
    RECURSIVE SUBROUTINE findCell(self, part, oldCell)
      USE moduleSpecies
      USE moduleErrors
      USE OMP_LIB
      IMPLICIT NONE

      CLASS(meshCell), INTENT(inout):: self
      CLASS(particle), INTENT(inout), TARGET:: part
      CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell
      REAL(8):: Xi(1:3)
      CLASS(meshElement), POINTER:: neighbourElement
      INTEGER:: sp

      Xi = self%phy2log(part%r)
      !Checks if particle is inside 'self' cell
      IF (self%inside(Xi)) THEN
        part%cell = self%n
        part%Xi   = Xi
        part%n_in = .TRUE.
        !Assign particle to listPart_in
        IF (listInCells) THEN
          CALL OMP_SET_LOCK(self%lock)
          sp = part%species%n
          CALL self%listPart_in(sp)%add(part)
          self%totalWeight(sp) = self%totalWeight(sp) + part%weight
          CALL OMP_UNSET_LOCK(self%lock)

        END IF

      ELSE
        !If not, searches for a neighbour and repeats the process.
        CALL self%neighbourElement(Xi, neighbourElement)
        !Defines the next step
        SELECT TYPE(neighbourElement)
        CLASS IS(meshCell)
          !Particle moved to new cell, repeat find procedure
          CALL neighbourElement%findCell(part, self)

        CLASS IS (meshEdge)
          !Particle encountered a surface, apply boundary
          CALL neighbourElement%boundariesParticle(part%species%n)%obj%apply(neighbourElement,part)

          !If particle is still inside the domain, call findCell
          IF (part%n_in) THEN
            IF(PRESENT(oldCell)) THEN
              CALL self%findCell(part, oldCell)

            ELSE
              CALL self%findCell(part)

            END IF
          END IF

        CLASS DEFAULT
          WRITE (*, "(A, I6)") "Element = ", self%n
          CALL criticalError("No connectivity found for element", "findCell")

        END SELECT

      END IF

    END SUBROUTINE findCell

    !If Coll and Particle are the same, simply copy the part%cell into part%cellColl
    SUBROUTINE findCellSameMesh(part)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part

      part%cellColl = part%cell

    END SUBROUTINE findCellSameMesh

    !TODO: try to combine this with the findCell for a regular mesh
    !Find the volume in which particle reside in the mesh for collisions
    !No boundary interaction taken into account
    SUBROUTINE findCellCollMesh(part)
      USE moduleSpecies
      IMPLICIT NONE

      TYPE(particle), INTENT(inout):: part
      LOGICAL:: found
      CLASS(meshCell), POINTER:: cell
      REAL(8), DIMENSION(1:3):: Xi
      CLASS(meshElement), POINTER:: neighbourElement
      INTEGER:: sp

      found = .FALSE.

      cell => meshColl%cells(part%cellColl)%obj
      DO WHILE(.NOT. found)
        Xi = cell%phy2log(part%r)
        IF (cell%inside(Xi)) THEN
          part%cellColl = cell%n
          IF (listInCells) THEN
            CALL OMP_SET_LOCK(cell%lock)
            sp = part%species%n
            CALL cell%listPart_in(sp)%add(part)
            cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight
            CALL OMP_UNSET_LOCK(cell%lock)

          END IF
          found = .TRUE.

        ELSE
          CALL cell%neighbourElement(Xi, neighbourElement)
          SELECT TYPE(neighbourElement)
          CLASS IS(meshCell)
            !Try next element
            cell => neighbourElement

          CLASS DEFAULT
            !Should never happend, but just in case, stops loops
            found = .TRUE.

          END SELECT

        END IF

      END DO

    END SUBROUTINE findCellCollMesh

    !Returns index of volume associated to a position (if any)
    !If no voulme is found, returns 0
    !WARNING: This function is slow and should only be used in initialization phase
    FUNCTION findCellBrute(self, r) RESULT(nVol)
      USE moduleSpecies
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(in):: self
      REAL(8), DIMENSION(1:3), INTENT(in):: r
      INTEGER:: nVol
      INTEGER:: e
      REAL(8), DIMENSION(1:3):: Xi

      !Inits RESULT
      nVol = 0

      DO e = 1, self%numCells
        Xi = self%cells(e)%obj%phy2log(r)
        IF(self%cells(e)%obj%inside(Xi)) THEN
          nVol = self%cells(e)%obj%n
          EXIT

        END IF

      END DO

    END FUNCTION findCellBrute

    !Computes collisions in element
    SUBROUTINE doCollisions(self)
      USE moduleCollisions
      USE moduleSpecies
      USE moduleList
      use moduleRefParam
      USE moduleRandom
      USE moduleOutput
      USE moduleMath
      USE moduleCaseParam, ONLY: timeStep
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(inout), TARGET:: self
      INTEGER:: e
      CLASS(meshCell), POINTER:: cell
      INTEGER:: k, i, j
      INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
      REAL(8):: pMax !Maximum probability of collision
      INTEGER:: nColl
      TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
      TYPE(particle), POINTER:: part_i, part_j
      INTEGER:: n, c
      REAL(8):: vRel, rMass, eRel
      REAL(8):: sigmaVrelTotal
      REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
      REAL(8):: rnd_real !Random number for collision
      INTEGER:: rnd_int  !Random number for collision

      IF (MOD(timeStep, everyColl) == 0) THEN
        !Collisions need to be performed in this iteration
        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
        DO e=1, self%numCells

          cell => self%cells(e)%obj

          !TODO: Simplify this, to many sublevels
          !Iterate over the number of pairs
          DO k = 1, nCollPairs
            !Reset tally of collisions
            IF (collOutput) THEN
              cell%tallyColl(k)%tally = 0

            END IF

            IF (interactionMatrix(k)%amount > 0) THEN
              !Select the species for the collision pair
              i = interactionMatrix(k)%sp_i%n
              j = interactionMatrix(k)%sp_j%n

              !Number of particles per species in the collision pair
              nPart_i = cell%listPart_in(i)%amount
              nPart_j = cell%listPart_in(j)%amount

              IF (nPart_i > 0 .AND. nPart_j > 0) THEN
                !Total number of particles for the collision pair
                nPart = nPart_i + nPart_j

                !Resets the number of collisions in the cell
                nColl = 0

                !Probability of collision for pair i-j
                pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume

                !Number of collisions in the cell
                nColl = NINT(REAL(nPart)*pMax*0.5D0)

                !Converts the list of particles to an array for easy access
                IF (nColl > 0) THEN
                  partTemp_i = cell%listPart_in(i)%convert2Array()
                  partTemp_j = cell%listPart_in(j)%convert2Array()

                END IF

                DO n = 1, nColl
                  !Select random particles
                  part_i => NULL()
                  part_j => NULL()
                  rnd_int = random(1, nPart_i)
                  part_i => partTemp_i(rnd_int)%part
                  rnd_int = random(1, nPart_j)
                  part_j => partTemp_j(rnd_int)%part
                  !If they are the same particle, skip
                  !TODO: Maybe try to improve this
                  IF (ASSOCIATED(part_i, part_j)) THEN
                    CYCLE

                  END IF

                  !If particles do not belong to the species, skip collision
                  !This can happen, for example, if particle has been previously ionized or removed
                  !TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision?
                  IF (part_i%species%n /= i .OR. &
                      part_j%species%n /= j) THEN
                      CYCLE

                  END IF
                  !Obtain the cross sections for the different processes
                  !TODO: From here it might be a procedure in interactionMatrix
                  vRel = NORM2(part_i%v-part_j%v)
                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
                  eRel = rMass*vRel**2
                  CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)

                  !Update maximum sigma*v_rel
                  IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN
                    cell%sigmaVrelMax(k) = sigmaVrelTotal

                  END IF

                  ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
                  probabilityColl = 0.0
                  DO c = 1, interactionMatrix(k)%amount
                    probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1))

                  END DO

                  !Selects random number between 0 and 1
                  rnd_real = random()

                  !If the random number is below the total probability of collision, collide particles
                  IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN

                    !Loop over collisions
                    DO c = 1, interactionMatrix(k)%amount
                      IF (rnd_real <= probabilityColl(c)) THEN
                        !$OMP CRITICAL
                        CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
                        !$OMP END CRITICAL

                        !If collisions are gonna be output, count the collision
                        IF (collOutput) THEN
                          cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1

                        END IF

                        !A collision has ocurred, exit the loop
                        EXIT

                      END IF

                    END DO

                  END IF

                  !Deallocate arrays for next collision
                  DEALLOCATE(sigmaVrel, probabilityColl)

                !End loop collisions in cell
                END DO

              END IF

            END IF

          !End loop collision pairs
          END DO

        !End loop volumes
        END DO
        !$OMP END DO

      END IF

    END SUBROUTINE doCollisions

    SUBROUTINE doCoulomb(self)
      USE moduleCoulomb
      USE moduleRandom
      USE moduleOutput
      USE moduleList
      USE moduleMath
      USE moduleRefParam
      USE moduleConstParam
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in), TARGET:: self
      CLASS(meshCell), POINTER:: cell
      TYPE(interactionsCoulomb):: pair
      INTEGER:: e
      INTEGER:: k
      INTEGER:: i, j
      INTEGER:: n
      INTEGER:: p
      TYPE(lNode), POINTER:: partTemp
      INTEGER(8), ALLOCATABLE:: cellNodes(:)
      CLASS(meshNode), POINTER:: node
      TYPE(outputFormat):: output
      REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node
      REAL(8):: density, velocity(1:3), temperature!values at particle position
      REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference
      REAL(8):: l, lW, l2
      REAL(8):: GlW, HlW
      REAL(8):: normC
      REAL(8):: cosThe, sinThe
      REAL(8):: cosPhi, sinPhi
      REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame
      REAL(8):: A, AW
      REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W
      REAL(8):: theta_per !Random angle for perpendicular direction
      REAL(8):: eps = 1.D-12
      REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij
      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij
      REAL(8):: massSum_ij
      REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji
      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji
      REAL(8):: massSum_ji
      REAL(8):: alpha_num, alpha_den, alpha, beta(1:3)


      !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
      DO e = 1, self%numCells
        cell => self%cells(e)%obj
        cellNodes = cell%getNodes(cell%nNodes)

        ALLOCATE(densityNodes(1:cell%nNodes),       &
                 velocityNodes(1:cell%nNodes, 1:3), &
                 temperatureNodes(1:cell%nNodes))

        DO k=1, nCoulombPairs
          pair = coulombMatrix(k)
          i = pair%sp_i%n
          j = pair%sp_j%n

          !Do scattering of particles from species_i due to species j
          !Compute background properties of species_j
          DO n = 1, cell%nNodes
            node => self%nodes(cellNodes(n))%obj
            CALL calculateOutput(node%output(j), output, node%v, pair%sp_j)
            densityNodes(n)      = output%density/n_ref
            velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
            temperatureNodes(n)  = output%temperature/T_ref

          END DO

          ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
          ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
          ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
          deltaV_ij = 0.D0
          p_ij      = 0.D0
          mass_ij   = 0.D0
          !Loop over particles of species_i
          partTemp => cell%listPart_in(i)%head
          p = 1
          DO WHILE(ASSOCIATED(partTemp))
            density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
            velocity    = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
            temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)

            !If cell temperature is too low, skip particle to avoid division by zero
            IF (temperature>eps) THEN
              l2 = pair%l2_j/temperature
              l  = SQRT(l2)

            ELSE
              partTemp => partTemp%next

              CYCLE

            END IF

            A  = pair%A_i*density

            C = partTemp%part%v - velocity
            normC = NORM2(C)

            !C_3 = z; C_1, C2 = x, y (per)
            C_per = NORM2(C(1:2))
            cosPhi = C(1)  / C_per
            sinPhi = C(2)  / C_per
            cosThe = C(3)  / normC
            sinThe = C_per / normC

            !Rotation matrix to go from W to C
            rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
                                        -sinPhi,        cosPhi,    0.D0,    & !Second column
                                  sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
                                  (/ 3, 3 /))

            !W at start is = (0, 0, normC), so normW = normC
            lW = l * normC
            GlW = G(lW)
            HlW = H(lW)
            AW = A / normC

            !Calculate changes in W due to collision process
            deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin
            deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
            deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()

            !Random angle to distribute perpendicular change in velocity
            theta_per = PI2*random()

            !Change W
            W(1) =         deltaW_per_square * COS(theta_per) 
            W(2) =         deltaW_per_square * SIN(theta_per)
            W(3) = normC + deltaW_par + deltaW_par_square

            !Compute changes in velocity for each particle
            deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
            mass_ij(p)       = pair%sp_i%m*partTemp%part%weight
            p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v

            !Move to the next particle in the list
            partTemp => partTemp%next
            p = p + 1

          END DO

          !Do corresponding collisions
          IF (i /= j) THEN
            !Do scattering of particles from species_j due to species i
            !Compute background properties of species_i
            DO n = 1, cell%nNodes
              node => self%nodes(cellNodes(n))%obj
              CALL calculateOutput(node%output(i), output, node%v, pair%sp_i)
              densityNodes(n)      = output%density/n_ref
              velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
              temperatureNodes(n)  = output%temperature/T_ref

            END DO

            ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
            ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
            ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
            deltaV_ji = 0.D0
            p_ji      = 0.D0
            mass_ji   = 0.D0
            !Loop over particles of species_j
            partTemp => cell%listPart_in(j)%head
            p = 1
            DO WHILE(ASSOCIATED(partTemp))
              density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
              velocity    = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
              temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)

              !If cell temperature is too low, skip particle to avoid division by zero
              IF (temperature>eps) THEN
                l2 = pair%l2_i/temperature
                l  = SQRT(l2)

              ELSE
                partTemp => partTemp%next

                CYCLE

              END IF
              A  = pair%A_j*density

              C = partTemp%part%v - velocity
              normC = NORM2(C)

              !C_3 = z; C_1, C2 = x, y (per)
              C_per = NORM2(C(1:2))
              cosPhi = C(1)  / C_per
              sinPhi = C(2)  / C_per
              cosThe = C(3)  / normC
              sinThe = C_per / normC

              !Rotation matrix to go from W to C
              rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
                                          -sinPhi,        cosPhi,    0.D0,    & !Second column
                                    sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
                                    (/ 3, 3 /))

              !W at start is = (0, 0, normC), so normW = normC
              lW = l * normC
              GlW = G(lW)
              HlW = H(lW)
              AW = A / normC

              !Calculate changes in W due to collision process
              deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin
              deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
              deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()

              !Random angle to distribute perpendicular change in velocity
              theta_per = PI2*random()

              !Change W
              W(1) =         deltaW_per_square * COS(theta_per) 
              W(2) =         deltaW_per_square * SIN(theta_per)
              W(3) = normC + deltaW_par + deltaW_par_square

              !Compute changes in velocity for each particle
              deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
              p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v

              !Move to the next particle in the list
              partTemp => partTemp%next
              p = p + 1

            END DO

          END IF

          !Calculate correction
          !Total mass
          massSum_ij = SUM(mass_ij)
          massSum_ji = 0.D0

          !Beta
          beta = 0.D0
          DO p = 1, cell%listPart_in(i)%amount
            beta = beta + mass_ij(p) * deltaV_ij(p,1:3)

          END DO
          
          IF (i /= j) THEN
            massSum_ji = SUM(mass_ji)
            DO p = 1, cell%listPart_in(j)%amount
              beta = beta + mass_ji(p) * deltaV_ji(p,1:3)

            END DO

          END IF

          beta = beta / (massSum_ij + massSum_ji)

          !Alpha
          alpha_num = 0.D0
          alpha_den = 0.D0
          DO p =1, cell%listPart_in(i)%amount
            alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3))
            alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2

          END DO

          IF (i /= j) THEN
            DO p = 1, cell%listPart_in(j)%amount
              alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3))
              alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2

            END DO

          END IF

          alpha = -2.D0*alpha_num / alpha_den

          !Apply correction to particles velocity
          partTemp => cell%listPart_in(i)%head
          p = 1
          DO WHILE(ASSOCIATED(partTemp))
            partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3))
            partTemp => partTemp%next
            p = p + 1

          END DO

          IF (i /= j) THEN
            partTemp => cell%listPart_in(j)%head
            p = 1
            DO WHILE(ASSOCIATED(partTemp))
              partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3))
              partTemp => partTemp%next
              p = p + 1

            END DO

          END IF

          DEALLOCATE(deltaV_ij, p_ij, mass_ij)

          IF (i /= j) THEN
            DEALLOCATE(deltaV_ji, p_ji, mass_ji)

          END IF

        END DO

        DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)

      END DO
      !$OMP END DO

    END SUBROUTINE doCoulomb

end submodule elements