Influence3DGeoGaussianRayCen Subroutine

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

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

! stent needs to be applied here already

! can generate an inexact result if the resulting integer is too big ! assumes uniform spacing in Pos%Rr ! sign is flipable inside abs $ __ $ $ ! 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

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~~influence3dgeogaussianraycen~~CallsGraph proc~influence3dgeogaussianraycen Influence3DGeoGaussianRayCen interface~cross_product cross_product proc~influence3dgeogaussianraycen->interface~cross_product proc~applycontribution ApplyContribution proc~influence3dgeogaussianraycen->proc~applycontribution proc~raynormal RayNormal proc~influence3dgeogaussianraycen->proc~raynormal proc~rtoir RToIR proc~influence3dgeogaussianraycen->proc~rtoir proc~scalebeam ScaleBeam proc~influence3dgeogaussianraycen->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~~influence3dgeogaussianraycen~~CalledByGraph proc~influence3dgeogaussianraycen Influence3DGeoGaussianRayCen proc~bellhopcore BellhopCore proc~bellhopcore->proc~influence3dgeogaussianraycen program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

  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