TraceRay3D Subroutine

subroutine TraceRay3D(alpha, beta, epsilon, Amp0)

Uses

  • proc~~traceray3d~~UsesGraph proc~traceray3d TraceRay3D module~reflect3dmod Reflect3DMod proc~traceray3d->module~reflect3dmod module~step3dmod Step3DMod proc~traceray3d->module~step3dmod module~bellhopmod bellhopMod module~reflect3dmod->module~bellhopmod module~cross_products cross_products module~reflect3dmod->module~cross_products module~raynormals RayNormals module~reflect3dmod->module~raynormals module~bdry3dmod bdry3Dmod module~step3dmod->module~bdry3dmod module~step3dmod->module~bellhopmod module~step3dmod->module~cross_products module~step3dmod->module~raynormals module~sspmod sspmod module~step3dmod->module~sspmod module~fatalerror FatalError module~bdry3dmod->module~fatalerror module~monotonicmod monotonicMod module~bdry3dmod->module~monotonicmod module~subtabulate SubTabulate module~bdry3dmod->module~subtabulate module~mathconstants MathConstants module~bellhopmod->module~mathconstants module~sspmod->module~fatalerror module~sspmod->module~monotonicmod module~splinec splinec module~sspmod->module~splinec

Traces the beam corresponding to a particular take off angle

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in) :: alpha
real(kind=8), intent(in) :: beta
complex(kind=8), intent(in) :: epsilon(2)
real(kind=8), intent(in) :: Amp0

Calls

proc~~traceray3d~~CallsGraph proc~traceray3d TraceRay3D proc~distances3d Distances3D proc~traceray3d->proc~distances3d proc~getbotseg3d GetBotSeg3D proc~traceray3d->proc~getbotseg3d proc~gettopseg3d GetTopSeg3D proc~traceray3d->proc~gettopseg3d proc~reflect3d Reflect3D proc~traceray3d->proc~reflect3d proc~step3d Step3D proc~traceray3d->proc~step3d proc~errout ERROUT proc~getbotseg3d->proc~errout proc~gettopseg3d->proc~errout interface~cross_product cross_product proc~reflect3d->interface~cross_product proc~evaluatessp3d EvaluateSSP3D proc~reflect3d->proc~evaluatessp3d proc~interpolatereflectioncoefficient InterpolateReflectionCoefficient proc~reflect3d->proc~interpolatereflectioncoefficient proc~raynormal RayNormal proc~reflect3d->proc~raynormal proc~step3d->interface~cross_product proc~computedeltapq ComputeDeltaPQ proc~step3d->proc~computedeltapq proc~step3d->proc~evaluatessp3d 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~evaluatessp3d->proc~errout proc~analytic3d Analytic3D proc~evaluatessp3d->proc~analytic3d proc~ccubic cCubic proc~evaluatessp3d->proc~ccubic proc~clinear cLinear proc~evaluatessp3d->proc~clinear 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~~traceray3d~~CalledByGraph proc~traceray3d TraceRay3D proc~bellhopcore BellhopCore proc~bellhopcore->proc~traceray3d program~bellhop3d BELLHOP3D program~bellhop3d->proc~bellhopcore

Source Code

SUBROUTINE TraceRay3D( alpha, beta, epsilon, Amp0 )
!! Traces the beam corresponding to a particular take off angle

  USE Step3DMod
  USE Reflect3DMod

  REAL     ( KIND=8 ), INTENT( IN ) :: Amp0  ! initial amplitude
  REAL     ( KIND=8 ), INTENT( IN ) :: alpha, beta    ! take-off angles of the ray
  COMPLEX  ( KIND=8 ), INTENT( IN ) :: epsilon( 2 )   ! beam initial conditions
  INTEGER             :: is, is1
  REAL     ( KIND=8 ) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot!, &
                         !c, cimag, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho   ! soundspeed derivatives
  REAL     (KIND=8) :: TopnInt( 3 ), BotnInt( 3 )
  REAL     (KIND=8) :: s1, s2
  REAL     (KIND=8) :: z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy
  REAL     (KIND=8) :: tinit( 3 )
  LOGICAL           :: topRefl, botRefl
  REAL     (KIND=8) :: term_minx, term_miny, term_maxx, term_maxy, term_x( 3 ), term_t( 3 )

  ! *** Initial conditions ***

  iSmallStepCtr = 0
  ray3D( 1 )%x    = xs_3D
  tinit = [ COS( alpha ) * COS( beta ), COS( alpha ) * SIN( beta ), SIN( alpha ) ]
  !CALL EvaluateSSP3D(  ray3D( 1 )%x, tinit, c, cimag, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, 'TAB' )

  ray3D( 1 )%t    = tinit / ray3D( 1 )%c ! LP: Set before PickEpsilon independent of ray angles; never overwritten
  !ray3D( 1 )%f    = epsilon( 2 )
  !ray3D( 1 )%g    = epsilon( 1 )
  !ray3D( 1 )%h    = 0.0
  !ray3D( 1 )%DetP = 1.0
  !ray3D( 1 )%DetQ = epsilon( 1 ) * epsilon( 2 )
  !ray3D( 1 )%c    = c   ! LP
  ray3D( 1 )%phi  = 0.0

  ray3D( 1 )%tau       = 0.0
  ray3D( 1 )%Amp       = Amp0
  ray3D( 1 )%Phase     = 0.0
  ray3D( 1 )%NumTopBnc = 0
  ray3D( 1 )%NumBotBnc = 0

  !ray3D( 1 )%p_tilde = [ 1.0, 0.0 ]   use these for complex Cerveny beams
  !ray3D( 1 )%q_tilde = [ epsilon( 1 ), CMPLX( 0.0, KIND=8 ) ]
  !ray3D( 1 )%p_hat   = [ 0.0, 1.0 ]
  !ray3D( 1 )%q_hat   = [ CMPLX( 0.0, KIND=8 ), epsilon( 2 ) ]

  ray3D( 1 )%p_tilde = [ 1.0, 0.0 ]
  ray3D( 1 )%q_tilde = [ 0.0, 0.0 ]
  ray3D( 1 )%p_hat   = [ 0.0, 1.0 ]
  ray3D( 1 )%q_hat   = [ 0.0, 0.0 ]

  ! dummy BotSeg info to force GetBotSeg to search for the active segment on first call
  xTopSeg = [ +big, -big ]
  yTopSeg = [ +big, -big ]
  xBotSeg = [ +big, -big ]
  yBotSeg = [ +big, -big ]

  CALL GetTopSeg3D( xs_3D, ray3D( 1 )%t, .TRUE. )   ! identify the top    segment above the source
  CALL GetBotSeg3D( xs_3D, ray3D( 1 )%t, .TRUE. )   ! identify the bottom segment below the source

  ! Trace the beam (note that Reflect alters the step index is)
  is = 0

  CALL Distances3D( ray3D( 1 )%x, Topx,  Botx, Topn, Botn, DistBegTop, DistBegBot )

  IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN
     Beam%Nsteps = 1
     WRITE( PRTFile, * ) 'Terminating the ray trace because the source is on or outside the boundaries'
     RETURN       ! source must be within the medium
  END IF

  Stepping: DO WHILE ( .TRUE. )
     is  = is + 1
     is1 = is + 1

     CALL Step3D( ray3D( is ), ray3D( is1 ), topRefl, botRefl )

     CALL GetTopSeg3D( ray3D( is1 )%x, ray3D( is1 )%t, .FALSE. )   ! identify the top    segment above the source
     CALL GetBotSeg3D( ray3D( is1 )%x, ray3D( is1 )%t, .FALSE. )   ! identify the bottom segment below the source

     IF ( IsegTopx == 0 .OR. IsegTopy == 0 .OR. IsegBotx == 0 .OR. IsegBoty == 0 ) THEN ! we escaped the box
        Beam%Nsteps = is
        EXIT Stepping
     END IF

     ! Reflections?
     ! Tests that ray at step is is inside, and ray at step is+1 is outside
     ! to detect only a crossing from inside to outside
     ! DistBeg is the distance at step is, which is saved
     ! DistEnd is the distance at step is+1, which needs to be calculated

     CALL Distances3D( ray3D( is1 )%x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     IF ( topRefl ) THEN
        IF ( STEP_DEBUGGING ) &
           WRITE( PRTFile, * ) 'Top reflecting'
        IF ( atiType == 'C' ) THEN
           s1 = ( ray3D( is1 )%x( 1 ) - Topx( 1 ) ) / ( xTopSeg( 2 ) - xTopSeg( 1 ) )   ! proportional distance along segment
           s2 = ( ray3D( is1 )%x( 2 ) - Topx( 2 ) ) / ( yTopSeg( 2 ) - yTopSeg( 1 ) )   ! proportional distance along segment

           TopnInt = Top( IsegTopx,     IsegTopy     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Top( IsegTopx + 1, IsegTopy + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Top( IsegTopx,     IsegTopy + 1 )%Noden * ( 1 - s1 ) * ( s2     )
           z_xx = Top( IsegTopx, IsegTopy )%z_xx
           z_xy = Top( IsegTopx, IsegTopy )%z_xy
           z_yy = Top( IsegTopx, IsegTopy )%z_yy

           kappa_xx = Top( IsegBotx, IsegBoty )%kappa_xx
           kappa_xy = Top( IsegBotx, IsegBoty )%kappa_xy
           kappa_yy = Top( IsegBotx, IsegBoty )%kappa_yy
        ELSE
           TopnInt = Topn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0

           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect3D( is, Bdry%Top%HS, 'TOP', TopnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RTop, NTopPTS )
        ray3D( is + 1 )%NumTopBnc = ray3D( is )%NumTopBnc + 1
        CALL Distances3D( ray3D( is + 1 )%x, Topx,  Botx, Topn, Botn, DistEndTop, DistEndBot )

     ELSE IF ( botRefl ) THEN
        IF ( STEP_DEBUGGING ) &
           WRITE( PRTFile, * ) 'Bottom reflecting'

        IF ( btyType == 'C' ) THEN
           s1 = ( ray3D( is1 )%x( 1 ) - Botx( 1 ) ) / ( xBotSeg( 2 ) - xBotSeg( 1 ) )   ! proportional distance along segment
           s2 = ( ray3D( is1 )%x( 2 ) - Botx( 2 ) ) / ( yBotSeg( 2 ) - yBotSeg( 1 ) )   ! proportional distance along segment

           BotnInt = Bot( IsegBotx,     IsegBoty     )%Noden * ( 1 - s1 ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty     )%Noden * ( s1     ) * ( 1 - s2 ) +  &
                     Bot( IsegBotx + 1, IsegBoty + 1 )%Noden * ( s1     ) * ( s2     ) +  &
                     Bot( IsegBotx,     IsegBoty + 1 )%Noden * ( 1 - s1 ) * ( s2     )
           z_xx = Bot( IsegBotx, IsegBoty )%z_xx
           z_xy = Bot( IsegBotx, IsegBoty )%z_xy
           z_yy = Bot( IsegBotx, IsegBoty )%z_yy

           kappa_xx = Bot( IsegBotx, IsegBoty )%kappa_xx
           kappa_xy = Bot( IsegBotx, IsegBoty )%kappa_xy
           kappa_yy = Bot( IsegBotx, IsegBoty )%kappa_yy
        ELSE
           BotnInt = Botn   ! normal is constant in a segment
           z_xx = 0
           z_xy = 0
           z_yy = 0

           kappa_xx = 0
           kappa_xy = 0
           kappa_yy = 0
        END IF

        CALL Reflect3D( is, Bdry%Bot%HS, 'BOT', BotnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RBot, NBotPTS )
        ray3D( is + 1 )%NumBotBnc = ray3D( is )%NumBotBnc + 1
        CALL Distances3D( ray3D( is + 1 )%x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot )

     END IF

     ! Has the ray exited the beam box, lost its energy, escaped the BTY, ATI boundaries, or exceeded storage limit?
     ! LP: See explanation for changes in bdry3DMod: GetTopSeg3D.
     term_x = ray3D( is + 1 )%x
     term_t = ray3D( is + 1 )%t
     term_minx = MAX( Bot( 1,            1 )%x( 1 ), Top( 1,            1 )%x( 1 ) )
     term_miny = MAX( Bot( 1,            1 )%x( 2 ), Top( 1,            1 )%x( 2 ) )
     term_maxx = MIN( Bot( NBTYPts( 1 ), 1 )%x( 1 ), Top( NATIPts( 1 ), 1 )%x( 1 ) )
     term_maxy = MIN( Bot( 1, NBTYPts( 2 ) )%x( 2 ), Top( 1, NATIPts( 2 ) )%x( 2 ) )
     IF ( &
          ABS( term_x( 1 ) - xs_3D( 1 ) ) >= Beam%Box%x .OR. &
          ABS( term_x( 2 ) - xs_3D( 2 ) ) >= Beam%Box%y .OR. &
          ABS( term_x( 3 )              ) >= Beam%Box%z .OR. & ! box is centered at z=0
          term_x( 1 ) < term_minx .OR. &
          term_x( 2 ) < term_miny .OR. &
          term_x( 1 ) > term_maxx .OR. &
          term_x( 2 ) > term_maxy .OR. &
          ( term_x( 1 ) == term_minx .AND. term_t( 1 ) < 0.0 ) .OR. &
          ( term_x( 2 ) == term_miny .AND. term_t( 2 ) < 0.0 ) .OR. &
          ( term_x( 1 ) == term_maxx .AND. term_t( 1 ) > 0.0 ) .OR. &
          ( term_x( 2 ) == term_maxy .AND. term_t( 2 ) > 0.0 ) .OR. &
          ray3D( is + 1 )%Amp < 0.005 .OR. &
          iSmallStepCtr > 50 ) THEN
        Beam%Nsteps = is + 1
        EXIT Stepping
     ELSE IF ( is >= MaxN - 3 ) THEN
        Beam%Nsteps = is
        WRITE( PRTFile, * ) 'Warning in TraceRay3D : Insufficient storage for ray trajectory'
        WRITE( PRTFile, * ) 'Angles are  alpha = ', alpha * RadDeg, '    beta = ', beta * RadDeg
        EXIT Stepping
     END IF

     DistBegTop = DistEndTop
     DistBegBot = DistEndBot

  END DO Stepping
END SUBROUTINE TraceRay3D