Coverage Report: influence.f90

Generated from GCOV analysis of Fortran source code

35.8%
Lines Executed
441 total lines
60.8%
Branches Executed
525 total branches
100.0%
Calls Executed
10 total calls
0
-
Source:influence.f90
0
-
Graph:influence.gcno
0
-
Data:influence.gcda
0
-
Runs:73
1
-
!! Beam influence computation for pressure field calculation
2
-
3
-
MODULE Influence
4
-
!! Computes beam contributions to complex pressure fields using various beam weighting approaches
5
-
6
-
! Compute the beam influence, i.e. the contribution of a single beam to the complex pressure
7
-
! mbp 12/2018, based on much older subroutines
8
-
9
-
USE bellhopMod
10
-
USE SourceReceiverPositions
11
-
USE ArrMod
12
-
USE sspMod ! used to construct image beams in the Cerveny style beam routines
13
-
USE WriteRay
14
-
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: ERROR_UNIT
15
-
16
-
IMPLICIT NONE
17
-
PUBLIC
18
-
19
-
INTEGER, PRIVATE :: iz, ir, iS
20
-
REAL (KIND=8), PRIVATE :: Ratio1 = 1.0D0 ! scale factor for a line source
21
-
REAL (KIND=8), PRIVATE :: W, s, n, Amp, phase, const, phaseInt, q0, q, qold, RcvrDeclAngle, rA, rB
22
-
COMPLEX (KIND=8), PRIVATE :: delay
23
-
24
-
CONTAINS
25
#####
SUBROUTINE InfluenceCervenyRayCen( U, eps, alpha, iBeamWindow2, RadiusMax )
26
-
!! Paraxial (Cerveny-style) beams in ray-centered coordinates
27
-
28
-
INTEGER, INTENT( IN ) :: IBeamWindow2
29
-
REAL (KIND=8), INTENT( IN ) :: alpha, RadiusMax ! take-off angle
30
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
31
-
COMPLEX (KIND=8), INTENT( IN ) :: eps ! LP: EPSILON is an intrinsic
32
-
INTEGER :: ir1, ir2, KMAHV( MaxN ), KMAH, image
33
-
REAL (KIND=8) :: nA, nB, nSq, c, zr, Polarity
34
#####
REAL (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps ) ! ray normal
35
-
COMPLEX (KIND=8) :: pVB( MaxN ), qVB( MaxN ), q, epsV( MaxN ), contri, gammaV( MaxN ), gamma, P_n, P_s
36
-
COMPLEX (KIND=8) :: tau
37
-
38
#####
SELECT CASE ( Beam%RunType( 1 : 1 ) )
39
-
CASE ( 'C', 'I', 'S' ) ! TL
40
-
CASE DEFAULT
41
#####
WRITE( PRTFile, * ) 'Cerveny influence does not support eigenrays or arrivals'
42
#####
CALL ERROUT( 'InfluenceCervenyRayCen', 'Cerveny influence does not support eigenrays or arrivals' )
43
-
END SELECT
44
-
45
-
!!! need to add logic related to NRz_per_range
46
-
47
-
! During reflection imag(q) is constant and adjacent normals cannot bracket a segment of the TL
48
-
! line, so no special treatment is necessary
49
-
50
#####
IF ( Beam%Type( 2 : 2 ) == 'C' ) THEN
51
#####
epsV( 1 : Beam%Nsteps ) = i * ABS( ray2D( 1 : Beam%Nsteps )%q( 1 ) / ray2D( 1 : Beam%Nsteps )%q( 2 ) )
52
-
ELSE
53
#####
epsV( 1 : Beam%Nsteps ) = eps
54
-
END IF
55
-
56
#####
pVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%p( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%p( 2 )
57
#####
qVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%q( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%q( 2 )
58
#####
gammaV( 1 : Beam%Nsteps ) = pVB( 1 : Beam%Nsteps ) / qVB( 1 : Beam%Nsteps )
59
-
60
-
! pre-calculate ray normal based on tangent with c(s) scaling
61
#####
znV = -ray2D( 1 : Beam%Nsteps )%t( 1 ) * ray2D( 1 : Beam%Nsteps )%c
62
#####
rnV = ray2D( 1 : Beam%Nsteps )%t( 2 ) * ray2D( 1 : Beam%Nsteps )%c
63
-
64
#####
IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) ) ! point source
65
-
66
-
! compute KMAH index
67
-
! Following is incorrect for 'Cerveny'-style beamwidth (narrow as possible)
68
#####
KMAHV( 1 ) = 1
69
-
70
#####
DO iS = 2, Beam%Nsteps
71
#####
KMAHV( iS ) = KMAHV( iS - 1 )
72
#####
CALL BranchCut( qVB( iS - 1 ), qVB( iS ), Beam%Type, KMAHV( iS ) )
73
-
END DO
74
-
75
#####
RcvrDepths: DO iz = 1, NRz_per_range
76
#####
zR = Pos%Rz( iz )
77
-
78
#####
Images: DO image = 1, Beam%Nimage
79
-
80
-
! LP: Previous code did rnV = -rnV for image 2 and 3. When Nimage = 2,
81
-
! this means rnV changes sign every step, which can't possibly be
82
-
! correct. This was fixed by mbp in InfluenceCervenyCart.
83
#####
Polarity = 1.0D0
84
#####
IF ( image == 2 ) Polarity = -1.0D0
85
-
86
-
!!! This logic means that the first step along the ray is skipped
87
-
!!! which is a problem if deltas is very large, e.g. isospeed problems
88
-
!!! I fixed this in InfluenceGeoHatRayCen
89
#####
ir1 = HUGE( ir1 )
90
-
91
#####
Stepping: DO iS = 2, Beam%Nsteps
92
-
93
-
! Compute ray-centered coordinates, (znV, rnV)
94
-
95
-
! If normal parallel to TL-line, skip to next step on ray
96
-
! LP: Changed from TINY( znV( iS ) ), see README.md.
97
#####
IF ( ABS( znV( iS ) ) < EPSILON( znV( iS ) ) ) CYCLE Stepping
98
-
99
#####
SELECT CASE ( image ) ! Images of beams
100
-
CASE ( 1 ) ! True beam
101
#####
nB = ( zR - ray2D( iS )%x( 2 ) ) / znV( iS )
102
-
CASE ( 2 ) ! Surface-reflected beam
103
#####
nB = ( zR - ( 2.0 * Bdry%Top%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
104
-
CASE ( 3 ) ! Bottom-reflected beam
105
#####
nB = ( zR - ( 2.0 * Bdry%Bot%HS%Depth - ray2D( iS )%x( 2 ) ) ) / znV( iS )
106
-
CASE DEFAULT
107
#####
WRITE( PRTFile, * ) 'Beam%Nimage = ', Beam%Nimage
108
#####
CALL ERROUT( 'InfluenceCervenyRayCen', 'Nimage must be 1, 2, or 3' )
109
-
END SELECT
110
-
111
#####
rB = ray2D( iS )%x( 1 ) + nB * rnV( iS ) * Polarity
112
#####
ir2 = RToIR( rB )
113
-
114
-
! detect and skip duplicate points (happens at boundary reflection)
115
#####
IF ( ir1 >= ir2 .OR. &
116
#####
& ABS( ray2D( iS )%x( 1 ) - ray2D( iS - 1 )%x( 1 ) ) < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) THEN
117
#####
rA = rB
118
#####
nA = nB
119
#####
ir1 = ir2
120
#####
CYCLE Stepping
121
-
END IF
122
-
123
#####
RcvrRanges: DO ir = ir1 + 1, ir2 ! Compute influence for each rcvr
124
#####
W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
125
#####
q = qVB( iS - 1 ) + W * ( qVB( iS ) - qVB( iS - 1 ) )
126
#####
gamma = gammaV( iS - 1 ) + W * ( gammaV( iS ) - gammaV( iS - 1 ) )
127
#####
n = nA + W * ( nB - nA )
128
#####
nSq = n * n
129
#####
IF ( AIMAG( gamma ) > 0 ) THEN
130
#####
WRITE( PRTFile, * ) 'Unbounded beam'
131
#####
CYCLE RcvrRanges
132
-
END IF
133
-
134
#####
IF ( -0.5 * omega * AIMAG( gamma ) * nSq < iBeamWindow2 ) THEN ! Within beam window?
135
#####
c = ray2D( iS - 1 )%c
136
#####
tau = ray2D( iS - 1 )%tau + W * ( ray2D( iS )%tau - ray2D( iS - 1 )%tau )
137
#####
contri = ratio1 * ray2D( iS )%Amp * SQRT( c * ABS( epsV( iS ) ) / q ) * &
138
#####
EXP( -i * ( omega * ( tau + 0.5 * gamma * nSq ) - ray2D( iS )%phase ) )
139
-
140
#####
SELECT CASE ( Beam%Component )
141
-
CASE ( 'P' ) ! pressure
142
-
CASE ( 'V' ) ! vertical component
143
#####
P_n = -i * omega * gamma * n * contri
144
#####
P_s = -i * omega / c * contri
145
#####
contri = c * DOT_PRODUCT( [ P_n, P_s ], ray2D( iS )%t )
146
-
CASE ( 'H' ) ! horizontal component
147
#####
P_n = -i * omega * gamma * n * contri
148
#####
P_s = -i * omega / c * contri
149
#####
contri = c * ( -P_n * ray2D( iS )%t( 2 ) + P_s * ray2D( iS )%t( 1 ) )
150
-
CASE DEFAULT
151
#####
WRITE( ERROR_UNIT, * ) 'InfluenceCervenyRayCen: Unknown component: ', Beam%Component
152
#####
CALL ERROUT( 'InfluenceCervenyRayCen', 'Unknown beam component' )
153
-
END SELECT
154
-
155
#####
KMAH = KMAHV( iS - 1 )
156
#####
CALL BranchCut( qVB( iS - 1 ), q, Beam%Type, KMAH ) ! Get correct branch of SQRT
157
-
158
#####
IF ( KMAH < 0 ) contri = -contri
159
#####
contri = Polarity * contri
160
-
161
#####
SELECT CASE ( Beam%RunType( 1 : 1 ) )
162
-
CASE ( 'I', 'S' ) ! Incoherent or Semi-coherent TL
163
#####
contri = ABS( contri ) ** 2
164
-
CASE DEFAULT
165
#####
WRITE( ERROR_UNIT, * ) 'InfluenceCervenyRayCen: Unknown RunType: ', Beam%RunType( 1 : 1 )
166
#####
CALL ERROUT( 'InfluenceCervenyRayCen', 'Unknown RunType' )
167
-
END SELECT
168
-
169
#####
U( iz, ir ) = U( iz, ir ) + CMPLX( Hermite( n, RadiusMax, 2 * RadiusMax ) * contri )
170
-
END IF
171
-
END DO RcvrRanges
172
#####
rA = rB
173
#####
nA = nB
174
#####
ir1 = ir2
175
-
END DO Stepping
176
-
END DO Images
177
-
END DO RcvrDepths
178
-
179
#####
END SUBROUTINE InfluenceCervenyRayCen
180
-
181
-
! **********************************************************************!
182
-
183
#####
SUBROUTINE InfluenceCervenyCart( U, eps, alpha, iBeamWindow2, RadiusMax )
184
-
!! Paraxial (Cerveny-style) beams in Cartesian coordinates
185
-
186
-
INTEGER, INTENT( IN ) :: IBeamWindow2
187
-
REAL (KIND=8), INTENT( IN ) :: alpha, RadiusMax ! take-off angle
188
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
189
-
COMPLEX (KIND=8), INTENT( IN ) :: eps ! LP: EPSILON is an intrinsic
190
-
INTEGER :: KMAHV( MaxN ), KMAH, irA, irB, Image
191
-
REAL (KIND=8), SAVE :: Polarity = 1
192
-
REAL (KIND=8) :: x( 2 ), rayt( 2 ), rayn( 2 ), Tr, Tz, zr, &
193
-
c, cimag, cs, cn, csq, gradc( 2 ), crr, crz, czz, rho, deltaz
194
-
COMPLEX (KIND=8) :: pVB( MaxN ), qVB( MaxN ), q, epsV( MaxN ), contri, gammaV( MaxN ), gamma, const
195
-
COMPLEX (KIND=8) :: tau
196
-
197
#####
SELECT CASE ( Beam%RunType( 1 : 1 ) )
198
-
CASE ( 'C', 'I', 'S' ) ! TL
199
-
CASE DEFAULT
200
#####
WRITE( PRTFile, * ) 'Cerveny influence does not support eigenrays or arrivals'
201
#####
CALL ERROUT( 'InfluenceCervenyRayCen', 'Cerveny influence does not support eigenrays or arrivals' )
202
-
END SELECT
203
-
204
-
! need to add logic related to NRz_per_range
205
-
206
-
! During reflection imag(q) is constant and adjacent normals cannot bracket a segment of the TL
207
-
! line, so no special treatment is necessary
208
-
209
#####
IF ( Beam%Type( 2 : 2 ) == 'C' ) THEN
210
#####
epsV( 1 : Beam%Nsteps ) = i * ABS( ray2D( 1 : Beam%Nsteps )%q( 1 ) / ray2D( 1 : Beam%Nsteps )%q( 2 ) )
211
-
ELSE
212
#####
epsV( 1 : Beam%Nsteps ) = eps
213
-
END IF
214
-
215
#####
pVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%p( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%p( 2 )
216
#####
qVB( 1 : Beam%Nsteps ) = ray2D( 1 : Beam%Nsteps )%q( 1 ) + epsV( 1 : Beam%Nsteps ) * ray2D( 1 : Beam%Nsteps )%q( 2 )
217
-
218
#####
IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) ) ! point source
219
-
220
-
! Form gamma and KMAH index
221
-
! Treatment of KMAH index is incorrect for 'Cerveny' style beam width BeamType
222
-
223
#####
Stepping0: DO iS = 1, Beam%Nsteps
224
-
225
#####
rayt = ray2D( iS )%c * ray2D( iS )%t ! unit tangent
226
#####
rayn = [ rayt( 2 ), -rayt( 1 ) ] ! unit normal
227
-
228
#####
CALL EvaluateSSP( ray2D( iS )%x, ray2D( iS )%t, c, cimag, gradc, crr, crz, czz, rho, Freq, 'TAB' )
229
-
230
#####
csq = c * c
231
#####
cS = DOT_PRODUCT( gradc, rayt )
232
#####
cN = DOT_PRODUCT( gradc, rayn )
233
-
234
#####
Tr = rayt( 1 )
235
#####
Tz = rayt( 2 )
236
-
237
#####
gammaV( iS ) = 0.0
238
#####
IF ( qVB( iS ) /= 0.0 ) THEN
239
#####
gammaV( iS ) = 0.5 * ( pVB( iS ) / qVB( iS ) * Tr**2 + 2.0 * cN / csq * Tz * Tr - cS / csq * Tz**2 )
240
-
END IF
241
-
242
#####
IF ( iS == 1 ) THEN
243
#####
KMAHV( 1 ) = 1
244
-
ELSE
245
#####
KMAHV( iS ) = KMAHV( iS - 1 )
246
#####
CALL BranchCut( qVB( iS - 1 ), qVB( iS ), Beam%Type, KMAHV( iS ) )
247
-
END IF
248
-
END DO Stepping0
249
-
250
#####
Stepping: DO iS = 3, Beam%Nsteps
251
-
! LP: BUG: Assumes rays may never travel left.
252
#####
IF ( ray2D( iS )%x( 1 ) > Pos%Rr( Pos%NRr ) ) RETURN
253
#####
rA = ray2D( iS - 1 )%x( 1 )
254
#####
rB = ray2D( iS )%x( 1 )
255
#####
IF ( ABS( rB - rA ) < 1.0D3 * SPACING( rB ) ) CYCLE Stepping ! don't process duplicate points
256
-
257
#####
irA = RToIR( rA )
258
#####
irB = RToIR( rB )
259
-
260
#####
IF ( irA >= irB ) CYCLE Stepping
261
-
262
#####
RcvrRanges: DO ir = irA + 1, irB
263
-
264
#####
W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
265
-
266
#####
x = ray2D( iS - 1 )%x + W * ( ray2D( iS )%x - ray2D( iS - 1 )%x )
267
#####
rayt = ray2D( iS - 1 )%t + W * ( ray2D( iS )%t - ray2D( iS - 1 )%t )
268
#####
c = ray2D( iS - 1 )%c + W * ( ray2D( iS )%c - ray2D( iS - 1 )%c )
269
#####
q = qVB( iS - 1 ) + W * ( qVB( iS ) - qVB( iS - 1 ) )
270
#####
tau = ray2D( iS - 1 )%tau + W * ( ray2D( iS )%tau - ray2D( iS - 1 )%tau )
271
#####
gamma = gammaV( iS - 1 ) + W * ( gammaV( iS ) - gammaV( iS - 1 ) )
272
-
273
#####
IF ( AIMAG( gamma ) > 0 ) THEN
274
#####
WRITE( PRTFile, * ) 'Unbounded beam'
275
#####
WRITE( PRTFile, * ) gammaV( iS - 1 ), gammaV( iS ), gamma
276
#####
CYCLE RcvrRanges
277
-
END IF
278
-
279
#####
const = Ratio1 * SQRT( c * ABS( epsV( iS - 1 ) ) / q )
280
-
281
-
! Get correct branch of SQRT
282
#####
KMAH = KMAHV( iS - 1 )
283
#####
CALL BranchCut( qVB( iS - 1 ), q, Beam%Type, KMAH )
284
#####
IF ( KMAH < 0 ) const = -const
285
-
286
#####
RcvrDepths: DO iz = 1, NRz_per_range
287
#####
zR = Pos%Rz( iz )
288
-
289
#####
contri = 0.0
290
#####
ImageLoop: DO Image = 1, Beam%Nimage
291
#####
SELECT CASE ( Image )
292
-
CASE ( 1 ) ! True beam
293
#####
deltaz = zR - x( 2 )
294
#####
Polarity = +1.0D0
295
-
CASE ( 2 ) ! Surface reflected beam
296
#####
deltaz = -zR + 2.0 * Bdry%Top%HS%Depth - x( 2 )
297
#####
Polarity = -1.0D0
298
-
CASE ( 3 ) ! Bottom reflected beam
299
#####
deltaz = -zR + 2.0 * Bdry%Bot%HS%Depth - x( 2 )
300
#####
Polarity = +1.0D0 ! assumes rigid bottom
301
-
CASE DEFAULT
302
#####
WRITE( ERROR_UNIT, * ) 'InfluenceCervenyCart: Unknown image index: ', Image
303
#####
CALL ERROUT( 'InfluenceCervenyCart', 'Unknown image index' )
304
-
END SELECT
305
-
306
#####
IF ( omega * AIMAG( gamma ) * deltaz ** 2 < iBeamWindow2 ) THEN
307
#####
contri = contri + Polarity * ray2D( iS )%Amp * Hermite( deltaz, RadiusMax, 2.0 * RadiusMax ) * &
308
#####
EXP( -i * ( omega * ( tau + rayt( 2 ) * deltaz + gamma * deltaz**2 ) - ray2D( iS )%Phase ) )
309
-
END IF
310
-
END DO ImageLoop
311
-
312
-
! contribution to field
313
#####
SELECT CASE( Beam%RunType( 1 : 1 ) )
314
-
CASE ( 'C' ) ! coherent
315
#####
contri = const * contri
316
-
CASE ( 'I', 'S' ) ! incoherent or semi-coherent
317
#####
contri = ABS( const * contri ) ** 2
318
-
CASE DEFAULT
319
#####
WRITE( ERROR_UNIT, * ) 'InfluenceCervenyCart: Unknown RunType: ', Beam%RunType( 1 : 1 )
320
#####
CALL ERROUT( 'InfluenceCervenyCart', 'Unknown RunType' )
321
-
END SELECT
322
#####
U( iz, ir ) = U( iz, ir ) + CMPLX( contri )
323
-
END DO RcvrDepths
324
-
END DO RcvrRanges
325
-
END DO Stepping
326
-
327
-
END SUBROUTINE InfluenceCervenyCart
328
-
329
-
! **********************************************************************!
330
-
331
#####
SUBROUTINE InfluenceGeoHatRayCen( U, alpha, dalpha )
332
-
!! Geometrically-spreading beams with a hat-shaped beam in ray-centered coordinates
333
-
334
-
REAL (KIND=8), INTENT( IN ) :: alpha, dalpha ! take-off angle
335
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
336
-
INTEGER :: irA, irB, II
337
#####
REAL (KIND=8) :: nA, nB, zr, L, dq( Beam%Nsteps - 1 )
338
#####
REAL (KIND=8) :: znV( Beam%Nsteps ), rnV( Beam%Nsteps ), RcvrDeclAngleV ( Beam%Nsteps )
339
#####
COMPLEX (KIND=8) :: dtau( Beam%Nsteps - 1 )
340
#####
REAL (KIND=8) :: KMAHphase( Beam%Nsteps )
341
-
342
-
! need to add logic related to NRz_per_range
343
-
344
-
! LP: See discussion of this change in the readme.
345
#####
qOld = ray2D( 1 )%q( 1 )
346
#####
phase = 0
347
#####
KMAHphase( 1 ) = 0
348
#####
DO is = 2, Beam%Nsteps
349
#####
q = ray2D( is )%q( 1 )
350
#####
CALL IncPhaseIfCaustic( .TRUE. )
351
#####
qOld = q
352
#####
KMAHphase( is ) = phase
353
-
END DO
354
-
355
#####
q0 = ray2D( 1 )%c / Dalpha ! Reference for J = q0 / q
356
#####
SrcDeclAngle = RadDeg * alpha ! take-off angle in degrees
357
-
358
#####
dq = ray2D( 2 : Beam%Nsteps )%q( 1 ) - ray2D( 1 : Beam%Nsteps - 1 )%q( 1 )
359
#####
dtau = ray2D( 2 : Beam%Nsteps )%tau - ray2D( 1 : Beam%Nsteps - 1 )%tau
360
-
361
-
! pre-calculate ray normal based on tangent with c(s) scaling
362
#####
znV = -ray2D( 1 : Beam%Nsteps )%t( 1 ) * ray2D( 1 : Beam%Nsteps )%c
363
#####
rnV = ray2D( 1 : Beam%Nsteps )%t( 2 ) * ray2D( 1 : Beam%Nsteps )%c
364
-
365
#####
RcvrDeclAngleV( 1 : Beam%Nsteps ) = RadDeg * ATAN2( ray2D( 1 : Beam%Nsteps )%t( 2 ), ray2D( 1 : Beam%Nsteps )%t( 1 ) )
366
-
367
-
! During reflection imag(q) is constant and adjacent normals cannot bracket a segment of the TL
368
-
! line, so no special treatment is necessary
369
-
370
#####
IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) ) ! point source
371
-
372
#####
ray2D( 1 : Beam%Nsteps )%Amp = Ratio1 * SQRT( ray2D( 1 : Beam%Nsteps )%c ) * ray2D( 1 : Beam%Nsteps )%Amp
373
-
! pre-apply some scaling
374
-
375
#####
RcvrDepths: DO iz = 1, NRz_per_range
376
#####
zR = Pos%Rz( iz )
377
-
378
#####
IF ( ABS( znV( 1 ) ) < 1D-6 ) THEN ! normal parallel to horizontal receiver line
379
#####
nA = 1D10
380
#####
rA = 1D10
381
#####
irA = 1
382
-
ELSE
383
#####
nA = ( zR - ray2D( 1 )%x( 2 ) ) / znV( 1 )
384
#####
rA = ray2D( 1 )%x( 1 ) + nA * rnV( 1 )
385
#####
irA = RToIR( rA )
386
-
END IF
387
-
388
#####
Stepping: DO iS = 2, Beam%Nsteps
389
-
390
-
! Compute ray-centered coordinates, (znV, rnV)
391
-
392
#####
IF ( ABS( znV( iS ) ) < 1D-10 ) CYCLE Stepping ! If normal parallel to TL-line, skip to next step on ray
393
#####
nB = ( zR - ray2D( iS )%x( 2 ) ) / znV( iS )
394
#####
rB = ray2D( iS )%x( 1 ) + nB * rnV( iS )
395
#####
irB = RToIR( rB )
396
-
397
-
! detect and skip duplicate points (happens at boundary reflection)
398
#####
IF ( ABS( ray2D( iS )%x( 1 ) - ray2D( iS - 1 )%x( 1 ) ) < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) .OR. irA == irB ) THEN
399
#####
rA = rB
400
#####
nA = nB
401
#####
irA = irB
402
#####
CYCLE Stepping
403
-
END IF
404
-
405
#####
RcvrDeclAngle = RcvrDeclAngleV( iS )
406
-
407
-
! *** Compute contributions to bracketed receivers ***
408
-
409
#####
II = 0
410
#####
IF ( irB <= irA ) II = 1 ! going backwards in range
411
-
412
#####
RcvrRanges: DO ir = irA + 1 - II, irB + II, SIGN( 1, irB - irA ) ! Compute influence for each rcvr
413
#####
W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
414
#####
n = ABS( nA + W * ( nB - nA ) )
415
#####
q = ray2D( iS - 1 )%q( 1 ) + W * dq( iS - 1 ) ! interpolated amplitude
416
#####
L = ABS( q ) / q0 ! beam radius
417
-
418
#####
IF ( n < L ) THEN ! in beamwindow?
419
#####
delay = ray2D( iS - 1 )%tau + W * dtau( iS - 1 ) ! interpolated delay
420
#####
const = ray2D( iS )%Amp / SQRT( ABS( q ) )
421
#####
W = ( L - n ) / L ! hat function: 1 on center, 0 on edge
422
#####
Amp = const * W
423
-
424
#####
phase = KMAHphase( iS - 1 )
425
#####
qOld = ray2D( is - 1 )%q( 1 )
426
#####
CALL FinalPhase( .FALSE. )
427
-
428
#####
CALL ApplyContribution( U( iz, ir ) )
429
-
END IF
430
-
END DO RcvrRanges
431
#####
rA = rB
432
#####
nA = nB
433
#####
irA = irB
434
-
END DO Stepping
435
-
END DO RcvrDepths
436
-
437
#####
END SUBROUTINE InfluenceGeoHatRayCen
438
-
439
-
! **********************************************************************!
440
-
441
260537
SUBROUTINE InfluenceGeoHatCart( U, alpha, Dalpha )
442
-
!! Geometric, hat-shaped beams in Cartesisan coordinates
443
-
444
-
REAL (KIND=8), INTENT( IN ) :: alpha, dalpha ! take-off angle, angular spacing
445
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
446
-
INTEGER :: irT( 1 ), irTT
447
260537*
REAL (KIND=8) :: x_ray( 2 ), rayt( 2 ), rayn( 2 ), x_rcvr( 2, NRz_per_range ), rLen, RadiusMax, zMin, zMax, dqds
448
-
COMPLEX (KIND=8) :: dtauds
449
-
450
260537
q0 = ray2D( 1 )%c / Dalpha ! Reference for J = q0 / q
451
260537
SrcDeclAngle = RadDeg * alpha ! take-off angle in degrees
452
260537
phase = 0.0
453
260537
qOld = ray2D( 1 )%q( 1 ) ! used to track KMAH index
454
260537
rA = ray2D( 1 )%x( 1 ) ! range at start of ray
455
-
456
-
! what if never satisfied?
457
-
! what if there is a single receiver (ir = 0 possible)
458
23494879*
irT = MINLOC( Pos%Rr( 1 : Pos%NRr ), MASK = Pos%Rr( 1 : Pos%NRr ) > rA ) ! find index of first receiver to the right of rA
459
260537
ir = irT( 1 )
460
260537
IF ( ray2D( 1 )%t( 1 ) < 0.0d0 .AND. ir > 1 ) ir = ir - 1
461
-
! if ray is left-traveling, get the first receiver to the left of rA
462
-
463
260537
IF ( Beam%RunType( 4 : 4 ) == 'R' ) Ratio1 = SQRT( ABS( COS( alpha ) ) ) ! point source
464
-
465
33151858*
Stepping: DO iS = 2, Beam%Nsteps
466
32891321*
rB = ray2D( iS )%x( 1 )
467
98673963*
x_ray = ray2D( iS - 1 )%x
468
-
469
-
! compute normalized tangent (compute it because we need to measure the step length)
470
98673963*
rayt = ray2D( iS )%x - ray2D( iS - 1 )%x
471
98673963
rlen = NORM2( rayt )
472
32891321*
IF ( rlen < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) CYCLE Stepping
473
-
! if duplicate point in ray, skip to next step along the ray
474
95141808
rayt = rayt / rlen ! unit tangent to ray
475
95141808
rayn = [ -rayt( 2 ), rayt( 1 ) ] ! unit normal to ray
476
31713936
RcvrDeclAngle = RadDeg * ATAN2( rayt( 2 ), rayt( 1 ) )
477
-
478
31713936*
dqds = ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 )
479
31713936*
dtauds = ray2D( iS )%tau - ray2D( iS - 1 )%tau
480
-
481
31713936*
q = ray2D( iS - 1 )%q( 1 )
482
31713936
CALL IncPhaseIfCaustic( .TRUE. )
483
31713936
qold = q
484
-
485
31713936*
RadiusMax = MAX( ABS( ray2D( iS - 1 )%q( 1 ) ), ABS( ray2D( iS )%q( 1 ) ) ) / q0 / ABS( rayt( 1 ) )
486
-
! beam radius projected onto vertical line
487
-
488
-
! depth limits of beam
489
31713936
IF ( ABS( rayt( 1 ) ) > 0.5 ) THEN ! shallow angle ray
490
27293053*
zmin = min( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) - RadiusMax
491
27293053*
zmax = max( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) + RadiusMax
492
-
ELSE ! steep angle ray
493
4420883
zmin = -HUGE( zmin )
494
4420883
zmax = +HUGE( zmax )
495
-
END IF
496
-
497
-
! compute beam influence for this segment of the ray
498
6923948
RcvrRanges: DO
499
-
! is Rr( ir ) contained in [ rA, rB )? Then compute beam influence
500
38637884*
IF ( Pos%Rr( ir ) >= MIN( rA, rB ) .AND. Pos%Rr( ir ) < MAX( rA, rB ) ) THEN
501
-
502
392346710*
x_rcvr( 1, 1 : NRz_per_range ) = Pos%Rr( ir )
503
6992325
IF ( Beam%RunType( 5 : 5 ) == 'I' ) THEN
504
#####
x_rcvr( 2, 1 ) = Pos%Rz( ir ) ! irregular grid
505
-
ELSE
506
392346710*
x_rcvr( 2, 1 : NRz_per_range ) = Pos%Rz( 1 : NRz_per_range ) ! rectilinear grid
507
-
END IF
508
-
509
392346710*
RcvrDepths: DO iz = 1, NRz_per_range
510
385354385*
IF ( x_rcvr( 2, iz ) < zmin .OR. x_rcvr( 2, iz ) > zmax ) CYCLE RcvrDepths
511
-
512
497242470*
s = DOT_PRODUCT( x_rcvr( :, iz ) - x_ray, rayt ) / rlen ! proportional distance along ray
513
497242470*
n = ABS( DOT_PRODUCT( x_rcvr( :, iz ) - x_ray, rayn ) ) ! normal distance to ray
514
165747490*
q = ray2D( iS - 1 )%q( 1 ) + s * dqds ! interpolated amplitude
515
165747490
RadiusMax = ABS( q / q0 ) ! beam radius
516
-
517
172739815
IF ( n < RadiusMax ) THEN
518
6373209*
delay = ray2D( iS - 1 )%tau + s * dtauds ! interpolated delay
519
6373209*
const = Ratio1 * SQRT( ray2D( iS )%c / ABS( q ) ) * ray2D( iS )%Amp
520
6373209
W = ( RadiusMax - n ) / RadiusMax ! hat function: 1 on center, 0 on edge
521
6373209
Amp = const * W
522
6373209
CALL FinalPhase( .FALSE. )
523
-
524
6373209*
CALL ApplyContribution( U( iz, ir ) )
525
-
END IF
526
-
END DO RcvrDepths
527
-
END IF
528
-
529
-
! bump receiver index, ir, towards rB
530
38637884*
IF ( Pos%Rr( ir ) < rB ) THEN
531
12488329
IF ( ir >= Pos%NRr ) EXIT RcvrRanges ! go to next step on ray
532
11068368
irTT = ir + 1 ! bump right
533
11068368*
IF ( Pos%Rr( irTT ) >= rB ) EXIT RcvrRanges
534
-
ELSE
535
26149555
IF ( ir <= 1 ) EXIT RcvrRanges ! go to next step on ray
536
685768
irTT = ir - 1 ! bump left
537
685768*
IF ( Pos%Rr( irTT ) <= rB ) EXIT RcvrRanges
538
-
END IF
539
6923948
ir = irTT
540
-
END DO RcvrRanges
541
-
542
31974473
rA = rB
543
-
END DO Stepping
544
-
545
260537
END SUBROUTINE InfluenceGeoHatCart
546
-
547
-
! **********************************************************************!
548
-
549
9200
SUBROUTINE InfluenceGeoGaussianCart( U, alpha, Dalpha )
550
-
!! Geometric, Gaussian beams in Cartesian coordintes
551
-
552
-
INTEGER, PARAMETER :: BeamWindow = 4 ! beam window: kills beams outside e**(-0.5 * ibwin**2 )
553
-
REAL (KIND=8), INTENT( IN ) :: alpha, dalpha ! take-off angle, angular spacing
554
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
555
-
INTEGER :: irT( 1 ), irTT
556
-
REAL (KIND=8) :: x_ray( 2 ), rayt( 2 ), rayn( 2 ), x_rcvr( 2 ), rLen, RadiusMax, zMin, zMax, sigma, lambda, A, dqds
557
-
COMPLEX (KIND=8) :: dtauds
558
-
559
9200
q0 = ray2D( 1 )%c / Dalpha ! Reference for J = q0 / q
560
9200
SrcDeclAngle = RadDeg * alpha ! take-off angle in degrees
561
9200
phase = 0
562
9200
qOld = ray2D( 1 )%q( 1 ) ! used to track KMAH index
563
9200
rA = ray2D( 1 )%x( 1 ) ! range at start of ray
564
-
565
-
! what if never satisfied?
566
-
! what if there is a single receiver (ir = 0 possible)
567
-
568
4627600*
irT = MINLOC( Pos%Rr( 1 : Pos%NRr ), MASK = Pos%Rr( 1 : Pos%NRr ) > rA )
569
-
! find index of first receiver to the right of rA
570
9200
ir = irT( 1 )
571
-
572
9200*
IF ( ray2D( 1 )%t( 1 ) < 0.0d0 .AND. ir > 1 ) ir = ir - 1
573
-
! if ray is left-traveling, get the first receiver to the left of rA
574
-
575
-
! sqrt( 2 * pi ) represents a sum of Gaussians in free space
576
9200
IF ( Beam%RunType( 4 : 4 ) == 'R' ) THEN
577
9200
Ratio1 = SQRT( ABS( COS( alpha ) ) ) / SQRT( 2. * pi ) ! point source
578
-
ELSE
579
#####
Ratio1 = 1 / SQRT( 2. * pi ) ! line source
580
-
END IF
581
-
582
2315871*
Stepping: DO iS = 2, Beam%Nsteps
583
-
584
2306671*
rB = ray2D( iS )%x( 1 )
585
6920013*
x_ray = ray2D( iS - 1 )%x
586
-
587
-
! compute normalized tangent (compute it because we need to measure the step length)
588
6920013*
rayt = ray2D( iS )%x - ray2D( iS - 1 )%x
589
6920013
rlen = NORM2( rayt )
590
2306671*
IF ( rlen < 1.0D3 * SPACING( ray2D( iS )%x( 1 ) ) ) CYCLE Stepping
591
-
! if duplicate point in ray, skip to next step along the ray
592
6697533
rayt = rayt / rlen
593
6697533
rayn = [ -rayt( 2 ), rayt( 1 ) ] ! unit normal to ray
594
2232511
RcvrDeclAngle = RadDeg * ATAN2( rayt( 2 ), rayt( 1 ) )
595
-
596
2232511*
dqds = ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 )
597
2232511*
dtauds = ray2D( iS )%tau - ray2D( iS - 1 )%tau
598
-
599
2232511*
q = ray2D( iS - 1 )%q( 1 )
600
2232511
CALL IncPhaseIfCaustic( .TRUE. )
601
2232511
qold = q
602
-
603
-
! calculate beam width
604
2232511*
lambda = ray2D( iS - 1 )%c / freq
605
2232511*
sigma = MAX( ABS( ray2D( iS - 1 )%q( 1 ) ), ABS( ray2D( iS )%q( 1 ) ) ) / q0 / ABS( rayt( 1 ) )
606
-
! beam radius projected onto vertical line
607
2232511*
sigma = MAX( sigma, MIN( 0.2 * freq * REAL( ray2D( iS )%tau ), pi * lambda ) )
608
2232511
RadiusMax = BeamWindow * sigma
609
-
610
-
! depth limits of beam
611
-
! LP: For rays shot at exactly 60 degrees, they will hit this edge case.
612
-
! This is a sharp edge--the handling on each side of this edge may be
613
-
! significantly different. So, moved the edge away from the round number.
614
2232511
IF ( ABS( rayt( 1 ) ) > 0.50001 ) THEN ! shallow angle ray
615
1995989*
zmin = min( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) - RadiusMax
616
1995989*
zmax = max( ray2D( iS - 1 )%x( 2 ), ray2D( iS )%x( 2 ) ) + RadiusMax
617
-
ELSE ! steep angle ray
618
236522
zmin = -HUGE( zmin )
619
236522
zmax = +HUGE( zmax )
620
-
END IF
621
-
622
-
! compute beam influence for this segment of the ray
623
1512903
RcvrRanges: DO
624
-
! WRITE( PRTFile, * ) 'iS', iS-2, 'ir', ir-1
625
-
626
-
! is Rr( ir ) contained in [ rA, rB )? Then compute beam influence
627
3745414*
IF ( Pos%Rr( ir ) >= MIN( rA, rB ) .AND. Pos%Rr( ir ) < MAX( rA, rB ) ) THEN
628
-
! WRITE( PRTFile, * ) ' rA', rA, 'rB', rB
629
-
630
916368614*
RcvrDepths: DO iz = 1, NRz_per_range
631
914846407
IF ( Beam%RunType( 5 : 5 ) == 'I' ) THEN
632
#####
x_rcvr = [ Pos%Rr( ir ), Pos%Rz( ir ) ] ! irregular grid
633
-
ELSE
634
2744539221*
x_rcvr = [ Pos%Rr( ir ), Pos%Rz( iz ) ] ! rectilinear grid
635
-
END IF
636
914846407
IF ( x_rcvr( 2 ) < zmin .OR. x_rcvr( 2 ) > zmax ) CYCLE RcvrDepths
637
-
638
339691308
s = DOT_PRODUCT( x_rcvr - x_ray, rayt ) / rlen ! proportional distance along ray
639
339691308
n = ABS( DOT_PRODUCT( x_rcvr - x_ray, rayn ) ) ! normal distance to ray
640
113230436*
q = ray2D( iS - 1 )%q( 1 ) + s * dqds ! interpolated amplitude
641
113230436
sigma = ABS( q / q0 ) ! beam radius
642
113230436*
sigma = MAX( sigma, MIN( 0.2 * freq * REAL( ray2D( iS )%tau ), pi * lambda ) ) ! min pi * lambda, unless near
643
-
644
114752643
IF ( n < BeamWindow * sigma ) THEN ! Within beam window?
645
-
! IF ( ir-1 >= 0 .AND. ir-1 <= 2 ) THEN
646
-
! WRITE( PRTFile, * ) ' iz n sigma', iz-1, n, sigma
647
-
! END IF
648
76021312
A = ABS( q0 / q )
649
76021312*
delay = ray2D( iS - 1 )%tau + s * dtauds ! interpolated delay
650
76021312*
const = Ratio1 * SQRT( ray2D( iS )%c / ABS( q ) ) * ray2D( iS )%Amp
651
76021312
W = EXP( -0.5 * ( n / sigma ) ** 2 ) / ( sigma * A ) ! Gaussian decay
652
76021312
Amp = const * W
653
76021312
CALL FinalPhase( .TRUE. )
654
-
655
76021312*
CALL ApplyContribution( U( iz, ir ) )
656
-
END IF
657
-
END DO RcvrDepths
658
-
END IF
659
-
660
-
! receiver not bracketed; bump receiver index, ir, towards rB
661
3745414*
IF ( rB > Pos%Rr( ir ) ) THEN
662
2917215
IF ( ir >= Pos%NRr ) EXIT RcvrRanges ! go to next step on ray
663
2904528
irTT = ir + 1 ! bump right
664
2904528*
IF ( Pos%Rr( irTT ) >= rB ) EXIT RcvrRanges ! go to next step on ray
665
-
ELSE
666
828199*
IF ( ir <= 1 ) EXIT RcvrRanges ! go to next step on ray
667
828199
irTT = ir - 1 ! bump left
668
828199*
IF ( Pos%Rr( irTT ) <= rB ) EXIT RcvrRanges ! go to next step on ray
669
-
END IF
670
1512903
ir = irTT
671
-
672
-
END DO RcvrRanges
673
-
674
2241711
rA = rB
675
-
END DO Stepping
676
-
677
9200
END SUBROUTINE InfluenceGeoGaussianCart
678
-
679
-
! **********************************************************************!
680
-
681
82394521
SUBROUTINE ApplyContribution( U )
682
-
!! Applies beam contribution to pressure field
683
-
684
-
COMPLEX, INTENT( INOUT ) :: U
685
-
COMPLEX ( KIND=4 ) :: dfield
686
-
687
343
SELECT CASE( Beam%RunType( 1 : 1 ) )
688
-
CASE ( 'E' ) ! eigenrays
689
686
IF ( Title( 1 : 9 ) == 'BELLHOP- ' ) THEN ! BELLHOP run
690
343
CALL WriteRay2D( SrcDeclAngle, iS )
691
-
ELSE ! BELLHOP3D run
692
#####
CALL WriteRay3D( DegRad * SrcDeclAngle, DegRad * SrcAzimAngle, is ) ! produces no output if NR=1
693
-
END IF
694
-
CASE ( 'A', 'a' ) ! arrivals
695
-
CALL AddArr( omega, iz, ir, Amp, phaseInt, delay, SrcDeclAngle, &
696
1430*
& RcvrDeclAngle, ray2D( iS )%NumTopBnc, ray2D( iS )%NumBotBnc )
697
-
CASE ( 'C' ) ! coherent TL
698
82392748
dfield = CMPLX( Amp * EXP( -i * ( omega * delay - phaseInt ) ) )
699
-
! WRITE( PRTFile, * ) 'ApplyContribution dfield', dfield
700
82392748*
U = U + dfield
701
-
! omega * n * n / ( 2 * ray2d( iS )%c**2 * delay ) ) ) ) ! curvature correction
702
-
CASE DEFAULT ! incoherent/semicoherent TL
703
82395951*
IF ( Beam%Type( 1 : 1 ) == 'B' ) THEN ! Gaussian beam
704
#####
U = U + SNGL( SQRT( 2. * pi ) * ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
705
-
ELSE
706
#####
U = U + SNGL( ( const * EXP( AIMAG( omega * delay ) ) ) ** 2 * W )
707
-
END IF
708
-
END SELECT
709
-
710
82394521
END SUBROUTINE ApplyContribution
711
-
712
-
! **********************************************************************!
713
-
714
#####
SUBROUTINE InfluenceSGB( U, alpha, Dalpha, RadiusMax )
715
-
!! Bucker's Simple Gaussian Beams in Cartesian coordinates
716
-
717
-
REAL (KIND=8), INTENT( IN ) :: alpha, dalpha, RadiusMax ! take-off angle, angular spacing
718
-
COMPLEX, INTENT( INOUT ) :: U( NRz_per_range, Pos%NRr ) ! complex pressure field
719
-
REAL (KIND=8) :: x( 2 ), rayt( 2 ), A, beta, cn, CPA, Adeltaz, deltaz, DS, sint, SX1, thet
720
-
COMPLEX (KIND=8) :: contri, tau
721
-
722
-
! LP: Added ABS to match other influence functions. Without it, this will
723
-
! crash (sqrt of negative real) for rays shot backwards.
724
#####
Ratio1 = SQRT( ABS( COS( alpha ) ) )
725
#####
phase = 0
726
#####
qOld = 1.0
727
#####
BETA = 0.98 ! Beam Factor
728
#####
A = -4.0 * LOG( BETA ) / Dalpha**2
729
#####
CN = Dalpha * SQRT( A / pi )
730
#####
rA = ray2D( 1 )%x( 1 )
731
#####
ir = 1
732
-
733
#####
Stepping: DO iS = 2, Beam%Nsteps
734
-
735
#####
RcvrDeclAngle = RadDeg * ATAN2( ray2D( iS )%t( 2 ), ray2D( iS )%t( 1 ) )
736
-
737
#####
rB = ray2D( iS )%x( 1 )
738
-
739
-
! phase shifts at caustics
740
#####
q = ray2D( iS - 1 )%q( 1 )
741
#####
CALL IncPhaseIfCaustic( .FALSE. )
742
#####
qold = q
743
-
744
-
! Loop over bracketed receiver ranges
745
-
! LP: BUG: This way of setting up the loop assumes the ray always travels
746
-
! towards positive R, which is not true for certain bathymetries (or for
747
-
! rays simply shot backwards, which previously would also crash during
748
-
! the setup, see above).
749
#####
RcvrRanges: DO WHILE ( ABS( rB - rA ) > 1.0D3 * SPACING( rA ) .AND. rB > Pos%Rr( ir ) )
750
-
751
#####
W = ( Pos%Rr( ir ) - rA ) / ( rB - rA )
752
#####
x = ray2D( iS - 1 )%x + W * ( ray2D( iS )%x - ray2D( iS - 1 )%x )
753
#####
rayt = ray2D( iS - 1 )%t + W * ( ray2D( iS )%t - ray2D( iS - 1 )%t )
754
#####
q = ray2D( iS - 1 )%q( 1 ) + W * ( ray2D( iS )%q( 1 ) - ray2D( iS - 1 )%q( 1 ) )
755
#####
tau = ray2D( iS - 1 )%tau + W * ( ray2D( iS )%tau - ray2D( iS - 1 )%tau )
756
-
757
-
! BUG: following is incorrect because ray doesn't always use a step of deltas
758
-
! LP: The do while ignores extremely small steps, but those small steps
759
-
! still increment iS, so the later ray segments still treat it as if
760
-
! all steps leading up to them were of size deltas.
761
#####
SINT = ( iS - 1 ) * Beam%deltas + W * Beam%deltas
762
-
763
#####
CALL IncPhaseIfCaustic( .FALSE. )
764
-
765
-
! WRITE( PRTFile, * ) 'is ir', is-2, ir-1
766
-
767
#####
RcvrDepths: DO iz = 1, NRz_per_range
768
#####
deltaz = Pos%Rz( iz ) - x( 2 ) ! ray to rcvr distance
769
-
! LP: Reinstated this condition for eigenrays and arrivals, as
770
-
! without it every ray would be an eigenray / arrival.
771
#####
Adeltaz = ABS( deltaz )
772
-
IF ( Adeltaz < RadiusMax .OR. Beam%RunType( 1 : 1 ) == 'C' &
773
#####
.OR. Beam%RunType( 1 : 1 ) == 'I' .OR. Beam%RunType( 1 : 1 ) == 'S' ) THEN
774
-
! LP: Changed to use ApplyContribution in order to support
775
-
! incoherent, semi-coherent, and arrivals.
776
#####
SrcDeclAngle = RadDeg * alpha ! take-off angle in degrees
777
-
CPA = ABS( deltaz * ( rB - rA ) ) / SQRT( ( rB - rA )**2 + &
778
#####
& ( ray2D( iS )%x( 2 ) - ray2D( iS - 1 )%x( 2 ) )**2 )
779
#####
DS = SQRT( deltaz **2 - CPA **2 )
780
#####
SX1 = SINT + DS
781
#####
thet = ATAN( CPA / SX1 )
782
#####
delay = tau + rayt( 2 ) * deltaz
783
#####
const = Ratio1 * CN * ray2D( iS )%Amp / SQRT( SX1 )
784
#####
W = EXP( -A * thet ** 2 )
785
#####
Amp = const * W
786
#####
phaseInt = ray2D( iS )%Phase + phase
787
#####
CALL ApplyContribution( U( iz, ir ) )
788
-
789
-
END IF
790
-
END DO RcvrDepths
791
-
792
#####
qOld = q
793
#####
ir = ir + 1
794
#####
IF ( ir > Pos%NRr ) RETURN
795
-
END DO RcvrRanges
796
-
797
#####
rA = rB
798
-
END DO Stepping
799
-
800
-
END SUBROUTINE InfluenceSGB
801
-
802
-
! **********************************************************************!
803
-
804
#####
SUBROUTINE BranchCut( q1C, q2C, BeamType, KMAH )
805
-
!! Checks for a branch cut crossing and updates KMAH accordingly
806
-
807
-
COMPLEX (KIND=8), INTENT( IN ) :: q1C, q2C
808
-
CHARACTER (LEN=4), INTENT( IN ) :: BeamType
809
-
INTEGER, INTENT( INOUT ) :: KMAH
810
-
REAL (KIND=8) :: q1, q2
811
-
812
#####
SELECT CASE ( BeamType( 2 : 2 ) )
813
-
CASE ( 'W' ) ! WKBeams
814
#####
q1 = REAL( q1C )
815
#####
q2 = REAL( q2C )
816
#####
IF ( ( q1 < 0.0 .AND. q2 >= 0.0 ) .OR. &
817
#####
( q1 > 0.0 .AND. q2 <= 0.0 ) ) KMAH = -KMAH
818
-
CASE DEFAULT
819
#####
IF ( REAL( q2C ) < 0.0 ) THEN
820
#####
q1 = AIMAG( q1C )
821
#####
q2 = AIMAG( q2C )
822
#####
IF ( ( q1 < 0.0 .AND. q2 >= 0.0 ) .OR. &
823
#####
( q1 > 0.0 .AND. q2 <= 0.0 ) ) KMAH = -KMAH
824
-
END IF
825
-
END SELECT
826
-
827
#####
END SUBROUTINE BranchCut
828
-
829
-
! **********************************************************************!
830
-
831
16*
SUBROUTINE ScalePressure( Dalpha, c, r, U, NRz, Nr, RunType, freq )
832
-
!! Scale the pressure field
833
-
834
-
REAL, PARAMETER :: pi = 3.14159265
835
-
INTEGER, INTENT( IN ) :: NRz, Nr
836
-
REAL, INTENT( IN ) :: r( Nr ) ! ranges
837
-
REAL (KIND=8), INTENT( IN ) :: Dalpha, freq, c ! angular spacing between rays, source frequency, nominal sound speed
838
-
COMPLEX, INTENT( INOUT ) :: U( NRz, Nr ) ! Pressure field
839
-
CHARACTER (LEN=5), INTENT( IN ) :: RunType
840
-
REAL (KIND=8) :: const, factor
841
-
842
-
! Compute scale factor for field
843
#####
SELECT CASE ( RunType( 2 : 2 ) )
844
-
CASE ( 'C' ) ! Cerveny Gaussian beams in Cartesian coordinates
845
#####
const = -Dalpha * SQRT( freq ) / c
846
-
CASE ( 'R' ) ! Cerveny Gaussian beams in Ray-centered coordinates
847
#####
const = -Dalpha * SQRT( freq ) / c
848
-
CASE DEFAULT
849
16
const = -1.0
850
-
END SELECT
851
-
852
16*
IF ( RunType( 1 : 1 ) /= 'C' ) U = SQRT( REAL( U ) ) ! For incoherent run, convert intensity to pressure
853
-
854
-
! scale and/or incorporate cylindrical spreading
855
27527*
Ranges: DO ir = 1, Nr
856
27511
IF ( RunType( 4 : 4 ) == 'X' ) THEN ! line source
857
501
factor = -4.0 * SQRT( pi ) * const
858
-
ELSE ! point source
859
27010*
IF ( r ( ir ) == 0 ) THEN
860
9
factor = 0.0D0 ! avoid /0 at origin, return pressure = 0
861
-
ELSE
862
27001*
factor = const / SQRT( ABS( r( ir ) ) )
863
-
END IF
864
-
END IF
865
2378438*
U( :, ir ) = SNGL( factor ) * U( :, ir )
866
-
END DO Ranges
867
-
868
16
END SUBROUTINE ScalePressure
869
-
870
-
! **********************************************************************!
871
-
872
#####
REAL (KIND=8 ) FUNCTION Hermite( x, x1, x2 )
873
-
874
-
! Calculates a smoothing function based on the h0 hermite cubic
875
-
! x is the point where the function is to be evaluated
876
-
! returns:
877
-
! [ 0, x1 ] = 1
878
-
! [ x1, x2 ] = cubic taper from 1 to 0
879
-
! [ x2, inf ] = 0
880
-
881
-
REAL (KIND=8 ), INTENT( IN ) :: x, x1, x2
882
-
REAL (KIND=8 ) :: Ax, u
883
-
884
#####
Ax = ABS( x )
885
-
886
#####
IF ( Ax <= x1 ) THEN
887
#####
Hermite = 1.0d0
888
#####
ELSE IF ( Ax >= x2 ) THEN
889
#####
Hermite = 0.0d0
890
-
ELSE
891
#####
u = ( Ax - x1 ) / ( x2 - x1 )
892
#####
Hermite = ( 1.0d0 + 2.0d0 * u ) * ( 1.0d0 - u ) ** 2
893
-
END IF
894
-
895
-
!hermit = hermit / ( 0.5 * ( x1 + x2 ) )
896
-
897
#####
END FUNCTION Hermite
898
-
899
-
! **********************************************************************!
900
-
901
82394521
SUBROUTINE FinalPhase( isGaussian )
902
-
LOGICAL, INTENT( IN ) :: isGaussian
903
-
INTEGER :: phaseStepNum
904
-
905
-
!! phase shifts at caustics
906
-
907
-
! this should be precomputed [LP: While IncPhaseIfCaustic can be
908
-
! precomputed, FinalPhase cannot, as it is dependent on the interpolated `q`
909
-
! value which is not known until the main run.]
910
-
! LP: All 2D functions discard the ray point phase if the condition is met,
911
-
! probably BUG
912
-
! LP: 2D Gaussian Cartesian reads the phase from the current point, all
913
-
! others (including 3D) read the phase from the previous point, probably BUG
914
82394521
IF ( isGaussian ) THEN
915
76021312
phaseStepNum = iS
916
-
ELSE
917
6373209
phaseStepNum = iS - 1
918
-
END IF
919
-
920
82394521*
phaseInt = ray2D( phaseStepNum )%Phase + phase
921
82394521
IF ( IsAtCaustic( .TRUE. ) ) THEN
922
89091
phaseInt = phase + pi / 2.
923
-
END IF
924
-
925
82394521
END SUBROUTINE FinalPhase
926
-
927
-
! **********************************************************************!
928
-
929
33946447
SUBROUTINE IncPhaseIfCaustic( qleq0 )
930
-
931
-
!! phase shifts at caustics
932
-
933
-
LOGICAL, INTENT( IN ) :: qleq0
934
-
935
33946447
IF ( IsAtCaustic( qleq0 ) ) THEN
936
4898
phase = phase + pi / 2.
937
-
END IF
938
-
939
33946447
END SUBROUTINE IncPhaseIfCaustic
940
-
941
-
! **********************************************************************!
942
-
943
116340968
LOGICAL FUNCTION IsAtCaustic( qleq0 )
944
-
945
-
! LP: There are two versions of the phase shift condition used in the
946
-
! BELLHOP code, with the equalities in opposite positions. qleq0 false is
947
-
! only used in SGB.
948
-
949
-
LOGICAL, INTENT( IN ) :: qleq0
950
-
951
116340968
IF ( qleq0 ) THEN
952
116340968
IsAtCaustic = q <= 0.0d0 .AND. qOld > 0.0d0 .OR. q >= 0.0d0 .AND. qOld < 0.0d0
953
-
ELSE
954
#####
IsAtCaustic = q < 0.0d0 .AND. qOld >= 0.0d0 .OR. q > 0.0d0 .AND. qOld <= 0.0d0
955
-
END IF
956
-
957
116340968
END FUNCTION IsAtCaustic
958
-
959
-
END MODULE Influence

Legend:

Execution Count: Number of times each line was executed. ##### indicates never executed, - indicates non-executable.