Influence3DGeoHatCart Subroutine

  SUBROUTINE Influence3DGeoHatCart( alpha, beta, Dalpha, Dbeta, P, x_rcvrMat, t_rcvr )
    !! Geometrically-spreading beams with a hat-shaped beam

    REAL ( KIND=8 ), INTENT( IN  ) :: alpha, beta, Dalpha, Dbeta         ! ray take-off angle
    REAL ( KIND=8 ), INTENT( IN  ) :: x_rcvrMat( 2, Pos%Ntheta, Pos%NRr ), t_rcvr( 2, Pos%Ntheta ) ! rcvr coordinates and tangent
    COMPLEX        , INTENT( OUT ) :: P( Pos%Ntheta, Pos%Nrz, Pos%NRr )  ! complex pressure
    INTEGER            :: irT( 1 ), irTT
    REAL    ( KIND=8 ) :: s, rlen, x_ray( 3 ), n_ray_z( 3 ), n_ray_theta(3 ), &
         e1( 3 ), e2( 3 ), x_rcvr( 3 ), x_rcvr_ray( 3 ), &
         L_z, L_diag, e_theta( 3 ), m_prime, a, b, zMin, zMax, &
         Det_Q, Det_Qold
    Ratio1 = SQRT( ABS( COS( alpha ) ) ) * SQRT( Dalpha * Dbeta ) / ray3D( 1 )%c
    CALL ScaleBeam( alpha, Dalpha, Dbeta )   ! scaling for geometric beams

    phase    = 0.0
    Det_QOld = ray3D( 1 )%DetQ   ! used to track phase changes at caustics (rotations of Det_Q)

    ! Compute nearest rcvr before normal
    ! LP: Silly / partly dead code, see README.
    rA  = NORM2( ray3D( 1 )%x( 1 : 2 ) - xs_3D( 1 : 2 ) )         ! range of ray point
    irT = MINLOC( Pos%Rr( 1 : Pos%NRr ), MASK = Pos%Rr( 1 : Pos%NRr ) > rA )        ! index of receiver
    ir  = irT( 1 )

    Stepping: DO is = 2, Beam%Nsteps
       ! Compute nearest rcvr before normal
       rB  = NORM2( ray3D( is )%x( 1 : 2 ) - xs_3D( 1 : 2 ) )         ! range of ray point

       IF ( ABS( rB - rA ) > 1.0D3 * SPACING( rA ) ) THEN   ! jump to next step if duplicate point
          ! initialize the index of the receiver range
          IF ( is == 2 ) THEN
             ! LP: Silly / partly dead code, see README.
             IF ( rB > rA ) THEN   ! ray is moving outward in range
                ir = 1             ! index all the way in
             ELSE                  ! ray is moving inward in range
                ir = Pos%NRr       ! index all the way out
             END IF
          END IF

          x_ray = ray3D( is - 1 )%x

          ! compute normalized tangent (we need it to measure the step length)
          rayt = ray3D( is )%x - ray3D( is - 1 )%x
          rlen = NORM2( rayt )

          IF ( rlen > 1.0D3 * SPACING( ray3D( is )%x( 1 ) ) ) THEN  ! Make sure this is not a duplicate point
              ! WRITE( PRTFile, * ) 'is rA rB', is-2, rA, rB

             rayt = rayt / rlen                                     ! unit tangent to ray
             CALL RayNormal_unit( rayt, ray3D( is )%phi, e1, e2 )   ! Get ray normals e1 and e2

             ! phase shifts at caustics
             Det_Q  = ray3D( is - 1 )%DetQ
             IF ( Det_Q <= 0.0d0 .AND. Det_QOld > 0.0d0 .OR. Det_Q >= 0.0d0 .AND. Det_QOld < 0.0d0 ) phase = phase + pi / 2.
             Det_Qold = Det_Q

             L1 = MAX( NORM2( ray3D( is - 1 )%q_tilde ), NORM2( ray3D( is )%q_tilde ) ) ! beamwidths
             L2 = MAX( NORM2( ray3D( is - 1 )%q_hat   ), NORM2( ray3D( is )%q_hat   ) )

             L_diag = SQRT( L1 ** 2 + L2 ** 2 )   ! worst case is when rectangle is rotated to catch the hypotenuse

             n_ray_theta = [ -rayt( 2 ), rayt( 1 ), 0.D0 ]  ! normal to the ray in the horizontal receiver plane

             ! *** Compute contributions to bracketed receivers ***
             dq_tilde = ray3D( is )%q_tilde - ray3D( is - 1 )%q_tilde
             dq_hat   = ray3D( is )%q_hat   - ray3D( is - 1 )%q_hat
             dtau     = ray3D( is )%tau     - ray3D( is - 1 )%tau

             Ranges: DO
                ! is r( ir ) contained in [ rA, rB ]? Then compute beam influence
                IF ( Pos%Rr( ir ) >= MIN( rA, rB ) .AND. Pos%Rr( ir ) < MAX( rA, rB ) ) THEN

                   Radials: DO itheta = 1, Pos%Ntheta   ! Loop over radials of receiver line

                      x_rcvr( 1 : 2 ) = x_rcvrMat( 1 : 2, itheta, ir )
                      m_prime = ABS( DOT_PRODUCT( x_rcvr( 1 : 2 ) - x_ray( 1 : 2 ), n_ray_theta( 1 : 2 ) ) )  ! normal distance from rcvr to ray segment

                      IF ( m_prime > L_diag ) THEN
                         ! IF ( itheta == 36 ) THEN
                         !    WRITE( PRTFile, * ) 'Skip theta b/c m_prime', m_prime, L_diag
                         ! END IF
                         CYCLE Radials
                      END IF

                      ! The set of possible receivers is a ring
                      ! However, extrapolating the beam backwards produces contributions with s negative and large
                      ! We do not want to accept these contributions--- they have the proper range but are 180 degrees
                      ! away from this segment of the ray
                      ! s = DOT_PRODUCT( x_rcvr( 1 : 2 ) - x_ray( 1 : 2 ), rayt( 1 : 2 ) / NORM2( rayt( 1 : 2 ) ) )   ! a distance along ray    (in x-y plane)
                      ! s = DOT_PRODUCT( x_rcvr( 1 : 2 ) - xs_3D( 1 : 2 ), t_rcvr( 1 : 2, itheta ) ) ! a distance along radial (in x-y plane)
                      s = DOT_PRODUCT( x_rcvr( 1 : 2 ) - xs_3D( 1 : 2 ), x_ray( 1 : 2 ) - xs_3D( 1 : 2 ) ) ! vector to rcvr dotted into vector to ray point

                      ! The real receivers have an s-value in [0, R_len]
                      ! IF ( s < 0D0 .OR. s > NORM2( ray3D( is )%x( 1 : 2 ) - ray3D( is - 1 )%x( 1 : 2 ) ) ) THEN
                      ! IF ( ABS( s ) > NORM2( ray3D( is )%x( 1 : 2 ) - ray3D( is - 1 )%x( 1 : 2 ) ) ) THEN

                      IF ( s < 0D0 ) THEN
                         ! WRITE( PRTFile, * ) 'Skip theta b/c s'
                         CYCLE Radials
                      END IF

                      ! calculate z-limits for the beam (could be pre-cacluated for each itheta)
                      ! LP: mbp seems to have realized only some of these components were being used
                      ! and is only calculating the needed ones, but still storing all of them
                      e_theta      = [ -t_rcvr( 2, itheta ), t_rcvr( 1, itheta ), 0.0D0 ]  ! normal to the vertical receiver plane
                      ! n_ray_z    = CROSS_PRODUCT( rayt, e_theta )                        ! normal to the ray in the vertical receiver plane
                      n_ray_z( 3 ) = rayt( 1 ) * e_theta( 2 ) - rayt( 2 ) * e_theta( 1 )   ! normal to the ray in the vertical receiver plane

                      IF ( ABS( n_ray_z( 3 ) ) < 1D-9 ) THEN
                         ! WRITE( PRTFile, * ) 'Skip theta b/c L_z divide by zero', n_ray_z( 3 )
                         CYCLE Radials   ! avoid divide by zero
                      END IF
                      L_z          = L_diag / ABS( n_ray_z( 3 ) )

                      zmin = MIN( ray3D( is - 1 )%x( 3 ), ray3D( is )%x( 3 ) ) - L_z  ! min depth of ray segment
                      zmax = MAX( ray3D( is - 1 )%x( 3 ), ray3D( is )%x( 3 ) ) + L_z  ! max depth of ray segment

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

                      ReceiverDepths: DO iz = 1, Pos%Nrz
                         x_rcvr( 3 ) = DBLE( Pos%rz( iz ) )   ! z coordinate of the receiver
                         IF ( x_rcvr( 3 ) < zmin .OR. x_rcvr( 3 ) > zmax ) CYCLE ReceiverDepths

                         x_rcvr_ray = x_rcvr - x_ray

!!! rlen factor could be built into rayt
                         ! linear interpolation of q's
                         s       = DOT_PRODUCT( x_rcvr_ray, rayt ) / rlen  ! proportional distance along ray
                         q_tilde = ray3D( is - 1 )%q_tilde + s * dq_tilde
                         q_hat   = ray3D( is - 1 )%q_hat   + s * dq_hat

                         L1 = NORM2( q_tilde )  ! beamwidth is length of the vector
                         L2 = NORM2( q_hat   )  ! beamwidth

                         n  = ABS( DOT_PRODUCT( x_rcvr_ray, e1 ) )         ! normal distance to ray
                         m  = ABS( DOT_PRODUCT( x_rcvr_ray, e2 ) )         ! normal distance to ray

                         ! represent (m, n) as a linear combination a q + b q
                         DetQint = q_tilde( 1 ) * q_hat( 2 ) - q_hat( 1 ) * q_tilde( 2 )   ! area of parallelogram formed by ray tube

                         IF ( DetQint == 0.0 ) CYCLE ReceiverDepths
                         a = ABS( ( -q_hat(   1 ) * m + q_hat(   2 ) * n ) / DetQint )
                         b = ABS( (  q_tilde( 1 ) * m - q_tilde( 2 ) * n ) / DetQint )

                         ! suppose qtilde( 2 ) = qhat( 1 ) = 0
                         !DetQint = q_tilde( 1 ) * q_hat( 2 )   ! area of parallelogram formed by ray tube
                         !a = ABS( n / q_tilde( 1 )  )
                         !b = ABS( m / q_hat(   2 )  )

                         IF ( MAX( a, b ) > 1 .OR. L1 == 0 .OR. L2 == 0 ) CYCLE ReceiverDepths   ! receiver is outside the beam

                         delay = ray3D( is - 1 )%tau + s * dtau

                         ! phase shift at caustics
                         phaseInt = ray3D( is - 1 )%Phase + phase
                         IF ( DetQint <= 0.0d0 .AND. Det_QOld > 0.0d0 .OR. &
                              DetQint >= 0.0d0 .AND. Det_QOld < 0.0d0 ) phaseInt = phaseInt + pi / 2.

                         const = ray3D( is )%Amp / SQRT( ABS( DetQint ) )
                         W     = ( 1 - a ) * ( 1 - b )   ! hat function: 1 on center, 0 on edge
                         Amp   = const * W

                         ! WRITE( PRTFile, * ) iz-1

                         CALL ApplyContribution( alpha, beta, P( itheta, iz, ir ) )
                      END DO ReceiverDepths
                   END DO Radials

                END IF
                ! bump receiver index, ir, towards rB
                IF ( Pos%Rr( ir ) < rB ) THEN
                   IF ( ir >= Pos%NRr ) EXIT Ranges   ! jump out of the search and go to next step on ray
                   irTT = ir + 1                      ! bump right
                   IF ( Pos%Rr( irTT ) >= rB ) EXIT Ranges
                ELSE
                   IF ( ir <= 1  ) EXIT Ranges        ! jump out of the search and go to next step on ray
                   irTT = ir - 1                      ! bump left
                   IF ( Pos%Rr( irTT ) <= rB ) EXIT Ranges
                END IF
                ir = irTT
             END DO Ranges
          END IF
       END IF
       rA = rB
    END DO Stepping

  END SUBROUTINE Influence3DGeoHatCart