!injection of particles
MODULE moduleInject

  TYPE:: injectGeneric
    INTEGER:: id
    CHARACTER(:), ALLOCATABLE:: name
    REAL(8):: v !Velocity (module)
    REAL(8):: T(1:3) !Temperature
    REAL(8):: n(1:3) !Direction of injection
    INTEGER:: nParticles !Number of particles to introduce each time step
    INTEGER:: pt !Species of injection
    INTEGER:: nEdges
    INTEGER, ALLOCATABLE:: edges(:) !Array with edges
    REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
    REAL(8):: sumWeight
    CONTAINS
      PROCEDURE, PASS:: init         => initInject
      PROCEDURE, PASS:: addParticles => addParticlesMaxwellian

  END TYPE injectGeneric

  INTEGER:: nInject
  TYPE(injectGeneric), ALLOCATABLE:: inject(:)

  CONTAINS
    !Injection of particles
    SUBROUTINE doInjects()
      IMPLICIT NONE

      INTEGER:: i

      DO i=1, nInject
        CALL inject(i)%addParticles()
      END DO

    END SUBROUTINE doInjects

    SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface)
      USE moduleMesh
      USE moduleMeshCyl
      USE moduleRefParam
      USE moduleConstParam
      USE moduleSpecies
      IMPLICIT NONE

      CLASS(injectGeneric), INTENT(inout):: self
      INTEGER, INTENT(in):: i
      REAL(8), INTENT(in):: v, n(1:3), T(1:3)
      INTEGER, INTENT(in):: pt, physicalSurface
      REAL(8), INTENT(in):: flow
      INTEGER:: e, et
      INTEGER:: phSurface(1:mesh%numEdges)

      self%id         = i
      self%v          = v/v_ref
      self%n          = n
      self%T          = T/T_ref
      self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight)
      self%pt         = pt

      !Gets the edge elements from which particles are injected
      !TODO: Improve this A LOT
      DO e = 1, mesh%numEdges
        phSurface(e) = mesh%edges(e)%obj%physicalSurface

      END DO

      self%nEdges = COUNT(phSurface == physicalSurface)
      ALLOCATE(inject(i)%edges(1:self%nEdges))
      ALLOCATE(inject(i)%weight(1:self%nEdges))
      et = 0
      DO e=1, mesh%numEdges
        IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
          et = et + 1
          self%edges(et)  = mesh%edges(e)%obj%n
          SELECT TYPE(edge => mesh%edges(e)%obj)
          CLASS IS (meshEdgeCyl)
            self%weight(et) = (edge%r(1)+edge%r(2))/2.D0

          END SELECT

        END IF

      END DO
      self%sumWeight = SUM(self%weight)
      nPartInj = nPartInj + self%nParticles

    END SUBROUTINE

    !Random velocity from maxwellian distribution
    REAL(8) FUNCTION vBC(u, vth)
      USE moduleConstParam, ONLY: PI
      REAL(8), INTENT (in):: u, vth
      REAL(8):: x, y
      vBC=0.D0
      x = 0.D0

      DO WHILE (x == 0.D0)
        x=RAND()
      END DO
      y=RAND()

      vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)

    END FUNCTION vBC

    SUBROUTINE addParticlesMaxwellian(self)
      USE moduleSpecies
      USE moduleSolver
      USE moduleMesh
      USE moduleMeshCyl
      IMPLICIT NONE

      CLASS(injectGeneric), INTENT(in):: self
      REAL(8):: randomX
      INTEGER:: j
      INTEGER:: randomEdge
      REAL(8):: randomPos
      REAL(8):: vVec(1:3), vTh(1:3)
      !Edge nodes
      INTEGER:: n1 = 0, n2 = 0
      !Edge nodes coordinates
      REAL(8):: p1(1:3) = 0.D0, p2(1:3) = 0.D0
      INTEGER:: nMin, nMax !Min and Max index in partInj array
      INTEGER:: n

      vVec = self%v * self%n
      vTh  = DSQRT(self%T/species(self%pt)%obj%m)

      !Insert particles
      !TODO: Adjust for multiple injections
      nMin = 1
      nMax = self%nParticles
      !Assign particle type
      partInj(nMin:nMax)%pt = self%pt
      !Assign weight to particle.
      partInj(nMin:nMax)%weight = species(self%pt)%obj%weight
      partInj(nMin:nMax)%n_in = .TRUE.
      DO n = nMin, nMax
        !Select edge randomly from which inject particle
        randomX = RAND()*self%sumWeight
        DO j = 1, self%nEdges
          IF (randomX < self%weight(j)) EXIT
          randomX = randomX - self%weight(j)

        END DO

        randomEdge = self%edges(j)

        !Get coordinates of edge nodes
        SELECT TYPE(edge => mesh%edges(randomEdge)%obj)
        CLASS IS (meshEdgeCyl)
          n1 = edge%n1%n
          n2 = edge%n2%n

        END SELECT
        SELECT TYPE(node => mesh%nodes(n1)%obj)
        TYPE IS (meshNodeCyl)
          p1(1) = node%z
          p1(2) = node%r

        END SELECT
        SELECT TYPE(node => mesh%nodes(n2)%obj)
        TYPE IS (meshNodeCyl)
          p2(1) = node%z
          p2(2) = node%r

        END SELECT

        !Volume associated to the edge:
        IF (DOT_PRODUCT(self%n, mesh%edges(randomEdge)%obj%normal) >= 0.D0) THEN
          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e1%n
        
        ELSE
          partInj(n)%e_p = mesh%edges(randomEdge)%obj%e2%n

        END IF

        partInj(n)%v = (/ vBC(vVec(1), vTh(1)), &
                           vBC(vVec(2), vTh(2)), &
                           vBC(vVec(3), vTh(3)) /)

        !Random position in edge
        !TODO: Use edge coordinates and transformations for this process
        randomPos = RAND()
        partInj(n)%r(1) = p1(1) + randomPos*(p2(1) - p1(1))
        partInj(n)%r(2) = p1(2) + randomPos*(p2(2) - p1(2))
        partInj(n)%r(3) = p1(3) + randomPos*(p2(3) - p1(3))


        !Push new particle
        CALL push(partInj(n))
        !Assign cell to new particle
        CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n))

      END DO

    END SUBROUTINE addParticlesMaxwellian

END MODULE moduleInject