Corrections to: “Accurate computation of gravitational field of a tesseroid” by Fukushima (2018) in J. Geod. 92(12):1371–1386

An accurate method with conditional split, double exponential quadrature rule, and numerical differentiation has been proposed in the paper “Accurate computation of gravitational field of a tesseroid” (Fukushima in J Geod 92(12):1371–1386, https://doi.org/10.1007/s00190-018-1126-2, 2018) to compute the gravitational field (i.e. gravitational potential, gravitational acceleration vector, and gravity gradient tensor) of a tesseroid. This study presents the corrections for some formulas in the main paper and electronic supplementary material of Fukushima (J Geod 92(12):1371–1386, https://doi.org/10.1007/s00190-018-1126-2, 2018). Moreover, the FORTRAN subroutines gtess (or qgtess) and ggtess (or qggtess) in the original codes xtess.txt (or xqtess.txt) in double (or quadrature) precision provided by Fukushima (J Geod 92(12):1371–1386, https://doi.org/10.1007/s00190-018-1126-2, 2018) are revised. The revised parts have impacts on the calculation of these components of the gravitational acceleration vector (gΦ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{\varPhi }$$\end{document} and gΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{\varLambda }$$\end{document}) and gravity gradient tensor (ΓΦΦ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{\varPhi \varPhi }$$\end{document}, ΓΦΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{\varPhi \varLambda }$$\end{document}, ΓΦH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{\varPhi H}$$\end{document}, ΓΛΛ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{\varLambda \varLambda }$$\end{document}, ΓΛH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{\varLambda H}$$\end{document}, and ΓHH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma _{H H}$$\end{document}). The revised FORTRAN codes xtess.f90 and xqtess.f90 in double and quadrature precision are presented at the GitHub website https://github.com/xiaoledeng/xtess-xqtess. These revised FORTRAN codes can accurately compute the gravitational field of a tesseroid in double and quadrature precision no matter the computation point is located outside, near the surface of, on the surface of, or inside the tesseroid. They can be applied to calculate the gravitational field of the different layers (e.g. atmosphere, topography, crust, and mantle) of the Earth or other celestial bodies, which helps investigate the various geoscience applications, e.g. geoid determination in geodesy and gravity interpretation in geophysics.


Introduction
In the paper "Accurate computation of gravitational field of a tesseroid" (Fukushima 2018), an accurate method with conditional split, double exponential quadrature rule, and numerical differentiation has been presented for calculating the gravitational potential, gravitational acceleration vector, and gravity gradient tensor of a tesseroid, whereas the expressions for the Γ ΦΛ and Γ ΛH components of the gravity gradient tensor in Fukushima (2018) contain formal errors and need to be corrected. Meanwhile, some typos of the negative sign need to be corrected. The two FORTRAN codes xtess.txt and xqtess.txt provided by Fukushima (2018) to calculate the gravitational acceleration vector and B Xiao-Le Deng xiaole.deng@gis.uni-stuttgart.de; xldeng@whu.edu.cn 1 Institute of Geodesy, University of Stuttgart, Stuttgart 70174, Germany gravity gradient tensor of a tesseroid need to be revised in double and quadruple precision.
This study presents the corrections for some formulas and two FORTRAN codes in Fukushima (2018). Furthermore, a numerical experiment is performed to reveal the consequences of the modified formulas and FORTRAN codes.
This study is organized as follows. Sections 2 and 3 present the corrections to some formulas and FORTRAN codes in Fukushima (2018), respectively. The corrected formulas and FORTRAN codes are numerically investigated against analytical solutions for a spherical shell in Sect. 4. Section 5 offers the revised FORTRAN codes xtess.f90 and xqtess.f90. Conclusions and consequences for Fukushima (2018) are summarized in Sect. 6. 8 Page 2 of 7 X.-L. Deng 2 Corrections to some formulas in Fukushima (2018) Regarding the expressions of the diagonal component (∂ 2 V /∂ H 2 ) Φ,Λ of the second-order partial derivatives in Fukushima (2018), the where the sign is wrong, and it is a typo. Equation (25) should be changed to: Similarly, for the expressions of the non-diagonal component [∂ 2 V /(∂Φ∂Λ)] H of the second-order partial derivatives in Eq. (26) of Fukushima (2018) (26) should be changed to: Regarding the expressions for the Γ ΦΛ component of the gravity gradient tensor in Eq. (31) of Fukushima (2018), the + g Φ tan Φ R 2 cos Φ term is wrong and should be corrected to the + tan Φ R g Λ term. Equation (31) should be changed to: Similarly, for the formula of the Γ ΛH component of the gravity gradient tensor in Fukushima (2018), the − g H R 2 cos Φ term in Eq. (33) is wrong and should be changed to the − g Λ R term. Then, Eq. (33) should be changed to: where the expressions in Eqs. (31) and (33) belong to the incorrect derivation. The correct derivation of the formulas in Eqs. (3) and (4) is provided in "Appendix 1". In the formula of the radial-radial component of the gravity gradient tensor Γ H H of the spherical shell when the computation point is located inside the spherical shell (H B < H < H T ) in Fukushima (2018) i.e. the negative sign '−' should be added, and it is a typo. Then, Eq. (48) should be changed to (Lin et al. 2020, Eq. (15c)): Regarding the expression for the relative error of the second-order central difference formula in Eq. (77), the i.e. the negative sign '−' should be removed, and it is a typo. Then, Eq. (77) should be changed to:

Corrections to FORTRAN codes of Fukushima (2018)
The two FORTRAN subroutines qgtess and qggtess to calculate the gravitational acceleration vector (g Φ and g Λ ) and gravity gradient tensor (Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H ) of the tesseroid in the code xqtess.txt in Tables 7-11 of the electronic supplementary material of Fukushima (2018) need to be revised in quadruple precision. The corrected subroutines qgtess and qggtess are presented in Tables 7-11 of this paper's electronic supplementary material, where the revised parts are in bold fonts. The same parts in the code xtess.txt should be revised in double precision. Note that the FORTRAN function vtess is revised to qvtess in the subroutine qggtess, which is not in bold font. Based on the revised contents in Tables 7-11 of this paper's electronic supplementary material, these revised codes will affect the evaluation of these components g Φ , g Λ , Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H , where the calculation of the g H will not be affected.

Numerical investigations
In this section, a simple numerical experiment is performed to reveal the magnitude of the error one makes when using the original codes of Fukushima (2018). The numerical experiment is carried out in a double-precision environment considering the computational efficiency, because quadruple precision generally requires more time than double precision. For instance, the ratio of quadruple precision to double precision is about 50.463/0.089 ≈ 567 per computation point with a single tesseroid on a desktop computer with an Intel i5-10400 CPU at 2.9 GHz and single-threaded operation in Sect. 5. Specifically, the values of the gravitational acceleration vector (g Φ and g Λ ) and gravity gradient tensor (Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H ) of a spherical shell are served as the analytical reference values for the calculated values of the discretized tesseroids forming the whole spherical shell with the original codes provided by Fukushima (2018) and revised codes.
For other components g Φ , g Λ , Γ ΦΛ , Γ Φ H , and Γ ΛH , the absolute errors in log 10 scale are presented as: When the computation point is located below (H < H B ) the spherical shell for the components Γ ΦΦ , Γ ΛΛ , and Γ H H , the absolute errors in log 10 scale are applied.
Regarding the components Γ ΦΦ , Γ ΛΛ , and Γ H H when the computation point is located inside (H B < H < H T ) and above (H T < H ) the spherical shell, their relative errors in log 10 scale are presented as: where F = Γ ΦΦ , Γ ΛΛ , or Γ H H . means the sum of the calculated values of the discretized tesseroids forming the whole spherical shell. F analytical represents the analytical reference values of the spherical shell. The absolute errors and relative errors in log 10 scale of the components g Φ , g Λ , Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H are shown in Fig. 1 using the revised FORTRAN codes and in Fig. 2 using the original FORTRAN codes provided by Fukushima (2018). The factors and units of the absolute errors of the gravitational acceleration vector (δg Φ and δg Λ ) and gravity gradient tensor (δΓ ΦΦ , δΓ ΦΛ , δΓ Φ H , δΓ ΛΛ , δΓ ΛH , and δΓ H H ) are m −1 and m −2 , respectively. When evaluating the practical results of the gravitational acceleration vector and gravity gradient tensor using the original or revised FORTRAN codes, the term Gρ R 2 0 (i.e. the units of G, ρ, and R 0 are m 3 kg −1 s −2 , kg m −3 , and m) should be multiplied.
In Fig. 1a, Fukushima (2018), other parameters are the same as in Fig. 1 gravitational potential. Specifically, the relative accuracy of the partially differentiated quantities is √ δ for the gravity gradient tensor (Fukushima 2018), where δ is the error tolerance to compute the gravitational potential.
In Fig. 1b, the absolute errors in log 10 scale of the δg Φ , δg Λ , δΓ ΦΛ , δΓ Φ H , and δΓ ΛH are mostly less than or equal to 10 −16 , which are the random errors in double precision. Finally, numerical results reveal that the evaluation of the absolute errors and relative errors in log 10 scale of the δΓ H H does not be affected by using the original codes, because one of the conditional splits (i.e. if Φ S ≤ Φ ≤ Φ N and in the revised code of the Γ H H is not triggered. When using the original codes for the δΓ ΦΦ and δΓ ΛΛ , their absolute and relative errors in log 10 scale were incorrectly improved by about 2-4 orders of magnitude. Regarding the δΓ ΛH , the precision of its absolute errors in log 10 scale by using the original codes is erroneously improved by approximately 3-6 orders of magnitude in the range of [−100 km, −40 km) and reduced by about 4-8 orders of magnitude in the range of [−40 km, +100 km]. The precision of the absolute errors in log 10 scale of the δΓ Φ H is reduced by about 9-11 orders of magnitude by using the original codes. Notably, the absolute errors in log 10 scale of the δg Φ , δg Λ , δΓ ΦΛ and the majority of δΓ Φ H are not presented in Fig. 2b because of no numerical values. This is due to the fact that these values equal negative infinity. In other words, the calculated values of the δg Φ , δg Λ , δΓ ΦΛ and the most of δΓ Φ H when using the original codes in Fukushima (2018) are equal to zero and their absolute errors in log 10 scale are log 10 0 = −∞.

Computer programs of the xtess.f90
and xqtess.f90 To make better use of the revised FORTRAN codes, the xtess.f90 and xqtess.f90 are presented at the GitHub website https://github.com/xiaoledeng/xtess-xqtess. The total CPU time of the xtess.f90 and xqtess.f90 to calculate all 10 components of the V , g Φ , g Λ , g H , Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H is 3.568/40 s ≈ 0.089 s per computation point when δ = 10 −16 in double precision and 50.463 s per computation point when δ = 10 −33 in quadrature precision with a single tesseroid mass body. These values are obtained on a desktop computer with an Intel i5-10400 CPU at 2.9 GHz using the single-threaded operation. When applying the gravitational potential, gravitational acceleration vector, and gravity gradient tensor using these FORTRAN codes in the practical applications, the term Gρ R 2 0 should be multiplied with the output numerical values to obtain the units m 2 s −2 for the gravitational potential, m s −2 for the gravitational acceleration vector, and s −2 for the gravity gradient tensor, where the units of G, ρ, and R 0 are m 3 kg −1 s −2 , kg m −3 , and m, respectively.

Conclusions and consequences for Fukushima (2018)
Theoretically, the revised parts in the original codes xtess.txt and xqtess.txt have impacts on the calculation of these components of the gravitational acceleration vector (i.e. g Φ and g Λ ) and gravity gradient tensor (i.e. Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , Γ ΛH , and Γ H H ). Regarding the Γ H H , the influence on the calculation results only occurs when one of the conditional splits (i.e. if Φ S ≤ Φ ≤ Φ N and The evaluation of the gravitational potential V and the radial component of the gravitational acceleration vector g H in the original codes xtess.txt and xqtess.txt will not be affected. Numerical results reveal that the absolute errors obtained by using the corrected codes are at a lower precision level by 2-4 orders of magnitude for the δΓ ΦΦ and δΓ ΛΛ components of the gravity gradient tensor. Regarding the δΓ ΛH component of the gravity gradient tensor, the precision level of the absolute errors by using the corrected codes is reduced by 3-6 orders of magnitude when the computation point is located below the spherical shell and improved by 4-8 orders of magnitude when the computation point is located in and above the spherical shell. When the computation point is located on the bottom and top boundaries of the spherical shell, the precision level of the absolute errors for the δΓ Φ H component of the gravity gradient tensor by using the corrected codes is improved by 9-11 orders of magnitude. If using the original codes for the g Φ and g Λ components of the gravitational acceleration vector and the δΓ ΦΛ and most of the δΓ Φ H components of the gravity gradient tensor, their calculated values of the tesseroids are equal to zero. When using the revised codes for these components, their absolute errors in log 10 scale are in ranges of about [−23, −21] for the δg Φ , [−21, −16] for the δg Λ , [−27, −25] for the δΓ ΦΛ , and [−27, −25] for the δΓ Φ H .
The previous study that quoted the original FORTRAN codes provided by Fukushima (2018) to calculate the gravitational field of a tesseroid should be carefully considered based on the above potential impacts. For example, Lin et al. (2020) calculated the gravitational potential V , radial component V z of the gravitational acceleration vector, and radial-radial component V zz of the gravity gradient tensor in Sect. 3.7 based on the use of the original FORTRAN codes in Fukushima (2018). Fortunately, these three components V , V z , and V zz are not affected, although the FORTRAN code of the V zz needs to be modified, whereas in Sect. 4.2 gravitational acceleration vector of tesseroid, Sect. 4.3 gravity gradient tensor of tesseroid, Sect. 4.4 polar tesseroid, and Sect. 4.5 polar cap slab of the electronic supplement material of Fukushima (2018), the evaluation of the gravity field quantities (e.g. the total gravitational acceleration g = g 2 Φ + g 2 Λ + g 2 H , magnitude of the vector representing the deflection of the vertical referred to the normal vector of the reference sphere θ = tan −1 angle of the vector representing the deflection of the vertical g Λ , and these components Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , and Γ ΛH of the gravity gradient tensor) may be wrong and needs to be carefully investigated especially when these components of the gravitational acceleration vector (i.e. g Φ and g Λ ) and gravity gradient tensor (i.e. Γ ΦΦ , Γ ΦΛ , Γ Φ H , Γ ΛΛ , and Γ ΛH ) are included.
Regarding the impacts of replacing the revised codes in this study with the original codes in Fukushima (2018) when considering practical applications for modeling the gravitational signals of mass distributions of the Earth or other planetary bodies, further empirical research is required to explore these in the future.
where p, q, s = 1, 2, 3. Note that the Einstein summation convention has been applied in Eq. (16). Γ s qp is the Christoffel's symbol of the second kind, which can be referred to in Eqs. (39)-(41) of Casotto and Fantino (2009) and Table 2 of Deng and Ran (2021).
The Γ ΦΛ component of the gravity gradient tensor is obtained with p = 2 and q = 1 in Eq. (16) as: where Γ 1 12 = − tan Φ and Γ 2 12 = Γ 3 12 = 0 are referred to in Table 2 of Deng and Ran (2021). The equation 1 (15) is applied. Similarly, the expression for the Γ ΛΦ can be obtained with p = 1 and q = 2 in Eq. (16) and is the same as the Γ ΦΛ .
Furthermore, the Γ ΛH component of the gravity gradient tensor is obtained with p = 1 and q = 3 in Eq. (16) as: where Γ 1 31 = 1/R and Γ 2 31 = Γ 3 31 = 0 are referred to in Table 2 of Deng and Ran (2021). The equation 1 R cos Φ ∂ V ∂Λ = g Λ in Eq. (15) is applied. Analogously, the expression for the Γ H Λ can be obtained with p = 3 and q = 1 in Eq. (16) and is the same as the Γ ΛH .
Note that when the subscripts H and Φ are added for the first terms in Eqs. (17) and (18), Eqs. (17) and (18) are the same as Eqs. (3) and (4).