Coverage Report: influence.f90

Generated from GCOV analysis of Fortran source code

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

Legend:

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