!moduleMeshBoundary: Boundary functions
MODULE moduleMeshBoundary
  USE moduleMesh

  CONTAINS
    SUBROUTINE reflection(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
      REAL(8):: tI
      REAL(8):: taup !time step for reflecting process

      !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(edge, part)
      USE moduleCaseParam
      USE moduleSpecies
      IMPLICIT NONE

      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(part)

      ELSE
        CALL edge%e2%scatter(part)

      END IF

    END SUBROUTINE absorption

    !Transparent boundary condition
    SUBROUTINE transparent(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 transparent

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

      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      INTEGER:: i

      !Modifies particle velocity according to wall temperature
      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
      TYPE IS(boundaryWallTemperature)
        DO i = 1, 3
          part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()

        END DO
          
      END SELECT

      CALL reflection(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(edge, part)
      USE moduleList
      USE moduleSpecies
      USE moduleMesh
      USE moduleRefParam
      USE moduleRandom
      IMPLICIT NONE

      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass
      REAL(8):: ionizationRate
      INTEGER:: ionizationPair, p
      REAL(8):: v0(1:3) !random velocity of neutral
      TYPE(particle), POINTER:: newElectron
      TYPE(particle), POINTER:: newIon

      SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
      TYPE IS(boundaryIonization)
        mRel = (bound%m0*species(part%sp)%obj%m)*(bound%m0+species(part%sp)%obj%m)
        vRel = SUM(DABS(part%v-bound%v0))
        eRel = mRel*vRel**2*5.D-1

        IF (eRel > bound%eThreshold) THEN
          ionizationRate = part%weight*bound%n0*bound%crossSection%get(eRel)*vRel

          !Rounds the number of particles up
          ionizationPair = NINT(ionizationRate*bound%effectiveTime/species(bound%sp)%obj%weight)

          !Create the new pair of particles
          DO p = 1, ionizationPair
            !Assign random velocity to the neutral
            v0(1) = bound%v0(1) + bound%vTh*randomMaxwellian()
            v0(2) = bound%v0(2) + bound%vTh*randomMaxwellian()
            v0(3) = bound%v0(3) + bound%vTh*randomMaxwellian()

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

            newElectron%sp = part%sp
            newIon%sp      = bound%sp

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

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

            newElectron%vol = part%vol
            newIon%vol      = part%vol

            newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
            newIon%xi      = newElectron%xi

            newElectron%qm = part%qm
            SELECT TYPE(spe => species(bound%sp)%obj)
            TYPE IS(speciesCharged)
              newIon%qm = spe%qm

            END SELECT

            newElectron%weight = species(bound%sp)%obj%weight
            newIon%weight      = newElectron%weight

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

            !Add particles to list
            CALL OMP_SET_LOCK(lockSurfaces)
            CALL partSurfaces%add(newElectron)
            CALL partSurfaces%add(newIon)
            CALL OMP_UNSET_LOCK(lockSurfaces)

          END DO

        END IF

      END SELECT

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

    END SUBROUTINE ionization

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

      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part

      CALL reflection(edge, part)

    END SUBROUTINE symmetryAxis

    !Points the boundary function to specific type
    SUBROUTINE pointBoundaryFunction(edge, s)
      USE moduleErrors
      IMPLICIT NONE

      CLASS(meshEdge), INTENT(inout):: edge
      INTEGER, INTENT(in):: s !Species index

      SELECT TYPE(obj => edge%boundary%bTypes(s)%obj)
      TYPE IS(boundaryAbsorption)
        edge%fBoundary(s)%apply => absorption

      TYPE IS(boundaryReflection)
        edge%fBoundary(s)%apply => reflection

      TYPE IS(boundaryTransparent)
        edge%fBoundary(s)%apply => transparent

      TYPE IS(boundaryWallTemperature)
        edge%fBoundary(s)%apply => wallTemperature

      TYPE IS(boundaryIonization)
        edge%fBoundary(s)%apply => ionization

      TYPE IS(boundaryAxis)
        edge%fBoundary(s)%apply => symmetryAxis

      CLASS DEFAULT
        CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')

      END SELECT

    END SUBROUTINE pointBoundaryFunction

END MODULE moduleMeshBoundary