SURF Report: High Accuracy Methods for Computing Gravitational Potential and Gravitational Force Fields Near the Surface of Irregularly Shaped 3-Dimensional Bodies¹¹1Copyright: © 2024. California Institute of Technology

Accurate gravity field calculations are necessary for landing on planets, moons, asteroids, Earth’s “minimoons,” or other irregularly shaped bodies, but current methods become increasingly inaccurate and slow near the surface. We present high accuracy, fast methods for computing gravitational potential and gravitational force fields, which are needed for future space missions. Notably, gravitational force and potential computations are simplified, with high accuracy enhanced by bringing the derivative inside the gravitational potential integral. In addition, we present a new gravitational field “calculus,” which lets us combine simpler potentials and force fields to create more complex ones without accuracy loss. Several examples are provided, for instance, where we subtract different shapes from a spherical body making a variety of craters. The calculus will also work well with volumetric octree methods. Additionally, we use new bounds in the gravitational potential integral, to avoid trying to fit smooth basis functions to non-smooth curves, and harness new computational tools where tasks can be migrated to GPUs. We also have found that cylindrical coordinates can have special advantages in tailoring shape models. We have created a series of algorithms and preliminary MATLAB and Mathematica toolboxes utilizing these methods and the gravitational calculus. These methods are newly customizable for necessary high-accuracy gravity computations in future missions planned by JPL and other space agencies to navigate near irregularly shaped bodies in the solar system.

Computing the gravitational potential field $V$ and gravitational force field $F$ are highly applicable and significant results for navigational mathematics, physics and engineering necessary for landing and modeling near-body dynamics. Throughout this paper we refer to ($V,F$) as the “gravity field.” While the reader may assume that fast, accurate computations of the gravity field near non-spherical objects are already available, this is not the case. Surprisingly, current faster state-of-the-art gravitational force computations exhibit orders of tens of percent error near asteroid and moon surfaces (Scheeres, 2016), (Scheeres et al., 2019), (Takahashi & Scheeres, 2014) and are insufficient for accurately representing the dynamics. Most of these surface gravity models create errors on the order of tens of percent for bodies that stray from spherical symmetry, and the more aspherical the shape is, the more inaccurate the gravitational force models near the surface are (Takahashi & Scheeres, 2014). We foresee that our research opens the door to developing strategies for faster and far more accurate navigation near the surface of bodies.

In particular, space mission design and flight design are increasingly in need of fast accurate surface gravity force computations for the sake of landing, as well as for general navigation near the surface of real bodies. As landing on or flying near asteroids, moons, and other non-spherical bodies in space becomes increasingly pertinent, our current improvement of accuracy in gravity computations near the surface of realistic, non-spherical bodies will be increasingly beneficial.

New Methods: We present new methods for fast, reliable navigation near realistic irregularly shaped bodies that can be used for determining the optimal navigation strategy in a given mission. We present a number of methods that can be used for highly accurate, fast and reliable calculations, with special focus on our “calculus” of gravitational potentials and forces. Additionally, our methodology avoids issues in the “conventional methods” by taking the derivative of the potential integrand, solving the definition of the gravitational force directly, and not working with smooth basis functions to fit to non-smooth curves. See Figure 1 for a non-smooth force curve that is not suitably fitted by smooth basis functions. This work has been accepted as California Institute of Technology Jet Propulsion Laboratory New Technology Report 53321 (Lo et al., 2024).

We address a current issue that lies in computing sufficiently accurate derivatives of standard potential functions to obtain accurate gravitational force fields near the surface (Werner, 2010), (Scheeres, 2016). In the example case of the Mars moon Phobos, computational methods to derive the gravity field with accuracy and speed are crucial to allow for reliable navigation, and the gravitational force calculations on Phobos are especially relevant where the atmosphere is thin and asymmetries make the body highly non-spherical (Willner et al., 2010), (Scheeres et al., 2019).

At the core of our new methods, we focus on matching solutions closely to the “true” gravitational potential field $V$ and gravitational force field $F$. Letting $\nabla$ denote the gradient, the “true gravity field definitions,” with no additional assumptions, are given by Equation 1 and Equation 2 for gravitational potential $V$ and gravitational force $F$, respectively, evaluated at any point $\overrightarrow{\mathbf{x}}\in\mathbf{R}^{3}$.

In Appendix B, Section 8, we present insights into the challenges that arise in the conventional methods, and we also present an implementable algorithm (B.1) derived from Takahashi & Scheeres (2014) which gives the previous state-of-the-art fast computation method.

Looking directly at the definition of the gravitational potential $V$ in real space induced by a shape $\mathbf{S}$ with density function $\rho(\overrightarrow{\mathbf{r}})$, given by Equation 1, we can compute the gravitational potential field $V$ induced by the shape $\mathbf{S}$. Integrating the integrand $\frac{G\rho(\overrightarrow{\mathbf{r}})d\mathbf{v}(\overrightarrow{\mathbf{r}})}{|\overrightarrow{\mathbf{x}}-\overrightarrow{\mathbf{r}}|}$ over $\overrightarrow{\mathbf{r}}\in\mathbf{R}^{3}$, where $G$ is the gravitational constant, $\rho(\overrightarrow{\mathbf{r}})$ is the density function which gives the density of the body at a given point in space, $d\mathbf{v}(\overrightarrow{\mathbf{r}})$ is the volume element, and $|\overrightarrow{\mathbf{x}}-\overrightarrow{\mathbf{r}}|$ is the distance between input vector $\overrightarrow{\mathbf{x}}\in\mathbf{R}^{3}$ and $\overrightarrow{\mathbf{r}}\in\mathbf{R}^{3}$, we derive the expression for the true gravitational potential field $V$ induced by the shape $\mathbf{S}$.

The integration process just described, where we integrate as in Equation 1 over the shape $\mathbf{S}$, is the most direct way of computing the gravitational potential field, and is derived directly from the gravitational potential definition. Staying “near” the direct definition presents a powerful way of limiting error in our solution; we can derive closed-form expressions for the gravity field for certain shapes $\mathbf{S}$ by conventional integration tricks.

The gravitational potential $V_{S}$ for a uniform density sphere is a continuous $C^{1}$ function (see Figure 1), but its derivative exhibits a fundamental issue near the surface. This restricts the standard methods from accurately deriving the gravitational force $F_{S}$ for the sphere, which is the derivative of the potential. The graph of gravitational force may seem to be smooth despite its underlying behavior as a $C^{1}$, not $C^{2}$, function, as demonstrated by Figure 1.

Closed-form expressions for the uniform density sphere’s gravity field can be derived by substituting the sphere into Equation 1, where over $\mathbf{R}^{3}$ the density $\rho({\overrightarrow{\mathbf{r}}})$ is 0 outside of the sphere and 1 inside the sphere.

See Equations 3 and 4 for closed-form gravity field expressions about the uniform mass sphere with $F_{S}$ denoting the gravitational force field about the sphere and $V_{S}$ denoting the gravitational potential field, where we take the positive direction to be pointing inwards and let $||\overrightarrow{\mathbf{x}}||$ denote the norm of vector $\overrightarrow{\mathbf{x}}\in\mathbf{R}^{3}$.

See Figure 1 for a plot of these expressions in any given direction from the center of the sphere, where the gravitational force has a cusp at the sphere boundary.

New Methods to Eliminate Conventional Errors: Real shapes of interest typically have other characteristics beyond spherical symmetry and constant density that may make dynamics near the surface additionally complex. We present a method whereby integrating Equation 1 directly, there is no computational noise which may come from unsuitable basis functions, and taking the gradient of Equation 1 prior to integration results in a highly accurate numerical or closed-form gravitational force field. See Algorithm A.1 in Section 7 for a new procedure for gravity field computation via integration which stays near the true gravity field definitions given by Equations 1 and 2.

Symbolic and Numeric Integration: The triple integrals for computing gravitational potential and force, as given by Equations 1 and 2 evaluated for a given shape $\mathbf{S}$, at a point $\overrightarrow{\mathbf{x}}=(x,y,z)$, may have one, two, or three symbolic integrations depending on the choice of coordinate system $(\phi_{1}(x),\phi_{2}(y),\phi_{3}(z))$, where choosing a convenient coordinate system can significantly speed up the computation while retaining accuracy. Each symbolic integration speeds up the computation because it provides a closed-form expression. This choice of coordinate system is motivated by selection of coordinates which exploit symmetries. We select the coordinate system $(\phi_{1}(x),\phi_{2}(y),\phi_{3}(z))$ which allows for a higher number of symbolic integrations in order to maintain speed in the computation and limit numerical integration.

We find that in cases where there is symmetry about the z-axis and the height can be written as a function of the radial distance from the center, the use of conventional cylindrical coordinates $(\rho,\theta,z)$—where $\rho$ is radial distance in the xy-plane from the origin, $\theta$ is the angle relative to the x-axis in the xy-plane, and $z$ is height—is convenient. This choice of coordinates lends way to two symbolic integrations, following the definitions in Equations 1 and 2. Spherical symmetry leads to a constant factor from integration in $\theta$, and we can write the bounds in $z$ as a function of planar coordinate $\rho$, which allows for symbolic integration in $z$. These properties hold for many notable volumes which may represent true surface features, such as craters. An example of the setup for computing the gravity field of a crater volume, where integration is simplified with two symbolic integrations via choice of cylindrical coordinates, is shown in Figure 2.

New “Gravitational Calculus:” We present uses of our new “gravitational calculus,” where computation of disjoint shapes’ gravity field can be translated or rotated, then added or subtracted to the overall gravity field. See Figures 3, 4, and 5, for examples of shapes with craters and mountains. The operations of translation and rotation via the gravitational calculus are given by Algorithm A.3.1 and A.3.2 in Appendix A, Section 7.

Certain closed-form solutions can be used in our new type of “gravitational calculus.” Closed-form solutions are possible when a convenient choice of coordinate system allows for three symbolic integrations. This is the case for the gravity field about the sphere. Additionally, highly accurate numerical computations can be used as a part of the calculus, where once a shape’s gravity field has been computed, use of the calculus allows for translation or rotation of the shape. In cases where one or more numerical integrations are computed to represent the gravity field of a shape $\mathbf{S}$, although not exact, that gravity field can also be utilized further via the gravitational calculus with rotation or translation.

Navigation: For navigation purposes, it is often necessary to investigate crater-filled shapes. Scientific classifications exist for craters, with an example of a simple and complex crater shown in Figure 6, while other craters may have even more complicated shapes, which carry other classifications (Kiefer, 2003). Impact craters may have specific morphological shapes including mounds, flat floors, or concentric shapes (Senft & Stewart, 2007).

A potential shape to investigate, in order to model a given crater, would have some bounds in radius, within a choice coordinate system (for instance, spherical coordinates for the following consideration), that are functions $c_{lower}(\phi)$ and $c_{upper}(\phi)$ of $\phi$ chosen in such a way of closely representing the “real crater.” This can be used to model real craters that form and lie in real bodies such as Phobos. See Figure 5 for a sphere with a “test crater” which algebraically simplifies since $c_{lower}(\phi)$ and $c_{upper}(\phi)$ are constants, where use of the gravitational calculus presents a gravity field computation for an asteroid-like shape.

We present additional convenient applications of the gravitational calculus in Figures 3 and 4, where a shape is rotated for use in computing the final shape’s desired gravity field, and the gravitational calculus makes the computations convenient. The gravitational calculus allows us to model mountains, craters, or other surface features.

Octree Methods: We present another new method, which utilizes octrees for shape and density representations to compute the gravity field in Algorithm A.2 of Appendix A, Section 7 (Meagher, 1982). Algorithm A.2 provides a guide for implementing octree representation. The use of quad tree algorithms has been examined in aerospace applications for shape representation of satellite mapping coverage (Lo et al., 2003), and many of the same benefits coming from this form of shape representation of completeness and efficiency may apply in our case of gravity field computations, given our gravitational calculus. The octree is a three-dimensional generalization of the quad tree representation, and this method allows for complete shape representation in three dimensions.

We have demonstrated that the nature of the gravitational force, with non-smooth behavior, encourages use of new fast methods for high accuracy, yet efficient, computations, which we propose. These methods are summarized by the algorithms that are stated in the following section.

While it has been previously demonstrated that within the minimum Brillouin sphere, where all mass is enclosed, the force function is not smooth by nature (Costin et al., 2020), (Šprlák & Han, 2021), we find that, in fact, anytime a body’s section crosses a mass boundary, such as in craters or mountains, there can be cusps arising in the true force function. See Figure 7 for an example of the gravity field about a crater modeled to fit the form of a “simple crater,” where we have a seemingly innocent potential function whose derivative has two cusps, which correspond to the points on the boundary of the crater volume.

We have developed tools for integrating over certain “crater sections” by directly integrating the potential definition in Mathematica and MATLAB code. See Appendix (Section 7).

Algorithms: We have constructed a set of algorithms which guide computation of gravitational potential field $V$ and gravitational force field $F$ for a given shape $\mathbf{S}$. See Appendix A, Section 7, for full set of new method algorithms. Algorithm A.1 states a workflow for computing the gravitational potential and force fields. The algorithms that follow, in Section 7, present direct methods involving integration in solving for $V$, $F$, the octree method (Algorithm A.2), gravitational calculus operator properties (Algorithms A.3.1, A.3.2), and tools for simplifying computations using the gravitational calculus and properties of the gravitational potential (Algorithms A.4, A.5). See Section 5 for additional future pathways for computing gravity fields. In Appendix B, we condense the previous conventional methods into a convenient algorithm; Algorithm B.1 re-states in brief, implementable fashion, a state-of-the-art method as presented in Takahashi & Scheeres (2014). This algorithm is also implemented in Mathematica notebook “Potential Reconstruction Algorithm,” created for demonstration of the previous methods versus the new methods. Figure 11 demonstrates the success of following the new proposed gravity field workflow in the case of the uniform mass sphere, where there is only numerical error, whereas in the conventional methods, there is additional error (See Appendix B, Section 8).

MATLAB Tools: We have a preliminary MATLAB package, GravToolbox, which computes, and plots in a 2-dimensional cross-section, the numerical potential of a specified shape. Inputting the bounds which represent the crater shape (as a spherical section which can be expressed in terms of numerical bounds for $(r,\theta,\phi)$ in spherical coordinates) into our functions, we can compute a (highly accurate) numerical cross-section of the potential and force field. Figure 7 shows the gravity field computation for a uniform density crater modelled to fit the shape of a “simple crater.”

Mathematica Tools: We have preliminary Mathematica notebooks which compute and plot in a 2-dimensional cross-section, the numerical potential of a specified shape, as done in the MATLAB package GravToolbox. We also have an example notebook “Potential Reconstruction Algorithm” which demonstrates Algorithm 1 and error measurements (see Appendix, Section 7).

As we have begun to build the calculus for the gravitational potential and gravitational forces, we have found that many of the potential integrations become algebraically complicated, leading to long run times and challenging computations. Taking the derivative of the potential’s integrand prior to integration to find the gravitational force field can lead to a more challenging integral than that of the potential integral. However, these run times can be sped up by some optimal combination of pre-computed octrees or pre-computed “definition integrals” in convenient coordinate systems and bounds for shapes, each of which need only be computed once. Our methods, at their finest, are highly accurate, but potentially slow in such a way that certain computations may be infeasibly slow without GPUs. This is the issue with the polyhedron model (Werner, 1997), (Werner, 2010). However, we anticipate that computation of shapes in parallel using GPUs creates a modern, sufficiently fast method which retains high accuracy up to specified tolerance. We can ensure a desired accuracy, up to some tolerance, in the case of computing numerically by computing the exact potential at a sufficiently large number of points and interpolating. Similarly, with use of octrees, we can continue the algorithm up to construction of a shape that is sufficiently close to the true shape of interest. We also anticipate that there may be special cases where a scaling property could be used with the gravitational calculus to speed up computations.

In the example case of Phobos, a highly accurate form, we predict, would require adaptability to a wide span of shapes that constitute the body in union, and this may involve numerical and octree computations in conjunction with our potential “calculus” (see (Werner, 1997) for Polyhedron approximation of non-spherical bodies). The approach would also extend to navigation near the Earth’s “minimoons.” These accurate computations model gravitational force fields which are very “cusp-like” in behavior, and are not limited to smooth functions.

Further work involves a deep investigation into how changing integration order, reformatting integrands, optimizing code, finding shapes that represent “true craters” in ways that are more conveniently derivable, changing coordinate systems, or using more powerful computing can be combined in such a way that creates a directly implementable algorithm for space agencies to use. Additionally, new algorithms must be able to model non-constant density distributions, which may be done using the octree method in conjunction with machine learning algorithms, to fit a power series density function to each closed-form box in the octree. Our GravToolbox package can be expanded into a re-usable set of code for navigational purposes and our algorithms in Section 7 can also be implemented for specific gravity field computations in missions requiring high accuracy navigation near bodies.

Additional Methods: It may be possible to use a different solution to Poisson’s Equation $\nabla^{2}U=-4\pi G\rho$ using complete basis functions such as wavelets (Iglewska-Nowak & Stefaniak, 2023) and a more complete representation of the problem itself as the single differential equation given by Equation 5 where for an ellipsoid with axes a, b, c, in the x, y, and z directions respectively, $f(x,y,z)=1-\frac{x^{2}}{a^{2}}-\frac{y^{2}}{b^{2}}-\frac{z^{2}}{c^{2}}$.

One solution method in the new formulation involving the Heaviside function representation could be to use a Monte Carlo algorithm which samples points randomly. Other investigations could involve a Fourier Transform Method. Others could investigate other distributions to represent the solution. Speed-up and parallel computations may utilize GPUs.

Thomas MacLean is supported as a William H. and Helen Lang Summer Undergraduate Research Fellow with funding support provided by William H. Lang and Helen Lang and the California Institute of Technology (Caltech) Student-Faculty Programs Office. Dr. Martin Lo at the NASA Jet Propulsion Laboratory has been a vital source of ideas and collaborator for the research. This research is directly funded and supported by the Caltech Barr Graphics Lab. The NASA Jet Propulsion Lab has provided funding, allocated via the Mission Design Section and made possible by Dr. Jon Sims and Dr. Martin Lo. The research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

© 2024. California Institute of Technology. Government sponsorship acknowledged.

Packages: GravToolbox. Presented in MATLAB, this directory gives a set of files which can compute specified crater or base sphere gravity fields in cylindrical coordinates. The package also generalizes as a convenient tool for computations of shapes which are a linear combination of craters and a base shape. Potential Reconstruction Algorithm. Mathematica Demonstration of Algorithm 1 versus direct integration. The package features a presentation of a current state-of-the-art implementation on a sphere, where there is some error as shown in Figure 11, while the numerical method used has extremely low error.

Appendix B analyzes the challenges that arise in conventional methods of gravitational potential field (here denoted by $U$ to match previous notations) and gravitational force field ($\nabla U$) calculations. We also present an implementable algorithm (B.1) giving the state-of-the-art fast computation method. Most methods use spherical harmonics expansions, though ellipsoidal expansions have also been used, and in certain cases improved accuracy (Romain & Jean-Pierre, 2001). These, however, are each derived in such a way that require assumptions in the posing of the problem.

Conventional Aerospace Method Challenges: Conventional state-of-the-art aerospace methods use series expansions of solutions to the Laplace Equation, given by Equation 10, which is solved by separation of variables to determine orthonormal basis functions which determine the spherical harmonics basis functions whose coefficients are conventionally given by Equation 11 (Takahashi & Scheeres, 2014), (Scheeres, 2016), (Zamaro & Biggs, 2015). This gives approximations for the gravitational potential $U^{e}$ outside the body given by Equation 12.

This solution method works best for the homogeneous differential equation (Equation 10) and is effective for solving for the gravity field of the body. To calculate the potential or force at points interior to the shape of interest, Poisson’s Equation is conventionally solved, as given by Equation 13.

Basis functions for the interior case have conventionally been found as Bessel function coefficients presented by Equation 14, though this work assumes that gravitational potential $V$ is proportional to the density distribution (Takahashi & Scheeres, 2014).

This state-of-the-art procedure for both the interior and exterior cases in this procedure are summarized by Algorithm B.1.

Main Sources of Error: There is a key region where the errors in current state-of-the-art methods arise in the calculation of the gravitational force field for any non-spherical body. Namely, assuming we have some mass $M$ with non-negligible volume $V$, a Brillouin sphere is any sphere centered at the center of $M$ which entirely encloses $M$. Within the Brillouin sphere, conventional gravitational force predictions are generally demonstrated to be inaccurate and divergent (Costin et al., 2020), (Šprlák & Han, 2021). Even for “nearly spherical” bodies, the divergence of the predication proves to be significant; historically, state-of-the-art methods give errors on the order of tens of percent for calculations within the minimum Brillouin sphere of the body (Takahashi & Scheeres, 2014), (Scheeres, 2016), (Scheeres et al., 2019). See Figure 11 for a comparison of these methods with our method of direct integration, about a uniform mass sphere.

Spherical harmonics expansions are a standard approximation to the gravitational field far from the body, but, in fact, convergence of spherical harmonics occurs with probability zero inside the minimum Brillouin sphere (Costin et al., 2020). Even for bodies such as the moon, where spherical harmonics expansions and other smooth approximations have been used throughout literature, direct surface measurements taken by GRAIL (Gravity Recovery and Interior Laboratory) indicate divergences and significant inaccuracy to conventional aerospace calculations within the minimum Brillouin sphere (Šprlák & Han, 2021). Our problem lies in that the “cusp-like” characteristics of the true force function make smooth approximations highly inaccurate near the surface; that is, the true force function is riddled with cusps and non-differentiable points as was shown in Figure 1 for the gravitational force field near a sphere. Algorithm B.1: Gravitational Potential $U(r,\lambda,\phi)$ Computation (State-of-the-art, Takahashi & Scheeres (2014))
———————————————————————————————————— 1. Eigenvalues $\alpha^{i}_{ln}$. Compute interior (Bessel) gravity field eigenvalues $\alpha^{i}_{ln}$ via boundary condition $j_{n-1}(\alpha^{i}_{ln})=0$ upon defining $j_{n}$ as the spherical Bessel function of the first kind, letting $n$ be the degree and $l$ be the power. Refer to Table 1 for pre-computed values. 2. Normalization factors $\mathscr{E}^{i}_{q}$, constants $G$, $M$, $R^{*}_{e}$. Working in spherical coordinates, define $\mathscr{E}_{q}^{i}(\alpha_{ln}^{i})=\alpha_{ln}^{i2}\left[j_{q}(\alpha_{ln}^{i})\right]^{2}$, latitude $\lambda$, longitude $\phi$, volume element $dv^{\prime}=r^{\prime 2}\sin\theta\,dr^{\prime}\,d\phi^{\prime}\,d\lambda^{\prime}$ via integration over points $(r^{\prime},\lambda^{\prime},\phi^{\prime})$ with nonzero density, and colatitude $\theta=\pi/2-\phi$. Define constants $M$ to be reference mass of the body, $R^{*}_{e}$ the reference radius of the body, and $G$ the gravitational constant. Define $P_{nm}$ to be the Legendre polynomials of order n and degree m. 3. Interior gravity field basis coefficients $\mathscr{A}^{i}_{lnm}$, $\mathscr{B}^{i}_{lnm}$. Compute interior (Bessel) gravity field non-dimensional coefficients $\mathscr{A}^{i}_{lnm}$, $\mathscr{B}^{i}_{lnm}$ via Equation (14) and scale them by the normalization factor $R^{*}/M^{*}G$. 4. Exterior gravity field basis coefficients $C^{e}_{nm}$, $S^{e}_{nm}$. Compute exterior gravity field non-dimensional coefficients $C^{e}_{nm}$, $S^{e}_{nm}$ via Equation (11). 5. Interior gravity field potential. Define interior gravity potential field $U^{i}$ as in Equation (15); using previous definitions and computations, this provides an analytical expression for $U^{i}$ as a function of any point $(r,\lambda,\phi)\epsilon\mathbf{R}^{3}$ up to $r<R^{*}_{e}$. 6. Exterior gravity field potential. Define exterior gravity potential field $U^{e}$ as in Equation (12) for an analytical expression of potential at any point $(r,\lambda,\phi)\epsilon\mathbf{R}^{3}$ such that $r\geq R^{*}_{e}$.