The lagrange padé parabolic equation (LP-PE) for the prediction of long range sound propagation in the atmosphere

Applied Acoustics, Vol. 49, No. 2, pp. 105-125, 1996 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0003-682X/96/$1 5.00 + 0.00 ELSEVIER

0003-682X(95)00064-X

The Lagrange Pad6 Parabolic Equation (LP-PE) for the Prediction of Long Range Sound Propagation in the Atmosphere R. A. Sack & M. West The Telford Institute of Acoustics, University of Salford, Salford, UK M5 4WT (Received 15 August 1995; revised version received 6 November 1995; accepted 15 November 1995)

ABSTRACT A new formulation of the parabolic equation (PE) method for the solution of the two-dimensional sound wave propagation in realistic atmospheres is derived. The new PE simultaneously reduces the errors implicit in the parabolic approximation and those due to finite step sizes in the range variable, thereby allowing large range steps to be used leading to much higher computational speeds than the standard Crank-Nicholson (CN) method. The new algorithm is based on an expansion of implicit functions due to Lagrange and can be interpreted as a sequence of CN-like steps with varying parameters. Copyright 0 1996 Elsevier Science Ltd Keywords: Parabolic equation, long range sound propagation.

1 INTRODUCTION In the numerical solution of the two-dimensional

Helmholtz equation

by means of a parabolic partial differential equation (PE) method, the wave function T$ is factorized into a rapidly varying exponential ‘carrier’ and a slowly varying ‘modulator’ cp Ii/(x,z) = exp (ikox)dx, z) 105

(2)

R. A. Sack.

106

The differential

equation

M.

Wesr

‘ g,,,P cp=o [ 1

for cp becomes

$9 j$ + 2&, g +

The last term in eqn (3) is assumed to be small so that k varies slowly with x and the term a2v/dx2 can be ignored to first approximation. This is equivalent to saying that the wave number k varies slowly with the co-ordinates, so that the constant kO is always a close approximation to k and k* - ko2 is small. In assessing the accuracy of any PE method it is useful to classify various errors which can arise and should be minimized. They are:

(A) Errors due to replacing

the original

elliptic

PDE’

the

by its parabolic

formulation. due to taking discrete steps in the marching direction2 (of increasing abscissa x). CC)Errors due to the discretization of the vertical variable, z.~ (D) Problems in dealing with infinite or semi-infinite domains, errors due to taking only a finite range in Z; this involves a careful choice of boundary conditions.4.5 (El Problems relating to the choice of the initial conditions at the starting value for x.4.’ problems in curvilinear co(F) When dealing with two-dimensional ordinates equations arise of the form

(W Errors

(4)

where P and R may vary with x, problems arise from their non-constancy.’ Not much gain is achieved if only one of the errors is minimized and others, of the same magnitude, left uncorrected. Thus, the treatment by Collins2 goes to considerable trouble to reduce errors of type (A) without reducing those of category (B). Likewise, it is clear that improvements in (C) and (D) beyond the standard Crank-Nicholson (CN) procedure can only yield marginally better results without efficient methods regarding (A) and (B). The present paper is aimed at a simultaneous improvement of methods dealing with (A) and (B) on the assumption that P and R in eqn (4) are independent of x, so that the complications under (F) do not arise. This means that on formally expressing Taylor’s series as an operator exponential cp(s+ h, -_) = exp

(

/r$ >

60(x, z)

Prediction of long range sound propagation

107

where @ax satisfies (6) with, according to eqns (3) and (4) (k* - k; + &)

b=-R/P=-

Pko) E = -l/P

(7)

= -1/(2iko)

the operator a/ax can be treated like a scalar variable X, defined by A = b + rh*

(8)

Improvements relating to (D) and (E) are left to a later paper, and (C) enters into consideration only in so far as the essential numerical treatment must be broken down into individual steps linear in b, which, according to eqn (7) involves the operator #/az2. In Section 2 a series expansion is derived for exp(hX) in powers of eb and hb jointly; this is achieved by means of an expansion theorem due to Lagrange.6 This power series is then replaced by a Pad& approximant. The polynomials of b occurring therein are in turn factorized into linear factors. The coefficients entering into these terms have to be determined numerically. In Section 3 we examine the core algorithm using a one-dimensional formulation, which provides a qualitative understanding of the errors implicit in the method. A representative selection of the resulting coefficients for a wide range of h/e or k,,h is tabulated (see Table 1). The new method is also applied to the one-dimensional problem corresponding to the choice of $(x,z) = exp (ikx) in eqn (2), and the approximate value of k obtained is compared to its value given at the outset. In Section 4 the new approach is applied to the original PE problem [eqns (l)-(3)] and the results compared with those obtained by the standard Crank-Nicholson (CN) method. The improvements achieved to date and an outline of future developments are discussed in Section 5. 2 DERIVATION

OF THE PADti APPROXIMANTS

The Lagrange series expansion Lagrange’s theorem6 states that for a variable X defined by I=b+@(A)

(9)

K. A. Suck, M. West

10X

an analytic

functionf’(X)

can be expanded

as a power series in

(10)

(Note that here 4 denotes an analytic function of the single variable b or X; it must not be confused with the modulator cp in eqns (2) and (3); the convergence criteria for the series ( 10) are very involved and will not be considered here.) In the given case we have for the exponential operator in eqn (5) with the use of eqns (6) and (8)

which yields for ( 10) II- I ehh = ehh + 2

7

&

[b2”ehb]

(12)

I

or on expanding

the exponential

where C”(, = 1; C,” = 0

C’,, =

(0 ’ 0)

(2n+mI)! n!(m - I)!@ + n)!

(14) (m ’ 0)

These coefficients are tabulated in Table 1. As far as the applications with h as defined in eqn (7) are concerned, all that matters are the powers of b, so that eqn (13) is best expressed as

exp(hh) = c

W,bN .\‘=(I

(15)

where the coefficients W, are polynomials in t and h, homogeneous of degree N. These polynomials can be evaluated using the coefficients along

Prediction of long range sound propagation

diagonals with constant recurrence relation c

nm

109

n +m = N, but as the coefficients satisfy a simple

=Wn+2m-3)C_

+ n ‘3m

n+m

Gm-2

(n+m)(n+m-

1)

(16)

the equivalent for the IV, becomes W

N

= 2(2n - 3)

N

cw,-1

+

h2 WN_2

N(N-

1)

(17)

Beginning with W,, = 1, W, = h, W, = ch + h2/2

(18)

the higher coefficients are best evaluated by means of eqn (17). A more concise method of deriving eqn (17) is given in Appendix 1. It should be pointed out that error terms proportional to the first power in horizontal step size Ax = h cannot be reduced by making h smaller, as the number of steps needed for a given range is inversely proportional to h. Alternatively it follows from the expansion exp hX = 1 + hX + h2X2/2 .. . that errors in the term proportional to h must involve errors in X itself and are thus independent of the PE approximation. Thus, the emphasis placed by Collins* to get terms linear in h correct to as high an order in E as possible is in full agreement with the aim of minimizing errors of type (A); the present approach refines this aim by taking account at the same time of higher powers of h, to reduce errors of type (B). Conversion to Pad6 approximants

The next step consists of finding Padi approximants (15)

for the power series

(19)

The degrees of the polynomials are best kept identical; for L < M the same number M of linear equations involving b has to be solved for a lower overall accuracy; for L> A4 the numerical process is likely to lead to instabilities.

K. A. Suck. M. West

I IO

The coefficients (~1and &,., are best determined by multiplying eqn (19) by the denominator V,db) leading to a set of 2M linear equations CV,+ W/&,Bl + W/_& + ...WO/Y/= u/ (I < M) w/~/-M+W/~l~/_M,,+...WO~I=O

(Mtld2M)

(20)

For L =-A4 = 1 this yields, in view of (18) h+/l,

;h’+ht+hflI

=cyi:

=0

(21)

with the solution B, = -h/2

- t:

aI =h/2-t

(22)

This is the basis for the standard ‘wide angle‘ approximation. For L= A4 > 1 a larger set of linear equations has to be solved; thus for LXMZ2 h+h f’f”,

+

w2/31

=al; +

W2

hB2 = 0:

+

hB,

+

82

=

a2

w, + w,@, + w2B2 = o

(23)

This set can be solved by Kramer’s rule, but is better tackled numerically. For larger values of L = M a numerical approach appears unavoidable. It should also be noted that in view of cxp (hh) x exp (-hi)

not only do the coefhcients

= 1

c’,,, in eqn (14) satisfy c,,,n(c. -h)

= (-_)“‘C,.m(c h)

(24)

but likewise

U/.M(C,-h) = VML(C,h)

(25)

In general this is not of practical use as the individual coefficients CQand &, though still homogeneous functions of E and h, are no longer polynomials in these variables. However there is an exceptional case when h and t are pure imaginary,

and h is real

(26)

Prediction of long range soundpropagation

111

which is the case when eqn (1) is turned into a PE. Not only does this make X also pure imaginary, but all terms in odd powers of hb are imaginary, and the others real. Hence changing the sign of h turns U and V into their complex conjugates which for M= L yields for eqn (25)

where the asterisk denotes the complex conjugate. This has the double advantage of simplifying the expressions and ensuring that U and V have the same modulus, and thus their ratio has modulus unity in agreement with the property of the exponential of the imaginary argument hX. Factorization of Pad&approximants Lastly the denominator polynomial VMdb) must be factorized into linear terms, analogous to the partial fraction expansion V,,(b)

= fic1 +

Pi4

(28)

i=l

again the roots pi must be determined numerically. In the general case one would have to get the complete decomposition of the Pad& approximant by a modified cover-up rule, but in the relevant case (eqn (26)) this simplifies to

M 1 -p;b

UMM

l-I i=l 1 + pib

-= V&U

Equation (29) is valid as long as also imaginary (29) becomes UMh4 -=

VMM

c/h

is pure imaginary, if, in addition,

M (1 + P&j*

l-I , 1 + pib

(2%

b

is

(30)

This is applicable for the PE, provided that there is no air absorption. In the more general case (Re(e/h)#O), the polynomials UMM and VMM must be factorized separately UMM -=

VMM

M 1 + sib

l-I-, 1 + Pib

(31)

K. A. Such. M. West

112

finding the correspondence between the (T; and p, may provide some difficulty. A further useful modification is to write the Padi approximants in the form (32)

valid under the conditions applying to eqns (29), (30) and (31) respectively. This has the advantage that the coefficients &and&i depend merely on the quotient c/h, and not on the absolute values E and h. (In the PE c/h= I /(2ik&).)

3 A ON E-D1 M ENSIONAL

APPLICATION

In contrast to the case t = 0 in eqns (8) or (9) for which the Padi approximants for ehh are analytically given as confluent hypergeometric functions (Laguerre polynomials) of the variable h (see Appendix 2), in the general case (the coellkients of) the polynomials Cl and V in eqn (19) have to be determined numerically, as must the coefficients p and c in eqns (28~(32). Table 1 lists the coefficients for a representative set of values of e, which in turn are determined by the ratios h/X0,where (33) However, this (or any more expanded) table should only serve as an indicator of how the coefficients behave; to use fixed values of h/X, would deprive the method of its flexibility. An exception is the case M= 1 for which, according to (7) and (22) p1 = - ;/I -- ih”/(47r);

p, = -;

- &J/(47rh)

(34)

In general one should proceed numerically as explained in Section 2. Table 1 does, however, give a good indication of the behaviour of the coefficients for A4 > 1. We find

c

Re(fi,)=-::

c

CRe(pi)=-ih

(354

Re (&) Im (6;) = ho/(8rrh)

(35b)

Prediction of long range sound propagation

113

Thus, while twice the real parts of the pi can be interpreted as sub-intervals of h, the imaginary parts fluctuate widely and may take either sign. To test the accuracy of the method described in Section 2, we have applied it to the one-dimensional problem resulting from eqns (l)-(3) if $I is taken as independent of z. The true solution for $ in eqn (1) is always * = C

exp (ikx)

(36)

regardless of how it is split up according to eqn (2) to which the numerical solution exp(ik,,x)xcp serves as an approximation (the alternative solution exp(-ikx) is excluded by virtue of the PE approach). Hence any solution v(h) obtained from the initial conditions ~(0) = 1 is an approximation to exp [i(k-k,#z] and consequently k0 +

Im IIn dM/h = kcak

(37)

is an approximation to the analytic value k entering into the solution, eqn (36). Some representative results, obtained by double precision calculations for the differences kcalc-k/k0 are shown in Fig. l(a)(c). The reason why comparison is made between approximations with constant h/M, rather than with constant h, is that in the former case the same computational effort is required to progress a given distance x( >> X). A number of conclusions can be drawn from these results. (A) Regardless of the values of h and A4, the calculated values of k are always smaller than the theoretical ones as long as k > ko, and larger than the true values of k for k < k,,. (This cannot be inferred from the logarithmic graphs, but follows from a print-out of the corresponding differences.) (B) As k tends to 0, all the approximations yield bad results. This is a consequence of the divergence of the Lagrange expansion at k = 0, where it is not clear whether exp (ikx) should be considered as an analytic continuation of exp(ikG) or exp(-ikox). (A lower limit for which the approximation is adequate has not been found to date.) (C) As k tends to 2k,, and beyond, there is little improvement in accuracy if A4 is increased at constant h/M= X2. (D) For M= 1 the curves for h = X/S and h = X/32 are practically identical on the scales given; this follows from the error in y itself, rather than the error in ehJ’,as explained in Section 2.

H. A. Suck. M. West

114

TABLE 1 Coefficients bL from eqn (32) M’2

M=3

M=4

0

-0.5

-0.2500 zt 1443;

~0.1423iO.13581’ 0.2153

-0.0916*0.1157i -0.1584 i 0.04741

0.5

~-o.s--o.o398l

-0.2506 + 0.1096; -0.2494-0.1899;

-0.1433 +O.l052i -0.1403-0.18751 0.2164-0.04001’

-0.0928 + 0.0885i -0.0878-O. 17391 -0.1607+0.0124i

I

-0.5-0.07961

-0.2543 + 0.082% -0.2457-0.24731’

-0.1489 + 0.0868i -0.1302-0.26491’ -0.2209-0.08251’

-0.0992 + 0.07581’ -0.0724-0.27221’ -0.1698-0.01901’ -0.158550.15381‘

2

--0.5.-0.1592,

-0.27 14 i- 0.0390; em0.228660.3944i

-0.23 17-O. 18OOi -0.1681 +O.O618i -0.1002-0.4816i

-0.1893-0.08481’ -0.1162 + 0.0604i -0.1477-0.30771’ -0.0468-0.53621

4

-0.5-0.31831

m0.2264-0.4032i ~0.2015+0.0187i --0.0722-0.99321

-0.2030-0.22731’ -0.1430+0.03391 -0.1205-0.6721i -0.0335-1.1008i

8

-O.S-0.63661

--0.3429-O. 1904i -0.1571-1.6109i

-m0.2392-0.0644i -0.2006-0.89931 -0.0602-2.03541’

-0.1923-0.54101’ -0.1762-0.0174i -0.1028-1.4319i -0.0287-2.23041’

16

-0.55

0.356550.4575; --0.1435. 3.3oooi

--0.2613-0.2152i em0.1830- I .919Oi -0.0557-4.1 l42i

-0.201330.1146i -0.1764- 1.2015i -0.0954-2.94621’ -0.0269-4.48501’

32

-0.5-2.54651

mm0.3604-0.9579i -0.1396-6.64921’

-0.2688-0.48481’ --0.1769-3.92641’ m0.0543-8.2576i

-0.2115-0.28571‘ -0.1694-2.5051i -0.0928-5.95331’ -0.0262-8.98751’

I .2732i

XC,= reference wavelength. h = range step. M = order of the Padi approximant

4 TEST

CASES

The Lagrange Padi: (LP) expansion described here has been applied to the numerical solution of the two-dimensional field* eqn (3) where the propagation is assumed to be predominantly in the x direction. According to eqns (6) *Although the method is two dimensional it may be applied to an axially symmetric dimensional case by scaling the numerical results with a factor Jx.’

three-

Prediction of long range sound propagation

H -

161

H *

LAMDAo

/

115

2

-65-,o-.--__ -76

_:..

-

-60

_.-..

.i_ I

I

4

12

6 K /

lEi

Fig. 1. Accuracy

H -

test for the one-dimensional eqn (19) is 1 for --; 2 for ...;

is

20

KO

M Z LAMDAo

iii1

,’

8

approximation in eqn (37). Parameter 3 for --; and 4 for -.-..

M in

and (8) the algorithm developed in Sections 2 and 3 for X may now be applied to the differential operator @iax. The z discretization remains the same as in the standard CN method. The new implementation of the field solution is best performed using a set of sub-ranges with increasing order of the Pad& approximation and increase

in horizontal step size. The LP-PE sub-range configuration chosen for use at all frequencies is shown in Table 2. The optimum configuration for maximum speed while maintaining a specified numerical accuracy (say within 1 dB of the CN-PE) has not yet been determined. The extraction of a suitable algorithm for setting up this optimum sub-range configuration will be the subject of a later paper. The above configuration was chosen so that the near-source zone, which includes the region where the PE breaks down,’ for I’ > X, is covered by a unity factor Pade which is equivalent to our original CN-PE. We have to be especially careful in this region since spurious reflections at the upper boundary can easily be generated if the number of Padi factors becomes large at low frequencies. The remaining sub-ranges were chosen so that the computational speed change from one sub-range to the next was roughly logarithmic.

Sub-range

configuration

No. 01’Pa& ,firctor.v

TABLE 2 the computational speed in each sub-range the CN speed

showing

Start

i mj

Finish i WI)

AxiJ,

compared

with

Speed (CNSJ

1 3

I0

800

0.2

1

; 4

1400 800 2500

2500 1400

I0.000

0.8 3.2 8.0

25 IO

Prediction of long range sound propagation

117

Generally the LP-PE, like the CN-PE, is very sensitive to even tiny errors in the configuration of the upper boundary condition. An elaborate ray tracing procedure was used whenever rays could be made to intersect the upper atmospheric boundary to allow the upper boundary condition to be correctly set up for upward bound rays. The LP-PE algorithm’s accuracy has been assessed by comparison of its predictions with corresponding predictions using the CN-PE. The test cases chosen use ground and atmospheric parameters specified in some of the ‘benchmark’ cases in Ref. 7, viz. Frequency Ground impedance

f zo

Source height Receiver height Range

=R

=s

10 and 100 Hz (38.79, 38.41) at 10 Hz (12.81, 11.62) at 100 Hz 5m lm 10 km

The sound speed gradients, g, recommended are 0, -0.1, 0.1 m/s/m and an alternating gradient case (0.1 (m/s)/m up to 100 m and -0.1 (m/s)/m up to 300 m and 0 above 300 m). Note that our calculations do not include the small atmospheric attenuation used.7 Our model has not yet been configured for operation at frequencies above 400 Hz. Comparing the LP and CN-PE predictions we can see at both 10 Hz (Fig. 2) and 100 Hz (Fig. 3) that excellent agreement is obtained. At the larger ranges, particularly beyond 2.5 km where the LP-PE’s horizontal steps are 8X, we obviously lose the detailed spatial behaviour, but the attenuations obtained by the LP-PE still agree closely with those calculated using the CNPE. In the ‘upwind’ cases (Figs 2(c) and 3(c), where g= -0.1, because of the rapid attenuation with distance, the LP-PE is not able to use its larger range steps. Figure 4 compares the contour plots obtained with the LP and CN PE methods. It shows that for all practical purposes the plots are indistinguishable.

5 DISCUSSION

OF RESULTS

It has been shown in eqn (34) that under the condition (26) the Lagrange-Pad& decomposition of the PE over a single range step Ax = h can be interpreted as a sequence of wide angle Crank-Nicholson (WAC-N) steps over a number

R. A. Sack. hf. West

Fig. 2. Ten Hz benchmark test. (a) g = 0. (b) g = 0. I, (c) g = -0.1 and (d) g = 0. I. -0.1, dashed line is the free field curve. PE step length is X/S.

0. The

119

Prediction of long range sound propagation Lagrange-Pa&

__ _.___-

c._.................,... ..-. Si_’

-55; k.....-

,. ...iao” ..,._...... 2’ij,

%\ ?.

~~~..

--.

.

Parabo?lc

__

,.

.z60

$.g6..

._

~~~

^

Distance

Equation

.._ _-

..-. _____l---F

7~~ .......

a6F

y6.6

m

..I i

10’

(c)

r %

-30

-

----________

-------____

-40 -50 0

100

200

300

400

Clis5tOsqrCe

600 II

(d) Fig. 2 - contd.

700

BOO

900

101

Lagrange-Pade

10”

,

I

,

.

-

.

Parabolic I

1

Equation

I

.

r

”

”

”

.

.

.

,

-1

H-

op,.

”

*

”

G

”

Wide

IO0

ZCC

”

Angle

Parabolac

”

1

Equation

I 300

400

500 Distance

600

700

BOO

900

101

n

Fig. 3. One hundred Hz benchmark test. (a) K==0. (b) g= 0.1, (c) g= -0.1 and (d) g=O. I. --0. I, 0. The dashed line is the free field curve. PE step length is X/5.

Prediction qf long runge sound propagation

Lagrange-Pads

Parsbolic

Eouation

-----mm_____ -4o-50 IOi oiatanca

kl

c 8 ”

(d) Fig. 3 -

conrd.

122

H. A. Suck. M.

Range m Wdc angle parabolic

loo

200

300

400

500 Renpc

Fig. 4. Attenuation

WCS/

equalon

600

(HIPE)

700

800

900

I”

,,I

contour plot at 50 Hz t‘or still air. Contours in 3 dB intervals. flow resistivity is 122 kRayls. Source height is 2 m.

Ground

A4 of sub-intervals of h with widely varying effective wave numbers. These sub-intervals, however, may only be taken as virtual as the theory does not prescribe in which order they are to be taken, and if, in addition, any of the parameters vary with x, the results would depend on the order selected. Even in the examples discussed in Section 4 in which all the parameters are independent of x, the Sommerfeld radiation condition used in determining the relation between the values of ;p(x, z) in the top-most positions of z, is affected by the value of x.’ Where waves with rapidly varying energy impinge on the top layer, the Sommerfeld matching has not always proved adequate and small but noticeable reflections can occur, leading to interference with upward radiation at larger ranges. It is these reflections which, so far, have vitiated the potentially most powerful improvement arising from the new approach. At large angles (zzx), the results obtained by the WAC-N method are qualitatively wrong, yielding contour lines that are convex towards the x-axis instead of concave. This is due to the fact that the factorization (2) is a bad separation in a region in which the wave propagation is essentially upwards rather than horizontal. Collins’ has shown that a Pad& decomposition reduces errors at

Prediction of long range sound propagation

123

large angles. However, when the new method is applied at short distances with values of M > 1, long spikes appear in the contour lines at very large angles. These spikes are of no consequence in the region where the results are qualitatively wrong anyway, but their reflections exacerbate interference effects in regions which give perfectly smooth results when M is initially taken as 1 (WAC-N approach) and gradually increased. The main improvement achieved to date by means of the new method lies in the speed of calculation at large distances. With M=4 step lengths of 8X,, or more can easily be accommodated leading to a speed-up factor of 10 or more compared with the WAC-N method. Where the sound profile can be calculated theoretically the errors arising by the use of the two methods are essentially the same. The sub-range configuration turns out not to have a critical effect on the solution accuracy. For a wide variation of sub-range parameters the values of the predicted sound pressure level minima were consistent to within 0.2 dB. Likewise complete consistency was obtained (to within 0.02 dB) when the height of source and receiver were interchanged. Further research is in progress to find methods of reducing the errors (C)(E) of Section 1 and to adapt the new method for case (F).

ACKNOWLEDGEMENTS The authors wish to thank the Directorate of Health and Safety, UK Ministry of Defence and the US Army Atmospheric Research Laboratory, White Sands for their financial support of this work. The authors also wish to thank the referees for their particularly helpful and constructive comments which have greatly improved this paper.

REFERENCES 1. Sack, R. A. & West, M., Representation of elliptic by parabolic partial differential equations with an application to axially symmetric sound propagation. Applied Acoustics, 37 (1992) 141-149. 2. Collins, M. D., Applications and time-domain solution of higher-order parabolic equations in underwater acoustics. Journal of the Acoustical Society of America, 86(3) (1989) 1097-l 102. 3. Gilbert, K. E. & White, M. J., Application of the parabolic equation to sound propagation in a refracting atmosphere. Journal of the Acoustical Society of America, 85(2) (1989) 63C637. 4. Gilbert, K. E. & Di, X., A fast Green’s function method for one-way sound propagation in the atmosphere. Journal of the Acoustical Society of America, 93 (1992) 712-714.

R. A. Sack. M.

124

West

of an algorithm for prediction of the sound field from a spherical acoustic source using the parabolic approximation. In Proceedings of the 5th International Symposium on Long Range Sound Propagation. The Open University, Milton Keynes, 1992, pp. 11S-127. 6. Lagrange, J. L., MPmoires de I’Academie de Berlin, 24 (1770) Oeuvres 2, p. 5.

West, M. & Sack, R. A., Development

25. [See also Whittaker, E. T. & Watson, G. N., A Course of Modern Analysis, 4th edn. Cambridge University Press, Cambridge, 1927, Section 7.32, p. 132.1 K., et al. Benchmark cases for outdoor sound propagation 7. Attenborough, models. Journal of the Acoustical Society of America, 97 (1) (1995) 173-191.

8. Abramowitz, M. & Stegun, I. A. (eds), Pocketbook Harri Deutsch,

Frankfurt,

APPENDIX

1984, Section

of Mathematical

Functions.

13.1.27.

1: ALTERNATIVE OF EQUATION

DERIVATION (17)

The recurrence relation (17) was not originally deduced from eqn (16). Instead, the polynomials W, were identified as multiples of confluent hypergeometric functions (, F,)

w,, =

,‘,‘“rlE!he”-’

,F,( 1 - n, 2 - 2n; h/e)

(AlI

l)...(cu+n-

(A31

where

and (a), =a(~+

1)

Now for y = 2~3, c):> 0 the functions ,Fl(a; 23 x) are closely related to Bessel functions, which, in turn, satisfy a three-term recurrence formula J/+I (x) + J/p I(x) = Wlx)J/W As the parameters in (Al) also satisfy y = 2a, a similar three-term to be expected for the polynomials W, of the form WN

=

Pn wn-

I +

qn

w,,-2

(A41 relation

(A3

is

Prediction of long range sound propagation

125

where pn and qn are linear and quadratic monomials in E and h, respectively. As we have throughout

p,,

and q,, cannot involve the first power of h, which suggests Pn = cl&;

(A71

qn = c;l12

where C, and C,’ are constants, to be determined by comparing C,_,,, with C,,i and CO,nwith CO,n-2, respectively. This leads to eqn (17); (16) merely served as a verification that (17) was valid in all its terms.

APPENDIX

2: PADF APPROXIMANTS

FOR THE CASE E=O

The use of confluent hypergeometric functions in the Pad& approximants of ex is well known but is explained here from first principles for completeness. The functions satisfy Kummer’s relation* IFI@; Y; x) =eX14(y-cr;

Y; -x)

W)

For positive integers L and M, and fractional r, this becomes ~FI(--L;

-L-M+r;

x)=eXIFl(-M+r;

-L-Mfr,

-x)

tA9)

where the series on the left reduces to a polynomial of degree L. As r tends to 0, the terms of the series on the right of degree (L + 1) to (L + M) inclusive also tend to 0, but resume for higher powers in x. Thus #1(-L;

-L - M; x)

1F,(-M;

-L - M; -x)

ex =

+

o(xL+“+‘)

=The series in (AlO) taken as polynomials agree with U,, and V,, (19). The first row in Table 1 corresponds to this special area.

(AlO) in eqn