!moduleMeshCylBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
  USE moduleMeshCyl

  CONTAINS
    SUBROUTINE reflection(edge, part)
      USE moduleSpecies
      IMPLICIT NONE

      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)

      !TODO: Try to do this without select
      SELECT TYPE(edge)
      TYPE IS(meshEdgeCyl)
        edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
        cosT = (edge%z(2)-edge%z(1))/edgeNorm
        sinT = DSQRT(1-cosT**2)

        rp(1) = part%r(1) - edge%z(1);
        rp(2) = part%r(2) - edge%r(1);

        rpp(1) = cosT*rp(1) - sinT*rp(2)
        rpp(2) = sinT*rp(1) + cosT*rp(2)
        rpp(2) = -rpp(2)

        vpp(1) = cosT*part%v(1) - sinT*part%v(2)
        vpp(2) = sinT*part%v(1) + cosT*part%v(2)
        vpp(2) = -vpp(2)

        part%r(1) =  cosT*rpp(1) + sinT*rpp(2) + edge%z(1);
        part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%r(1);
        part%v(1) =  cosT*vpp(1) + sinT*vpp(2)
        part%v(2) = -sinT*vpp(1) + cosT*vpp(2)

      END SELECT

      part%n_in = .TRUE.

    END SUBROUTINE reflection

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

      CLASS(meshEdge), INTENT(inout):: edge
      CLASS(particle), INTENT(inout):: part
      REAL(8):: rEdge(1:2) !Position of particle projected to the edge
      REAL(8):: a, b, c
      REAL(8):: a2b2
      REAL(8):: d !Distance from particle to edge

      SELECT TYPE(edge)
      TYPE IS(meshEdgeCyl)
        a = (edge%z(1) - edge%z(2))
        b = (edge%r(1) - edge%r(2))
        c =  edge%z(1)*edge%r(2) - edge%z(2)*edge%r(1)

        a2b2 = a**2 + b**2

        rEdge(1) = (b*( b*part%r(1) - a*part%r(2)) - a*c)/a2b2
        rEdge(2) = (a*(-b*part%r(1) + a*part%r(2)) - b*c)/a2b2

        d = NORM2(rEdge - part%r(1:2))
        !Reduce weight of particle by the distance to the edge and move it to the edge
        IF (d > 0.D0) THEN
          part%weight = part%weight / d
          part%r(1:2) = rEdge

        END IF

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

        ELSE
          CALL edge%e2%scatter(part)

        END IF

      END SELECT

      !Remove particle from the domain
      part%n_in = .FALSE.

    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
      REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
      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

      !Reflects particle in the edge
      SELECT TYPE(edge)
      TYPE IS(meshEdgeCyl)
        edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
        cosT = (edge%z(2)-edge%z(1))/edgeNorm
        sinT = DSQRT(1-cosT**2)

        rp(1) = part%r(1) - edge%z(1);
        rp(2) = part%r(2) - edge%r(1);

        rpp(1) = cosT*rp(1) - sinT*rp(2)
        rpp(2) = sinT*rp(1) + cosT*rp(2)
        rpp(2) = -rpp(2)

        vpp(1) = cosT*part%v(1) - sinT*part%v(2)
        vpp(2) = sinT*part%v(1) + cosT*part%v(2)
        vpp(2) = -vpp(2)

        part%r(1) =  cosT*rpp(1) + sinT*rpp(2) + edge%z(1);
        part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%r(1);
        part%v(1) =  cosT*vpp(1) + sinT*vpp(2)
        part%v(2) = -sinT*vpp(1) + cosT*vpp(2)

      END SELECT

      part%n_in = .TRUE.

    END SUBROUTINE wallTemperature

    !Symmetry axis. Dummy function
    SUBROUTINE symmetryAxis(edge, part)
      USE moduleSpecies
      IMPLICIT NONE

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

    END SUBROUTINE symmetryAxis

END SUBMODULE moduleMeshCylBoundary