The recently published set of programs H_FUN for calculation of the Chandrasekhar function [Comput. Phys. Commun. 235 (2019) 489-501] is replaced with a new set H_FUN_v2 which has the following improvements: (i) increased accuracy from 14 decimals to 15–16 decimals, (ii) the upper limit of albedo values increased from $\omega =0.85$ to $\omega =1$, (iii) a reasonably short execution time, and (iv) brevity and simplicity of the code. The previously edited program HFUNELE dedicated to applications in electron transport theory is replaced with a short universal program HFUNIV that allows to reach the desired accuracy for small values of arguments, for albedo values close to unity, and also for conservative scattering. This program turns out to be very compact with listing close to one page (96 Fortran lines). Additionally, a new test program is added that illustrates the procedure for determination of accuracy by comparison with the enclosed extensive database containing the Chandrasekhar function reference values determined with accuracy of 21 decimals.

New version program summary

Program Title: H_FUN_v2

Program Files doi: http://dx.doi.org/10.17632/xf4k25cnpk.2

Licensing provisions: GNU General Public License 3

Programming language: Fortran 90

Journal reference of previous version: Comput. Phys. Commun. 235 (2019) 489-501.

Does the new version supersede the previous version?: Yes

Reasons for the new version:

In recent literature, a number of publications have addressed methods for reaching high accuracy of Chandrasekhar function (H-function) values for isotropic scattering [1-4]. Algorithms leading to an accuracy of 14 or 15 decimals were proposed in these reports. It was proved that implementation of the double-exponential rule (DE-rule) in these algorithms considerably improves the execution time. The recently published program HFUNELE based partially on the DE-rule was dedicated to applications in theoretical models of electron transport in condensed matter [4]. This algorithm ensured an accuracy of at least 14 decimals for arguments varying in ranges needed in these applications, i.e., from zero to unity for the directional variable, x, and from zero to 0.85 for the albedo, $\omega$.

The integral representations of the H-function exhibit an integrable singularity for conservative scattering ( $\omega =1$) which typically requires an additional procedure for its removal. It was presently found that a suitable transformation of the relevant integrands for application of the DE-rule also removes this singularity; no additional procedures are needed. It effect a short universal program for the H-function was edited that allows to reach the desired accuracy of at least 15 decimals: (i) for small values of arguments, (ii) for albedo values close to unity, and (iii) for conservative scattering. The presently developed universal program HFUNIV turned out to be very compact with listing of 122 lines (96 lines excluding comments), in contrast with 868 lines of previously published program HFUNELE. Nonetheless, the H-function values are calculated with accuracy of at least 15 decimals for full range of arguments, i.e., $0\le \mu \le 1$ and $0\le \omega \le 1$.

Summary of revisions:

(1) The DE-rule was applied to integrals in different integral representations of the H-function. Among these representations, the expression reported by Ivanov (ref. [5], p. 128), (1)$H\left(\mu ,\omega \right)=exp\left[\frac{-\mu }{\pi }{\int }_0^{\infty }\frac{1}{1+{\mu }^2{x}^2}ln\left(1-\omega \frac{\arctan x}{x}\right)dx\right],$ turned out to be the most convenient to use in terms of computing time, accuracy, and ranges of arguments. In calculations of the integrand in Eq. (1) a modification was added now which is active in the case of small integration arguments x and large $\omega$ values. The expression under logarithm can be rewritten as follows (2)$1-\omega \frac{\arctan x}{x}=\left(1-\omega \right)+\omega \sum _{k=1}^{\infty }\frac{{\left(-1\right)}^{k+1}}{2k+1}{x}^{2k}.$ For integration variable, $x<0.15$, it is recommended to use the series representation given by Eq. (2) truncated to eleven terms, i.e., for $1\le k\le 11$.

(2) Accuracy of 15 or 16 decimals is close to the number of decimals in data representation in double precision arithmetic. A major problem is then to precisely evaluate an accuracy of the H-function values. For this purpose, much effort was devoted to create a reference set of high precision H-function values (21 decimals) for arguments linearly or logarithmically distributed in presently considered ranges (in total, 49562 values). These values were calculated in the quadruple precision arithmetic, i.e., 128 bit representation of real variables, functions and subroutines. This is equivalent to precision of 33 significant digits, and the exponent close to $\pm 4000$. The procedure used and the exemplary values were presented in earlier publications [4, 6]. A subset of reference data (10000 values) is presently included into the test programs as it may be of interest to use for evaluation of other algorithms providing the H-function.

(3) For small albedo values, the H-function can be calculated with accuracy of 15 or 16 decimals from an approximate expression [4, 6] (3)$H\left(\mu ,\omega \right)\cong {\left(\frac{1+\mu }{\mu }\right)}^{\frac{\mu \omega }{2}}.$ To evaluate the accuracy of Eq. (3) and determine the range of validity, the set of reference H-function values for dense grid of logarithmically distributed arguments was created with arguments varying from ${10}^{-15}$ to unity (151 × 151=22801 values). The pairs of arguments for which the ultimate accuracy of 16 decimals was reached are shown in Fig. 1. The border line (value of the albedo at the border, ${\omega }_b,$ as a function of directional variable, x) is also shown in this figure; thus, it is recommended to use Eq. (3) for albedo values smaller than ${\omega }_b$ at a given value of $\mu$. It should be noted that the number of decimals referred to here as an accuracy is taking into account all digits of a given H-function value, including the integer part. Details are discussed in earlier works [4, 6].

(4) A new Fortran procedure HFDE for calculating the H-function based on Eq. (1) was presently developed. It is implementing the DE-rule based on the variable transformation recommended in Ref. [4]. As mentioned, Eq. (2) was used in calculations of the integrand for integration variable values $x<0.15$. Furthermore, care was taken to prevent the loss of accuracy in the region of albedo values close to unity. This action is realized within a few lines of the Fortran code:

 IF (OMEGA .GT. 0.999D0) THEN

 ARGUM=DNINT((1.0D0-OMEGA)*1.0D15)*1.0D-15+OMEGA*ARGUM

 ELSE

 ARGUM=1.0D0-OMEGA+OMEGA*ARGUM

 ENDIF

 B=DLOG(ARGUM)

where the variable ARGUM initially contains a particular value of the ratio $\arctan \left(x\right)∕x$. Accuracy of the H-function values calculated from subroutine HFDE for 22801 pairs of arguments is shown in Fig. 2.

Nature of problem: In some applications of the theory of electron transport in solids to problems of surface analysis by electron spectroscopies, we need values of the H-function calculated with high accuracy. A further need is a reasonably short execution time of a given program calculating the H-function due to the fact that this function may be frequently called in computational algorithms implementing the theory of electron transport. Execution times of the order of single microseconds accompanied by accuracy of 15 decimals would improve performance of these algorithms. Since the albedo values for some materials may in general exceed 0.85, the high accuracy of the H-function should be reached in a full range of arguments, i.e., directional variable, $0<\mu \le 1$, and albedo, $0<\omega \le 1$.

Solution method: Two different integral representations of the H-function were selected for analysis. On suitable transformation of the corresponding integrands, the DE-rule was applied in both cases. As recommended for this approach, the trapezoidal rule was used for integration. This simple quadrature was sufficiently accurate to ensure the desired performance of both considered integral representations. However, the algorithm implementing Eq. (1) turned out to have shorter execution time and shorter listing. For this reason it is recommended for use. As shown in Fig. 1, the precision of at least 15 decimals is also reached by the approximate analytical formula in the range of small albedo values. Consequently, a mixed algorithm, HFUNIV, was created to ensure additional decrease of the execution time. In this algorithm, the approximate formula is used for $\omega <{\omega }_b$ while the algorithm HFDE is used in remaining cases. This algorithm is finally recommended for general use.

The procedure for evaluation of accuracy is the same as in earlier publications [4, 6]. To illustrate this procedure, a short test program HFUNIV_NO_DECIMALS is presently added. The evaluation is performed for arguments varying in ranges $0.01\le \mu \le 1$ and $0.01\le \omega \le 1$ in steps of 0.01 (10000 H-function values). The second program, TEST_PROGRAM_v2, is the same as in ref. [4], except the calculations are performed using subroutine HFUNIV.

Additional comments:

(i) Restrictions: The acceptable arguments of the Chandrasekhar function are restricted to ranges of frequent interest: $0<\mu \le 1$ and $0<\omega \le 1$. No additional restrictions are needed. The program HFUNIV provides accurate values of the H-function also in the range of small arguments down to $\mu ={10}^{-15}$ and $\omega ={10}^{-15}$.

(ii) Unusual features: In the integral representations of the H-function, a singularity appears when the albedo approaches unity. In the past, different procedures of singularity removal were proposed to obtain the nonsingular integrand. However, it was presently found that the transformation of integrands recommended for application of the DE-rule does not require any further transformation to remove singularity. The DE-rule accompanied with trapezoidal quadrature is also convergent for the conservative scattering.

Furthermore, the program HFDE is convergent in the vicinity of conservative scattering. It was found that accuracy of 15 or more decimals was reached for albedo values varying from $\omega =1-{10}^{-3}$ to $\omega =1-{10}^{-14}$. These features resulted in brevity of the code HFDE implementing the DE-rule; the listing of this program has only 74 lines. Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] N. Mohankumar and A. Rawat, Some comments on the evaluation of the Chandrasekhar H function, Comput. Phys. Commun. 184 (2013) 579-581. K. Kawabata, 15-digit accuracy calculations of Chandrasekhar’s H-function for isotropic scattering by means of the double exponential formula, Astrophys. Space Sci. 361 (2016) 373. K. Kawabata, 15-digit accuracy calculations of Ambartsumian-Chandrasekhar’s H-functions for four-term phase functions with the double-exponential formula, Astrophys. Space Sci. 363 (2018) 1. A. Jablonski, The Chandrasekhar function for modeling photoelectron transport in solids, Comput. Phys. Commun. 235 (2019) 489-501. V. V. Ivanov, Transfer of Radiation in Spectral Lines, National Bureau of Standards Special Publication 385, National Bureau of Standards, Washington D.C., 1973 (English edition of Radiative Transfer and the Spectra of Celestial Bodies, Nauka, Moskva, 1969). [6] A. Jablonski, The Chandrasekhar function revisited, Comput. Phys. Commun. 196 (2015) 416-428.

中文翻译：

---

用于计算Chandrasekhar函数的高精度值的改进算法

最近发布的程序集H_FUN，用于计算Chandrasekhar函数[Comput。物理 公社 235（2019）489-501]替换为具有以下改进的新集合H_FUN_v2：（i）将精度从14位小数提高到15-16位小数，（ii）反照率值的上限从$\omega =0。85$ 至 $\omega =\mathrm{1个}$，（iii）相当短的执行时间，以及（iv）代码的简洁和简洁。先前编辑的专门用于电子传输理论的程序HFUNELE被替换为简短的通用程序HFUNIV，该程序可以在较小的自变量值，接近于1的反照率值以及保守散射方面达到所需的精度。该程序非常紧凑，只列出了一页（96个Fortran行）。此外，添加了一个新的测试程序，该程序说明了与随附的包含Chandrasekhar函数参考值（精确度为21位小数）的广泛数据库进行比较而确定精度的过程。

新版本程序摘要

程式名称： H_FUN_v2

程序文件doi： http : //dx.doi.org/10.17632/xf4k25cnpk.2

许可条款： GNU通用公共许可证3

编程语言： Fortran 90

先前版本的期刊参考：计算。物理 社区 235（2019）489-501。

新版本会取代旧版本吗？：是

新版本的原因：

在最近的文献中，许多出版物已经提出了用于达到各向同性散射的Chandrasekhar函数（H函数）值的高精度方法[1-4]。这些报告中提出了导致小数点后14位或15位精度的算法。事实证明，在这些算法中双指数规则（DE-rule）的实现大大提高了执行时间。最近发布的部分基于DE规则的程序HFUNELE致力于在凝聚态电子传输理论模型中的应用[4]。对于在这些应用程序所需的范围内变化的参数，该算法确保了至少14位小数的精度，即方向变量x从零到1，反照率从零到0.85，$\omega$。

H函数的积分表示对于保守散射显示出可积分的奇点（$\omega =\mathrm{1个}$），通常需要执行其他删除程序。目前发现，对适用于DE规则的相关积分进行适当的转换也消除了这种奇异性。无需其他程序。它编辑了一个简短的H函数通用程序，该程序允许达到至少15个小数的所需精度：（i）小参数值，（ii）接近1的反照率值，（iii）保守散射。与以前发布的HFUNELE程序的868行相比，当前开发的通用程序HFUNIV的结构非常紧凑，列出了122行（不含注释的96行）。尽管如此，H-函数值对于所有参数范围的精度至少为15位小数，即 $0\le \mu \le \mathrm{1个}$ 和 $0\le \omega \le \mathrm{1个}$。

修订摘要：

（1）DE规则应用于H函数的不同积分表示形式的积分。在这些表示中，伊万诺夫（Ivanov）报告了该表达（参考[5]，第128页），（1）$H\left(\mu ，\omega \right)=经验值\left[\frac{-\mu }{\pi }{\int }_0^{\infty }\frac{\mathrm{1个}}{\mathrm{1个}+{\mu }^2{X}^2}ln\left(\mathrm{1个}-\omega \frac{\mathrm{Arctan}X}{X}\right)dX\right]，$ 在计算时间，准确性和参数范围方面，事实证明是最方便使用的。在等式中被积分的计算中。（1）现在添加了一个修改，该修改在小积分参数x和大积分参数的情况下有效$\omega$价值观。对数下的表达式可以重写如下（2）$\mathrm{1个}-\omega \frac{\mathrm{Arctan}X}{X}=\left(\mathrm{1个}-\omega \right)+\omega \sum _{ķ=\mathrm{1个}}^{\infty }\frac{{\left(-\mathrm{1个}\right)}^{ķ+\mathrm{1个}}}{2ķ+\mathrm{1个}}{X}^{2ķ}。$ 对于积分变量， $X<0。15$，建议使用方程式给出的序列表示。（2）截短为11个词，即$\mathrm{1个}\le ķ\le 11$。

（2）15位或16位小数的精度接近双精度算术中数据表示中的小数位数。于是，主要的问题是精确评估H函数值的准确性。为此，为在当前考虑的范围内线性或对数分布的参数（总共49562个值）创建了高精度H函数值（21个十进制小数）的参考集。这些值以四倍精度算术（即实数变量，函数和子例程的128位表示）进行计算。这相当于33个有效数字的精度，指数接近$\pm 4000$。先前的出版物[4，6]中介绍了所使用的过程和示例性值。参考数据的子集（10000个值）目前包含在测试程序中，因为可能有兴趣将其用于评估提供H函数的其他算法。

（3）为小反照率值，则ħ -function可以具有15或16位小数精度计算从近似表达[4,6] （3）$H\left(\mu ，\omega \right)\cong {\left(\frac{\mathrm{1个}+\mu }{\mu }\right)}^{\frac{\mu \omega }{2}}。$ 要评估方程式的准确性。（3）并确定有效范围，创建了对数分布自变量的密集网格的参考H函数值集，其中自变量为${10}^{-15}$统一（151×151 = 22801值）。图1显示了达到16位小数的最终精度的参数对。边界线（边界处反照率的值，${\omega }_b，$作为方向变量的函数，x）也显示在该图中；因此，建议使用公式。（3）对于小于${\omega }_b$ 给定值 $\mu$。应当注意，此处称为精度的小数位数是考虑到给定H函数值的所有数字，包括整数部分。在较早的著作中讨论了细节[4，6]。

（4）一种新的Fortran程序HFDE，用于基于等式计算H函数。（1）目前已开发。它正在基于参考文献中建议的变量转换来实施DE规则。[4]。如前所述，（2）用于积分变量值的积分计算$X<0。15$。此外，要小心以防止在接近单位的反照率值区域中损失精度。在Fortran代码的几行中即可实现此操作：

 如果（OMEGA .GT.0.999D0）然后

 ARGUM = DNINT（（1.0D0-OMEGA）* 1.0D15）* 1.0D-15 + OMEGA * ARGUM

 其他

 ARGUM = 1.0D0-OMEGA + OMEGA * ARGUM

 万一

 B = DLOG（ARGUM）

其中变量ARGUM最初包含比率的特定值$\mathrm{Arctan}（X）∕X$。由子例程HFDE计算的22801对自变量的H函数值的精度如图2所示。

问题的性质：在固体电子传输理论对电子光谱学表面分析问题的某些应用中，我们需要以高精度计算出的H函数值。进一步的需求是给定程序计算H函数的执行时间相当短，这是由于在实现电子传输理论的计算算法中可能经常调用该函数。执行时间只有几微秒，并伴有15位小数的精度，可以提高这些算法的性能。由于某些材料的反照率值通常可能超过0.85，因此H的高精度-function应该在所有参数（即方向变量， $0<\mu \le \mathrm{1个}$和反照率， $0<\omega \le \mathrm{1个}$。

求解方法：H的两种不同的积分表示选择功能进行分析。在对相应积分进行合适的转换后，在两种情况下均采用了DE规则。按照此方法的建议，梯形法则用于积分。这个简单的正交足以保证两个考虑的积分表示的期望性能。但是，实现等式的算法。（1）证明执行时间更短，上市时间更短。因此，建议使用。如图1所示，在小的反照率值范围内，近似解析公式也可以达到至少15位小数的精度。因此，创建了混合算法HFUNIV，以确保进一步减少执行时间。在此算法中，近似公式用于$\omega <{\omega }_b$其余情况则使用HFDE算法。最后，建议将该算法用于一般用途。

准确性评估的程序与早期出版物[4，6]相同。为了说明此过程，目前添加了一个简短的测试程序HFUNIV_NO_DECIMALS。对范围不同的参数执行评估$0。01\le \mu \le \mathrm{1个}$ 和 $0。01\le \omega \le \mathrm{1个}$以0.01（10000 H-函数值）为单位。第二个程序TEST_PROGRAM_v2与ref中的相同。[4]，除了计算是使用子例程HFUNIV执行的。

附加评论：

（i）限制： Chandrasekhar函数的可接受参数仅限于经常引起关注的范围：$0<\mu \le \mathrm{1个}$ 和 $0<\omega \le \mathrm{1个}$。无需其他限制。该方案HFUNIV提供的精确的值ħ也是在小参数下降到范围-function$\mu ={10}^{-15}$ 和 $\omega ={10}^{-15}$。

（ii）异常特征：在H函数的积分表示中，当反照率趋于统一时出现奇异性。过去，提出了不同的去除奇异点的方法来获得非奇异的被积数。但是，目前发现，推荐用于DE规则应用的整物转换不需要任何进一步的转换即可消除奇异性。伴随梯形正交的DE规则也收敛于保守散射。

此外，程序HFDE在保守散射附近收敛。结果发现，反照率值从$\omega =\mathrm{1个}-{10}^{-3}$ 至 $\omega =\mathrm{1个}-{10}^{-14}$。这些功能使实施DE-规则的HFDE代码更加简洁。该程序的列表只有74行。竞争利益声明

作者宣称，他们没有已知的竞争财务利益或个人关系，这些关系或个人关系似乎会影响本文报道的工作。参考文献[1] N. Mohankumar和A. Rawat，对Chandrasekhar H函数的评估，计算机。物理 公社 184（2013）579-581。K. Kawabata，借助双指数公式Astrophys对Chandrasekhar的H函数进行各向同性散射的15位精度计算。太空科学 361（2016）373. K. Kawabata，Ambartsumian-Chandrasekhar的H的15位精度计算-具有双指数公式Astrophys的四项相位函数。太空科学 363（2018）1. A. Jablonski，Chandrasekhar函数，用于建模固体中的光电子传输，Comput。物理 社区 235（2019）489-501。伊万诺夫（VV Ivanov），《光谱线中辐射的转移》，国家标准局特别出版物385，国家标准局，华盛顿特区，1973年（《辐射转移和天体光谱》，莫斯科，瑙卡，英文版，1969年）。[6] A. Jablonski，Chandrasekhar函数，Comput。物理 社区 196（2015）416-428。