Generalization of the parabolic equation for EM waves in a dielectric layer of nonuniform thickness

Wave Motion 17 (1993) 337-347 Elsevier

337

Generalization of the parabolic equation for EM waves in a dielectric layer of nonuniform thickness V.A. Baranov and A.V. Popov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation, Academy of Science, IZMIRAN, 142092 Troitsk, Moscow Region, Russia

Received 3 July 1992

The parabolic equationmethod by Leontovichand Fock is usuallyapplied to scalar wave equationsor to such vector problems that admitintroductionof scalarpotentials.Here, we seek an approximatedescriptionof three-dimensionalwavepropagationin a nonuniform dielectric layer. It is supposed that the characteristicthicknessB is large compared with the wavelengthA= 21r/k, and small compared with the characteristicscale of horizontalnonuniformityL. In manycases of practical interest,the Fresnelparameter o-= kB2/L is of the order of unity.Underthis assumption,we constructan asymptoticsolutionof Maxwell'sequationsas a series in powers of the smoothness parameter v= B/L << 1. Its leadingtermcan be expressed via solutionsof two parabolicequationsfor slowlychangingscalaramplitudes.

1. Introduction The Leontovich-Fock wave parabolic equation ( W P E ) [1,2] has been used widely for the calculation of narrowangle-band wave packets (elg. radio wave propagation along the earth surface, microwaves in quasioptic devices etc). The transition from rigorous wave theory to the parabolic equation approximation can be justified by the physical conception of transverse diffusion [3] as well as by formally expanding the solution to be found or the wave operator in asymptotic series in negative powers of the wavenumber [ 4 - 6 ] . It must be noted that most of these papers deal with scalar waves or with such vector problems that can be reduced to scalar potentials. Our purpose is to extend the parabolic equation method to the case of EM propagation in a three-dimensional nonuniform environment when one has to solve the complete set of Maxwell's equations. More exactly, we seek an approximate description of EM wave packets gliding along a dielectric layer of nonuniform thickness b ( X , Y). It is supposed that the characteristic thickness B ~ b ( X , Y) is great compared with the wavelength A = 2~r/k and small compared with the characteristic scale of the horizontal nonuniformityL. It means that v = B / L and p = (kB) - 1 are small parameters of the problem we treat. In many cases of practical interest (for instance, VLF propagation in the earth-ionosphere waveguide) the order of magnitude of these parameters are such that the ratio o-= v / p = k B 2 / L is of the order of unity. Under this condition, making use of the well-known technique of the ray method [ 7 ], we construct an asymptotic solution of the Maxwell equations as a series in powers of the smoothness parameter v = B~ L << 1 (a similar approach has been developed in [8] for underwater acoustic propagation). This solution has the features of geometric optics in the horizontal plane (X, Y) cf. [9] and depends on the vertical coordinate Z via two series of scalar amplitude functions governed by partial differential equations that can be reduced to the LeontovichFock WPE. 0165-2125/93/$06.00 © 1993- ElsevierSciencePublishersB.V. All rights reserved

V.A. Baranov, A. V. Popov / Parabolic equations for EM waves

338

2. P r o b l e m

formulation

and asymptotic

solution

Consider harmonic ( ~ e-i,o,) EM waves propagating in a nonuniform dielectric layer 0 < Z < b(X, Y) and satisfying the Maxwell equations rot E = ikH,

rot H = - ikeE

( 1)

where k = o9/c is the free-space wavenumber, and the dielectric permittivity e, like the layer thickness b, is supposed to be a smoothly varying function of horizontal coordinates (X, Y). At the curvilinear upper boundary Z= b(X, Y), let the tangential components of the electric and magnetic fields be connected by the Leontovich condition E, - E - n ( n E )

= -_~ n ×H

(2)

with a normal impedance ~ ( X , Y) depending smoothly on (X, Y). The same boundary condition E, = ~ o no × H

(3)

is supposed to describe the lower boundary Z= 0 (we allow it to be flat just for simplicity). Here,

bx n= ( -~/l+bZx+bZ v, ~ / l ++b b~ z ' - b v

v / l + b Z1 x +) b 2

,

n o = ( 0 , 0 , 1)

(4)

are the normal to the surfaces Z = b(X, Y) and Z= O. Assuming that the characteristic value of the layer thickness b(X, Y) is B and the characteristic scale of its variations along the horizontal plane is L >> B, we introduce dimensionless variables x = X/L, y = Y/L, z = Z/B so that b = Bh(x, y), E= e(x, y), z,= B/L << 1 and the Fresnel number o-= kB2/L ~ 1. In order to simplify calculations we use a 6-component vector of the electromagnetic field

E, ll=

Ez

(5)

nx Hv

Hence, the Maxwell eqs. ( 1 ) can be rewritten in the following form

2[ all all) all =io'Rll vtL--~x +M---~ + v N a-----Z

(6)

where L, M, N, R are block matrices

l=

(i ° i) (°°i) 0

l

-

,

m=

0

0

- 1

0

,

n=

{0, i) {i°i) 1

0

~0

0

,

I=

1

.

(7)

0

We search for asymptotic solutions of the singularly perturbed eq. (6) having derivatives with respect to the horizontal variables (x, y) of the order O( v -2) and of the order O( 1) with respect to z variable. These solutions

V.A. Baranov, A. V. Popov /Parabolic equations for EM waves

339

correspond to the lowest-order modes of the dielectric layer propagating along the (X, Y) plane via successful, nearly glancing, reflections between the layer boundaries. After separating the rapidly oscillating exponent:

lI(x, y, z) =U(x, y, z) exp[kr/v24(x, y)]

(8)

with the phase function 4(x, y) not determined yet, one obtains for the slowly varying amplitude vector U(x, y, z) the following equation

io'(R-4xL-4yM)U=

vN-~z + v 2 L-~x + M

.

(9)

It can be solved in the form of an asymptotic expansion in the small parameter v:

U(x, y, z) = Uo +

v U i 3I- /22U2 -4- ...

(10)

Substitution of the series ( 1O) into eq. ( 9 ) yields a recurrent chain of linear algebraic equations for the coefficients

E: SUo=O

SU~ = 1 N OUo ,~ -Cz - F I ,

'

1 [.. OU,

sv2=

SU 3 = -1 10r

+'

OUo

OUo'~ _

+M Oy ) - , 2 ,

(NOU2 "~-L ~U1 + M ~ ) = _ F 3 -~Z

OX

•..

(11)

Here, the matrix S defined as

S = R - 4xL - 4yM

(12)

depends on the unknown phase function 4(x, y). It is obvious that nontrivial solutions Uo(x, y, z) of the zeroth-order eq. ( 11 ) exist only if det S = ~ ( 4 2 + 4 2 - E) 2 = 0.

(13)

This means that the phase function 4(x, y) must satisfy the following two-dimensional eikonal equation: (V4)~=~(x'Y)

~x'

"

(14)

One can check easily that under the condition (13), the matrix S has two annulling eigenvectors 0 0 1

a=4y

'--4y 0

,

b--

0

(15)

o ff

so that the general solution of the zeroth eq. ( 11 ) has the following form

Uo(x, y, z) =u(x, y, z)a(x, y) +v(x, y, z)b(x, y)

(16)

V.A. Baranov, A. V. Popov / Parabolic equations)or EM waves

340

where u and v are scalar functions to be determined. In order to obtain them, one has to refer to eqs. ( 11 ), extended for the higher-order corrections Uj, j = !, 2 . . . . . All of them have the same algebraic structure

SUj =Fj

(17)

where S is the matrix (12) and F:, at each step, is an already known vector function depending on the preceding solutions U:-i, Uj 2. By virtue of eq. (13), the matrix S is twice degenerated, so the inhomogeneous eq. (17) has a solution only under the condition [ 10] that the right-hand part Fj is orthogonal to the eigenvectors of the transposed matrix s'r: q~,.

0 0 E

~=

0

'

/3=

_q~

o

q,x

1

0

(18)

It is easy to check that the solution of eq. (17), under the condition (aUg'j) = ( f l F j ) = 0 ,

(19)

can be represented in the form

U: = uja + %b + TFj

(20)

where T is the block matrix

t=

1 0

,

(21)

and uj(x, y, z) and vj(x, y, z) are arbitrary scalar functions. After the preceding amplitudes U:_ l, Uj_2 inherent in the right-hand side of eq. (17) have been found except for the scalar functions uj _ i, c) l, u~_ 2, c) _ 2, the orthogonality conditions (19) yield a pair of equations connecting these coefficients not determined so far. Here, we will obtain the equations governing the coefficients of the zeroth-order approximation (16). After writing down explicitly the right-hand side F~:

F~ = - -

lo"

Na+

7z

Nb

(22)

and taking into account the following identities (a, N a ) = (a, Nb ) = ( fl, N a ) = ( fl, Nb) = 0

(23)

one can see that, f o r j = 1, the orthogonality defined by eq. (19) is ensured unconditionally. The amplitudes u(x, y, z) and v(x, y, z) remain, at this step, undetermined. One could be worried by the fact that the expression for the first-order correction

UI = u l a + v l b +

--

10"

TNa+

-~Z

TNb

(24)

ensued from eq. (20) contains already four arbitrary functions: u, v, Ul, c~. Nevertheless, by constructing the following corrections, this ambiguity will remove.

341

V.A. Baranov, A. V. Popov /Parabolic equationsfor EM waves

In fact, let us calculate the right-hand side of eq. (17) for U2: F2 =

k--~z Na +

+u L~x

Nb +

Oy]

L+

L+

a+

~ ~x

- --~I~z2NTNa+

-NTN2boz

(25)

.

Employing the identities (23) and explicit formulas for the remaining scalar products (a, NTNa) = 1, a , L Tx = 4 , ~ ,

(a, La) = 2q~x,

(a, m a ) = 2q~y,

a, g

(or, N r N b ) = ( a , Lb) = ( a , g b ) = 0

(a, L ~ x ) = ( a , M ~ ) = 0 ,

(fl, N T N a ) = ( f l , L a ) = ( f l , M a ) = O ,

(/K L-~x) =(fl, M ~ ) = 0 ,

(fl, NTNb) = e,

(fl, L a~--~)=,~xx+ G~bx,

(fi, Lb) = 2eqbx,

(fl, Mb) = 2eq~y, (26)

([~,M~)=~t~yy-}-~yCI) v ,

appearing in eq. (19), we see that the solvability condition of eq. (17) for j = 2 can be reduced to a pair of independent equations i o" 2 ~X -~x + dp~ -

i~r 2E qJX~x +qJy

"~-( ~xx "[- (I)yy! U q - - OZ2 ~- 0 '

+[e(qG+qbyy)+Gq~x+Eyq~y]v

+ -0z - 2 =0,

(27)

for the scalar potentials u(x, y, z) and v(x, y, z). They can be rewritten in a more compact form as io'(2V~Vu+uA4))+

02u = 0

3Z2

'

io" 2 V ~ V v +

Aq)+

v + -

0Z2

=0

(28)

by using the notations V= (0/0x, O/Oy), A = 02/Ox 2 + 02/Oy 2 for the horizontal gradient and the Laplace operator in dimensionless the variables (x, y). Quite similarly, the solvability condition of eq. (17) forj = 3 yields partial differential equations governing the first-order corrections u~ (x, y, z), Vl(x, y, z). Omitting the rather complicated derivation, we just write down the result:

02Ul = Ex~y--~v~1~x OU io'(2 V~VUl "}-U1 ~7~)) -{- - Oz2

[

(

io" 2Vcl, Vv, + V ~ +

~TEV(J)~ ]

•

02UI

-~ ] U l J ''[- ~OZ

Oz

--

~';Y~X--~x~yOU if2 OZ"

(29)

In summary, we can state that to determine the asymptotic solution of Maxwell's equations one has to solve the two-dimensional eikonal equation (14) and the parabolic equations (28). Their solutions ~(x, y), u(x, y, z), v(x, y, z) inserted into eqs. (8), ( 15)-(16) yield the leading term of the 6-component EM wave propagating along the

342

V.A. Baranov, A. V. Popov / Parabolic equations for EM waves

nonuniform dielectric layer. The first-order corrections (with respect to the smoothness parameter v) are given by the solutions of eqs. (29) that are to be substituted into eq. (24).

3. Derivation of approximate boundary conditions The approximate boundary conditions for the parabolic equations (28)-29) can be deduced from equations ( 2 ) (3) by expanding them in the small parameter 1,. Their concrete form depends on the magnitude of the impedances S~ and -9~0. We shall confine ourselves to the typical case when these quantities are small (this is always true, for instance, for non-magnetic boundaries if the Leontovich condition is valid [2] ). For simplicity of the formal derivation, we will consider the impedances to be proportional to the first power of the small parameter v:

_9~= vg(x, y) ,

~q~o= vgo(x, y) .

(30)

Let us insert the asymptotic solution of Maxwell's equations obtained above into the boundary condition (2). After truncating the series (10) at the quantities of the order v 2, one obtains an approximate representation of the EM field 6-vector

//=eiO./~,2~ (u+ vul + v2u2)a+(v+ vv, + v2v2)b+ V__,0- 0Ulox + v-~]Hva~-~ Ox + v

+ --10- L ~ + ~ 122 0 -2

a+ 3x L+

b+u L-~x +M 3y]

[32U 32V "~'~ rNrN~-~z~ a + -0Z- 2 bj f +O(,,3).

TNb

L-~x +M

(31)

From here, easily ensue formulas for the electric and magnetic field vectors E and H. It is convenient to write them down in a local cartesian coordinate set formed by the basis vectors

Or 0)

'

q=(

Or 0)

'

,=,001,

It follows from eqs. (2), (30) that, intending to obtain boundary conditions for u, v, u~, v~, one should have an approximation of the electric field vector accurate up to the quantities of the order of ~:

E=

ioq/-~ -~z + ~--~z ] -['- ~

U+I~

1 + /'F2U2

10-EqV(ev) p

0 "21~-~aZ

) + ~pV(Ev)

q

+ (u+ VUl + v2u2)k}ei'~/~2a~+O(v)3

(33)

and of the magnetic vector - up to O( v 2) :

H=(

vlf~E az p - V/-e(u+ vu~)q+'(v+ w l ) k } ei'~/~2a~+O(v2) .

(34)

343

V.A. Baranov, A. V. Popov /Parabolic equations for EM waves

b(X, Y) = Bh(X/L, Y/L):

Let us calculate the EM field components orthogonal to the surface z =

E, - E - (nE)n =

e ~ / ~ 2 ~ [ 1 + v2(Vh) 2] - l

- v2~/-~)(q~)v

- .---qg(~v) p lore

"~-[~fEE(U'+" /~, Jr- 1/2U2) "1- v(qVh)(u+ vuj)

1,2

p2 a2v]

+ m/r~(qVh)(v+ ~ , ) + v2(Vh)2u - ~

(pVh) -~z..I . ) '

(35)

n X H - ~ e~'/~2~[ 1 + v2(17h) 2] - ~/2 X{[~/'~E(u+

vul)- v,(qVh)v]p+[ ve(pVh)v

vv/-~~z q+ vV/~(pVh)uk}

Substituting the vectors (35) into eq. (2) and projecting ontop and q yield two scalar identities

_Tz

)-

,,,

= --g[~/~E(//'Jt - 1,q./l)- pE(q~Th)v] , I/2 ~ ( ~ + w , + v~v~) + v ( q ~ ) + v ~ ( p v h ) ~ v + r--p~z(Ev)

1.,2 ~2 v

~o-~

or2.¢/~ 0Z 2 (36)

After arranging them in powers of v, one gets the following boundary conditions for the potentials u, v

O~ - i ~ ~z

+

u=O ,

v---O

(37)

and for the first-order corrections u j, vz

__ ~ul

-io'~

g+

pVh

at the curved upper surface z = 0u ~z +i~regou= 0 ,

uj

=-~e(qVv)

h(x, y). Quite u=O

for the leading term coefficients and

-~ - ~ ( q V h )

~z ,

v,

=

u

(38)

similarly, at the fiat bottom z = 0 analogous boundary conditions arise: (39)

V.A. Baranov, A. V. Popov/Parabolic equationsfor EM waves

344

Ou~ +i~regoul = - V / - e ( q V v ) = O , ~Z

L'l = 0 .

(40)

for the first corrections. Formulation of the initial condition in the vicinity of a source will not be discussed here. It has been considered in detail in [2,4].

4. Analysis of the results In the previous sections, an asymptotic solution of Maxwell's equations has been constructed that describes propagation of EM waves along a uniform dielectric layer. The layer is characterized by a slowly varying thickness b( X, Y ) = Bh( X / L, Y/ L ) , as well as a permittivity e( X / L, Y/ L ) and boundary impedances £U ( X / L, Y/ L ) , 2 o ( X / L, Y/L). The assumption of slow variation of the layer parameters in the (X, Y) plane ( v = B / L << 1 ) made us able to reduce approximately the exact three-dimensional vector problem ( 1 ) to the scalar parabolic eqs. (28). Now we will demonstrate that these equations are actually two-dimensional and correspond to an approximate decomposition of the total EM field in TM and TE waves. Let us consider a solution q~(x, y) of the two-dimensional eikonal equation (14) determined by an initial wave front. We introduce curvilinear coordinates ( ~', ~) associated with the wave fronts @(x, y) = Const and horizontal rays - characteristic lines of eq. (14). Suppose that the variable ~-varies along the rays and is normalized so that (ds) 2 = e(x, y) (dT)2, and ~:is constant on the rays so that the vectors p and q introduced above are the basis vectors of this coordinate set. It is seen from eqs. ( 3 3 ) - ( 3 4 ) that the vectors E and H are almost orthogonal to the average propagation directionp, and small longitudinal projections are proportional to ~u/Oz and ~v/Oz being associated with vertical nonuniformity of the EM field. Since the potentials u(x, y, z), v(x, y, z) satisfy the independent linear differential eqs. (28) and boundary conditions (37), (39), there exist two approximately independent classes of solutions. The first one

(

ETM = uk

io_gt-e Oz p e ~'~/"2~,

HTM = --UV/--eq e ~'/'2~t"

(41)

may be called the TM mode while the second one

E T E - ~ C ~ q e i~T/~2~,

live ~-

c,k- io.vr~~ ~zp e i~/~T24~

(42)

is approximate the TE mode because the vector q is transverse to the propagation direction. It must be emphasized that variations of the layer thickness b = Bh(x, y) and permittivity e(x, y) are not supposed to be small - TM and TE waves separation is possible due to the adiabatic character of these variations. As it can be seen from the differential equations (29) and the boundary conditions (38), (40), small coupling between TM and TE modes arises in the first-order approximation (with respect to v) when transversal gradients of the permittivity q Ve or the thickness q Vh are encountered. By employing the ray coordinates (~-, ~), the parabolic eqs. (28) governing the potentials u(x, y, z), v(x, y, z) can be rewritten in equivalent form as

(uvS)

2io'--

0~-

02(u¢5) =0,

+ -

0z 2

2 i o 3t

+

3z 2

=0

(43)

V.A. Baranov, A. V. Popov / Parabolic equations for EM waves

345

where J = D(x, y ) / D ( ~', ~) is the ray divergency [ 11 ]. Since E and J on the ray sc= Const are known functions of the variable ~-, substitution of

u=

A(x, y, z) f),

v=

B(x, y, z) ~

(44)

will reduce the problem to integration of the Leontovich-Fock WPE OW a 2 W

2i~ aT az 2 = 0 ,

(45)

which governs both functions A and B. Equation (45) has constant coefficients. Nevertheless, in most nontrivial cases, the nonuniformity of the layer boundary and of the impedances inhibits separation of variables. Therefore, application of numerical methods for the EM field calculation in a nonuniform dielectric layer is practically inevitable. The asymptotic approach proposed above can be very useful for the numerical implementation because it reduces the three-dimensional vector problem to integration of two-dimensional marching-type equations for slowly varying scalar amplitude functions. Equation (42) shows that the asymptotic solution constructed above diverges in transition regions near caustics and points of focusing of the horizontal rays where J( ~-, ~) = 0. In these regions the amplitudes u and v vary rapidly and, for a uniformly valid approximation, the Ansatz (8) must be modified - cf. [ 11,12 ]. To finish our analysis, we will write down the integral relation reflecting energy conservation for waves approximately described by the parabolic equations ( 2 8 ) - ( 2 9 ) . Let us introduce a normalized density of the horizontal energy flow h(xv)

W(x, y ) = R e

f

[E×H*]lldz

(46)

0

Here, subscript II denotes projection of a vector onto the (x, y) plane: aii = a - (ka)k. By making use of eqs. ( 3 3 ) (34) we get h

W=p~

(lu+ vu, 12+Ely+ ~v, 12) dz+

~qf~[uv*-u*vl~:~+O(v2).

(47)

0

As a consequence of the boundary conditions (37) and (39), the terms outside the integral are equal to zero. From the parabolic eqs. (28) the following energetic identities follow

'[

h(x.z)

div( f / 0

ou'..ouT

2i--~ u az

-~zJz= o '

~ [ av*_v, a~]z=h

h(x,z)

elvl2dz =~

c - - ~z-

0

Analogous relations for the corrections u~ and v~

~zdz=o"

(48)

346

V.A. Baranov, A. V. Popov / Parabolic equations for EM waves

h(x,v) div [ Vq~

f(UU*I-}-U:~Ul)dZ]=(V~I)~Th)(RbI~q' -U*UI)'z=h o [ 1

+2i~

au* u ~z

u* ~ul +uj

-

--~-

-

OU :~

OU] z=h

-u*--

~

0ZJ,= 0

h(x,y)

+

Oz - u ~ - z ] dz,

2io-E

o

(49) h(x,y)

div [ V~

I

e(vv* +v*t'l) dz] =e( V~Vh)(vv* +c*v~)I~=,~

o

e [ Or*-v* + -2iov - Oz -

av -Oz - +vl

av* OvT'=h ~ - - t ' * - - OZj.=o

h(x.y)

+ exep~-e,,~x . . 2io-e

I(

.

0

.U*

au* . Oz

U~

az

dz,

can be derived from eqs. (28)-(29) considered jointly. After combining eqs. (46) and (47), we obtain an approximate expression for the divergence of the horizontal energy flow vector div W= - e ( R e gtu+ vu~ Izlz=h~x,y~+ R e go lu + uu1121:-o) + O ( p 2) •

(50)

It makes it possible to estimate the attenuation of TM waves due to absorption in the waveguide boundaries. TE waves have lower attenuation. To actually calculate the attenuation, one should use a more accurate approximation of the energy flow (46) - up to quantities of the order O(v3).

5. Conclusion

The asymptotic solution of the Maxwell equations constructed in this paper provides the simplest description of the EM waves propagating along a nonuniform dielectric layer. The leading term of the phase ~b(x,y) is determined by the inhomogeneous dielectric permittivity E(x, y). The nonuniform thickness h(x, y) affects only the global distribution of the complex wave amplitude U(x, y, z) calculated along the horizontal rays via numerical integration of two-dimensional parabolic equations for scalar potentials. In contrast to the adiabatic approximation [ 9 ], grazing wave diffraction effects disturb the separation of individual vertical modes. The corrections connected with global curvature of the layer and weak vertical inhomogeneity of the dielectric permittivity e= e-o(X,y) + u~ (x, y, z) can be found within the framework of the method described above. Field calculations in the regions of horizontal ray focusing and at great distances (D >> L) from the source requires another approach.

References

[ 1 ] M.A. Leontovich and V.A. Fock, "Solution of the problem of propagation of electromagnetic waves along the earth's surface by the method of parabolic equations", J. Phys. USSR 10 ( 1 ), 13 (1946).

V.A. Baranov, A. V. Popov / Parabolic equations for EM waves

347

[ 2 ] V.A. Fock, Electromagnetic Diffraction and Propagation Problems, Pergamon, London (1965). [ 3 ] G.D. Maluzhinets, "Development of conceptions on diffraction phenomena", ( in Russian), Soviet Physics (Uspekhi) 69 (2), 321 (1959). [4] F.D. Tappert, in: Wat,e Propagation and Underwater Acoustics, Ed. by J.B. Keller and J.S. Papadakis, Springer, Berlin-Heidelberg-New York, Ch. 5, (1977). [5] E.A. Polansky, "On the relation between solutions of Helmhotz and Schr6dinger type equations", (in Russian), J. Comp. Math. and Math. Phys. 12 ( 1), 241-249 (1972). [6] A.V. Popov and S.A. Hoziosky, "On a generalization of the diffraction parabolic equation", (in Russian), J. Comp. Math. and Math. Phys. 17 (2), 527-533 (1977). [7] V.M. Babi~ and V.S. Buldyrev, Short-Wauelength Diffraction Theory (Asymptotic Methods), Springer, Berlin-Heidelberg-New York, ( 1991 ). [8] G. Kriegsman, " A multiscale derivation of a new parabolic equation which includes density variations", Computers and Math. with Applications 11 (7/8), 817-821 (1985). [9] R. Burridge and H. Weinberg, in: Wave propagation and UnderwaterAcoustices, Ed. by J.B. Keller and J.S. Papadakis, Springer, BerlinHeidelberg-New York, Ch. 3, (1977). [ 10] R. Courant and D. Hilbert, Methods of Mathematical Physics, v. 1. Wiley, Interscience, New York (1953). [ I 1] Yu.A. Kravtsov and Yu.I. Orlov, Geometrical Optics oflnhomogeneous Media, Springer, Berlin-Heidelberg-New York (1990). [ 12] V.A. Borovikov and A.V. Popov, "Short-wave modes for wave ducts in quasistratified media", Wat,e Motion (6), 525-546 (1984).