InfluenceGeoHatRayCen Subroutine

public subroutine InfluenceGeoHatRayCen(U, alpha, dalpha)

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

Arguments

Type IntentOptional Attributes Name
complex, intent(inout) :: U(NRz_per_range,Pos%NRr)
real(kind=8), intent(in) :: alpha
real(kind=8), intent(in) :: dalpha

Calls

proc~~influencegeohatraycen~~CallsGraph proc~influencegeohatraycen InfluenceGeoHatRayCen proc~applycontribution~2 ApplyContribution proc~influencegeohatraycen->proc~applycontribution~2 proc~finalphase FinalPhase proc~influencegeohatraycen->proc~finalphase proc~incphaseifcaustic IncPhaseIfCaustic proc~influencegeohatraycen->proc~incphaseifcaustic proc~rtoir RToIR proc~influencegeohatraycen->proc~rtoir proc~addarr AddArr proc~applycontribution~2->proc~addarr proc~writeray2d WriteRay2D proc~applycontribution~2->proc~writeray2d proc~writeray3d WriteRay3D proc~applycontribution~2->proc~writeray3d sngl sngl proc~applycontribution~2->sngl proc~isatcaustic IsAtCaustic proc~finalphase->proc~isatcaustic proc~incphaseifcaustic->proc~isatcaustic proc~addarr->sngl

Called by

proc~~influencegeohatraycen~~CalledByGraph proc~influencegeohatraycen InfluenceGeoHatRayCen proc~bellhopcore BellhopCore proc~bellhopcore->proc~influencegeohatraycen proc~bellhopcore~2 BellhopCore proc~bellhopcore~2->proc~influencegeohatraycen program~bellhop BELLHOP program~bellhop->proc~bellhopcore~2 program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

  SUBROUTINE InfluenceGeoHatRayCen( U, alpha, dalpha )
  !! Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates

    REAL (KIND=8), INTENT( IN    ) :: alpha, dalpha                 ! take-off angle
    COMPLEX,       INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )   ! complex pressure field
    INTEGER          :: irA, irB, II
    REAL    (KIND=8) :: nA, nB, zr, L, dq( Beam%Nsteps - 1 )
    REAL    (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps ),  RcvrDeclAngleV ( Beam%Nsteps )
    COMPLEX (KIND=8) :: dtau( Beam%Nsteps - 1 )
    REAL    (KIND=8) :: KMAHphase( Beam%Nsteps )

    ! need to add logic related to NRz_per_range

    ! LP: See discussion of this change in the readme.
    qOld = ray2D( 1 )%q( 1 )
    phase = 0
    KMAHphase( 1 ) = 0
    DO is = 2, Beam%Nsteps
       q = ray2D( is )%q( 1 )
       CALL IncPhaseIfCaustic( .TRUE. )
       qOld = q
       KMAHphase( is ) = phase
    END DO

    q0           = ray2D( 1 )%c / Dalpha   ! Reference for J = q0 / q
    SrcDeclAngle = RadDeg * alpha          ! take-off angle in degrees

    dq   = ray2D( 2 : Beam%Nsteps )%q( 1 ) - ray2D( 1 : Beam%Nsteps - 1 )%q( 1 )
    dtau = ray2D( 2 : Beam%Nsteps )%tau    - ray2D( 1 : Beam%Nsteps - 1 )%tau

    ! pre-calculate ray normal based on tangent with c(s) scaling
    znV = -ray2D( 1 : Beam%Nsteps )%t( 1 ) * ray2D( 1 : Beam%Nsteps )%c
    rnV =  ray2D( 1 : Beam%Nsteps )%t( 2 ) * ray2D( 1 : Beam%Nsteps )%c

    RcvrDeclAngleV( 1 : Beam%Nsteps ) = RadDeg * ATAN2( ray2D( 1 : Beam%Nsteps )%t( 2 ), ray2D( 1 : Beam%Nsteps )%t( 1 ) )

    ! During reflection imag(q) is constant and adjacent normals cannot bracket a segment of the TL
    ! line, so no special treatment is necessary

    IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) )  ! point source

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

    RcvrDepths: DO iz = 1, NRz_per_range
       zR = Pos%Rz( iz )

       IF ( ABS( znV( 1 ) ) < 1D-6 ) THEN   ! normal parallel to horizontal receiver line
          nA  = 1D10
          rA  = 1D10
          irA = 1
       ELSE
          nA  = ( zR - ray2D( 1 )%x( 2 )   ) / znV( 1 )
          rA  = ray2D( 1 )%x( 1 ) + nA * rnV( 1 )
          irA = RToIR( rA )
       END IF

       Stepping: DO iS = 2, Beam%Nsteps

          ! Compute ray-centered coordinates, (znV, rnV)

          IF ( ABS( znV( iS ) ) < 1D-10 ) CYCLE Stepping   ! If normal parallel to TL-line, skip to next step on ray
          nB  = ( zR - ray2D( iS )%x( 2 )   ) / znV( iS )
          rB  = ray2D( iS )%x( 1 ) + nB * rnV( iS )
          irB = RToIR( rB )

          ! detect and skip duplicate points (happens at boundary reflection)
          IF ( ABS( ray2D( iS )%x( 1 ) - ray2D( iS - 1 )%x( 1 ) ) < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) .OR. irA == irB ) THEN
             rA  = rB
             nA  = nB
             irA = irB
             CYCLE Stepping
          END IF

          RcvrDeclAngle = RcvrDeclAngleV( iS )

          ! *** Compute contributions to bracketed receivers ***

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

          RcvrRanges: DO ir = irA + 1 - II, irB + II, SIGN( 1, irB - irA )  ! Compute influence for each rcvr
             W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
             n = ABS( nA                + W * ( nB - nA ) )
             q = ray2D( iS - 1 )%q( 1 ) + W * dq( iS - 1 )     ! interpolated amplitude
             L = ABS( q ) / q0   ! beam radius

             IF ( n < L ) THEN   ! in beamwindow?
                delay    = ray2D( iS - 1 )%tau + W * dtau( iS - 1 )   ! interpolated delay
                const    = ray2D( iS )%Amp / SQRT( ABS( q ) )
                W        = ( L - n ) / L   ! hat function: 1 on center, 0 on edge
                Amp      = const * W

                phase = KMAHphase( iS - 1 )
                qOld = ray2D( is - 1 )%q( 1 )
                CALL FinalPhase( .FALSE. )

                CALL ApplyContribution( U( iz, ir ) )
             END IF
          END DO RcvrRanges
          rA  = rB
          nA  = nB
          irA = irB
       END DO Stepping
    END DO RcvrDepths

  END SUBROUTINE InfluenceGeoHatRayCen