SUBROUTINE TraceRay2D( alpha, beta, Amp0 )
! Traces the beam corresponding to a particular take-off angle
USE ReflectMod
REAL (KIND=8), INTENT( IN ) :: alpha, beta, Amp0 ! initial angles, amplitude
INTEGER :: is, is1 ! index for a step along the ray
REAL (KIND=8) :: x( 3 ) ! ray coordinate
REAL (KIND=8) :: t_o( 3 ) ! tangent in ocean space
!REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
REAL (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom
REAL (KIND=8) :: tinit( 2 ), tradial( 2 ), BotnInt( 3 ), TopnInt( 3 ), s1, s2
REAL (KIND=8) :: z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy
REAL (KIND=8) :: term_minx, term_miny, term_maxx, term_maxy
LOGICAL :: topRefl, botRefl
! *** Initial conditions ***
iSmallStepCtr = 0
tinit = [ COS( alpha ), SIN( alpha ) ]
tradial = [ COS( beta ), SIN( beta ) ]
ray2D( 1 )%x = [ 0.0D0, xs_3D( 3 ) ]
! CALL EvaluateSSP2D( ray2D( 1 )%x, tinit, c, cimag, gradc, crr, crz, czz, rho, xs_3D, tradial, freq )
ray2D( 1 )%t = tinit / ray2D( 1 )%c ! LP: Set before PickEpsilon independent of ray angles; never overwritten
ray2D( 1 )%p = [ 1.0, 0.0 ]
ray2D( 1 )%q = [ 0.0, 1.0 ]
ray2D( 1 )%tau = 0.0
ray2D( 1 )%Amp = Amp0
ray2D( 1 )%Phase = 0.0
ray2D( 1 )%NumTopBnc = 0
ray2D( 1 )%NumBotBnc = 0
IsegTopx = 1
IsegTopy = 1
IsegBotx = 1
IsegBoty = 1
t_o = RayToOceanT( ray2D( 1 )%t, tradial )
CALL GetTopSeg3D( xs_3D, t_o, .TRUE. ) ! identify the top segment above the source
CALL GetBotSeg3D( xs_3D, t_o, .TRUE. ) ! identify the bottom segment below the source
! Trace the beam (note that Reflect alters the step index is)
is = 0
CALL Distances3D( xs_3D, Topx, Botx, Topn, Botn, DistBegTop, DistBegBot )
IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
Beam%Nsteps = 1
RETURN ! source must be within the medium
END IF
Stepping: DO WHILE ( .TRUE. )
is = is + 1
is1 = is + 1
CALL Step2D( ray2D( is ), ray2D( is1 ), tradial, topRefl, botRefl )
! convert polar coordinate of ray to x-y coordinate
x = RayToOceanX( ray2D( is1 )%x, xs_3D, tradial )
t_o = RayToOceanT( ray2D( is1 )%t, tradial )
CALL GetTopSeg3D( x, t_o, .FALSE. ) ! identify the top segment above the source
CALL GetBotSeg3D( x, t_o, .FALSE. ) ! identify the bottom segment below the source
IF ( IsegTopx == 0 .OR. IsegTopy == 0 .OR. IsegBotx == 0 .OR. IsegBoty == 0 ) THEN ! we escaped the box
Beam%Nsteps = is
EXIT Stepping
END IF
! Reflections?
! Tests that ray at step is is inside, and ray at step is+1 is outside
! to detect only a crossing from inside to outside
! DistBeg is the distance at step is, which is saved
! DistEnd is the distance at step is+1, which needs to be calculated
CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )
IF ( topRefl ) THEN
IF ( atiType == 'C' ) THEN
! LP: This is superfluous, it's the same x calculated above.
x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )
s1 = ( x( 1 ) - Topx( 1 ) ) / ( xTopSeg( 2 ) - xTopSeg( 1 ) ) ! proportional distance along segment
s2 = ( x( 2 ) - Topx( 2 ) ) / ( yTopSeg( 2 ) - yTopSeg( 1 ) ) ! proportional distance along segment
TopnInt = Top( IsegTopx, IsegTopy )%Noden * ( 1 - s1 ) * ( 1 - s2 ) + &
Top( IsegTopx + 1, IsegTopy )%Noden * ( s1 ) * ( 1 - s2 ) + &
Top( IsegTopx + 1, IsegTopy + 1 )%Noden * ( s1 ) * ( s2 ) + &
Top( IsegTopx, IsegTopy + 1 )%Noden * ( 1 - s1 ) * ( s2 )
z_xx = Top( IsegTopx, IsegTopy )%z_xx
z_xy = Top( IsegTopx, IsegTopy )%z_xy
z_yy = Top( IsegTopx, IsegTopy )%z_yy
kappa_xx = Top( IsegTopx, IsegTopy )%kappa_xx
kappa_xy = Top( IsegTopx, IsegTopy )%kappa_xy
kappa_yy = Top( IsegTopx, IsegTopy )%kappa_yy
ELSE
TopnInt = Topn ! normal is constant in a segment
z_xx = 0
z_xy = 0
z_yy = 0
kappa_xx = 0
kappa_xy = 0
kappa_yy = 0
END IF
CALL Reflect2D( is, Bdry%Top%HS, 'TOP', TopnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RTop, NTopPTS, tradial)
ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1
x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )
CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )
ELSE IF ( botRefl ) THEN
! write( *, * ) 'Reflecting', x, Botx
! write( *, * ) 'Botn', Botn
! write( *, * ) 'Distances', DistEndTop, DistEndBot
IF ( btyType == 'C' ) THEN
x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )
s1 = ( x( 1 ) - Botx( 1 ) ) / ( xBotSeg( 2 ) - xBotSeg( 1 ) ) ! proportional distance along segment
s2 = ( x( 2 ) - Botx( 2 ) ) / ( yBotSeg( 2 ) - yBotSeg( 1 ) ) ! proportional distance along segment
BotnInt = Bot( IsegBotx, IsegBoty )%Noden * ( 1 - s1 ) * ( 1 - s2 ) + &
Bot( IsegBotx + 1, IsegBoty )%Noden * ( s1 ) * ( 1 - s2 ) + &
Bot( IsegBotx + 1, IsegBoty + 1 )%Noden * ( s1 ) * ( s2 ) + &
Bot( IsegBotx, IsegBoty + 1 )%Noden * ( 1 - s1 ) * ( s2 )
z_xx = Bot( IsegBotx, IsegBoty )%z_xx
z_xy = Bot( IsegBotx, IsegBoty )%z_xy
z_yy = Bot( IsegBotx, IsegBoty )%z_yy
kappa_xx = Bot( IsegBotx, IsegBoty )%kappa_xx
kappa_xy = Bot( IsegBotx, IsegBoty )%kappa_xy
kappa_yy = Bot( IsegBotx, IsegBoty )%kappa_yy
ELSE
BotnInt = Botn ! normal is constant in a segment
z_xx = 0
z_xy = 0
z_yy = 0
kappa_xx = 0
kappa_xy = 0
kappa_yy = 0
END IF
CALL Reflect2D( is, Bdry%Bot%HS, 'BOT', BotnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RBot, NBotPTS, tradial)
ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1
x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial )
CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )
END IF
! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit?
! LP: See explanation for changes in bdry3DMod: GetTopSeg3D.
t_o = RayToOceanT( ray2D( is + 1 )%t, tradial )
term_minx = MAX( Bot( 1, 1 )%x( 1 ), Top( 1, 1 )%x( 1 ) )
term_miny = MAX( Bot( 1, 1 )%x( 2 ), Top( 1, 1 )%x( 2 ) )
term_maxx = MIN( Bot( NBTYPts( 1 ), 1 )%x( 1 ), Top( NATIPts( 1 ), 1 )%x( 1 ) )
term_maxy = MIN( Bot( 1, NBTYPts( 2 ) )%x( 2 ), Top( 1, NATIPts( 2 ) )%x( 2 ) )
! LP: The Beam%Box conditions were inexplicably commented out in 2022 revision of Nx2D, see README.
IF ( ABS( x( 1 ) - xs_3D( 1 ) ) >= Beam%Box%x .OR. &
ABS( x( 2 ) - xs_3D( 2 ) ) >= Beam%Box%y .OR. &
ABS( x( 3 ) ) >= Beam%Box%z .OR. & ! LP: Removed xs( 3 ) for consistency
x( 1 ) < term_minx .OR. &
x( 2 ) < term_miny .OR. &
x( 1 ) > term_maxx .OR. &
x( 2 ) > term_maxy .OR. &
( x( 1 ) == term_minx .AND. t_o( 1 ) < 0.0 ) .OR. &
( x( 2 ) == term_miny .AND. t_o( 2 ) < 0.0 ) .OR. &
( x( 1 ) == term_maxx .AND. t_o( 1 ) > 0.0 ) .OR. &
( x( 2 ) == term_maxy .AND. t_o( 2 ) > 0.0 ) .OR. &
ray2D( is + 1 )%Amp < 0.005 .OR. &
! ray2D( is + 1 )%t( 1 ) < 0 .OR. & ! kills off a backward traveling ray
iSmallStepCtr > 50 ) THEN
Beam%Nsteps = is + 1
EXIT Stepping
ELSE IF ( is >= MaxN - 3 ) THEN
WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory'
WRITE( PRTFile, * ) 'Angles are alpha = ', alpha * RadDeg, ' beta = ', beta * RadDeg
Beam%Nsteps = is
EXIT Stepping
END IF
DistBegTop = DistEndTop
DistBegBot = DistEndBot
END DO Stepping
END SUBROUTINE TraceRay2D
