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