Accurate analytic approximation for the Chapman grazing incidence function

A new analytical approximation for the Chapman mapping integral, Ch\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {Ch}}$$\end{document}, for exponential atmospheres is proposed. This formulation is based on the derived relation of the Chapman function to several classes of the incomplete Bessel functions. Application of the uniform asymptotic expansion to the incomplete Bessel functions allowed us to establish the precise analytical approximation to Ch\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {Ch}}$$\end{document}, which outperforms established analytical results. In this way the resource consuming numerical integration can be replaced by the derived approximation with higher accuracy. The obtained results are useful for various branches of atmospheric physics such as the calculations of optical depths in exponential atmospheres at large grazing angles, physical and chemical aeronomy, atmospheric optics, ionospheric modeling, and radiative transfer theory.


Introduction
In modeling of atmospheres their chemical constituents are often assumed to be distributed according to the hydrostatic equilibrium within the column extending from a point of interest towards the sun. This columnar content can be conveniently mapped to the reference value via the Chapman grazing incidence function Chapman (1931). The reference value corresponds to the vertical column, i.e., to the situation when the sun is at the zenith of the point of interest. As an approximation to the hydrostatic equation the distribution of atmospheric constituents is assumed to follow the exponential distribution within the vertical column. Being originally introduced for the description of energy deposition in the ionosphere due to the absorption of solar radiation, the Chapman function is now an important mathematical object in other branches of the atmospheric physics, such as radiative transfer theory (Dahlback and Stamnes 1991), chemical and physical aeronomy (Brasseur and Solomon 2005;Schunk and Nagy 2009), description of light scattering in planetary atmospheres (Wallach and Hapke 1985). The function is also used in modeling of ionospheric D-layer (Friedrich and Torkar 2001), in calculations of atmospheric transparency (Green et al. 1964), optical air masses (Kocifaj 1996;Rapp-Arrarás and Domingo-Santos 2011), cosmic ray ionization rates (Velinov 1974), or for air shower simulations (Extensive and Showers 2010), to name just a few.
In general case the evaluation of the Chapman function requires numerical integration, which, in turn, might require considerable computational resources in complex physical models of exponential atmospheres. This is the reason why many attempts have been made in order to provide simple analytic approximative expressions for the Chapman integral. The most notorious approximations of the Chapman function have been collected in the D. Huestis work (Huestis 2001). This work is also an excellent source on numerical and asymptotical evaluation of the various integral representations of this mapping function.
The aim of the present article is to draw the reader attention to the less known interrelation of the Chapman function and the incomplete cylindrical functions (Agrest and Maximov 1971). Based on this observation we obtain new asymptotic representations of the Chapman function, which approximate the exact integral expressions with the high order of accuracy. One of the derived approximations yields better approximation of the exact Chapman integral than other analytic approximations found in the literature. We hope that the obtained result would be useful for physical models of ionospheres, theories of optical propagation in the atmospheres, radiative transfer theories, and in remote sensing applications The article is organized as follows. In section "Definitions and current analytic approximations" we summarize the basic definitions and properties of the Chapman function and give an overview of its existing analytic approximations. In section "Derivation of new analytic approximations" we derive the relationship of the Chapman integral with several classes of the incomplete cylindrical functions. Based on these relationships and using the uniform asymptotic expansions of the incomplete Macdonald integrals we derive analytic approximations of the Chapman mapping function. The methodological aspects of comparison of different analytic approximations are discussed in "Methodology of comparison". The performance of the obtained approximations is discussed in section "Results of comparison and discussion" and is compared with other approximations found in the literature. The obtained results are summarized in "Conclusions".

Definitions and current analytic approximations
The derivation of Chapman mapping function relies on several simplifications concerning the properties and geometry of atmospheric layers. Firstly, the spherical symmetry of the atmospheric shells is assumed. Secondly, the hydrostatic equilibrium within the atmospheric shells is expected, meaning that in a state of rest, the atmospheric pressure at specified height must support the weight of the gas above it. Additionally to this the temperature and the gravitational acceleration are supposed not to vary with the altitude. This restriction results in exponential distribution of an ideal gas in hydrostatic equilibrium. Finally, the bending of solar rays due to atmospheric refraction is neglected.
Under specified restrictions the classical representation of the Chapman function as derived in Chapman (1931) reads as where ζ is the solar zenith angle at the point of interest and (1) is the dimensionless ratio of the spatial coordinate L measured from the center of the Earth to the scale height H of particular atmospheric gas. The scale height characterizes the exponential decay of the atmospheric constituent density with the height. For an ideal gas this quantity is expressed through the Boltzman's constant k, the temperature T, the mean constituent mass m, and the acceleration due to gravity g as If the spherical earth model is under consideration, we have L = R ⊕ + z with R ⊕ being the earth radius and z being the hight of the point of interest. The derivation of Eq. (1), based on geometrical considerations, is summarized in Swider (1964).
To avoid numerical integration of Eq. (1) several analytic approximations have been proposed in the literature. These approximations are based either on asymptotic evaluation of the integral or on semi-empirical considerations. Table 1 summarizes the most widely used approximative Chapman mapping functions. Most of the expressions are given in terms of the complimentary error function which in contrast to Eq. (1) being a single argument function is properly tabulated and approximated, cf. Abramowitz and Stegun (1972). The formulas in Table 1 are given for the solar zenith angle 0 ≤ ζ ≤ π/2 and can be extended to larger zenith angles via the relationship (Chapman 1953;Titheridge 1988): Here and K 1 is the modified Bessel function of the second kind of the first order (Macdonald function), cf. Abramowitz and Stegun (1972). At the limiting case of ζ = π/2 we obtain the Chapman function at the horizon for the point of interest: Due to Eq. (3) we restrict our attention in the following to zenith angles in the range ζ ∈ [0, π/2]. (2) Ch(X, ζ ) =2 exp [X(1 − sin ζ )]Ch(X sin ζ , π/2) − Ch(X, π − ζ ), ζ > π/2.

Relation to incomplete cylindrical functions
Here we show the relationship of the Chapman function to several families of the incomplete cylindrical functions. To start with, we change the integration variable t in Eq. (1) according to s = csc(t) . In this way we arrive at the alternative integral representation of the Chapman function where Integrating Eq. (6) by parts we obtain From this formula it is easy to derive the result (5) for ζ = π/2 , i.e., cosh β = 1 , by using the integral representation of the Macdonald function (Abramowitz and Stegun 1972, Eq. 9.6.23) where Ŵ(ν) is the Gamma function.
In the view of Eq. (9) the Chapman function (8) can be rewritten as with being the incomplete Macdonald function in the form introduced in Agrest and Maximov (1971). The nomenclature "incomplete cylindrical function" or "incomplete Bessel function" becomes clear from this definition as only in the limit y → ∞ the function (11) becomes the "complete" Macdonald function (9). Table 1 Analytical approximations of the Chapman function (in chronological order). The solar zenith angle ζ is given in radian units

Index Reference Approximation of the Chapman function
(a) Green et al. (1964); Green and Barnum (1963) Ch(X , ζ ) ≈ exp (d) Titheridge (1988Titheridge ( , 2000 Ch Kocifaj (1996) In order to relate the Chapman function to other types of incomplete cylindrical functions we perform the variable change in (8) according to s = cosh θ and obtain This expression allows one to express the Chapman function through the incomplete Macdonald functions as defined by Jones (2007) or as defined by Terras (1981). Namely, using the definition for K n (x, y) the Chapman function can be represented as The corresponding relationship between the Chapman function and the incomplete cylindrical function (14) is easily obtained from (15) by observing the obvious relation The recurrence relations also allow to establish the relationship between the incomplete Bessel functions in form of Agrest and Maximov, Jones, and Terras. For example, we have in our case of interest the following relations and (12) We also note that the limiting case (5) for ζ = π/2 is obtained from (15) by observing that K n (x, 0) = K n (x) 1 and with the use of the well-known recurrence relations for the Bessel functions.
The recently discovered relationship of the Chapman function to the solution of an inhomogeneous modified Bessel's equation (Huestis 2001) yields the following representation of this mapping function: where I n (x) are the modified Bessel functions of nth order. The advantage of this integral representation is that the Champan function in this formulation is stable for all values of X and ζ as has been shown by Huestis (2001). The integrals in Eq. (19) cannot be analytically integrated in general case except for X = 0 . One can show that these integrals possess a weak singularity for t → ∞, ζ = π/2 while I 2 possesses an integrable singularity for t → 0, ζ → 0.
Equation (19) allows one to establish the relation of the Chapman function and another class of the incomplete cylindrical functions known as the incomplete Lipschitz-Hankel integrals (Agrest and Maximov 1971). Namely, the integrals (20) can be expressed through the incomplete Lipschitz-Hankel integrals (Agrest and Maximov 1971) Alternatively the property is applied if the Terras' representation for incomplete Bessel functions is used. Planets and Space (2021) 73:112 as Here the parameter β is given by Eq. (7). Substituting Eqs. (22), (23) in Eq. (19) we obtain This representation reduces Ch(X, π/2) directly to the limiting case (5).
Using Eq. (19) we now derive the relationship of the Chapman function to the incomplete Bessel functions (13). We start with the result obtained in Ref. Agrest and Maximov (1971): where the function L 0 is defined in Eq. (11). Differentiating both sides of Eq. (25) with respect to x and using the relations dI 0 (x)/dx = I 1 (x) and dK 0 (x)/dx = −K 1 (x) we obtain This equation allows one to rewrite the Chapman function (19) as where the property K −n (x, y) = K n (x, y) has been used. This is alternative to Eq. (15) representation of the Chapman function in terms of the incomplete Macdonald functions (13).
To conclude this section we relate the Chapman function with another class of the incomplete cylindrical (21)

functions, namely the incomplete Weber integrals (Agrest and Maximov 1971)
These integrals appear in the theories of optical diffraction, pulsed radar theories, and mathematical statistics, see e.g. Vasylyev et al. (2013). Reference (Agrest and Maximov 1971) derives the relationship between the incomplete Weber integral Q 0 (x, y) and the incomplete Lipschitz-Hankel integral Ie 0 (x, y) as given in Eq. (21). In our notations this relationship for 0 ≤ ζ ≤ π/2 reads as: The incomplete integral Ke 0 (x, y) can be expressed in terms of Ie 0 (x, y) as shown in Ref. Agrest and Maximov (1971) as where L n (x, y) is given by (11). Combining Eq. (29) with (30) and substituting them in (24) we obtain where we have used Eq. (17) for the incomplete Bessel function L 0 (x, y) . The obtained expression is rather complex to be applied in practice. Its derivation serves merely for sake of completeness of the relations between the Chapman function and the incomplete cylindrical functions of different kinds. However, it appears that (31) becomes useful for the analytic approximation of the Chapman function when the scale parameter is small.

Analytic approximations for the mapping function
In the previous section we have obtained several formulas that relate the Chapman mapping function (1) to the incomplete cylindrical functions of various types. Specifically, Eq. (10) relates the Chapman integral to the incomplete Macdonald function L n (x, y) of the Agrest and Maksimov type, whereas Eqs. (15) and (27) express the former in terms of the incomplete Macdonald function K n (x, y) of the Jones type or alternatively in terms of the Terras incomplete Macdonald function K n (x, v) . Finally, Eq. (24) is the representation in terms of the incomplete Lipschitz-Hankel integrals Ie n (x, y), Ke n (x, y) and can be also rewritten in terms of the incomplete Weber integral Q n (x, y) and the incomplete Macdonald function as given in Eq. (31). Since there are relationships among these various types of incomplete Bessel functions as discussed in the previous section, we consider the analytic approximation of the Chapman function using Eqs. (15) and (27) written in terms of K n (x, y) . The special case of small (29) Ie 0 (csc ζ , X sin ζ ) = tan ζ 2e −2X sin 2 ζ 2 Q 0 X sin 2 ζ 2 , X sin ζ − 1 − e −X I 0 (X sin ζ ) .
(30) Ke 0 (csc ζ , X sin ζ ) = K 0 (X sin ζ ) I 0 (X sin ζ ) Ie 0 (csc ζ , X sin ζ ) Ch(X, ζ ) =X sin ζ e X K 1 (X sin ζ ) + I 1 (X sin ζ ) values of the scale parameter X, which is not of interest in atmospheric physics applications, is based on Eq. (31) and is considered in the Appendix A.
dictated by suitability of its integral representation (13) for uniform asymptotic expansions. We refer to the work (Jones 2007) for details in derivation of the expansion formula and give here the lowest terms in the asymptotic expansion of K n (x, y) Here β is given in Eq. (7) and (32) Inspection of Eq. (36) shows that the first term on the r.h.s resembles the asymptotic expansion of the Macdonald function K 1 (X sin ζ ) with the complementary error function. This observation shows the similarity of the obtained asymptotic expansion formula (36) with the largeζ expansion of Chapman,cf. Eq.(38) in his original article (Chapman 1931). If one performs the expansion of the Macdonald function up to the second term in the first summand of the Chapman asymptotic expression, one obtains exactly the first term of Eq. (36). With this respect it is worth to note that Chapman considered the asymptotic of small and large zenith angles separately. Moreover, Chapman's original asymptotic expansions are complex polynomial series with poor convergence, the issue corrected in Ref. Huestis (2001). Equation (36) originates from uniform (35) asymptotic expansions and is valid for all values of zenith angle and is simpler to implement in practice.
To conclude this section we provide an explicit expression for approximation (36) in the form useful for numeric calculations Here the zenith angle enters as an argument of sine function only and the expression sin ζ is optimal to be evaluated prior the substitution in (37). Finally we note that for sufficiently large values of X as ζ → 0 the formula (37) approaches the expected limit sec ζ.

Methodology of comparison
Fortran code was written to evaluate Ch(X, ζ ) through the analytical approximations (35), (36) and (43). These equations has adopted in the form optimal for numerical calculations, cf. Eq. (37). For evaluation of the numerical integral (1) and some of the previous analytical approximations of Table 1 the code written by Huestis (2001) has been adopted. We have chosen the double precision calculations for both approximative formulas and numerical integral (1) to minimize the effect of rounding errors.
The code was compiled and tested using the Intel Fortran Compiler under Windows operating system. Based on the Fortran program the Python package has been created using F2Py Fortran-Python interface generator and has been implemented as the corresponding wrapper program.
In order to show the precision of various analytical approximations we plotted the logarithm of the relative error as a function of the solar zenith angle in Figs. 1 and 2 for scale parameters X = 50 and X = 500 , respectively. The choice of the scale parameters is arbitrary and chosen to satisfy X ≫ 1 and to represent the values different in one order of magnitude. The relative error is defined as where Ch (e) (X, ζ ) is the exact value of the Chapman function taken from the double precision numerical integral (37) Ch (e) (X, ζ ) , (1) and Ch (a) (X, ζ ) is the analytical approximation of this function calculated with the double precision. We use the relative error instead of the absolute error ε a = |Ch (e) − Ch (a) | due to the following reasoning. For applications in atmospheric physics the values of the scale parameter that satisfy X ≫ 1 are of interest. For ζ ∈ [0, π/2] in this asymptotic one can show that is the order symbol (Erdélyi 2003). These relations follows from Eq. (5) and from the asymptotic expansion of the Macdonald function. In order to remove the O( √ X) trend in the error analysis the relative error has been chosen. Similar approach has been chosen by D. L. Huestis in his thorough analysis of Chapman function evaluation (Huestis 2001). The current error analysis is thus fully in agreement with those performed in Ref. Huestis (2001) apart from the fact that D. L. Huestis chosen the scale parameters X = 50 and X = 800 and we use X = 50 and X = 500 in the current study. We have also checked the correspondence of several approximative formulas from Table 1 with the error analysis of corresponding expressions in Ref. Swider and Gardner (1967) for X = 50 and Fig. 3 Relative error as a function of the scale parameter for small (dashed) and large (solid) solar zenith angles. Three analytical formulas are compared: a is Eq. (36) obtained in this work, b is the Kocifaj formula (Kocifaj 1996) that shows the best performance among analytical approximations found in the literature, and c is the empirical formula of Green and Barnum (1963); Green et al. (1964) X = 500 . A good agreement has been found despite of difference in some computational aspects.
As has been already mentioned, the majority of analytic approximations for the Chapman function involve the complimentary error function. For all such expressions we have used the specific rational approximations for this function due to Smith and Smith (1972) In this way we reduced the approximation of Chapman function due to Smith & Smith to the Fitzmaurice approximation (Fitzmaurice 1964). We note that the approximative formula (39) originates from the tables of approximations collected in Ref. Hart (1968) and yields the absolute error in double precision calculations of orders 10 −5 and 10 −3 for x ∈ [0, 8] and x ∈ (8, 100] , respectively.

Results of comparison and discussion
The plotted in Figs. 1, 2 and 3 relative errors (38) for the analytical approximations of Chapman function [see Table 1 and Eqs. (35), (36)] exhibit different performance. For example, the Green & Barnum empirical formula shows the pure performance for zenith angles ζ → π/2 as has also been noted in Huestis (2001), see Figs. 2(a) and 3(c). The Fitmaurice formula reaches the discrepancy of several percents from the numerically obtained values at moderate vales of X, e.g. around 1% at X = 50 , see Fig. 1(b). Due to its analytical simplicity this function gained the popularity in modeling of the atmosphere. For example, the Thermosphere Ionosphere Electrodynamics General Circulation Model (Dickinson et al. 1981) or the Gauss-Seidel limb scattering radiative transfer model (Loughman et al. 2015) use the Fitzmaurice formula and hence it also incorporates these discrepancies from the exact formula. With this respect it is worth to point that the Kocifaj formula (curve (e) in 1, 2) has an attractive property of being quite accurate and analytically simple. As one can see this formula, remaining rather unknown to atmospheric physics community, performs even better than the highly-accurate formula of Swider [curve (c)] or the approximation due to Huestis [curve (f )].
Another accurate analytic approximation of the Chapman integral is the semi-empirical formula of Titheridge introduced in Titheridge (1988) and corrected in Titheridge (2000), see curve (d). The discrepancy of this approximation from the Chapman function grows however for large values of X and when ζ → π/2 as can be seen in Fig. 2.
and (27). As has been shown in the previous section the asymptotic expansion of Eq. (15), cf. curve (g), reduces to the zeroth-order term in uniform asymptotic expansion of Eq. (27), cf. curve (h). Still, being only the first term given by (35), it outperforms the analytical approximations of Fitzmaurice [curve(b)] and the first-term of asymptotic expansion of Huestis [curve (f )]. The analytical approximation (36), curve (h), shows the best performance among the analytical approximations. Oscillatory behavior of this curve arises from variation in accuracy of polynomial approximations to the mathematical function erfc(x) and of the expression √ π/2x(1 + 3/8x) exp(−x) being the asymptotic expansion of Bessel function K 1 (x).
In order to illustrate the dependence of the obtained approximation formula (36) on the scale parameter, we compare it with the analytical approximation of Kocifaj and Green & Barnum in Fig. 3. The curve (a) calculated with Eq. (36) performs better than the Kocifaj formula [curve (b)] or Green & Barnum formula [curve (c)] for all relevant scale parameters values in atmospheric applications. It is worth to note that in some applications (e.g., radiative transfer and radiometric sounding), where the Chapman mapping function is used as a core function inside multiple integrals, the complexity of derived analytical approximations (35), (36) will require additional computational resources and analytically simpler and accurate Kocifaj's approximation is preferable. Finally, we remark that for small values of the scale parameter the approximation (43) of Appendix A performs better than other analytical approximations, however this case is less interesting for known application cases. Fig. 4 Relative error as a function of the scale parameter in the asymptotic X ≤ 1 . The curve a corresponds to Eq. (36), the curve b stems from the Kocifaj result (Kocifaj 1996), the curve c is due to Green and Barnum (1963), and the curve d is obtained from Eq. (43). The cases of small (dashed) and large (solid) solar zenith angles are shown