InfluenceGeoGaussianCart Subroutine

  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