Precise and Fast Computation of the Gravitational Field of a General Finite Body and Its Application to the Gravitational Study of Asteroid Eros

Let us consider the gravitational potential of a general finite-body-like solar system object. It is mathematically described as a volume integral in Equation (1). Without losing generality, we assume that the body is composed of multiple layers and all the layers are well defined in a certain spherical polar coordinate system associated with the object. If this assumption is not satisfied in a single spherical polar coordinate system, such as in the case of comet 67P Churyumov–Gerasimenko (Jorda et al. 2016), one may separate the finite body into some pieces so that the radial uniqueness condition is satisfied in each piece after taking different spherical polar coordinate systems with appropriate shift of the coordinate centers.

We denote by ${R}_{{\rm{T}}}(\phi ,\lambda )$ and ${R}_{{\rm{B}}}(\phi ,\lambda )$ the radii of the top and the bottom surfaces of the layer, respectively. They are uniquely defined as functions of the latitude ϕ and the longitude λ. We also write the volume mass density of the layer as $\rho (r,\phi ,\lambda )$, where r is the radius vector. Hereafter, we focus on the gravitational field of a single layer. Within the layer, we assume that $\rho (r,\phi ,\lambda )$ is properly expanded as a low-degree polynomial of r as

$$\begin{eqnarray}&&\rho (r,\phi ,\lambda )=\displaystyle \sum _{n=0}^{N}{\rho }_{n}(\phi ,\lambda ){\left(\displaystyle \frac{r}{{R}_{0}}\right)}^{n},\end{eqnarray} \tag{ 5 }$$

where R₀ is a certain reference value of the radius and N is a certain integer, say, in the range of 0–5. In fact, the density profile of the Earth is approximated by piecewise linear or cubic polynomials (Dziewonski & Anderson 1981; Kennett 1998; Ganguly et al. 2009). The density profile of the tropospheric air mass can also be represented by a quintic polynomial (Karcol 2011). The separation of rⁿ from the full expression of the density profile, $\rho (r,\phi ,\lambda )$, is a central trick of the new method. Indeed, as we show below, the kernel function can be analytically integrable even after the inclusion of rⁿ.

The gravitational potential of the layer at an arbitrary point is written as a similar finite series as

$$\begin{eqnarray}&&V(R,{\rm{\Phi }},{\rm{\Lambda }})=\displaystyle \sum _{n=0}^{N}{V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }}),\end{eqnarray} \tag{ 6 }$$

where $(R,{\rm{\Phi }},{\rm{\Lambda }})$ are the spherical polar coordinates of the evaluation point. Each term of the potential is expressed as a volume integral in the spherical polar coordinate system as

$$\begin{eqnarray}{V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \equiv & \displaystyle \frac{G}{{R}_{0}}{\displaystyle \int }_{-\pi /2}^{\pi /2}\left[{\displaystyle \int }_{0}^{2\pi }{\rho }_{n}(\phi ,\lambda )\left({\displaystyle \int }_{{R}_{{\rm{B}}}(\phi ,\lambda )}^{{R}_{{\rm{T}}}(\phi ,\lambda )}\right.\right.\\ & & \times \left.\left.\displaystyle \frac{{r}^{2}}{S(r,\phi ,\lambda )}{\left(\displaystyle \frac{r}{{R}_{0}}\right)}^{n}{dr}\right)d\lambda \right]\cos \phi \,d\phi .\end{eqnarray} \tag{ 7 }$$

Here $S(r,\phi ,\lambda )$ is a square root representing the normalized distance between the internal and evaluation points as

$$\begin{eqnarray}&&S(r,\phi ,\lambda )\equiv \sqrt{A(\phi ,\lambda )+{[B(\phi ,\lambda )+\zeta ]}^{2}},\end{eqnarray} \tag{ 8 }$$

where $A(\phi ,\lambda )$ and $B(\phi ,\lambda )$ are a pair of auxiliary functions defined as

$$\begin{eqnarray}&&A(\phi ,\lambda )\equiv B(\phi ,\lambda )[2\alpha -B(\phi ,\lambda )],\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&B(\phi ,\lambda )\equiv 2\alpha \left[{\sin }^{2}\left(\displaystyle \frac{\xi }{2}\right)+\cos \phi \cos {\rm{\Phi }}{\sin }^{2}\left(\displaystyle \frac{\eta }{2}\right)\right],\end{eqnarray} \tag{ 10 }$$

and α, ξ, η, and ζ are non-dimensional quantities defined as

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{R}{{R}_{0}},\,\,\xi \equiv \phi -{\rm{\Phi }},\,\,\eta \equiv \lambda -{\rm{\Lambda }},\,\,\zeta \equiv \displaystyle \frac{r-R}{{R}_{0}}.\end{eqnarray} \tag{ 11 }$$

The innermost line integral of ${V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }})$, namely that with respect to r, can be analytically evaluated in a closed form as explained in Appendix A. As a result, ${V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }})$ is rewritten as a surface integral as

$$\begin{eqnarray}&&{V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }})={{GR}}_{0}^{2}{\int }_{-\pi /2}^{\pi /2}\left({\int }_{0}^{2\pi }{\rho }_{n}(\phi ,\lambda ){K}_{n}(\phi ,\lambda )d\lambda \right)\cos \phi d\phi ,\end{eqnarray} \tag{ 12 }$$

where ${K}_{n}(\phi ,\lambda )$ is the kernel function of degree n defined as a difference of integrals as

$$\begin{eqnarray}&&{K}_{n}(\phi ,\lambda )\\ &&\equiv \,{I}_{n}({R}_{{\rm{T}}}(\phi ,\lambda ),\phi ,\lambda )-{I}_{n}({R}_{{\rm{B}}}(\phi ,\lambda ),\phi ,\lambda ),\end{eqnarray} \tag{ 13 }$$

while ${I}_{n}(r,\phi ,\lambda )$ is an indefinite line integral defined and analytically expressed as

$$\begin{eqnarray}{I}_{n}(r,\phi ,\lambda ) & \equiv & \displaystyle \frac{1}{{R}_{0}^{n+3}}{\displaystyle \int }^{r}\displaystyle \frac{{(r^{\prime} )}^{n+2}{dr}^{\prime} }{S(r^{\prime} ,\phi ,\lambda )}={C}_{n}(\phi ,\lambda )L(r,\phi ,\lambda )\\ & & +{D}_{n}(r,\phi ,\lambda )S(r,\phi ,\lambda ).\end{eqnarray} \tag{ 14 }$$

Here $L(r,\phi ,\lambda )$ is a logarithmic function conditionally defined as

$$\begin{eqnarray}&&L(r,\phi ,\lambda )\\ \,\equiv \left\{\begin{array}{ll}\mathrm{ln}[S(r,\phi ,\lambda )+B(\phi ,\lambda )+\zeta ], & (B(\phi ,\lambda )+\zeta \geqslant 0),\\ \mathrm{ln}[A(\phi ,\lambda )/\{S(r,\phi ,\lambda ) & \\ \quad -B(\phi ,\lambda )-\zeta \}], & (B(\phi ,\lambda )+\zeta \lt 0),\end{array}\right.\end{eqnarray} \tag{ 15 }$$

and ${C}_{n}(\phi ,\lambda )$ and ${D}_{n}(r,\phi ,\lambda )$ are the coefficient functions of degree n defined as polynomials of $B(\phi ,\lambda )$, α, and ζ. The first four pairs are explicitly displayed here as

$$\begin{eqnarray}&&{C}_{0}\equiv {\alpha }^{2}-3\alpha B+\displaystyle \frac{3}{2}{B}^{2},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{D}_{0}\equiv \left(2\alpha -\displaystyle \frac{3}{2}B\right)+\displaystyle \frac{1}{2}\zeta ,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{C}_{1}\equiv (\alpha -B)\left({\alpha }^{2}-5\alpha B+\displaystyle \frac{5}{2}{B}^{2}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{D}_{1}\equiv \left(3{\alpha }^{2}-\displaystyle \frac{35}{6}\alpha B+\displaystyle \frac{5}{2}{B}^{2}\right)+\left(\displaystyle \frac{3}{2}\alpha -\displaystyle \frac{5}{6}B\right)\zeta +\displaystyle \frac{1}{3}{\zeta }^{2},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{C}_{2}\equiv {\alpha }^{4}-10{\alpha }^{3}B+\displaystyle \frac{45}{2}{\alpha }^{2}{B}^{2}-\displaystyle \frac{35}{2}\alpha {B}^{3}+\displaystyle \frac{35}{8}{B}^{4},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{D}_{2}\equiv \left(4{\alpha }^{3}-\displaystyle \frac{43}{3}{\alpha }^{2}B+\displaystyle \frac{175}{12}\alpha {B}^{2}-\displaystyle \frac{35}{8}{B}^{3}\right)\\ &&\quad +\,\left(3{\alpha }^{2}-\displaystyle \frac{49}{12}\alpha B+\displaystyle \frac{35}{24}{B}^{2}\right)\zeta +\left(\displaystyle \frac{4}{3}\alpha -\displaystyle \frac{7}{12}B\right){\zeta }^{2}+\displaystyle \frac{1}{4}{\zeta }^{3},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{C}_{3}\equiv (\alpha -B)\left({\alpha }^{4}-14{\alpha }^{3}B+\displaystyle \frac{77}{2}{\alpha }^{2}{B}^{2}-\displaystyle \frac{63}{2}\alpha {B}^{3}+\displaystyle \frac{63}{8}{B}^{4}\right),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}{D}_{3} & \equiv & \left(5{\alpha }^{4}-\displaystyle \frac{85}{3}{\alpha }^{3}B+\displaystyle \frac{1001}{20}{\alpha }^{2}{B}^{2}-\displaystyle \frac{273}{8}\alpha {B}^{3}+\displaystyle \frac{63}{8}{B}^{4}\right)\\ & & +\left(5{\alpha }^{3}-\displaystyle \frac{145}{12}{\alpha }^{2}B+\displaystyle \frac{399}{40}\alpha {B}^{2}-\displaystyle \frac{21}{8}{B}^{3}\right)\zeta \\ & & +\left(\displaystyle \frac{10}{3}{\alpha }^{2}-\displaystyle \frac{69}{20}\alpha B+\displaystyle \frac{21}{20}{B}^{2}\right){\zeta }^{2}\\ & & +\left(\displaystyle \frac{5}{4}\alpha -\displaystyle \frac{9}{20}{B}^{2}\right){\zeta }^{3}+\displaystyle \frac{1}{5}{\zeta }^{4},\end{eqnarray} \tag{ 23 }$$

where we dropped the argument dependence for simplicity. Refer to Appendix A for a recursive algorithm to obtain the coefficient functions of arbitrarily high degrees. The correctness of the analytical expression of ${I}_{n}(r,\phi ,\lambda )$, Equation (14), is easily confirmed by its direct differentiation with respect to r.

Originally, the functions $S(r,\phi ,\lambda )$, $A(\phi ,\lambda )$, and $B(\phi ,\lambda )$ are expressed as

$$\begin{eqnarray}&&S(r,\phi ,\lambda )=\displaystyle \frac{\sqrt{{R}^{2}-2{rR}\cos \psi +{r}^{2}}}{{R}_{0}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&A(\phi ,\lambda )={\alpha }^{2}{\sin }^{2}\psi ,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&B(\phi ,\lambda )=2\alpha {\sin }^{2}\displaystyle \frac{\psi }{2},\end{eqnarray} \tag{ 26 }$$

where ψ is the spherical angle between the direction vectors toward the internal and the evaluation points from the coordinate center, such that

$$\begin{eqnarray}&&\cos \psi =\sin \phi \sin {\rm{\Phi }}+\cos \phi \cos {\rm{\Phi }}\cos (\lambda -{\rm{\Lambda }}).\end{eqnarray} \tag{ 27 }$$

Nevertheless, the computation of ψ from this conventional definition suffers from a cancellation problem when $\phi \approx {\rm{\Phi }}$ and $\lambda \approx {\rm{\Lambda }}$, and therefore $\psi \approx 0$. This is the reason why we use the rewritten expressions of $A(\phi ,\lambda )$ and $B(\phi ,\lambda )$, Equations (9) and (10), which are robust against the smallness of $\xi =\phi -{\rm{\Phi }}$ and $\eta =\lambda -{\rm{\Lambda }}$. Refer to Appendix A for a more detailed explanation.

The analytical expression of ${I}_{n}(r,\phi ,\lambda )$, Equation (14), is not found in the literature. Although the lowest degree expression, ${I}_{0}(r,\phi ,\lambda )$, is mathematically equivalent with the existing form (Heck & Seitz 2007; Equation (20)), the rewritten expression here is numerically more stable in the sense that no cancellation occurs when ξ and η are small, then $A(\phi ,\lambda )$ and $B(\phi ,\lambda )$ are also small, thus $S(r,\phi ,\lambda )\approx | \zeta |$, and therefore $L(r,\phi ,\lambda )\approx \mathrm{ln}(\zeta +| \zeta | )$ becomes large if $\zeta \lt 0$.

Note that if the layer is thin, the above expression of the kernel functions in a differenced form, Equation (13), suffers from another kind of information loss. Refer to Appendix B for its cancellation-free rewriting. The robustness of the rewritten kernel functions is the key factor to ensure an accurate integration of the gravitational potential, as we show below.

Now, the problem is reduced to the evaluation of a weakly singular two-dimensional integral. Since the integral cannot be analytically evaluated in a closed form, we evaluate it by a method of numerical quadrature.

When the evaluation point is inside or on the surface of the layer, the kernel functions, ${K}_{n}(\phi ,\lambda )$, contain a logarithmically blowing-up singularity at the point of angular coincidence, namely when $\xi =\eta =0$. This is because, in this case, the argument of the logarithm in $L(r,\phi ,\lambda )$ reduces to the form of $\zeta +| \zeta |$, which can completely vanish if $\zeta \leqslant 0$. In addition, even if the evaluation point is outside the layer, such that $\zeta \gt 0$, the main functions, $L(r,\phi ,\lambda )$ and $S(r,\phi ,\lambda )$, have a sharp peak or a steep kink at the point of angular coincidence when $| \zeta |$ is sufficiently small. In such cases, the behavior of the integrand is difficult approximate as a low-degree polynomial, which is the commonly assumed condition on which most of the ordinary quadrature rules are based. As a result, practical difficulties appear in the numerical integration by the standard numerical integrators like the Gauss–Legendre quadrature or the Romberg rule (Press et al. 2007; Chapter 4).

In order to circumvent this difficulty, we adopt the DE rule (Takahashi & Mori 1973, 1974) as a method of numerical integration. In short, the DE rule is the Romberg rule applied to the integrand transformed by the so-called TANH–SINH transformation and other integral variable transformations of DE nature (Press et al. 2007; Section 4.5). It is known as one of the best quadrature rules (Bailey et al. 2005). This is because the integration error typically decreases exponentially with respect to the number of integrand evaluations.

Moreover, the DE rule is easy to use since its algorithm is adaptive. In fact, one only has to specify δ, the input relative error tolerance. Then, the algorithm automatically takes care of test integrations in the Romberg rule so as to minimize the number of integrand evaluations while satisfying the error requirement. Experimentally, we learn that the algorithm is so reliable in the sense that the achieved relative errors are usually smaller than the input requirement if the length of the integration interval after the transformation is appropriately chosen. For the readers' convenience, we added Appendix C, where we describe ${\mathtt{dqde}}$, a Fortran subroutine of the DE rule in the double-precision environment with an optimal setting of the integration interval.

The DE rule is known to be able to integrate weakly singular and/or nearly singular integrals properly if the singularities and/or the point of peaks are located at the endpoints of the integration interval (Mori 1985). In the present case, the singularities and/or the sharp peaks are located at the point of angular coincidence, namely when $\xi =\eta =0$. Therefore, we first change the integration variables from ϕ and λ to ξ and η, respectively. Next, we split the integration interval of ξ into positive and negative parts. In general, one must also split the integration interval of η similarly. However, we only have to shift the transformed integration interval as $[0,2\pi )$ in terms of not λ, but η. This is because the integrand is of $2\pi$ periodicity with respect to λ.

The final form of the split quadrature becomes

$$\begin{eqnarray}&&{V}_{n}(R,{\rm{\Phi }},{\rm{\Lambda }})={\int }_{-\pi /2-{\rm{\Phi }}}^{0}{P}_{n}(\xi )d\xi +{\int }_{0}^{\pi /2-{\rm{\Phi }}}{P}_{n}(\xi )d\xi ,\end{eqnarray} \tag{ 28 }$$

where the new integrand ${P}_{n}(\xi )$ is a definite line integral defined as

$$\begin{eqnarray}&&{P}_{n}(\xi )\equiv {\int }_{0}^{2\pi }{Q}_{n}(\xi ,\eta )d\eta ,\end{eqnarray} \tag{ 29 }$$

while ${Q}_{n}(\xi ,\eta )$ is an abbreviation expressed as

$$\begin{eqnarray}&&{Q}_{n}(\xi ,\eta )\equiv {{GR}}_{0}^{2}{\rho }_{n}({\rm{\Phi }}+\xi ,{\rm{\Lambda }}+\eta ){K}_{n}({\rm{\Phi }}+\xi ,{\rm{\Lambda }}+\eta ).\end{eqnarray} \tag{ 30 }$$

Of course, if the object is homogeneous, then we can set N = 0 and simplify the integrand Q₀ as ${Q}_{0}(\xi ,\eta )\,={{GR}}_{0}^{2}{\rho }_{00}{K}_{0}({\rm{\Phi }}+\xi ,{\rm{\Lambda }}+\eta )$, where ${\rho }_{00}$ is the constant density value of the object. We evaluate all the line integrals in Equations (28) and (29) by the DE rule.

Let us move on to the computation of the gravitational force. In this case, the situation is significantly different. Even after the radial integration is analytically conducted, the numerical integration of its surface integral form suffers from the algebraic singularities (Klees & Lehmann 1998). Therefore, we adopt another approach: the numerical partial differentiation of the numerically integrated potential (Fukushima 2016a, 2016d, 2017). Namely, assuming the availability of $V({\boldsymbol{X}})$ for arbitrary values of ${\boldsymbol{X}}$, we obtain ${\boldsymbol{a}}({\boldsymbol{X}})$, and other quantities defined as the partial derivatives in general, by numerically evaluating their definitions as partial derivatives of $V({\boldsymbol{X}})$.

As the method of numerical differentiation, we previously used Ridders' method (Ridders 1982). Of course, this is a sure method. Nevertheless, it is also true that its computational cost is significantly higher, say a few tens of times more, than that of the potential. Thus, in order to reduce the computational labor while maintaining as high an accuracy as possible, we replace it with a conditional switch between the central difference formula and the single-sided difference formulas (Fukushima 2017). More specifically speaking, as long as the following test displacement of the evaluation point does not cross the surface of the layer, we approximate the gravitational acceleration components in the spherical polar coordinate system by the central difference formula as

$$\begin{eqnarray}{a}_{R}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \approx & [V(R+{\rm{\Delta }}R,{\rm{\Phi }},{\rm{\Lambda }})\\ & & -V(R-{\rm{\Delta }}R,{\rm{\Phi }},{\rm{\Lambda }})]/(2{\rm{\Delta }}R),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}{a}_{{\rm{\Phi }}}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \approx & [V(R,{\rm{\Phi }}+{\rm{\Delta }}{\rm{\Phi }},{\rm{\Lambda }})\\ & & -V(R,{\rm{\Phi }}-{\rm{\Delta }}{\rm{\Phi }},{\rm{\Lambda }})]/(2R{\rm{\Delta }}{\rm{\Phi }}),\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}{a}_{{\rm{\Lambda }}}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \approx & [V(R,{\rm{\Phi }},{\rm{\Lambda }}+{\rm{\Delta }}{\rm{\Lambda }})\\ & & -V(R,{\rm{\Phi }},{\rm{\Lambda }}-{\rm{\Delta }}{\rm{\Lambda }})]/[2(R\cos {\rm{\Phi }}){\rm{\Delta }}{\rm{\Lambda }}],\end{eqnarray} \tag{ 33 }$$

by setting the test displacements in the spherical polar coordinates as

$$\begin{eqnarray}&&{\rm{\Delta }}R/R={\rm{\Delta }}{\rm{\Phi }}={\rm{\Delta }}{\rm{\Lambda }}={\delta }_{2},\end{eqnarray} \tag{ 34 }$$

where we set the relative displacement as

$$\begin{eqnarray}&&{\delta }_{2}=\sqrt[3]{\delta }.\end{eqnarray} \tag{ 35 }$$

This expression of the test displacement is obtained by equating the approximation error of the central difference formula with the input relative error of the potential integration (Press et al. 2007; Section 5.7) after approximating the behavior of the potential and its third-order partial derivatives as those of a masspoint.

On the other hand, if any of the above test displacements crosses the surface of the layer, we compute the gravitational acceleration vector by appropriately chosen single-sided difference formulas such as

$$\begin{eqnarray}&&{a}_{R}(R,{\rm{\Phi }},{\rm{\Lambda }})\approx [V(R+{\rm{\Delta }}R,{\rm{\Phi }},{\rm{\Lambda }})-V(R,{\rm{\Phi }},{\rm{\Lambda }})]/({\rm{\Delta }}R),\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{a}_{{\rm{\Phi }}}(R,{\rm{\Phi }},{\rm{\Lambda }})\approx [V(R,{\rm{\Phi }}+{\rm{\Delta }}{\rm{\Phi }},{\rm{\Lambda }})-V(R,{\rm{\Phi }},{\rm{\Lambda }})]/(R{\rm{\Delta }}{\rm{\Phi }}),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}{a}_{{\rm{\Lambda }}}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \approx & [V(R,{\rm{\Phi }},{\rm{\Lambda }}+{\rm{\Delta }}{\rm{\Lambda }})\\ & & -V(R,{\rm{\Phi }},{\rm{\Lambda }})]/[(R\cos {\rm{\Phi }}){\rm{\Delta }}{\rm{\Lambda }}],\end{eqnarray} \tag{ 38 }$$

or

$$\begin{eqnarray}&&{a}_{R}(R,{\rm{\Phi }},{\rm{\Lambda }})\approx [V(R,{\rm{\Phi }},{\rm{\Lambda }})-V(R-{\rm{\Delta }}R,{\rm{\Phi }},{\rm{\Lambda }})]/({\rm{\Delta }}R),\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{a}_{{\rm{\Phi }}}(R,{\rm{\Phi }},{\rm{\Lambda }})\approx [V(R,{\rm{\Phi }},{\rm{\Lambda }})-V(R,{\rm{\Phi }}-{\rm{\Delta }}{\rm{\Phi }},{\rm{\Lambda }})]/(R{\rm{\Delta }}{\rm{\Phi }}),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}{a}_{{\rm{\Lambda }}}(R,{\rm{\Phi }},{\rm{\Lambda }}) & \approx & [V(R,{\rm{\Phi }},{\rm{\Lambda }})\\ & & -V(R,{\rm{\Phi }},{\rm{\Lambda }}-{\rm{\Delta }}{\rm{\Lambda }})]/[(R\cos {\rm{\Phi }}){\rm{\Delta }}{\rm{\Lambda }}].\end{eqnarray} \tag{ 41 }$$

This time, we set the test displacements in the spherical polar coordinates as

$$\begin{eqnarray}&&{\rm{\Delta }}R/R={\rm{\Delta }}{\rm{\Phi }}={\rm{\Delta }}{\rm{\Lambda }}={\delta }_{1},\end{eqnarray} \tag{ 42 }$$

where a new relative displacement is chosen as

$$\begin{eqnarray}&&{\delta }_{1}=\sqrt{\delta }.\end{eqnarray} \tag{ 43 }$$

It is determined by equating the approximation error of the single-sided difference formulas with the input relative error of the potential integration (Press et al. 2007; Section 5.7) after approximating the behavior of the potential and its second-order partial derivatives as those of a masspoint. Of course, the input relative error tolerance of the potential integration, δ, is different from the achieved accuracy of the potential integrated. However, as we show below, we confirm that δ is a reliable indicator of the actual integration error. Therefore, we used the former as an appropriate estimate of the latter.

The expected relative accuracy of the thus approximated gravitational acceleration vector is ${\delta }_{2}^{2}=\sqrt[3]{{\delta }^{2}}$ for the central difference formulas and ${\delta }_{1}=\sqrt{\delta }$ for the single-sided difference formulas, respectively. For example, if we set $\delta ={10}^{-12}$, then ${\delta }_{2}={10}^{-4}$ and ${\delta }_{1}={10}^{-6}$, and we may therefore anticipate the eight-digit accuracy of the gravitational acceleration components if the evaluation point is far from the surface and the six-digit accuracy if near the surface.

The computational cost for each component of the acceleration vector is twice higher than that of the gravitational potential when the central difference formula is employed. Thus, the total CPU time to obtain the full acceleration vector becomes six times longer. On the other hand, when the single-sided difference formula is used, the extra computational cost per component is the same as that of the potential. Since the base potential value can be shared, the total computational labor is only four times greater.

Below, we show the computational performance of the new method. Let us begin with the case where the analytical solution is known: a homogeneous spherical layer. We denote the reference radius and the top and bottom surface radii of the layer by R₀, ${R}_{{\rm{T}}}$, and ${R}_{{\rm{B}}}$, respectively. Then, the gravitational potential of the layer at a given radius R is expressed as

$$\begin{eqnarray}&&V(R)\\ \equiv \,\left\{\begin{array}{ll}2\pi G\rho ({R}_{{\rm{T}}}+{R}_{{\rm{B}}})H, & (0\leqslant R\lt {R}_{{\rm{B}}}),\\ (2\pi G\rho /3)(3{R}_{{\rm{T}}}^{2}-{R}^{2}-2{R}_{{\rm{B}}}^{3}/R), & ({R}_{{\rm{B}}}\leqslant R\lt {R}_{{\rm{T}}}),\\ {GM}/R, & ({R}_{{\rm{T}}}\leqslant R),\end{array}\right.\end{eqnarray} \tag{ 44 }$$

where ρ is the constant density of the layer, $H\equiv {R}_{{\rm{T}}}-{R}_{{\rm{B}}}$ is the thickness of the layer, and M is the mass of the layer, expressed as

$$\begin{eqnarray}&&M\equiv (4\pi /3)\rho ({R}_{{\rm{T}}}^{2}+{R}_{{\rm{T}}}{R}_{{\rm{B}}}+{R}_{{\rm{B}}}^{2})H.\end{eqnarray} \tag{ 45 }$$

The radial component of the gravitational force is expressed as

$$\begin{eqnarray}&&a(R)\equiv \displaystyle \frac{{dV}(R)}{{dR}}\\ \,=\,\left\{\begin{array}{ll}0, & (0\leqslant R\lt {R}_{{\rm{B}}}),\\ (4\pi G\rho /3)(-R+{R}_{{\rm{B}}}^{3}/{R}^{2}), & ({R}_{{\rm{B}}}\leqslant R\lt {R}_{{\rm{T}}}),\\ -{GM}/{R}^{2}, & ({R}_{{\rm{T}}}\leqslant R).\end{array}\right.\end{eqnarray} \tag{ 46 }$$

As a test case, we set ${R}_{0}=6380\,\mathrm{km}$, ${R}_{{\rm{T}}}={R}_{0}+10\,\mathrm{km}$, and ${R}_{{\rm{B}}}={R}_{0}-40\,\mathrm{km}$ so as to model a spherical crust of the Earth. Then, we measured the relative errors of the gravitational potential integrated by the new method by comparing with the exact solution as

$$\begin{eqnarray}&&\delta V\equiv ({V}_{\mathrm{integrated}}-{V}_{\mathrm{exact}})/{V}_{\mathrm{exact}}.\end{eqnarray} \tag{ 47 }$$

Figure 3 illustrates the measured errors as functions of R for various values of δ, the input relative error tolerance, in a logarithmic manner. The obtained errors are almost constant with respect to R, namely they are independent of the location of the evaluation point. The relative errors of the integrated potential is roughly the same as δ if $\delta \gt {10}^{-10}$. Otherwise, the errors increase by about a factor 30.

Zoom In Zoom Out Reset image size

Let us move to the error of the computed acceleration vector. This time, the relative errors are not well defined since the exact value vanishes when $R\lt {R}_{{\rm{B}}}$. Instead, we introduce the normalized error of acceleration defined as

$$\begin{eqnarray}&&\delta a\equiv ({a}_{\mathrm{integrated}}-{a}_{\mathrm{exact}})/{a}_{\max },\end{eqnarray} \tag{ 48 }$$

where ${a}_{\max }\equiv a({R}_{{\rm{T}}})$ is the maximum value of a(R). Figure 4 depicts the normalized errors of the gravitational acceleration computed by the new method as functions of R. This time, they are mostly unchanged if outside the layer. In fact, the measured errors are roughly equal to $\sqrt[3]{{\delta }^{2}}$ as we expected. Inside the layer, however, the errors increase roughly by a factor of 10–30. Also, the errors become large near the surface, say $6\times {10}^{-5}$ when $\delta ={10}^{-12}$.

Zoom In Zoom Out Reset image size

If we add the contribution of the main part of the Earth, namely that of its mantle and core, the relative magnitude of the computational errors of the crustal component become much smaller, say, by a factor of 10⁻³–10⁻⁴. Therefore, we think that the new method can be sufficiently accurate in the computation of the gravitational field. This is the case if the considered body is nearly spherical in the sense that the contribution of the nonspherical part is relatively small such as the cases of Ceres, Vesta, Deimos, Phobos, and especially the asteroid Ryugu, the target object of the HAYABUSA2mission (Müller et al. 2016).

Finally, we measure the computational labor of the new method. Figure 5 depicts ${N}_{\mathrm{func}}$, the number of transcendental function evaluations required in computing the gravitational potential of the spherical layer by the new method. The figure illustrates the R-dependence of ${N}_{\mathrm{func}}$. We omit the similar graphs showing the Φ- and Λ-dependences, respectively, since the results are practically the same as those depicted in Figure 5. At this scale, the computational cost seems to be independent on the location of the evaluation point. This situation is the same as in the polyhedron model. However, as will be seen later, the computational cost in actual cases reduces when the distance from the object increases significantly.

Zoom In Zoom Out Reset image size

As a more realistic example, let us investigate the gravitational field of 433 Eros. This is one of the most difficult cases of the gravitational field computation because the shape of the asteroid is significantly irregular. As its surface function, ${R}_{{\rm{T}}}(\phi ,\lambda )$, we used the 24 × 24 spherical harmonic expansion (Miller et al. 2002). Its bird's eye view as well as the associated contour map are shown in Figure 6. By setting ${R}_{{\rm{B}}}(\phi ,\lambda )=0$ for simplicity and assuming the homogeneity of the asteroid, we integrated the gravitational potential, $V(R,{\rm{\Phi }},{\rm{\Lambda }})$, by setting $\delta ={10}^{-12}$, and computed the acceleration vector, ${\boldsymbol{a}}(R,{\rm{\Phi }},{\rm{\Lambda }})$, by setting ${\delta }_{2}={10}^{-4}$ far from the surface and ${\delta }_{1}={10}^{-6}$ near the surface.

Zoom In Zoom Out Reset image size

Figures 7 and 8 depict a bird's eye view and the associated contour map of the integrated potential and the magnitude of the computed acceleration vector at a constant radius, respectively. They are drawn as two-dimensional maps with respect to the angular coordinates. The chosen radius, $R=2{R}_{0}$, represents a case of the intermediate zone, namely inside the Brillouin sphere but above the bottom of the surface of the object.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The displayed view/map of the potential are much smoother than those of the surface function (Figure 6). This is due to the averaging nature of the gravitational field computation as we had experienced with that of the spiral galaxy M74 (Fukushima 2016b). On the other hand, those of the acceleration magnitude have a sharp kink on two oval-shaped loci. They roughly correspond to the cross-section of the surface of Eros cut by the sphere of constant radius, $R=2\,{R}_{0}$. This is natural because the density is discontinuous across the loci: finite inside the loci and zero outside them.

In order to show the nonanalyticity of the gravitational field more clearly, we prepared Figures 9 and 10, displaying the contour map of the gravitational potential and the magnitude of the gravitational acceleration evaluated on the equator, namely when ${\rm{\Phi }}=0$, respectively. The contours of the potential value become dense near the surface of the object. This implicitly shows the nonanalyticity there. The situation is similar to that of the potential map of some finite axisymmetric bodies (Fukushima 2016d, Appendices A and C). Also, the contours of the acceleration magnitude exhibit a sharp boundary, which exactly coincides with the surface of the object. This is the same feature we experienced with that of the spherical layer in the previous subsection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With an aim to clarify the nonanalytic feature of the gravitational field, we examined the radial variation of the gravitational field more closely. Refer to Figure 11, plotting the potential and the three acceleration components along the direction passing through the lowest point of the surface of Eros. Again, we confirm that the acceleration components are continuous but not differentiable at the surface. This is because the potential curves are not analytic there.

Zoom In Zoom Out Reset image size

We noticed that, even outside the object, all the acceleration components are not monotonically decreasing when R increases. The phenomenon is prominent near the surface, say, when $0.5\lt R/{R}_{0}\lt 0.7$. Indeed, their peak values are achieved there. This interesting feature is caused by the gravitational attraction of the mountainous area surrounding the lowest point.

On the other hand, Figure 12 illustrates the radial variation of the computational cost of the gravitational potential. When the distance from the asteroid is large, the potential value is almost inversely proportional to R, the radial coordinate. Therefore, for an efficient computation, one may increase the relative error tolerance so as to be in proportion to R. This policy decreases the computational cost as seen in Figure 12. This feature is different from the polyhedron model, which requries the same computational cost independently of the location of the evaluation point. Having said so, we must say that the spherical or spheroidal harmonic expansions are much more efficient if outside of the Brillouin sphere or spheroid.

Zoom In Zoom Out Reset image size

Let us focus on the gravitational feature on the surface, which is most difficult to obtain by the existing methods including the polyhedron models. Figure 13 illustrates the bird's eye view and the associated contour map of the surface gravitational potential of the asteroid. By comparing it with Figure 6, we notice that the magnitude of the surface potential is inversely correlated with the shape function. Namely, it becomes small in the mountainous area and large in the basins and valleys.

Zoom In Zoom Out Reset image size

Next, we prepared Figure 14, showing the similar view and map of the magnitude of the surface gravitational acceleration vector of the asteroid. This time, the variation of the magnitude of the acceleration vector is significantly smaller than those of the potential. Also, the magnitude of the acceleration is smaller not only in the mountainous area but also in the basins and valleys. Third, we calculated the deflection angle of the surface acceleration vector,

$$\begin{eqnarray}&&\theta \equiv {\tan }^{-1}(\sqrt{{a}_{{\rm{\Phi }}}^{2}+{a}_{{\rm{\Lambda }}}^{2}}/{a}_{R}).\end{eqnarray} \tag{ 49 }$$

Its bird's eye view and the associated contour map are provided in Figure 15. The deflection angle is large in most areas, especially in the steep slopes of the mountainous area. These facts will be important in designing the descending orbit of a space probe to land on irregular shaped objects (Liu et al. 2017b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We begin with the expression of the normalized mutual distance, S. Using the spherical polar coordinates of the internal and evaluation points, ${\boldsymbol{x}}$ expressed as $(r,\phi ,\lambda )$ and ${\boldsymbol{X}}$ expresed as $(R,{\rm{\Phi }},{\rm{\Lambda }})$, we write the mutual distance as

$$\begin{eqnarray}&&| {\boldsymbol{x}}-{\boldsymbol{X}}| =\sqrt{{r}^{2}-2{rR}\cos \psi +{R}^{2}},\end{eqnarray} \tag{ 50 }$$

where ψ is the angular distance between the direction vectors, ${\boldsymbol{x}}/r$ and ${\boldsymbol{X}}/R$, such that its cosine is expressed as

$$\begin{eqnarray}&&\cos \psi =({\boldsymbol{x}}\cdot {\boldsymbol{X}})/({rR})=\sin \phi \sin {\rm{\Phi }}\\ &&+\,\cos \phi \cos {\rm{\Phi }}\cos (\lambda -{\rm{\Lambda }}).\end{eqnarray} \tag{ 51 }$$

This expression is fragile against the proximity of ${\boldsymbol{x}}$ and ${\boldsymbol{X}}$, namely when $\phi \approx {\rm{\Phi }}$, $\lambda \approx {\rm{\Lambda }}$, and/or $r\approx R$. In order to overcome the numerical problems caused by the proximity, we introduce new variables and an auxiliary quantity defined as

$$\begin{eqnarray}&&\xi \equiv \phi -{\rm{\Phi }},\,\,\eta \equiv \lambda -{\rm{\Lambda }},\,\,\zeta \equiv (r-R)/{R}_{0},\,\,\alpha \equiv R/{R}_{0},\end{eqnarray} \tag{ 52 }$$

where R₀ is a certain scale constant. Using them, we first rewrite $\cos \psi$ as

$$\begin{eqnarray}\cos \psi & = & \sin \phi \sin {\boldsymbol{\Phi }}+\cos \phi \cos {\boldsymbol{\Phi }}-\cos \phi \cos {\boldsymbol{\Phi }}\\ & & \times [1-\cos (\lambda -{\rm{\Lambda }})]\\ & = & \cos \xi -\cos \phi \cos {\boldsymbol{\Phi }}(1-\cos \eta ).\end{eqnarray} \tag{ 53 }$$

Next, we obtain a more robust form representing ψ as

$$\begin{eqnarray}{\sin }^{2}(\psi /2) & = & (1-\cos \psi )/2=(1-\cos \xi )/2\,+\,\cos \phi \cos {\boldsymbol{\Phi }}\\ & & \times (1-\cos \eta )/2\\ & = & {\sin }^{2}(\xi /2)+\cos \phi \cos {\rm{\Phi }}{\sin }^{2}(\eta /2).\end{eqnarray} \tag{ 54 }$$

After this preparation, we rewrite the normalized mutual distance, S, as

$$\begin{eqnarray}S & \equiv & | {\boldsymbol{x}}-{\boldsymbol{X}}| /{R}_{0}=\sqrt{{(r-R\cos \psi )}^{2}+{R}^{2}{\sin }^{2}\psi }/{R}_{0}\\ & = & \sqrt{{(\zeta +B)}^{2}+A},\end{eqnarray} \tag{ 55 }$$

where B and A are expressed as

$$\begin{eqnarray}B & = & (R/{R}_{0})(1-\cos \psi )=2\alpha {\sin }^{2}(\psi /2)\\ & = & 2\alpha [{\sin }^{2}(\xi /2)+\cos \phi \cos {\boldsymbol{\Phi }}{\sin }^{2}(\eta /2)],\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}A & = & {(R/{R}_{0})}^{2}{\sin }^{2}\psi ={\alpha }^{2}{[2\sin (\psi /2)\cos (\psi /2)]}^{2}\\ & = & 2\alpha {\sin }^{2}(\psi /2)[2\alpha -2\alpha {\sin }^{2}(\psi /2)]=B(2\alpha -B).\end{eqnarray} \tag{ 57 }$$

Thus, we arrived at the expression of S in terms of A and B in the main text.

Next, we consider the radial integral I_n. It is compactly expressed as

$$\begin{eqnarray}{I}_{n} & = & \displaystyle \frac{1}{{R}_{0}^{n+3}}{\displaystyle \int }^{r}\displaystyle \frac{{(r^{\prime} )}^{n+2}{dr}^{\prime} }{S}={\displaystyle \int }^{\zeta }\displaystyle \frac{{(\zeta ^{\prime} +\alpha )}^{n+2}d\zeta ^{\prime} }{\sqrt{{(\zeta ^{\prime} +B)}^{2}+A}}\\ & = & {\displaystyle \int }^{\zeta +B}\displaystyle \frac{{(t+\alpha -B)}^{n+2}{dt}}{\sqrt{{t}^{2}+A}}\\ & = & \displaystyle \sum _{j=0}^{n+2}\left(\begin{array}{c}n+2\\ j\end{array}\right){(\alpha -B)}^{n+2-j}{L}_{j}(\zeta +B),\end{eqnarray} \tag{ 58 }$$

where L_j(t) is an indefinite integral defined as

$$\begin{eqnarray}&&{L}_{j}(t)\equiv {\int }^{t}\displaystyle \frac{{\tau }^{j}d\tau }{\sqrt{{\tau }^{2}+A}}.\end{eqnarray} \tag{ 59 }$$

The integrand of L_j(t) is of the form of a monomial of τ divided by a square root of a simple quadratic polynomial of τ. Thus, it can be analytically integrated in terms of elementary functions. For example, the first two of them are explicitly obtained as

$$\begin{eqnarray}&&{L}_{0}(t)=\mathrm{ln}(t+\sqrt{{t}^{2}+A})=L,\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{L}_{1}(t)=\sqrt{{t}^{2}+A}=S.\end{eqnarray} \tag{ 61 }$$

Namely, they reduce to L and S introduced in the main text. For the computation of the higher order integrals, L_j(t), we employ its recurrence relation as

$$\begin{eqnarray}&&{{jL}}_{j}(t)={t}^{j-1}S-(j-1){{AL}}_{j-2},\,\,(j=2,\mathrm{3,...}),\end{eqnarray} \tag{ 62 }$$

by starting from a pair of initial values, ${L}_{0}(t)$ and ${L}_{1}(t)$. The correctness of the recurrence relation is easily confirmed by the direct differentiation. At any rate, L_j(t) are mechanically obtained from the recurrence formula such as

$$\begin{eqnarray}&&2{L}_{2}={tS}-{AL},\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&3{L}_{3}=({t}^{2}-2A)S,\end{eqnarray} \tag{ 64 }$$

where we omitted the dependency on t for further simplicity. Refer to Karcol (2011) presenting a slightly different derivation. In any case, it is trivial to see that these integration formulas result in the expression of ${C}_{j}(\phi ,\lambda )$ and ${D}_{j}(r,\phi ,\lambda )$ for $j=0,1,2,3$ provided in the main text.

Here we provide a couple of Fortran programs to evaluate a line integral expressed as

$$\begin{eqnarray}&&s\equiv {\int }_{a}^{b}f(x){dx},\end{eqnarray} \tag{ 76 }$$

where f(x) is either (i) regular in the whole integration interval, $a\leqslant x\leqslant b$, or (ii) regular inside the interval, $a\lt x\lt b$, and has integrable singularities or sharp peaks at one or multiple endpoints, namely when x = a and/or x = b. As a basic method of its numerical evaluation, we adopt the DE rule using the TANH–SINH rule (Takahashi & Mori 1974).

The most problematic issue in its practical implementation is the truncation of the integration interval after the variable transformation, which is infinite by definition. More specifically speaking, we restrict the infinite interval, $(-\infty ,+\infty )$, into a finite one, $(-{h}_{\max },+{h}_{\max })$. Thus, the question is the choice of ${h}_{\max }$. If it is too small, the computational precision is degraded. Meanwhile, if it is too large, the computational labor becomes too much. Therefore, there is an optimal value of ${h}_{\max }$, which may depend on the general nature of the integrand before the transformation and the input relative error tolerance, δ.

After some preliminary experiments, we fixed its functional form of ${h}_{\max }$ to be a piecewise polynomial of ${\mathrm{log}}_{10}\delta$ as

$$\begin{eqnarray}&&{h}_{\max }={p}_{j}+{q}_{j}(-{\mathrm{log}}_{10}\delta ),\\ &&\qquad ({t}_{j-1}\lt -{\mathrm{log}}_{10}\delta \lt {t}_{j};j=1,2,\,...,\,J),\end{eqnarray} \tag{ 77 }$$

where p_j, q_j, and t_j and the total number of pieces J are free parameters to be tuned. Then, by restricting the type of integrand to $\mathrm{log}(x-a)$ or $\mathrm{log}(b-x)$, we experimentally determined an optimal set of the parameters as

$$\begin{eqnarray}&&{p}_{1}=2.75,\,\,{q}_{1}=0.75,\,\,{t}_{1}=4,\end{eqnarray} \tag{ 78 }$$

$$\begin{eqnarray}&&{p}_{2}=3.75,\,\,{q}_{2}=0.50,\,\,{t}_{2}=6,\end{eqnarray} \tag{ 79 }$$

$$\begin{eqnarray}&&{p}_{3}=5.25,\,\,{q}_{3}=0.25,\,\,{t}_{3}=10,\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{p}_{4}=6.50,\,\,{q}_{4}=0.125,\,\,{t}_{4}=\infty ,\end{eqnarray} \tag{ 81 }$$

while ${t}_{0}=-\infty$.

Consequently, we implemented the algorithm into ${\mathtt{dqde}}$, a Fortran subroutine in the double-precision environment as listed in Table 1. Its quadruple precision extension, ${\mathtt{qqde}}$, is obtained from ${\mathtt{dqde}}$ by converting all the double-precision constants into the corresponding quadruple precision ones such as ${\mathtt{TWOPI}}\,={\mathtt{6}}.{\mathtt{283185307179586476925286766559}}{\mathtt{q}}{\mathtt{0}}$ and modifying the parameters as ${\mathtt{MMAX}}={\mathtt{4092}}$.

In order to show the performance of the implemented Fortran programs, we integrate a test integral defined as

$$\begin{eqnarray}&&F(\zeta )\equiv {\int }_{0}^{2\pi }f(\eta ,\zeta )d\eta ,\end{eqnarray} \tag{ 82 }$$

where the integrand is expressed as

$$\begin{eqnarray}f(\eta ,\zeta )\equiv \left\{\begin{array}{ll}\mathrm{ln}(\zeta +\sqrt{{\eta }^{2}+{\zeta }^{2}}), & (\zeta \geqslant 0),\\ \mathrm{ln}[{\eta }^{2}/(-\zeta +\sqrt{{\eta }^{2}+{\zeta }^{2}})], & (\zeta \lt 0).\end{array}\right.\end{eqnarray} \tag{ 83 }$$

This integrand is obtained by simplifying $L(r,\phi ,\lambda )$ in the main text such that the integral is analytically computable. Indeed, $F(\zeta )$ is explicitly written in a closed form as

$$\begin{eqnarray}&&F(\zeta )=2\pi [f(2\pi ,\zeta )-1]+\zeta \mathrm{ln}[(2\pi +\sqrt{4{\pi }^{2}+{\zeta }^{2}})/| \zeta | ].\end{eqnarray} \tag{ 84 }$$

By using the analytical expression, Equation (84), we examined the achieved accuracy and the computational cost of the Fortran programs. First, Figure 16 illustrates the relation between the achieved integration errors, $| \delta F|$, and the input relative error tolerance, δ. On the other hand, Figure 17 shows the relation between $\delta F$ and the computational cost ${N}_{\mathrm{eval}}$, namely the number of integrand evaluations. These figures indicate a credibility and a good cost performance of the DE rule. Indeed, they are (i) reliable in the sense that the achieved relative error $| \delta F|$ is almost the same as δ, and (ii) efficient because $| \delta F|$ decreases exponentially when ${N}_{\mathrm{eval}}$, the computational cost, increases linearly as ${\mathrm{log}}_{10}\,| \delta F| =-3.5-0.05{N}_{\mathrm{eval}}$. From these two figures, we learn that the programs ${\mathtt{dqde}}$ and ${\mathtt{qqde}}$ are trustable and highly effective.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size