ComputeBdryTangentNormal Subroutine

  SUBROUTINE ComputeBdryTangentNormal( Bdry, BotTop )

    ! Does some pre-processing on the boundary points to pre-compute segment
    ! lengths  (%Len),
    ! tangents (%t, %nodet),
    ! normals  (%n, %noden), and
    ! curvatures (%kappa)
    !
    ! The boundary is also extended with a constant depth to infinity to cover cases where the ray
    ! exits the domain defined by the user

    INTEGER, SAVE                    :: NPts( 2 ) = [ 0, 0 ]
    REAL      (KIND=8)               :: p1( 3 ), p2( 3 ), p3( 3 ), p4( 3 ), U( 3 ), V( 3 )
    REAL      (KIND=8)               :: n1( 3 ), n2( 3 )      ! normal vectors to the pair of triangles
    REAL      (KIND=8)               :: tvec( 3 ), Len
    TYPE(BdryPt), INTENT( INOUT )    :: Bdry( :, : )
    CHARACTER (LEN=3),  INTENT( IN ) :: BotTop           ! Flag indicating bottom or top reflection
    CHARACTER (LEN=2), SAVE          :: CurvilinearFlag = '-'
    REAL      (KIND=8)               :: mx, my, n( 3 )

    SELECT CASE ( BotTop )
    CASE ( 'Bot' )
       NPts = NbtyPts
       CurvilinearFlag = btyType
    CASE ( 'Top' )
       NPts = NatiPts
       CurvilinearFlag = atiType
    END SELECT

    ! normals on triangle faces
    DO ix = 1, NPts( 1 ) - 1
       DO iy = 1, NPts( 2 ) - 1
          ! coordinates of corner nodes, moving counter-clockwise around the rectangle
          p1 = Bdry( ix,     iy     )%x
          p2 = Bdry( ix + 1, iy     )%x
          p3 = Bdry( ix + 1, iy + 1 )%x
          p4 = Bdry( ix,     iy + 1 )%x

          ! edges for triangle 1
          U = p2 - p1   ! tangent along one edge
          V = p3 - p1   ! tangent along another edge

          ! normal vector is the cross-product of the edge tangents
          n1( 1 ) = U( 2 ) * V( 3 ) - U( 3 ) * V( 2 )
          n1( 2 ) = U( 3 ) * V( 1 ) - U( 1 ) * V( 3 )
          n1( 3 ) = U( 1 ) * V( 2 ) - U( 2 ) * V( 1 )
          IF ( BotTop == 'Top' ) n1 = -n1

          Bdry( ix, iy )%n1 = n1 / NORM2( n1 )   ! scale to make it a unit normal

          ! edges for triangle 2
          U = p3 - p1   ! tangent along one edge
          V = p4 - p1   ! tangent along another edge

          ! normal vector is the cross-product of the edge tangents
          n2( 1 ) = U( 2 ) * V( 3 ) - U( 3 ) * V( 2 )
          n2( 2 ) = U( 3 ) * V( 1 ) - U( 1 ) * V( 3 )
          n2( 3 ) = U( 1 ) * V( 2 ) - U( 2 ) * V( 1 )
          IF ( BotTop == 'Top' ) n2 = -n2

          Bdry( ix, iy )%n2 = n2 / NORM2( n2 )   ! scale to make it a unit normal

       END DO
    END DO

    ! normals at nodes
    ! use forward, centered, or backward difference formulas
    DO ix = 1, NPts( 1 )
       DO iy = 1, NPts( 2 )
          IF ( ix == 1 ) THEN
             mx = ( Bdry( ix + 1, iy     )%x( 3 ) - Bdry( ix    , iy     )%x( 3 ) ) / &
                  ( Bdry( ix + 1, iy     )%x( 1 ) - Bdry( ix    , iy     )%x( 1 ) )
          ELSE IF ( ix == Npts( 1 ) ) THEN
             mx = ( Bdry( ix    , iy     )%x( 3 ) - Bdry( ix - 1, iy     )%x( 3 ) ) / &
                  ( Bdry( ix    , iy     )%x( 1 ) - Bdry( ix - 1, iy     )%x( 1 ) )
          ELSE
             mx = ( Bdry( ix + 1, iy     )%x( 3 ) - Bdry( ix - 1, iy     )%x( 3 ) ) / &
                  ( Bdry( ix + 1, iy     )%x( 1 ) - Bdry( ix - 1, iy     )%x( 1 ) )
          END IF

          IF ( iy == 1 ) THEN
             my = ( Bdry( ix    , iy + 1 )%x( 3 ) - Bdry( ix    , iy     )%x( 3 ) ) / &
                  ( Bdry( ix    , iy + 1 )%x( 2 ) - Bdry( ix    , iy     )%x( 2 ) )
          ELSE IF ( iy == Npts( 2 ) ) THEN
             my = ( Bdry( ix    , iy     )%x( 3 ) - Bdry( ix    , iy - 1 )%x( 3 ) ) / &
                  ( Bdry( ix    , iy     )%x( 2 ) - Bdry( ix    , iy - 1 )%x( 2 ) )
          ELSE
             my = ( Bdry( ix    , iy + 1 )%x( 3 ) - Bdry( ix    , iy - 1 )%x( 3 ) ) / &
                  ( Bdry( ix    , iy + 1 )%x( 2 ) - Bdry( ix    , iy - 1 )%x( 2 ) )
          END IF

          n = [ -mx, -my, 1.0D0 ]   ! this a normal to the surface

          IF ( ix < NPts( 1 ) .AND. iy < NPts( 2 ) ) THEN
             ! xx term
             Bdry( ix, iy )%phi_xx = atan2( n( 3 ), n( 1 ) )   ! this is the angle at each node

             ! xy term
             tvec = Bdry( ix + 1, iy + 1 )%x - Bdry( ix, iy )%x
             Len  = SQRT( tvec( 1 ) ** 2 + tvec( 2 ) ** 2 )
             tvec = tvec / Len
             Bdry( ix, iy )%phi_xy = atan2( n( 3 ), n( 1 ) * tvec( 1 ) + n( 2 ) * tvec( 2 ) )   ! this is the angle at each node

             ! yy term
             Bdry( ix, iy )%phi_yy = atan2( n( 3 ), n( 2 ) )   ! this is the angle at each node
          END IF

          Bdry( ix, iy )%Noden_unscaled = n
          Bdry( ix, iy )%Noden = n / NORM2( n )
       END DO
    END DO

    IF ( CurvilinearFlag( 1 : 1 ) == 'C' ) THEN ! curvilinear option: compute derivative as centered difference between two nodes

       ! compute curvatures in each segment
       ! ALLOCATE( phi( NPts ), Stat = IAllocStat )

       ! - sign below because the node normal = ( -mx, -my, 1 )
       DO ix = 1, NPts( 1 ) - 1
          DO iy = 1, NPts( 2 ) - 1
             ! z_xx (difference in x of z_x)

             Bdry( ix, iy )%z_xx = -( Bdry( ix + 1, iy     )%Noden_unscaled( 1 ) - Bdry( ix, iy )%Noden_unscaled( 1 ) ) / &
                                    ( Bdry( ix + 1, iy     )%x(              1 ) - Bdry( ix, iy )%x(              1 ) )

             tvec = Bdry( ix + 1, iy )%x - Bdry( ix, iy )%x
             Len  = SQRT( tvec( 1 ) ** 2 + tvec( 3 ) ** 2 )
             Bdry( ix, iy )%kappa_xx = ( Bdry( ix + 1, iy )%phi_xx - Bdry( ix, iy )%phi_xx ) / Len ! this is curvature = dphi/ds

             ! z_xy (difference in y of z_x)

             Bdry( ix, iy )%z_xy = -( Bdry( ix    , iy + 1 )%Noden_unscaled( 1 ) - Bdry( ix, iy )%Noden_unscaled( 1 ) ) / &
                                    ( Bdry( ix    , iy + 1 )%x(              2 ) - Bdry( ix, iy )%x(              2 ) )

             tvec = Bdry( ix + 1, iy + 1 )%x - Bdry( ix, iy )%x
             Len  = SQRT( tvec( 1 ) ** 2 + tvec( 2 ) ** 2 + tvec( 3 ) ** 2 )
             Bdry( ix, iy )%kappa_xy = ( Bdry( ix + 1, iy + 1 )%phi_xy - Bdry( ix, iy )%phi_xy ) / Len ! this is curvature = dphi/ds

             ! new
             tvec = Bdry( ix, iy + 1 )%x - Bdry( ix, iy )%x
             Len  = SQRT( tvec( 2 ) ** 2 + tvec( 3 ) ** 2 )
             Bdry( ix, iy )%kappa_xy = ( Bdry( ix, iy + 1 )%phi_xx - Bdry( ix, iy )%phi_xx ) / Len ! this is curvature = dphi/ds

             ! z_yy (difference in y of z_y)

             Bdry( ix, iy )%z_yy = -( Bdry( ix    , iy + 1 )%Noden_unscaled( 2 ) - Bdry( ix, iy )%Noden_unscaled( 2 ) ) / &
                                    ( Bdry( ix    , iy + 1 )%x(              2 ) - Bdry( ix, iy )%x(              2 ) )

             tvec = Bdry( ix, iy + 1 )%x - Bdry( ix, iy )%x
             Len  = SQRT( tvec( 2 ) ** 2 + tvec( 3 ) ** 2 )
             Bdry( ix, iy )%kappa_yy = ( Bdry( ix, iy + 1 )%phi_yy - Bdry( ix, iy )%phi_yy ) / Len ! this is curvature = dphi/ds

             ! introduce Len factor per Eq. 4.4.18 in Cerveny's book
             Len = NORM2( Bdry( ix, iy )%Noden_unscaled )
             Bdry( ix, iy )%z_xx = Bdry( ix, iy )%z_xx / Len
             Bdry( ix, iy )%z_xy = Bdry( ix, iy )%z_xy / Len
             Bdry( ix, iy )%z_yy = Bdry( ix, iy )%z_yy / Len
          END DO
       END DO
    ELSE
       Bdry%z_xx = 0
       Bdry%z_xy = 0
       Bdry%z_yy = 0

       Bdry%kappa_xx = 0
       Bdry%kappa_xy = 0
       Bdry%kappa_yy = 0
    END IF
!!$    write( *, * ) 'ix=1', Bdry( 1, : )%kappa_xx
!!$    write( *, * ) 'ix=1', Bdry( 1, : )%kappa_xy
!!$    write( *, * ) 'ix=1', Bdry( 1, : )%kappa_yy
!!$    write( *, * ) 'iy=1', Bdry( :, 1 )%kappa_xx
!!$    write( *, * ) 'iy=1', Bdry( :, 1 )%kappa_xy
!!$    write( *, * ) 'iy=1', Bdry( :, 1 )%kappa_yy
!!$    write( *, * ) 'D'
!!$    write( *, * ) Bdry( :, : )%z_xx
!!$    write( *, * ) Bdry( :, : )%z_xy
!!$    write( *, * ) Bdry( :, : )%z_yy
!!$
!!$    write( *, * ) 'kappa'
!!$    write( *, * ) Bdry( :, : )%kappa_xx
!!$    write( *, * ) Bdry( :, : )%kappa_xy
!!$    write( *, * ) Bdry( :, : )%kappa_yy

  END SUBROUTINE ComputeBdryTangentNormal