Coverage Report: influence.f90

Generated from GCOV analysis of Fortran source code

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

Legend:

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