!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
      use moduleRefParam, only: v_ref
      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)
        part%v = part%v + 40e3/v_ref
        if (dot_product(part%v, edge%normal) <= 0.d0) then
          call transparent(edge, part)

        end if

        ! 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