- Email: [email protected]
Accurate evaluation of the Chapman function for atmospheric attenuation
Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
Accurate evaluation of the Chapman function for atmospheric attenuation David L. Huestis* Molecular Physics Laboratory, SRI International, Menlo Park, CA 94025, USA Received 10 March 2000; accepted 21 May 2000
Abstract We report development of the "rst practical procedures for arbitrarily accurate analytical evaluation of the Chapman function, Ch(X, ), which describes the attenuation of solar radiation by an exponential atmosphere with respect to altitude (z), solar zenith angle (), atmospheric scale height (H), and planetary radius (R), with X"(R#z)/H. A new asymptotic expansion in negative powers of X is derived with good accuracy that is nearly independent of the value of . An entirely new approach, with unlimited potential accuracy for all values of X and , is developed based on the newly discovered inhomogeneous Bessel's di!erential equation satis"ed by the Chapman function. We also provide e$cient procedures for accurate numerical integration and review the performance of formulas available in the literature. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Atmosphere; Solar; Radiation; Absorption; Column depth; Scale height; Chapman
1. Introduction In pioneering papers [1}4], Sidney Chapman investigated the altitude dependence of energy deposition by absorption of solar radiation including the apparent thickening of the atmosphere with the sun near the horizon. This work led to a sequence of subsequent publications [5}10] attempting to construct suitable analytical approximations to the `Chapman integrala or `Chapman functiona for an exponential atmosphere. Wilkes [11] provided a table of values evaluated by numerical integration, accurate to three places after the decimal point for ranging from 20 to 1003
* Corresponding author. Fax: #1-650-859-6196. E-mail address: [email protected] (D.L. Huestis). 0022-4073/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 0 ) 0 0 1 0 7 - 2
710
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
and X from 50 to 1000. Here is the solar zenith angle and X"(R#z)/H is a dimensionless ratio involving the planetary radius (R), altitude (z), and atmospheric scale height (H). Typical values of X for the terrestrial atmosphere range from 200 to 1300. The past work continues to be cited in more modern literature, [12,13], in spite of the fact that the scale height of the actual terrestrial atmosphere varies signi"cantly with altitude, due to variation in temperature, composition, and mechanisms of mixing, and thus the Chapman function for an exponential atmosphere is not quantitatively appropriate. Indeed, Pressly [14] used the "rst rocket determinations of the altitude dependence of the Earth's atmospheric density to compute the e!ective path-integrated air mass using the numerical Chapman Integral for a `reala (i.e., non-exponential) atmosphere. The work reported here was developed as part of testing of a numerical integration routine for arbitrary atmospheric pro"les. The idea was to use an exponential atmosphere as an accurately known reference case. We were somewhat surprised to discover the rather unsatisfactory computational performance of the formulas available in the literature. All the published analytical approaches are valid only asymptotically for large values of X. Some compact semiempirical expressions ([7], Eq. (3); [8], Eq. (52) corrected by [10]) are accurate to better than 0.1%, independent of , for large X. Others [6] have poor performance for certain ranges of . None of these is well-motivated mathematically or suitable for extension to higher numerical accuracy. The other analytic formulas correspond to asymptotic expansions in negative powers of X cos ([2], Eq. (22)) or X sin [5; 2, Eq. (38); 8, Eq. (42)], and thus are computationally useful only for small and large values of , respectively, and of uncertain potential for extension to higher accuracy, especially for intermediate values of . Here, we report development of the "rst practical procedures for arbitrarily accurate analytical evaluation of the Chapman function, Ch(X, ). A new asymptotic expansion in negative powers of X, for moderately large values of X ('36), yields good accuracy that is nearly independent of the value of . An entirely new approach, with unlimited potential accuracy for all values of X and , is based on the newly discovered inhomogeneous Bessel's di!erential equation satis"ed by the Chapman function. This is the "rst analytical method to achieve good accuracy for small X. We also provide e$cient procedures for accurate numerical integration. Source code is provided online at the web address http://www-mpl.sri.com/software/chapman/chapman.html.
2. Changing the variable of integration Following Rees [15, Section 2.2] we write the atmospheric column depth of a species of density n(z) from a point at an altitude of z to the sun at a solar zenith angle of as
N(, z )"
n(z) 1!
X
for 4/2 and
R#z \ sin dz, R#z
N(, z )"2N , z !N(!, z ) 2 Q
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
711
for '/2, where z "(R#z )sin !R, and R is the radius of the solid planet, e.g., Earth. For an exponential atmosphere, with n(z)"n(z )exp[!(z!z )/H], with a constant scale height, H, the column depth can be represented by the Chapman function, Ch(X, ) through the relation N(, z )"Hn(z ) Ch(X, ), where X"(R#z )/H. The Chapman function relates the actual column depth with that for a vertical sun. In what follows, we will usually limit our study to values of 4/2, because we can calculate Ch(X, ) for '/2 through the relation Ch(X, )"2Ch(X, /2)!Ch(X, !), if 5/2, where we already know [2, Eq. (16)] that Ch(X, /2)"Xe6K (X) with K being a well-known modi"ed Bessel function. We write the Chapman function as 1 R#z \ sin Ch(X, )" exp[!(z!z )/H] 1! dz. H R#z X If we make the substitutions y"(R#z)/(R#z )!1 and X"(R#z )/H we obtain sin \ Ch(X, )"X e\6W 1! dy. 1#y The further substitutions s"(1#y)/sin and s "1/sin give Ch(X, )"X sin e6 exp[!(X sin )s]s(s!1)\ ds. Q From this equation we can establish [16, Eqs. (9.6.23), (9.6.27)] that
(1)
(2)
(3)
Ch(X, /2)"Xe6K (X). If we now substitute s"1/sin we obtain
Q exp[X(1!sin csc )] csc d (4) which is the same as Chapman's original formulation [2, Eq. (10), p. 486]. If we "rst expand the square root in (2) Ch(X, )"X sin
1/((1!q)"1#q/((1!q)[1#((1!q)]
Note the sign change in the de"nition of K ( ).
712
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
with q"[sin /(1#y)]"sin, we obtain a new form of (4) without numerical singularities and simple convergence to Ch(0, )"Ch(X, 0)"1
Ch(X, )"1#Xsin
Q
exp[X(1!sin csc )]/(1#cos ) d. Alternatively, we may substitute t"y#t , with t "1!sin , into Eq. (2), "nding Ch(X, )"2X exp[!X(t!t )] f (t, ) dt, R where
(5)
(6)
f (t, )"(t#sin )/(t#2sin ) "(t!t #1)/(t!t #1#sin ) . Below we will use Eq. (6) to develop a compact asymptotic expansion of good accuracy for moderately large values of X ('36) and all values of , Eq. (5) to develop an e$cient numerical integral for Ch(X, ) for all values of X and , and Eq. (3) to derive an explicit analytical representation, valid for all values of X and , based on the di!erential equation satis"ed by Ch(X, ).
3. Asymptotic expansion At altitudes in the lower terrestrial thermosphere and below the values of X are rather large, in the range 200}1300. This means that only a limited range of t-values contribute signi"cantly to the integral (6) above, i.e. t!t (0.1. This suggests that we construct an asymptotic power series based on the expansion f (t,)" C ()(t!t )L L L that is nominally convergent only for (t!t )/(1#sin )(1. We "nd that C "(1#sin )\ , C "(1#2sin )/[2(1#sin ) ], C "!(1#4sin )/[8(1#sin ) ], C "(1#6sin )/[16(1#sin ) ]. Thus, we represent the Chapman function for 4/2 as Ch(X, )+2X C I , L L L where
I " L
exp[!X(t!t )](t!t )L dt. R
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
713
Evaluating these integrals in terms of the incomplete gamma functions [16, Section 6.5], we write
L n! I "1/(XL> ) exp(y ) (!y )H e\W yL\H dy L j!(n!j)! W H L n! (!y )H (n!j#1/2, y ) "1/(2XL> ) exp(y ) j!(n!j)! H with y "Xt "X(1!sin ). For small values of y , i.e. as P/2, we write [16, Eq. (6.5.17)] (1/2, y )"( erfc(y ) and calculate (k#1/2, y ) by upward recursion [16, Eqs. (6.5.3), (6.5.22)] (k#3/2, y )"(k#1/2)(k#1/2, y )#yI>exp(!y ). For large values of y we use the asymptotic expansion of the incomplete gamma function [16, Eq. (6.5.31)] to write I "[1/(2X y )][1!(1/2)/y #(3/4)/y !(15/8)/y #(105/16)/y !2] I "[1/(2X y )][1!1/y #(9/4)/y !(15/2)/y #(525/16)/y !2] I "[1/(2X y )][2!3/y #9/y !(75/2)/y #(1575/8)/y !2] I "[1/(2X y )][6!12/y #45/y !225/y #(11025/8)/y !2] For large values of X, now consider the special case of "sin "0. This gives y "(X, C "1, C "1/2, C "!1/8, and C "1/16. Using the above asymptotic expansions for I above, we "nd L Ch(X, 0)+2X C I L L L "1#O(X\). Similarly, the special case "/2, sin "1 gives y "0, C "1/(2, C "3/(4(2), C "!5/(32(2), C "7/(128(2), I "(/(2X ), I "(/(4X ), I "3(/(8X ), I "15(/(16X ), and Ch(X, /2)+2X C I L L L "((X/2)[1#3/(8X)!15/(128X)#105/(1024X)#O(X\)] which we recognize [2, Eq. (19), p. 487; 16, Eq. (9.7.2)] as Ch(X, /2)"X e6K (X)#O(X\ ). The above analysis shows that for su$ciently large values of X, for example X'36, the four-term asymptotic expansion developed here should have good relative accuracy for all values of . An interesting case is provided by considering only the "rst term in the asymptotic expansion: C "(1#sin )\ , I "(/4X) exp(Xt ) erfc((Xt ).
714
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
This results in the relatively simple expression Ch(X, )+[X/(1#sin )] exp[X(1!sin )] erfc[X(1!sin )] which is slightly more accurate than the semiempirical formula of Fitzmaurice [7, Eq. (3)] and slightly less accurate than that of Swider ([8, Eq. (52)]; corrected by [10]).
4. Numerical integration We are integrating Eq. (5) below:
Ch(X, )"1#Xsin
Q
exp[X(1!sin csc )]/(1#cos ) d. (5) The integrand is numerically very smooth, and rapidly varying only near "0. For X'0, we choose the lower limit of numerical integration such that the integrand is exponentially small, 7;10\ (3;10\ for double precision). The domain of integration is divided into 64 equal intervals (6000 for double precision), and integrated numerically using the 9-point closed Newton}Cotes formula from Hildebrand [17, Eq. (3.5.17), p. 75]. The results agree with the earlier numerical integrals of Wilkes [11] to within his stated accuracy, 0.001 absolute.
5. Inhomogeneous di4erential equation Motivated by the form of Eq. (3) above, consider the function de"ned by the integral
Z(t)"
e\RQ(s!1) ds
Q
Ch(X, )"!te6
d Z(t) dt
evaluated for t"Xsin and s "1/sin . Di!erentiating under the integral sign we have Z(t)"! e\RQs(s!1)\ ds Q and Z(t)" e\RQs(s!1)\ ds Q Integrating by parts we write Z as
Z(t)"exp(!ts )(s !1) !t e\RQ(s!1) ds. Q Similarly, we rewrite Z as
Z(t)"
Q
e\RQ(s!1)\ #(s!1) ds.
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
715
Combining the equations for Z and Z we obtain 1 1 Z(t)# Z(t)!Z(t)" exp(!ts )(s !1) . t t We recognize this as an inhomogeneous modi"ed Bessel's equation of order zero ([16], Eq. (9.6.1); [18], Eq. (5.62)). Following the prescription of `variation of parametersa given by Rabenstein [8, Section 1.12] we can write the general solution of this equation in the form
Z(t)"AI (t)#BK (t)!(s !1)
exp(!qs )[I (t)K (q)!K (t)I (q)] dq R with the coe$cients A and B to be determined by matching boundary conditions. Di!erentiating with respect to t, we obtain Ch(X, )"X sin e6!AI (X sin )#BK (Xsin ) #cos e\W[I (X sin )K (y sin )#K (Xsin )I (y sin )] dy . 6 Applying the boundary condition Ch(X, 0),1 requires that B"0. The requirement [2, Eq. (12), p. 486] that Ch(X, ) approach the "nite value of sec as XPR implies that A"0. Thus, we have the relatively compact equation
Ch(X, )"X sin cos e6
e\W[I (Xsin )K (y sin ) 6 #K (Xsin )I (y sin )] dy. (7) The integrand is positive everywhere, with an integrable logarithmic singularity in K (y sin ) as yP0 or P0, which is e!ectively suppressed because Xsin I (Xsin )PXsin for small X or small . A similarly weak singularity occurs as yPR if P/2. We examine
e\WI (ysin ) dy 6 which is always less than or equal to cos
e\WI (y sin ) dy which we can evaluate by substituting an integral representation [16, Eq. (9.6.16)] of I (y sin ) 1 exp(y sin cos) d I (y sin )" which gives [19, Eq. (858.524)] cos
1 e\WI (y sin ) dy" cos /(1!sin cos ) d"1 establishing that the above formula (7) for Ch(X, ) is stable for all values of X and , the "rst such analytical formulation. cos
716
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
The remaining task is thus to evaluate the integrals
6
e\W I (y sin ) dy
and
6
e\WK (y sin ) dy
for general values of X and . Given that accurate approximations for I (z) and K (z) are available [16, Eqs. (9.8.1}2), (9.8.5}6); 20; 20, pp. 178}179 (BESSI0.FOR and BESSK0.FOR)], we can evaluate the above by term-by-term integration of the approximating polynomials in ascending and descending powers of z and estimate the error in the analytical integration from the stated maximum approximation errors.
6. Performance FORTRAN code was written to evaluate Ch(X, x) through Eqs. (5), (6), and (7) above, as developed in previous sections. The source code is available on-line at the Web address http://www-mpl.sri.com/software/chapman/chapman.html. The methods are brie#y summarized below. Eq. (5): Numerical integral Eq. (6): Asymptotic expansion Eq. (7): Inhomogeneous di!erential equation The code was compiled, tested, and timed using Microsoft FORTRAN 5.1, Digital Equipment VAX FORTRAN, and IBM RS/6000 aix XL FORTRAN. Executed in MS-DOS mode under Windows 95 on a 160 MHz PC the following maximum relative errors were calculated for 044903,14X4800
Error
Method (comments)
6;10\ 1;10\ 6;10\ 1.5;10\ 7;10\
Eq. Eq. Eq. Eq. Eq.
(5), (5), (6), (6), (7)
real*4 real*8 (convergence test) X536 X'60
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
717
The following `typicala execution times, in microseconds, were observed: Time
Method (comments)
500 5000 25 120
Eq. Eq. Eq. Eq.
(5), real*4 (5), real*8 (6) (7)
The `recommendeda approach uses Eq. (6), asymptotic expansion, for X536, and Eq. (7), inhomogeneous di!erential equation, for X(36. No claims are made that the code is optimized for speed, accuracy, or compactness. The principal objectives were (1) (2) (3) (4)
Robustness with respect to argument values. Rigorous mathematical derivation and error control. Maximal use of `well-knowna mathematical functions. Ease of readability and mapping of theory to coding.
All internal formulas are calculated in double precision (real*8) and reduced to single precision (real*4) in the function return, as appropriate. The real*8 accuracy could be improved with more accurate representations of E1(), erfc(), I0(), I1(), K0(), and K1() than the polynomial approximations given by Abramowitz and Stegun [16] and Press et al. [20]. In the course of development, many additional representations and approximations of the Chapman function were attempted that failed to be robustly extendable to machine precision.
7. Comparison with formulas in the literature 7.1. Reanalysis of the work of Chapman The best analytical work in the literature is that of Chapman [2] himself, based on asymptotic expansions for large X. His small- expansion [2, Eq. (22)] of the form Ch(X, )"sec # b ()/XL L L where b ()"!tan sec, b ()"3tan sec, b ()"!(15 tan#12 tan) sec, b ()"(105 tan#60 tan) sec, b ()"!(945 tan#1260 tan#360 tan) sec
718
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
is quite accurate for small to moderate values of , and moderate to large values of X. On the other hand, his large- expansion performs poorly if implemented as originally written [2, Eq. (38)]: Ch(X, )"ye6K (y)erfc[((2y<)]#0.5 S (<)/(4y)K K K where y"X sin , <"[0.5(1/sin !1)] , S (<)" (2n!1)2(2n!2m#1)c <L\K\, K L LK> with coe$cients c taken from the expansion L (1#2<)/(1#<)\ " c <L. L L Although this expansion is nominally convergent for all <41, its convergence can be quite slow. Convergence of the series for S (<) is very poor for larger values of m. For m"3, a 4-term K polynomial (through c <) is usable only for <(0.4 ('503), while a 7-term sum (through c <) is only slightly better, requiring <(0.5 ('423). We have overcome these convergence di$culties by analytically summing the series for S (<), K using 4-term polynomial representations only for <(0.1 ('783). From the series that de"nes the coe$cients c we can write L S (<)"[(1#2<)/(1#<)\ !1]/< and derive expressions for S (<) for m'0 by analytical di!erentiation: K S (<)"[(d/d<)S (<)!c ]/<, S (<)"[(d/d<)S (<)!3c ]/<, S (<)"[(d/d<)S (<)!15c ]/<. With this formulation, for small <, the S (<) vanish as <, while for large <, S (<) approaches the K constant 2, and the higher S (<) vanish as 1/<. K 7.2. Relative accuracy versus Fig. 1 illustrates the relative accuracy for various approximations to the Chapman function, calculated versus for the "xed value of X"800. The exact value is taken from the double precision numerical integral, Eq. (5). The relative error is then expressed as relative error" approximation-exact /exact.
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
719
Fig. 1. Accuracy of Chapman function approximations vs. solar zenith Angle for X"800 (a) long-dash, Green et al., semiempirical rational fraction, [6, Eq. (11)], (b) long-dash, Fitzmaurice, semiempirical error function [7, Eq. (3)], (c) solid, present work, asymptotic, 1 term Eq. (6), (d) long-dash, Swider, semiempirical error function, [8, Eq. (52)]; corrected by [10], (e) short-dash, Chapman, large , 1 term ([2, Eq. (38)]; improved in present work), (f) solid, present work, inhomogeneous di!erential equation (Eq. (7)), (g) X points, present work, numerical (Eq. (5)), real*4, (h) short-dash, Chapman, large ([2, Eq. (38)]; improved in present work), (i) solid, present work, asymptotic (Eq. (6)), ( j) short-dash, Chapman, small [2, Eq. (22)].
The three short-dashed curves in Fig. 1 come from the work of Chapman as described in this section. The curve labeled ( j) comes from Chapman's small- approximation [2, Eq. (22)]. It is very accurate for small values of . The curve labeled (h) results from Chapman's large- approximation [2, Eq. (38)] as improved and implemented above [2, Eq. (22)]; improved in present work). Oscillatory behavior in this curve as well as in curve (i), discussed below, for large values of results from variations in accuracy of the polynomial approximations to the mathematical functions erfc() and K (). We see that very good accuracy is achieved even down to "13. The curve labeled (e) is from the "rst term in this asymptotic expansion (m"0 only from above). The three solid curves in Fig. 1 come from the present work. The curve labeled (i) comes from the asymptotic expansion based on Eq. (6). Curve (c) is from the "rst term. The curve labeled (f) comes from Eq. (7), which is based on the new inhomogeneous di!erential equation for the Chapman function. The ;-symbols, labeled (g) come from the single-precision numerical integral based on Eq. (5). The solid curves (present work) with the short-dashed curves (Chapman) show that good accuracy can be achieved, for large X, with either set of formulas. The three long-dashed curves in Fig. 1 come from oft cited analytical formulas from the literature. The curve labeled (a) comes from Green et al. [6, Eq. (11)]. Its poor performance was a surprise. Banks and Kockarts [12] attribute an accuracy of better than 2%, apparently based on the four values reported by Swider [8, p. 780] at X"50 and 1000 and "70 and 903. Our
720
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
investigation of this formula showed large errors for between 80 and 893 for all values of X. For example, at X"1000 and "863 we "nd a relative error of 16%. The other two long-dashed curves in Fig. 1 come from the empirical formulas of Fitzmaurice [7, Eq. (3)] and Swider ([8, Eq. (52)]; corrected by [10]), also cited by Banks and Kockarts [12]. Both of these formulas are comparable in accuracy to the "rst terms in the asymptotic expansions developed here and previously by Chapman (large-), curves (c) and (e), respectively.
7.3. Relative accuracy versus X Fig. 2 illustrates the relative worst-case accuracy for various approximations to the Chapman function. For each value of X, the exact and approximate values were calculated for values of ranging from 0 to 903, with the largest relative error reported. For all the curves shown, the relative error continues to decrease as X is increased above the largest value shown. As before, the short-dashed curve, labeled (b), comes from the (improved) Chapman formulas, with the small- formula [2, Eq. (22)] or large- formula ([2, Eq. (38)]; improved in present work) chosen based on which has the smallest `last terma. The solid curves come from the present work, with curve (a) from the Eq. (7), inhomogeneous di!erential equation, curve (c) from Eq. (6), asymptotic expansion, and curve (d) from Eq. (5), single precision numerical integration. By choosing Eq. (7) for X(36 and Eq. (6) for X536, we clearly have an accurate representation of Ch(X, ) for all values of X and .
Fig. 2. Worst-case accuracy of Chapman function approximations vs. X (a) present work, inhomogeneous di!erential equation (Eq. (7)), (b) Chapman, small- and large- combined [2, Eq. (22)]; [8, Eq. (52)]; corrected by [10], (c) present work, asymptotic (Eq. (6)), (d) present work, numerical (Eq. (5) real*4).
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
721
8. Conclusions We have developed the "rst analytical formulations for the Chapman function, Ch(X, ), for atmospheric attenuation in an exponential atmosphere, that are extendable to arbitrary accuracy for all values of the reduced radius of curvature, X"(R#z)/H, and the solar zenith angle, . The newly discovered inhomogeneous di!erential equation satis"ed by Ch(X, ) provides the basis for accurate calculation for small values of X. For any "xed value of X, the present formulas provide a better worst-case accuracy versus than any available in the literature. Acknowledgements This work was supported by grants from the NASA Sun}Earth Connection and NSF Aeronomy programs.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
Chapman S. Proc Phys Soc (London) 1931;43:26. Chapman S. Proc Phys Soc (London) 1931;43:483. Chapman S. Proc Phys Soc (London) 1939;51:93. Chapman S. Proc Phys Soc B (London) 1953;66:510. Hulburt EO. Phys Rev 1939;55:639. Green AES, Lindemeyer CS, Griggs M. J Geophys Res 1964;69:493. Fitzmaurice JA. Appl Opt 1964;3:640. Swider Jr. W. Planet Space Sci 1964;12:761. Swider W Jr. Gardner ME. On the accuracy of certain approximations for the Chapman function, Environmental Research Papers No. 272, AFCRL-67-0468 Air Force Cambridge Research Laboratory, Bedford, MA, 1967. Swider Jr. W, Gardner ME. Appl Opt 1969;8:725. Wilkes MV. Proc Phys Soc B (London) 1954;67:304. Banks PM, Kockarts G. Aeronomy. Part B. New York: Academic, 1973, p. 107. Brasseur G, Solomon S. Aeronomy of the middle atmosphere. 2nd ed. Dortrecht: Reidel, 1986, p. 106. Pressly EC. Phys Rev 1953;89:654. Rees MH. Physics and chemistry of the upper atmosphere. Cambridge: Cambridge University Press, 1989. Abramowitz M, Stegun IA. Handbook of mathematical functions, Vol. NBS AMS 55. Washington, DC: USGPO, 1964. Hildebrand FB. Introduction to numerical analysis. New York: McGraw-Hill, 1956. Rabenstein AL. Introduction to ordinary di!erential equations. New York: Academic Press, 1966. Dwight HB. Tables of integrals and other mathematical data, 4th ed. New York: Macmillan, 1961. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes: the art of scienti"c programming. London: Cambridge University Press, 1986.
Accurate evaluation of the Chapman function for atmospheric attenuation David L. Huestis* Molecular Physics Laboratory, SRI International, Menlo Park, CA 94025, USA Received 10 March 2000; accepted 21 May 2000
Abstract We report development of the "rst practical procedures for arbitrarily accurate analytical evaluation of the Chapman function, Ch(X, ), which describes the attenuation of solar radiation by an exponential atmosphere with respect to altitude (z), solar zenith angle (), atmospheric scale height (H), and planetary radius (R), with X"(R#z)/H. A new asymptotic expansion in negative powers of X is derived with good accuracy that is nearly independent of the value of . An entirely new approach, with unlimited potential accuracy for all values of X and , is developed based on the newly discovered inhomogeneous Bessel's di!erential equation satis"ed by the Chapman function. We also provide e$cient procedures for accurate numerical integration and review the performance of formulas available in the literature. 2001 Elsevier Science Ltd. All rights reserved. Keywords: Atmosphere; Solar; Radiation; Absorption; Column depth; Scale height; Chapman
1. Introduction In pioneering papers [1}4], Sidney Chapman investigated the altitude dependence of energy deposition by absorption of solar radiation including the apparent thickening of the atmosphere with the sun near the horizon. This work led to a sequence of subsequent publications [5}10] attempting to construct suitable analytical approximations to the `Chapman integrala or `Chapman functiona for an exponential atmosphere. Wilkes [11] provided a table of values evaluated by numerical integration, accurate to three places after the decimal point for ranging from 20 to 1003
* Corresponding author. Fax: #1-650-859-6196. E-mail address: [email protected] (D.L. Huestis). 0022-4073/01/$ - see front matter 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 0 ) 0 0 1 0 7 - 2
710
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
and X from 50 to 1000. Here is the solar zenith angle and X"(R#z)/H is a dimensionless ratio involving the planetary radius (R), altitude (z), and atmospheric scale height (H). Typical values of X for the terrestrial atmosphere range from 200 to 1300. The past work continues to be cited in more modern literature, [12,13], in spite of the fact that the scale height of the actual terrestrial atmosphere varies signi"cantly with altitude, due to variation in temperature, composition, and mechanisms of mixing, and thus the Chapman function for an exponential atmosphere is not quantitatively appropriate. Indeed, Pressly [14] used the "rst rocket determinations of the altitude dependence of the Earth's atmospheric density to compute the e!ective path-integrated air mass using the numerical Chapman Integral for a `reala (i.e., non-exponential) atmosphere. The work reported here was developed as part of testing of a numerical integration routine for arbitrary atmospheric pro"les. The idea was to use an exponential atmosphere as an accurately known reference case. We were somewhat surprised to discover the rather unsatisfactory computational performance of the formulas available in the literature. All the published analytical approaches are valid only asymptotically for large values of X. Some compact semiempirical expressions ([7], Eq. (3); [8], Eq. (52) corrected by [10]) are accurate to better than 0.1%, independent of , for large X. Others [6] have poor performance for certain ranges of . None of these is well-motivated mathematically or suitable for extension to higher numerical accuracy. The other analytic formulas correspond to asymptotic expansions in negative powers of X cos ([2], Eq. (22)) or X sin [5; 2, Eq. (38); 8, Eq. (42)], and thus are computationally useful only for small and large values of , respectively, and of uncertain potential for extension to higher accuracy, especially for intermediate values of . Here, we report development of the "rst practical procedures for arbitrarily accurate analytical evaluation of the Chapman function, Ch(X, ). A new asymptotic expansion in negative powers of X, for moderately large values of X ('36), yields good accuracy that is nearly independent of the value of . An entirely new approach, with unlimited potential accuracy for all values of X and , is based on the newly discovered inhomogeneous Bessel's di!erential equation satis"ed by the Chapman function. This is the "rst analytical method to achieve good accuracy for small X. We also provide e$cient procedures for accurate numerical integration. Source code is provided online at the web address http://www-mpl.sri.com/software/chapman/chapman.html.
2. Changing the variable of integration Following Rees [15, Section 2.2] we write the atmospheric column depth of a species of density n(z) from a point at an altitude of z to the sun at a solar zenith angle of as
N(, z )"
n(z) 1!
X
for 4/2 and
R#z \ sin dz, R#z
N(, z )"2N , z !N(!, z ) 2 Q
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
711
for '/2, where z "(R#z )sin !R, and R is the radius of the solid planet, e.g., Earth. For an exponential atmosphere, with n(z)"n(z )exp[!(z!z )/H], with a constant scale height, H, the column depth can be represented by the Chapman function, Ch(X, ) through the relation N(, z )"Hn(z ) Ch(X, ), where X"(R#z )/H. The Chapman function relates the actual column depth with that for a vertical sun. In what follows, we will usually limit our study to values of 4/2, because we can calculate Ch(X, ) for '/2 through the relation Ch(X, )"2Ch(X, /2)!Ch(X, !), if 5/2, where we already know [2, Eq. (16)] that Ch(X, /2)"Xe6K (X) with K being a well-known modi"ed Bessel function. We write the Chapman function as 1 R#z \ sin Ch(X, )" exp[!(z!z )/H] 1! dz. H R#z X If we make the substitutions y"(R#z)/(R#z )!1 and X"(R#z )/H we obtain sin \ Ch(X, )"X e\6W 1! dy. 1#y The further substitutions s"(1#y)/sin and s "1/sin give Ch(X, )"X sin e6 exp[!(X sin )s]s(s!1)\ ds. Q From this equation we can establish [16, Eqs. (9.6.23), (9.6.27)] that
(1)
(2)
(3)
Ch(X, /2)"Xe6K (X). If we now substitute s"1/sin we obtain
Q exp[X(1!sin csc )] csc d (4) which is the same as Chapman's original formulation [2, Eq. (10), p. 486]. If we "rst expand the square root in (2) Ch(X, )"X sin
1/((1!q)"1#q/((1!q)[1#((1!q)]
Note the sign change in the de"nition of K ( ).
712
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
with q"[sin /(1#y)]"sin, we obtain a new form of (4) without numerical singularities and simple convergence to Ch(0, )"Ch(X, 0)"1
Ch(X, )"1#Xsin
Q
exp[X(1!sin csc )]/(1#cos ) d. Alternatively, we may substitute t"y#t , with t "1!sin , into Eq. (2), "nding Ch(X, )"2X exp[!X(t!t )] f (t, ) dt, R where
(5)
(6)
f (t, )"(t#sin )/(t#2sin ) "(t!t #1)/(t!t #1#sin ) . Below we will use Eq. (6) to develop a compact asymptotic expansion of good accuracy for moderately large values of X ('36) and all values of , Eq. (5) to develop an e$cient numerical integral for Ch(X, ) for all values of X and , and Eq. (3) to derive an explicit analytical representation, valid for all values of X and , based on the di!erential equation satis"ed by Ch(X, ).
3. Asymptotic expansion At altitudes in the lower terrestrial thermosphere and below the values of X are rather large, in the range 200}1300. This means that only a limited range of t-values contribute signi"cantly to the integral (6) above, i.e. t!t (0.1. This suggests that we construct an asymptotic power series based on the expansion f (t,)" C ()(t!t )L L L that is nominally convergent only for (t!t )/(1#sin )(1. We "nd that C "(1#sin )\ , C "(1#2sin )/[2(1#sin ) ], C "!(1#4sin )/[8(1#sin ) ], C "(1#6sin )/[16(1#sin ) ]. Thus, we represent the Chapman function for 4/2 as Ch(X, )+2X C I , L L L where
I " L
exp[!X(t!t )](t!t )L dt. R
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
713
Evaluating these integrals in terms of the incomplete gamma functions [16, Section 6.5], we write
L n! I "1/(XL> ) exp(y ) (!y )H e\W yL\H dy L j!(n!j)! W H L n! (!y )H (n!j#1/2, y ) "1/(2XL> ) exp(y ) j!(n!j)! H with y "Xt "X(1!sin ). For small values of y , i.e. as P/2, we write [16, Eq. (6.5.17)] (1/2, y )"( erfc(y ) and calculate (k#1/2, y ) by upward recursion [16, Eqs. (6.5.3), (6.5.22)] (k#3/2, y )"(k#1/2)(k#1/2, y )#yI>exp(!y ). For large values of y we use the asymptotic expansion of the incomplete gamma function [16, Eq. (6.5.31)] to write I "[1/(2X y )][1!(1/2)/y #(3/4)/y !(15/8)/y #(105/16)/y !2] I "[1/(2X y )][1!1/y #(9/4)/y !(15/2)/y #(525/16)/y !2] I "[1/(2X y )][2!3/y #9/y !(75/2)/y #(1575/8)/y !2] I "[1/(2X y )][6!12/y #45/y !225/y #(11025/8)/y !2] For large values of X, now consider the special case of "sin "0. This gives y "(X, C "1, C "1/2, C "!1/8, and C "1/16. Using the above asymptotic expansions for I above, we "nd L Ch(X, 0)+2X C I L L L "1#O(X\). Similarly, the special case "/2, sin "1 gives y "0, C "1/(2, C "3/(4(2), C "!5/(32(2), C "7/(128(2), I "(/(2X ), I "(/(4X ), I "3(/(8X ), I "15(/(16X ), and Ch(X, /2)+2X C I L L L "((X/2)[1#3/(8X)!15/(128X)#105/(1024X)#O(X\)] which we recognize [2, Eq. (19), p. 487; 16, Eq. (9.7.2)] as Ch(X, /2)"X e6K (X)#O(X\ ). The above analysis shows that for su$ciently large values of X, for example X'36, the four-term asymptotic expansion developed here should have good relative accuracy for all values of . An interesting case is provided by considering only the "rst term in the asymptotic expansion: C "(1#sin )\ , I "(/4X) exp(Xt ) erfc((Xt ).
714
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
This results in the relatively simple expression Ch(X, )+[X/(1#sin )] exp[X(1!sin )] erfc[X(1!sin )] which is slightly more accurate than the semiempirical formula of Fitzmaurice [7, Eq. (3)] and slightly less accurate than that of Swider ([8, Eq. (52)]; corrected by [10]).
4. Numerical integration We are integrating Eq. (5) below:
Ch(X, )"1#Xsin
Q
exp[X(1!sin csc )]/(1#cos ) d. (5) The integrand is numerically very smooth, and rapidly varying only near "0. For X'0, we choose the lower limit of numerical integration such that the integrand is exponentially small, 7;10\ (3;10\ for double precision). The domain of integration is divided into 64 equal intervals (6000 for double precision), and integrated numerically using the 9-point closed Newton}Cotes formula from Hildebrand [17, Eq. (3.5.17), p. 75]. The results agree with the earlier numerical integrals of Wilkes [11] to within his stated accuracy, 0.001 absolute.
5. Inhomogeneous di4erential equation Motivated by the form of Eq. (3) above, consider the function de"ned by the integral
Z(t)"
e\RQ(s!1) ds
Q
Ch(X, )"!te6
d Z(t) dt
evaluated for t"Xsin and s "1/sin . Di!erentiating under the integral sign we have Z(t)"! e\RQs(s!1)\ ds Q and Z(t)" e\RQs(s!1)\ ds Q Integrating by parts we write Z as
Z(t)"exp(!ts )(s !1) !t e\RQ(s!1) ds. Q Similarly, we rewrite Z as
Z(t)"
Q
e\RQ(s!1)\ #(s!1) ds.
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
715
Combining the equations for Z and Z we obtain 1 1 Z(t)# Z(t)!Z(t)" exp(!ts )(s !1) . t t We recognize this as an inhomogeneous modi"ed Bessel's equation of order zero ([16], Eq. (9.6.1); [18], Eq. (5.62)). Following the prescription of `variation of parametersa given by Rabenstein [8, Section 1.12] we can write the general solution of this equation in the form
Z(t)"AI (t)#BK (t)!(s !1)
exp(!qs )[I (t)K (q)!K (t)I (q)] dq R with the coe$cients A and B to be determined by matching boundary conditions. Di!erentiating with respect to t, we obtain Ch(X, )"X sin e6!AI (X sin )#BK (Xsin ) #cos e\W[I (X sin )K (y sin )#K (Xsin )I (y sin )] dy . 6 Applying the boundary condition Ch(X, 0),1 requires that B"0. The requirement [2, Eq. (12), p. 486] that Ch(X, ) approach the "nite value of sec as XPR implies that A"0. Thus, we have the relatively compact equation
Ch(X, )"X sin cos e6
e\W[I (Xsin )K (y sin ) 6 #K (Xsin )I (y sin )] dy. (7) The integrand is positive everywhere, with an integrable logarithmic singularity in K (y sin ) as yP0 or P0, which is e!ectively suppressed because Xsin I (Xsin )PXsin for small X or small . A similarly weak singularity occurs as yPR if P/2. We examine
e\WI (ysin ) dy 6 which is always less than or equal to cos
e\WI (y sin ) dy which we can evaluate by substituting an integral representation [16, Eq. (9.6.16)] of I (y sin ) 1 exp(y sin cos) d I (y sin )" which gives [19, Eq. (858.524)] cos
1 e\WI (y sin ) dy" cos /(1!sin cos ) d"1 establishing that the above formula (7) for Ch(X, ) is stable for all values of X and , the "rst such analytical formulation. cos
716
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
The remaining task is thus to evaluate the integrals
6
e\W I (y sin ) dy
and
6
e\WK (y sin ) dy
for general values of X and . Given that accurate approximations for I (z) and K (z) are available [16, Eqs. (9.8.1}2), (9.8.5}6); 20; 20, pp. 178}179 (BESSI0.FOR and BESSK0.FOR)], we can evaluate the above by term-by-term integration of the approximating polynomials in ascending and descending powers of z and estimate the error in the analytical integration from the stated maximum approximation errors.
6. Performance FORTRAN code was written to evaluate Ch(X, x) through Eqs. (5), (6), and (7) above, as developed in previous sections. The source code is available on-line at the Web address http://www-mpl.sri.com/software/chapman/chapman.html. The methods are brie#y summarized below. Eq. (5): Numerical integral Eq. (6): Asymptotic expansion Eq. (7): Inhomogeneous di!erential equation The code was compiled, tested, and timed using Microsoft FORTRAN 5.1, Digital Equipment VAX FORTRAN, and IBM RS/6000 aix XL FORTRAN. Executed in MS-DOS mode under Windows 95 on a 160 MHz PC the following maximum relative errors were calculated for 044903,14X4800
Error
Method (comments)
6;10\ 1;10\ 6;10\ 1.5;10\ 7;10\
Eq. Eq. Eq. Eq. Eq.
(5), (5), (6), (6), (7)
real*4 real*8 (convergence test) X536 X'60
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
717
The following `typicala execution times, in microseconds, were observed: Time
Method (comments)
500 5000 25 120
Eq. Eq. Eq. Eq.
(5), real*4 (5), real*8 (6) (7)
The `recommendeda approach uses Eq. (6), asymptotic expansion, for X536, and Eq. (7), inhomogeneous di!erential equation, for X(36. No claims are made that the code is optimized for speed, accuracy, or compactness. The principal objectives were (1) (2) (3) (4)
Robustness with respect to argument values. Rigorous mathematical derivation and error control. Maximal use of `well-knowna mathematical functions. Ease of readability and mapping of theory to coding.
All internal formulas are calculated in double precision (real*8) and reduced to single precision (real*4) in the function return, as appropriate. The real*8 accuracy could be improved with more accurate representations of E1(), erfc(), I0(), I1(), K0(), and K1() than the polynomial approximations given by Abramowitz and Stegun [16] and Press et al. [20]. In the course of development, many additional representations and approximations of the Chapman function were attempted that failed to be robustly extendable to machine precision.
7. Comparison with formulas in the literature 7.1. Reanalysis of the work of Chapman The best analytical work in the literature is that of Chapman [2] himself, based on asymptotic expansions for large X. His small- expansion [2, Eq. (22)] of the form Ch(X, )"sec # b ()/XL L L where b ()"!tan sec, b ()"3tan sec, b ()"!(15 tan#12 tan) sec, b ()"(105 tan#60 tan) sec, b ()"!(945 tan#1260 tan#360 tan) sec
718
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
is quite accurate for small to moderate values of , and moderate to large values of X. On the other hand, his large- expansion performs poorly if implemented as originally written [2, Eq. (38)]: Ch(X, )"ye6K (y)erfc[((2y<)]#0.5 S (<)/(4y)K K K where y"X sin , <"[0.5(1/sin !1)] , S (<)" (2n!1)2(2n!2m#1)c <L\K\, K L LK> with coe$cients c taken from the expansion L (1#2<)/(1#<)\ " c <L. L L Although this expansion is nominally convergent for all <41, its convergence can be quite slow. Convergence of the series for S (<) is very poor for larger values of m. For m"3, a 4-term K polynomial (through c <) is usable only for <(0.4 ('503), while a 7-term sum (through c <) is only slightly better, requiring <(0.5 ('423). We have overcome these convergence di$culties by analytically summing the series for S (<), K using 4-term polynomial representations only for <(0.1 ('783). From the series that de"nes the coe$cients c we can write L S (<)"[(1#2<)/(1#<)\ !1]/< and derive expressions for S (<) for m'0 by analytical di!erentiation: K S (<)"[(d/d<)S (<)!c ]/<, S (<)"[(d/d<)S (<)!3c ]/<, S (<)"[(d/d<)S (<)!15c ]/<. With this formulation, for small <, the S (<) vanish as <, while for large <, S (<) approaches the K constant 2, and the higher S (<) vanish as 1/<. K 7.2. Relative accuracy versus Fig. 1 illustrates the relative accuracy for various approximations to the Chapman function, calculated versus for the "xed value of X"800. The exact value is taken from the double precision numerical integral, Eq. (5). The relative error is then expressed as relative error" approximation-exact /exact.
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
719
Fig. 1. Accuracy of Chapman function approximations vs. solar zenith Angle for X"800 (a) long-dash, Green et al., semiempirical rational fraction, [6, Eq. (11)], (b) long-dash, Fitzmaurice, semiempirical error function [7, Eq. (3)], (c) solid, present work, asymptotic, 1 term Eq. (6), (d) long-dash, Swider, semiempirical error function, [8, Eq. (52)]; corrected by [10], (e) short-dash, Chapman, large , 1 term ([2, Eq. (38)]; improved in present work), (f) solid, present work, inhomogeneous di!erential equation (Eq. (7)), (g) X points, present work, numerical (Eq. (5)), real*4, (h) short-dash, Chapman, large ([2, Eq. (38)]; improved in present work), (i) solid, present work, asymptotic (Eq. (6)), ( j) short-dash, Chapman, small [2, Eq. (22)].
The three short-dashed curves in Fig. 1 come from the work of Chapman as described in this section. The curve labeled ( j) comes from Chapman's small- approximation [2, Eq. (22)]. It is very accurate for small values of . The curve labeled (h) results from Chapman's large- approximation [2, Eq. (38)] as improved and implemented above [2, Eq. (22)]; improved in present work). Oscillatory behavior in this curve as well as in curve (i), discussed below, for large values of results from variations in accuracy of the polynomial approximations to the mathematical functions erfc() and K (). We see that very good accuracy is achieved even down to "13. The curve labeled (e) is from the "rst term in this asymptotic expansion (m"0 only from above). The three solid curves in Fig. 1 come from the present work. The curve labeled (i) comes from the asymptotic expansion based on Eq. (6). Curve (c) is from the "rst term. The curve labeled (f) comes from Eq. (7), which is based on the new inhomogeneous di!erential equation for the Chapman function. The ;-symbols, labeled (g) come from the single-precision numerical integral based on Eq. (5). The solid curves (present work) with the short-dashed curves (Chapman) show that good accuracy can be achieved, for large X, with either set of formulas. The three long-dashed curves in Fig. 1 come from oft cited analytical formulas from the literature. The curve labeled (a) comes from Green et al. [6, Eq. (11)]. Its poor performance was a surprise. Banks and Kockarts [12] attribute an accuracy of better than 2%, apparently based on the four values reported by Swider [8, p. 780] at X"50 and 1000 and "70 and 903. Our
720
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
investigation of this formula showed large errors for between 80 and 893 for all values of X. For example, at X"1000 and "863 we "nd a relative error of 16%. The other two long-dashed curves in Fig. 1 come from the empirical formulas of Fitzmaurice [7, Eq. (3)] and Swider ([8, Eq. (52)]; corrected by [10]), also cited by Banks and Kockarts [12]. Both of these formulas are comparable in accuracy to the "rst terms in the asymptotic expansions developed here and previously by Chapman (large-), curves (c) and (e), respectively.
7.3. Relative accuracy versus X Fig. 2 illustrates the relative worst-case accuracy for various approximations to the Chapman function. For each value of X, the exact and approximate values were calculated for values of ranging from 0 to 903, with the largest relative error reported. For all the curves shown, the relative error continues to decrease as X is increased above the largest value shown. As before, the short-dashed curve, labeled (b), comes from the (improved) Chapman formulas, with the small- formula [2, Eq. (22)] or large- formula ([2, Eq. (38)]; improved in present work) chosen based on which has the smallest `last terma. The solid curves come from the present work, with curve (a) from the Eq. (7), inhomogeneous di!erential equation, curve (c) from Eq. (6), asymptotic expansion, and curve (d) from Eq. (5), single precision numerical integration. By choosing Eq. (7) for X(36 and Eq. (6) for X536, we clearly have an accurate representation of Ch(X, ) for all values of X and .
Fig. 2. Worst-case accuracy of Chapman function approximations vs. X (a) present work, inhomogeneous di!erential equation (Eq. (7)), (b) Chapman, small- and large- combined [2, Eq. (22)]; [8, Eq. (52)]; corrected by [10], (c) present work, asymptotic (Eq. (6)), (d) present work, numerical (Eq. (5) real*4).
D.L. Huestis / Journal of Quantitative Spectroscopy & Radiative Transfer 69 (2001) 709}721
721
8. Conclusions We have developed the "rst analytical formulations for the Chapman function, Ch(X, ), for atmospheric attenuation in an exponential atmosphere, that are extendable to arbitrary accuracy for all values of the reduced radius of curvature, X"(R#z)/H, and the solar zenith angle, . The newly discovered inhomogeneous di!erential equation satis"ed by Ch(X, ) provides the basis for accurate calculation for small values of X. For any "xed value of X, the present formulas provide a better worst-case accuracy versus than any available in the literature. Acknowledgements This work was supported by grants from the NASA Sun}Earth Connection and NSF Aeronomy programs.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]
Chapman S. Proc Phys Soc (London) 1931;43:26. Chapman S. Proc Phys Soc (London) 1931;43:483. Chapman S. Proc Phys Soc (London) 1939;51:93. Chapman S. Proc Phys Soc B (London) 1953;66:510. Hulburt EO. Phys Rev 1939;55:639. Green AES, Lindemeyer CS, Griggs M. J Geophys Res 1964;69:493. Fitzmaurice JA. Appl Opt 1964;3:640. Swider Jr. W. Planet Space Sci 1964;12:761. Swider W Jr. Gardner ME. On the accuracy of certain approximations for the Chapman function, Environmental Research Papers No. 272, AFCRL-67-0468 Air Force Cambridge Research Laboratory, Bedford, MA, 1967. Swider Jr. W, Gardner ME. Appl Opt 1969;8:725. Wilkes MV. Proc Phys Soc B (London) 1954;67:304. Banks PM, Kockarts G. Aeronomy. Part B. New York: Academic, 1973, p. 107. Brasseur G, Solomon S. Aeronomy of the middle atmosphere. 2nd ed. Dortrecht: Reidel, 1986, p. 106. Pressly EC. Phys Rev 1953;89:654. Rees MH. Physics and chemistry of the upper atmosphere. Cambridge: Cambridge University Press, 1989. Abramowitz M, Stegun IA. Handbook of mathematical functions, Vol. NBS AMS 55. Washington, DC: USGPO, 1964. Hildebrand FB. Introduction to numerical analysis. New York: McGraw-Hill, 1956. Rabenstein AL. Introduction to ordinary di!erential equations. New York: Academic Press, 1966. Dwight HB. Tables of integrals and other mathematical data, 4th ed. New York: Macmillan, 1961. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical recipes: the art of scienti"c programming. London: Cambridge University Press, 1986.