Step3D Subroutine

public subroutine Step3D(ray0, ray2, topRefl, botRefl)

! what if we cross isegx, isegy, or isegz at the same time?

Arguments

Type IntentOptional Attributes Name
type(ray3DPt), intent(inout) :: ray0
type(ray3DPt), intent(inout) :: ray2
logical, intent(out) :: topRefl
logical, intent(out) :: botRefl

Calls

proc~~step3d~~CallsGraph proc~step3d Step3D interface~cross_product cross_product proc~step3d->interface~cross_product proc~computedeltapq ComputeDeltaPQ proc~step3d->proc~computedeltapq proc~evaluatessp3d EvaluateSSP3D proc~step3d->proc~evaluatessp3d proc~raynormal RayNormal proc~step3d->proc~raynormal proc~raynormal_unit RayNormal_unit proc~step3d->proc~raynormal_unit proc~reducestep3d ReduceStep3D proc~step3d->proc~reducestep3d proc~steptobdry3d StepToBdry3D proc~step3d->proc~steptobdry3d proc~updateraypq UpdateRayPQ proc~step3d->proc~updateraypq 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~analytic3d Analytic3D proc~evaluatessp3d->proc~analytic3d proc~ccubic cCubic proc~evaluatessp3d->proc~ccubic proc~clinear cLinear proc~evaluatessp3d->proc~clinear proc~errout ERROUT proc~evaluatessp3d->proc~errout proc~hexahedral Hexahedral proc~evaluatessp3d->proc~hexahedral proc~n2linear n2Linear proc~evaluatessp3d->proc~n2linear proc~readssp ReadSSP proc~ccubic->proc~readssp proc~splineall SPLINEALL proc~ccubic->proc~splineall proc~updatedepthsegmentt UpdateDepthSegmentT proc~ccubic->proc~updatedepthsegmentt proc~clinear->proc~readssp proc~clinear->proc~updatedepthsegmentt proc~hexahedral->proc~errout interface~monotonic monotonic proc~hexahedral->interface~monotonic proc~hexahedral->proc~readssp proc~update3dxsegmentt Update3DXSegmentT proc~hexahedral->proc~update3dxsegmentt proc~update3dysegmentt Update3DYSegmentT proc~hexahedral->proc~update3dysegmentt proc~update3dzsegmentt Update3DZSegmentT proc~hexahedral->proc~update3dzsegmentt proc~n2linear->proc~readssp proc~n2linear->proc~updatedepthsegmentt proc~monotonic_dble monotonic_dble interface~monotonic->proc~monotonic_dble proc~monotonic_sngl monotonic_sngl interface~monotonic->proc~monotonic_sngl proc~readssp->proc~errout proc~crci CRCI proc~readssp->proc~crci proc~crci->proc~errout proc~franc_garr Franc_Garr proc~crci->proc~franc_garr

Called by

proc~~step3d~~CalledByGraph proc~step3d Step3D proc~traceray3d TraceRay3D proc~traceray3d->proc~step3d proc~bellhopcore BellhopCore proc~bellhopcore->proc~traceray3d program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

  SUBROUTINE Step3D( ray0, ray2, topRefl, botRefl )

    ! Does a single step along the ray
    ! x denotes the ray coordinate, ( x, y, z )
    ! t denotes the scaled tangent to the ray (previously (xi, eta, zeta) )
    ! c * t would be the unit tangent

    ! rays
    TYPE( ray3DPt ), INTENT( INOUT ) :: ray0, ray2
    TYPE( ray3DPt ) :: ray1

    LOGICAL, INTENT( OUT ) :: topRefl, botRefl
    INTEGER         :: iSegx0, iSegy0, iSegz0, snapDim
    REAL  (KIND=8 ) :: gradc0( 3 ), gradc1( 3 ), gradc2( 3 ), &
         c0, cimag0, csq0, cxx0, cyy0, czz0, cxy0, cxz0, cyz0, cnn0, cmn0, cmm0, &
         c1, cimag1, csq1, cxx1, cyy1, czz1, cxy1, cxz1, cyz1, cnn1, cmn1, cmm1, &
         c2, cimag2,       cxx2, cyy2, czz2, cxy2, cxz2, cyz2, c_mat1( 2, 2 ), &
         urayt0( 3 ), urayt1( 3 ), urayt2( 3 ), h, halfh, hw0, hw1, w0, w1, rho
    COMPLEX( KIND=8 ) :: d_phi0, d_phi1
    REAL  (KIND=8 ) :: d_p_tilde0( 2 ), d_p_hat0( 2 ), d_q_tilde0( 2 ), d_q_hat0( 2 ), &
                       d_p_tilde1( 2 ), d_p_hat1( 2 ), d_q_tilde1( 2 ), d_q_hat1( 2 )
    REAL  (KIND=8 ) :: gradcjump( 3 ), csjump, cn1jump, cn2jump, tBdry( 3 ), nBdry( 3 ), RM, R1, R2, Tg, Th, &
         rayt( 3 ), rayn1( 3 ), rayn2( 3 )
    REAL  (KIND=8 ) :: e1( 3 ), e2( 3 )                              ! ray normals for ray-centered coordinates
    REAL  (KIND=8 ) :: p_tilde_in(  2 ), p_hat_in(  2 ), q_tilde_in(  2 ), q_hat_in(  2 ), &
         p_tilde_out( 2 ), p_hat_out( 2 ), RotMat( 2, 2 )

    IF ( STEP_DEBUGGING ) THEN
       WRITE( PRTFile, * )
       WRITE( PRTFile, * ) 'ray0 x t', ray0%x, ray0%t
       WRITE( PRTFile, * ) 'ray0 p q', ray0%p_tilde, ray0%q_tilde
       WRITE( PRTFile, * ) '        ', ray0%p_hat,   ray0%q_hat
       WRITE( PRTFile, * ) 'iSegx iSegy iSegz Top_td_side Bot_td_side', iSegx, iSegy, iSegz, Top_td_side, Bot_td_side
    END IF

    ! The numerical integrator used here is a version of the polygon (a.k.a. midpoint, leapfrog, or Box method), and similar
    ! to the Heun (second order Runge-Kutta method).
    ! However, it's modified to allow for a dynamic step change, while preserving the second-order accuracy).

    ! *** Phase 1 (an Euler step)

    !write( *, * ) '______________________________________________________________'
    !write( *, * ) 'in segment ', ISegBoty
    !write( *, * ) 'current coord ', ray0%x

    CALL EvaluateSSP3D( ray0%x, ray0%t, c0, cimag0, gradc0, cxx0, cyy0, czz0, cxy0, cxz0, cyz0, rho, freq, 'TAB' )
    CALL RayNormal( ray0%t, ray0%phi, c0, e1, e2 ) ! Compute ray normals e1 and e2
    CALL Get_c_partials( cxx0, cxy0, cxz0, cyy0, cyz0, czz0, e1, e2, cnn0, cmn0, cmm0 ) ! Compute second partials of c along ray normals
    CALL ComputeDeltaPQ( ray0, c0, gradc0, cnn0, cmn0, cmm0, d_phi0, d_p_tilde0, d_p_hat0, d_q_tilde0, d_q_hat0)


    iSegx0 = iSegx     ! make note of current layer
    iSegy0 = iSegy
    iSegz0 = iSegz

    csq0   = c0 * c0
    urayt0 = c0 * ray0%t  ! unit tangent
    h = Beam%deltas       ! initially set the step h, to the basic one, deltas

    CALL ReduceStep3D( ray0%x, urayt0, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary
    halfh = 0.5 * h       ! first step of the modified polygon method is a half step

    ray1%x = ray0%x    + halfh * urayt0
    ray1%t = ray0%t    - halfh * gradc0 / csq0
    CALL UpdateRayPQ ( ray1, ray0, halfh, d_phi0, d_p_tilde0, d_p_hat0, d_q_tilde0, d_q_hat0 )

    ! *** Phase 2

    CALL EvaluateSSP3D( ray1%x, ray1%t, c1, cimag1, gradc1, cxx1, cyy1, czz1, cxy1, cxz1, cyz1, rho, freq, 'TAB' )
    ! LP: Fixed; should be ray1%phi; ray2%phi would be uninitialized memory or
    ! left over from the previous ray
    CALL RayNormal( ray1%t, ray1%phi, c1, e1, e2 )
    CALL Get_c_partials( cxx1, cxy1, cxz1, cyy1, cyz1, czz1, e1, e2, cnn1, cmn1, cmm1 ) ! Compute second partials of c along ray normals
    CALL ComputeDeltaPQ( ray1, c1, gradc1, cnn1, cmn1, cmm1, d_phi1, d_p_tilde1, d_p_hat1, d_q_tilde1, d_q_hat1)

    csq1   = c1 * c1
    urayt1 = c1 * ray1%t   ! unit tangent

    CALL ReduceStep3D( ray0%x, urayt1, iSegx0, iSegy0, iSegz0, h ) ! reduce h to land on boundary

    ! use blend of f' based on proportion of a full step used.
    w1  = h / ( 2.0d0 * halfh )
    w0  = 1.0d0 - w1
    urayt2 = w0 * urayt0 + w1 * urayt1
    ! Take the blended ray tangent ( urayt2 ) and find the minimum step size ( h )
    ! to put this on a boundary, and ensure that the resulting position
    ! ( ray2%x ) gets put precisely on the boundary.
    CALL StepToBdry3D( ray0%x, ray2%x, urayt2, iSegx0, iSegy0, iSegz0, h, &
       topRefl, botRefl, snapDim )
    !write( *, * ) 'final coord ', ray2%x, w0, w1

    ! Update other variables with this new h
    ! LP: Fixed: ray2%phi now depends on hw0 & hw1 like the other parameters,
    ! originally only depended on h and ray1 vars
    hw0 = h * w0
    hw1 = h * w1
    ray2%t   = ray0%t   - hw0 * gradc0 / csq0 - hw1 * gradc1 / csq1
    ray2%tau = ray0%tau + hw0 / CMPLX( c0, cimag0, KIND=8 ) + hw1 / CMPLX( c1, cimag1, KIND=8 )
    CALL UpdateRayPQ ( ray2, ray0, hw0, d_phi0, d_p_tilde0, d_p_hat0, d_q_tilde0, d_q_hat0 )
    CALL UpdateRayPQ ( ray2, ray2, hw1, d_phi1, d_p_tilde1, d_p_hat1, d_q_tilde1, d_q_hat1 ) ! Not a typo, accumulating into 2

    IF ( STEP_DEBUGGING ) &
       WRITE( PRTFile, * ) 'ray2%t', ray2%t

    ray2%Amp       = ray0%Amp
    ray2%Phase     = ray0%Phase
    ray2%NumTopBnc = ray0%NumTopBnc
    ray2%NumBotBnc = ray0%NumBotBnc

    ! *** If we crossed an interface, apply jump condition ***

    CALL EvaluateSSP3D( ray2%x, ray2%t, c2, cimag2, gradc2, cxx2, cyy2, czz2, cxy2, cxz2, cyz2, rho, freq, 'TAB' )

    ray2%c = c2

    IF ( iSegx /= iSegx0 .OR. &
         iSegy /= iSegy0 .OR. &
         iSegz /= iSegz0 ) THEN

       gradcjump =  gradc2 - gradc0  ! this is precise only for c-linear layers

!!! what if we cross isegx, isegy, or isegz at the same time?
       IF ( iSegz /= iSegz0 ) THEN
          nBdry = [ 0D0, 0D0, -SIGN( 1.0D0, ray2%t( 3 ) ) ]   ! inward normal to layer
       ELSE IF ( iSegx /= iSegx0 ) THEN
          nBdry = [ -SIGN( 1.0D0, ray2%t( 1 ) ), 0D0, 0D0 ]   ! inward normal to x-segment
       ELSE
          nBdry = [ 0D0, -SIGN( 1.0D0, ray2%t( 2 ) ), 0D0 ]   ! inward normal to y-segment
       END IF
       CALL CurvatureCorrection

    END IF

    ! WRITE( PRTFile, * ) 'ray2 p q', ray2%p_tilde, ray2%q_tilde
    ! WRITE( PRTFile, * ) '        ', ray2%p_hat,   ray2%q_hat

  CONTAINS
    SUBROUTINE CurvatureCorrection

      ! correct p-q due to jumps in the gradient of the sound speed

      Th    = DOT_PRODUCT( ray2%t, nBdry )   ! component of ray tangent, normal to boundary
      tBdry = ray2%t - Th * nBdry            ! tangent, along the boundary, in the reflection plane
      tBdry = tBdry / NORM2( tBdry )         ! unit boundary tangent
      Tg    = DOT_PRODUCT( ray2%t, tBdry )   ! component of ray tangent, along the boundary

      rayt = c2 * ray2%t                     ! unit tangent to ray

      rayn2 = cross_product( rayt, nBdry )   ! ray tangent x boundary normal gives refl. plane normal
      rayn2 = rayn2 / NORM2( rayn2 )         ! unit normal
      rayn1 = -cross_product( rayt, rayn2 )  ! ray tangent x refl. plane normal is first ray normal

      ! normal and tangential derivatives of the sound speed
      cn1jump = DOT_PRODUCT( gradcjump, rayn1 )
      cn2jump = DOT_PRODUCT( gradcjump, rayn2 )
      csjump  = DOT_PRODUCT( gradcjump, rayt  )

      ! WRITE( PRTFile, * ) 'cn1 cn2 cs jumps', cn1jump, cn2jump, csjump

      RM = Tg / Th   ! this is tan( alpha ) where alpha is the angle of incidence
      R1 = RM * ( 2 * cn1jump - RM * csjump ) / c2 ** 2
      R2 = RM * cn2jump / c2 ** 2

      ! *** curvature correction ***

      CALL RayNormal_unit( rayt, ray2%phi, e1, e2 )   ! Compute ray normals e1 and e2

      RotMat( 1, 1 ) = DOT_PRODUCT( rayn1, e1 )
      RotMat( 1, 2 ) = DOT_PRODUCT( rayn1, e2 )
      RotMat( 2, 1 ) = -RotMat( 1, 2 )             ! DOT_PRODUCT( rayn2, e1 )
      RotMat( 2, 2 ) = DOT_PRODUCT( rayn2, e2 )

      ! rotate p-q values in e1, e2 system, onto rayn1, rayn2 system

      p_tilde_in = RotMat( 1, 1 ) * ray2%p_tilde + RotMat( 1, 2 ) * ray2%p_hat
      p_hat_in   = RotMat( 2, 1 ) * ray2%p_tilde + RotMat( 2, 2 ) * ray2%p_hat

      q_tilde_in = RotMat( 1, 1 ) * ray2%q_tilde + RotMat( 1, 2 ) * ray2%q_hat
      q_hat_in   = RotMat( 2, 1 ) * ray2%q_tilde + RotMat( 2, 2 ) * ray2%q_hat

      ! here's the actual curvature change

      p_tilde_out = p_tilde_in - q_tilde_in * R1 - q_hat_in * R2
      p_hat_out   = p_hat_in   - q_tilde_in * R2

      ! rotate p back to e1, e2 system, q does not change
      ! Note RotMat^(-1) = RotMat^T

      ray2%p_tilde = RotMat( 1, 1 ) * p_tilde_out + RotMat( 2, 1 ) * p_hat_out
      ray2%p_hat   = RotMat( 1, 2 ) * p_tilde_out + RotMat( 2, 2 ) * p_hat_out

    END SUBROUTINE CurvatureCorrection

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

    SUBROUTINE Get_c_partials( cxx, cxy, cxz, cyy, cyz, czz, e1, e2, cnn, cmn, cmm )

      ! Computes the second partials of c along ray normals

      REAL (KIND=8), INTENT( IN  ) :: cxx, cxy, cxz, cyy, cyz, czz  ! curvature of sound speed (cartesian)
      REAL (KIND=8), INTENT( IN  ) :: e1( 3 ), e2( 3 )              ! principal normals
      REAL (KIND=8), INTENT( OUT ) :: cnn, cmn, cmm                 ! curvature of sound speed (ray-centered)

      cnn = cxx * e1( 1 )**2 + cyy * e1( 2 )**2 + czz * e1( 3 )**2 + 2.0 * cxy * e1( 1 ) * e1( 2 ) + &
           2.0 * cxz * e1( 1 ) * e1( 3 ) + 2.0 * cyz * e1( 2 ) * e1( 3 )

      cmn = cxx * e1( 1 ) * e2( 1 ) + cyy * e1( 2 ) * e2( 2 ) + czz * e1( 3 ) * e2( 3 ) + &
           cxy * ( e1( 1 ) * e2( 2 ) + e2( 1 ) * e1( 2 ) ) + cxz * ( e1( 1 ) * e2( 3 ) + e2( 1 ) * e1( 3 ) ) +  &
           cyz * ( e1( 2 ) * e2( 3 ) + e2( 2 ) * e1( 3 ) )

      cmm = cxx * e2( 1 )**2 + cyy * e2( 2 )**2 + czz * e2( 3 )**2 + 2.0 * cxy * e2( 1 ) * e2( 2 ) + &
           2.0 * cxz * e2( 1 ) * e2( 3 ) + 2.0 * cyz * e2( 2 ) * e2( 3 )

      RETURN
    END SUBROUTINE Get_c_partials
  END SUBROUTINE Step3D