!moduleMeshBoundary: Boundary functions for the mesh edges
submodule(moduleMesh) boundaryParticle
  contains
    ! Returns the array index of the boundary based on the name
    module function boundaryParticleName_to_Index(boundaryName) result(bp)
      use moduleErrors
      implicit none

      character(:), allocatable:: boundaryName
      integer:: bp
      integer:: b

      bp = 0
      do b = 1, nBoundariesParticle
        if (boundaryName == boundariesParticle(b)%obj%name) then
          bp = boundariesParticle(b)%obj%n

        end if

      end do

      if (bp == 0) then
        call criticalError('Boundary ' // boundaryName // ' not found', 'boundaryParticleName_to_Index')
      
      end if
    
    end function boundaryParticleName_to_Index

    module SUBROUTINE reflection_apply(self, edge, part)
      USE moduleCaseParam
      USE moduleSpecies
      IMPLICIT NONE

      class(boundaryReflection), intent(inout):: self
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      !rp  = intersection between particle and edge
      !rpp = final position of particle
      !vpp = final velocity of particle
      REAL(8), DIMENSION(1:3):: rp, vpp

      !Reflect particle velocity
      vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
      part%v = vpp

      rp = edge%intersection(part%r)

      part%r = 2.D0*(rp - part%r) + part%r

      !particle is assumed to be inside
      part%n_in = .TRUE.

    END SUBROUTINE reflection_apply

    !Absoption in a surface
    module SUBROUTINE absorption_apply(self, edge, part)
      USE moduleCaseParam
      USE moduleSpecies
      IMPLICIT NONE

      class(boundaryAbsorption), intent(inout):: self
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      REAL(8):: rpp(1:3) !Position of particle projected to the edge
      REAL(8):: d !Distance from particle to edge

      rpp = edge%intersection(part%r)

      d = NORM2(rpp - part%r)

      IF (d >= 0.D0) THEN
        part%weight = part%weight/d

      END IF

      !Assign new position to particle
      part%r = rpp
      !Remove particle from the domain
      part%n_in = .FALSE.

      !Scatter particle in associated volume
      IF (ASSOCIATED(edge%e1)) THEN
        CALL edge%e1%scatter(edge%e1%nNodes, part)

      ELSE
        CALL edge%e2%scatter(edge%e2%nNodes, part)

      END IF

    END SUBROUTINE absorption_apply

    !Transparent boundary condition
    module SUBROUTINE transparent_apply(self, edge, part)
      USE moduleSpecies
      IMPLICIT NONE

      class(boundaryTransparent), intent(inout):: self
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part

      !Removes particle from domain
      part%n_in = .FALSE.

    END SUBROUTINE transparent_apply

    !Symmetry axis. Reflects particles.
    !Although this function should never be called, it is set as a reflective boundary
    !to properly deal with possible particles reaching a corner and selecting this boundary.
    module SUBROUTINE symmetryAxis_apply(self, edge, part)
      USE moduleSpecies
      IMPLICIT NONE

      class(boundaryAxis), intent(inout):: self
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part

      CALL genericReflection(edge, part)

    END SUBROUTINE symmetryAxis_apply

    ! wallTemperature
    ! During the collision, the boundary transfers part of the energy to the incident particle
    ! Init
    module SUBROUTINE wallTemperature_init(boundary, T, c)
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
      REAL(8):: vTh

      vTh = DSQRT(c * T) / v_ref
      allocate(boundaryWallTemperature:: boundary)

      select type(boundary)
      type is(boundaryWallTemperature)
        boundary%vTh = vTh

      end select

    END SUBROUTINE wallTemperature_init

    ! Apply
    module SUBROUTINE wallTemperature_apply(self, edge, part)
      USE moduleSpecies
      USE moduleRandom
      IMPLICIT NONE

      class(boundaryWallTemperature), intent(inout):: self
      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      INTEGER:: i

      !Modifies particle velocity according to wall temperature
      DO i = 1, 3
        part%v(i) = part%v(i) + self%vTh*randomMaxwellian()

      END DO
          
      CALL genericReflection(edge, part)

    END SUBROUTINE wallTemperature_apply

    ! ionization
    ! An electron will pass through the surface and create a series of ion-electron pairs based on a neutral background
    !init
    module SUBROUTINE ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
      use moduleRefParam
      use moduleSpecies
      use moduleCaseParam
      use moduleConstParam
      use moduleErrors, only: criticalError
      use omp_lib, only: omp_init_lock
      implicit none

      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      real(8), intent(in):: mImpact
      real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
      integer, intent(in):: ion
      real(8), intent(in):: effTime
      character(:), allocatable, intent(in):: crossSection
      real(8), intent(in):: eThreshold
      integer, optional, intent(in):: electronSecondary

      ALLOCATE(boundaryIonization:: boundary)

      SELECT TYPE(boundary)
      TYPE IS(boundaryIonization)
        boundary%m0 = m0 / m_ref
        boundary%n0 = n0 * Vol_ref
        boundary%v0 = v0 / v_ref
        boundary%vTh = DSQRT(kb*T0/m0)/v_ref
        boundary%species => species(ion)%obj
        IF (PRESENT(electronSecondary)) THEN
          SELECT TYPE(sp => species(electronSecondary)%obj)
          TYPE IS(speciesCharged)
            boundary%electronSecondary => sp

          CLASS DEFAULT
            CALL criticalError("Species " // sp%name // " chosen for " // &
                               "secondary electron is not a charged species", 'initIonization')

          END SELECT

        ELSE
          boundary%electronSecondary => NULL()

        END IF
        boundary%effectiveTime = effTime / ti_ref
        CALL boundary%crossSection%init(crossSection)
        CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
        boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
        boundary%deltaV = DSQRT(boundary%eThreshold/mImpact)

        boundary%update => ionization_update
        boundary%print  => ionization_print

        boundary%tally = 0
        call omp_init_lock(boundary%tally_lock)

      END SELECT

    END SUBROUTINE ionization_init

    ! Apply
    module SUBROUTINE ionization_apply(self, edge, part)
      USE moduleList
      USE moduleSpecies
      USE moduleMesh
      USE moduleRefParam
      USE moduleRandom
      USE moduleMath
      use omp_lib, only: omp_set_lock, omp_unset_lock
      IMPLICIT NONE

      class(boundaryIonization), intent(inout):: self
      CLASS(meshEdge), INTENT(inout)::        edge
      CLASS(particle), INTENT(inout)::        part
      REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
      INTEGER:: nIonizations !Number of ionizations based on eRel
      REAL(8):: pIonization !Probability of ionization of each event
      INTEGER:: p
      REAL(8):: v0(1:3) !random velocity of neutral
      TYPE(particle), POINTER:: newElectron
      TYPE(particle), POINTER:: newIon

      mRel = reducedMass(self%m0, part%species%m)
      vRel = SUM(DABS(part%v-self%v0))
      eRel = mRel*vRel**2*5.D-1

      !Maximum number of possible ionizations based on relative energy
      nIonizations = FLOOR(eRel/self%eThreshold)

      DO p = 1, nIonizations
        !Get probability of ionization
        pIonization = 1.D0 - DEXP(-self%n0*self%crossSection%get(eRel)*vRel*self%effectiveTime/REAL(nIonizations))

        !If a random number is below the probability of ionization, create new pair of ion-electron
        IF (random() < pIonization) THEN
          !Assign random velocity to the neutral
          v0(1) = self%v0(1) + self%vTh*randomMaxwellian()
          v0(2) = self%v0(2) + self%vTh*randomMaxwellian()
          v0(3) = self%v0(3) + self%vTh*randomMaxwellian()

          !Allocates the new particles
          ALLOCATE(newElectron)
          ALLOCATE(newIon)

          IF (ASSOCIATED(self%electronSecondary)) THEN
            newElectron%species => self%electronSecondary

          ELSE
            newElectron%species => part%species

          END IF
          newIon%species      => self%species

          newElectron%v = v0 + (1.D0 + self%deltaV*v0/NORM2(v0))
          newIon%v      = v0

          newElectron%r = edge%randPos()
          newIon%r      = newElectron%r

          IF (ASSOCIATED(edge%e1)) THEN
            newElectron%cell = edge%e1%n

          ELSEIF (ASSOCIATED(edge%e2)) THEN
            newElectron%cell = edge%e2%n

          END IF
          newIon%cell      = newElectron%cell

          newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r)
          newIon%Xi      = newElectron%Xi

          newElectron%weight = part%weight
          newIon%weight      = newElectron%weight

          newElectron%n_in = .TRUE.
          newIon%n_in      = .TRUE.

          !Add particles to list
          CALL partSurfaces%setLock()
          CALL partSurfaces%add(newElectron)
          CALL partSurfaces%add(newIon)
          CALL partSurfaces%unsetLock()

          !Electron loses energy due to ionization
          eRel = eRel - self%eThreshold
          vRel = 2.D0*DSQRT(eRel)/mRel

          !Reduce number of possible ionizations
          nIonizations = nIonizations - 1

          ! Increase the tally
          call omp_set_lock(self%tally_lock)
          self%tally = self%tally + 1
          call omp_unset_lock(self%tally_lock)

        END IF

      END DO

      !Removes ionizing electron regardless the number of pair created
      part%n_in = .FALSE.

    END SUBROUTINE ionization_apply

    ! Update
    subroutine ionization_update(self)
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self

      select type(self)
      type is(boundaryIonization)
        self%tally = 0

      end select

    end subroutine ionization_update

    ! Print
    subroutine ionization_print(self, fileID)
      use moduleOutput, only: fmtInt
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self
      integer, intent(in):: fileID

      write(fileID, '(A)') self%name
      select type(self)
      type is(boundaryIonization)
        write(fileID, '(A)') 'Number of successful ionization events per iteration'
        write(fileID, '('//fmtInt//')') self%tally

      end select

    end subroutine ionization_print

    ! quasiNeutrality
    ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge.
    ! Init
    module subroutine quasiNeutrality_init(boundary, s_incident)
      implicit none

      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      integer, intent(in):: s_incident

      allocate(boundaryQuasiNeutrality:: boundary)

      select type(boundary)
      type is(boundaryQuasiNeutrality)
        allocate(boundary%edges(0))

        boundary%s_incident = s_incident

        boundary%update => quasiNeutrality_update
        boundary%print  => quasiNeutrality_print

      end select

    end subroutine quasiNeutrality_init

    ! Apply
    module subroutine quasiNeutrality_apply(self, edge, part)
      use moduleRandom
      implicit none

      class(boundaryQuasiNeutrality), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

      if (random() <= self%alpha(edge%n)) then
        call genericReflection(edge, part)

      else
        call genericTransparent(edge, part)

      end if

    end subroutine quasiNeutrality_apply

    ! Update
    subroutine quasiNeutrality_update(self)
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self
      integer:: e, n, s
      integer, allocatable:: nodes(:)
      class(meshEdge), pointer:: edge
      real(8), allocatable:: density_nodes(:)
      class(meshNode), pointer:: node
      real(8):: den_center
      real(8):: density_incident, density_rest
      real(8):: alpha

      select type(self)
      type is(boundaryQuasiNeutrality)
        do e = 1, size(self%edges)
          edge => self%edges(e)%obj

          density_incident = 0.d0
          density_rest     = 0.d0

          nodes = edge%getNodes(edge%nNodes)
          allocate(density_nodes(1:edge%nNodes))
          do s = 1, nSpecies
            do n = 1, edge%nNodes
              node => mesh%nodes(nodes(n))%obj

              density_nodes(n) = node%output(s)%den

            end do

            ! Gather the density at the edge center and multiply by the species charge
            den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)

            if (s == self%s_incident) then
              density_incident = den_center

            else
              density_rest = density_rest + den_center

            end if

          end do

          ! Correction for this time step
          if (density_rest > 1.0d-10) then
            ! If there is a rest population, correct
            alpha = 1.d0 + density_incident/density_rest

          else
            ! If not, alpha is assumed unchaged
            alpha = 0.d0

          end if

          ! Apply correction with a factor of 0.1 to avoid fast changes
          self%alpha(edge%n) = self%alpha(edge%n) + 1.0d-2 * alpha

          ! Limit alpha between 0 and 1
          self%alpha(edge%n) = min(self%alpha(edge%n), 1.d0)
          self%alpha(edge%n) = max(self%alpha(edge%n), 0.d0)

          deallocate(density_nodes)

        end do

      end select

    end subroutine quasiNeutrality_update

    ! Print output
    subroutine quasiNeutrality_print(self, fileID)
      use moduleOutput, only: fmtColInt, fmtColReal
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self
      integer, intent(in):: fileID
      integer:: e

      write(fileID, '(A)') self%name
      select type(self)
      type is(boundaryQuasiNeutrality)
        write(fileID, '(A,",",A)') '"Edge id"', '"alpha"'
        do e = 1, size(self%edges)
          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)

        end do

      end select

    end subroutine quasiNeutrality_print

    ! outflowAdaptive
    ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow
    ! Init
    subroutine outflowAdaptive_init(boundary, s_incident)
      implicit none

      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
      integer, intent(in):: s_incident

      allocate(boundaryOutflowAdaptive:: boundary)
      
      select type(boundary)
      type is(boundaryOutflowAdaptive)
        allocate(boundary%edges(0))

        boundary%s_incident = s_incident

        boundary%update => outflowAdaptive_update
        boundary%print  => outflowAdaptive_print

      end select

    end subroutine outflowAdaptive_init

    ! Apply
    module subroutine outflowAdaptive_apply(self, edge, part)
      use moduleSpecies
      implicit none

      class(boundaryOutflowAdaptive), intent(inout):: self
      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part

      call genericReflection(edge, part)

      part%v(1) = part%v(1) - self%velocity_shift(edge%n)

      if (dot_product(part%v,edge%normal) <= 0.d0) then

        call genericTransparent(edge, part)

      end if

    end subroutine outflowAdaptive_apply

    ! Update
    subroutine outflowAdaptive_update(self)
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self
      integer:: e, n, s
      integer, allocatable:: nodes(:)
      class(meshEdge), pointer:: edge
      real(8), allocatable:: den_nodes(:)
      real(8), allocatable:: mom_nodes(:)
      class(meshNode), pointer:: node
      real(8):: den_center, mom_center, vel_center
      real(8):: vel_incident, vel_rest
      real(8):: alpha
      
      select type(self)
      type is(boundaryOutflowAdaptive)
        do e = 1, size(self%edges)
          edge => self%edges(e)%obj

          vel_incident = 0.d0
          vel_rest     = 0.d0

          nodes = edge%getNodes(edge%nNodes)
          allocate(den_nodes(1:edge%nNodes))
          allocate(mom_nodes(1:edge%nNodes))
          do s = 1, nSpecies
            den_center = 0.0d0
            mom_center = 0.0d0
            vel_center = 0.0d0
            select type(sp => species(s)%obj)
            type is(speciesCharged)
              do n = 1, edge%nNodes
                node => mesh%nodes(nodes(n))%obj

                den_nodes(n) = node%output(s)%den
                mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)

              end do

              den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
              if (den_center > 1.0d-10) then
                mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
                vel_center = mom_center / den_center

              else
                vel_center = 0.0d0
              
              end if

              if (s == self%s_incident) then
                vel_incident = vel_center
              else
                vel_rest = vel_rest + vel_center
              end if

            end select

          end do

          alpha = vel_rest - vel_incident

          ! Apply correction with a factor of 0.1 to avoid fast changes
          self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha

          deallocate(den_nodes)
          deallocate(mom_nodes)

        end do

      end select

    end subroutine outflowAdaptive_update

    ! Print
    subroutine outflowAdaptive_print(self, fileID)
      use moduleOutput, only: fmtColInt, fmtColReal
      implicit none

      class(boundaryParticleGeneric), intent(inout):: self
      integer, intent(in):: fileID
      integer:: e

      write(fileID, '(A)') self%name
      select type(self)
      type is(boundaryOutflowAdaptive)
        write(fileID, '(A,",",A)') '"Edge id"', '"vel shift non-dimensional"'
        do e = 1, size(self%edges)
          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%velocity_shift(self%edges(e)%obj%n)
        end do

      end select

    end subroutine outflowAdaptive_print


    ! Generic boundary conditions for internal use
    ! reflection
    module subroutine genericReflection(edge, part)
      use moduleCaseParam
      use moduleSpecies
      implicit none

      class(meshEdge), intent(inout):: edge
      class(particle), intent(inout):: part
      !rp  = intersection between particle and edge
      !rpp = final position of particle
      !vpp = final velocity of particle
      real(8), dimension(1:3):: rp, vpp

      !Reflect particle velocity
      vpp = part%v - 2.D0*dot_product(part%v, edge%normal)*edge%normal
      part%v = vpp

      rp = edge%intersection(part%r)

      part%r = 2.D0*(rp - part%r) + part%r

      !particle is assumed to be inside
      part%n_in = .TRUE.

    end subroutine genericReflection

    ! transparent
    module subroutine genericTransparent(edge, part)
      use moduleSpecies
      implicit none

      class(meshEdge), intent(inout)::         edge
      class(particle), intent(inout)::         part

      !Removes particle from domain
      part%n_in = .FALSE.

    end subroutine genericTransparent

    ! Updates the boundary condition values when %update procedure is associated
    module subroutine boundariesParticle_update()
      implicit none

      integer:: b

      do b = 1, nBoundariesParticle
        if (associated(boundariesParticle(b)%obj%update)) then
          call boundariesParticle(b)%obj%update

        end if

      end do

    end subroutine boundariesParticle_update

    ! Writes the output into the Step_XXXXX_boundaryParticles file when the %print procedure is associated
    module subroutine boundariesParticle_write()
      use moduleCaseparam, only: timeStep
      use moduleOutput, only: fileID_boundaryParticle, &
                              formatFileName, &
                              informFileCreation,&
                              boundaryParticleOutput,&
                              generateFilePath, &
                              prefix
      implicit none

      integer:: b
      character(:), allocatable:: fileName

      if (boundaryParticleOutput) then
        fileName = formatFileName(prefix, 'boundariesParticle', 'csv', timeStep)
        call informFileCreation(fileName)
        open(fileID_boundaryParticle, file = generateFilePath(fileName))

        do b = 1, nBoundariesParticle
          if (associated(boundariesParticle(b)%obj%print)) then
            call boundariesParticle(b)%obj%print(fileID_boundaryParticle)

          end if

        end do

        close(fileID_boundaryParticle)

      end if
        
    end subroutine boundariesParticle_write

end submodule boundaryParticle
