StepToBdry2D Subroutine

  SUBROUTINE StepToBdry2D( x0, x2, urayt, h, topRefl, botRefl, &
       iSegz0, iSegr0, Topx, Topn, Botx, Botn )
    USE BdryMod
    REAL (KIND=8), INTENT( IN    ) :: x0( 2 ), urayt( 2 )
    REAL (KIND=8), INTENT( INOUT ) :: x2( 2 ), h
    INTEGER,       INTENT( IN    ) :: iSegz0, iSegr0                             ! SSP layer the ray is in
    REAL (KIND=8), INTENT( IN    ) :: Topx( 2 ), Topn( 2 ), Botx( 2 ), Botn( 2 ) ! Top, bottom coordinate and normal
    LOGICAL,       INTENT( OUT   ) :: topRefl, botRefl
    REAL (KIND=8)                  :: d( 2 ), d0( 2 ), rSeg( 2 )

    ! Original step due to maximum step size
    h = Beam%deltas
    x2 = x0 + h * urayt

    ! interface crossing in depth
    IF ( ABS( urayt( 2 ) ) > EPSILON( h ) ) THEN
       IF      ( SSP%z( iSegz0     ) > x2( 2 ) ) THEN
          h  = ( SSP%z( iSegz0     ) - x0( 2 ) ) / urayt( 2 )
          x2( 1 ) = x0( 1 ) + h * urayt( 1 )
          x2( 2 ) = SSP%z( iSegz0 )
          IF ( STEP_DEBUGGING ) THEN
             WRITE( PRTFile, * ) 'StepToBdry2D lower depth h to', h, x2
          END IF
       ELSE IF ( SSP%z( iSegz0 + 1 ) < x2( 2 ) ) THEN
          h  = ( SSP%z( iSegz0 + 1 ) - x0( 2 ) ) / urayt( 2 )
          x2( 1 ) = x0( 1 ) + h * urayt( 1 )
          x2( 2 ) = SSP%z( iSegz0 + 1 )
          IF ( STEP_DEBUGGING ) THEN
             WRITE( PRTFile, * ) 'StepToBdry2D upper depth h to', h, x2
          END IF
       END IF
    END IF

    ! ray mask using a box centered at ( 0, 0 )
    IF ( ABS( x2( 1 ) ) > Beam%Box%r ) THEN
       h = ( Beam%Box%r - ABS( x0( 1 ) ) ) / ABS( urayt( 1 ) )
       x2( 2 ) = x0( 2 ) + h * urayt( 2 )
       x2( 1 ) = SIGN( Beam%Box%r, x0( 1 ) )
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D beam box crossing R h to', h, x2
       END IF
    END IF
    IF ( ABS( x2( 2 ) ) > Beam%Box%z ) THEN
       h = ( Beam%Box%z - ABS( x0( 2 ) ) ) / ABS( urayt( 2 ) )
       x2( 1 ) = x0( 1 ) + h * urayt( 1 )
       x2( 2 ) = SIGN( Beam%Box%z, x0( 2 ) )
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D beam box crossing Z h to', h, x2
       END IF
    END IF

    ! top crossing
    d = x2 - Topx             ! vector from top to ray
    ! LP: Originally, this value had to be > a small positive number, meaning the
    ! new step really had to be outside the boundary, not just to the boundary.
    ! Also, this is not missing a normalization factor, Topn is normalized so
    ! this is actually the distance above the top in meters.
    IF ( DOT_PRODUCT( Topn, d ) > -INFINITESIMAL_STEP_SIZE ) THEN
       d0 = x0 - Topx         ! vector from top node to ray origin
       h = -DOT_PRODUCT( d0, Topn ) / DOT_PRODUCT( urayt, Topn )
       x2 = x0 + h * urayt
       ! Snap to exact top depth value if it's flat
       IF ( ABS( Topn( 1 ) ) < EPSILON( Topn( 1 ) ) ) THEN
          x2( 2 ) = Topx( 2 )
       END IF
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D top crossing h to', h, x2
       END IF
       topRefl = .TRUE.
    ELSE
       topRefl = .FALSE.
    END IF

    ! bottom crossing
    d = x2 - Botx              ! vector from bottom to ray
    ! See comment above for top case.
    IF ( DOT_PRODUCT( Botn, d ) > -INFINITESIMAL_STEP_SIZE ) THEN
       d0  = x0 - Botx         ! vector from bottom node to ray origin
       h = -DOT_PRODUCT( d0, Botn ) / DOT_PRODUCT( urayt, Botn )
       x2 = x0 + h * urayt
       ! Snap to exact bottom depth value if it's flat
       IF ( ABS( Botn( 1 ) ) < EPSILON( Botn( 1 ) ) ) THEN
          x2( 2 ) = Botx( 2 )
       END IF
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D bottom crossing h to', h, x2
       END IF
       botRefl = .TRUE.
       ! Should not ever be able to cross both, but in case it does, make sure
       ! only the crossing we exactly landed on is active
       topRefl = .FALSE.
    ELSE
       botRefl = .FALSE.
    END IF

    ! top or bottom segment crossing in range
    rSeg( 1 ) = MAX( rTopSeg( 1 ), rBotSeg( 1 ) ) ! LP: lower range bound (not an x value)
    rSeg( 2 ) = MIN( rTopSeg( 2 ), rBotSeg( 2 ) ) ! LP: upper range bound (not a y value)

    IF ( SSP%Type == 'Q' ) THEN
       rSeg( 1 ) = MAX( rSeg( 1 ), SSP%Seg%r( iSegr0     ) )
       rSeg( 2 ) = MIN( rSeg( 2 ), SSP%Seg%r( iSegr0 + 1 ) )
    END IF

    IF ( ABS( urayt( 1 ) )  > EPSILON( h ) ) THEN
       IF      ( x2( 1 ) < rSeg( 1 ) ) THEN
          h = -( x0( 1 ) - rSeg( 1 ) ) / urayt( 1 )
          x2( 1 ) = rSeg( 1 )
          x2( 2 ) = x0( 2 ) + h * urayt( 2 )
          IF ( STEP_DEBUGGING ) THEN
             WRITE( PRTFile, * ) 'StepToBdry2D lower range h to', h, x2
          END IF
          topRefl = .FALSE.
          botRefl = .FALSE.
       ELSE IF ( x2( 1 ) > rSeg( 2 ) ) THEN
          h = -( x0( 1 ) - rSeg( 2 ) ) / urayt( 1 )
          x2( 1 ) = rSeg( 2 )
          x2( 2 ) = x0( 2 ) + h * urayt( 2 )
          IF ( STEP_DEBUGGING ) THEN
             WRITE( PRTFile, * ) 'StepToBdry2D upper range h to', h, x2
          END IF
          topRefl = .FALSE.
          botRefl = .FALSE.
       END IF
    END IF

    IF ( h < INFINITESIMAL_STEP_SIZE * Beam%deltas ) THEN   ! is it taking an infinitesimal step?
       h = INFINITESIMAL_STEP_SIZE * Beam%deltas            ! make sure we make some motion
       x2 = x0 + h * urayt
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D small step forced h to', h, x2
       END IF
       ! Recheck reflection conditions
       d = x2 - Topx             ! vector from top to ray
       IF ( DOT_PRODUCT( Topn, d ) > EPSILON( h ) ) THEN
          topRefl = .TRUE.
       ELSE
          topRefl = .FALSE.
       END IF
       d = x2 - Botx              ! vector from bottom to ray
       IF ( DOT_PRODUCT( Botn, d ) > EPSILON( h ) ) THEN
          botRefl = .TRUE.
          topRefl = .FALSE.
       ELSE
          botRefl = .FALSE.
       END IF
    ELSE
       IF ( STEP_DEBUGGING ) THEN
          WRITE( PRTFile, * ) 'StepToBdry2D normal h to', h, x2
       END IF
    END IF

  END SUBROUTINE StepToBdry2D