InfluenceCervenyCart Subroutine

  SUBROUTINE InfluenceCervenyCart( U, eps, alpha, iBeamWindow2, RadiusMax )
  !! Paraxial (Cerveny-style) beams in Cartesian coordinates

    INTEGER,          INTENT( IN    ) :: IBeamWindow2
    REAL    (KIND=8), INTENT( IN    ) :: alpha, RadiusMax                ! take-off angle
    COMPLEX,          INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )  ! complex pressure field
    COMPLEX (KIND=8), INTENT( IN    ) :: eps                          ! LP: EPSILON is an intrinsic
    INTEGER          :: KMAHV( MaxN ), KMAH, irA, irB, Image
    REAL    (KIND=8), SAVE :: Polarity = 1
    REAL    (KIND=8) :: x( 2 ), rayt( 2 ), rayn( 2 ), Tr, Tz, zr, &
         c, cimag, cs, cn, csq, gradc( 2 ), crr, crz, czz, rho, deltaz
    COMPLEX (KIND=8) :: pVB( MaxN ), qVB( MaxN ), q, epsV( MaxN ), contri, gammaV( MaxN ), gamma, const
    COMPLEX (KIND=8) :: tau

    SELECT CASE ( Beam%RunType( 1 : 1 ) )
    CASE ( 'C', 'I', 'S' )   ! TL
    CASE DEFAULT
       WRITE( PRTFile, * ) 'Cerveny influence does not support eigenrays or arrivals'
       CALL ERROUT( 'InfluenceCervenyRayCen', 'Cerveny influence does not support eigenrays or arrivals' )
    END SELECT

    ! need to add logic related to NRz_per_range

    ! During reflection imag(q) is constant and adjacent normals cannot bracket a segment of the TL
    ! line, so no special treatment is necessary

    IF ( Beam%Type( 2 : 2 ) == 'C' ) THEN
       epsV( 1 : Beam%Nsteps ) = i * ABS( ray2D( 1 : Beam%Nsteps )%q( 1 ) / ray2D( 1 : Beam%Nsteps )%q( 2 ) )
    ELSE
       epsV( 1 : Beam%Nsteps ) = eps
    END IF

    pVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%p( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%p( 2 )
    qVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%q( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%q( 2 )

    IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) )  ! point source

    ! Form gamma and KMAH index
    ! Treatment of KMAH index is incorrect for 'Cerveny' style beam width BeamType

    Stepping0: DO iS = 1, Beam%Nsteps

       rayt = ray2D( iS )%c * ray2D( iS )%t   ! unit tangent
       rayn = [ rayt( 2 ), -rayt( 1 ) ]       ! unit normal

       CALL EvaluateSSP( ray2D( iS )%x, ray2D( iS )%t, c, cimag, gradc, crr, crz, czz, rho, Freq, 'TAB' )

       csq = c * c
       cS  = DOT_PRODUCT( gradc, rayt )
       cN  = DOT_PRODUCT( gradc, rayn )

       Tr  = rayt(  1 )
       Tz  = rayt(  2 )

       gammaV( iS ) = 0.0
       IF ( qVB( iS ) /= 0.0 ) THEN
          gammaV( iS ) = 0.5 * ( pVB( iS ) / qVB( iS ) * Tr**2 + 2.0 * cN / csq * Tz * Tr - cS / csq * Tz**2 )
       END IF

       IF ( iS == 1 ) THEN
          KMAHV( 1 ) = 1
       ELSE
          KMAHV( iS ) = KMAHV( iS - 1 )
          CALL BranchCut( qVB( iS - 1 ), qVB( iS ), Beam%Type, KMAHV( iS ) )
       END IF
    END DO Stepping0

    Stepping: DO iS = 3, Beam%Nsteps
       ! LP: BUG: Assumes rays may never travel left.
       IF ( ray2D( iS     )%x( 1 ) > Pos%Rr( Pos%NRr ) ) RETURN
       rA = ray2D( iS - 1 )%x( 1 )
       rB = ray2D( iS     )%x( 1 )
       IF ( ABS( rB - rA ) < 1.0D3 * SPACING( rB ) ) CYCLE Stepping   ! don't process duplicate points

       irA = RToIR( rA )
       irB = RToIR( rB )

       IF ( irA >= irB ) CYCLE Stepping

       RcvrRanges: DO ir = irA + 1, irB

          W     = ( Pos%Rr( ir ) - rA ) / ( rB - rA )

          x     = ray2D(  iS - 1 )%x    + W * ( ray2D(  iS )%x   -  ray2D(  iS - 1 )%x   )
          rayt  = ray2D(  iS - 1 )%t    + W * ( ray2D(  iS )%t   -  ray2D(  iS - 1 )%t   )
          c     = ray2D(  iS - 1 )%c    + W * ( ray2D(  iS )%c   -  ray2D(  iS - 1 )%c   )
          q     = qVB(    iS - 1 )      + W * ( qVB(    iS )     -  qVB(    iS - 1 )     )
          tau   = ray2D(  iS - 1 )%tau  + W * ( ray2D(  iS )%tau -  ray2D(  iS - 1 )%tau )
          gamma = gammaV( iS - 1 )      + W * ( gammaV( iS )     -  gammaV( iS - 1 )     )

          IF ( AIMAG( gamma ) > 0 ) THEN
             WRITE( PRTFile, * ) 'Unbounded beam'
             WRITE( PRTFile, * ) gammaV( iS - 1 ), gammaV( iS ), gamma
             CYCLE RcvrRanges
          END IF

          const = Ratio1 * SQRT( c * ABS( epsV( iS - 1 ) ) / q )

          ! Get correct branch of SQRT
          KMAH = KMAHV( iS - 1 )
          CALL BranchCut( qVB( iS - 1 ), q, Beam%Type, KMAH )
          IF ( KMAH < 0 ) const = -const

          RcvrDepths: DO iz = 1, NRz_per_range
             zR = Pos%Rz( iz )

             contri = 0.0
             ImageLoop: DO Image = 1, Beam%Nimage
                SELECT CASE ( Image )
                CASE ( 1 )   ! True beam
                   deltaz = zR - x( 2 )
                   Polarity = +1.0D0
                CASE ( 2 )   ! Surface reflected beam
                   deltaz = -zR + 2.0 * Bdry%Top%HS%Depth - x( 2 )
                   Polarity = -1.0D0
                CASE ( 3 )   ! Bottom  reflected beam
                   deltaz = -zR + 2.0 * Bdry%Bot%HS%Depth - x( 2 )
                   Polarity = +1.0D0   ! assumes rigid bottom
                END SELECT

                IF ( omega * AIMAG( gamma ) * deltaz ** 2 < iBeamWindow2 ) &
                     contri =  contri + Polarity * ray2D( iS )%Amp * Hermite( deltaz, RadiusMax, 2.0 * RadiusMax ) * &
                     EXP( -i * ( omega * ( tau + rayt( 2 ) * deltaz + gamma * deltaz**2 ) - ray2D( iS )%Phase ) )
             END DO ImageLoop

             ! contribution to field
             SELECT CASE( Beam%RunType( 1 : 1 ) )
             CASE ( 'C' )        ! coherent
                contri = const * contri
             CASE ( 'I', 'S' )   ! incoherent or semi-coherent
                contri = ABS( const * contri ) ** 2
             END SELECT
             U( iz, ir ) = U( iz, ir ) + CMPLX( contri )
          END DO RcvrDepths
       END DO RcvrRanges
    END DO Stepping

  END SUBROUTINE InfluenceCervenyCart