bellhop3D.f90 Source File
Provides BELLHOP3D program definition
This file depends on
Source Code
!! Provides BELLHOP3D program definition PROGRAM BELLHOP3D !! Gaussian beam tracing for ocean acoustics in three dimensions ! Gaussian beam tracing in three dimensions ! Michael B. Porter ! Copyright (C) 2009 Michael B. Porter ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! You should have received a copy of the GNU General Public License ! along with this program. If not, see <http://www.gnu.org/licenses/>. ! Original version done at NRL in 1986 ! Published at SACLANT Undersea Research Center 1989 ! ! Converted to Fortran 2003 and greatly enhanced for outside users in 2010 ! with the support of the Agency for Defense Development, Republic of Korea ! Supported also by the U.S. Office of Naval Research ! ! Changes included: ! Use of standard BELLHOP input files ! Inclusion of multiple sources and receivers ! Integrated option for Nx2D or full 3D ! Implementation of several standard beam options for both the Nx2D and 3D cases ! Reading in: ! oceanography from a SSP file ! bathymetry or altimetry from a BTY or ATI file ! a reflection coefficient from a BRC or TRC file ! a source beam pattern from a SBP file ! Loose ends: ! Nx2D version does not handle jumps in cx, cy (routine step2d) ! cannot specify isingle( 2 ) for alpha and beta ! Trilinear interpolation (hex) ignores cxy values. ! If the water depth for the lower half space is much larger than that for the SSP ! in the water column, it will use a step that is too large and exit the ray box. ! Cerveny beams (rarely used): ! influenceC calls SSP; need to select EvaluateSSP2D or EvaluateSSP3D for that to work in BELLHOP3D ! efficiency changes for Cerveny beams as well (and in BELLHOP) ! space filling or minimum width options are applied to both alpha and beta--- often only want space filling in azimuth ! Influend3DGeoHat writes no eigenray info if number of receiver ranges NR=1 ! r( 1 ) = 1 m in BELLHOP plus logic for reversing ! geogaussian should pre-calculate dtauds and dqds like geohat? ! fix automatic deltas selection ! detect Sz below bottom with bty file ! Should probably use a Structure of Arrays instead of an Array of Structures for speed ! Desired additional features: ! Terrain following option for receiver depth ! Variable bottom type vs.lat/long USE bellhopMod USE ReadEnvironmentBell USE RefCoef USE bdry3DMod USE BeamPattern USE sspMod USE influence USE Influence3D USE FatalError IMPLICIT NONE CHARACTER ( LEN=80 ) :: FileRoot ThreeD = .TRUE. ! get the file root for naming all input and output files ! should add some checks here ... CALL GET_COMMAND_ARGUMENT( 1, FileRoot ) ! Open the print file OPEN( UNIT = PRTFile, FILE = TRIM( FileRoot ) // '.prt', STATUS = 'UNKNOWN', IOSTAT = iostat ) ! Read in control data CALL ReadEnvironment( FileRoot, ThreeD ) CALL ReadATI3D( FileRoot, Bdry%Top%HS%Opt( 5 : 5 ), Bdry%Top%HS%Depth, PRTFile ) ! AlTImetry CALL ReadBTY3D( FileRoot, Bdry%Bot%HS%Opt( 2 : 2 ), Bdry%Bot%HS%Depth, PRTFile ) ! BaThYmetry CALL ReadReflectionCoefficient( FileRoot, Bdry%Bot%HS%Opt( 1 : 1 ), Bdry%Top%HS%Opt( 2 : 2 ), PRTFile ) ! (top and bottom) SBPFlag = Beam%RunType( 3 : 3 ) CALL ReadPAT( FileRoot, PRTFile ) ! Source Beam Pattern CALL OpenOutputFiles( FileRoot, ThreeD ) CALL BellhopCore CONTAINS ! **********************************************************************! SUBROUTINE BellhopCore !! Core subroutine to run Bellhop algorithm USE angleMod USE SourceReceiverPositions USE ArrMod USE WriteRay ! @note "LP" ! Was 400M; this translates to 16 GB of arrivals, which is half my computer's ! memory, so I can't run this and the C++/CUDA version in parallel. This large ! a size should not be necessary for most runs anyway, and it's not great for ! a program to allocate this kind of memory just as a default (rather than in ! response to the user specifying a very large problem size). So, cut in half. ! @endnote INTEGER, PARAMETER :: ArrivalsStorage = 200000000 INTEGER :: IBPvec( 1 ), ibp, iBeamWindow2, irz, itheta, isx, isy, isz, iRec, ir REAL ( KIND=8 ) :: Tstart, Tstop REAL ( KIND=8 ) :: Amp0, RadMax, S REAL ( KIND=8 ) :: c0, cimag0, gradc( 3 ), cxx, cyy, czz, cxy, cxz, cyz, rho REAL ( KIND=8 ) :: tdummy( 3 ) REAL ( KIND=8 ), ALLOCATABLE :: x_rcvrMat( :, :, : ), t_rcvr( :, : ) COMPLEX ( KIND=8 ) :: epsilon( 2 ) COMPLEX, ALLOCATABLE :: P( :, :, : ), U( :, : ) CALL CPU_TIME( Tstart ) omega = 2.0 * pi * freq IF ( Beam%deltas == 0.0 ) Beam%deltas = ( Bdry%Bot%HS%Depth - Bdry%Top%HS%Depth ) / 10.0 ! Automatic step size selection Angles%alpha = DegRad * Angles%alpha ! convert to radians Angles%Dalpha = 0.0 IF ( Angles%Nalpha /= 1 ) & Angles%Dalpha = ( Angles%alpha( Angles%Nalpha ) - Angles%alpha( 1 ) ) / ( Angles%Nalpha - 1 ) ! angular spacing between beams SELECT CASE ( Beam%RunType( 5 : 5 ) ) CASE ( 'I' ) NRz_per_range = 1 ! irregular grid CASE ( 'R' ) NRz_per_range = Pos%NRz ! rectilinear grid END SELECT ! for a TL calculation, allocate space for the pressure matrix SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'C', 'S', 'I' ) ! TL calculation ALLOCATE ( P( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat ) ALLOCATE ( U( NRz_per_range, Pos%NRr ), Stat = IAllocStat ) ! used for 2D option IF ( IAllocStat /= 0 ) & CALL ERROUT( 'BELLHOP', 'Insufficient memory for TL matrix: reduce Nr * NRz' ) CASE ( 'A', 'a', 'R', 'E' ) ! Arrivals calculation ALLOCATE ( P( 1, 1, 1 ), Stat = IAllocStat ) ! open a dummy variable ALLOCATE ( U( 1, 1 ), Stat = IAllocStat ) ! open a dummy variable END SELECT ! for an arrivals run, allocate space for arrivals matrices SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'A', 'a' ) MaxNArr = MAX( ArrivalsStorage / ( Pos%Ntheta * NRz_per_range * Pos%NRr ), 10 ) ! allow space for at least 10 arrivals WRITE( PRTFile, * ) WRITE( PRTFile, * ) '( Maximum # of arrivals = ', MaxNArr, ')' ALLOCATE ( Arr3D( Pos%Ntheta, NRz_per_range, Pos%NRr, MaxNArr ), & NArr3D( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat ) IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', & 'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' ) ! For a 2D Arrivals run, also need to allocate a structure for that IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN ALLOCATE ( Arr( NRz_per_range, Pos%NRr, MaxNArr ), & NArr( NRz_per_range, Pos%NRr ), Stat = IAllocStat ) IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', & 'Insufficient memory to allocate arrivals matrix; reduce parameter ArrivalsStorage' ) NArr( 1 : NRz_per_range, 1 : Pos%NRr ) = 0 END IF CASE DEFAULT MaxNArr = 1 ALLOCATE ( Arr3D( Pos%Ntheta, NRz_per_range, Pos%NRr, 1 ), NArr3D( Pos%Ntheta, NRz_per_range, Pos%NRr ), Stat = IAllocStat ) END SELECT ALLOCATE( x_rcvrMat( 2, Pos%Ntheta, Pos%NRr ), t_rcvr( 2, Pos%Ntheta ), Stat = IAllocStat ) IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP', & 'Insufficient memory to x_rcvrMat; reduce Nr or Ntheta' ) ! 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 ) ) WRITE( PRTFile, * ) Source_z: DO isz = 1, Pos%NSz ! loop over source z-coordinate Source_x: DO isx = 1, Pos%Nsx ! loop over source x-coordinate Source_y: DO isy = 1, Pos%Nsy ! loop over source y-coordinate P = 0.0 ! Zero out field matrix U = 0.0 NArr3D = 0 ! LP: Zero out 3D arrival matrix ! IF ( r( 1 ) == 0.0 ) r( 1 ) = 1.0 xs_3D = [ Pos%sx( isx ), Pos%sy( isy ), Pos%sz( isz ) ] WRITE( PRTFile, * ) WRITE( PRTFile, "( 'xs = ', G11.3, 2X, G11.3, 2X, G11.3 )" ) xs_3D ! positions of rcvrs in the x-y plane; this is pre-calculated for InfluenceGeoHatCart ! It is not clear that the pre-calculation saves time ... DO ir = 1, Pos%NRr DO itheta = 1, Pos%Ntheta x_rcvrMat( 1 : 2, itheta, ir ) = xs_3D( 1 : 2 ) + Pos%Rr( ir ) * t_rcvr( :, itheta ) ! x-y coordinate of the receiver END DO END DO ! *** Compute 'optimal' beam constant *** tdummy = [ 0.0, 0.0, 1.0 ] CALL EvaluateSSP3D( xs_3D, tdummy, c0, cimag0, gradc, cxx, cyy, czz, cxy, cxz, cyz, rho, freq, 'TAB' ) ray2D( 1 )%c = c0 ray3D( 1 )%c = c0 CALL PickEpsilon( Beam%Type( 1 : 2 ), omega, c0, Angles%Dalpha, Angles%Dbeta, Beam%rLoop, Beam%epsMultiplier, epsilon ) ! beam constant ! *** Trace successive beams *** AzimuthalAngle: DO ibeta = 1, Angles%Nbeta ! this is also the receiver bearing angle for a 2D run SrcAzimAngle = RadDeg * Angles%beta( ibeta ) ! take-off azimuthal angle in degrees !if ( ibeta /= 134 ) cycle AzimuthalAngle IF ( Angles%iSingle_beta == 0 .OR. ibeta == Angles%iSingle_beta ) THEN ! Single beam run? !IF ( ibeta == 2 ) THEN ! Single beam run? !IF ( mod( ibeta+1, 2 ) == 0 ) THEN ! Single beam run? WRITE( PRTFile, FMT = "( 'Tracing azimuthal beam ', I4, F10.2 )" ) ibeta, SrcAzimAngle FLUSH( PRTFile ) ! WRITE( *, FMT = "( 'Tracing beam ', I4, F10.2 )" ) ibeta, RadDeg * Angles%beta( ibeta ) DeclinationAngle: DO ialpha = 1, Angles%Nalpha SrcDeclAngle = RadDeg *Angles%alpha( ialpha ) ! take-off declination angle in degrees IF ( Angles%iSingle_alpha == 0 .OR. ialpha == Angles%iSingle_alpha ) THEN ! Single beam run? !IF ( ialpha == 11 .or. ialpha == 12 ) THEN ! Single beam run? IF ( Angles%iSingle_alpha > 0 .OR. Angles%iSingle_beta > 0 .OR. & Angles%Nalpha * Angles%Nbeta <= 20000 .OR. MOD( ialpha, 20 ) == 1 ) THEN WRITE( PRTFile, FMT = "( ' Tracing declination beam ', I4, F10.2 )" ) ialpha, SrcDeclAngle END IF !flush( prtfile ) IBPvec = maxloc( SrcBmPat( :, 1 ), mask = SrcBmPat( :, 1 ) < SrcDeclAngle ) ! index of ray angle in beam pattern IBP = IBPvec( 1 ) IBP = MAX( IBP, 1 ) ! don't go before beginning of table IBP = MIN( IBP, NSBPPts - 1 ) ! don't go past end of table ! linear interpolation to get amplitude s = ( SrcDeclAngle - SrcBmPat( IBP, 1 ) ) / ( SrcBmPat( IBP + 1, 1 ) - SrcBmPat( IBP, 1 ) ) Amp0 = ( 1 - s ) * SrcBmPat( IBP, 2 ) + s * SrcBmPat( IBP + 1, 2 ) ! Lloyd mirror pattern for semi-coherent option IF ( Beam%RunType( 1 : 1 ) == 'S' ) & Amp0 = Amp0 * SQRT( 2.0 ) * ABS( SIN( omega / c0 * xs_3D( 3 ) * SIN( Angles%alpha( ialpha ) ) ) ) SELECT CASE ( Beam%RunType( 6 : 6 ) ) ! flag for 2D or 3D calculation CASE ( '2' ) ! Nx2D calculation, neglecting horizontal refraction CALL TraceRay2D( Angles%alpha( ialpha ), Angles%beta( ibeta ), Amp0 ) IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN ! If not a ray trace run, calculate the field SELECT CASE ( Beam%Type( 1 : 1 ) ) CASE ( 'R' ) iBeamWindow2 = Beam%iBeamWindow ** 2 RadMax = 50 * c0 / freq ! 50 wavelength max radius CALL InfluenceCervenyRayCen( U, epsilon( 1 ), Angles%alpha( ialpha ), IBeamWindow2, RadMax ) CASE ( 'C' ) CALL ERROUT( 'BELLHOP3D', 'Run Type ''C'' not supported at this time' ) iBeamWindow2 = Beam%iBeamWindow ** 2 RadMax = 50 * c0 / freq ! 50 wavelength max radius CALL InfluenceCervenyCart( U, epsilon( 1 ), Angles%alpha( ialpha ), IBeamWindow2, RadMax ) CASE ( 'g' ) CALL InfluenceGeoHatRayCen( U, Angles%alpha( ialpha ), Angles%Dalpha ) CASE ( 'S' ) RadMax = 50 * c0 / freq ! LP: Was missing / uninitialized. CALL InfluenceSGB( U, Angles%alpha( ialpha ), Angles%Dalpha, RadMax ) CASE ( 'B' ) CALL InfluenceGeoGaussianCart( U, Angles%alpha( ialpha ), Angles%Dalpha ) CASE ( 'G', '^', ' ' ) CALL InfluenceGeoHatCart( U, Angles%alpha( ialpha ), Angles%Dalpha ) CASE DEFAULT CALL ERROUT( 'BELLHOP3D', 'Invalid Run Type' ) END SELECT END IF CASE ( '3' ) ! full 3D calculation CALL TraceRay3D( Angles%alpha( ialpha ), Angles%beta( ibeta ), epsilon, Amp0 ) IF ( Beam%RunType( 1 : 1 ) /= 'R' ) THEN ! If not a ray trace run, calculate the field SELECT CASE ( Beam%RunType( 2 : 2 ) ) CASE ( 'C' ) ! Cerveny style beams CALL ERROUT( 'BELLHOP3D', 'Run Type ''C'' not supported at this time' ) ! option: assemble f, g, h from p-q !ray3D( 1 : Beam%Nsteps )%f = ray3D( 1 : Beam%Nsteps )%p_tilde( 1 ) * & ! ray3D( 1 : Beam%Nsteps )%q_hat( 2 ) - & ! ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%p_hat( 1 ) !ray3D( 1 : Beam%Nsteps )%g = -( ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%q_hat( 1 ) - & ! ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) * & ! ray3D( 1 : Beam%Nsteps )%p_hat( 2 ) ) !ray3D( 1 : Beam%Nsteps )%h = ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%q_hat( 2 ) - & ! ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%p_hat( 2 ) !ray3D( 1 : Beam%Nsteps )%DetP = ray3D( 1 : Beam%Nsteps )%p_tilde( 1 ) * & ! ray3D( 1 : Beam%Nsteps )%p_hat( 2 ) - & ! ray3D( 1 : Beam%Nsteps )%p_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%p_hat( 1 ) !ray3D( 1 : Beam%Nsteps )%DetQ = ray3D( 1 : Beam%Nsteps )%q_tilde( 1 ) * & ! ray3D( 1 : Beam%Nsteps )%q_hat( 2 ) - & ! ray3D( 1 : Beam%Nsteps )%q_tilde( 2 ) * & ! ray3D( 1 : Beam%Nsteps )%q_hat( 1 ) !CALL Influence3D_3D( xs_3D, Angles%alpha( ialpha ), iBeamWindow, P ) CASE ( 'g' ) ! Geometric-beams with hat-shape CALL Influence3DGeoHatRayCen( Angles%alpha( ialpha ), Angles%beta( ibeta ), & Angles%Dalpha, Angles%Dbeta, P ) CASE ( 'G', '^', ' ' ) ! Geometric-beams with hat-shape in Cartesian coordinates CALL Influence3DGeoHatCart( Angles%alpha( ialpha ), Angles%beta( ibeta ), & Angles%Dalpha, Angles%Dbeta, P, x_rcvrMat, t_rcvr ) CASE ( 'b' ) ! Geometric-beams with Gaussian-shape CALL Influence3DGeoGaussianRayCen( Angles%alpha( ialpha ), Angles%beta( ibeta ), & Angles%Dalpha, Angles%Dbeta, P ) CASE ( 'B' ) ! Geometric-beams with Gaussian-shape in Cartesian coordiantes CALL Influence3DGeoGaussianCart( Angles%alpha( ialpha ), Angles%beta( ibeta ), & Angles%Dalpha, Angles%Dbeta, P, x_rcvrMat, t_rcvr ) CASE DEFAULT CALL ERROUT( 'BELLHOP3D', 'Invalid Run Type' ) END SELECT END IF END SELECT ! Optionally dump rays to a disk file IF ( Beam%RunType( 1 : 1 ) == 'R' ) THEN CALL WriteRay3D( Angles%alpha( ialpha ), Angles%beta( ibeta ), Beam%Nsteps ) ENDIF ENDIF ! closes iSingle test END DO DeclinationAngle ! for a 2D TL run, scale the pressure and copy the 2D slice into the radial of the 3D field IF ( Beam%RunType( 6 : 6 ) == '2' ) THEN ! 2D calculation SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'C', 'S', 'I' ) ! TL calculation CALL ScalePressure( Angles%Dalpha, ray2D( 1 )%c, Pos%Rr, U, NRz_per_range, Pos%NRr, Beam%RunType, freq ) P( ibeta, :, : ) = U U = 0 ! clear out the pressure field on the radial in prep for next radial CASE ( 'A' ) ! arrivals calculation, ascii NArr3D( ibeta, :, : ) = NArr( :, : ) Arr3D( ibeta, :, :, : ) = Arr( :, :, : ) Arr3D( ibeta, :, :, : )%SrcAzimAngle = SNGL( SrcAzimAngle ) ! angle Arr3D( ibeta, :, :, : )%RcvrAzimAngle = SNGL( SrcAzimAngle ) ! angle (rcvr angle is same as source angle) Narr = 0 ! this clears out the 2D arrival structure CASE ( 'a' ) ! arrivals calculation, binary NArr3D( ibeta, :, : ) = NArr( :, : ) Arr3D( ibeta, :, :, : ) = Arr( :, :, : ) Arr3D( ibeta, :, :, : )%SrcAzimAngle = SNGL( SrcAzimAngle ) ! angle Arr3D( ibeta, :, :, : )%RcvrAzimAngle = SNGL( SrcAzimAngle ) ! angle (rcvr angle is same as source angle) Narr = 0 ! this clears out the 2D arrival structure END SELECT END IF ENDIF ! closes iSingle test END DO AzimuthalAngle ! *** Scale the complex pressure field *** IF ( Beam%RunType( 6 : 6 ) == '3' ) THEN ! 3D calculation SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'C', 'S', 'I' ) ! TL calculation CALL ScalePressure3D( Angles%Dalpha, Angles%Dbeta, ray2D( 1 )%c, epsilon, P, & Pos%Ntheta, NRz_per_range, Pos%NRr, Beam%RunType, freq ) END SELECT END IF ! Write out the field SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'C', 'S', 'I' ) ! TL calculation DO irz = 1, Pos%NRz RcvrBearing: DO itheta = 1, Pos%Ntheta iRec = 10 + ( isx - 1 ) * Pos%Nsy * Pos%Ntheta * Pos%NSz * Pos%NRz + & ( isy - 1 ) * Pos%Ntheta * Pos%NSz * Pos%NRz + & ( itheta - 1 ) * Pos%NSz * Pos%NRz + & ( isz - 1 ) * Pos%NRz + irz WRITE( SHDFile, REC = IRec ) P( itheta, irz, 1 : Pos%NRr ) END DO RcvrBearing END DO CASE ( 'A' ) ! arrivals calculation, ascii CALL WriteArrivalsASCII3D( Pos%Rr, Pos%Ntheta, NRz_per_range, Pos%NRr ) CASE ( 'a' ) ! arrivals calculation, binary CALL WriteArrivalsBinary3D( Pos%Rr, Pos%Ntheta, NRz_per_range, Pos%NRr ) END SELECT END DO Source_y END DO Source_x END DO Source_z ! close all files SELECT CASE ( Beam%RunType( 1 : 1 ) ) CASE ( 'C', 'S', 'I' ) ! TL calculation CLOSE( SHDFile ) CASE ( 'A', 'a' ) ! arrivals calculation CLOSE( ARRFile ) CASE ( 'R' ) ! ray trace CLOSE( RAYFile ) END SELECT ! Display run time CALL CPU_TIME( Tstop ) WRITE( PRTFile, "( /, ' CPU Time = ', G15.3, 's' )" ) Tstop - Tstart END SUBROUTINE BellhopCore ! **********************************************************************! SUBROUTINE PickEpsilon( BeamType, omega, c, Dalpha, Dbeta, rLoop, EpsMultiplier, epsilon ) ! Picks the optimum value for epsilon REAL (KIND=8), INTENT( IN ) :: omega, c ! angular frequency, sound speed REAL (KIND=8), INTENT( IN ) :: Dalpha, Dbeta ! angular spacing for ray fan REAL (KIND=8), INTENT( IN ) :: epsMultiplier, Rloop ! multiplier, loop range COMPLEX (KIND=8), INTENT( OUT ) :: epsilon( 2 ) ! beam initial conditions CHARACTER (LEN= 2), INTENT( IN ) :: BeamType LOGICAL, SAVE :: INIFlag = .TRUE. REAL (KIND=8), SAVE :: HalfWidth( 2 ) = [ 0.0, 0.0 ] COMPLEX (KIND=8) :: epsilonOpt( 2 ) CHARACTER (LEN=80) :: TAG ! LP: BUG: Multiple codepaths do not set epsilonOpt, leads to UB SELECT CASE ( BeamType( 1 : 1 ) ) CASE ( 'C', 'R' ) ! Cerveny beams TAG = 'Cerveny style beam' SELECT CASE ( BeamType( 2 : 2 ) ) CASE ( 'F' ) TAG = 'Space filling beams' HalfWidth( 1 ) = 2.0 / ( ( omega / c ) * Dalpha ) HalfWidth( 2 ) = 0.0 IF ( Dbeta /= 0.0 ) HalfWidth( 2 ) = 2.0 / ( ( omega / c ) * Dbeta ) epsilonOpt = i * 0.5 * omega * HalfWidth ** 2 CASE ( 'M' ) TAG = 'Minimum width beams' HalfWidth( 1 ) = SQRT( 2.0 * c * 1000.0 * rLoop / omega ) HalfWidth( 2 ) = HalfWidth( 1 ) epsilonOPT = i * 0.5 * omega * HalfWidth **2 CASE ( 'C' ) TAG = 'Cerveny style beam' END SELECT CASE ( 'g' ) TAG = 'Geometric beam, hat-shaped, Ray coord.' epsilonOPT = 0.0 CASE ( 'G', '^' ) TAG = 'Geometric beam, hat-shaped, Cart. coord.' epsilonOPT = 0.0 CASE ( 'b' ) TAG = 'Geometric beam, Gaussian-shaped, Ray coord.' epsilonOPT = 0.0 CASE ( 'B' ) TAG = 'Geometric beam, Gaussian-shaped, Cart. coord.' epsilonOPT = 0.0 CASE ( 'S' ) TAG = 'Simple Gaussian beams' halfwidth = 2.0 / ( ( omega / c ) * Dalpha ) epsilonOpt = i * 0.5 * omega * halfwidth ** 2 END SELECT epsilon = epsMultiplier * epsilonOPT ! On first call write info to prt file IF ( INIFlag ) THEN WRITE( PRTFile, * ) WRITE( PRTFile, * ) TAG WRITE( PRTFile, * ) 'HalfWidth1 = ', HalfWidth( 1 ) WRITE( PRTFile, * ) 'HalfWidth2 = ', HalfWidth( 2 ) WRITE( PRTFile, * ) 'epsilonOPT1 = ', epsilonOPT( 1 ) WRITE( PRTFile, * ) 'epsilonOPT2 = ', epsilonOPT( 2 ) WRITE( PRTFile, * ) 'EpsMult = ', EpsMultiplier WRITE( PRTFile, * ) INIFlag = .FALSE. END IF END SUBROUTINE PickEpsilon ! **********************************************************************! SUBROUTINE TraceRay2D( alpha, beta, Amp0 ) ! Traces the beam corresponding to a particular take-off angle USE ReflectMod REAL (KIND=8), INTENT( IN ) :: alpha, beta, Amp0 ! initial angles, amplitude INTEGER :: is, is1 ! index for a step along the ray REAL (KIND=8) :: x( 3 ) ! ray coordinate REAL (KIND=8) :: t_o( 3 ) ! tangent in ocean space !REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho REAL (KIND=8) :: DistBegTop, DistEndTop, DistBegBot, DistEndBot ! Distances from ray beginning, end to top and bottom REAL (KIND=8) :: tinit( 2 ), tradial( 2 ), BotnInt( 3 ), TopnInt( 3 ), s1, s2 REAL (KIND=8) :: z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy REAL (KIND=8) :: term_minx, term_miny, term_maxx, term_maxy LOGICAL :: topRefl, botRefl ! *** Initial conditions *** iSmallStepCtr = 0 tinit = [ COS( alpha ), SIN( alpha ) ] tradial = [ COS( beta ), SIN( beta ) ] ray2D( 1 )%x = [ 0.0D0, xs_3D( 3 ) ] ! CALL EvaluateSSP2D( ray2D( 1 )%x, tinit, c, cimag, gradc, crr, crz, czz, rho, xs_3D, tradial, freq ) ray2D( 1 )%t = tinit / ray2D( 1 )%c ! LP: Set before PickEpsilon independent of ray angles; never overwritten ray2D( 1 )%p = [ 1.0, 0.0 ] ray2D( 1 )%q = [ 0.0, 1.0 ] ray2D( 1 )%tau = 0.0 ray2D( 1 )%Amp = Amp0 ray2D( 1 )%Phase = 0.0 ray2D( 1 )%NumTopBnc = 0 ray2D( 1 )%NumBotBnc = 0 IsegTopx = 1 IsegTopy = 1 IsegBotx = 1 IsegBoty = 1 t_o = RayToOceanT( ray2D( 1 )%t, tradial ) CALL GetTopSeg3D( xs_3D, t_o, .TRUE. ) ! identify the top segment above the source CALL GetBotSeg3D( xs_3D, t_o, .TRUE. ) ! identify the bottom segment below the source ! Trace the beam (note that Reflect alters the step index is) is = 0 CALL Distances3D( xs_3D, Topx, Botx, Topn, Botn, DistBegTop, DistBegBot ) IF ( DistBegTop <= 0 .OR. DistBegBot <= 0 ) THEN Beam%Nsteps = 1 RETURN ! source must be within the medium END IF Stepping: DO WHILE ( .TRUE. ) is = is + 1 is1 = is + 1 CALL Step2D( ray2D( is ), ray2D( is1 ), tradial, topRefl, botRefl ) ! convert polar coordinate of ray to x-y coordinate x = RayToOceanX( ray2D( is1 )%x, xs_3D, tradial ) t_o = RayToOceanT( ray2D( is1 )%t, tradial ) CALL GetTopSeg3D( x, t_o, .FALSE. ) ! identify the top segment above the source CALL GetBotSeg3D( x, t_o, .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( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot ) IF ( topRefl ) THEN IF ( atiType == 'C' ) THEN ! LP: This is superfluous, it's the same x calculated above. x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial ) s1 = ( x( 1 ) - Topx( 1 ) ) / ( xTopSeg( 2 ) - xTopSeg( 1 ) ) ! proportional distance along segment s2 = ( 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( IsegTopx, IsegTopy )%kappa_xx kappa_xy = Top( IsegTopx, IsegTopy )%kappa_xy kappa_yy = Top( IsegTopx, IsegTopy )%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 Reflect2D( is, Bdry%Top%HS, 'TOP', TopnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RTop, NTopPTS, tradial) ray2D( is + 1 )%NumTopBnc = ray2D( is )%NumTopBnc + 1 x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial ) CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot ) ELSE IF ( botRefl ) THEN ! write( *, * ) 'Reflecting', x, Botx ! write( *, * ) 'Botn', Botn ! write( *, * ) 'Distances', DistEndTop, DistEndBot IF ( btyType == 'C' ) THEN x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial ) s1 = ( x( 1 ) - Botx( 1 ) ) / ( xBotSeg( 2 ) - xBotSeg( 1 ) ) ! proportional distance along segment s2 = ( 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 Reflect2D( is, Bdry%Bot%HS, 'BOT', BotnInt, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RBot, NBotPTS, tradial) ray2D( is + 1 )%NumBotBnc = ray2D( is )%NumBotBnc + 1 x = RayToOceanX( ray2D( is + 1 )%x, xs_3D, tradial ) CALL Distances3D( x, Topx, Botx, Topn, Botn, DistEndTop, DistEndBot ) END IF ! Has the ray left the box, lost its energy, escaped the boundaries, or exceeded storage limit? ! LP: See explanation for changes in bdry3DMod: GetTopSeg3D. t_o = RayToOceanT( ray2D( is + 1 )%t, tradial ) 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 ) ) ! LP: The Beam%Box conditions were inexplicably commented out in 2022 revision of Nx2D, see README. IF ( ABS( x( 1 ) - xs_3D( 1 ) ) >= Beam%Box%x .OR. & ABS( x( 2 ) - xs_3D( 2 ) ) >= Beam%Box%y .OR. & ABS( x( 3 ) ) >= Beam%Box%z .OR. & ! LP: Removed xs( 3 ) for consistency x( 1 ) < term_minx .OR. & x( 2 ) < term_miny .OR. & x( 1 ) > term_maxx .OR. & x( 2 ) > term_maxy .OR. & ( x( 1 ) == term_minx .AND. t_o( 1 ) < 0.0 ) .OR. & ( x( 2 ) == term_miny .AND. t_o( 2 ) < 0.0 ) .OR. & ( x( 1 ) == term_maxx .AND. t_o( 1 ) > 0.0 ) .OR. & ( x( 2 ) == term_maxy .AND. t_o( 2 ) > 0.0 ) .OR. & ray2D( is + 1 )%Amp < 0.005 .OR. & ! ray2D( is + 1 )%t( 1 ) < 0 .OR. & ! kills off a backward traveling ray iSmallStepCtr > 50 ) THEN Beam%Nsteps = is + 1 EXIT Stepping ELSE IF ( is >= MaxN - 3 ) THEN WRITE( PRTFile, * ) 'Warning in TraceRay2D : Insufficient storage for ray trajectory' WRITE( PRTFile, * ) 'Angles are alpha = ', alpha * RadDeg, ' beta = ', beta * RadDeg Beam%Nsteps = is EXIT Stepping END IF DistBegTop = DistEndTop DistBegBot = DistEndBot END DO Stepping END SUBROUTINE TraceRay2D ! **********************************************************************! SUBROUTINE Step2D( ray0, ray2, tradial, topRefl, botRefl ) ! Does a single step along the ray ! x denotes the ray coordinate, (r,z) ! t denotes the scaled tangent to the ray (previously (rho, zeta)) ! c * t would be the unit tangent USE Step3DMod REAL (KIND=8), INTENT( IN ) :: tradial( 2 ) ! coordinate of source and ray bearing angle LOGICAL, INTENT( OUT ) :: topRefl, botRefl TYPE( ray2DPt ) :: ray1 TYPE( ray2DPt ), INTENT( IN ) :: ray0 TYPE( ray2DPt ), INTENT( INOUT ) :: ray2 INTEGER :: iSegx0, iSegy0, iSegz0, snapDim REAL (KIND=8 ) :: gradc0( 2 ), gradc1( 2 ), gradc2( 2 ), rho, & c0, cimag0, crr0, crz0, czz0, csq0, cnn0_csq0, & c1, cimag1, crr1, crz1, czz1, csq1, cnn1_csq1, & c2, cimag2, crr2, crz2, czz2, urayt0( 2 ), urayt1( 2 ), urayt2( 2 ), & h, halfh, hw0, hw1, ray2n( 2 ), RM, RN, gradcjump( 2 ), cnjump, csjump, w0, w1, & rayx3D( 3 ), rayt3D( 3 ), x2_o( 3 ), x2_o_out( 3 ) IF ( STEP_DEBUGGING ) THEN WRITE( PRTFile, * ) WRITE( PRTFile, * ) 'ray0 x t', ray0%x, ray0%t WRITE( PRTFile, * ) 'iSegr iSegz', iSegr, iSegz 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) CALL EvaluateSSP2D( ray0%x, ray0%t, c0, cimag0, gradc0, crr0, crz0, czz0, rho, xs_3D, tradial, freq ) csq0 = c0 * c0 cnn0_csq0 = crr0 * ray0%t( 2 )**2 - 2.0 * crz0 * ray0%t( 1 ) * ray0%t( 2 ) + czz0 * ray0%t( 1 )**2 iSegx0 = iSegx ! make note of current layer iSegy0 = iSegy iSegz0 = iSegz h = Beam%deltas ! initially set the step h, to the basic one, deltas urayt0 = c0 * ray0%t ! unit tangent rayx3D = RayToOceanX( ray0%x, xs_3D, tradial ) rayt3D = RayToOceanT( urayt0, tradial ) CALL ReduceStep3D( rayx3D, rayt3D, 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 ray1%p = ray0%p - halfh * cnn0_csq0 * ray0%q ray1%q = ray0%q + halfh * c0 * ray0%p ! *** Phase 2 CALL EvaluateSSP2D( ray1%x, ray1%t, c1, cimag1, gradc1, crr1, crz1, czz1, rho, xs_3D, tradial, freq ) csq1 = c1 * c1 cnn1_csq1 = crr1 * ray1%t( 2 )**2 - 2.0 * crz1 * ray1%t( 1 ) * ray1%t( 2 ) + czz1 * ray1%t( 1 )**2 ! The Munk test case with a horizontally launched ray caused problems. ! The ray vertexes on an interface and can ping-pong around that interface. ! Have to be careful in that case about big changes to the stepsize (that invalidate the leap-frog scheme) in phase II. ! A modified Heun or Box method could also work. urayt1 = c1 * ray1%t ! unit tangent rayt3D = RayToOceanT( urayt1, tradial ) CALL ReduceStep3D( rayx3D, rayt3D, 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 rayt3D = RayToOceanT( urayt2, tradial ) ! 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( rayx3D, x2_o, rayt3D, iSegx0, iSegy0, iSegz0, h, & topRefl, botRefl, snapDim ) ray2%x = OceanToRayX( x2_o, xs_3D, tradial, urayt2, snapDim ) IF ( STEP_DEBUGGING ) THEN x2_o_out = RayToOceanX( ray2%x, xs_3D, tradial ) WRITE( PRTFile, * ) 'OceanToRayX in ', x2_o WRITE( PRTFile, * ) ' ==> ', ray2%x WRITE( PRTFile, * ) ' ==> ', x2_o_out END IF !write( *, * ) 'final coord ', x2_o, ray2%x, w0, w1 hw0 = h * w0 hw1 = h * w1 ray2%t = ray0%t - hw0 * gradc0 / csq0 - hw1 * gradc1 / csq1 ray2%p = ray0%p - hw0 * cnn0_csq0 * ray0%q - hw1 * cnn1_csq1 * ray1%q ray2%q = ray0%q + hw0 * c0 * ray0%p + hw1 * c1 * ray1%p ray2%tau = ray0%tau + hw0 / CMPLX( c0, cimag0, KIND=8 ) + hw1 / CMPLX( c1, cimag1, KIND=8 ) 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 EvaluateSSP2D( ray2%x, ray2%t, c2, cimag2, gradc2, crr2, crz2, czz2, rho, xs_3D, tradial, freq ) ray2%c = c2 !!! this needs modifying like the full 3D version to handle jumps in the x-y direction IF ( iSegz /= iSegz0 ) THEN gradcjump = gradc2 - gradc0 ! this is precise only for c-linear layers ray2n = [ -ray2%t( 2 ), ray2%t( 1 ) ] cnjump = DOT_PRODUCT( gradcjump, ray2n ) csjump = DOT_PRODUCT( gradcjump, ray2%t ) RM = ray2%t( 1 ) / ray2%t( 2 ) RN = RM * ( 2 * cnjump - RM * csjump ) / c2 RN = -RN ray2%p = ray2%p + ray2%q * RN END IF END SUBROUTINE Step2D FUNCTION RayToOceanX( x, xs, tradial ) !! Transform ray coordinates to ocean coordinates REAL ( KIND=8 ) :: RayToOceanX( 3 ) REAL ( KIND=8 ), INTENT( IN ) :: x( 2 ), xs( 3 ), tradial( 2 ) RayToOceanX = [ xs( 1 ) + x( 1 ) * tradial( 1 ), & xs( 2 ) + x( 1 ) * tradial( 2 ), & x( 2 ) ] END FUNCTION RayToOceanX FUNCTION RayToOceanT( t, tradial ) !! Transform ray tangent to ocean tangent coordinates REAL ( KIND=8 ) :: RayToOceanT( 3 ) REAL ( KIND=8 ), INTENT( IN ) :: t( 2 ), tradial( 2 ) RayToOceanT = [ t( 1 ) * tradial( 1 ), & t( 1 ) * tradial( 2 ), & t( 2 ) ] END FUNCTION RayToOceanT FUNCTION OceanToRayX( x, xs, tradial, t, snapDim ) !! Transform ocean coordinates to ray coordinates ! LP: Going back and forth through the coordinate transform won't ! always keep the precise value, so we may have to finesse the floats. ! Valid values of snapDim: ! -2: Snap to X or Y unspecified ! -1: No snap ! 0: Snap to X ! 1: Snap to Y ! 2: Snap to Z REAL ( KIND=8 ) :: OceanToRayX( 2 ) REAL ( KIND=8 ), INTENT( IN ) :: x( 3 ), xs( 3 ), tradial( 2 ), t( 2 ) INTEGER, INTENT( IN ) :: snapDim REAL ( KIND=8 ) :: ret( 2 ), x_back( 3 ), wantdir( 2 ), errdir( 2 ) LOGICAL :: correctdir( 2 ) INTEGER :: i ! Depth always transfers perfectly--not changed. ret( 2 ) = x( 3 ) ! For range, use larger dimension--this avoids divide-by-zero or divide ! by a small number causing accuracy problems. IF ( ABS( tradial( 1 ) ) >= ABS( tradial( 2 ) ) ) THEN ret( 1 ) = ( x( 1 ) - xs( 1 ) ) / tradial( 1 ) ELSE ret( 1 ) = ( x( 2 ) - xs( 2 ) ) / tradial( 2 ) END IF IF ( snapDim < -2 .OR. snapDim == -1 .OR. snapDim >= 2 ) THEN ! Either: ! snapDim out of range (won't happen, but want to help compiler) ! No snap selected--this is the best estimate ! Snap to Z--Z already perfect, this is the best estimate OceanToRayX = ret RETURN END IF ! Only do this iteration a few times, then give up. DO i = 1, 4 ! Go back from 2D to 3D, compare to original x. x_back = RayToOceanX( ret, xs, tradial ) ! If we can't be on the boundary, we want to be slightly forward of ! the boundary, measured in terms of the ray tangent range. This also ! encompasses cases where one component exactly matches (errdir.x_or_y ! == RL(0.0)). For both of these values, only the sign matters. wantdir = tradial * t( 1 ) errdir = x_back( 1 : 2 ) - x( 1 : 2 ) correctdir = ( wantdir * errdir ) >= 0.0D0 IF ( ( snapDim == 0 .AND. correctdir( 1 ) ) .OR. & ( snapDim == 1 .AND. correctdir( 2 ) ) .OR. & ( correctdir( 1 ) .AND. correctdir( 2 ) ) ) THEN OceanToRayX = ret RETURN END IF ! Move to the next floating point value for ret, in the direction ! of the ray tangent. How do we know this will actually produce a ! new value for x_back? If the scale of the floating point steps ! around x_back is larger than that of ret (several values of ret ! map to the same value of x_back after the transform), there ! should be no problem in finding a value that maps exactly, so ! we would not get here. ret( 1 ) = NEAREST( ret( 1 ), t( 1 ) ) END DO WRITE( PRTFile, * ) 'Warning in OceanToRayX: Failed to transform 3D -> 2D -> 3D consistently' OceanToRayX = ret END FUNCTION OceanToRayX ! **********************************************************************! 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 ! **********************************************************************! SUBROUTINE Distances3D( rayx, Topx, Botx, Topn, Botn, DistTop, DistBot ) !! Computes distances from ray to boundaries ! Formula differs from JKPS because code uses outward pointing normals ! Note that Topx, Botx, just need to be any node on the diagonal that divides each square into triangles ! In bdry3DMod, the node is selected as the one at the lowest x, y, index and that defines the triangles REAL (KIND=8), INTENT( IN ) :: rayx( 3 ) ! ray coordinate REAL (KIND=8), INTENT( IN ) :: Topx( 3 ), Botx( 3 ) ! top, bottom boundary coordinate for node REAL (KIND=8), INTENT( IN ) :: Topn( 3 ), Botn( 3 ) ! top, bottom boundary normal REAL (KIND=8), INTENT( OUT ) :: DistTop, DistBot ! distance from the ray to top, bottom boundaries REAL (KIND=8) :: dTop( 3 ), dBot( 3 ) dTop = rayx - Topx ! vector pointing from top to ray dBot = rayx - Botx ! vector pointing from bottom to ray DistTop = -DOT_PRODUCT( Topn, dTop ) DistBot = -DOT_PRODUCT( Botn, dBot ) !!$ write( *, * ) !!$ write( *, * ) 'Distances3D', DistBot !!$ write( *, * ) 'rayx', rayx !!$ write( *, * ) 'Botx', Botx !!$ write( *, * ) 'dBot', dBot !!$ write( *, * ) 'Normal', Botn !!$ write( *, * ) END SUBROUTINE Distances3D END PROGRAM BELLHOP3D