TraceRay2D Subroutine

SUBROUTINE TraceRay2D( xs, alpha, Amp0 )
!! Traces the beam corresponding to a particular take-off angle

  USE Step
  USE WriteRay

  REAL     (KIND=8), INTENT( IN ) :: xs( 2 )      ! x-y coordinate of the source
  REAL     (KIND=8), INTENT( IN ) :: alpha, Amp0  ! initial angle, amplitude
  INTEGER           :: is, is1                    ! index for a step along the ray
  REAL     (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho
  REAL     (KIND=8) :: dEndTop( 2 ), dEndBot( 2 ), TopnInt( 2 ), BotnInt( 2 ), ToptInt( 2 ), BottInt( 2 )
  REAL     (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom
  REAL     (KIND=8) :: sss
  REAL     (KIND=8) :: tinit( 2 )
  LOGICAL           :: topRefl, botRefl

  ! Initial conditions

  ! LP: This initialization was missing (only done once globally). In some cases
  ! (a source on a boundary), in conjunction with other subtle issues, the
  ! missing initialization caused the initial state of one ray to be dependent
  ! on the final state of the previous ray, which is obviously non-physical.
  iSegz = 1
  iSegr = 1

  iSmallStepCtr = 0
  tinit = [ COS( alpha ), SIN( alpha ) ]
  CALL EvaluateSSP( xs, tinit, c, cimag, gradc, crr, crz, czz, rho, freq, 'TAB' )
  ray2D( 1 )%c         = c
  ray2D( 1 )%x         = xs
  ray2D( 1 )%t         = tinit / c
  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

  ! second component of qv is not used in geometric beam tracing
  ! set I.C. to 0 in hopes of saving run time
  IF ( Beam%RunType( 2 : 2 ) == 'G' ) ray2D( 1 )%q = [ 0.0, 0.0 ]

  CALL GetTopSeg( xs( 1 ), ray2D( 1 )%t( 1 ) )   ! identify the top    segment above the source
  CALL GetBotSeg( xs( 1 ), ray2D( 1 )%t( 1 ) )   ! identify the bottom segment below the source

  ! convert range-dependent geoacoustic parameters from user to program units
  ! compiler is not accepting the copy of the whole structure at once ...
  IF ( atiType( 2 : 2 ) == 'L' ) THEN
     Bdry%Top%HS%cp  = Top( IsegTop )%HS%cp   ! grab the geoacoustic info for the new segment
     Bdry%Top%HS%cs  = Top( IsegTop )%HS%cs
     Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
  END IF

  IF ( btyType( 2 : 2 ) == 'L' ) THEN
     Bdry%Bot%HS%cp  = Bot( IsegBot )%HS%cp
     Bdry%Bot%HS%cs  = Bot( IsegBot )%HS%cs
     Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
  END IF

  ! Trace the beam (note that Reflect alters the step index, is)
  is = 0
  CALL Distances2D( ray2D( 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop,    dEndBot,  &
       Top( IsegTop )%n, Bot( IsegBot )%n, DistBegTop, DistBegBot )

  IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
     Beam%Nsteps = 1
     WRITE( PRTFile, * ) 'Terminating the ray trace because the source is on or outside the boundaries'
     RETURN       ! source must be within the medium
  END IF

  ! LP: Changed from loop istep to MaxN-1, as is will hit MaxN-3 first.
  Stepping: DO WHILE ( .TRUE. )
     is  = is + 1
     is1 = is + 1

     CALL Step2D( ray2D( is ), ray2D( is1 ),  &
          topRefl, botRefl, &
          Top( IsegTop )%x, Top( IsegTop )%n, &
          Bot( IsegBot )%x, Bot( IsegBot )%n )

     ! New altimetry segment?
     IF (        ray2D( is1 )%x( 1 )  < rTopSeg( 1 ) &
          .OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) <  0.0 ) &
          .OR.   ray2D( is1 )%x( 1 )  > rTopSeg( 2 ) &
          .OR. ( ray2D( is1 )%x( 1 ) == rTopSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
        CALL GetTopSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
        IF ( atiType( 2 : 2 ) == 'L' ) THEN
           Bdry%Top%HS%cp  = Top( IsegTop )%HS%cp   ! grab the geoacoustic info for the new segment
           Bdry%Top%HS%cs  = Top( IsegTop )%HS%cs
           Bdry%Top%HS%rho = Top( IsegTop )%HS%rho
        END IF
     END IF

     ! New bathymetry segment?
     IF (        ray2D( is1 )%x( 1 )  < rBotSeg( 1 ) &
          .OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 1 ) .AND. ray2D( is1 )%t( 1 ) <  0.0 ) &
          .OR.   ray2D( is1 )%x( 1 )  > rBotSeg( 2 ) &
          .OR. ( ray2D( is1 )%x( 1 ) == rBotSeg( 2 ) .AND. ray2D( is1 )%t( 1 ) >= 0.0 )) THEN
        CALL GetBotSeg( ray2D( is1 )%x( 1 ), ray2D( is1 )%t( 1 ) )
        IF ( btyType( 2 : 2 ) == 'L' ) THEN
           Bdry%Bot%HS%cp  = Bot( IsegBot )%HS%cp   ! grab the geoacoustic info for the new segment
           Bdry%Bot%HS%cs  = Bot( IsegBot )%HS%cs
           Bdry%Bot%HS%rho = Bot( IsegBot )%HS%rho
        END IF
     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 Distances2D( ray2D( is1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop,    dEndBot,  &
          Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )

     IF      ( topRefl ) THEN
        ! WRITE( PRTFile, * ) 'Top reflecting'
        IF ( atiType == 'C' ) THEN
           sss     = DOT_PRODUCT( dEndTop, Top( IsegTop )%t ) / Top( IsegTop )%Len   ! proportional distance along segment
           TopnInt = ( 1 - sss ) * Top( IsegTop )%Noden + sss * Top( 1 + IsegTop )%Noden
           ToptInt = ( 1 - sss ) * Top( IsegTop )%Nodet + sss * Top( 1 + IsegTop )%Nodet
        ELSE
           TopnInt = Top( IsegTop )%n   ! normal is constant in a segment
           ToptInt = Top( IsegTop )%t
        END IF

        CALL Reflect2D( is, Bdry%Top%HS, 'TOP', ToptInt, TopnInt, Top( IsegTop )%kappa, RTop, NTopPTS )
        ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1

        CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop,    dEndBot,  &
             Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )

     ELSE IF ( botRefl ) THEN
        ! WRITE( PRTFile, * ) 'Bottom reflecting'
        IF ( btyType == 'C' ) THEN
           sss     = DOT_PRODUCT( dEndBot, Bot( IsegBot )%t ) / Bot( IsegBot )%Len   ! proportional distance along segment
           BotnInt = ( 1 - sss ) * Bot( IsegBot )%Noden + sss * Bot( 1 + IsegBot )%Noden
           BottInt = ( 1 - sss ) * Bot( IsegBot )%Nodet + sss * Bot( 1 + IsegBot )%Nodet
        ELSE
           BotnInt = Bot( IsegBot )%n   ! normal is constant in a segment
           BottInt = Bot( IsegBot )%t
        END IF

        CALL Reflect2D( is, Bdry%Bot%HS, 'BOT', BottInt, BotnInt, Bot( IsegBot )%kappa, RBot, NBotPTS )
        ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1
        CALL Distances2D( ray2D( is + 1 )%x, Top( IsegTop )%x, Bot( IsegBot )%x, dEndTop,    dEndBot, &
             Top( IsegTop )%n, Bot( IsegBot )%n, DistEndTop, DistEndBot )

     END IF

     ! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit?
     ! There is a test here for when two adjacent ray points are outside the boundaries
     ! BELLHOP3D handles this differently using a counter-limit for really small steps.
     IF ( ABS( ray2D( is + 1 )%x( 1 ) ) >= Beam%Box%r .OR. &
          ABS( ray2D( is + 1 )%x( 2 ) ) >= Beam%Box%z .OR. ray2D( is + 1 )%Amp < 0.005 .OR. &
          ( DistBegTop < 0.0 .AND. DistEndTop < 0.0 ) .OR. &
          ( DistBegBot < 0.0 .AND. DistEndBot < 0.0 ) ) THEN
          ! ray2D( is + 1 )%t( 1 ) < 0 ) THEN ! this last test kills off a backward traveling ray
        Beam%Nsteps = is + 1
        EXIT Stepping
     ELSE IF ( is >= MaxN - 3 ) THEN
        WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory'
        Beam%Nsteps = is
        EXIT Stepping
     END IF

     DistBegTop = DistEndTop
     DistBegBot = DistEndBot

  END DO Stepping
END SUBROUTINE TraceRay2D