SUBROUTINE StepToBdry3D( x0, x2, urayt, iSegx0, iSegy0, iSegz0, h, &
topRefl, botRefl, snapDim )
INTEGER, INTENT( IN ) :: iSegx0, iSegy0, iSegz0
REAL (KIND=8), INTENT( IN ) :: x0( 3 ), urayt( 3 ) ! ray coordinate and tangent
REAL (KIND=8), INTENT( INOUT ) :: x2( 3 ), h ! output coord, reduced step size
LOGICAL, INTENT( OUT ) :: topRefl, botRefl
INTEGER, INTENT( OUT ) :: snapDim
REAL (KIND=8) :: d( 3 ), d0( 3 ), tri_n( 3 )
REAL (KIND=8) :: xSeg( 2 ), ySeg( 2 ) ! boundary limits
INTEGER :: k
REAL (KIND=8) :: hnew, over_diag_amount
LOGICAL :: newSide, newOnEdge
! Original step due to maximum step size
h = Beam%deltas
x2 = x0 + h * urayt
snapDim = -1
! interface crossing in depth
IF ( ABS( urayt( 3 ) ) > EPSILON( h ) ) THEN
IF ( SSP%z( iSegz0 ) > x2( 3 ) .AND. iSegz0 > 1 ) THEN
h = ( SSP%z( iSegz0 ) - x0( 3 ) ) / urayt( 3 )
x2 = x0 + h * urayt
x2( 3 ) = SSP%z( iSegz0 )
snapDim = 2
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D shallower h to', h, x2
ELSE IF ( SSP%z( iSegz0 + 1 ) < x2( 3 ) .AND. iSegz0 + 1 < SSP%Nz ) THEN
h = ( SSP%z( iSegz0 + 1 ) - x0( 3 ) ) / urayt( 3 )
x2 = x0 + h * urayt
x2( 3 ) = SSP%z( iSegz0 + 1 )
snapDim = 2
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D deeper h to', h, x2
END IF
END IF
! ray mask using a box centered at ( xs_3D( 1 ), xs_3D( 2 ), 0 )
IF ( ABS( x2( 1 ) - xs_3D( 1 ) ) > Beam%Box%x ) THEN
h = ( Beam%Box%x - ABS( ( x0( 1 ) - xs_3D( 1 ) ) ) ) / ABS( urayt( 1 ) )
x2 = x0 + h * urayt
x2( 1 ) = xs_3D( 1 ) + SIGN( Beam%Box%x, x0( 1 ) - xs_3D( 1 ) )
snapDim = 0
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D beam box crossing X h to', h, x2
END IF
IF ( ABS( x2( 2 ) - xs_3D( 2 ) ) > Beam%Box%y ) THEN
h = ( Beam%Box%y - ABS( ( x0( 2 ) - xs_3D( 2 ) ) ) ) / ABS( urayt( 2 ) )
x2 = x0 + h * urayt
x2( 2 ) = xs_3D( 2 ) + SIGN( Beam%Box%y, x0( 2 ) - xs_3D( 2 ) )
snapDim = 1
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D beam box crossing Y h to', h, x2
END IF
IF ( ABS( x2( 3 ) ) > Beam%Box%z ) THEN
h = ( Beam%Box%z - ABS( x0( 3 ) ) ) / ABS( urayt( 3 ) )
x2 = x0 + h * urayt
x2( 3 ) = SIGN( Beam%Box%z, x0( 3 ) )
snapDim = 2
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D beam box crossing Y h to', h, x2
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 ) ) .AND. ABS( Topn( 2 ) ) < EPSILON( Topn( 2 ) ) ) THEN
x2( 3 ) = Topx( 3 )
END IF
snapDim = 2 ! Even if not flat, exactness in Z most important
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D top crossing h to', h, x2
topRefl = .TRUE.
ELSE
topRefl = .FALSE.
END IF
! bottom crossing
d = x2 - Botx ! vector from bottom to ray
! LP: 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 ) ) .AND. ABS( Botn( 2 ) ) < EPSILON( Botn( 2 ) ) ) THEN
x2( 3 ) = Botx( 3 )
END IF
snapDim = 2 ! Even if not flat, exactness in Z most important
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D bottom crossing h to', h, x2
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/bottom segment crossing in x
xSeg( 1 ) = MAX( xTopSeg( 1 ), xBotSeg( 1 ) ) ! LP: lower range bound (not an x value)
xSeg( 2 ) = MIN( xTopSeg( 2 ), xBotSeg( 2 ) ) ! LP: upper range bound (not a y value)
IF ( SSP%Type == 'H' ) THEN ! ocean segment
xSeg( 1 ) = MAX( xSeg( 1 ), SSP%Seg%x( iSegx0 ) )
xSeg( 2 ) = MIN( xSeg( 2 ), SSP%Seg%x( iSegx0 + 1 ) )
END IF
IF ( ABS( urayt( 1 ) ) > EPSILON( h ) ) THEN
IF ( x2( 1 ) < xSeg( 1 ) ) THEN
h = -( x0( 1 ) - xSeg( 1 ) ) / urayt( 1 )
x2 = x0 + h * urayt
x2( 1 ) = xSeg( 1 )
snapDim = 0
topRefl = .FALSE.
botRefl = .FALSE.
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D X min bound h to', h, x2
ELSE IF ( x2( 1 ) > xSeg( 2 ) ) THEN
h = -( x0( 1 ) - xSeg( 2 ) ) / urayt( 1 )
x2 = x0 + h * urayt
x2( 1 ) = xSeg( 2 )
snapDim = 0
topRefl = .FALSE.
botRefl = .FALSE.
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D X max bound h to', h, x2
END IF
END IF
! top/bottom segment crossing in y
ySeg( 1 ) = MAX( yTopSeg( 1 ), yBotSeg( 1 ) )
ySeg( 2 ) = MIN( yTopSeg( 2 ), yBotSeg( 2 ) )
IF ( SSP%Type == 'H' ) THEN ! ocean segment
ySeg( 1 ) = MAX( ySeg( 1 ), SSP%Seg%y( iSegy0 ) )
ySeg( 2 ) = MIN( ySeg( 2 ), SSP%Seg%y( iSegy0 + 1 ) )
END IF
! LP: removed 1000 * epsilon which mbp had comment questioning also
IF ( ABS( urayt( 2 ) ) > EPSILON( h ) ) THEN
IF ( x2( 2 ) < ySeg( 1 ) ) THEN
h = -( x0( 2 ) - ySeg( 1 ) ) / urayt( 2 )
x2 = x0 + h * urayt
x2( 2 ) = ySeg( 1 )
snapDim = 1
topRefl = .FALSE.
botRefl = .FALSE.
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Y min bound h to', h, x2
ELSE IF ( x2( 2 ) > ySeg( 2 ) ) THEN
h = -( x0( 2 ) - ySeg( 2 ) ) / urayt( 2 )
x2 = x0 + h * urayt
x2( 2 ) = ySeg( 2 )
snapDim = 1
topRefl = .FALSE.
botRefl = .FALSE.
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Y max bound h to', h, x2
END IF
END IF
! triangle crossing within a top segment
IF ( .NOT. Top_td_onEdge ) THEN
d = x2 - Topxmid ! vector from top center to ray end
d0 = x0 - Topxmid ! vector from top center to ray origin
tri_n = [ -( yTopSeg( 2 ) - yTopSeg( 1 ) ), xTopSeg( 2 ) - xTopSeg( 1 ), 0.0d0 ]
tri_n = tri_n / NORM2( tri_n )
over_diag_amount = DOT_PRODUCT( d, tri_n )
newSide = ( over_diag_amount >= 0.0D0 )
IF ( newSide .NEQV. Top_td_side ) THEN
newOnEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
IF ( newOnEdge ) THEN
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Top naturally stepped to tri diag edge, h over_diag_amount', &
over_diag_amount
ELSE
hnew = -DOT_PRODUCT( d0, tri_n ) / DOT_PRODUCT( urayt, tri_n )
IF ( hnew < 0.0D0 ) THEN
WRITE( PRTFile, * ) 'StepToBdry3D Top BHC_WARN_TRIDIAG_H_NEGATIVE'
h = 0.0D0
ELSE
IF ( hnew >= h ) &
WRITE( PRTFile, * ) 'StepToBdry3D Top BHC_WARN_TRIDIAG_H_GROWING'
h = hnew
END IF
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Top tri diag crossing hnew dot(n, d0) dot(n, d)', &
hnew, DOT_PRODUCT( tri_n, d0 ), over_diag_amount
x2 = x0 + h * urayt
Top_td_justSteppedTo = .TRUE.
Top_td_outgoingSide = newSide
snapDim = -2
topRefl = .FALSE.
botRefl = .FALSE.
END IF
END IF
END IF
! triangle crossing within a bottom segment
IF ( .NOT. Bot_td_onEdge ) THEN
d = x2 - Botxmid ! vector from bottom center to ray end
d0 = x0 - Botxmid ! vector from bottom center to ray origin
tri_n = [ -( yBotSeg( 2 ) - yBotSeg( 1 ) ), xBotSeg( 2 ) - xBotSeg( 1 ), 0.0d0 ]
tri_n = tri_n / NORM2( tri_n )
over_diag_amount = DOT_PRODUCT( d, tri_n )
newSide = ( over_diag_amount >= 0.0D0 )
IF ( newSide .NEQV. Bot_td_side ) THEN
newOnEdge = ( ABS( over_diag_amount ) < TRIDIAG_THRESH )
IF ( newOnEdge ) THEN
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Bot naturally stepped to tri diag edge, h over_diag_amount', &
over_diag_amount
ELSE
hnew = -DOT_PRODUCT( d0, tri_n ) / DOT_PRODUCT( urayt, tri_n )
IF ( hnew < 0.0D0 ) THEN
WRITE( PRTFile, * ) 'StepToBdry3D Bot BHC_WARN_TRIDIAG_H_NEGATIVE'
h = 0.0D0
ELSE
IF ( hnew >= h ) &
WRITE( PRTFile, * ) 'StepToBdry3D Bot BHC_WARN_TRIDIAG_H_GROWING'
h = hnew
END IF
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D Bot tri diag crossing hnew dot(n, d0) dot(n, d)', &
hnew, DOT_PRODUCT( tri_n, d0 ), over_diag_amount
x2 = x0 + h * urayt
Bot_td_justSteppedTo = .TRUE.
Bot_td_outgoingSide = newSide
snapDim = -2
topRefl = .FALSE.
botRefl = .FALSE.
END IF
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
snapDim = -1
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D small step forced h to ', h, x2
! Recheck reflection conditions
d = x2 - Topx ! vector from top to ray
IF ( DOT_PRODUCT( Topn, d ) > EPSILON( d( 1 ) ) ) THEN
topRefl = .TRUE.
Top_td_justSteppedTo = .FALSE.
Bot_td_justSteppedTo = .FALSE.
snapDim = 2
ELSE
topRefl = .FALSE.
END IF
d = x2 - Botx ! vector from bottom to ray
IF ( DOT_PRODUCT( Botn, d ) > EPSILON( d( 1 ) ) ) THEN
botRefl = .TRUE.
topRefl = .FALSE.
Top_td_justSteppedTo = .FALSE.
Bot_td_justSteppedTo = .FALSE.
snapDim = 2
ELSE
botRefl = .FALSE.
END IF
ELSE
IF ( STEP_DEBUGGING ) &
WRITE( PRTFile, * ) 'StepToBdry3D regular h to ', h, x2
END IF
END SUBROUTINE StepToBdry3D
