!moduleMeshBoundary: Boundary functions for the mesh edges
submodule(moduleMesh) boundary
  contains
    module subroutine initBoundary(self, config, object)
      use json_module
      use moduleRefParam, only: m_ref
      use moduleConstParam, only: me
      implicit none

      class(boundaryGeneric), allocatable, intent(out):: self
      type(json_file), intent(inout):: config
      character(:), allocatable, intent(in):: object
      character(:), allocatable:: bType
      logical:: found
      real(8):: Tw, cw !Wall temperature and specific heat
      !neutral Properties
      real(8):: m0, n0, T0
      real(8), dimension(:), allocatable:: v0
      real(8):: effTime
      real(8):: eThreshold !Energy threshold
      integer:: speciesID, electronSecondaryID
      character(:), allocatable:: speciesName, crossSection, electronSecondary

      CALL config%get(object // '.name', self%name, found)
      CALL config%get(object // '.type', bType,     found)
      SELECT CASE(bType)
      CASE('reflection')
        ALLOCATE(boundaryReflection:: self)

      CASE('absorption')
        ALLOCATE(boundaryAbsorption:: self)

      CASE('transparent')
        ALLOCATE(boundaryTransparent:: self)

      CASE('axis')
        ALLOCATE(boundaryAxis:: self)

      CASE('wallTemperature')
        CALL config%get(object // '.temperature', Tw, found)
        IF (.NOT. found) CALL criticalError("temperature not found for wallTemperature boundary type", 'readBoundary')
        CALL config%get(object // '.specificHeat', cw, found)
        IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')

        CALL initWallTemperature(self, Tw, cw)

      CASE('ionization')
        !Neutral parameters
        CALL config%get(object // '.neutral.ion', speciesName, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'ion' for neutrals in ionization", 'readBoundary')
        speciesID = speciesName2Index(speciesName)
        CALL config%get(object // '.neutral.mass', m0, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'mass' for neutrals in ionization", 'readBoundary')
        CALL config%get(object // '.neutral.density', n0, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'density' for neutrals in ionization", 'readBoundary')
        CALL config%get(object // '.neutral.velocity', v0, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'velocity' for neutrals in ionization", 'readBoundary')
        CALL config%get(object // '.neutral.temperature', T0, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'temperature' for neutrals in ionization", 'readBoundary')

        CALL config%get(object // '.effectiveTime', effTime, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'effectiveTime' for ionization", 'readBoundary')

        CALL config%get(object // '.energyThreshold', eThreshold, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'eThreshold' in ionization", 'readBoundary')

        CALL config%get(object // '.crossSection', crossSection, found)
        IF (.NOT. found) CALL criticalError("missing parameter 'crossSection' for neutrals in ionization", 'readBoundary')

        CALL config%get(object // '.electronSecondary', electronSecondary, found)
        electronSecondaryID = speciesName2Index(electronSecondary)
        IF (found) THEN
          CALL initIonization(self, me/m_ref, m0, n0, v0, T0, &
                                     speciesID, effTime, crossSection, eThreshold,electronSecondaryID)


        ELSE
          CALL initIonization(self, me/m_ref, m0, n0, v0, T0, &
                                     speciesID, effTime, crossSection, eThreshold)

        END IF

      case('quasiNeutrality')
        call initQuasiNeutrality(self)

      CASE DEFAULT
        CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')

      END SELECT

    end subroutine initBoundary

    SUBROUTINE initWallTemperature(boundary, T, c)
      USE moduleRefParam
      IMPLICIT NONE

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

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

    END SUBROUTINE initWallTemperature

    SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
      USE moduleRefParam
      USE moduleSpecies
      USE moduleCaseParam
      USE moduleConstParam
      USE moduleErrors
      IMPLICIT NONE

      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
      real(8), intent(in):: mImpact
      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
      INTEGER, INTENT(in):: ion
      INTEGER, OPTIONAL, INTENT(in):: electronSecondary
      REAL(8):: effTime
      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
      REAL(8), INTENT(in):: eThreshold

      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)

      END SELECT

    END SUBROUTINE initIonization

    subroutine initQuasiNeutrality(boundary)
      implicit none

      class(boundaryGeneric), allocatable, intent(out):: boundary
      integer:: e, et

      allocate(boundaryQuasiNeutrality:: boundary)

      select type(boundary)
      type is(boundaryQuasiNeutrality)
        boundary%alpha = 0.d0

      end select

    end subroutine initQuasiNeutrality

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

      class(boundaryReflection), intent(in):: 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

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

      class(boundaryAbsorption), intent(in):: 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

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

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

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

    END SUBROUTINE transparent

    !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.
    SUBROUTINE symmetryAxis(self, edge, part)
      USE moduleSpecies
      IMPLICIT NONE

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

      CALL genericReflection(edge, part)

    END SUBROUTINE symmetryAxis

    !Wall with temperature
    SUBROUTINE wallTemperature(self, edge, part)
      USE moduleSpecies
      USE moduleRandom
      IMPLICIT NONE

      class(boundaryWallTemperature), intent(in):: 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

    !Ionization surface: an electron will pass through the surface
    ! and create an ion-electron pair based on a neutral background
    SUBROUTINE ionization(self, edge, part)
      USE moduleList
      USE moduleSpecies
      USE moduleMesh
      USE moduleRefParam
      USE moduleRandom
      USE moduleMath
      IMPLICIT NONE

      class(boundaryIonization), intent(in):: 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

        END IF

      END DO

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

    END SUBROUTINE ionization

    subroutine quasiNeutrality(self, edge, part)
      use moduleRandom
      implicit none

      class(boundaryQuasiNeutrality), intent(in):: self
      class(meshEdge), intent(inout)::             edge
      class(particle), intent(inout)::             part
      real(8), allocatable:: density(:)
      class(meshCell), pointer:: cell
      real(8):: EF_dir
      real(8):: alpha


      if (associated(edge%e1)) then
        cell => edge%e1

      else
        cell => edge%e2

      end if
      
      if (random() <= alpha) then
        call genericReflection(edge, part)

      else
        call genericTransparent(edge, part)

      end if

    end subroutine quasiNeutrality

    ! Generic boundary conditions for internal use
    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

    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

end submodule boundary