InfluenceSGB Subroutine

  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