Influence3DGeoHatRayCen Subroutine

public subroutine Influence3DGeoHatRayCen(alpha, beta, Dalpha, Dbeta, P)

Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates

! sign is flippable inside abs

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in) :: alpha
real(kind=8), intent(in) :: beta
real(kind=8), intent(in) :: Dalpha
real(kind=8), intent(in) :: Dbeta
complex, intent(out) :: P(Pos%Ntheta,Pos%Nrz,Pos%NRr)

Calls

proc~~influence3dgeohatraycen~~CallsGraph proc~influence3dgeohatraycen Influence3DGeoHatRayCen interface~cross_product cross_product proc~influence3dgeohatraycen->interface~cross_product proc~applycontribution ApplyContribution proc~influence3dgeohatraycen->proc~applycontribution proc~raynormal RayNormal proc~influence3dgeohatraycen->proc~raynormal proc~rtoir RToIR proc~influence3dgeohatraycen->proc~rtoir proc~scalebeam ScaleBeam proc~influence3dgeohatraycen->proc~scalebeam proc~cross_product_dble cross_product_dble interface~cross_product->proc~cross_product_dble proc~cross_product_sngl cross_product_sngl interface~cross_product->proc~cross_product_sngl proc~addarr3d AddArr3D proc~applycontribution->proc~addarr3d proc~writeray3d WriteRay3D proc~applycontribution->proc~writeray3d sngl sngl proc~applycontribution->sngl proc~addarr3d->sngl

Called by

proc~~influence3dgeohatraycen~~CalledByGraph proc~influence3dgeohatraycen Influence3DGeoHatRayCen proc~bellhopcore BellhopCore proc~bellhopcore->proc~influence3dgeohatraycen program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

  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