SUBROUTINE Step2D( ray0, ray2, topRefl, botRefl, Topx, Topn, Botx, Botn )
! 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 BdryMod
TYPE( ray2DPt ), INTENT( INOUT ) :: ray0, ray2
TYPE( ray2DPt ) :: ray1
REAL (KIND=8 ), INTENT( IN ) :: Topx( 2 ), Topn( 2 ), Botx( 2 ), Botn( 2 )
LOGICAL, INTENT( OUT ) :: topRefl, botRefl
INTEGER :: iSegz0, iSegr0
REAL (KIND=8 ) :: gradc0( 2 ), gradc1( 2 ), gradc2( 2 ), &
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 ), &
h, halfh, ray2n( 2 ), RM, RN, gradcjump( 2 ), cnjump, csjump, w0, w1, rho
REAL (KIND=8 ) :: urayt2( 2 ), unitdt( 2 ), unitdp( 2 ), unitdq( 2 )
COMPLEX (KIND=8 ) :: unitdtau
IF ( STEP_DEBUGGING ) THEN
WRITE( PRTFile, * )
WRITE( PRTFile, * ) 'ray0 x t p q tau amp', ray0%x, ray0%t, ray0%p, ray0%q, ray0%tau, ray0%Amp
WRITE( PRTFile, * ) 'iSegr iSegz', iSegr, iSegz
END IF
! IF ( ray0%x( 1 ) > 10.0 ) THEN
! STOP 'Enough'
! 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 EvaluateSSP( ray0%x, ray0%t, c0, cimag0, gradc0, crr0, crz0, czz0, rho, freq, 'TAB' )
csq0 = c0 * c0
cnn0_csq0 = crr0 * ray0%t( 2 )**2 - 2.0 * crz0 * ray0%t( 1 ) * ray0%t( 2 ) + czz0 * ray0%t( 1 )**2
iSegz0 = iSegz ! make note of current layer
iSegr0 = iSegr
h = Beam%deltas ! initially set the step h, to the basic one, deltas
urayt0 = c0 * ray0%t ! unit tangent
CALL ReduceStep2D( ray0%x, urayt0, iSegz0, iSegr0, Topx, Topn, Botx, Botn, h ) ! reduce h to land on boundary
! WRITE( PRTFile, * ) 'out h, urayt0', h, urayt0
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
! WRITE( PRTFile, * ) 'ray1 x t p q', ray1%x, ray1%t, ray1%p, ray1%q
! *** Phase 2
CALL EvaluateSSP( ray1%x, ray1%t, c1, cimag1, gradc1, crr1, crz1, czz1, rho, freq, 'TAB' )
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
CALL ReduceStep2D( ray0%x, urayt1, iSegz0, iSegr0, Topx, Topn, Botx, Botn, h ) ! reduce h to land on boundary
! WRITE( PRTFile, * ) 'urayt1', urayt1, 'out h 2', h
! use blend of f' based on proportion of a full step used.
w1 = h / ( 2.0d0 * halfh )
w0 = 1.0d0 - w1
! WRITE( PRTFile, * ) 'w1 w0', w1, w0
urayt2 = w0 * urayt0 + w1 * urayt1
unitdt = -w0 * gradc0 / csq0 - w1 * gradc1 / csq1
unitdp = -w0 * cnn0_csq0 * ray0%q - w1 * cnn1_csq1 * ray1%q
unitdq = w0 * c0 * ray0%p + w1 * c1 * ray1%p
unitdtau = w0 / CMPLX( c0, cimag0, KIND=8 ) + w1 / CMPLX( c1, cimag1, KIND=8 )
! 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 StepToBdry2D(ray0%x, ray2%x, urayt2, h, topRefl, botRefl, &
iSegz0, iSegr0, Topx, Topn, Botx, Botn)
ray2%t = ray0%t + h * unitdt
ray2%p = ray0%p + h * unitdp
ray2%q = ray0%q + h * unitdq
ray2%tau = ray0%tau + h * unitdtau
! WRITE( PRTFile, * ) 'ray2 x t p q tau', ray2%x, ray2%t, ray2%p, ray2%q, ray2%tau
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 EvaluateSSP( ray2%x, ray2%t, c2, cimag2, gradc2, crr2, crz2, czz2, rho, freq, 'TAB' )
ray2%c = c2
IF ( iSegz /= iSegz0 .OR. iSegr /= iSegr0 ) THEN
gradcjump = gradc2 - gradc0
ray2n = [ -ray2%t( 2 ), ray2%t( 1 ) ] ! ray normal
cnjump = DOT_PRODUCT( gradcjump, ray2n )
csjump = DOT_PRODUCT( gradcjump, ray2%t )
IF ( iSegz /= iSegz0 ) THEN ! crossing in depth
RM = +ray2%t( 1 ) / ray2%t( 2 ) ! this is tan( alpha ) where alpha is the angle of incidence
ELSE ! crossing in range
RM = -ray2%t( 2 ) / ray2%t( 1 ) ! this is tan( alpha ) where alpha is the angle of incidence
END IF
RN = RM * ( 2 * cnjump - RM * csjump ) / c2
ray2%p = ray2%p - ray2%q * RN
END IF
END SUBROUTINE Step2D
