InfluenceCervenyRayCen – The BELLHOP underwater acoustics ray tracing tool

public subroutine InfluenceCervenyRayCen(U, eps, alpha, IBeamWindow2, RadiusMax)

Paraxial (Cerveny-style) beams in ray-centered coordinates

! need to add logic related to NRz_per_range

! This logic means that the first step along the ray is skipped ! which is a problem if deltas is very large, e.g. isospeed problems ! I fixed this in InfluenceGeoHatRayCen

Arguments

Type	Intent		Name
complex,	intent(inout)	::	U(NRz_per_range,Pos%NRr)
complex(kind=8),	intent(in)	::	eps
real(kind=8),	intent(in)	::	alpha
integer,	intent(in)	::	IBeamWindow2
real(kind=8),	intent(in)	::	RadiusMax

Calls

Help

Called by

Help

Source Code

  SUBROUTINE InfluenceCervenyRayCen( U, eps, alpha, iBeamWindow2, RadiusMax )
  !! Paraxial (Cerveny-style) beams in ray-centered coordinates

    INTEGER,          INTENT( IN    ) :: IBeamWindow2
    REAL    (KIND=8), INTENT( IN    ) :: alpha, RadiusMax                ! take-off angle
    COMPLEX,          INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr )  ! complex pressure field
    COMPLEX (KIND=8), INTENT( IN    ) :: eps                          ! LP: EPSILON is an intrinsic
    INTEGER          :: ir1, ir2, KMAHV( MaxN ), KMAH, image
    REAL    (KIND=8) :: nA, nB, nSq, c, zr, Polarity
    REAL    (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps )   ! ray normal
    COMPLEX (KIND=8) :: pVB( MaxN ), qVB( MaxN ), q, epsV( MaxN ), contri, gammaV( MaxN ), gamma, P_n, P_s
    COMPLEX (KIND=8) :: tau

    SELECT CASE ( Beam%RunType( 1 : 1 ) )
    CASE ( 'C', 'I', 'S' )   ! TL
    CASE DEFAULT
       WRITE( PRTFile, * ) 'Cerveny influence does not support eigenrays or arrivals'
       CALL ERROUT( 'InfluenceCervenyRayCen', 'Cerveny influence does not support eigenrays or arrivals' )
    END SELECT

!!! need to add logic related to NRz_per_range

    ! 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%Type( 2 : 2 ) == 'C' ) THEN
       epsV( 1 : Beam%Nsteps ) = i * ABS( ray2D( 1 : Beam%Nsteps )%q( 1 ) / ray2D( 1 : Beam%Nsteps )%q( 2 ) )
    ELSE
       epsV( 1 : Beam%Nsteps ) = eps
    END IF

    pVB(    1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%p( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%p( 2 )
    qVB(    1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%q( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%q( 2 )
    gammaV( 1 : Beam%Nsteps ) = pVB(   1 : Beam%Nsteps ) / qVB( 1 : Beam%Nsteps )

    ! 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

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

    ! compute KMAH index
    ! Following is incorrect for 'Cerveny'-style beamwidth (narrow as possible)
    KMAHV(  1 ) = 1

    DO iS = 2, Beam%Nsteps
       KMAHV(  iS ) = KMAHV( iS - 1 )
       CALL BranchCut( qVB( iS - 1 ), qVB( iS ), Beam%Type, KMAHV( iS ) )
    END DO

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

       Images: DO image = 1, Beam%Nimage

          ! LP: Previous code did rnV = -rnV for image 2 and 3. When Nimage = 2,
          ! this means rnV changes sign every step, which can't possibly be
          ! correct. This was fixed by mbp in InfluenceCervenyCart.
          Polarity = 1.0D0
          IF ( image == 2 ) Polarity = -1.0D0

          !!! This logic means that the first step along the ray is skipped
          !!! which is a problem if deltas is very large, e.g. isospeed problems
          !!! I fixed this in InfluenceGeoHatRayCen
          ir1 = HUGE( ir1 )

          Stepping: DO iS = 2, Beam%Nsteps

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

             ! If normal parallel to TL-line, skip to next step on ray
             ! LP: Changed from TINY( znV( iS ) ), see README.md.
             IF ( ABS( znV( iS ) ) < EPSILON( znV( iS ) ) ) CYCLE Stepping

             SELECT CASE ( image )     ! Images of beams
             CASE ( 1 )                ! True beam
                nB  = ( zR -                             ray2D( iS )%x( 2 )   ) / znV( iS )
             CASE ( 2 )                ! Surface-reflected beam
                nB  = ( zR - ( 2.0 * Bdry%Top%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
             CASE ( 3 )                ! Bottom-reflected beam
                nB  = ( zR - ( 2.0 * Bdry%Bot%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
             CASE DEFAULT
                WRITE( PRTFile, * ) 'Beam%Nimage = ', Beam%Nimage
                CALL ERROUT( 'InfluenceCervenyRayCen', 'Nimage must be 1, 2, or 3' )
             END SELECT

             rB  = ray2D( iS )%x( 1 ) + nB * rnV( iS ) * Polarity
             ir2 = RToIR( rB )

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

             RcvrRanges: DO ir = ir1 + 1, ir2    ! Compute influence for each rcvr
                W     = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
                q     =    qVB( iS - 1 ) + W * (    qVB( iS ) -    qVB( iS - 1 ) )
                gamma = gammaV( iS - 1 ) + W * ( gammaV( iS ) - gammaV( iS - 1 ) )
                n     = nA + W * ( nB - nA )
                nSq   = n * n
                IF ( AIMAG( gamma ) > 0 ) THEN
                   WRITE( PRTFile, * ) 'Unbounded beam'
                   CYCLE RcvrRanges
                END IF

                IF ( -0.5 * omega * AIMAG( gamma ) * nSq < iBeamWindow2 ) THEN   ! Within beam window?
                   c      = ray2D( iS - 1 )%c
                   tau    = ray2D( iS - 1 )%tau + W * ( ray2D( iS )%tau - ray2D( iS - 1 )%tau )
                   contri = ratio1 * ray2D( iS )%Amp * SQRT( c * ABS( epsV( iS ) ) / q ) * &
                        EXP( -i * ( omega * ( tau + 0.5 * gamma * nSq ) - ray2D( iS )%phase ) )

                   SELECT CASE ( Beam%Component )
                   CASE ( 'P' )   ! pressure
                   CASE ( 'V' )   ! vertical component
                      P_n    = -i * omega * gamma * n * contri
                      P_s    = -i * omega / c         * contri
                      contri = c * DOT_PRODUCT( [ P_n, P_s ], ray2D( iS )%t )
                   CASE ( 'H' )   ! horizontal component
                      P_n    = -i * omega * gamma * n * contri
                      P_s    = -i * omega / c         * contri
                      contri = c * ( -P_n * ray2D( iS )%t( 2 ) + P_s * ray2D( iS )%t( 1 ) )
                   END SELECT

                   KMAH = KMAHV( iS - 1 )
                   CALL BranchCut( qVB( iS - 1 ), q, Beam%Type, KMAH ) ! Get correct branch of SQRT

                   IF ( KMAH  < 0  ) contri = -contri
                   contri = Polarity * contri

                   SELECT CASE ( Beam%RunType( 1 : 1 ) )
                   CASE ( 'I', 'S' )   ! Incoherent or Semi-coherent TL
                      contri = ABS( contri ) ** 2
                   END SELECT

                   U( iz, ir ) = U( iz, ir ) + CMPLX( Hermite( n, RadiusMax, 2 * RadiusMax ) * contri )
                END IF
             END DO RcvrRanges
             rA  = rB
             nA  = nB
             ir1 = ir2
          END DO Stepping
       END DO Images
    END DO RcvrDepths

  END SUBROUTINE InfluenceCervenyRayCen

InfluenceCervenyRayCen Subroutine

Contents

Source Code

public subroutine InfluenceCervenyRayCen(U, eps, alpha, IBeamWindow2, RadiusMax)

Arguments

Calls

Called by

Source Code