Coverage Report: ReflectMod.f90

Generated from GCOV analysis of Fortran source code

0.0%
Lines Executed
80 total lines
0.0%
Branches Executed
0 total branches
0.0%
Calls Executed
0 total calls
0
-
Source:ReflectMod.f90
0
-
Graph:ReflectMod.gcno
0
-
Data:ReflectMod.gcda
0
-
Runs:29
1
-
!! Reflection calculations for acoustic boundaries
2
-
MODULE ReflectMod
3
-
!! Ray reflection computations at acoustic boundaries with loss and phase calculations
4
-
5
-
USE bellhopMod
6
-
7
-
IMPLICIT NONE
8
-
PUBLIC
9
-
10
-
CONTAINS
11
-
12
-
!! Computes reflection of ray at acoustic boundary
13
#####
SUBROUTINE Reflect2D( is, HS, BotTop, nBdry3d, z_xx, z_xy, z_yy, kappa_xx, kappa_xy, kappa_yy, RefC, Npts, tradial )
14
-
15
-
!USE norms
16
-
USE RefCoef
17
-
USE sspMod
18
-
USE Cone ! if using analytic formulas for the cone (conical seamount)
19
-
20
-
INTEGER, INTENT( INOUT ) :: is
21
-
TYPE( HSInfo ), INTENT( IN ) :: HS ! half-space properties
22
-
CHARACTER (LEN=3), INTENT( IN ) :: BotTop ! Flag indicating bottom or top reflection
23
-
REAL (KIND=8), INTENT( IN ) :: nBdry3d( 3 ) ! Normal to the boundary
24
-
REAL (KIND=8), INTENT( IN ) :: z_xx, z_xy, z_yy
25
-
REAL (KIND=8), INTENT( IN ) :: kappa_xx, kappa_xy, kappa_yy
26
-
TYPE(ReflectionCoef), INTENT( IN ) :: RefC( NPts ) ! reflection coefficient
27
-
INTEGER, INTENT( IN ) :: Npts
28
-
REAL (KIND=8), INTENT( IN ) :: tradial( 2 )
29
-
INTEGER :: is1
30
-
REAL (KIND=8) :: c, cimag, gradc( 2 ), crr, crz, czz, rho ! derivatives of sound speed
31
-
REAL (KIND=8) :: RM, RN, Tg, Th, rayt( 2 ), rayn( 2 ), rayt_tilde( 2 ), rayn_tilde( 2 ), cnjump, csjump ! for curvature change
32
-
REAL (KIND=8) :: ck, co, si, cco, ssi, pdelta, rddelta, sddelta ! for beam shift
33
-
COMPLEX (KIND=8) :: gamma1, gamma2, gamma1Sq, gamma2Sq, GK, Refl, ch, a, b, d, sb, delta, ddelta
34
-
REAL (KIND=8) :: kappa ! Boundary curvature
35
-
REAL (KIND=8) :: tBdry( 2 ), nBdry( 2 ) ! tangent, normal to boundary
36
-
TYPE(ReflectionCoef) :: RInt
37
-
REAL (KIND=8) :: t_rot( 2 ), n_rot( 2 ), RotMat( 2, 2 ), kappaMat( 2, 2 ), DMat( 2, 2 ) ! for cone reflection
38
-
39
#####
is = is + 1
40
#####
is1 = is + 1
41
-
42
-
! following should perhaps be pre-calculated in Bdry3DMod
43
#####
nBdry( 1 ) = DOT_PRODUCT( nBdry3d( 1 : 2 ), tradial )
44
#####
nBdry( 2 ) = nBdry3d( 3 )
45
-
46
-
!CALL ConeFormulas( z_xx, z_xy, z_yy, nBdry, xs_3D, tradial, ray2D( is )%x, BotTop ) ! special case of a conical seamount
47
-
48
-
!!!! use kappa_xx or z_xx?
49
-
! kappa = ( z_xx * tradial( 1 ) ** 2 + 2 * z_xy * tradial( 1 ) * tradial( 2 ) + z_yy * tradial( 2 ) ** 2 )
50
#####
kappa = ( kappa_xx * tradial( 1 ) ** 2 + 2 * kappa_xy * tradial( 1 ) * tradial( 2 ) + kappa_yy * tradial( 2 ) ** 2 )
51
-
52
#####
IF ( BotTop == 'TOP' ) kappa = -kappa
53
-
54
#####
nBdry = nBdry / NORM2( nBdry )
55
-
56
#####
Th = DOT_PRODUCT( ray2D( is )%t, nBdry ) ! component of ray tangent normal to boundary
57
-
58
#####
tBdry = ray2D( is )%t - Th * nBdry ! component of ray tangent along the boundary
59
#####
tBdry = tBdry / NORM2( tBdry )
60
-
! could also calculate tbdry as +/- of [ nbdry( 2), -nbdry( 1 ) ], but need sign
61
-
62
#####
Tg = DOT_PRODUCT( ray2D( is )%t, tBdry ) ! component of ray tangent along boundary
63
-
64
#####
ray2D( is1 )%NumTopBnc = ray2D( is )%NumTopBnc
65
#####
ray2D( is1 )%NumBotBnc = ray2D( is )%NumBotBnc
66
#####
ray2D( is1 )%x = ray2D( is )%x
67
#####
ray2D( is1 )%t = ray2D( is )%t - 2.0 * Th * nBdry ! changing the ray direction
68
-
69
-
! Calculate the change in curvature
70
-
! Based on formulas given by Muller, Geoph. J. R.A.S., 79 (1984).
71
-
72
#####
CALL EvaluateSSP2D( ray2D( is1 )%x, ray2D( is1 )%t, c, cimag, gradc, crr, crz, czz, rho, xs_3D, tradial, freq )
73
-
74
-
! incident and reflected (tilde) unit ray tangent and normal
75
#####
rayt = c * ray2D( is )%t ! unit tangent to ray
76
#####
rayt_tilde = c * ray2D( is1 )%t ! unit tangent to ray
77
#####
rayn = +[ -rayt( 2 ), rayt( 1 ) ] ! unit normal to ray
78
#####
rayn_tilde = -[ -rayt_tilde( 2 ), rayt_tilde( 1 ) ] ! unit normal to ray
79
-
80
#####
RN = 2 * kappa / c ** 2 / Th ! boundary curvature correction
81
-
82
-
! get the jumps (this could be simplified, e.g. jump in rayt is roughly 2 * Th * nbdry
83
#####
cnjump = -DOT_PRODUCT( gradc, rayn_tilde - rayn )
84
#####
csjump = -DOT_PRODUCT( gradc, rayt_tilde - rayt )
85
-
86
#####
IF ( BotTop == 'TOP' ) THEN
87
#####
cnjump = -cnjump ! flip sign for top reflection
88
#####
RN = -RN
89
-
END IF
90
-
91
#####
RM = Tg / Th ! this is tan( alpha ) where alpha is the angle of incidence
92
#####
RN = RN + RM * ( 2 * cnjump - RM * csjump ) / c ** 2
93
-
94
#####
SELECT CASE ( Beam%Type( 2 : 2 ) )
95
-
CASE ( 'D' )
96
#####
RN = 2.0 * RN
97
-
CASE ( 'Z' )
98
#####
RN = 0.0
99
-
END SELECT
100
-
101
#####
ray2D( is1 )%c = c
102
#####
ray2D( is1 )%tau = ray2D( is )%tau
103
#####
ray2D( is1 )%p = ray2D( is )%p + ray2D( is )%q * RN
104
#####
ray2D( is1 )%q = ray2D( is )%q
105
-
106
-
! amplitude and phase change
107
-
108
#####
SELECT CASE ( HS%BC )
109
-
CASE ( 'R' ) ! rigid
110
#####
ray2D( is1 )%Amp = ray2D( is )%Amp
111
#####
ray2D( is1 )%Phase = ray2D( is )%Phase
112
-
CASE ( 'V' ) ! vacuum
113
#####
ray2D( is1 )%Amp = ray2D( is )%Amp
114
#####
ray2D( is1 )%Phase = ray2D( is )%Phase + pi
115
-
CASE ( 'F' ) ! file
116
#####
RInt%theta = RadDeg * ABS( ATAN2( Th, Tg ) ) ! angle of incidence (relative to normal to bathymetry)
117
#####
IF ( RInt%theta > 90 ) RInt%theta = 180. - RInt%theta ! reflection coefficient is symmetric about 90 degrees
118
#####
CALL InterpolateReflectionCoefficient( RInt, RefC, Npts, PRTFile )
119
#####
ray2D( is1 )%Amp = ray2D( is )%Amp * RInt%R
120
#####
ray2D( is1 )%Phase = ray2D( is )%Phase + RInt%phi
121
-
CASE ( 'A', 'G' ) ! half-space
122
#####
GK = omega * Tg ! wavenumber in direction parallel to bathymetry
123
#####
gamma1Sq = ( omega / c ) ** 2 - GK ** 2 - i * tiny( omega ) ! tiny prevents g95 giving -zero, and wrong branch cut
124
#####
gamma2Sq = ( omega / HS%cP ) ** 2 - GK ** 2 - i * tiny( omega )
125
#####
gamma1 = SQRT( -gamma1Sq )
126
#####
gamma2 = SQRT( -gamma2Sq )
127
-
128
#####
Refl = ( HS%rho * gamma1 - rho * gamma2 ) / ( HS%rho * gamma1 + rho * gamma2 )
129
-
130
-
! write( *, * ) abs( Refl ), c, HS%cp, rho, HS%rho
131
#####
IF ( ABS( Refl ) < 1.0E-5 ) THEN ! kill a ray that has lost its energy in reflection
132
#####
ray2D( is1 )%Amp = 0.0
133
#####
ray2D( is1 )%Phase = ray2D( is )%Phase
134
-
ELSE
135
#####
ray2D( is1 )%Amp = ABS( Refl ) * ray2D( is )%Amp
136
#####
ray2D( is1 )%Phase = ray2D( is )%Phase + ATAN2( AIMAG( Refl ), REAL( Refl ) )
137
-
138
-
! compute beam-displacement Tindle, Eq. (14)
139
-
! needs a correction to beam-width as well ...
140
-
! IF ( REAL( gamma2Sq ) < 0.0 ) THEN
141
-
! rhoW = 1.0 ! density of water
142
-
! rhoWSq = rhoW * rhoW
143
-
! rhoHSSq = rhoHS * rhoHS
144
-
! DELTA = 2 * GK * rhoW * rhoHS * ( gamma1Sq - gamma2Sq ) /
145
-
! &( gamma1 * i * gamma2 *
146
-
! &( -rhoWSq * gamma2Sq + rhoHSSq * gamma1Sq ) )
147
-
! RV( is + 1 ) = RV( is + 1 ) + DELTA
148
-
! END IF
149
-
150
#####
if ( Beam%Type( 4 : 4 ) == 'S' ) then ! beam displacement & width change (Seongil's version)
151
-
152
#####
ch = ray2D( is )%c / conjg( HS%cP )
153
#####
co = ray2D( is )%t( 1 ) * ray2D( is )%c
154
#####
si = ray2D( is )%t( 2 ) * ray2D( is )%c
155
#####
ck = omega / ray2D( is )%c
156
-
157
#####
a = 2 * HS%rho * ( 1 - ch * ch )
158
#####
b = co * co - ch * ch
159
#####
d = HS%rho * HS%rho * si * si + b
160
#####
sb = sqrt( b )
161
#####
cco = co * co
162
#####
ssi = si * si
163
-
164
-
! LP: Missing divide by 0 check on si from 2D version re-added.
165
#####
IF ( si /= 0.0 ) THEN
166
#####
delta = a * co / si / ( ck * sb * d ) ! Do we need an abs() on this???
167
-
ELSE
168
#####
delta = 0.0
169
-
END IF
170
-
171
#####
pdelta = real( delta ) / ( ray2D( is )%c / co)
172
-
ddelta = -a / ( ck*sb*d ) - a*cco / ssi / (ck*sb*d) + a*cco / (ck*b*sb*d) &
173
#####
-a*co / si / (ck*sb*d*d) * (2* HS%rho * HS%rho *si*co-2*co*si)
174
#####
rddelta = -real( ddelta )
175
#####
sddelta = rddelta / abs( rddelta )
176
-
177
#####
print *, 'beam displacing ...'
178
#####
ray2D( is1 )%x( 1 ) = ray2D( is1 )%x( 1 ) + real( delta ) ! displacement
179
#####
ray2D( is1 )%tau = ray2D( is1 )%tau + pdelta ! phase change
180
#####
ray2D( is1 )%q = ray2D( is1 )%q + sddelta * rddelta * si * c * ray2D( is )%p ! beam-width change
181
-
endif
182
-
183
-
ENDIF
184
-
END SELECT
185
#####
END SUBROUTINE Reflect2D
186
-
187
-
END MODULE ReflectMod

Legend:

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