Step2D Subroutine

SUBROUTINE Step2D( ray0, ray2, tradial, topRefl, botRefl )

  ! Does a single step along the ray
  ! x denotes the ray coordinate, (r,z)
  ! t denotes the scaled tangent to the ray (previously (rho, zeta))
  ! c * t would be the unit tangent

  USE Step3DMod

  REAL (KIND=8), INTENT( IN ) :: tradial( 2 )   ! coordinate of source and ray bearing angle
  LOGICAL, INTENT( OUT ) :: topRefl, botRefl
  TYPE( ray2DPt )    :: ray1
  TYPE( ray2DPt ), INTENT( IN )    :: ray0
  TYPE( ray2DPt ), INTENT( INOUT ) :: ray2
  INTEGER            :: iSegx0, iSegy0, iSegz0, snapDim
  REAL     (KIND=8 ) :: gradc0( 2 ), gradc1( 2 ), gradc2( 2 ), rho, &
                        c0, cimag0, crr0, crz0, czz0, csq0, cnn0_csq0, &
                        c1, cimag1, crr1, crz1, czz1, csq1, cnn1_csq1, &
                        c2, cimag2, crr2, crz2, czz2, urayt0( 2 ), urayt1( 2 ), urayt2( 2 ), &
                        h, halfh, hw0, hw1, ray2n( 2 ), RM, RN, gradcjump( 2 ), cnjump, csjump, w0, w1, &
                        rayx3D( 3 ), rayt3D( 3 ), x2_o( 3 ), x2_o_out( 3 )

  IF ( STEP_DEBUGGING ) THEN
     WRITE( PRTFile, * )
     WRITE( PRTFile, * ) 'ray0 x t', ray0%x, ray0%t
     WRITE( PRTFile, * ) 'iSegr iSegz', iSegr, iSegz
  END IF

  ! The numerical integrator used here is a version of the polygon (a.k.a. midpoint, leapfrog, or Box method), and similar
  ! to the Heun (second order Runge-Kutta method).
  ! However, it's modified to allow for a dynamic step change, while preserving the second-order accuracy).

  ! *** Phase 1 (an Euler step)

  CALL EvaluateSSP2D( ray0%x, ray0%t, c0, cimag0, gradc0, crr0, crz0, czz0, rho, xs_3D, tradial, freq )

  csq0      = c0 * c0
  cnn0_csq0 = crr0 * ray0%t( 2 )**2 - 2.0 * crz0 * ray0%t( 1 ) * ray0%t( 2 ) + czz0 * ray0%t( 1 )**2

  iSegx0 = iSegx     ! make note of current layer
  iSegy0 = iSegy
  iSegz0 = iSegz

  h = Beam%deltas            ! initially set the step h, to the basic one, deltas

  urayt0 = c0 * ray0%t  ! unit tangent
  rayx3D = RayToOceanX( ray0%x, xs_3D, tradial )
  rayt3D = RayToOceanT( urayt0, tradial )

  CALL ReduceStep3D( rayx3D, rayt3D, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary

  halfh = 0.5 * h   ! first step of the modified polygon method is a half step

  ray1%x = ray0%x + halfh * urayt0
  ray1%t = ray0%t - halfh * gradc0 / csq0
  ray1%p = ray0%p - halfh * cnn0_csq0 * ray0%q
  ray1%q = ray0%q + halfh * c0        * ray0%p

  ! *** Phase 2

  CALL EvaluateSSP2D( ray1%x, ray1%t, c1, cimag1, gradc1, crr1, crz1, czz1, rho, xs_3D, tradial, freq )
  csq1      = c1 * c1
  cnn1_csq1 = crr1 * ray1%t( 2 )**2 - 2.0 * crz1 * ray1%t( 1 ) * ray1%t( 2 ) + czz1 * ray1%t( 1 )**2

  ! The Munk test case with a horizontally launched ray caused problems.
  ! The ray vertexes on an interface and can ping-pong around that interface.
  ! Have to be careful in that case about big changes to the stepsize (that invalidate the leap-frog scheme) in phase II.
  ! A modified Heun or Box method could also work.

  urayt1 = c1 * ray1%t   ! unit tangent
  rayt3D = RayToOceanT( urayt1, tradial )

  CALL ReduceStep3D( rayx3D, rayt3D, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary

  ! use blend of f' based on proportion of a full step used.
  w1  = h / ( 2.0d0 * halfh )
  w0  = 1.0d0 - w1
  urayt2 = w0 * urayt0 + w1 * urayt1
  rayt3D = RayToOceanT( urayt2, tradial )
  ! Take the blended ray tangent ( urayt2 ) and find the minimum step size ( h )
  ! to put this on a boundary, and ensure that the resulting position
  ! ( ray2%x ) gets put precisely on the boundary.
  CALL StepToBdry3D( rayx3D, x2_o, rayt3D, iSegx0, iSegy0, iSegz0, h, &
     topRefl, botRefl, snapDim )
  ray2%x = OceanToRayX( x2_o, xs_3D, tradial, urayt2, snapDim )
  IF ( STEP_DEBUGGING ) THEN
     x2_o_out = RayToOceanX( ray2%x, xs_3D, tradial )
     WRITE( PRTFile, * ) 'OceanToRayX in ', x2_o
     WRITE( PRTFile, * ) '           ==> ', ray2%x
     WRITE( PRTFile, * ) '           ==> ', x2_o_out
  END IF
  !write( *, * ) 'final coord ', x2_o, ray2%x, w0, w1

  hw0 = h * w0
  hw1 = h * w1
  ray2%t   = ray0%t   - hw0 * gradc0 / csq0      - hw1 * gradc1 / csq1
  ray2%p   = ray0%p   - hw0 * cnn0_csq0 * ray0%q - hw1 * cnn1_csq1 * ray1%q
  ray2%q   = ray0%q   + hw0 * c0        * ray0%p + hw1 * c1        * ray1%p
  ray2%tau = ray0%tau + hw0 / CMPLX( c0, cimag0, KIND=8 ) + hw1 / CMPLX( c1, cimag1, KIND=8 )

  ray2%Amp       = ray0%Amp
  ray2%Phase     = ray0%Phase
  ray2%NumTopBnc = ray0%NumTopBnc
  ray2%NumBotBnc = ray0%NumBotBnc

  ! If we crossed an interface, apply jump condition

  CALL EvaluateSSP2D( ray2%x, ray2%t, c2, cimag2, gradc2, crr2, crz2, czz2, rho, xs_3D, tradial, freq )
  ray2%c = c2

  !!! this needs modifying like the full 3D version to handle jumps in the x-y direction
  IF ( iSegz /= iSegz0 ) THEN
     gradcjump =  gradc2 - gradc0  ! this is precise only for c-linear layers
     ray2n     = [ -ray2%t( 2 ), ray2%t( 1 ) ]

     cnjump    = DOT_PRODUCT( gradcjump, ray2n  )
     csjump    = DOT_PRODUCT( gradcjump, ray2%t )

     RM        = ray2%t( 1 ) / ray2%t( 2 )
     RN        = RM * ( 2 * cnjump - RM * csjump ) / c2
     RN        = -RN
     ray2%p    = ray2%p + ray2%q * RN

  END IF

END SUBROUTINE Step2D