!moduleMeshBoundary: Boundary functions for the mesh edges
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

      !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(edge%e1%nNodes, part)

      ELSE
        CALL edge%e2%scatter(edge%e2%nNodes, 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

    !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

    !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%species%n)%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
      USE moduleMath
      IMPLICIT NONE

      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

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

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

        DO p = 1, nIonizations
          !Get probability of ionization
          pIonization = 1.D0 - DEXP(-bound%n0*bound%crossSection%get(eRel)*vRel*bound%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) = 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)

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

            ELSE
              newElectron%species => part%species

            END IF
            newIon%species      => bound%species

            newElectron%v = v0 + (1.D0 + bound%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 - bound%eThreshold
            vRel = 2.D0*DSQRT(eRel)/mRel

            !Reduce number of possible ionizations
            nIonizations = nIonizations - 1

          END IF

        END DO

      END SELECT

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

    END SUBROUTINE ionization

    subroutine outflowAdaptive(edge, part)
      use moduleRandom
      implicit none

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

      select type(bound => edge%boundary%bTypes(part%species%n)%obj)
      type is(boundaryOutflowAdaptive)

        if (random() < 0.844d0) then
          call reflection(edge, part)

        else
          call transparent(edge, part)

        end if

      end select

    end subroutine outflowAdaptive

    !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(boundaryAxis)
        edge%fBoundary(s)%apply => symmetryAxis

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

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

      type is(boundaryOutflowAdaptive)
        edge%fBoundary(s)%apply => outflowAdaptive

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

      END SELECT

    END SUBROUTINE pointBoundaryFunction

END MODULE moduleMeshBoundary