当前位置: X-MOL 学术Comput. Phys. Commun. › 论文详情
Improved algorithm for calculating high accuracy values of the Chandrasekhar function
Computer Physics Communications ( IF 3.309 ) Pub Date : 2020-03-03 , DOI: 10.1016/j.cpc.2020.107237
A. Jablonski

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 ω=0.85 to ω=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. 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, ω. The integral representations of the H-function exhibit an integrable singularity for conservative scattering ( ω=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) in the vicinity of albedo equal 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≤μ≤1 and 0≤ω≤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μ,ω=exp−μπ∫0∞11+μ2x2ln1−ωarctanxxdx,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 ω values. The expression under logarithm can be rewritten as follows (2)1−ωarctanxx=1−ω+ω∑k=1∞−1k+12k+1x2k.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≤k≤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 ±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μ,ω≅1+μμμω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 is shown in Fig. 1. The border line (value of the albedo at the border, ω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 ωb at a given value of μ. 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. 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 material 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<μ≤1, and albedo, 0<ω≤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 ω<ω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≤μ≤1 and 0.01≤ω≤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<μ≤1 and 0<ω≤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 μ=10−15 and ω=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 ω=1−10−3 to ω=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.
更新日期:2020-03-03

 

全部期刊列表>>
智控未来
聚焦商业经济政治法律
跟Nature、Science文章学绘图
控制与机器人
招募海内外科研人才,上自然官网
隐藏1h前已浏览文章
课题组网站
新版X-MOL期刊搜索和高级搜索功能介绍
ACS材料视界
x-mol收录
湖南大学化学化工学院刘松
上海有机所
李旸
南方科技大学
西湖大学
X-MOL
支志明
中山大学化学工程与技术学院
试剂库存
天合科研
down
wechat
bug