!! Three-dimensional beam influence computation for pressure fields

MODULE Influence3D
  !! 3D beam influence calculations with complex pressure field contributions and spatial weighting

  ! 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 WriteRay
  USE RayNormals
  USE SourceReceiverPositions
  USE ArrMod
  USE sspMod   ! used to construct image beams in the Cerveny style beam routines
  USE cross_products

  IMPLICIT NONE
  PUBLIC

  INTEGER,       PRIVATE :: itheta, iz, ir, is
  REAL (KIND=8), PRIVATE :: W, s, m, n, Amp, phase, const, phaseInt, Ratio1, &
       L1, L2, rayt( 3 ), &
       RcvrDeclAngle, RcvrAzimAngle, &
       rA, rB, lambda
  REAL     (KIND=8) :: q_tilde( 2 ), q_hat( 2 ), dq_tilde( 2 ), dq_hat( 2 ), DetQint
  COMPLEX ( KIND=8 ) :: delay, dtau
  LOGICAL, PARAMETER, PRIVATE :: INFL_DEBUGGING = .FALSE.
  INTEGER, PARAMETER, PRIVATE :: INFL_DEBUGGING_ITHETA = 0, INFL_DEBUGGING_IR = 299, &
       INFL_DEBUGGING_IZ = 0

CONTAINS

  SUBROUTINE Influence3DGeoHatRayCen( alpha, beta, Dalpha, Dbeta, P )
    !! Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates

    REAL ( KIND=8 ), INTENT( IN  ) :: alpha, beta, Dalpha, Dbeta         ! ray take-off angle
    COMPLEX        , INTENT( OUT ) :: P( Pos%Ntheta, Pos%Nrz, Pos%NRr )  ! complex pressure
    INTEGER          :: irA, irB, II
    REAL    (KIND=8) :: nA, nB, mA, mB, deltaA, deltaB, t_rcvr( 2, Pos%Ntheta ), &
         a, b, KMAHphase( Beam%Nsteps )
    REAL    (KIND=8) :: e1G( 3, Beam%Nsteps ), e2G( 3, Beam%Nsteps ), e1xe2( 3, Beam%Nsteps ), &
         xt( 3, Beam%Nsteps ),  xtxe1( 3, Beam%Nsteps ),  xtxe2( 3, Beam%Nsteps )

    Ratio1 = SQRT( ABS( COS( alpha ) ) ) * SQRT( Dalpha * Dbeta ) / ray3D( 1 )%c
    CALL ScaleBeam( alpha, Dalpha, Dbeta )   ! scaling for geometric beams

    ! phase shifts at caustics (rotations of Det_Q)
    KMAHphase( 1 ) = 0
    DO is = 2, Beam%Nsteps
       KMAHphase( is ) = KMAHphase( is - 1 )
       IF ( ray3D( is )%DetQ <= 0.0d0 .AND. ray3D( is - 1 )%DetQ > 0.0d0 .OR. &
            ray3D( is )%DetQ >= 0.0d0 .AND. ray3D( is - 1 )%DetQ < 0.0d0 ) KMAHphase( is ) = KMAHphase( is - 1 ) + pi / 2.
    END DO

    ! pre-calculate tangents, normals (e1, e2), etc.
    DO is = 1, Beam%Nsteps
       xt( :, is ) = ray3D( is )%x - xs_3D   ! vector from the origin of the receiver plane to this point on the ray
       CALL RayNormal( ray3D( is )%t, ray3D( is )%phi, ray3D( is )%c, e1G( :, is ), e2G( :, is ) )

       ! e1xe2( :, is ) = cross_product( e1G( :, is ), e2G( :, is ) )
       e1xe2( :, is ) = ray3D( is )%c * ray3D( is )%t   ! unit tangent to ray
    END DO

    ! tangent along receiver bearing line
    t_rcvr( 1, : ) = COS( DegRad * Pos%theta( 1 : Pos%Ntheta ) )
    t_rcvr( 2, : ) = SIN( DegRad * Pos%theta( 1 : Pos%Ntheta ) )

    ReceiverDepths: DO iz = 1, Pos%Nrz
       ! precalculate vector from receiver to each step of ray
       xt( 3, : ) = ray3D( 1 : Beam%Nsteps )%x( 3 ) - Pos%rz( iz )
       DO is = 1, Beam%Nsteps
          xtxe1( :, is ) = cross_product( xt( :, is ), e1G( :, is ) )
          xtxe2( :, is ) = cross_product( xt( :, is ), e2G( :, is ) )
       END DO

       Radials: DO itheta = 1, Pos%Ntheta
          ! step along the beam ...
          ! Most of the time the beam makes no contribution to a receiver
          ! Therefore we try to test that quickly and move on to the next receiver

          Stepping: DO is = 2, Beam%Nsteps

             ! *** Compute coordinates of intercept: nB, mB, rB ***
             deltaA = -DOT_PRODUCT( t_rcvr( :, itheta ), e1xe2( 1 : 2, is - 1 ) )
             deltaB = -DOT_PRODUCT( t_rcvr( :, itheta ), e1xe2( 1 : 2, is     ) )

             IF ( ABS( deltaA ) < 1e-3 .OR. ABS( deltaB ) < 1e-3 ) CYCLE Stepping   ! ray normal || radial of rcvr line

             mA =  DOT_PRODUCT( t_rcvr( :, itheta ), xtxe1( 1 : 2, is - 1 ) ) / deltaA
             mB =  DOT_PRODUCT( t_rcvr( :, itheta ), xtxe1( 1 : 2, is     ) ) / deltaB

             nA = -DOT_PRODUCT( t_rcvr( :, itheta ), xtxe2( 1 : 2, is - 1 ) ) / deltaA
             nB = -DOT_PRODUCT( t_rcvr( :, itheta ), xtxe2( 1 : 2, is     ) ) / deltaB

             rA  = -DOT_PRODUCT( xt( :, is - 1 ), e1xe2( :, is - 1 ) ) / deltaA
             rB  = -DOT_PRODUCT( xt( :, is     ), e1xe2( :, is     ) ) / deltaB

             irA = RToIR( rA )
             irB = RToIR( rB )
             ! detect and skip duplicate points (happens at boundary reflection)
             ! IF ( irA /= irB .AND. NORM2( ray3D( is )%x - ray3D( is - 1 )%x ) > 1.0D3 * SPACING( ray3D( is )%x( 1 ) ) ) THEN
             ! ^^ too slow ^^
             IF ( irA /= irB .AND. NORM2( ray3D( is )%x - ray3D( is - 1 )%x ) > 1.0e-4 ) THEN

                ! *** 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

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

                Ranges: DO ir = irA + 1 - II, irB + II, SIGN( 1, irB - irA )
                   s = ( Pos%Rr( ir ) - rA ) / ( rB - rA )

                   ! linear interpolation of q's
                   q_tilde = ray3D( is - 1 )%q_tilde + s * dq_tilde
                   q_hat   = ray3D( is - 1 )%q_hat   + s * dq_hat
                   n       = ABS( nA                 + s * ( nB - nA ) )   ! normal distance to ray
                   m       = ABS( mA                 + s * ( mB - mA ) )   ! 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    ) CYCLE Ranges  ! receiver is outside the beam

                   a = ABS( ( -q_hat(   1 ) * m + q_hat(   2 ) * n ) / DetQint ) !!! sign is flippable inside abs
                   b = ABS( (  q_tilde( 1 ) * m - q_tilde( 2 ) * n ) / DetQint )
                   IF ( MAX( a, b ) > 1 ) CYCLE Ranges  ! receiver is outside the beam

                   W = ( 1 - a ) * ( 1 - b )   ! beamshape is a hat function: 1 on center, 0 on edge
                   delay = ray3D( is - 1 )%tau + s * dtau

                   const = ray3D( is )%Amp / SQRT( ABS( DetQint ) )
                   Amp   = const * W   ! hat function

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

                   CALL ApplyContribution( alpha, beta, P( itheta, iz, ir ) )
                END DO Ranges
             END IF
             ! LP: These values are overwritten on the next loop; this does nothing.
             mA     = mB
             deltaA = deltaB
          END DO Stepping
       END DO Radials
    END DO ReceiverDepths

  END SUBROUTINE Influence3DGeoHatRayCen

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

  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 )
                      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

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

  SUBROUTINE Influence3DGeoGaussianRayCen( alpha, beta, Dalpha, Dbeta, P )
    !! Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates

    REAL, PARAMETER                :: BeamWindow = 4.0                   ! kills beams outside e**(-0.5 * BeamWindow**2 )
    REAL ( KIND=8 ), INTENT( IN  ) :: alpha, beta, Dalpha, Dbeta         ! ray take-off angle
    COMPLEX        , INTENT( OUT ) :: P( Pos%Ntheta, Pos%Nrz, Pos%NRr )  ! complex pressure
    INTEGER          :: irA, irB, II
    REAL    (KIND=8) :: nA, nB, mA, mB, a, b, deltaA, deltaB, t_rcvr( 2, Pos%Ntheta ), &
         L1_stent, L2_stent, KMAHphase( Beam%Nsteps )
    REAL    (KIND=8) :: e1G( 3, Beam%Nsteps ), e2G( 3, Beam%Nsteps ), e1xe2( 3, Beam%Nsteps ), &
         xt( 3, Beam%Nsteps ),  xtxe1( 3, Beam%Nsteps ),  xtxe2( 3, Beam%Nsteps )
    REAL    (KIND=8) :: MaxRadius_m( Beam%Nsteps - 1 ), MaxRadius_n( Beam%Nsteps - 1 )

    Ratio1 = SQRT( ABS( COS( alpha ) ) ) * SQRT( Dalpha * Dbeta ) / ray3D( 1 )%c / ( 2.0 * pi )
    CALL ScaleBeam( alpha, Dalpha, Dbeta )   ! scaling for geometric beams

    ! phase shifts at caustics (rotations of Det_Q)
    KMAHphase( 1 ) = 0
    DO is = 2, Beam%Nsteps
       KMAHphase( is ) = KMAHphase( is - 1 )
       IF ( ray3D( is )%DetQ <= 0.0d0 .AND. ray3D( is - 1 )%DetQ > 0.0d0 .OR. &
            ray3D( is )%DetQ >= 0.0d0 .AND. ray3D( is - 1 )%DetQ < 0.0d0 ) KMAHphase( is ) = KMAHphase( is - 1 ) + pi / 2.
    END DO

    ! pre-calculate tangents, normals (e1, e2), etc.
    DO is = 1, Beam%Nsteps
       xt( :, is ) = ray3D( is )%x - xs_3D   ! vector from the origin of the receiver plane to this point on the ray
       CALL RayNormal( ray3D( is )%t, ray3D( is )%phi, ray3D( is )%c, e1G( :, is ), e2G( :, is ) )

       ! e1xe2( :, is ) = cross_product( e1G( :, is ), e2G( :, is ) )
       e1xe2( :, is ) = ray3D( is )%c * ray3D( is )%t   ! unit tangent to ray
    END DO

    DO is = 1, Beam%Nsteps - 1
       MaxRadius_m( is ) = BeamWindow * MAX( ABS( ray3D( is )%q_hat(   2 ) ), ABS( ray3D( is + 1 )%q_hat(   2 ) ) )
       MaxRadius_n( is ) = BeamWindow * MAX( ABS( ray3D( is )%q_tilde( 1 ) ), ABS( ray3D( is + 1 )%q_tilde( 1 ) ) )
    END DO

    ! tangent along receiver bearing line
    t_rcvr( 1, : ) = COS( DegRad * Pos%theta( 1 : Pos%Ntheta ) )
    t_rcvr( 2, : ) = SIN( DegRad * Pos%theta( 1 : Pos%Ntheta ) )

    ReceiverDepths: DO iz = 1, Pos%Nrz
       ! precalculate vector from receiver to each step of ray
       xt( 3, : ) = ray3D( 1 : Beam%Nsteps )%x( 3 ) - Pos%rz( iz )
       DO is = 1, Beam%Nsteps
          xtxe1( :, is ) = cross_product( xt( :, is ), e1G( :, is ) )
          xtxe2( :, is ) = cross_product( xt( :, is ), e2G( :, is ) )
       END DO

       Radials: DO itheta = 1, Pos%Ntheta
          ! *** Compute coordinates of intercept: nA, mA, rA ***
          is = 1

          ! LP: These values are immediately overwritten inside the loop. This
          ! looks like an initialization of these values, and they're saved from
          ! B to A at the end of the loop, but the computation of deltaA, mA,
          ! and irA would have to be removed below.
          deltaA = -DOT_PRODUCT( t_rcvr( :, itheta ), e1xe2( 1 : 2, is ) )

          ! Check for ray normal || radial of rcvr line
          IF ( ABS( deltaA ) < 1D3 * SPACING( deltaA ) ) THEN
             irA = 0   ! serves as a flag that this normal can't be used
          ELSE
             mA  =  DOT_PRODUCT( t_rcvr( :, itheta ), xtxe1( 1 : 2, is ) ) / deltaA
          END IF

          ! step along the beam ...
          ! Most of the time the beam makes no contribution to a receiver
          ! Therefore we try to test that quickly and move on to the next receiver

          Stepping: DO is = 2, Beam%Nsteps
             lambda = ray3D( is - 1 )%c / freq   ! local wavelength

             ! *** Compute coordinates of intercept: nB, mB, rB ***
             deltaA = -DOT_PRODUCT( t_rcvr( :, itheta ), e1xe2( 1 : 2, is - 1 ) )
             deltaB = -DOT_PRODUCT( t_rcvr( :, itheta ), e1xe2( 1 : 2, is     ) )

             IF ( ABS( deltaA ) < 1e-3 .OR. ABS( deltaB ) < 1e-3 ) THEN
                IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                   .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                   WRITE( PRTFile, * ) 'Skipping b/c deltaA deltaB', deltaA, deltaB
                END IF
                CYCLE Stepping   ! ray normal || radial of rcvr line
             END IF

             mA  =  DOT_PRODUCT( t_rcvr( :, itheta ), xtxe1( 1 : 2, is - 1 ) ) / deltaA
             mB  =  DOT_PRODUCT( t_rcvr( :, itheta ), xtxe1( 1 : 2, is     ) ) / deltaB
!!! stent needs to be applied here already

             ! Possible contribution if max possible beamwidth > min possible distance to receiver
             IF ( MaxRadius_m( is - 1 ) > MIN( ABS( mA ), ABS( mB ) ) .OR. ( mA * mB < 0 ) ) THEN
                nA = -DOT_PRODUCT( t_rcvr( :, itheta ), xtxe2( 1 : 2, is - 1 ) ) / deltaA
                nB = -DOT_PRODUCT( t_rcvr( :, itheta ), xtxe2( 1 : 2, is     ) ) / deltaB

                ! Possible contribution if max possible beamwidth > min possible distance to receiver
                IF ( MaxRadius_n( is - 1 ) > MIN( ABS( nA ), ABS( nB ) ) .OR. ( nA * nB < 0 ) ) THEN
!!! can generate an inexact result if the resulting integer is too big
!!! assumes uniform spacing in Pos%Rr
                   rA  = -DOT_PRODUCT( xt( :, is - 1 ), e1xe2( :, is - 1 ) ) / deltaA
                   rB  = -DOT_PRODUCT( xt( :, is     ), e1xe2( :, is     ) ) / deltaB

                   irA = RToIR( rA )
                   irB = RToIR( rB )

                   IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                      .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                      WRITE( PRTFile, * ) 'rA rB irA irB', rA, rB, irA, irB
                   END IF

                   ! detect and skip duplicate points (happens at boundary reflection)
                   IF ( irA /= irB .AND. NORM2( ray3D( is )%x - ray3D( is - 1 )%x ) &
                      & > 1.0D3 * SPACING( ray3D( is )%x( 1 ) ) ) THEN

                      ! *** 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

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

                      Ranges: DO ir = irA + 1 - II, irB + II, SIGN( 1, irB - irA )

                         s = ( Pos%Rr( ir ) - rA ) / ( rB - rA )

                         ! linear interpolation of q's
                         q_tilde = ray3D( is - 1 )%q_tilde + s * dq_tilde
                         q_hat   = ray3D( is - 1 )%q_hat   + s * dq_hat
                         n       = ABS( nA                 + s * ( nB - nA ) )   ! normal distance to ray
                         m       = ABS( mA                 + s * ( mB - mA ) )   ! 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    ) CYCLE Ranges  ! receiver is outside the beam

                         a = ABS( ( -q_hat(   1 ) * m + q_hat(   2 ) * n ) / DetQint ) !!! sign is flipable inside abs
                         b = ABS( (  q_tilde( 1 ) * m - q_tilde( 2 ) * n ) / DetQint )

                         !IF ( a + b > BeamWindow ) CYCLE Ranges   ! receiver is outside the beam .OR. L1 == 0 .OR. L2 == 0

                         ! Beam shape function (Gaussian)
                         ! The factor involving L1, L2 compensates for the change in intensity when widening the beam
                         !W   = EXP( -.5 * ( ( n / L1_stent ) ** 2 + ( m / L2_stent ) ** 2 ) ) * L1 / L1_stent * L2 / L2_stent
                         W   = EXP( -.5 * ( a ** 2 + b ** 2 ) )
                         delay = ray3D( is - 1 )%tau + s * dtau

                         const = ray3D( is )%Amp / SQRT( ABS( DetQint ) )
                         Amp   = const * W   ! hat function
                         ! phase shift at caustics
                         phaseInt = ray3D( is - 1 )%Phase + KMAHphase( is - 1 )
                         IF ( DetQint <= 0.0d0 .AND. ray3D( is - 1 )%DetQ > 0.0d0 .OR. &
                              DetQint >= 0.0d0 .AND. ray3D( is - 1 )%DetQ < 0.0d0 ) phaseInt = phaseInt + pi / 2.

                         IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                            .AND. ir-1 == INFL_DEBUGGING_IR .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                            WRITE( PRTFile, * ) 'is itheta iz ir: const W delay phaseInt', &
                               is-2, itheta-1, iz-1, ir-1, const, W, delay, phaseInt
                         END IF

                         CALL ApplyContribution( alpha, beta, P( itheta, iz, ir ) )
!!$                         ______________
!!$
!!$                         ! Within beam window?
!!$                         L1 = ABS( q_tilde( is - 1 ) + W * dq_tilde( is - 1 ) )   ! beamwidth
!!$                         IF ( n > BeamWindow * L1 ) CYCLE Ranges                  ! in beamwindow?
!!$
!!$                         L2 = ABS( q_hat(   is - 1 ) + W * dq_hat(   is - 1 ) )   ! beamwidth
!!$                         IF ( m > BeamWindow * L2 ) CYCLE Ranges                  ! in beamwwindow?
!!$
!!$                         ! calculate the beamwidth (must be at least lambda, except in the nearfield)
!!$                         ! comment out the following 2 lines to turn the stent off
!!$                         L1_stent = MAX( L1, MIN( 0.2 * freq * REAL( ray3D( is - 1 )%tau ), pi * lambda ) )
!!$                         L2_stent = MAX( L2, MIN( 0.2 * freq * REAL( ray3D( is - 1 )%tau ), pi * lambda ) )
!!$
!!$                         DetQint = L1 * L2 * ray3D( 1 )%c ** 2 / ( Dalpha * Dbeta )  ! based on actual beamwidth
                      END DO Ranges
                   ELSE
                      IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                         .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                         IF ( irA == irB ) THEN
                            WRITE( PRTFile, * ) 'Skipping b/c irA == irB'
                         ELSE
                            WRITE( PRTFile, * ) 'Skipping b/c dupl'
                         END IF
                      END IF
                   END IF
                ELSE
                   IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                      .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                      WRITE( PRTFile, * ) 'Skipping b/c MaxRadius_n nA nB', MaxRadius_n( is - 1 ), nA, nB
                   END IF
                END IF
             ELSE
                IF ( INFL_DEBUGGING .AND. itheta-1 == INFL_DEBUGGING_ITHETA &
                   .AND. iz-1 == INFL_DEBUGGING_IZ ) THEN
                   WRITE( PRTFile, * ) 'Skipping b/c MaxRadius_m mA mB', MaxRadius_m( is - 1 ), mA, mB
                END IF
             END IF
             ! LP: These values are overwritten on the next loop; this does nothing.
             mA     = mB
             deltaA = deltaB
          END DO Stepping
       END DO Radials
    END DO ReceiverDepths

  END SUBROUTINE Influence3DGeoGaussianRayCen


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

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

    REAL, PARAMETER                :: BeamWindow = 4.0                   ! kills beams outside e**(-0.5 * BeamWindow**2 )
    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 ) :: rlen, a, b, x_ray( 3 ), n_ray_z( 3 ), n_ray_theta( 3 ), &
         e1( 3 ), e2( 3 ), x_rcvr( 3 ), x_rcvr_ray( 3 ), &
         L1_stent, L2_stent, L_z, L_diag, e_theta( 3 ), m_prime, zMin, zMax, &
         Det_Q, Det_Qold

    Ratio1 = SQRT( ABS( COS( alpha ) ) ) * SQRT( Dalpha * Dbeta ) / ray3D( 1 )%c / ( 2.0 * pi )
    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: Same issues as hat cart.
    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
       lambda = ray3D( is - 1 )%c / freq   ! local wavelength
       ! 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
             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
             rayt = rayt / rlen                                     ! unit tangent to ray
             CALL RayNormal_unit( rayt, ray3D( is )%phi, e1, e2 )   ! Get ray normals e1 and e2
             ! WRITE( PRTFile, * ) 'rayt phi', rayt, ray3D( is )%phi
             ! WRITE( PRTFile, * ) 'e1 e2', e1, 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 > BeamWindow * L_diag ) CYCLE Radials

                      ! 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
!!! pre-calculate unit ray tangent
                      ! 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 ) CYCLE Radials

                      ! calculate z-limits for the beam (could be pre-cacluated for each itheta)
                      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 )
                      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 ) CYCLE Radials   ! avoid divide by zero
                      L_z          = BeamWindow * 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
                         ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                         !    WRITE( PRTFile, * ) 'x_rcvr_ray e1, e2', x_rcvr_ray, e1, e2
                         ! END IF

!!! 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

                         IF ( L1 == 0 .OR. L2 == 0 ) THEN
                            ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                            !    WRITE( PRTFile, * ) 'Skipping z b/c L1/L2', iz
                            ! END IF
                            CYCLE ReceiverDepths
                         END IF

                         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 ) THEN
                            ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                            !    WRITE( PRTFile, * ) 'Skipping z b/c DetQint', iz
                            ! END IF
                            CYCLE ReceiverDepths
                         END IF

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

                         IF ( a + b > BeamWindow ) THEN
                            ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                            !    WRITE( PRTFile, * ) 'Skipping z b/c outside beam (a+b) BeamWindow', &
                            !      iz, a+b, BeamWindow
                            ! END IF
                            CYCLE ReceiverDepths   ! receiver is outside the beam
                         END IF
                         ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                         !   WRITE( PRTFile, * ) 'q_tilde', q_tilde
                         !   WRITE( PRTFile, * ) 'q_hat  ', q_hat
                         !   WRITE( PRTFile, * ) 'DetQint m n', DetQint, n, m
                         !   WRITE( PRTFile, * ) 'a b', a, b
                         ! END IF

                         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 ) )

                         ! calculate the beamwidth (must be at least lambda, except in the nearfield)
                         ! comment out the following 2 lines to turn the stent off
                         !L1_stent  = MAX( L1, MIN( 0.2 * freq * REAL( ray3D( is )%tau ), pi * lambda ) )
                         !L2_stent  = MAX( L2, MIN( 0.2 * freq * REAL( ray3D( is )%tau ), pi * lambda ) )

                         ! Beam shape function (Gaussian)
                         ! The factor involving L1, L2 compensates for the change in intensity when widening the beam

                         !W   = EXP( -.5 * ( n ** 2 + m ** 2 ) ) * L1 / L1_stent * L2 / L2_stent
                         W   = EXP( -.5 * ( a ** 2 + b ** 2 ) )   ! Gaussian
                         Amp = const * W

                         ! IF ( ir-1 == 7 .AND. itheta-1 == 59 ) THEN
                         !    WRITE( PRTFile, * ) 'iz const W delay phaseInt', iz, const, W, delay, phaseInt
                         ! END IF

                         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 Influence3DGeoGaussianCart

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

  SUBROUTINE ScaleBeam( alpha, Dalpha, Dbeta )
    !! Scaling for geometric beams

    REAL ( KIND = 8 ), INTENT( IN ) :: alpha, Dalpha, Dbeta

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

    ray3D( 1 : Beam%Nsteps )%DetQ = ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) * ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) - &
                                    ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * ray3D( 1 : Beam%Nsteps )%q_hat(   1 )

    ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) =                       Dalpha * ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) / ray3D( 1 )%c
    ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) =                       Dalpha * ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) / ray3D( 1 )%c
    ray3D( 1 : Beam%Nsteps )%q_hat(   1 ) = ABS( COS( alpha ) ) * Dbeta  * ray3D( 1 : Beam%Nsteps )%q_hat(   1 ) / ray3D( 1 )%c
    ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) = ABS( COS( alpha ) ) * Dbeta  * ray3D( 1 : Beam%Nsteps )%q_hat(   2 ) / ray3D( 1 )%c

  END SUBROUTINE ScaleBeam

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

  SUBROUTINE ApplyContribution( alpha, beta, U )

    REAL ( KIND=8 ), INTENT( IN    ) :: alpha, beta
      !! ray take-off angle
    COMPLEX,         INTENT( INOUT ) :: U
    REAL ( KIND=8 )                  :: rayt2( 3 )
    ! LP: Don't want to clobber rayt!

    SELECT CASE( Beam%RunType( 1 : 1 ) )
    CASE ( 'E' )      ! eigenrays
       CALL WriteRay3D( alpha, beta, is )   ! produces no output if NR=1
    CASE ( 'A', 'a' ) ! arrivals
       rayt2 = ray3D( is )%x - ray3D( is - 1 )%x ! ray tangent !!! does this always need to be done???
       RcvrDeclAngle = RadDeg * ATAN2( rayt2( 3 ), NORM2( rayt2( 1 : 2 ) ) )
       RcvrAzimAngle = RadDeg * ATAN2( rayt2( 2 ), rayt2( 1 ) )

       CALL AddArr3D( omega, itheta, iz, ir, Amp, phaseInt, delay, &
            SrcDeclAngle, SrcAzimAngle, RcvrDeclAngle, RcvrAzimAngle, &
            ray3D( is )%NumTopBnc, ray3D( is )%NumBotBnc )
    CASE ( 'C'  )     ! coherent TL
       U = U + CMPLX( Amp * EXP( -i * ( omega * delay - phaseInt ) ) )
    CASE DEFAULT      ! incoherent/semi-coherent TL
       SELECT CASE( Beam%Type( 1 : 1 ) )
       CASE ( 'B', 'b' )   ! Gaussian beam
          U = U + SNGL( 2. * pi * ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
       CASE DEFAULT
          U = U + SNGL(           ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
       END SELECT
    END SELECT

  END SUBROUTINE ApplyContribution

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

  SUBROUTINE ScalePressure3D( Dalpha, Dbeta, c, epsilon, P, Ntheta, Nrz, Nr, RunType, freq )
    !! Scale the pressure field

    INTEGER,            INTENT( IN    ) :: Ntheta, Nrz, Nr
    REAL    ( KIND=8 ), INTENT( IN    ) :: Dalpha, Dbeta        ! angular spacing between rays
    REAL    ( KIND=8 ), INTENT( IN    ) :: freq, c              ! source frequency, nominal sound speed
    COMPLEX,            INTENT( INOUT ) :: P( Ntheta, Nrz, Nr ) ! Pressure field
    COMPLEX ( KIND=8 ), INTENT( IN    ) :: epsilon( 2 )
    CHARACTER (LEN=5 ), INTENT( IN    ) :: RunType
    COMPLEX ( KIND=8 )                  :: const
!!!! this routine should be eliminated

    ! Compute scale factor for field
    SELECT CASE ( RunType( 2 : 2 ) )
    CASE ( 'C' )   ! Cerveny Gaussian beams in Cartesian coordinates
       ! epsilon is normally imaginary here, so const is complex
       const = SQRT( epsilon( 1 ) * epsilon( 2 ) ) * freq * Dbeta * Dalpha / ( SQRT( c ) ) **3
           ! put this factor into the beam instead?
       P( :, :, : ) = CMPLX( const, KIND=4 ) * P( :, :, : )
    CASE DEFAULT
       const = 1.0
    END SELECT

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

  END SUBROUTINE ScalePressure3D

END MODULE Influence3D
