1
-
!! Ray stepping module with adaptive step size control
3
-
!! Ray propagation with adaptive step size control and boundary interaction handling
11
-
REAL (KIND=8), PARAMETER, PRIVATE :: INFINITESIMAL_STEP_SIZE = 1.0d-6
15
13986421
SUBROUTINE Step2D( ray0, ray2, topRefl, botRefl, Topx, Topn, Botx, Botn )
17
-
! Does a single step along the ray
18
-
! x denotes the ray coordinate, (r,z)
19
-
! t denotes the scaled tangent to the ray (previously (rho, zeta))
20
-
! c * t would be the unit tangent
23
-
TYPE( ray2DPt ), INTENT( INOUT ) :: ray0, ray2
24
-
TYPE( ray2DPt ) :: ray1
25
-
REAL (KIND=8 ), INTENT( IN ) :: Topx( 2 ), Topn( 2 ), Botx( 2 ), Botn( 2 )
26
-
LOGICAL, INTENT( OUT ) :: topRefl, botRefl
27
-
INTEGER :: iSegz0, iSegr0
28
-
REAL (KIND=8 ) :: gradc0( 2 ), gradc1( 2 ), gradc2( 2 ), &
29
-
c0, cimag0, crr0, crz0, czz0, csq0, cnn0_csq0, &
30
-
c1, cimag1, crr1, crz1, czz1, csq1, cnn1_csq1, &
31
-
c2, cimag2, crr2, crz2, czz2, urayt0( 2 ), urayt1( 2 ), &
32
-
h, halfh, ray2n( 2 ), RM, RN, gradcjump( 2 ), cnjump, csjump, w0, w1, rho
33
-
REAL (KIND=8 ) :: urayt2( 2 ), unitdt( 2 ), unitdp( 2 ), unitdq( 2 )
34
-
COMPLEX (KIND=8 ) :: unitdtau
36
-
IF ( STEP_DEBUGGING ) THEN
38
-
WRITE( PRTFile, * ) 'ray0 x t p q tau amp', ray0%x, ray0%t, ray0%p, ray0%q, ray0%tau, ray0%Amp
39
-
WRITE( PRTFile, * ) 'iSegr iSegz', iSegr, iSegz
42
-
! IF ( ray0%x( 1 ) > 10.0 ) THEN
46
-
! The numerical integrator used here is a version of the polygon (a.k.a. midpoint, leapfrog, or Box method), and similar
47
-
! to the Heun (second order Runge-Kutta method).
48
-
! However, it's modified to allow for a dynamic step change, while preserving the second-order accuracy).
50
-
! *** Phase 1 (an Euler step)
52
13986421
CALL EvaluateSSP( ray0%x, ray0%t, c0, cimag0, gradc0, crr0, crz0, czz0, rho, freq, 'TAB' )
54
13986421
csq0 = c0 * c0
55
13986421
cnn0_csq0 = crr0 * ray0%t( 2 )**2 - 2.0 * crz0 * ray0%t( 1 ) * ray0%t( 2 ) + czz0 * ray0%t( 1 )**2
56
13986421
iSegz0 = iSegz ! make note of current layer
57
13986421
iSegr0 = iSegr
59
13986421
h = Beam%deltas ! initially set the step h, to the basic one, deltas
60
41959263
urayt0 = c0 * ray0%t ! unit tangent
62
13986421
CALL ReduceStep2D( ray0%x, urayt0, iSegz0, iSegr0, Topx, Topn, Botx, Botn, h ) ! reduce h to land on boundary
63
-
! WRITE( PRTFile, * ) 'out h, urayt0', h, urayt0
64
13986421
halfh = 0.5 * h ! first step of the modified polygon method is a half step
66
41959263
ray1%x = ray0%x + halfh * urayt0
67
41959263
ray1%t = ray0%t - halfh * gradc0 / csq0
68
41959263
ray1%p = ray0%p - halfh * cnn0_csq0 * ray0%q
69
41959263
ray1%q = ray0%q + halfh * c0 * ray0%p
71
-
! WRITE( PRTFile, * ) 'ray1 x t p q', ray1%x, ray1%t, ray1%p, ray1%q
75
13986421
CALL EvaluateSSP( ray1%x, ray1%t, c1, cimag1, gradc1, crr1, crz1, czz1, rho, freq, 'TAB' )
76
13986421
csq1 = c1 * c1
77
13986421
cnn1_csq1 = crr1 * ray1%t( 2 )**2 - 2.0 * crz1 * ray1%t( 1 ) * ray1%t( 2 ) + czz1 * ray1%t( 1 )**2
79
-
! The Munk test case with a horizontally launched ray caused problems.
80
-
! The ray vertexes on an interface and can ping-pong around that interface.
81
-
! Have to be careful in that case about big changes to the stepsize (that invalidate the leap-frog scheme) in phase II.
82
-
! A modified Heun or Box method could also work.
84
41959263
urayt1 = c1 * ray1%t ! unit tangent
86
13986421
CALL ReduceStep2D( ray0%x, urayt1, iSegz0, iSegr0, Topx, Topn, Botx, Botn, h ) ! reduce h to land on boundary
88
-
! WRITE( PRTFile, * ) 'urayt1', urayt1, 'out h 2', h
90
-
! use blend of f' based on proportion of a full step used.
91
13986421
w1 = h / ( 2.0d0 * halfh )
92
13986421
w0 = 1.0d0 - w1
93
-
! WRITE( PRTFile, * ) 'w1 w0', w1, w0
94
41959263
urayt2 = w0 * urayt0 + w1 * urayt1
95
41959263
unitdt = -w0 * gradc0 / csq0 - w1 * gradc1 / csq1
96
41959263
unitdp = -w0 * cnn0_csq0 * ray0%q - w1 * cnn1_csq1 * ray1%q
97
41959263
unitdq = w0 * c0 * ray0%p + w1 * c1 * ray1%p
98
13986421
unitdtau = w0 / CMPLX( c0, cimag0, KIND=8 ) + w1 / CMPLX( c1, cimag1, KIND=8 )
100
-
! Take the blended ray tangent (urayt2) and find the minimum step size (h)
101
-
! to put this on a boundary, and ensure that the resulting position
102
-
! (ray2.x) gets put precisely on the boundary.
103
-
CALL StepToBdry2D(ray0%x, ray2%x, urayt2, h, topRefl, botRefl, &
104
13986421
iSegz0, iSegr0, Topx, Topn, Botx, Botn)
105
41959263
ray2%t = ray0%t + h * unitdt
106
41959263
ray2%p = ray0%p + h * unitdp
107
41959263
ray2%q = ray0%q + h * unitdq
108
13986421
ray2%tau = ray0%tau + h * unitdtau
110
-
! WRITE( PRTFile, * ) 'ray2 x t p q tau', ray2%x, ray2%t, ray2%p, ray2%q, ray2%tau
112
13986421
ray2%Amp = ray0%Amp
113
13986421
ray2%Phase = ray0%Phase
114
13986421
ray2%NumTopBnc = ray0%NumTopBnc
115
13986421
ray2%NumBotBnc = ray0%NumBotBnc
117
-
! If we crossed an interface, apply jump condition
119
13986421
CALL EvaluateSSP( ray2%x, ray2%t, c2, cimag2, gradc2, crr2, crz2, czz2, rho, freq, 'TAB' )
122
13986421
IF ( iSegz /= iSegz0 .OR. iSegr /= iSegr0 ) THEN
123
2650527
gradcjump = gradc2 - gradc0
124
2650527
ray2n = [ -ray2%t( 2 ), ray2%t( 1 ) ] ! ray normal
126
2650527
cnjump = DOT_PRODUCT( gradcjump, ray2n )
127
2650527
csjump = DOT_PRODUCT( gradcjump, ray2%t )
129
883509
IF ( iSegz /= iSegz0 ) THEN ! crossing in depth
130
868780
RM = +ray2%t( 1 ) / ray2%t( 2 ) ! this is tan( alpha ) where alpha is the angle of incidence
131
-
ELSE ! crossing in range
132
14729
RM = -ray2%t( 2 ) / ray2%t( 1 ) ! this is tan( alpha ) where alpha is the angle of incidence
135
883509
RN = RM * ( 2 * cnjump - RM * csjump ) / c2
136
2650527
ray2%p = ray2%p - ray2%q * RN
140
13986421
END SUBROUTINE Step2D
142
-
! **********************************************************************!
144
27972842
SUBROUTINE ReduceStep2D( x0, urayt, iSegz0, iSegr0, Topx, Topn, Botx, Botn, h )
146
-
! calculate a reduced step size, h, that lands on any points where the environment changes
149
-
INTEGER, INTENT( IN ) :: iSegz0, iSegr0 ! SSP layer the ray is in
150
-
REAL (KIND=8), INTENT( IN ) :: x0( 2 ), urayt( 2 ) ! ray coordinate and tangent
151
-
REAL (KIND=8), INTENT( IN ) :: Topx( 2 ), Topn( 2 ), Botx( 2 ), Botn( 2 ) ! Top, bottom coordinate and normal
152
-
REAL (KIND=8), INTENT( INOUT ) :: h ! reduced step size
153
-
REAL (KIND=8) :: hInt, hBoxr, hBoxz, hTop, hBot, hSeg ! step sizes
154
-
REAL (KIND=8) :: x( 2 ), d( 2 ), d0( 2 ), rSeg( 2 )
156
-
! Detect interface or boundary crossing and reduce step, if necessary, to land on that crossing.
157
-
! Keep in mind possibility that user put source right on an interface
158
-
! and that multiple events can occur (crossing interface, top, and bottom in a single step).
160
83918526
x = x0 + h * urayt ! make a trial step
162
-
! interface crossing in depth
163
27972842
hInt = huge( hInt )
164
27972842
IF ( ABS( urayt( 2 ) ) > EPSILON( hInt ) ) THEN
165
27972842*
IF ( SSP%z( iSegz0 ) > x( 2 ) ) THEN
166
975361*
hInt = ( SSP%z( iSegz0 ) - x0( 2 ) ) / urayt( 2 )
167
26997481*
ELSE IF ( SSP%z( iSegz0 + 1 ) < x( 2 ) ) THEN
168
962189*
hInt = ( SSP%z( iSegz0 + 1 ) - x0( 2 ) ) / urayt( 2 )
172
-
! ray mask using a box centered at ( 0, 0 )
173
27972842
hBoxr = huge( hBoxr )
174
27972842
IF ( ABS( x( 1 ) ) > Beam%Box%r ) THEN
175
763736
hBoxr = ( Beam%Box%r - ABS( x0( 1 ) ) ) / ABS( urayt( 1 ) )
178
27972842
hBoxz = huge( hBoxz )
179
27972842
IF ( ABS( x( 2 ) ) > Beam%Box%z ) THEN
180
250632
hBoxz = ( Beam%Box%z - ABS( x0( 2 ) ) ) / ABS( urayt( 2 ) )
184
27972842
hTop = huge( hTop )
185
83918526
d = x - Topx ! vector from top to ray
186
-
! LP: Changed from EPSILON( hTop ) to 0 in conjunction with StepToBdry2D change here.
187
83918526
IF ( DOT_PRODUCT( Topn, d ) >= 0.0D0 ) THEN
188
2117493
d0 = x0 - Topx ! vector from top node to ray origin
189
3529155
hTop = -DOT_PRODUCT( d0, Topn ) / DOT_PRODUCT( urayt, Topn )
193
27972842
hBot = huge( hBot )
194
83918526
d = x - Botx ! vector from bottom to ray
195
-
! LP: Changed from EPSILON( hBot ) to 0 in conjunction with StepToBdry2D change here.
196
83918526
IF ( DOT_PRODUCT( Botn, d ) >= 0.0D0 ) THEN
197
1905915
d0 = x0 - Botx ! vector from bottom node to ray origin
198
3176525
hBot = -DOT_PRODUCT( d0, Botn ) / DOT_PRODUCT( urayt, Botn )
201
-
! top or bottom segment crossing in range
202
27972842
rSeg( 1 ) = MAX( rTopSeg( 1 ), rBotSeg( 1 ) )
203
27972842
rSeg( 2 ) = MIN( rTopSeg( 2 ), rBotSeg( 2 ) )
205
27972842
IF ( SSP%Type == 'Q' ) THEN
206
1525908*
rSeg( 1 ) = MAX( rSeg( 1 ), SSP%Seg%r( iSegr0 ) )
207
1525908*
rSeg( 2 ) = MIN( rSeg( 2 ), SSP%Seg%r( iSegr0 + 1 ) )
210
27972842
hSeg = huge( hSeg )
211
27972842
IF ( ABS( urayt( 1 ) ) > EPSILON( hSeg ) ) THEN
212
27972842
IF ( x( 1 ) < rSeg( 1 ) ) THEN
213
10443
hSeg = -( x0( 1 ) - rSeg( 1 ) ) / urayt( 1 )
214
27962399
ELSE IF ( x( 1 ) > rSeg( 2 ) ) THEN
215
1529313
hSeg = -( x0( 1 ) - rSeg( 2 ) ) / urayt( 1 )
219
27972842
h = MIN( h, hInt, hBoxr, hBoxz, hTop, hBot, hSeg ) ! take limit set by shortest distance to a crossing
220
27972842
IF ( h < INFINITESIMAL_STEP_SIZE * Beam%deltas ) THEN ! is it taking an infinitesimal step?
221
331
h = INFINITESIMAL_STEP_SIZE * Beam%deltas ! make sure we make some motion
222
331
iSmallStepCtr = iSmallStepCtr + 1 ! keep a count of the number of sequential small steps
223
-
IF ( STEP_DEBUGGING ) THEN
224
-
WRITE( PRTFile, * ) 'Small step forced to', h
227
27972511
iSmallStepCtr = 0 ! didn't do a small step so reset the counter
230
27972842
END SUBROUTINE ReduceStep2D
232
-
! **********************************************************************!
234
13986421
SUBROUTINE StepToBdry2D( x0, x2, urayt, h, topRefl, botRefl, &
235
-
iSegz0, iSegr0, Topx, Topn, Botx, Botn )
237
-
REAL (KIND=8), INTENT( IN ) :: x0( 2 ), urayt( 2 )
238
-
REAL (KIND=8), INTENT( INOUT ) :: x2( 2 ), h
239
-
INTEGER, INTENT( IN ) :: iSegz0, iSegr0 ! SSP layer the ray is in
240
-
REAL (KIND=8), INTENT( IN ) :: Topx( 2 ), Topn( 2 ), Botx( 2 ), Botn( 2 ) ! Top, bottom coordinate and normal
241
-
LOGICAL, INTENT( OUT ) :: topRefl, botRefl
242
-
REAL (KIND=8) :: d( 2 ), d0( 2 ), rSeg( 2 )
244
-
! Original step due to maximum step size
245
13986421
h = Beam%deltas
246
41959263
x2 = x0 + h * urayt
248
-
! interface crossing in depth
249
13986421
IF ( ABS( urayt( 2 ) ) > EPSILON( h ) ) THEN
250
13986421*
IF ( SSP%z( iSegz0 ) > x2( 2 ) ) THEN
251
714408*
h = ( SSP%z( iSegz0 ) - x0( 2 ) ) / urayt( 2 )
252
714408
x2( 1 ) = x0( 1 ) + h * urayt( 1 )
253
714408*
x2( 2 ) = SSP%z( iSegz0 )
254
-
IF ( STEP_DEBUGGING ) THEN
255
-
WRITE( PRTFile, * ) 'StepToBdry2D lower depth h to', h, x2
257
13272013*
ELSE IF ( SSP%z( iSegz0 + 1 ) < x2( 2 ) ) THEN
258
802206*
h = ( SSP%z( iSegz0 + 1 ) - x0( 2 ) ) / urayt( 2 )
259
802206
x2( 1 ) = x0( 1 ) + h * urayt( 1 )
260
802206*
x2( 2 ) = SSP%z( iSegz0 + 1 )
261
-
IF ( STEP_DEBUGGING ) THEN
262
-
WRITE( PRTFile, * ) 'StepToBdry2D upper depth h to', h, x2
267
-
! ray mask using a box centered at ( 0, 0 )
268
13986421
IF ( ABS( x2( 1 ) ) > Beam%Box%r ) THEN
269
767463
h = ( Beam%Box%r - ABS( x0( 1 ) ) ) / ABS( urayt( 1 ) )
270
767463
x2( 2 ) = x0( 2 ) + h * urayt( 2 )
271
767463
x2( 1 ) = SIGN( Beam%Box%r, x0( 1 ) )
272
-
IF ( STEP_DEBUGGING ) THEN
273
-
WRITE( PRTFile, * ) 'StepToBdry2D beam box crossing R h to', h, x2
276
13986421
IF ( ABS( x2( 2 ) ) > Beam%Box%z ) THEN
277
726
h = ( Beam%Box%z - ABS( x0( 2 ) ) ) / ABS( urayt( 2 ) )
278
726
x2( 1 ) = x0( 1 ) + h * urayt( 1 )
279
726
x2( 2 ) = SIGN( Beam%Box%z, x0( 2 ) )
280
-
IF ( STEP_DEBUGGING ) THEN
281
-
WRITE( PRTFile, * ) 'StepToBdry2D beam box crossing Z h to', h, x2
286
41959263
d = x2 - Topx ! vector from top to ray
287
-
! LP: Originally, this value had to be > a small positive number, meaning the
288
-
! new step really had to be outside the boundary, not just to the boundary.
289
-
! Also, this is not missing a normalization factor, Topn is normalized so
290
-
! this is actually the distance above the top in meters.
291
41959263
IF ( DOT_PRODUCT( Topn, d ) > -INFINITESIMAL_STEP_SIZE ) THEN
292
988533
d0 = x0 - Topx ! vector from top node to ray origin
293
1647555
h = -DOT_PRODUCT( d0, Topn ) / DOT_PRODUCT( urayt, Topn )
294
988533
x2 = x0 + h * urayt
295
-
! Snap to exact top depth value if it's flat
296
329511
IF ( ABS( Topn( 1 ) ) < EPSILON( Topn( 1 ) ) ) THEN
297
329511
x2( 2 ) = Topx( 2 )
299
-
IF ( STEP_DEBUGGING ) THEN
300
-
WRITE( PRTFile, * ) 'StepToBdry2D top crossing h to', h, x2
302
329511
topRefl = .TRUE.
304
13656910
topRefl = .FALSE.
308
41959263
d = x2 - Botx ! vector from bottom to ray
309
-
! See comment above for top case.
310
41959263
IF ( DOT_PRODUCT( Botn, d ) > -INFINITESIMAL_STEP_SIZE ) THEN
311
1042665
d0 = x0 - Botx ! vector from bottom node to ray origin
312
1737775
h = -DOT_PRODUCT( d0, Botn ) / DOT_PRODUCT( urayt, Botn )
313
1042665
x2 = x0 + h * urayt
314
-
! Snap to exact bottom depth value if it's flat
315
347555
IF ( ABS( Botn( 1 ) ) < EPSILON( Botn( 1 ) ) ) THEN
316
307350
x2( 2 ) = Botx( 2 )
318
-
IF ( STEP_DEBUGGING ) THEN
319
-
WRITE( PRTFile, * ) 'StepToBdry2D bottom crossing h to', h, x2
321
347555
botRefl = .TRUE.
322
-
! Should not ever be able to cross both, but in case it does, make sure
323
-
! only the crossing we exactly landed on is active
324
347555
topRefl = .FALSE.
326
13638866
botRefl = .FALSE.
329
-
! top or bottom segment crossing in range
330
13986421
rSeg( 1 ) = MAX( rTopSeg( 1 ), rBotSeg( 1 ) ) ! LP: lower range bound (not an x value)
331
13986421
rSeg( 2 ) = MIN( rTopSeg( 2 ), rBotSeg( 2 ) ) ! LP: upper range bound (not a y value)
333
13986421
IF ( SSP%Type == 'Q' ) THEN
334
762954*
rSeg( 1 ) = MAX( rSeg( 1 ), SSP%Seg%r( iSegr0 ) )
335
762954*
rSeg( 2 ) = MIN( rSeg( 2 ), SSP%Seg%r( iSegr0 + 1 ) )
338
13986421
IF ( ABS( urayt( 1 ) ) > EPSILON( h ) ) THEN
339
13986421
IF ( x2( 1 ) < rSeg( 1 ) ) THEN
340
753261
h = -( x0( 1 ) - rSeg( 1 ) ) / urayt( 1 )
341
753261
x2( 1 ) = rSeg( 1 )
342
753261
x2( 2 ) = x0( 2 ) + h * urayt( 2 )
343
-
IF ( STEP_DEBUGGING ) THEN
344
-
WRITE( PRTFile, * ) 'StepToBdry2D lower range h to', h, x2
346
753261
topRefl = .FALSE.
347
753261
botRefl = .FALSE.
348
13233160
ELSE IF ( x2( 1 ) > rSeg( 2 ) ) THEN
349
18654
h = -( x0( 1 ) - rSeg( 2 ) ) / urayt( 1 )
350
18654
x2( 1 ) = rSeg( 2 )
351
18654
x2( 2 ) = x0( 2 ) + h * urayt( 2 )
352
-
IF ( STEP_DEBUGGING ) THEN
353
-
WRITE( PRTFile, * ) 'StepToBdry2D upper range h to', h, x2
355
18654
topRefl = .FALSE.
356
18654
botRefl = .FALSE.
360
13986421
IF ( h < INFINITESIMAL_STEP_SIZE * Beam%deltas ) THEN ! is it taking an infinitesimal step?
361
747216
h = INFINITESIMAL_STEP_SIZE * Beam%deltas ! make sure we make some motion
362
2241648
x2 = x0 + h * urayt
363
-
IF ( STEP_DEBUGGING ) THEN
364
-
WRITE( PRTFile, * ) 'StepToBdry2D small step forced h to', h, x2
366
-
! Recheck reflection conditions
367
2241648
d = x2 - Topx ! vector from top to ray
368
2241648
IF ( DOT_PRODUCT( Topn, d ) > EPSILON( h ) ) THEN
371
747215
topRefl = .FALSE.
373
2241648
d = x2 - Botx ! vector from bottom to ray
374
2241648
IF ( DOT_PRODUCT( Botn, d ) > EPSILON( h ) ) THEN
375
#####
botRefl = .TRUE.
376
#####
topRefl = .FALSE.
378
747216
botRefl = .FALSE.
381
-
IF ( STEP_DEBUGGING ) THEN
382
-
WRITE( PRTFile, * ) 'StepToBdry2D normal h to', h, x2
386
13986421
END SUBROUTINE StepToBdry2D