InfluenceGeoHatCart Subroutine

  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