A fast and accurate DE-formula algorithm to evaluate Ambartsumian-Chandrarasekhar H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$H$\end{document}-function for isotropic scattering

In the present work, a compact, fast, and yet accurate algorithm is developed to calculate the numerical values of the Ambartsumian-Chandrasekhar’s H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$H$\end{document}-function for isotropic scattering and its moments on the basis of the double exponential (DE) formula of Takahashi and Mori (RIMS, Kyoto Univ., 9:721, 1974). The main improvement made in the new method is an elimination of the iterative procedure for automatic adjustment of the step-size of integrations carried out with the DE-formula. Instead, a set of optimal values for the upper limit of integration Tmax\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$T_{\text{max}}$\end{document} and the number of division points NT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$N_{\text{T}}$\end{document} to specify the step-size of the quadrature is predetermined for calculations of the H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$H$\end{document}-function with a 15-digit accuracy, and also another for the evaluations of the moments with an accuracy of 14-digits or better. FORTRAN90 subroutines HFISCA for the H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$H$\end{document}-function and HFMOMENT for the moments of arbitrary degrees are subsequently constructed (their source codes and a driver together with a sample set of output are shown in Appendix of this paper). Tables of sample calculations of the H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$H$\end{document}-function and its moments of degree −1 through 6 carried out by these programs are also presented. The routines HFISCA and HFMOMENT should prove useful not only in astrophysical applications but also in other disciplines of science such as the electron transports in condensed matter and remote-sensing data analyses. A request for a copy of the Fortran 90 source code of the program can be made by writing to kawabata@rs.kagu.tus.ac.jp.


Introduction
The intensity of radiation emerging from a plane-parallel, semi-infinite, vertically homogeneous absorbing-scattering layer illuminated by mono-directional light can be expressed using the Ambartsumian-Chandrasekhar's H -function (see, e.g., Ambartsumian 1942Ambartsumian , 1943aAmbartsumian ,b, 1958Chandrasekhar 1950Chandrasekhar , 1960Kourganoff 1963;Ivanov 1973;Sobolev 1975;Irvine 1975;van de Hulst 1980;Kelley 1980;Lenoble 1985;Chamberlain and Hunten 1987;Yanovitskij 1997;Nagirner and Ivanov 2020;Frisch 2022). The expressions in terms of the H -functions are useful particularly for modeling the angular distributions of such emergent intensities when the media comprising the layer are giving rise to mildly anisotropic scattering. They can also be used to produce ref-erence values for these quantities to check the results from computer codes to handle more realistic cases of anisotropic scattering. Or they can be employed to conveniently approximate the emergent intensities from anisotropically scattering layers of large but finite optical thicknesses (see, e.g., Sobolev 1975;Irvine 1975). Of particular interest is the work of Viik (1986), who succeeded in calculating with 14digit accuracy the values of the H -functions for isotropic scattering, Rayleigh scattering and linearly anisotropic scattering cases by approximating Sobolev's resolvent function using exponential series. More commonly used method is to solve iteratively the integral equations satisfied by the Hfunctions (see, e.g., Chandrasekhar 1950Chandrasekhar , 1960Kolesov and Smoktii 1972;Bosma and de Rooij 1983;Lenoble 1985;Kawabata 2015Kawabata , 2016Kawabata , 2018. The most basic is the H -function for isotropic scattering, for which some explicit integral representations of the solution of the defining integral equation are available depending on the choice of the integration variable (e.g., Stibbs and Weir 1959;Ivanov 1973;van de Hulst 1980;Rutily and Bergeat 1987). With the DE-formula for numerical integrations, the values of the H -function can be obtained rather straightforwardly with more than 14 or 15-digit accuracy (Mohankumar and Natarajan 2007;Kawabata 2016;Jablonski 2019Jablonski , 2020. The isotropic scattering H -function still plays a significant role not only in astrophysical applications, but also in other disciplines of science such as remotesensing (see, e.g., Heylen et al. 2014) and photoelectron transports in condensed matters (see, e.g., Jablonski 2012Jablonski , 2015Jablonski , 2019Jablonski , 2020, where numerical values of high accuracy (more than 10 significant digits) are required for this function.
Lately, Jablonski (2020) has published a set of Fortran 90 program H_FUN_v2 to calculate the values of the Hfunction for isotropic scattering for full parameter ranges presumably with accuracy of 15 or 16 decimals: actual calculations are performed in a subroutine subprogram HFU-NIV using either an integral representation of the H -function given in Ivanov (1973) applying the DE-formula or an approximate formula for the H -function for high-speed executions depending on the values of the parameters.
In view of the importance of the H -function for isotropic scattering, it must be of great interest to make a design of an algorithm for numerical evaluation of this function based solely on the DE-formula and assess to what extent it can compete the capability of HFUNIV, which we would like to undertake in the present work. If successful, we shall apply our method to numerical calculations of the higher order moments of the function.

Formulation for the DE formula integration
Suppose we wish to carry out the following definite integral numerically for a given function f (x): (1) where the integrand function f (x) may have a singularity or singularities at either x = a or x = b or both. Under such circumstances, the double exponential formula developed by Takahashi and Mori (1974) (see also Mori and Sugihara 2001;Mori 2005;Kyoya and Tanaka 2019;Bailey et al. 2000), which we shall hereafter refer to as the DE-formula, is known to be highly effective (its optimality has been rigorously proven by Sugihara 1997), so that we shall employ this DE-formula to evaluate the Ambartsumian-Chandrasekhar's H -function for isotropic scattering, since its integral representation possesses exactly this type of singularities. Although the general recipe for numerical integrations based on the DE-formula are available in various literatures such as Mori (1990Mori ( , 2005, Watanabe (1990), and Mori and Sugihara (2001), we present the DE-formula in this section in the form suitable for numerical calculations of the H -function.
Transforming the independent variable x to The central procedure of the DE-formula is to employ the following variable transformation ξ = φ(t) that maps x = a to t = −∞ and x = b to t = +∞: so that where we have defined the weight function w(t) in Eq. (4) as Then, Eq.
(2) takes the following form: which can be further rewritten as with the function η(t) being defined by Noting that for t 1, the weight function given in Eq. (5) decreases double-exponentially as we can truncate the integration in Eq. (7) at a sufficiently large upper limit T max such that Discretizing Eq. (10) with a step size h for the variable t and applying the trapezoidal rule for the integral, we get the following DE-formula for definite integrals on the interval [a, b]: where A notable computational advantage associated with the DEalgorithm is the ease with which an adaptive automatic integration (hereafter referred to as AAI) procedure can be implemented. In fact, the computer programs of Kawabata (2016Kawabata ( , 2018 and Jablonski (2020) possess such capability. However, since the AAI procedure is likely to be a principal cause for slowing down the DE calculations of the Hfunction, we would like to investigate if it is possible to carry out numerical evaluation of the H -function for isotropic scattering without the AAI still maintaining the same degree of accuracy as that realized by the DE-method equipped with the AAI procedure. Along this objective, we shall instead search for an appropriate set of the values for the principal parameters N T and T max that can minimize the largest absolute error involved in our values of H ( 0 , μ) in comparison with reference values.

Application of the DE-formula to H -function evaluation
The following is one of a few integral forms available to express the Ambartsumian-Chandrasekhar H -function H ( 0 , μ) for isotropic scattering (see, e.g., Stibbs and Weir 1959;Kourganoff 1963;Ivanov 1973;van de Hulst 1980;Rutily and Bergeat 1987;Limaye 2011, 2013;Kawabata 2016; Jablonski 2020 and the references cited therein): The parameter 0 involved in Eq. (13) is the single scattering albedo (0 ≤ 0 ≤ 1) of the isotropically scattering agents comprising a semi-infinite and vertically homogeneous layer, while the parameter μ (0 ≤ μ ≤ 1) is the cosine of the zenith angle θ (0 ≤ θ ≤ π/2) measured from the outward normal of the upper surface of the layer, and specifies the direction of a ray of radiation incident on or emergent from that surface. Eq. (13) is of the following form: where I indicates the integral with respect to x. Applying to this quantity the procedure employed to obtain Eq. (11) from Eq. (1) with a = 0 and b = π/2, we have with w i specified by Eq. (12d). The function f (x; 0 , μ) has singularities at x = 0 and x = π/2 in addition to a sharp spike-like variation depending on the values of 0 and μ. This is why an accurate numerical integration of Eq. (13) with respect to x is hard to achieve by means of a standard quadrature system such as the Simpson's rule as employed in Stibbs and Weir (1959), or even the Gauss-Legendre formula as used in .
Considering that the function f (x; 0 , μ) needs to be evaluated whenever the values of 0 and/or μ are changed, it is more convenient to isolate the parts of Eq. (15c) that depend solely on the running variable x and express it in the following form: where we have defined The quantities Q(x), U(x), and V (x) are just functions of x only, and are independent of both 0 and μ, so that they can be precomputed and stored in tabular form for repeated use in evaluations of H ( 0 , μ) for arbitrary values of 0 and μ. It should however be pointed out that the values of V (x) tend to become numerically unreliable particularly for x < 10 −7 , if an intrinsic function for tan x supplied by any ordinary Fortran compiler is utilized: 0 would be returned to V (x) well before we reach x = 0. Taking into account this fact, it is safer to make an appropriate approximation for V (x) for small values of x lest a logarithmic divergence should happen to f when 0 = 1. 1 A recourse is therefor made to Eq. (7) of Kawabata (2016), i.e., a power series approximation for V (x), for x ≤ x trans = 0.1238798627279953: where we have designated s = x 2 . We shall therefore employ the following: (19) Figure 1 shows the weight function w(t) given by Eq. (5)(solid line curve) and the integrand function Eq. (15c) for 0 = 1 and μ = 10 −3 (dotted line curve) as an example. We notice the presence of two steep ravines in the latter, viz., Fig. 1 The solid line curve shows the weight function w(t) (see Eq. (5)) for the DE-formula (to be referred to the bottom abscissa and the left-hand side ordinate), whereas the dotted line curve is the integrand f (x; 0 , μ) (see Eq. (15c)) with 0 = 1 and μ = 10 −3 (to be referred to the top abscissa and the right-hand side ordinate)

Fig. 2
New integrand f (x(t); 0 , μ)w(t) with 0 = 1 as a function of the DE-formula variable t . The numbers 1 through 7 signify seven values employed for μ, viz., 1: μ = 10 −15 , 2: μ = 10 −10 , 3: μ = 10 −5 , 4: μ = 10 −3 , 5: μ = 0.1, 6: μ = 0.3, and 7: μ = 1 one at x = 0 and the other at x = 1.5843775. Such characteristic features of f (x; 0 , μ) are the principal reason why a careful treatment is required in evaluation of H ( 0 , μ) as has already been mentioned. Figure 2 shows seven samples of f (x(t); 0 , μ)w(t) keeping 0 = 1 unchanged. No contribution to the integral is made from the integrand at x = 0 nor at x = π/2 due to the fact that their functional values are damped to a sufficient degree by the presence of the weight function w(t), thereby allowing us to bypass the procedure of the semianalytical integration Eq. (8) introduced in  to evaluate I . In addition, as a comparison between the integrand f in Fig. 1 and the curve No. 4 in Fig. 2 tells us, the behavior of f (x(t); 0 , μ)w(t) is significantly smoother than that of f (x; 0 , μ) itself, which helps increase the accuracy of numerical integration for H ( 0 , μ).

Numerical evaluation of H ( 0 , μ)
Based on Eqs. (15a)-(15c), and Eq. (16), we have produced a FORTRAN90 function subprogram (see Appendix for the source code), which we call HFISCA (H -Function for Isotropic SCAttering), involving two essential parameters N T and T max for calculating the values of H ( 0 , μ). Our principal task is then to determine a single set of numerical values of these parameters in such a way that the values of the H -function be generated with the highest accuracy possible. In order to actually assess their accuracies, we make calculations of the H -function on a 222 × 222 square mesh of ( 0 , μ) specified in the following manner (hereafter referred to as Domain A): 2 13), All the calculations in the present work are performed in double precision arithmetics and the results are compared with reference values produced independently by another Fortran90 function subprogram HFITER employed in Kawabata (2018) to iteratively solve the integral equation(shown later by Eq. (24a)) defining the H -function for isotropic scattering: this routine makes use of a RATFOR program DEAUTO published by Watanabe (1990) for performing numerical integrations employing the DE-formula and an AAI procedure.
In order to determine the most appropriate set of values for (N T , T max ), we first look into the maximum absolute error ε max (N T , T max ) involved in the values of the H -function calculated by HFISCA on the Domain A for a given set of values for (N T , T max ), from which we locate the smallest values ε min (N T ). Finally, the minimum value ε of those of ε min (N T ) should be found, and the values of N T and T max giving rise to this ε will be adopted as what we desire: where the notation ε min (N T ) tacitly presumes that the relevant value of T max is to be separately registered as a function of N T , and a similar remark also applies to the notation ε with respect to N T . Based on some preliminary experimentation, we have found it sufficient to make the above search for numerical candidates for N T and T max in the following ranges: enabling us to cover 60 ≤ N T ≤ 125 and 3 ≤ T max ≤ 3.25. In Fig. 3, the heavy filled circles with a label "a" show in logarithmic scale the values of ε max (N T , T max ) thus found as functions of N T in the range between 65 and 120 (l = 6, 7, . . . , 61). The solid lines indicated by a label "b" connect the data points possessing ε min (N T ), while the horizontal right arrow indicates the ordinate location of those with ε = 1.776 × 10 −15 , the minimum absolute error of our current values of the H -function attainable on the Domain A.
Obviously from the diagram, we find the minimum absolute error occurs not at a single data point but at 17 values of N T , viz., 99 through 115, and 528 values of T max , although the number of values of the latter parameter associated with each of the 17 values of N T differs: in the case of, e.g., N T = 107, we have 68. In other words, we need to impose some additional discriminatory criteria on these 528 sets of (N T , T max ) in order to arrive at a unique best pair.
Let us therefore inspect the maximum values of the root mean square errors δ RMS (N T , T max ) of these 528 sets of the values of H -function calculated again with HFISCA on the Domain A: where the values of T max in Eqs. (23a), (23b) and (23c)  The values of N T , T max , and δ min thus derived for the 17 cases are summarized in columns 2, 3, and 4 of Table 1 with column 1 giving their numbering. The smallest value δ of δ min (N T , T max ) is found to be 1.702 × 10 −16 given by the case No. 17, where we have N T = 115 and T max = 3.2053. However, we also notice comparable values 1.703 × 10 −16 for the case No. 15 and 1.705 × 10 −16 for the case No. 13. Hence, we should postpone our final decision on the values of N T and T max until a few more checks have been made.
Another check worth considering would be to see how well the values of H ( 0 , μ) HFISCA obtained for the 17 cases of (N T , T max ) satisfy the integral equation defining the Hfunction (see, e.g., Chandrasekhar 1960;van de Hulst 1980;Bosma and de Rooij 1983): The numerical integrations required to calculate the righthand side of Eq. (24a) are carried out using the routine DEAUTO of Watanabe (1990) except that the value of the parameter EPS is set to 10 −15 instead of 10 −12 for a stricter convergence criterion for the AAI procedure. The values of ε IE and δ IE obtained for the 17 cases of (N T , T max ) are shown in columns 5 and 6 of Table 1, respectively. The smallest value 1.776 × 10 −15 for ε IE is found at No. 9, 11, and 13, while the case No. 13 has the smallest value 1.720 × 10 −16 for δ IE . So far, we have checked the factors that tend to give an edge to larger values of N T . Another crucial factor from a practical point of view must be the CPU time used by the routine HFISCA for execution. Let us therefore make a timing comparison by measuring the CPU-time t CPU (No.) required to perform calculations of the H -function on the Domain A: the arithmetic mean of 100 measurements of the CPU-time is taken for each case and the results of measurements are presented in column 7 of Table 1 under the heading "t * " in the form of a fractional ratio t * ≡  Table 1 with the header "Product" shows the values of the product of the four quantities given in columns 5 through 7 and column 9 gives the rankings of the 17 cases according to these products such that the smallest comes first. We therefore select the case No. 13 as the best, which in turn implies we adopt N T = 111 and T max = 3.2215 for evaluation of the H -function. Table 2  The bottom row for ε F μ ( 0 ) (see Eq. (26) below) shows, for a given value of 0 , the maximum absolute error involved in 222 values H ( 0 , μ j ) (j = 1, 2, . . . , 222) with μ j specified by Eq. (20a): 1, 2, . . . , 22).
The maximum of ε F μ ( 0 ) is found to be 1.776 × 10 −15 taking place at 0 = 1 − 10 −7 , which coincides with ε of the values of H ( 0 , μ) HFISCA obtained on the Domain A. The values of the H -function produced by HFISCA should therefore be correct at least to the 14-th decimal place, and have the accuracy of 15 significant digits including the integer part.

Comparison with a competitive program for H -function
Whether or not the routine HFISCA is of some value must be assessed by comparing with a leading computer program or programs made publicly available for numerical calculations of the H -function. To the best of our knowledge, the For-tran90 program HFUNIV developed by Jablonski (2020) is best suited for this purpose. This program is consisted of two parts: one is HFDE to calculate the values of the H -function for isotropic scattering by means of the DE-formula, and the other is an approximation formula for the H -function to be employed for small values of 0 less than or equal to some demarcation values specified by the value of μ. The expression adopted for H ( 0 , μ) can be derived by putting x = arctan(s) in Eq. (13) of the present work and involves an integral from 0 to ∞ with respect to the new variable s. Applying the variable transformation required for the DEformula to this integral and truncating the resulting infinite integral at ±T max = ±6, we arrive at the DE-formula equation Jablonski (2020) employs. It should also be pointed out that the routine HFDE performs AAI as in Kawabata (2016Kawabata ( , 2018. It is claimed that the program HFUNIV is capable of producing the numerical values of the H -function with    (3) δ: the same as (1)  an accuracy of 15 to 16 significant digits for the parameter ranges 0 ≤ 0 ≤ 1 and 0 ≤ μ ≤ 1. First, we would like to make comparisons between HFISCA and HFUNIV in two crucial aspects, viz., 1) the CPU timing for calculation of the 222 × 222 values of H ( 0 , μ) on the Domain A (0 ≤ 0 ≤ 1; 0 ≤ μ ≤ 1), 2) numerical accuracy in terms of the maximum absolute error ε and the root mean square error δ: where the reference values are those produced by the routine HFITER of Kawabata (2018) as in the foregoing sections.
The results for (a) HFUNIV and (b) HFISCA are shown without parentheses in columns 2 through 6 of Table 3. The figure given in column 2 of (a) is the mean of the 100 measurements of the CPU-time taken by HFUNIV relative to that for HFISCA. The execution speed of HFISCA is found to be faster than that of HFUNIV by a factor of slightly less than five. Elimination of the AAI procedure from the algorithm of HFISCA is likely to be the primary cause for this speed-up.
Furthermore, the maximum absolute error ε of H ( 0 , μ) HFUNIV is 27.53 × 10 −15 , which takes place, e.g., at 0 = 1 − 10 −13 and μ = 0.998, where we have H ( 0 , μ) HFUNIV = 2.9042700256165590 in contrast to H ( 0 , μ) HFISCA = 2.9042700256165865. This fact indicates that some of the values of the H -function generated by HFUNIV are correct only to the 13-th decimal place, possessing the accuracy of merely 14 significant digits. The value of δ exhibited by H HFUNIV is also somewhat larger than that for H HFISCA , and seems to attest the higher accuracy of our calculations.
Lastly, let us make the same consistency check on H ( 0 , μ) HFUNIV that we have applied to the results of the code HFISCA using Eqs. (24a), (24b). The maximum deviation ε IE and the root mean square deviation δ IE given by Eqs. (25a), (25b) and (25c) with H HFISCA replaced with H HFUNIV are shown in columns 5 and 6 respectively in the row (a) together with the corresponding results for H HFISCA in the row (b). It is interesting to see that they are virtually the same as those for ε and δ, and that the values of H HFUNIV are significantly larger than those for H HFISCA .
Calculations of ε, δ, ε IE and δ IE and measurements of the CPU time for HFUNIV and HFISCA have therefore been repeated on a new rectangular domain (Domain C) formed by (0.85 0,i , μ j ) (i,j = 1, 2, . . . , 222) with 0,i and μ j being specified by Eqs. (20a), (20b). The results are shown parenthesized in Table 3. We see the values of all the quantities are reduced in comparison with those obtained for the case with the full parameter ranges of 0 , viz., 0 ≤ 0 ≤ 1. Of particular interest is the fact that the values of ε and ε IE Fig. 5 Maximum deviation of H ( 0 , μ) from its iterate H ( 0 , μ) new as a function of 0 (see Eq. (25a)). The dotted line "a" is for H ( 0 , μ) HFUNIV and the solid line "b" is for H ( 0 , μ) HFISCA for HFUNIV are now both 1.332 × 10 −15 , a reduction by a factor of slightly more than 20. We seem to be justified to conclude that the numerical accuracy of H ( 0 , μ) HFUNIV for the parameter range 0.85 ≤ 0 ≤ 1 is not as high as that of H ( 0 , μ) HFISCA .

Calculations of the moments of the H -function
Since we have established an efficient scheme for evaluating the H -function for isotropic scattering with intended accuracy and executional speed, it must be of interest to investigate if we can calculate in a similar manner with a comparable accuracy the numerical values of the moments of the H -function that are defined as (see, e.g., Ivanov 1973;van de Hulst 1980): where we have employed a new variable s = μ n for n ≥ 1. The quantity α * −1 ( 0 ), whose exact value is equal to 2 ln H ( 0 , 1) (Chamberlain and McElroy 1966;Ivanov 1973;van de Hulst 1980), is the convergent part of the minus first moment of H ( 0 , μ). Our principal task here is to find the values for the parameters N T and T max that can minimize the maximum absolute error ε max α 0 (N T , T max ) of the 0-th moment α 0 calculated for the 222 values of 0 specified by Eq. (20b): where α 0 ( 0 ) exa is the exact value of α 0 ( 0 ) given by which can be reduced to for small values of 0 . We switch from Eq. (30a) to Eq. (30b) for 0 ≤ 3.3 × 10 −3 to increase numerical accuracy.
We have made a search for the best combination of N T and T max by varying their values in the following way: This time, we need to check the relative errors instead of absolute errors in the calculated values of the moments to see how many significant figures they are correct to, since their magnitudes vary over wide ranges: in the case of, e.g., α * −1 ( 0 ), its value goes from 0 to 2.1348 · · · . Hence, we have made calculations of the maximum relative errors of α * −1 ( 0 ) as well as α n ( 0 ) (0 ≤ n ≤ 2000) taking the 222 values of 0 given by Eq. (20b).
The maximum relative error of α * −1 is found to be 8.22 × 10 −15 and that of α 0 is 4.66 × 10 −15 . This fact indicates that the values of α * −1 should be correct to at least 13 and most probably 14 significant digits, whereas those of α 0 should be correct to at least 14 and most probably 15 significant digits. On the other hand, the maximum relative errors of the values of the higher order moments α n ( 0 ) (n ≥ 1) have turned out to be roughly an order of magnitude smaller and fall in between 4.44 × 10 −16 and 9.99 × 10 −16 , which implies that also the values of α n are to be correct at least to 14 and most probably 15 significant digits. To further confirm the findings described above, we have extended the moment calculations up to n = 5000. All of the differences in the presently calculated values of the moments in comparison with the corresponding reference values are found to be confined exclusively in the 15-th significant digits: of 5002 × 222 values of the moments calculated by HFMOMENT, 1018654(91.7%) of them are in agreement with the reference values to entire 15 digits( 15 = 0), 91784(8.3%) of them have a one-unit deviation ( 15 = ±1), three of them have a four-unit deviation ( 15 = ±4), and the remaining three of them have a five-unit deviation ( 15 = ±5). Furthermore, both 15 = ±5 and ±4 are caused solely by α * −1 , and the ordinary moments α n ( 0 ) (n ≥ 0) have at most 15 = ±1.
In Table 4, we present the values of α * −1 as well as α 0 through α 6 calculated by HFMOMENT for 40 argument values of 0 between 0 and 1. The maximum relative error (not in %) found for each degree of the moments obtained for the 222 values of 0 generated by Eq. (20b)(not the 40 argument values) is given in the bottom row of Table 4. It should be reminded again that (28b) instead of Eq. (28a) is employed in HFMOMENT to calculate the values of α n ( 0 ) for n ≥ 1 for better numerical accuracy, whereas Eqs. (28c) and (28d) are used for α 0 and α * −1 , respectively. We presume that the tabulated values of α n (n ≥ 0) should be correct to at least 14 and mostly 15 significant digits with an uncertainty of ±1 in the 15-th digits. The values of α * −1 ( 0 ) shown in column 2 of Table 4 may be viewed as an extension of Table 9 of Ivanov (1973) (see also Sect. 8.3.3 of van de Hulst 1980), although their accuracies are somewhat lower having a maximum uncertainty of ±5 in the 15-th digits.

Conclusion
In the present work, we have developed a fast and sufficiently compact double exponential formula (DE) algorithm to calculate with 15-digit accuracy (correct to the 14-th decimal place) the values of the μ) for isotropic scattering for the full parameter ranges, viz., 0 ≤ 0 ≤ 1 and 0 ≤ μ ≤ 1, where 0 and μ are respectively the single scattering albedo and an angular variable. A major improvement over the similar DE algorithm employed by Kawabata (2016) is twofold: (1) a new formulation of the DE-method (see Eqs. (10) and (11)) so as to avoid the use of the semi-analytical integration approximation proposed by  and later by Kawabata (2016) to take care of a singular behavior exhibited by the integrand function at the origin of μ when 0 = 1, and (2) elimination of the procedure for automatic adjustment of the step-size (AAI) for numerical integration of the integrand involved in Eq. (13) of the H -function. The Fortran90 subroutine subprogram HFISCA developed based on this algorithm has turned out to be faster than those so far available.
To achieve our objective (2), we have searched for numerical values of the two essential parameters, viz., T max , the upper bound or a cut-off point to apply to an infinite range integral, and N T , the number of divisions of a closed interval [0, T max ], in such a way that the maximum absolute error, the root mean square error of the resulting values of the H HFISCA and similar quantities of the newly obtained iterate H new arising when the values of H HFISCA are substituted to a defining integral equation Eq. (24a) of the H -function as an approximation H old . In so doing, we have found T max = 3.2215 and N T = 111 as the best combination, with which the maximum absolute error and the root mean square error of the values of the H -function are 1.78×10 −15 and 1.71 × 10 −16 , respectively.
As a supplementary test of the accuracy of the values of the H -function obtained by HFISCA, calculations of α * −1 ( 0 ), the convergent part of the minus first moment, as well as the n-th degree moments α n ( 0 ) (n = 0, 1, . . . , 2000) of the H -function have been performed, where the required integrals of the H -function are calculated by means of the same algorithm as that employed for the H -function evaluations. The values of T max and N T are however determined separately by inspecting the maximum absolute errors of α 0 ( 0 ), to find the combination N T = 109 and T max = 3.3834 to be the best corresponding to the maximum absolute error of 4.7 × 10 −15 . A newly developed For- tran90 subroutine subprogram HFMOMENT (also given in Appendix) implements these parameter values.
Extended calculations of the maximum relative errors of the moments for the degrees up to 5000 and a detailed digit-

Appendix: a program to calculate H ( 0 , μ) and its moments
Given below is a Fortran90 program set consisted of a main program MAIN and three function subprograms CALRSPI0, HFISCA, and HFMOMENTS. The function subprogram HFISCA enables us to evaluate, by means of the double-exponential (DE) formula developed by Takahashi and Mori (1974), the numerical values of the μ) for isotropic scattering, with a 15-digit accuracy, for a given set of values for 0 , 1 − 0 , and μ. The value of 1 − 0 may better be calculated as shown in MAIN by using the function subprogram CALRSPI0, which should help minimize the effect of the loss of significant digits in the case the value of 0 is close to unity. The function subprogram HFMOMENT produces again by employing the DE-formula the numerical values of the n-th degree moment of the H -function (n ≥ −1) defined by Eqs. (28a)-(28d) for given values of 0 , 1 − 0 , and n. It should be reminded that, for n ≥ 1, we use the integral expression Eq. (28b) to improve numerical accuracy of integration: α n ( 0 ) = 1 n 1 0 H ( 0 , s 1/n )s 1/n ds (n ≥ 1), although the usual forms of defining equations are employed for n = −1 and n = 0: (Ivanov 1973: van de Hulst 1980. The source code of the program is as shown below: It should be noted that actual output from the program MAIN shown above produces the values of H -function as well as those of its moments α * −1 ( 0 ) and α n ( 0 ) (n = 0, 1, . . . , 4) to the 15-th decimal place, although the results displayed here are rounded to the 7-th decimal place. The header "Alpha_ − 1" should be interpreted as being "α * −1 ".