The computation of ELF radio wave fields in the Earth-ionosphere duct

Journal of Atmospheric and Terrestrial Physics, Vol. 51, No. 3, pp. 233 239, 1989. Printed in Great Britain.

0021-9169/89 $3.00+ .00 Pergamon Press plc

The computation of ELF radio wave fields in the Earth-ionosphere duct D. LLANWYN JONES and G. S. JOYCE Physics Department, King's College, Strand, l,ondon WC2R 2LS (Received in final./orm 7 December 1988)

A~traet--Evaluation of ELF radio wave fields using the mode theory requires the computation of Legendre functions. Four different representations of the Legendre functions in terms of hypergeometric series have been programmed for computation on a microcomputer. By an appropriate choice of series for given values of the propagation constant and distance from the source, only a few (typically eight) terms are needed to calculate the Legendre functions with a precision of six decimal digits.

field. In (3) and (4) ~, M and h are the wave angular frequency, the current moment of the source and the The electromagnetic fields radiated by a dipole source effective height of the waveguide, respectively. The complex quantity v is the residue eigenvalue, located in the Earth-ionosphere duct may be represented in terms of a modal (or residue) series (WAIT, the calculation of which is not considered here. v is 1970; GALEJS, 1972). In the ELF (extremely low fre- simply related to the propagation constant of the 0quency; 3 Hz-3 kHz) band only the leading term of travelling waves in the Earth-ionosphere spherical this mode expansion is required in the calculation of shell. Essentially, the real part of v determines the the fields provided that the source-observer sep- phase velocity, and the imaginary part the attenuation aration is greater than some 500 km, i.e. several constant of the propagating fields, v embodies the waveguide heights. Here the 'waveguide' is the spheri- boundary conditions imposed on the fields by the Earth and the ionosphere and may be calculated cal dielectric shell bounded by the surface of the Earth numerically for any (isotropic or anisotropic) speciand the lower (D-E-region) of the ionosphere. Using geocentric spherical polar co-ordinates (r, 0, fied ionospheric model--see, for example, GALEJS ~b) we consider the field at the observer (receiver) point (1972), MtJRAKAMI (1976), TRAN (1980), JONES (a, 0, q~) radiated by a vertical electric dipole source (1984). Ivl is of order ka = 2~a/2 = f / 7 . 5 , where f is located at (a, 0, 0) where a is the Earth's radius and the wave frequency in hertz. In the ELF band Ivl varies between 0.5 and 400 and is thus a small quantity 0 is the angle subtended at the centre of the Earth for the lower part of this band (3 Hz-100 Hz), which by the source-receiver arc of length d = aO. The principal field components are the radial (ver- is currently of most interest. tical) electric field Er and the azimuthal (horizontal) The two Legendre functions P ° ( - c o s 0 ) and magnetic field H~. The leading term of the modal P~ ( - c o s 0) also appear in the field formulas for the three other primary radiation sources--the horizontal series defines these components as : electric and magnetic dipoles and the vertical magnetic Er = Eo P ° ( - cos 0)/(sin v~) (1) dipole (see, for example, GALEJS, 1972, Chapter 4). Assuming v to be known, the central problem in the H~ = H o e ~ ( - cos 0)/(sin w ) (2) evaluation of the field components from (1) and (2) in which : (or from the formulas for the other three primary sources) is that of computing the Legendre function Eo = [iMv(v + 1)]/(4a2eooh) (3) pro, with m = 0,1. At the higher ELF frequencies, and Ho = - M / ( 4 a h ) (4) in the VLF band, Ivl is large enough for the widely used asymptotic representation of the Legendre func(see WAIT, 1970; JONES, 1974a). In (1) and (2) P~ is tions to be highly accurate [JoN~s, 1970, equations the Legendre function of degree v and order m with (6) and (7)]. In the lower ELF band (3-100 Hz) an m = 0 for the electric field and m = 1 for the magnetic alternative method of computation is required. 233 l. I N T R O D U C T I O N

234

D. LLANWYNJONES and G. S. JOYCE 2. THE ZONAL HARMONIC SERIES---

3. THE HYPERGEOMETRIC SERIES

AN EARLY COMPUTATION SOLUTION

SCHUMANN 0952) pointed out that is is possible to expand the Legendre function P7 with m = 0 in terms of the Legendre polynomials po, where n = 0, 1, 2 . . . . , JONES and KEMP (1970) also considered the case m = l and gave a series expansion for pm with m = 0,1 in terms of the zonal harmonic functions P~. The detailed form of this expansion is (ERDf~LYIet al., 1953) : P , 7 ' ( - c o s 0 ) _ ( - l ) "+l ~ (2n+l)P~(cosO) sinvrc x ,=0 ~ n(n+ l ) - v ( v + 1)"

Many of the special functions of mathematical physics, including the Legendre functions, have representations in terms of hypergeometric series. The formulas used in the present work are given or easily derived from results in ERDI~LYI et al. (1953), and can also be obtained from AnRAMOWlTZand STEGUN (1968) o r MAGNUS, OBERHETTINGERand SONI 0966). The particular hypergeometric series appropriate for the Legendre functions is 2F~(a, b; c; z), in which a, b, c and z are real or complex parameters. This function is defined as (ERDI~LYIet al., 1953, p. 56)

(a),(b),, z" 2Fl(a'b;c;z)=,~_ o (c),, n!'

(5) This series has been applied widely for computing ELF fields in past years (see JONES, 1974b for a review of this early work). Equation (5) has the physical significance that it represents the field as a summation of the electromagnetic resonance modes of the Earthionosphere spherical shell cavity (known as 'Schumann resonances'), the condition n(n+ 1) = v(v+ 1) giving the (complex) frequencies of the resonator modes. For large n and 0 < 0 < ~z, the terms in the zonal harmonic series (5) are of order nm-3/z. It is clear, therefore, that this series will be slowly convergent for m = 0,1. For example, to compute the electric field some 100 terms are required for graphical accuracy and considerably more are needed for the magnetic field. NICKOLAENKOand RABINOWITZ(1974) showed that the convergence of equation (5) could be greatly increased by an ingenious 'partial extraction of the source singularity'. Their published formula has been corrected by CONNOR and MACKAY (1978) who have shown that six decimal digit accuracy can be achieved using a CDC 7600 main frame computer with 'up to 100 terms' in the expansion. Connor and Mackay present values of p o ( _ cos 0) for v = (0.1 + i0. l) and (3+i3) and for 0 in the range 20-160 °. Our current research programme involves the interpretation of experimentally acquired ELF data and requires the computation of electromagnetic field components in the frequency band 5 Hz-75 Hz by means of a microcomputer. Since this calculation is inside a program iteration loop it should be rapid. An alternative method of computing Legendre functions accurately, without restriction on v or 0 and using a machine with a lower precision than a main frame was thus sought. A solution to this requirement has been found in the hypergeometric series representation of the Legendre functions.

(6)

where (a)n = F(n + a)/F(a) and F(a) is the gamma function. The terms of the series (6) can be computed recursively (without any need to compute gamma functions) as follows. Since : F ( a + 1) = aF(a) it follows that : (a),+, = F ( n + a + 1)/F(a) = (n + a) F(n + a)/F(a) = (n + a) (a),. Hence, writing :

2Fl(a,b;c;z) = ~ (1,

(7)

,-0

it is seen that :

F(n+a)(n+b)z] d,+, = L ~ l ) ja,,, ,

n >i o.

(8)

From the relation (8) the nth term of the series can be evaluated by upward recursion with the starting value d o = 1. The hypergeometric series (6) is absolutely convergent within the unit cicle Izl < i and divergent for Izl > 1. For rapid convergence it is clear that [zl should be as small as possible. In practice, the number of terms required to evaluate the series within a specified error bound depends both on Izl and on the values of the other parameters. If the parameter values are too large the summation is affected by large rounding errors and few, if any, of the computed digits in the result will be correct. In the computational method to be developed here we have used four different series of the hypergeometric type to compute ELF radio fields without

235

The computation of ELF radio wave fields restriction on v or 0 and with a precision of 5-7 digits using a 24-bit mantissa microcomputer (RM Nimbus). The calculation has been programmed in (single precision) Fortran. Early tests of the first two series used were also made using a (32-bit mantissa) BBC-B machine with a BASIC interpreter. An estimate of the precision of the calculation is obtained by associating an error of + 1/2 LSB (least significant bit) with each term contributing to the sum and accumulating an error-sum through the computation. When the computation scheme detailed below is followed, only 1-12 terms are required to achieve the stated precision of 5-7 decimal digits in the result for any value of 0 or of v (or wave frequency). 4. C O M P U T A T I O N LEGENDRE

OF

values of Ivl the terms become large leading to a loss of significant digits in the summation. Series 1 thus becomes unsuitable for computation at some v (or frequency) dependent distance from the antipode of the source. We thus consider an alternative representation which is useful in the vicinity of the source. 4.2. The source series (series 2) A large number of linear transformation formulas for the hypergeometric function can be found from the analytic continuation of (6) through an integral representation of 2Fl. Of particular utility is the transformation z=cos2(O/2) to 1-z=sin2(O/2) (see ERDI~LYIet al., 1953, p. 110), the corresponding equation for the Legendre function pO being :

THE

(~/sinvrOP°(-cosO) = ~ D,Hn,

FUNCTIONS

(10)

n.-O

4.1. The antipode series (series 1) A fundamental representation of the Legendre function in terms of the hypergeometric series is (ERD~L¥1 et al., 1953, p. 148) P T ( - c o s 0 ) = Am(O)2F~[1+m+ v , m - v ;

where :

Hn = In (s)--2~k(n+ 1) +~O(l +n+v)+~O(n--v)

l+m;cos2(O/2)]

(9)

with :

A°(O) = 1 and A~(O) = -(1/2)v(v+ l) sin0. If we compare (7) and (9) it is readily seen that :

s = sin z (8/2) and ~O(x) denotes the digamma function. The digamma function is defined as the logarithmic derivative of the gamma function ~k(z) = (d/dz) In F(z).

~c

Pm(--cosO) = A,.~(0) ~ d,, n

It follows from the relation F ( z + 1) = zF(z) that

0

~,(z+ 1) = ~ ( z ) + 1/z.

where d, satisfies the recurrence relation :

~(n+m)(n+m+ l ) - v ( v +

]

×cos 2(8/2) d.,

(11)

Using (11) and the reflection formula for the digamma function : @ ( - z ) = rtcot (~z)+~O(1 + z )

n>~0

with the initial condition do = 1. We shall refer to this result as series 1. In the neighbourhood of the source antipode (8 = 180 ~) the series 1 converges rapidly because the convergence factor z = cos 2 (8/2) is small. At the antipode we have P°(1) = 1 and P ~ ( 1 ) = 0. However, near the source (8 = 0) the convergence factor approaches unity and series 1 becomes slowly convergent. For small values of ]vl and 0 some 2000 terms are needed to compute P T ( - c o s 8) and, as Ivl (corresponding to the wave frequency) is increased, a point will be reached where there is a severe loss of precision until, ultimately, none of the computed digits is correct. This problem arises because for n + m < Re(v) the terms of the series alternate in sign and for large

(12)

a recurrence relation for the digamma term in (10) can be developed leading to the following equations, listed in a form suitable for computation. (a) Start by setting Do = 1 and calculate

Ho=ln(s)+27+7~cot(ltv)+2~b(l+v )

(13)

in which ? = -@(1) = 0.577 215 6 5 5 . . . is the Euler constant. This calculation requires the computation of just one digamma function for each value of v and is easily programmed (see Appendix). (b) The terms in the summation can then be calculated for increasing values of n using : Dn+l = [n(n+ 1) --v(v+ 1)]sD,/(n+ 1) 2 and

236

D. LLANWYNJONES and G. S. JOYCE

H,+I = H , +

[(n + 1) + 2v(v + 1)] n~>0. ( n + 1)[n(n + 1) - v ( v + 1)]' (14)

p0 has been computed using this technique. To calculate P~ we use the fundamental relation P,! ( - c o s 0) = - (d/dO)P°(-cos 0). By differentiating (10), the following result is obtained : ac

[rt/sin (w)]P,! ( - cos 0) = -- cot (0/2) ~ D,(nH, + 1). n--

0

(15) The computation of P~ can thus be carried out along with that of po in the same program loop. Series 2 is highly convergent in the vicinity of the source. Very close to the source only the leading term is required and this term of (10) or (15) demonstrates the nature of the source singularity. It should be remembered, however, that the field equations (1) or (2) are only valid beyond a distance of some 500 km (0 = 5 °) since only one mode is implicit in these equations. In the vicinity of the antipode of the source (0 = 180 °) the order of 1000 terms are required to compute the fields and, for the larger values of Ivl (Ivl > 3), precision is reduced progressively so that at Ivl = 5 only 1 or 2 digit precision is obtained. If only series 1 and 2 were available, series 1 would be used for 0 >~ 90 ° and series 2 for 0 < 90 °. The maximum value of z [in (6)] would then be 0.5 on 0 = 90 °. With this scheme the lowest precision is obtained on 0 = 90 C where there is a noticeable reduction in the n u m b e r of significant digits for Ivl > 3. A further transformation convergent at the equator of the source (0 = 90 °) was therefore implemented. 4.3. The equatorial series (series 3) The transformation used changes z = cos 2 (0/2) in (9) to cos20 and yields an equation in which z = 0 on 0 = 9ff~. This enables precise calculations to be made near the equator of the source. The appropriate result is (ERDI~LYIet al., 1953, p. 144)

4.4. The asymptotic series (series 4) The so-called asymptotic series for the Legendre function is (ERDI~LYIe t al., 1953, p. 147) : /z 1/2F(v + 3/2) P~ ( - cos 0) = (2 cosec 0) 1/2F(v 4- m + 1) x ~ G, sin [ ( v + n + 1/2)(rt-0) n--0

+(m+ 1/2)Oz/2)+mz/2],

(17)

where F ( m + 1 / 2 1 , ( - m 4 - l/2),q[ (v+3/2),,n! j -(cosec0)/2]".

G,, = ]

The factor G, is identical to that which appears in the hypergeometric series. The formula (17) is actually the sum of two hypergeometric series with : z = - (1/2) cosec (0) exp [_+ i(O- 3~z/2)]

2 "n '/2sin"OPT(--cosO) = [C(u+ 1/2)C(v+ 1/2)1 ' x 2F,(u, v ; 1/2 ; cos 2 0) + 2 cos 0 [F(u)F(v)]- '

x 2F,(u+l/2, v4-1/2;3/2;cos20),

(16)

where :

u= -(v+m)/2,

using the recursion relation (8) with the appropriate values of a, b, c and z = cos 2 (0). To calculate the various gamma functions in (16) a short, general purpose routine to evaluate In F(z) with a precision of 7 digits was written (see Appendix). Implementation of (16) greatly increases the computation precision near the equator. It is apparent that in (16) v/2, rather than v appears in the argument list of the hypergeometric functions. This is of advantage as the problems with precision are exacerbated by large values of Ivl. Figure l shows the dependence on 0 of the trigonometric convergence factor z in the three series 1, 2 and 3. By using all three series the maximum value ofz is 1/4, reached at 0 = 60 ° and 0 = 120 °. An interim computation scheme was to choose series 2 for 0 < 60°; series 3 for 60 ° ~< 0 ~< 120 ° and series 1 for 120 ° < 0 ~< 180 ~'. By invoking all three series 5 7 digit precision can be obtained for all 0 for values of Ivl up to about 5 (corresponding to a wave frequency of 40 Hz). To extend the frequency range of the program the well-known asymptotic series has been included.

v=(l+v-m)/2.

Equation (16) has been programmed for application near the source's equator. In this computation the two hypergeometric functions were determined

and is a convergent series for 30 ° < 0 < 150 ° (corresponding to ]z] < 1). Although the convergence factor z for this series is generally greater than that for the other series discussed (see Fig. i) having a minm u m value of 0.5 on 0 = 90 °, a great advantage of series 4 is that v appears only in the denominator so the terms decrease in magnitude rapidly for large Iv]. In fact, the leading term of (17) (and only the leading term) is frequently used to compute fields in the upperELF and VLF frequency bands. Equation (17) is an asymptotic series for v--+ oo.

The computation of ELF radio wave fields

237

1,0

u ta-

0.5 z taJ ~z ta~ z

0 0

30

60

90 AN6LE

120

150

180

O°

Fig. 1. Variation of the convergence factor z with angle 0 for the four hypergeometric series.

This has the significance that if (17) is terminated at the N t h term, the remainder (error) is of the order of (iv+3/21) N-~ for the large Ivl. Hence precise results can be obtained from (17) even outside the range of convergence provided that 0 is not too near 0 ° or 180 ° . When calculating using (17) outside its range of convergence, the computation is terminated when either the terms start to increase or when the prescribed error bound is reached. In the former case the error of computation is estimated as the last term calculated. Series 4 has been included in our computation scheme, being invoked when Ivl > 5 and 0 is not too near 0 ° or 180 °. F o r Ivl > 5 numerical experiments have shown that (17) usually gives the greatest precision (in comparison with the other three series) when 0 is within the range 3 5 ° < 0 < 145 °. Outside this range we could compute P 7 using both (17) and either series 2 or series 1 (series 2, if near the source and series 1 if near the source antipode), the final result of the calculation being taken to be that with the lower predicted error. However, the need to calculate using two different series can be eliminated since numerical tests show that for the (complex) values of v appropriate to propagation in the Earth-ionosphere duct, in the region of the source series 2 gives greater precision than series 4 if (Re(v)+4.5)0 < 310 °. In the region of the antipode series 1 is more precise than series 4 if (Re(v) + 2) (180 - 0) < 320 °.

4.5. Summary of computation scheme In summary, the final computation scheme is as follows : (a) For I v l < 5 . If 0 < 6 0 °, use series 2; if 60 ° ~< 0 ~< 120 °, use series 3; i f 0 > 120 °, use series 1. (b) For Ivl/>5. If 0 < 3 5 °, use series 2 if (Re(v)+4.5)O<310 ° and series 4 otherwise; if 35 ° ~< 0 ~< 145 ° use series 4 ; if 0 > 145 °, use series 1 if (Re(v) + 2 ) ( 1 8 0 - 0) < 320 ° and series 4 otherwise. The error associated with the real and imaginary parts of each computed term of a series is assessed as 3 x 10 6 of the absolute value corresponding to + (1/2) LSB (least significant bit) on a 24-bit mantissa machine. Separate error sums are accumulated for the real and imaginary parts as the series are summed. As an estimate of the accuracy of calculation, the number of correct significant figures in the final result is calculated as log~0 (Isuml/error-sum). 5. EXAMPLES

OF

NUMERICAL

RESULTS

The validity of the computational technique has been tested by calculating the magnitude and argument of the Legendre function po for the values of v and 0 for which data were tabulated by CONNOR and MACKAY (1978). Our computations are in agreement with these earlier results. To check the calculation of P,! (for which no data were presented by Connor and Mackay) we have used our calculated values of po to evaluate PJ by means of the formula :

238

D. LLANWYN JONES and G. S. JOYCE

sin OP,! ( - - c o s 0) = v ( v + 1)[P°+ i ( - c o s 0) _ po

( _ cos O)]/(2v + 1).

A n e x a m p l e o f n u m e r i c a l d a t a o b t a i n e d using the f o r m a l i s m p r e s e n t e d in this p a p e r is given in T a b l e 1. D a t a are p r e s e n t e d for two ( c o m p l e x ) values o f v a p p r o p r i a t e to p r o p a g a t i o n in the E a r t h - i o n o s p h e r e d u c t for wave f r e q u e n c i e s o f 10 a n d 100 H z (JONES, 1970). T h e electric a n d m a g n e t i c fields o f a vertical

electric d i p o l e source o f m o m e n t M c a n be o b t a i n e d f r o m (1) to (4) o n c e the values o f t h e L e g e n d r e functions h a v e b e e n calculated. 6. CONCLUSION A m e t h o d o f c o m p u t i n g the L e g e n d r e f u n c t i o n s w h i c h arise in the t h e o r y o f p r o p a g a t i o n o f e x t r e m e l y low f r e q u e n c y ( E L F ) r a d i o waves in the E a r t h - i o n o sphere duct has been described. T h e t e c h n i q u e pre-

Table 1. Computed values of P T ' ( - cos 0) for m = 0,1. The values of v used for computation are appropriate to wave frequencies of I0 and 100 Hz. The angle 0 is expressed in degrees. SF is an estimate of the number of correct digits in the result, NT is the number of terms used in the calculation and NS is the series number (as given in the text) Re v = 1.342

Im v = -0.194

(Frequency = 10 Hz)

m

0/deg

Re P

R SF

Im P

1 SF

NT

NS

0 0 0 0 0 0 0 0 0

5.00 22.50 45.00 67.50 90.00 112.50 135.00 157.50 180.00

0.730764E+ 00 -0.299349E+00 -0.661169E+00 -0.619184E+00 -0.316389E+00 0.121763E+00 0.561914E+00 0.882885E+00 0.100000E+01

7.3 7.4 7.4 7.4 7.5 6.7 7.1 7.4

-0.871410E+00 -0.480867E + 00 -0.175793E+00 0.351350E-01 0.137970E+00 0.142504E+ 00 0.879068E-01 0.261014E-01 0.000000E + 00

7.5 7.3 6.7 6.7 7.5 7.3 7.4 7.5

3 5 8 8 2 8 8 5

2 2 2 3 3 3 1 l

1 1 1 1 1 1 I 1

5.00 22.50 45.00 67.50 90.00 112.50 135.00 157.50

7.5 7.5 6.9 7.0 7.5 7.4 7.4 7.5

3 6 9 8 2 8 8 5

2 2 2 3 3 3 l 1

180.00

-0.241914E+01 - 0.906193E+00 -0.662919E+00 -0.402531E+00 -0.125396E+00 0.848998E-01 0.170144E+00 0.125637E+00 0.000000E + 00

7.5 7.5 7.3 7.4 7.5 6.7 7.2 7.5

1

0.776751E+01 0.168122E+01 0.327878E+00 -0.489614E+00 -0.999426E +00 -0.117447E+01 -0.101468E+01 -0.583288E+00 0.000000E + 00

Re v = 15.640

Im v = -0.879

(Frequency = I00 Hz)

m

0/deg

Re P

R SF

Im P

I SF

NT

NS

0 0 0 0 0 0 0 0 0

5.00 22.50 45.00 67.50 90.00 112.50 135.00 157.50 180.00

-0.394273E+00 0.169390E+01 0.867022E+00 0.529946E+ 00 0.359496E+ 00 0.269687E +00 0.231312E+00 0.250810E+00 0.100000 E + 01

6.2 7.5 7.5 7.5 7.5 7.5 7.5 7.5

-0.476713E+01 -0.617592E+00 -0.375185E+00 -0.258399E+00 -0.188385E+00 -0.141128E+00 - 0 . I04996E+00 -0.702448E-01 0.000000E + 00

7.1 7.4 7.5 7.5 7.5 7.5 7.5 7.4

7 11 7 6 6 6 7 11

2 4 4 4 4 4 4 4

1 1 1 1 1 I 1 1

5.00 22.50 45.00 67.50 90.00 112.50 135.00 157.50

0.785080E+ 02 0.136771E+02 0.745182E+01 0.502261E+01 0.377136E+01 0.312587E+01 0.295065E+ 01 0.348378E+01

7.4 7.5 7.5 7.5 7.5 7.5 7.5 7.5

6.9 7.5 7.5 7.5 7.5 7.5 7.5 7.5

7 11 7 6 6 6 7 10

2 4 4 4 4 4 4 4

1

180.00

0.000000E + 00

-0.361260E+02 0.257073E+02 0.130509E+02 0.775047E+01 0.494698E+ 01 0.328003E +01 0.218132E+01 0.132596E+01 0.000000E + 00

239

The computation of E L F radio wave fields

sented here is capable of giving high precision using a much smaller number of terms than that needed by alternative techniques and is well suited to implementation by a small microcomputer.

Acknowledgements The work described here was undertaken in support of SERC research grant S G D 10351. We have also benefited from discussions with C. P. BURKE who independently programmed some of the series.

REFERENCES

ABRAMOWITZA. and STEGUN I, A.

1968

Handbook of Mathematical Functions. Dover Publi-

CONNOR J. N. L. and MACKAY D. C. ERDI~LY1A., MAGNUS W., OBERHETT1NGERF. and TRICOMI F. G. GALEJS J.

1978 1953

J. atmos, terr. Phys. 40, 977. Higher Transcendental Functions. McGraw-Hill, New

t 972

Terrestrial Propagation of Long Electromagnetic Waves.

JONES D. L. JONES D. L.

1970 1974a

Radio Sci. 5, 803. ELF-VLF Radio Wave Propagation, p. 207. D. Reidel,

JONES D. L. JONES D. L. and KEMP D. T. MAGNUS W., OBERHETTINGERF. and SONI R. P.

1974b 1970 1966

MURAKAMI Y. N1CKOLAENKOA. P. and RABINOWITZL. M. SCHUMANN W. O. TRAN A. WAIT J. R.

1976 1974 1952 1980 1970

IEEE Trans. Comms. 22, 477. J. atmos, terr. Phys. 32, 1095. Formulas and Theorems for the Special Functions of Mathematical Physics. Springer, New York. J. atmos, terr. Phys. 38, 371. J. atmos, terr. Phys. 36, 979. Z. Naturforsch. 7A, 149. J. atmos, terr. Phys. 42, 261. Electromagnetic Waves in Stratified Media (2nd Edn).

cations, New York.

York. Pergamon Press, New York.

Holland.

Pergamon Press, New York.

Reference is also made to the following unpublished material." JONES D . L .

1984

A P P E N D I X - - C O M P U T A T I O N OF GAMMA AND DIGAMMA FUNCTIONS

The g a m m a and d i g a m m a functions required for calculation of the Legendre functions have been obtained by using standard asymptotic expansions presented ~ by ABRAMOWXTZand STEGUN (1968). Specifically : In F(z) ~ (z - 1/2)" In (z) - z + (1/2)" In (2~) + 1/(12z) - I/(360z 3) + 1/(1260z 5) - 1/(1680z 7) + .-. ~b(z) ~ In (z) - l/(2z) - 1/(12z 2) + 1/(120z 4) - 1/(252z 6) + --.. In both these expansions z is an arbitrary complex quantity except that it m u s t not be real and negative (both functions are singular if z is zero or a negative integer). The number of terms given above are adequate to calculate the g a m m a or d i g a m m a functions to seven significant digits provided that Izl > I0 and not too near a negative integer. If Izl > 10 and the real part of z is positive, the g a m m a

The Interpretation of E L F Propagation Data using a Full-Wave Matrix Method. Paper presented at the XXIst URSI General Assembly (Available as Technical Report RP-8408).

and digamma functions are computed using these equations directly. If [z[ < l0 the two functions are computed by adding a suitable integer n to z to ensure that Iz+nl > 10, i.e. we compute In F ( z + n ) or ~(z+n). In F(z) and ~(z) can then be computed by repeated application of the recurrence formulas :

r ( ~ - 1) = r ( ~ ) / ( z -

l)

~b(z- 1) = ~k(z) - 1 / ( z - 1). Further, if the real part o f z is negative, we compute In F ( - z ) using the prescription above and then use the reflection formula, F(z) = - {~/[zF(-z)]} cosec gz, to find In F(z). The same trick can be employed to compute the d i g a m m a function when Re(z) is negative [using (12)] but for the calculations discussed here [Section 4.2] the argument of the digamma function has a positive real part.