!! Beam influence computation for pressure field calculation

MODULE Influence
  !! Computes beam contributions to complex pressure fields using various beam weighting approaches

  ! Compute the beam influence, i.e. the contribution of a single beam to the complex pressure
  ! mbp 12/2018, based on much older subroutines

  USE bellhopMod
  USE SourceReceiverPositions
  USE ArrMod
  USE sspMod   ! used to construct image beams in the Cerveny style beam routines
  USE WriteRay
  USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT

  IMPLICIT NONE
  PUBLIC

  INTEGER,          PRIVATE :: iz, ir, iS
  REAL    (KIND=8), PRIVATE :: Ratio1 = 1.0D0   ! scale factor for a line source
  REAL    (KIND=8), PRIVATE :: W, s, n, Amp, phase, const, phaseInt, q0, q, qold, RcvrDeclAngle, rA, rB
  COMPLEX (KIND=8), PRIVATE :: delay

CONTAINS
  SUBROUTINE InfluenceCervenyRayCen( U, eps, alpha, iBeamWindow2, RadiusMax )
  !! Paraxial (Cerveny-style) beams in ray-centered 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          :: ir1, ir2, KMAHV( MaxN ), KMAH, image
    REAL    (KIND=8) :: nA, nB, nSq, c, zr, Polarity
    REAL    (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps )   ! ray normal
    COMPLEX (KIND=8) :: pVB( MaxN ), qVB( MaxN ), q, epsV( MaxN ), contri, gammaV( MaxN ), gamma, P_n, P_s
    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 )
    gammaV( 1 : Beam%Nsteps ) = pVB(   1 : Beam%Nsteps ) / qVB( 1 : Beam%Nsteps )

    ! pre-calculate ray normal based on tangent with c(s) scaling
    znV = -ray2D( 1 : Beam%Nsteps )%t( 1 ) * ray2D( 1 : Beam%Nsteps )%c
    rnV =  ray2D( 1 : Beam%Nsteps )%t( 2 ) * ray2D( 1 : Beam%Nsteps )%c

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

    ! compute KMAH index
    ! Following is incorrect for 'Cerveny'-style beamwidth (narrow as possible)
    KMAHV(  1 ) = 1

    DO iS = 2, Beam%Nsteps
       KMAHV(  iS ) = KMAHV( iS - 1 )
       CALL BranchCut( qVB( iS - 1 ), qVB( iS ), Beam%Type, KMAHV( iS ) )
    END DO

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

       Images: DO image = 1, Beam%Nimage

          ! LP: Previous code did rnV = -rnV for image 2 and 3. When Nimage = 2,
          ! this means rnV changes sign every step, which can't possibly be
          ! correct. This was fixed by mbp in InfluenceCervenyCart.
          Polarity = 1.0D0
          IF ( image == 2 ) Polarity = -1.0D0

          !!! This logic means that the first step along the ray is skipped
          !!! which is a problem if deltas is very large, e.g. isospeed problems
          !!! I fixed this in InfluenceGeoHatRayCen
          ir1 = HUGE( ir1 )

          Stepping: DO iS = 2, Beam%Nsteps

             ! Compute ray-centered coordinates, (znV, rnV)

             ! If normal parallel to TL-line, skip to next step on ray
             ! LP: Changed from TINY( znV( iS ) ), see README.md.
             IF ( ABS( znV( iS ) ) < EPSILON( znV( iS ) ) ) CYCLE Stepping

             SELECT CASE ( image )     ! Images of beams
             CASE ( 1 )                ! True beam
                nB  = ( zR -                             ray2D( iS )%x( 2 )   ) / znV( iS )
             CASE ( 2 )                ! Surface-reflected beam
                nB  = ( zR - ( 2.0 * Bdry%Top%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
             CASE ( 3 )                ! Bottom-reflected beam
                nB  = ( zR - ( 2.0 * Bdry%Bot%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
             CASE DEFAULT
                WRITE( PRTFile, * ) 'Beam%Nimage = ', Beam%Nimage
                CALL ERROUT( 'InfluenceCervenyRayCen', 'Nimage must be 1, 2, or 3' )
             END SELECT

             rB  = ray2D( iS )%x( 1 ) + nB * rnV( iS ) * Polarity
             ir2 = RToIR( rB )

             ! detect and skip duplicate points (happens at boundary reflection)
             IF ( ir1 >= ir2 .OR. &
                & ABS( ray2D( iS )%x( 1 ) - ray2D( iS - 1 )%x( 1 ) ) < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) THEN
                rA  = rB
                nA  = nB
                ir1 = ir2
                CYCLE Stepping
             END IF

             RcvrRanges: DO ir = ir1 + 1, ir2    ! Compute influence for each rcvr
                W     = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
                q     =    qVB( iS - 1 ) + W * (    qVB( iS ) -    qVB( iS - 1 ) )
                gamma = gammaV( iS - 1 ) + W * ( gammaV( iS ) - gammaV( iS - 1 ) )
                n     = nA + W * ( nB - nA )
                nSq   = n * n
                IF ( AIMAG( gamma ) > 0 ) THEN
                   WRITE( PRTFile, * ) 'Unbounded beam'
                   CYCLE RcvrRanges
                END IF

                IF ( -0.5 * omega * AIMAG( gamma ) * nSq < iBeamWindow2 ) THEN   ! Within beam window?
                   c      = ray2D( iS - 1 )%c
                   tau    = ray2D( iS - 1 )%tau + W * ( ray2D( iS )%tau - ray2D( iS - 1 )%tau )
                   contri = ratio1 * ray2D( iS )%Amp * SQRT( c * ABS( epsV( iS ) ) / q ) * &
                        EXP( -i * ( omega * ( tau + 0.5 * gamma * nSq ) - ray2D( iS )%phase ) )

                   SELECT CASE ( Beam%Component )
                   CASE ( 'P' )   ! pressure
                   CASE ( 'V' )   ! vertical component
                      P_n    = -i * omega * gamma * n * contri
                      P_s    = -i * omega / c         * contri
                      contri = c * DOT_PRODUCT( [ P_n, P_s ], ray2D( iS )%t )
                   CASE ( 'H' )   ! horizontal component
                      P_n    = -i * omega * gamma * n * contri
                      P_s    = -i * omega / c         * contri
                      contri = c * ( -P_n * ray2D( iS )%t( 2 ) + P_s * ray2D( iS )%t( 1 ) )
                   CASE DEFAULT
                      WRITE( ERROR_UNIT, * ) 'InfluenceCervenyRayCen: Unknown component: ', Beam%Component
                      CALL ERROUT( 'InfluenceCervenyRayCen', 'Unknown beam component' )
                   END SELECT

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

                   IF ( KMAH  < 0  ) contri = -contri
                   contri = Polarity * contri

                   SELECT CASE ( Beam%RunType( 1 : 1 ) )
                   CASE ( 'I', 'S' )   ! Incoherent or Semi-coherent TL
                      contri = ABS( contri ) ** 2
                   CASE DEFAULT
                      WRITE( ERROR_UNIT, * ) 'InfluenceCervenyRayCen: Unknown RunType: ', Beam%RunType( 1 : 1 )
                      CALL ERROUT( 'InfluenceCervenyRayCen', 'Unknown RunType' )
                   END SELECT

                   U( iz, ir ) = U( iz, ir ) + CMPLX( Hermite( n, RadiusMax, 2 * RadiusMax ) * contri )
                END IF
             END DO RcvrRanges
             rA  = rB
             nA  = nB
             ir1 = ir2
          END DO Stepping
       END DO Images
    END DO RcvrDepths

  END SUBROUTINE InfluenceCervenyRayCen

  ! **********************************************************************!

  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
               CASE DEFAULT
                  WRITE( ERROR_UNIT, * ) 'InfluenceCervenyCart: Unknown image index: ', Image
                  CALL ERROUT( 'InfluenceCervenyCart', 'Unknown image index' )
               END SELECT

                IF ( omega * AIMAG( gamma ) * deltaz ** 2 < iBeamWindow2 ) THEN
                     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 IF
             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
            CASE DEFAULT
               WRITE( ERROR_UNIT, * ) 'InfluenceCervenyCart: Unknown RunType: ', Beam%RunType( 1 : 1 )
               CALL ERROUT( 'InfluenceCervenyCart', 'Unknown RunType' )
            END SELECT
             U( iz, ir ) = U( iz, ir ) + CMPLX( contri )
          END DO RcvrDepths
       END DO RcvrRanges
    END DO Stepping

  END SUBROUTINE InfluenceCervenyCart

  ! **********************************************************************!

  SUBROUTINE InfluenceGeoHatRayCen( U, alpha, dalpha )
  !! Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates

    REAL (KIND=8), INTENT( IN    ) :: alpha, dalpha                 ! take-off angle
    COMPLEX,       INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )   ! complex pressure field
    INTEGER          :: irA, irB, II
    REAL    (KIND=8) :: nA, nB, zr, L, dq( Beam%Nsteps - 1 )
    REAL    (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps ),  RcvrDeclAngleV ( Beam%Nsteps )
    COMPLEX (KIND=8) :: dtau( Beam%Nsteps - 1 )
    REAL    (KIND=8) :: KMAHphase( Beam%Nsteps )

    ! need to add logic related to NRz_per_range

    ! LP: See discussion of this change in the readme.
    qOld = ray2D( 1 )%q( 1 )
    phase = 0
    KMAHphase( 1 ) = 0
    DO is = 2, Beam%Nsteps
       q = ray2D( is )%q( 1 )
       CALL IncPhaseIfCaustic( .TRUE. )
       qOld = q
       KMAHphase( is ) = phase
    END DO

    q0           = ray2D( 1 )%c / Dalpha   ! Reference for J = q0 / q
    SrcDeclAngle = RadDeg * alpha          ! take-off angle in degrees

    dq   = ray2D( 2 : Beam%Nsteps )%q( 1 ) - ray2D( 1 : Beam%Nsteps - 1 )%q( 1 )
    dtau = ray2D( 2 : Beam%Nsteps )%tau    - ray2D( 1 : Beam%Nsteps - 1 )%tau

    ! pre-calculate ray normal based on tangent with c(s) scaling
    znV = -ray2D( 1 : Beam%Nsteps )%t( 1 ) * ray2D( 1 : Beam%Nsteps )%c
    rnV =  ray2D( 1 : Beam%Nsteps )%t( 2 ) * ray2D( 1 : Beam%Nsteps )%c

    RcvrDeclAngleV( 1 : Beam%Nsteps ) = RadDeg * ATAN2( ray2D( 1 : Beam%Nsteps )%t( 2 ), ray2D( 1 : Beam%Nsteps )%t( 1 ) )

    ! 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%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) )  ! point source

    ray2D( 1 : Beam%Nsteps )%Amp = Ratio1 * SQRT( ray2D( 1 : Beam%Nsteps )%c ) * ray2D( 1 : Beam%Nsteps )%Amp
        ! pre-apply some scaling

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

       IF ( ABS( znV( 1 ) ) < 1D-6 ) THEN   ! normal parallel to horizontal receiver line
          nA  = 1D10
          rA  = 1D10
          irA = 1
       ELSE
          nA  = ( zR - ray2D( 1 )%x( 2 )   ) / znV( 1 )
          rA  = ray2D( 1 )%x( 1 ) + nA * rnV( 1 )
          irA = RToIR( rA )
       END IF

       Stepping: DO iS = 2, Beam%Nsteps

          ! Compute ray-centered coordinates, (znV, rnV)

          IF ( ABS( znV( iS ) ) < 1D-10 ) CYCLE Stepping   ! If normal parallel to TL-line, skip to next step on ray
          nB  = ( zR - ray2D( iS )%x( 2 )   ) / znV( iS )
          rB  = ray2D( iS )%x( 1 ) + nB * rnV( iS )
          irB = RToIR( rB )

          ! detect and skip duplicate points (happens at boundary reflection)
          IF ( ABS( ray2D( iS )%x( 1 ) - ray2D( iS - 1 )%x( 1 ) ) < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) .OR. irA == irB ) THEN
             rA  = rB
             nA  = nB
             irA = irB
             CYCLE Stepping
          END IF

          RcvrDeclAngle = RcvrDeclAngleV( iS )

          ! *** Compute contributions to bracketed receivers ***

          II = 0
          IF ( irB <= irA ) II = 1   ! going backwards in range

          RcvrRanges: DO ir = irA + 1 - II, irB + II, SIGN( 1, irB - irA )  ! Compute influence for each rcvr
             W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
             n = ABS( nA                + W * ( nB - nA ) )
             q = ray2D( iS - 1 )%q( 1 ) + W * dq( iS - 1 )     ! interpolated amplitude
             L = ABS( q ) / q0   ! beam radius

             IF ( n < L ) THEN   ! in beamwindow?
                delay    = ray2D( iS - 1 )%tau + W * dtau( iS - 1 )   ! interpolated delay
                const    = ray2D( iS )%Amp / SQRT( ABS( q ) )
                W        = ( L - n ) / L   ! hat function: 1 on center, 0 on edge
                Amp      = const * W

                phase = KMAHphase( iS - 1 )
                qOld = ray2D( is - 1 )%q( 1 )
                CALL FinalPhase( .FALSE. )

                CALL ApplyContribution( U( iz, ir ) )
             END IF
          END DO RcvrRanges
          rA  = rB
          nA  = nB
          irA = irB
       END DO Stepping
    END DO RcvrDepths

  END SUBROUTINE InfluenceGeoHatRayCen

  ! **********************************************************************!

  SUBROUTINE InfluenceGeoHatCart( U, alpha, Dalpha )
  !! Geometric, hat-shaped beams in Cartesisan coordinates

    REAL (KIND=8), INTENT( IN    ) :: alpha, dalpha                 ! take-off angle, angular spacing
    COMPLEX,       INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )   ! complex pressure field
    INTEGER          :: irT( 1 ), irTT
    REAL    (KIND=8) :: x_ray( 2 ), rayt( 2 ), rayn( 2 ), x_rcvr( 2, NRz_per_range ), rLen, RadiusMax, zMin, zMax, dqds
    COMPLEX (KIND=8) :: dtauds

    q0           = ray2D( 1 )%c / Dalpha   ! Reference for J = q0 / q
    SrcDeclAngle = RadDeg * alpha          ! take-off angle in degrees
    phase        = 0.0
    qOld         = ray2D( 1 )%q( 1 )       ! used to track KMAH index
    rA           = ray2D( 1 )%x( 1 )       ! range at start of ray

    ! what if never satisfied?
    ! what if there is a single receiver (ir = 0 possible)
    irT = MINLOC( Pos%Rr( 1 : Pos%NRr ), MASK = Pos%Rr( 1 : Pos%NRr ) > rA )   ! find index of first receiver to the right of rA
    ir  = irT( 1 )
    IF ( ray2D( 1 )%t( 1 ) < 0.0d0 .AND. ir > 1 ) ir = ir - 1
    ! if ray is left-traveling, get the first receiver to the left of rA

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

    Stepping: DO iS = 2, Beam%Nsteps
       rB     = ray2D( iS     )%x( 1 )
       x_ray  = ray2D( iS - 1 )%x

       ! compute normalized tangent (compute it because we need to measure the step length)
       rayt = ray2D( iS )%x - ray2D( iS - 1 )%x
       rlen = NORM2( rayt )
       IF ( rlen < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) CYCLE Stepping
           ! if duplicate point in ray, skip to next step along the ray
       rayt = rayt / rlen                    ! unit tangent to ray
       rayn = [ -rayt( 2 ), rayt( 1 ) ]      ! unit normal  to ray
       RcvrDeclAngle = RadDeg * ATAN2( rayt( 2 ), rayt( 1 ) )

       dqds   = ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 )
       dtauds = ray2D( iS )%tau    - ray2D( iS - 1 )%tau

       q  = ray2D( iS - 1 )%q( 1 )
       CALL IncPhaseIfCaustic( .TRUE. )
       qold = q

       RadiusMax = MAX( ABS( ray2D( iS - 1 )%q( 1 ) ), ABS( ray2D( iS )%q( 1 ) ) ) / q0 / ABS( rayt( 1 ) )
           ! beam radius projected onto vertical line

       ! depth limits of beam
       IF ( ABS( rayt( 1 ) ) > 0.5 ) THEN   ! shallow angle ray
          zmin   = min( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) - RadiusMax
          zmax   = max( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) + RadiusMax
       ELSE                                 ! steep angle ray
          zmin = -HUGE( zmin )
          zmax = +HUGE( zmax )
       END IF

       ! compute beam influence for this segment of the ray
       RcvrRanges: DO
          ! is Rr( ir ) contained in [ rA, rB )? Then compute beam influence
          IF ( Pos%Rr( ir ) >= MIN( rA, rB ) .AND. Pos%Rr( ir ) < MAX( rA, rB ) ) THEN

             x_rcvr( 1, 1 : NRz_per_range ) = Pos%Rr( ir )
             IF ( Beam%RunType( 5 : 5 ) == 'I' ) THEN
                x_rcvr( 2, 1 ) = Pos%Rz( ir )                  ! irregular   grid
             ELSE
                x_rcvr( 2, 1 : NRz_per_range ) = Pos%Rz( 1 : NRz_per_range )   ! rectilinear grid
             END IF

             RcvrDepths: DO iz = 1, NRz_per_range
                IF ( x_rcvr( 2, iz ) < zmin .OR. x_rcvr( 2, iz ) > zmax ) CYCLE RcvrDepths

                s         =      DOT_PRODUCT( x_rcvr( :, iz ) - x_ray, rayt ) / rlen ! proportional distance along ray
                n         = ABS( DOT_PRODUCT( x_rcvr( :, iz ) - x_ray, rayn ) )      ! normal distance to ray
                q         = ray2D( iS - 1 )%q( 1 ) + s * dqds               ! interpolated amplitude
                RadiusMax = ABS( q / q0 )                                   ! beam radius

                IF ( n < RadiusMax ) THEN
                   delay    = ray2D( iS - 1 )%tau + s * dtauds              ! interpolated delay
                   const    = Ratio1 * SQRT( ray2D( iS )%c / ABS( q ) ) * ray2D( iS )%Amp
                   W        = ( RadiusMax - n ) / RadiusMax   ! hat function: 1 on center, 0 on edge
                   Amp      = const * W
                   CALL FinalPhase( .FALSE. )

                   CALL ApplyContribution( U( iz, ir ) )
                END IF
             END DO RcvrDepths
          END IF

          ! bump receiver index, ir, towards rB
          IF ( Pos%Rr( ir ) < rB ) THEN
             IF ( ir >= Pos%NRr        ) EXIT RcvrRanges  ! go to next step on ray
             irTT = ir + 1                                ! bump right
             IF ( Pos%Rr( irTT ) >= rB ) EXIT RcvrRanges
          ELSE
             IF ( ir <= 1              ) EXIT RcvrRanges  ! go to next step on ray
             irTT = ir - 1                                ! bump left
             IF ( Pos%Rr( irTT ) <= rB ) EXIT RcvrRanges
          END IF
          ir = irTT
       END DO RcvrRanges

       rA = rB
    END DO Stepping

  END SUBROUTINE InfluenceGeoHatCart

  ! **********************************************************************!

  SUBROUTINE InfluenceGeoGaussianCart( U, alpha, Dalpha )
  !! Geometric, Gaussian beams in Cartesian coordintes

    INTEGER,       PARAMETER       :: BeamWindow = 4               ! beam window: kills beams outside e**(-0.5 * ibwin**2 )
    REAL (KIND=8), INTENT( IN    ) :: alpha, dalpha                ! take-off angle, angular spacing
    COMPLEX,       INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )  ! complex pressure field
    INTEGER          :: irT( 1 ), irTT
    REAL    (KIND=8) :: x_ray( 2 ), rayt( 2 ), rayn( 2 ), x_rcvr( 2 ), rLen, RadiusMax, zMin, zMax, sigma, lambda, A, dqds
    COMPLEX (KIND=8) :: dtauds

    q0           = ray2D( 1 )%c / Dalpha   ! Reference for J = q0 / q
    SrcDeclAngle = RadDeg * alpha          ! take-off angle in degrees
    phase        = 0
    qOld         = ray2D( 1 )%q( 1 )       ! used to track KMAH index
    rA           = ray2D( 1 )%x( 1 )       ! range at start of ray

    ! what if never satisfied?
    ! what if there is a single receiver (ir = 0 possible)

    irT = MINLOC( Pos%Rr( 1 : Pos%NRr ), MASK = Pos%Rr( 1 : Pos%NRr ) > rA )
    ! find index of first receiver to the right of rA
    ir  = irT( 1 )

    IF ( ray2D( 1 )%t( 1 ) < 0.0d0 .AND. ir > 1 ) ir = ir - 1
    ! if ray is left-traveling, get the first receiver to the left of rA

    ! sqrt( 2 * pi ) represents a sum of Gaussians in free space
    IF ( Beam%RunType( 4 : 4 ) == 'R' ) THEN
       Ratio1 = SQRT( ABS( COS( alpha ) ) ) / SQRT( 2. * pi )   ! point source
    ELSE
       Ratio1 = 1 / SQRT( 2. * pi )                             ! line  source
    END IF

    Stepping: DO iS = 2, Beam%Nsteps

       rB    = ray2D( iS     )%x( 1 )
       x_ray = ray2D( iS - 1 )%x

       ! compute normalized tangent (compute it because we need to measure the step length)
       rayt = ray2D( iS )%x - ray2D( iS - 1 )%x
       rlen = NORM2( rayt )
       IF ( rlen < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) CYCLE Stepping
           ! if duplicate point in ray, skip to next step along the ray
       rayt = rayt / rlen
       rayn = [ -rayt( 2 ), rayt( 1 ) ]      ! unit normal to ray
       RcvrDeclAngle = RadDeg * ATAN2( rayt( 2 ), rayt( 1 ) )

       dqds   = ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 )
       dtauds = ray2D( iS )%tau    - ray2D( iS - 1 )%tau

       q  = ray2D( iS - 1 )%q( 1 )
       CALL IncPhaseIfCaustic( .TRUE. )
       qold = q

       ! calculate beam width
       lambda    = ray2D( iS - 1 )%c / freq
       sigma     = MAX( ABS( ray2D( iS - 1 )%q( 1 ) ), ABS( ray2D( iS )%q( 1 ) ) ) / q0 / ABS( rayt( 1 ) )
                   ! beam radius projected onto vertical line
       sigma     = MAX( sigma, MIN( 0.2 * freq * REAL( ray2D( iS )%tau ), pi * lambda ) )
       RadiusMax = BeamWindow * sigma

       ! depth limits of beam
       ! LP: For rays shot at exactly 60 degrees, they will hit this edge case.
       ! This is a sharp edge--the handling on each side of this edge may be
       ! significantly different. So, moved the edge away from the round number.
       IF ( ABS( rayt( 1 ) ) > 0.50001 ) THEN   ! shallow angle ray
          zmin   = min( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) - RadiusMax
          zmax   = max( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) + RadiusMax
       ELSE                                 ! steep angle ray
          zmin = -HUGE( zmin )
          zmax = +HUGE( zmax )
       END IF

       ! compute beam influence for this segment of the ray
       RcvrRanges: DO
          ! WRITE( PRTFile, * ) 'iS', iS-2, 'ir', ir-1

          ! is Rr( ir ) contained in [ rA, rB )? Then compute beam influence
          IF ( Pos%Rr( ir ) >= MIN( rA, rB ) .AND. Pos%Rr( ir ) < MAX( rA, rB ) ) THEN
             ! WRITE( PRTFile, * ) '  rA', rA, 'rB', rB

             RcvrDepths: DO iz = 1, NRz_per_range
                IF ( Beam%RunType( 5 : 5 ) == 'I' ) THEN
                   x_rcvr = [ Pos%Rr( ir ), Pos%Rz( ir ) ]   ! irregular   grid
                ELSE
                   x_rcvr = [ Pos%Rr( ir ), Pos%Rz( iz ) ]   ! rectilinear grid
                END IF
                IF ( x_rcvr( 2 ) < zmin .OR. x_rcvr( 2 ) > zmax ) CYCLE RcvrDepths

                s      =      DOT_PRODUCT( x_rcvr - x_ray, rayt ) / rlen  ! proportional distance along ray
                n      = ABS( DOT_PRODUCT( x_rcvr - x_ray, rayn ) )       ! normal distance to ray
                q      = ray2D( iS - 1 )%q( 1 ) + s * dqds                ! interpolated amplitude
                sigma  = ABS( q / q0 )                                    ! beam radius
                sigma  = MAX( sigma, MIN( 0.2 * freq * REAL( ray2D( iS )%tau ), pi * lambda ) )  ! min pi * lambda, unless near

                IF ( n < BeamWindow * sigma ) THEN   ! Within beam window?
                   ! IF ( ir-1 >= 0 .AND. ir-1 <= 2 ) THEN
                   !    WRITE( PRTFile, * ) '    iz n sigma', iz-1, n, sigma
                   ! END IF
                   A        = ABS( q0 / q )
                   delay    = ray2D( iS - 1 )%tau + s * dtauds     ! interpolated delay
                   const    = Ratio1 * SQRT( ray2D( iS )%c / ABS( q ) ) * ray2D( iS )%Amp
                   W        = EXP( -0.5 * ( n / sigma ) ** 2 ) / ( sigma * A )   ! Gaussian decay
                   Amp      = const * W
                   CALL FinalPhase( .TRUE. )

                   CALL ApplyContribution( U( iz, ir ) )
                END IF
             END DO RcvrDepths
          END IF

          ! receiver not bracketed; bump receiver index, ir, towards rB
          IF ( rB > Pos%Rr( ir ) ) THEN
             IF ( ir >= Pos%NRr        ) EXIT RcvrRanges  ! go to next step on ray
             irTT = ir + 1                                ! bump right
             IF ( Pos%Rr( irTT ) >= rB ) EXIT RcvrRanges  ! go to next step on ray
          ELSE
             IF ( ir <= 1              ) EXIT RcvrRanges  ! go to next step on ray
             irTT = ir - 1                                ! bump left
             IF ( Pos%Rr( irTT ) <= rB ) EXIT RcvrRanges  ! go to next step on ray
          END IF
          ir = irTT

       END DO RcvrRanges

       rA = rB
    END DO Stepping

  END SUBROUTINE InfluenceGeoGaussianCart

  ! **********************************************************************!

  SUBROUTINE ApplyContribution( U )
  !! Applies beam contribution to pressure field

    COMPLEX, INTENT( INOUT ) :: U
    COMPLEX ( KIND=4 ) :: dfield

    SELECT CASE( Beam%RunType( 1 : 1 ) )
    CASE ( 'E' )                ! eigenrays
       IF ( Title( 1 :  9 ) == 'BELLHOP- ' ) THEN   ! BELLHOP run
          CALL WriteRay2D( SrcDeclAngle, iS )
       ELSE                                         ! BELLHOP3D run
          CALL WriteRay3D( DegRad * SrcDeclAngle, DegRad * SrcAzimAngle, is )   ! produces no output if NR=1
       END IF
    CASE ( 'A', 'a' )           ! arrivals
       CALL AddArr( omega, iz, ir, Amp, phaseInt, delay, SrcDeclAngle, &
                  & RcvrDeclAngle, ray2D( iS )%NumTopBnc, ray2D( iS )%NumBotBnc )
    CASE ( 'C' )                ! coherent TL
       dfield = CMPLX( Amp * EXP( -i * ( omega * delay - phaseInt ) ) )
       ! WRITE( PRTFile, * ) 'ApplyContribution dfield', dfield
       U = U + dfield
                     ! omega * n * n / ( 2 * ray2d( iS )%c**2 * delay ) ) ) )   ! curvature correction
    CASE DEFAULT                ! incoherent/semicoherent TL
       IF ( Beam%Type( 1 : 1 ) == 'B' ) THEN   ! Gaussian beam
          U = U + SNGL( SQRT( 2. * pi ) * ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
       ELSE
          U = U + SNGL(                   ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
       END IF
    END SELECT

  END SUBROUTINE ApplyContribution

  ! **********************************************************************!

  SUBROUTINE InfluenceSGB( U, alpha, Dalpha, RadiusMax )
  !! Bucker's Simple Gaussian Beams in Cartesian coordinates

    REAL (KIND=8), INTENT( IN    ) :: alpha, dalpha, RadiusMax      ! take-off angle, angular spacing
    COMPLEX,       INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )   ! complex pressure field
    REAL    (KIND=8)  :: x( 2 ), rayt( 2 ), A, beta, cn, CPA, Adeltaz, deltaz, DS, sint, SX1, thet
    COMPLEX (KIND=8)  :: contri, tau

    ! LP: Added ABS to match other influence functions. Without it, this will
    ! crash (sqrt of negative real) for rays shot backwards.
    Ratio1 = SQRT( ABS( COS( alpha ) ) )
    phase  = 0
    qOld   = 1.0
    BETA   = 0.98  ! Beam Factor
    A      = -4.0 * LOG( BETA ) / Dalpha**2
    CN     = Dalpha * SQRT( A / pi )
    rA     = ray2D( 1 )%x( 1 )
    ir     = 1

    Stepping: DO iS = 2, Beam%Nsteps

       RcvrDeclAngle = RadDeg * ATAN2( ray2D( iS )%t( 2 ), ray2D( iS )%t( 1 ) )

       rB = ray2D( iS )%x( 1 )

       ! phase shifts at caustics
       q  = ray2D( iS - 1 )%q( 1 )
       CALL IncPhaseIfCaustic( .FALSE. )
       qold = q

       ! Loop over bracketed receiver ranges
       ! LP: BUG: This way of setting up the loop assumes the ray always travels
       ! towards positive R, which is not true for certain bathymetries (or for
       ! rays simply shot backwards, which previously would also crash during
       ! the setup, see above).
       RcvrRanges: DO WHILE ( ABS( rB - rA ) > 1.0D3 * SPACING( rA ) .AND. rB > Pos%Rr( ir ) )

          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 )
          q     = ray2D( iS - 1 )%q( 1 ) + W * ( ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 ) )
          tau   = ray2D( iS - 1 )%tau    + W * ( ray2D( iS )%tau    - ray2D( iS - 1 )%tau )

          ! BUG: following is incorrect because ray doesn't always use a step of deltas
          ! LP: The do while ignores extremely small steps, but those small steps
          ! still increment iS, so the later ray segments still treat it as if
          ! all steps leading up to them were of size deltas.
          SINT  = ( iS - 1 ) * Beam%deltas + W * Beam%deltas

          CALL IncPhaseIfCaustic( .FALSE. )

          ! WRITE( PRTFile, * ) 'is ir', is-2, ir-1

          RcvrDepths: DO iz = 1, NRz_per_range
             deltaz =  Pos%Rz( iz ) - x( 2 )   ! ray to rcvr distance
             ! LP: Reinstated this condition for eigenrays and arrivals, as
             ! without it every ray would be an eigenray / arrival.
             Adeltaz    = ABS( deltaz )
             IF ( Adeltaz < RadiusMax .OR. Beam%RunType( 1 : 1 ) == 'C' &
                   .OR. Beam%RunType( 1 : 1 ) == 'I' .OR. Beam%RunType( 1 : 1 ) == 'S' ) THEN
                ! LP: Changed to use ApplyContribution in order to support
                ! incoherent, semi-coherent, and arrivals.
                SrcDeclAngle = RadDeg * alpha   ! take-off angle in degrees
                CPA    = ABS( deltaz * ( rB - rA ) ) / SQRT( ( rB - rA )**2 + &
                            & ( ray2D( iS )%x( 2 ) - ray2D( iS - 1 )%x( 2 ) )**2  )
                DS     = SQRT( deltaz **2 - CPA **2 )
                SX1    = SINT + DS
                thet   = ATAN( CPA / SX1 )
                delay  = tau + rayt( 2 ) * deltaz
                const  = Ratio1 * CN * ray2D( iS )%Amp / SQRT( SX1 )
                W      = EXP( -A * thet ** 2 )
                Amp    = const * W
                phaseInt = ray2D( iS )%Phase + phase
                CALL ApplyContribution( U( iz, ir ) )

             END IF
          END DO RcvrDepths

          qOld = q
          ir   = ir + 1
          IF ( ir > Pos%NRr ) RETURN
       END DO RcvrRanges

       rA = rB
    END DO Stepping

  END SUBROUTINE InfluenceSGB

  ! **********************************************************************!

  SUBROUTINE BranchCut( q1C, q2C, BeamType, KMAH )
  !! Checks for a branch cut crossing and updates KMAH accordingly

    COMPLEX  (KIND=8), INTENT( IN )    :: q1C, q2C
    CHARACTER (LEN=4), INTENT( IN )    :: BeamType
    INTEGER,           INTENT( INOUT ) :: KMAH
    REAL     (KIND=8)                  :: q1, q2

    SELECT CASE ( BeamType( 2 : 2 ) )
    CASE ( 'W' )   ! WKBeams
       q1 = REAL( q1C )
       q2 = REAL( q2C )
       IF ( ( q1 < 0.0 .AND. q2 >= 0.0 ) .OR. &
            ( q1 > 0.0 .AND. q2 <= 0.0 ) ) KMAH = -KMAH
    CASE DEFAULT
       IF ( REAL( q2C ) < 0.0 ) THEN
          q1 = AIMAG( q1C )
          q2 = AIMAG( q2C )
          IF ( ( q1 < 0.0 .AND. q2 >= 0.0 ) .OR. &
               ( q1 > 0.0 .AND. q2 <= 0.0 ) ) KMAH = -KMAH
       END IF
    END SELECT

  END SUBROUTINE BranchCut

  ! **********************************************************************!

  SUBROUTINE ScalePressure( Dalpha, c, r, U, NRz, Nr, RunType, freq )
  !! Scale the pressure field

    REAL,              PARAMETER       :: pi = 3.14159265
    INTEGER,           INTENT( IN    ) :: NRz, Nr
    REAL,              INTENT( IN    ) :: r( Nr )         ! ranges
    REAL     (KIND=8), INTENT( IN    ) :: Dalpha, freq, c ! angular spacing between rays, source frequency, nominal sound speed
    COMPLEX,           INTENT( INOUT ) :: U( NRz, Nr )    ! Pressure field
    CHARACTER (LEN=5), INTENT( IN    ) :: RunType
    REAL     (KIND=8)                  :: const, factor

    ! Compute scale factor for field
    SELECT CASE ( RunType( 2 : 2 ) )
    CASE ( 'C' )   ! Cerveny Gaussian beams in Cartesian coordinates
       const = -Dalpha * SQRT( freq ) / c
    CASE ( 'R' )   ! Cerveny Gaussian beams in Ray-centered coordinates
       const = -Dalpha * SQRT( freq ) / c
    CASE DEFAULT
       const = -1.0
    END SELECT

    IF ( RunType( 1 : 1 ) /= 'C' ) U = SQRT( REAL( U ) ) ! For incoherent run, convert intensity to pressure

    ! scale and/or incorporate cylindrical spreading
    Ranges: DO ir = 1, Nr
       IF ( RunType( 4 : 4 ) == 'X' ) THEN   ! line source
          factor = -4.0 * SQRT( pi ) * const
       ELSE                                  ! point source
          IF ( r ( ir ) == 0 ) THEN
             factor = 0.0D0                  ! avoid /0 at origin, return pressure = 0
          ELSE
             factor = const / SQRT( ABS( r( ir ) ) )
          END IF
       END IF
       U( :, ir ) = SNGL( factor ) * U( :, ir )
    END DO Ranges

  END SUBROUTINE ScalePressure

  ! **********************************************************************!

  REAL (KIND=8 ) FUNCTION Hermite( x, x1, x2 )

    ! Calculates a smoothing function based on the h0 hermite cubic
    ! x is the point where the function is to be evaluated
    ! returns:
    ! [  0, x1  ] = 1
    ! [ x1, x2  ] = cubic taper from 1 to 0
    ! [ x2, inf ] = 0

    REAL (KIND=8 ), INTENT( IN  ) :: x, x1, x2
    REAL (KIND=8 )                :: Ax, u

    Ax  = ABS( x  )

    IF ( Ax <= x1 ) THEN
       Hermite = 1.0d0
    ELSE IF ( Ax >= x2 ) THEN
       Hermite = 0.0d0
    ELSE
       u       = ( Ax - x1 ) / ( x2 - x1 )
       Hermite = ( 1.0d0 + 2.0d0 * u ) * ( 1.0d0 - u ) ** 2
    END IF

    !hermit = hermit / ( 0.5 * ( x1 + x2 ) )

  END FUNCTION Hermite

  ! **********************************************************************!

  SUBROUTINE FinalPhase( isGaussian )
    LOGICAL, INTENT( IN ) :: isGaussian
    INTEGER :: phaseStepNum

    !! phase shifts at caustics

    ! this should be precomputed [LP: While IncPhaseIfCaustic can be
    ! precomputed, FinalPhase cannot, as it is dependent on the interpolated `q`
    ! value which is not known until the main run.]
    ! LP: All 2D functions discard the ray point phase if the condition is met,
    ! probably BUG
    ! LP: 2D Gaussian Cartesian reads the phase from the current point, all
    ! others (including 3D) read the phase from the previous point, probably BUG
    IF ( isGaussian ) THEN
       phaseStepNum = iS
    ELSE
       phaseStepNum = iS - 1
    END IF

    phaseInt = ray2D( phaseStepNum )%Phase + phase
    IF ( IsAtCaustic( .TRUE. ) ) THEN
       phaseInt = phase + pi / 2.
    END IF

  END SUBROUTINE FinalPhase

  ! **********************************************************************!

  SUBROUTINE IncPhaseIfCaustic( qleq0 )

    !! phase shifts at caustics

    LOGICAL, INTENT( IN ) :: qleq0

    IF ( IsAtCaustic( qleq0 ) ) THEN
       phase = phase + pi / 2.
    END IF

  END SUBROUTINE IncPhaseIfCaustic

  ! **********************************************************************!

  LOGICAL FUNCTION IsAtCaustic( qleq0 )

    ! LP: There are two versions of the phase shift condition used in the
    ! BELLHOP code, with the equalities in opposite positions. qleq0 false is
    ! only used in SGB.

    LOGICAL, INTENT( IN ) :: qleq0

    IF ( qleq0 ) THEN
       IsAtCaustic = q <= 0.0d0 .AND. qOld >  0.0d0 .OR. q >= 0.0d0 .AND. qOld <  0.0d0
    ELSE
       IsAtCaustic = q <  0.0d0 .AND. qOld >= 0.0d0 .OR. q >  0.0d0 .AND. qOld <= 0.0d0
    END IF

  END FUNCTION IsAtCaustic

END MODULE Influence
