Fast and Efficient Lunar Finite Element Gravity Model

Abstract

In this paper, the finite element method (FEM) is integrated with orthogonal polynomial approximation in high-dimensional spaces to innovatively model the Moon’s surface gravity anomaly. The aim is to approximate solutions to Laplace’s classical differential equations of gravity, employing classical Chebyshev polynomials as basis functions. Using classical Chebyshev polynomials as basis functions, the least-squares approximation was used to approximate discrete samples of the approximation function. These test functions provide an understanding of errors in approximation and corresponding errors due to differentiation and integration. These test functions provide an understanding of errors in approximation and corresponding errors due to differentiation and integration. The first application of this project is to substitute the globally valid classical spherical harmonic series of approximations with locally valid series of orthogonal polynomial approximations (i.e., using the FEM approach). With an error tolerance set at ${10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$, this method is used to adapt the gravity model radially upwards from the lunar surface. The results showcase a need for a higher degree of approximation on and near the lunar surface, with the necessity decreasing as the radius increases. Notably, this method achieves a computational speedup of five orders of magnitude when applying the method to radial adaptation. More intrinsically, the second application involves using the methodology as an effective tool in solving boundary value problems. Specifically, this approach is implemented to solve classical differential equations involved with high-precision, long-term orbit propagation. This application provides a four-order-of-magnitude speedup in computational time while maintaining an error within the ${10}^{-10}$$\mathrm{m}{\mathrm{s}}^{-2}$ error range for various orbit propagation tests. Alongside the advancements in orthogonal approximation theory, the FEM enables revolutionary speedups in orbit propagation without compromising accuracy.

1. Introduction

On 4 March 2022, at approximately 4:25 a.m. PST, a rocket component was reported to have crashed into the far side of the Moon [1]. Although the specific rocket component is still not verified, losing control over spacecraft during lunar orbit missions has been an ongoing issue since the launch of the Luna 1 rover in 1959. Moreover, scientists have noticed a peculiar phenomenon about the lunar surface gravity; as satellites and probes orbited the Moon over certain craters and impact basins, they periodically veered off course and sometimes even plummeted toward the lunar surface. This is a particular problem when observing these effects on the far side of the Moon (See Figure 1). It was not until the Apollo 15 and 16 missions that scientists described the Moon’s surface gravity as “lumpy” due to the differing mass concentrations caused by the geology. These fluctuations in the lunar surface gravity have always puzzled scientists, which led to a high demand for precise lunar gravity data in the hopes of avoiding these crashes. Additionally, from the Lunar Reconnaissance Orbiter (LRO) [2] to the recent Artemis program [3], the field of lunar-based orbit propagation research has been increasingly growing over the decades. As a result of this, the importance of modeling the lunar surface gravity anomaly to aid in future Moon missions has become crucial in predicting the trajectories of spacecraft orbiting around the lunar surface. Using high-resolution gravity data, many lunar surface gravity models, such as NASA’s Gravity Recovery and Interior Laboratory’s (GRAIL) GRGM1200B model [4] and the GSFC LP-100k lunar potential Model [5], have been created as an attempt to constrain lunar surface gravity anomalies in response to this demand. Although sufficient enough to model the lunar surface gravity, these existing models rely on the global spherical harmonic (SH) series to solve Laplace’s equation for gravity.

The SH method for solving the equation of gravity was initially implemented as a way to model the geopotential through the use of Earth-based satellites recognized in the late 1950s [7]. Recent projects that implement gravitational acceleration modeling about the Moon can be viewed in the following literature: [8,9,10]. However, the SH approach proves to be computationally demanding, with three main issues arising [11,12,13]:

• At polynomial orders $n>2$, convergence proves to be inefficient and slow, which requires tens of thousands of terms to obtain accuracy for higher-order gravity representation.
• In any astronomical surface gravity modeling, the north and south poles represent singularities in SH, which result in poor approximations at the poles. For further examinations on obstacles in SH gravity approximations, see [14,15,16]
• SH is only accurate with spherical-shaped astronomical bodies. This means that the gravity of objects such as comets and oblate bodies cannot be modeled using this method.

The primary aim is to introduce a Finite Element Model (FEM) local gravity representation of the higher-order perturbation to offset the issue of slow convergence of the gravity models at polynomial orders greater than 2. Specifically, lower-order functions are used to sufficiently model the gravity perturbations about the lunar surface. This is highly motivated by the literature, where the projects propose similar advancements to the classical SH approach, such as in Junkins [17] and, more importantly, in recent investigations by the following: [18,19,20].

To establish an efficient resolution to the above computational challenges, this paper details a new approach to modeling surface gravity anomaly, increasing computational speed while maintaining a conservative error. Specifically, the paper proposes replacing the current LP-100k lunar SH gravity model with a series of local orthogonal polynomial approximations implemented by the finite element method of modeling lunar surface gravity. In Bani Younes’s dissertation [21], a FEM surface gravity model was implemented for the geopotential, which holds significant relevance to the current applications on the lunar FEM model. With this paper, steps have been taken to produce a cost-effective and accurate FEM lunar surface gravity model that mirrors that of the paper. Furthermore, parallelization of the lunar FEM surface model is implemented, building upon the research of [22], where a parallel evaluation of the Chebyshev approximation was studied. Parallelization of the original FEM lunar model and astrodynamics applications using the FEM model further encourages computational speedup in both regards. With this, creating a FEM lunar surface gravity model elevates the surface gravity modeling field with these main contributions:

• Computational efficiency;
• Memory allocation;
• Precision;
• Parallelization.

The following sections present the key background related to the proposed method and applications. First, the FEM is reviewed based on the Chebyshev polynomial interpolation. Then, the orbit propagation is adopted to integrate the approximated force function (i.e., FEM) to generate the satellite trajectory.

1.1. Finite Element Method

This section presents an overview on the research attempts based on finite element methods corresponding to gravitational modeling. There have been a variety of approximation techniques and basis functions implemented for gravity representation, which includes weighted functions [12,23,24], wavelets [25], splines [26], octrees [27], pseudocenters [28] and 3D digital modeling [29]. FEM can be used in a wide variety of applications that are both interesting and important in the fields of mathematical modeling, engineering, and computational analysis. Some examples include using FEM to analyze electric sail static structures [30], determination of the gravity potential and tidal stresses of asteroids [31], and using FEM to analyze a lunar lander during a soft landing [32]. Similar to our approach, each method balances accuracy with efforts to encourage computational speed. FEM was first proposed by Junkins [17] as a way to interpolate the geopotential of Earth-based gravitational modeling. The work in Bani Younes [21] provided the computational framework to implement FEM to model the geopotential, acceleration due to gravity, and radial adaption about the Earth’s surface. With this descent, the purpose of this paper is to interpolate the lunar gravitational acceleration about the surface and above the surface to help aid future lunar-based spacecraft missions that require precise orbital propagation.

Using the lunar FEM model, truncating the classical expansion at $n=2$ allows for local gravity representation at the higher-order perturbations. The FEM modeling approach is applicable to any spherical model and even to irregularly shaped astronomical bodies, while SH is purely for spherical bodies. Furthermore, this method increases computation by trading computer memory for run-time speed.

For our purposes, FEM is used on the lunar surface, where the method subdivides the entire lunar surface as a large system and subdivides it into smaller finite elements. The gravitational differential equations will then be solved in each shell and will then be connected by nodes. However, determining the correct shell iteration for these higher-order approximations is key with this method.

1.2. Chebyshev Polynomial Interpolation

Developing a more efficient theory to solve higher-order Laplace’s equation for gravity requires a structure that is both radially adaptive and efficient. If implemented correctly, interpolation at any radial or point on the surface can be conducted. One approach involves using the Chebyshev polynomial method. Created by mathematician Pafnuty Chebyshev in the nineteenth century, Chebyshev polynomials are two sequences of polynomials related to the cosine and sine functions denoted as $T_n\left(x\right)$ and $U_n\left(x\right)$, respectively. The polynomials can be defined starting with trigonometric functions where the Chebyshev polynomials of the first kind ($T_n$) are given by:

$$T_n\left(cos\theta \right)=cos\left(n\theta \right)$$

(1)

and the polynomials of the second kind ($U_n$) are:

$$U_n\left(cos\theta \right)sin\left(\theta \right)=sin\left((n+1)\theta \right)$$

(2)

in which n is any integer number that can be used to convert these definitions into a polynomial form. Chebyshev polynomials are orthogonal polynomials that are solutions of a special kind of Sturm–Liouville differential equation called a Chebyshev differential equation. Most importantly, the roots of $T_n\left(x\right)$, i.e., the Chebyshev nodes, are used as estimation points. Approximating at this point produces efficient polynomial interpolation, which is useful for precise radial adaptation and other applications which require higher-order approximation.

Although typically used as a numerical method for solving initial value problems (IVPs) and certain classes of boundary value problems (BVPs), this method can also be used to approximate solutions of linear or non-linear ordinary differential equations (ODEs) like that of Laplace’s equations for gravity. Using these Chebyshev polynomials, an orthogonal approximation can be placed for the FEM modeling process in which the path approximations would encourage computational speedup.

1.3. Orbit Propagation

The last application in this project includes using the FEM model to approximate a satellite’s orbital trajectory over a provided time span. Similar works which inspired this method of orbital propagation can be found in the following literature: [33,34,35,36,37]. This application will be achieved using the resulting Chebyshev polynomial coefficients from the original FEM routine at specific radii above the lunar surface to approximate the acceleration a satellite with negligible mass would experience during the course of a short- and long-term orbit around the Moon. This propagation considers the two-body problem, which consists of a satellite with negligible mass orbiting a larger astronomical body. This problem allows us to use Newton’s second law force equations to then come up with a finalized two-body equation of motion for the theoretical satellite orbiting the Moon. In this paper, first the orthogonal approximation techniques used alongside Chebyshev polynomials are covered in detail in Section 2. Then, in Section 3, the process of utilizing FEM on the LP-100k model is discussed and compared to the original SH model. In Section 4, a discussion is presented on how the coefficients produced by the FEM model on a two-body orbit propagation problem about the lunar surface are implemented. Lastly, in Section 5, a summary of the findings and concluding remarks is provided.

2. Orthogonal Approximation

The following are some examples of treatments of discrete approximation using Chebyshev polynomials: [24,38,39,40,41,42,43,44,45]. Although [41,42] hold more significance in the field of orthogonal approximation, in [44], the orthogonal approximation is used to show that multi-resolution approximation of input and output maps can be linear and nonlinear. Typically, a linear least-squares approximation is used for an orthogonal approximation when dealing with experimental data or an arbitrary function generating polynomials. In this paper, a least-squares approximation with classical Chebyshev polynomials is used to analyze discrete samples of the “to-be” approximated function. Furthermore, arranging the regression matrix to be “Kronecker factorable” (a technique used in [21]) allows the use of array algebra identities to generate the multidimensional orthogonal least-squares operators. This approach can be seamlessly adapted to higher-dimensional approximations. This section summarizes how using the Kronecker product technique extends a single classical Chebyshev polynomial approximation result to higher dimensions. This method avoids dimensional challenges related to how many attributes a dataset has and paves the way for efficient computation. The importance of this method allows considering the advantages of orthogonal approximation in gravity modeling. Most mathematical derivations are based on the work presented in [21].

First, let us briefly consider the following approximation of a single-valued function, $f\left(x\right)$, with x as one independent variable:

$$g\left(x\right),\left\{x_{min}\le x\le x_{max}\right\}$$

(3)

For simplification in multi-dimensional analysis, a new independent variable is introduced $\xi$, where $\{-1\le \xi \le 1\}$. Thus, through the process of normalization, the forward and inverse transformation for $\xi$ can be produced:

$$\begin{array}{c}Forward\left\{\begin{array}{c}\xi \left(x\right)=2(x-x_{min})/(x_{max}-x_{min})-1\end{array}\right.\\ \\ Inverse\left\{\begin{array}{c}x\left(\xi \right)=x_{min}+(\xi +1)(x_{max}-x_{min})/2\end{array}\right.\end{array}$$

(4)

The final function to be approximated can be created by substituting the forward and inverse functions into Equation (3):

$$f\left(\xi \right)=g\left(x_{min}+(\xi +1)(x_{max}-x_{min})/2\right)$$

(5)

In a general case, for basis functions, ${\varphi }_n\left(\xi \right)$, we seek to approximate $f\left(\xi \right)$ as a linear combination of a set of $N+1$ linearly independent basis functions $\{{\varphi }_0\left(\xi \right),{\varphi }_1\left(\xi \right),\dots ,{\varphi }_N\left(\xi \right\}$ as the following:

$$f\left(\xi \right)\approx \sum _{n=0}^N{a}_n{\varphi }_n\left(\xi \right)$$

(6)

For discrete cases, a set of sample points is introduced in which nodes is set as $\{{\xi }_0,{\xi }_i,\dots ,{\xi }_M;M\ge N\}$. With this, the residual approximation error at each node is considered to be:

$$r_j=f\left({\xi }_j\right)-\sum _{n=0}^N{a}_n{\varphi }_n\left(\xi \right);j=0,1,\dots ,M$$

(7)

In general, the decision to use a least-squares approximation routine is motivated by the advantages of applying the one-dimensional analysis to higher dimensions. Using least-squares approximation, the objective is to find the coefficient vector a that minimizes the weighted sum square of the residuals:

$$J={\displaystyle \frac{1}{2}}{(f-\Phi a)}^{T}W(f-\Phi a);W=W^T$$

(8)

where $\Phi ={\sum }_{n=0}^N{a}_n{\varphi }_n\left(\xi \right)$, and W is a positive definite weighted matrix determined by classical orthogonal conditions for Chebyshev polynomials.

Using a weighted least-squares solution and orthogonal conditions, this result can be displayed as normal equations:

$$a={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}W\mathbf{f}=C_x\mathbf{f}$$

(9)

This finalized relationship can be displayed in matrix form as the following:

$$f=\left[\begin{array}{c}f\left({\xi }_0\right)\\ f\left({\xi }_i\right)\\ ⋮\\ f\left({\xi }_N\right)\end{array}\right],C_x={\displaystyle \frac{1}{M_x}}\left[\begin{array}{ccccc}{\displaystyle \frac{T_0\left({\xi }_0\right)}{2}}& T_0\left({\xi }_1\right)& \cdots & T_0\left({\xi }_{M_{x-1}}\right)& {\displaystyle \frac{T_0\left({\xi }_{M_x}\right)}{2}}\\ \\ T_1\left({\xi }_0\right)& 2T_1\left({\xi }_1\right)& \cdots & 2T_1\left({\xi }_{M_{x-1}}\right)& T_1\left({\xi }_{M_x}\right)\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \\ T_{N_x-1}\left({\xi }_0\right)& 2T_{N_x-1}\left({\xi }_1\right)& \cdots & 2T_{N_x-1}\left({\xi }_{M_{x-1}}\right)& T_{N_x-1}\left({\xi }_{M_x}\right)\\ \\ T_{N_x}\left({\xi }_0\right)& 2T_{N_x}\left({\xi }_1\right)& \cdots & 2T_{N_x}\left({\xi }_{M_{x-1}}\right)& T_{N_x}\left({\xi }_{M_x}\right)\end{array}\right],a=\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_N\end{array}\right]$$

Similarly, for more than one independent variable, let us consider the function $f(x,y)$, with the two independent variables x and y:

$$g(x,y),\left\{x_{min}<x<x_{max}\right\},\left\{y_{min}<y<y_{max}\right\}$$

(10)

As in the one-variable approximation example, let us use a set of new independent variables, $\xi$ and $\eta$, to reduce the computational burden of dimensionality which leads to the normalized forward and inverse functions:

$$\begin{array}{c}Forward\left\{\begin{array}{c}\xi \left(x\right)=-1+2(x-x_{min})/(x_{max}-x_{min})\\ \\ \eta \left(x\right)=-1+2(y-y_{min})/(y_{max}-y_{min})\end{array}\right.\\ \\ Inverse\left\{\begin{array}{c}x\left(\xi \right)=x_{min}+(\xi +1)(x_{max}-x_{min})/2\\ \\ y\left(\eta \right)=y_{min}+(\xi +1)(y_{max}-y_{min})/2\end{array}\right.\end{array}$$

(11)

Substituting these equations into Equation (11) provides us with the function to be approximated:

$$f(\xi ,\eta )=g\left(x\left(\xi \right),y\left(\eta \right)\right)=g\left(x_{\min }+(\xi +1){\displaystyle \frac{(x_{\max }-x_{\min })}{2}},y_{\min }+(\eta +1){\displaystyle \frac{(y_{\max }-y_{\min })}{2}}\right)$$

(12)

and in the discrete measurement case tested at sample nodes, $\{{\xi }_0,{\xi }_1,\dots ,{\xi }_{M_x};M_x\ge N_x\}$, $\{{\eta }_0,{\eta }_1,\dots ,{\eta }_{M_y};M_y\ge N_y\}$, the residual approximation error at each node is:

$$r_{ij}=f({\xi }_i,{\eta }_j)-\sum _{\alpha =0}^{N_x}\sum _{\beta =0}^{N_y}{a}_{\alpha \beta }{\varphi }_{\alpha \beta }({\xi }_i,{\eta }_j);i=0,1,\dots ,M_x,j=0,1,\dots ,M_y$$

(13)

which can be displayed in vector-matrix notation: $r=f-\Phi a$ where $\Phi =\sum a_{ij}$. When in matrix form, a common issue that needs to be considered is dimensionality. To avoid this problem, efficient choices of basis functions and associated sample points must be made. This technique is examplified in [46]. Referring to section 1.6.2, $\Phi$ is formed as the Kronecker product of two basic matrices:

$$\Phi =\left[\begin{array}{ccc}1& {\xi }_0& {\xi }_0^2\\ \\ 1& {\xi }_1& {\xi }_1^2\\ \\ 1& {\xi }_2& {\xi }_2^2\\ \\ 1& {\xi }_3& {\xi }_3^2\end{array}\right]\otimes \left[\begin{array}{cc}1& {\eta }_0\\ \\ 1& {\eta }_1\\ \\ 1& {\eta }_2\end{array}\right]={\Phi }_x\otimes {\Phi }_y$$

(14)

For our purposes, the Kronecker product, $C=A\otimes B$, is implemented in MATLAB [47] via the “kron” command. Subsequently, the weighted least-square solution matrix becomes:

$$a={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}f=\left\{\left[{\left({\Phi }_x^T{\Phi }_x\right)}^{-1}{\Phi }_x^T\right]\otimes \left[{\left({\Phi }_y^T{\Phi }_y\right)}^{-1}{\Phi }_y^T\right]\right\}f$$

(15)

As Equation (15) is “Kronecker favorable”, it can be further represented as $\Phi ={\Phi }_x\otimes {\Phi }_y$. Furthermore, this result can be generalized to higher dimensions as the following: $\Phi ={\Phi }_x\otimes {\Phi }_y\otimes {\Phi }_z$, which allows for us to generalize Equation (15) into a three-dimensional approximation space $(x,y,z)$ or $(\xi ,\eta ,\zeta )$:

$$a={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}f=\left\{\left[{\left({\Phi }_x^T{\Phi }_x\right)}^{-1}{\Phi }_x^T\right]\otimes \left[{\left({\Phi }_y^T{\Phi }_y\right)}^{-1}{\Phi }_y^T\right]\otimes \left[{\left({\Phi }_z^T{\Phi }_z\right)}^{-1}{\Phi }_z^T\right]\right\}f$$

(16)

Using the Kronecker product of three small least-squares operators to produce a larger operator allows us to approximate in higher dimensions by taking the cube root of the size of the three inverted matrices. Going up dimensions would require us only to introduce a fourth variable, which is the beauty of this Kronecker product trick. for more details about the properties of the Kronecker product and how lower-dimensional analysis can be easily “upgraded” to higher dimensions. With this, the above equation is generalized to n-dimensions by taking the nth root of the inverted matrix dimension:

$$a={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}Wf\equiv Cf\equiv {\left({\overline{\Phi }}^T\overline{\Phi }\right)}^{-1}{\overline{\Phi }}^T\overline{f}\equiv \overline{C}\overline{f}$$

(17)

where C is the Cholesky square root of $W^{1/2}$ of the weighted matrix $W=W^{1/2}{W}^{1/2}$. and $\overline{\Phi }\equiv W^{1/2}\Phi$, $\overline{f}=W^{1/2}f$ which makes $\Phi$ Kronecker factorable. This final result allows the same advantages to be used in the identity-weighted matrix case. With this, Equation (17) can be expressed as:

$$a={\left({\overline{\Phi }}^T\overline{\Phi }\right)}^{-1}{\overline{\Phi }}^T\overline{f}=\left\{\left[{\left({\overline{\Phi }}_x^T{\overline{\Phi }}_x\right)}^{-1}{\overline{\Phi }}_x^T\right]\otimes \left[{\left({\overline{\Phi }}_y^T{\overline{\Phi }}_y\right)}^{-1}{\overline{\Phi }}_y^T\right]\otimes \left[{\left({\overline{\Phi }}_z^T{\overline{\Phi }}_z\right)}^{-1}{\overline{\Phi }}_z^T\right]\right\}\overline{f}$$

(18)

where $\overline{\Phi }\equiv W^{1/2}\Phi$, ${\overline{\Phi }}_x\equiv W^{1/2}{\Phi }_x$, ${\overline{\Phi }}_y\equiv W^{1/2}{\Phi }_y$, ${\overline{\Phi }}_z\equiv W^{1/2}{\Phi }_z$, and $\overline{f}\equiv W^{1/2}f$.

It is useful to note that $W,W_x,W_y,$ and $W_z$ have a Kronecker relationship: $W=W_x\otimes W_y\otimes W_z$, which allows for us to sample locations corresponding to the Chebyshev polynomials. The Kronecker product trick for constructing the two-dimensional Chebyshev least-squares matrix can be provided by:

$$\overline{C}\equiv {\left({\overline{\Phi }}^T\overline{\Phi }\right)}^{-1}{\overline{\Phi }}^T={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^T{W}^{1/2}=C{W}^{-1/2}\equiv {\overline{C}}_x\otimes {\overline{C}}_y$$

(19)

where $\overline{C_x}=\left({\overline{{\Phi }_x}}^T\overline{{\Phi }_x}\right)$. The most important takeaway from this step is that C is also Kronecker factorable, just as weighted matrix W is. From the above conclusions, this formulation is applied to n-dimensions. Take the following example for approximating a function in three dimensions:

$$\overline{C}={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}W=\overline{C_x}\otimes \overline{C_y}\otimes \overline{C_z},W\equiv W_x\otimes W_y\otimes W_z,\overline{f}=W^{1/2}f$$

(20)

For final generalized function, consider a function of n independent variables:

$$g(x_1,x_2,\dots ,x_n),\left\{x_{i,min}\le x_i\le x_{i,max}\right\},i=1,2,\dots ,n$$

(21)

with the following forward and inverse functions:

$$\begin{array}{c}Forward\left\{\begin{array}{c}{\xi }_i\left(x_i\right)=-1+2(x_i-x_{i,min})/(x_{i,max}-x_{i,min})\end{array}\right.\\ \\ Inverse\left\{\begin{array}{c}{x}_i\left({\xi }_i\right)=x_{i,min}+({\xi }_i+1)(x_{i,max}-x_{i,min})/2\end{array}\right.\end{array}$$

(22)

The function to be approximated is written by substituting Equations (22) into (21) using an independent variable (as performed before):

$$f({\xi }_1,{\xi }_2,\dots ,{\xi }_n)=f\left({\xi }_i\right)=g\left(x_{i,min}+({\xi }_i+1)(x_{i,max}-x_{i,min})/2\right)$$

(23)

Similarly, for discrete measurements, a set of sample points (nodes) is introduced: $\{{\xi }_{i,0},{\xi }_{i,i},\dots ,{\xi }_{i,Mx,i};M_{xi}\ge N_{xi}\}$. This produces the following residual approximation error at each measurement node:

$$r_{i_k}=f({\xi }_{1_k},{\xi }_{2_k},\dots ,{\xi }_{n_k})-\sum _{{\alpha }_1=0}^{N_{x_1}}\sum _{{\alpha }_2=0}^{N_{x_2}}\cdots \sum _{{\alpha }_n=0}^{N_{x_n}}{a}_{{\alpha }_1{\alpha }_1\dots {\alpha }_n}{\varphi }_{{\alpha }_1{\alpha }_1\dots {\alpha }_n}({\xi }_{1_k},{\xi }_{2_k},\dots ,{\xi }_{n_k})$$

(24)

where the k indices represent the measurement for the ith variable. Equation (24) can written in vector-matrix notation, $r=f-\Phi a$. Finally, as mentioned in the two- and three-dimensional examples, the weighted least-squares solution can be written as the Kronecker product of n smaller least-squares matrices:

$$a={({\Phi }^{T}W\Phi )}^{-1}{\Phi }^{T}Wf\equiv \overline{C_{x1}}\otimes \overline{C_{x2}}\otimes \dots \otimes \overline{C_{xn}}{W}^{1/2}f$$

(25)

where $\overline{C_{xi}}=C_{xi}{W}_{xi}^{-1/2}$, $W=W_{x1}\otimes W_{x2}\otimes \dots \otimes W_{xn}$, and $a=vec\left\{a_{1k,2k,\dots nk}\right\}$, and $f=vec\left\{f_{1k,2k,\dots nk}\right\}$.

To conclude this section, it is good to note that the previously discussed methods of orthogonal approximation, i.e., the use of least-squares approximation with Chebyshev polynomials alongside the Kronecker product trick, is implemented in the FEM approximation routine, specifically for the three-dimensional case in approximating Laplace’s gravitational differential equations. Figure 2 and Figure 3 summarize the steps taken in the approximation routine for the one-, two-, and n-dimension approximation routine.

3. Representations of Lunar Potential Using Orthogonal Finite Elements

In the following discussion, we consider a basic orthogonal FEM approximation to a gravitational acceleration model determined from NASA’s Lunar Prospector (LP) mission in 100 spherical harmonic degrees [48]. Released in 2000, the LP Gravity model allows considerations of an analogous FEM approximation of the ‘LP100K’ gravity model. For high-precision orbit propagation, all terms used in the $(m,n)=(100,100)$ LP model are considered. Using this model, Laplace’s equation is solved for gravity using FEM and compared to the classical spherical harmonic approximations. Initially, a test was conducted to produce a contour plot of the disturbance gravity acceleration on a small degree order using the FEM approach. This test ensured a reasonable error value on top of a reasonable computation time, which encouraged us to continue the procedure with a higher degree order for an accurate representation. Figure 4 shows the radial disturbance acceleration on the lunar surface from the FEM representation with $(4\times 4)$ degree square using the N = 10 approximation.

As mentioned in the introduction, the three main challenges faced when using SH for orthogonal approximations are (1) an upper limit discrepancy, (2) inefficiency at orders higher than 2, and (3) non-free singularities at the poles of the astronomical bodies. With these three challenges in mind, a finite element method model of surface gravity representations is considered in the hopes that lower degree functions produced can be implemented to effectively and sufficiently model and compute lunar gravity. Furthermore, this method can be parallelized to increase computational efficiency. With this, the FEM model is structured to be both radially adaptive and efficient, which will result in better interpolation at higher degrees (i.e., at different radial distances above the surface of the Moon).

The following gravitational potential model is made up of the reference and disturbance gravity terms where the reference term includes two-body effects on the lunar gravity ($U_{2B}$), while the disturbance term ($\Delta U(r,\lambda ,\varphi )$) includes any other effects on gravity for the higher degree and order terms to take place. With these two terms, the finalized potential, $U(r,\lambda ,\varphi )$, can be given as:

$$U(r,\lambda ,\varphi )=\underset{ReferenceTerm}{\underbrace{U_{2B}(r,\lambda ,\varphi )}}+\underset{DisturbanceTerm}{\underbrace{\Delta U(r,\lambda ,\varphi )}}$$

(26)

Furthermore, the acceleration can be derived from the potential and is represented as:

$$\ddot{r}=-{\displaystyle \frac{\partial U}{\partial r}}=\underset{ReferenceTerms}{\underbrace{{\displaystyle \frac{\partial }{\partial r}}\left\{{\displaystyle \frac{GM}{r}}\right\}}}-\underset{DisturbanceTerm}{\underbrace{{\displaystyle \frac{\partial }{\partial r}}\Delta U}}$$

(27)

It is important to note that adjusting the model to correct for the two-body interference is performed before the reference potential compact and efficiency are considered. This means only the perturbing acceleration terms are considered in the model.

3.1. Sampling: Radii and Shells

To help illustrate a sample FEM grid, consider a spherical grid with radius $R_m$ (i.e., a sphere with the radius of the Moon at the lunar equator). This sphere is then covered by a 2-D mesh from $0\le \lambda \le {360}^{\circ }$ and $-{88}^{\circ }\le \varphi \le {88}^{\circ }$. It is good to note that the poles are neglected in this routine due to the limitations of the original SH model. As discussed, the poles of astronomical bodies are treated as singularities in SH modeling. However, the FEM model can be used to approximate the surface gravity at the poles using the Cartesian coordinates. This is implemented in later approximation models. The following section allows us to discuss the process of finding the most efficient initial conditions like shell size and radius sample.

3.1.1. Radial Sampling

At an arbitrary radius, r in the range of $R_m\le r\le 1.11R_m$, let the gravity data $U(r,\lambda ,\varphi )$ be provided at a constant radius, r, on a specific spherical shell. These data on this shell should be transformed into $U(\xi ,\zeta ,\eta )$, where $-1\le \zeta \le +1$ and $-1\le \eta \le +1$; we can retrieve an accurate gravity modeling representation of $U(\zeta ,\eta )$ on “sufficiently-dense” sets of concentric spherical surfaces. However, sampling radial coordinates is key when considering this unique variation.

To determine a proper radial sample, two different sampling methods are considered: uniform sampling, approximating the gravity at every radius between the lunar surface to 200 km above in 10 km increments, and cosine sampling. For cosine sampling, consider the transformed position r in which the position follows a “cosine-like” transformation as a function of the previously mentioned variable $\xi$ as the following:

$$r=r_{min}+(r_{max}-r_{min})\left[1-cos\left({\displaystyle \frac{\pi }{4}}(1+\xi )\right)\right]$$

(28)

where

$$\xi ={\displaystyle \frac{4}{\pi }}co{s}^{-1}\left(1-{\displaystyle \frac{r-r_{min}}{r_{max}-r_{min}}}\right)-1$$

(29)

where $-1\le \xi \le +1$ and $r_{min}\le r\le r_{max}$, $r_{min}=R_m$ and $r_{max}=1.11R_m$. At higher radii above the lunar surface, gravity anomalies decline rapidly. With this fact, we seek to format an efficient way to sample surface gravity that allows memory allocation. With cosine sampling, this transformation “naturally” constrains dense measurements at smaller radii, ensuring that gravity perturbations are maximized and mimic how gravity would naturally act. This means that no matter the number of points N at lower radii, these nodes then bunch up in comparison to larger radii. This idea is illustrated in Figure 5 for 10 and 21 sample points; compare this method to uniform sampling. To compact this progression, the transformed radius variable $\xi$ was created to help treat r as an independent variable. It is again important to note that $\xi$ is sampled non-uniformly using another cosine transformation in order to satisfy the Chebyshev orthogonality conditions:

$${\xi }_i=-cos(i\pi /2M),i=0,1,2,\dots ,M$$

(30)

Not only does cosine sampling simulate a natural progression in gravity, but this sampling method also produces the error magnitude on the order of (${10}^{-9}$) while taking the same amount of computational power as uniform sampling. For these reasons, we continue using cosine sampling in this analysis at a latitude and longitude of $\lambda$ = $\varphi$ = $0^{\circ }$ degrees on the Moon.

3.1.2. Shells and Polynomial Order

After determining a proper radial sample, determining the required polynomial order as a function of radius is the next step. The objective is to determine the order that adaptively maintains an efficient error tolerance. The required acceleration error tolerance is determined from the polynomial order N at every radius. For example, if a 4-by-4-degree square area on the lunar surface is considered, the convergence error is the maximum absolute error between the LP SH model’s acceleration values and the approximated acceleration from the FEM method.

Consider a shell size of $(10\times 10)$ square degrees compared to a smaller size of $(2\times 2)$ square degrees. The $(2\times 2)$ shell size would require less computational power/time to compute all the coefficients over a small area, while the $(10\times 10)$ shell size would require a longer computational time to reach a similar error tolerance because of the larger number of coefficient values calculated. The compromise of time would be better justified with the resulting error produced. Although the larger shell size may take a longer period of time to produce results, if these results display a small error tolerance (${10}^{-9}\mathrm{m}{\mathrm{s}}^{-1}$), it would be best to use this. Thus, when determining the correct shell size for this model, it follows to pursue a shell size that maintains an efficient error tolerance while not compromising computational power.

With a maximum degree order of 100 and shell sizes of ($(4\times 4)$ and $(6\times 6)$, both sizes are elevated and checked for accuracy and computational time across a portion of the near side of the lunar surface (from $0\le \lambda \le {180}^{\circ }$ and $0^{\circ }\le \varphi \le {88}^{\circ }$). As discussed previously, the smaller shell size of $(4\times 4)$ would be more desirable if the computational time is reasonable. However, if the accuracy for the $(6\times 6)$ case is reasonable compared to the $(4\times 4)$ shell case, this shell size is considered. Figure 6 displays a graph of computational time and calculation error as a function of shell size. For a shell size of $(4\times 4)$, with maximum degree order of 100, it took the machine (see machine configurations in Section 4) approximately 3627.43 s to produce the relevant FEM approximation coefficients with an error of $6.16\times {10}^{-10}\mathrm{m}{\mathrm{s}}^{-1}$. For the $(6\times 6)$ case, it took the same machine approximately 10,653.11 s to produce the FEM approximation coefficients with an error of $1.40\times {10}^{-16}\mathrm{m}{\mathrm{s}}^{-1}$. With the shell size of $(6\times 6)$ providing such a small error compared to typical errors found in other FEM surface gravity models produced in [21], this error is insignificantly smaller given the increased computational time. With this, it is concluded that the $(4\times 4)$ case provides the best balance between computational time and reasonable accuracy.

3.2. Orthogonal Approximation of the Gravitational Acceleration

Using orthogonal approximation, the lunar gravitational acceleration is estimated by transforming gravity into sample points. These sample points are located in the distribution in the least-squares approximation. In this subsection, a brief discussion is given.

Let gravity data $U(r,\lambda ,\varphi )$ transform to $U({\xi }_i,{\zeta }_j,{\eta }_k)$ with sample points located as a cosine distribution. For the three components of gravity, the approximation is as follows:

$$\begin{array}{c}{G}_x({\xi }_i,{\zeta }_j,{\eta }_k)\cong \sum _{\alpha =0}^{N_x}\sum _{\beta =0}^{N_y}{a}_{\alpha \beta }\left({\xi }_i\right){\varphi }_{\alpha \beta }({\zeta }_j,{\eta }_k)={\Psi }_x^T({\zeta }_j,{\eta }_k)A_x\left({\xi }_i\right)\\ \\ G_y({\xi }_i,{\zeta }_j,{\eta }_k)\cong \sum _{\alpha =0}^{N_x}\sum _{\beta =0}^{N_y}{A}_{\alpha \beta }\left({\xi }_i\right){\varphi }_{\alpha \beta }({\zeta }_j,{\eta }_k)={\Psi }_y^T({\zeta }_j,{\eta }_k)a_y\left({\xi }_i\right)\\ \\ G_z({\xi }_i,{\zeta }_j,{\eta }_k)\cong \sum _{\alpha =0}^{N_x}\sum _{\beta =0}^{N_y}{a}_{\alpha \beta }\left({\xi }_i\right){\varphi }_{\alpha \beta }({\zeta }_j,{\eta }_k)={\Psi }_z^T({\zeta }_j,{\eta }_k)A_z\left({\xi }_i\right)\end{array}$$

(31)

where the indices $(i,j,k)$ range from 0 to $(M_r,M_{\lambda },M_{\varphi })$, and ${\Psi }_{x,y,z}$ are the basis functions (which are transposed) of the acceleration components.

Furthermore, Equation (32) can be expressed more compactly if the same basis function for the acceleration components uses $\Psi$ = ${\Psi }_x$ = ${\Psi }_y$ = ${\Psi }_z$, which provides us with the following relation: $\Psi ={\Psi }_x\otimes {\Psi }_y\otimes {\Psi }_z$. We then utilize the tensor product of $\Psi ({\eta }_j,{\zeta }_k)$ and the identity matrix, $I_n$, twice to “undo” the transpose acceleration components. This creates the finalized gravitational acceleration approximation equation:

$$G({\xi }_i,{\zeta }_j,{\eta }_k)=\left[\begin{array}{c}{G}_x({\xi }_i,{\zeta }_j,{\eta }_k)\\ \\ G_y({\xi }_i,{\zeta }_j,{\eta }_k)\\ \\ G_z({\xi }_i,{\zeta }_j,{\eta }_k)\end{array}\right]=[(\Psi ({\zeta }_j,{\eta }_k)\otimes I_n)\otimes I_n]\left[\begin{array}{c}{A}_x\left({\xi }_i\right)\\ A_y\left({\xi }_i\right)\\ A_z\left({\xi }_i\right)\end{array}\right]$$

(32)

where $\left[\begin{array}{c}{A}_x\left({\xi }_i\right)\\ A_y\left({\xi }_i\right)\\ A_z\left({\xi }_i\right)\end{array}\right]$ are the gravitational acceleration components in Cartesian representation. To use the FEM for computing gravitational acceleration, the maximum polynomial order is employed, ensuring consistency levels as compared to the overall model. In particular, at larger radii above the lunar surface, the computational speed decreases by an order of magnitude when only non-zero terms are included.

Parallelization

Recent developments have seen the parallelization of gravitational FEM models, as implemented in [22,49], utilizing the Open Multi-Processing (OpenMP) application interface, which supports shared-memory parallel programming in C/C++ and Fortran. In this project, the inputs are the spherical coordinates ($r,\varphi ,\lambda$) with corresponding normalized parameters $\xi$, $\zeta$, and $\eta$, which are used to evaluate the Chebyshev polynomials. The master thread is then split into slave threads, each evaluating the multivariate Chebyshev polynomials in parallel. These threads are then joined back together to be evaluated for the FEM acceleration by adding spherical and $J_2$ terms. A three-magnitude speedup was reported using OpenMP on a geopotential FEM model with eight threads. Even with a more efficient computation time, the parallelized FEM model still produced less than 10 digits of accuracy compared to the original SH model. With this in mind, OpenMP implementation to the FEM lunar model promises higher speedup gains if more threads are used. This study indicates potential efficiency gains in both memory allocation and computation time for future FEM lunar model developments.

3.3. Results of the FEM LP Analysis

The following discussion delves into the results derived from the above-mentioned FEM routine applying NASA’s LP100K SH model. The primary focus is on determining the best parameters for achieving an optimal balance between computational efficiency and error minimization.

For cosine sampling, the number of CGL nodes was highest near the surface and would decrease as the radii increased. The number of sample nodes should decrease as nodes are related to the number of coefficients used in the polynomial approximation at a specific radius. In order to find the best node quantity for each radius, the routine is iterated through uniform node values at each radius until a node value produces an error greater than a tolerance of ${10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$. Because the number of CGL nodes is related to the polynomial order, the goal is to further investigate by observing how the polynomial order affects computational time. With this, an analogous relationship is observed between CPU time and polynomial order (as expected). Figure 7 illustrates a decrease in time cost as a function of radius alongside the polynomial order decreasing in the same manner, confirming the idea that it takes a longer period of time to approximate at higher polynomial orders than at lower orders.

For higher-precision orbit propagation, all terms used in the $(m,n)=(100,100)$ LP spherical harmonic model are considered. As mentioned in the approximation theory, an error tolerance of ${10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$ is selected while substituting the high-degree order gravity model with a FEM approximation. Figure 8 showcases the maximum absolute error of the approximated norm of the “disturbance” acceleration as a function of the Chebyshev polynomial order N. Here, it is noted that as the radius increases, the error declines faster; a smaller polynomial order is sufficient enough to reach the placed error tolerance of ${10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$. However, in this plot, it is clear that there are two error tolerance values displayed. In the previous analysis, a polynomial order of 21 is identified that yielded an error of ${10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$. However, after reviewing the coefficient values in this polynomial order, it is found that the portion of corresponding errors was much smaller than allowed error tolerance. This suggests that this polynomial order was significantly implementing a high precision that was computationally expensive. With this, a “fetching” routine is implemented where any error values less than our error tolerance had their corresponding coefficient value to be approximated to 0.

As an outcome, most of the coefficient values were approximated to 0, further suggesting an excessive polynomial order. After revisiting the current analysis, and after the polynomial order of 18 was found for the surface, the same “fetching” routine was performed with $N=18$, finding that none of the error values strayed too far past the error tolerance, which was apparent in the weight of non-zero coefficient values found in that sample. This led us to believe that a maximum polynomial order of 18 was the most reasonable (computationally supportive and precise) option for the surface of the Moon. An error tolerance of ${10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$ proved to be too conservative. With this, a more liberal error tolerance value of ${10}^{-8}\mathrm{m}{\mathrm{s}}^{-2}$ was placed to ease the amount of computational burden while also still maintaining reasonable error.

As we get closer to the surface of the Moon, a higher polynomial order is required to reach this error tolerance. The (100, 100) LP model’s gravitational acceleration values are used as the truth and the errors between the FEM approximation are reduced by adjusting the degree value of the Chebyshev polynomials until it reaches the error tolerance. For the proposed lunar model, it is found that the most efficient maximum approximation error is ${10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$, which provides a fair balance between computational time and machine precision at a maximum polynomial order of N = 18 on the lunar surface with a $(4\times 4)$ square degree FEM cell size. However, when a polynomial order of 18 was applied to a radius 100 km above the lunar surface, the error produced was negligible. This confirmed the previous discovery about a decreasing radial polynomial order, i.e., as you go further from the surface, the polynomial order required to reach the error tolerance decreases as well. Furthermore, a maximum degree value (M) of 100 and a shell value of $(4\times 4)$ degree maintains an efficient error tolerance while constraining sensible computational power.

The determination of the correct polynomial order for each radial adaptation allowed for us to determine the number of coefficients required in approximation at each radius. At the lunar surface, the polynomial order of $N=18$ indicates that the gravity acceleration at the lunar surface can be approximated by ${(18+1)}^2=361$ orthogonal polynomial terms (or coefficients). At 200 km above the surface, a polynomial order of $N=6$ was found to be most efficient, which means it requires 49 coefficient terms. This suggests that as the polynomial order increases, the number of coefficients required in the approximation also increases as well. This relationship can also be supported when observing how an error is related to the node number at each radius. Figure 9 displays two plots: (a) polynomial order as a function of radius and (b) number of coefficients as a function of radius; these plots show a nearly similar result.

On the same machine, the FEM model displays three orders of magnitude of computational speedup at the lunar surface compared to the (100, 100) LP SH model. This can be seen in Figure 10. The FEM approach continues to be computationally efficient compared to the original (100, 100) SH model, even going as far as to decrease computational time as it reaches higher radii. At radii R near (150 km $+R_m)<R<$ (200 km $+R_m)$, it is found that the FEM approaches five orders of magnitude of computational efficiency compared to the SH model. If implemented for orbit propagation, the FEM gravity model should not only provide nine to ten digits of accuracy, but it also produces a three to five order of magnitude reduction in CPU time.

The Chebyshev coefficients were saved for each radial iteration. These first 49 coefficient terms are displayed for the final cosine sampling, as shown in Figure 11. With cosine sampling, higher-order non-zero coefficients have been found to be comparable to easy-to-fit smooth functions of similar nature. With this, it is observed that the norm of gravity decays rapidly, as it should. Figure 12 qualitatively displays a prompt decay of gravity anomalies in a finite element where the contours are created by the gravitational acceleration perturbation at three radii located at the surface, 100 km, and 200 km above the surface. These coefficients can then be used to approximate orbital propagation at different radii on or above the lunar surface. When considering the coefficients to be a function of radius, the two-dimensional approximations on each surface are used to promote three dimensions by representing each coefficient as an orthogonal Chebyshev function of radius. As a result, an average surface gravity acceleration value of $1.624\pm 7.626\times {10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$ was found for a polynomial order of 18 and at 200 km above the surface of the Moon, an average value of $0.5017\pm 1.36595\times {10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$. These results are illustrated in Figure 13 for all computed radii.

4. FEM Application in Astrodynamics: Satellite Trajectory Propagation

In this section, we discuss how the FEM can be applied to orbit trajectory propagation, and then we examine a case study to demonstrate how the coefficients produced from the FEM lunar model can aid in satellite trajectory propagation. Finally, we compare this result to an SH routine to analyze computational time and precision.

Modeling spacecraft trajectories has been one of the most important fields in astrodynamics and space guidance. Satellite orbit propagation predicts how an object moves based on its current state. The satellite’s position is defined by a group of approximate equations of motion where the degree of approximation is dependent on orbital information. A trajectory is the path a spacecraft follows through space as it is influenced by the gravitational force of one or more bodies, creating acceleration. The two-body model describes of a satellite with negligible mass orbiting a larger astronomical body. The equations of motion for this model can be solved, resulting in different conic sections depending on their distance away from the larger body. Although not represented in the two-body model, it is worth noting that when considering a satellite orbiting the Moon, we would need to consider the gravity of Earth. Introducing an “n-th” body example is to be investigated in future works.

In the equation of motion, the larger astronomical body is treated as fixed within the two-body frame. We then consider the body as a uniform point mass to further simplify the equation of motion. As stated previously, the spacecraft has negligible mass compared to the main body, which means its acceleration is only due to the gravity of the main body.

The following equation represents the two-body problem implemented using the FEM coefficients (represented by the disturbance acceleration):

$$\ddot{\overrightarrow{r}}=-{\displaystyle \frac{\mu }{r^3}}\overrightarrow{r}+\sum _{i=0}^k{G}_d^i$$

(33)

where $\mu$ is the gravitational constant of the Moon provided by the original ‘LP100K’ model: $\mu =4.902800238\times {10}^{12}\mathrm{m}{\mathrm{s}}^{-2}$ alongside the disturbance acceleration where $G_d^i$ is the ith coefficient of k gravity spherical harmonic perturbation terms. These terms were previously saved from the original FEM model.

As mentioned previously, the position of a spacecraft at any point in its orbit is approximated as long as we have at least one reference point to use as a starting approximation. Numerical methods of orbit propagation, also termed special perturbations, compute and approximate solutions of general equations of motion [50,51,52]. Figure 14 illustrates a general approach to satellite trajectory propagation provided initial conditions. In the following sections, it follows to implement the coefficients acquired from the FEM routine performed on a small patch of $(4\times 4)$ square degrees on various radii on and from the surface of the Moon to interpolate the trajectory of a point source satellite orbiting the Moon. Finally, the accuracy and computational time are checked with respect to that of the original LP (100 × 100) SH model. Orbital propagation using FEM has been implemented in a series of projects, some of which can be seen here: [53,54,55].

Trajectory Propagation Using FEM Versus LP (100 × 100) SH Gravity Model

In this section, the trajectory propagation results are compared using the coefficient values found from the FEM model to the classical SH method. Computations for the following examples were performed using a machine with the following configurations:

• Intel(R) Xeon(R) CPU 3.70 GHz, 16.0 GB of RAM;
• Windows 64-bit operating system, x64-based processor;
• MATLAB R2019a.

The final goal in this experiment is to justify that the coefficients produced by the original FEM routine can be used in various satellite orbit propagation problems. The steps needed for the orbit propagation experiment are the following:

• Produce a set of initial values: ${\varphi }_o$, ${\lambda }_o$, and R. These values represent the spacecraft’s azimuth, elevation, and radius from the lunar surface, respectively, at $t=0$.
• Set an “ending” position or time. This can be achieved by providing ${\varphi }_{final}$, ${\lambda }_{final}$ or providing a time span.
• Produce a small trajectory across a small shell. Using the coefficient values saved at a particular radius from the previous FEM routine to calculate the acceleration of the spacecraft along this initial trajectory using both the FEM and SH approaches. At this point, the resulting computational time and error between the two methods are then compared.
• Generate satellite positions. Using MATLAB’s “ode45” function to solve the differential two-body problem represented as Equation (33), we feed in the acceleration values calculated from step 3, and a time span to find the spacecraft’s position as a function of time.
• Generate satellite velocities. To achieve this, the MATLAB function “fsolve” is adopted to solve the systems of nonlinear equations of several values to solve for the satellite’s velocity over the course of the small shell provided initial velocities of $V_{x_0}$, $V_{y_0}$, $V_{z_0}$, the spacecraft’s position, and the time span.
• Generate a long-span satellite orbit trajectory. This is achieved, again, by utilizing the “ode45” function to generate satellite positions at a longer time span and the resulting final velocities.

Before determining the entire trajectory of a spacecraft, a smaller portion of its trajectory based on three different radii above the lunar surface (100 km, 150 km, and 200 km) is simulated first. This small trajectory spans over the (4 × 4) square degree shell at this radius. In this analysis, three tests on orbit propagation are considered. These tests differ in how the initial values are produced:

• Test 1: Randomized position changes. Using MATLAB’s “rand” function, we randomize a uniformly distributed pseudorandom set of values for $d\varphi$ and $d\lambda$, i.e., the changes in azimuth and elevation, respectively. Using these randomized values, two ranges for $\varphi$ and $\lambda$ are created:

  $$\begin{array}{c}[{\varphi }_{min}=0^o<\varphi +d\varphi <{\varphi }_{max}=4^o]\\ \\ [{\lambda }_{min}=0^o<\lambda +d\lambda <{\lambda }_{max}=4^o]\end{array}$$

  (34)
  With these simulated ranges and a radius 100 km above the lunar surface, these values are then converted into Cartesian coordinates to obtain starting x, y, and z positions of the simulated spacecraft.
• Test 2: A Monte Carlo simulation of initial conditions. With this approach, randomized values are repeated for $d\varphi$ and $d\lambda$ under the conditions that ${\varphi }_{min}<{\varphi }_{max}$ and ${\lambda }_{min}<{\lambda }_{max}$. This produces various amounts of ranges for $\varphi$ and $\lambda$ in which a randomly picked set of ranges is sought to achieve the initial positions of the simulated spacecraft.
• Test 3: Set starting and final positions. This approach is the most simple and straightforward, as it is reflective of the most basic problems in orbit propagation. It follows to start off with converting ${\varphi }_o=0^o$, ${\lambda }_o=0^o$, and R into the satellite’s initial Cartesian coordinate: $x_0$, $y_0$, and $z_0$. Unlike the previous two approaches, in this test, the satellite’s final position is found by converting ${\varphi }_{final}=4^o$, ${\lambda }_{final}=4^o$, and R into the satellite’s final Cartesian coordinate: $x_{final}$, $y_{final}$, and $z_{final}$. All that is left to do is to find the satellite’s positions during $t=0s$ and $t=t_{final}$.

For all three approaches, the initial ranges of $\varphi$ and $\lambda$ and R are used to generate satellite accelerations using the FEM routine and the LP $(100\times 100)$ SH routine. In the FEM analysis, the coefficients saved from the original model at $R=100\mathrm{km}+R_m$ were multiplied by the Chebyshev polynomials approximated at each elevation and azimuth. These results are then compared to SH results for computational time and error. Figure 15 displays the computational time for both approaches for each of the three tests for each of the three radii. For all three tests, FEM timing is three to four orders of magnitude smaller than that of SH.

For all three tests, the FEM model’s error is calculated by subtracting it from the “truth” value, specified before, to be the LP $(100\times 100)$ model. The resulting errors are illustrated in Figure 16. In Tests 1 and 2, the resulting FEM error from all three radii is in the order of ${10}^{-10}\mathrm{m}{\mathrm{s}}^{-2}$. With the previously mentioned increase in computational speedup, these two approaches produced above-adequate results for approximating the acceleration produced in the $(4\times 4)$ square degree shell 100 km above the lunar surface. However, for the Monte Carlo simulation (Test 2), there was a higher order of points in the iteration compared to Test 1, which is reflected in the computational time required for each of the tests. Figure 17 displays the initial trajectories created from these three tests over the $(4\times 4)$ square degree section. In Test 1, the trajectory is smooth and laminar, which is expected from a small spacecraft trajectory, while Test 2 displays a randomized path, which is representative of a typical Monte Carlo simulation. It is good to note that the trajectory found in Test 2 does not suggest that the model failed to produce a correct spacecraft orbit. This proves that the model was able to approximate a disorderly trajectory, as apparent in the small error.

The error found in Test 3 is in the order of ${10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$ across all three radii. It is also good to note that the error is lowest at the “end” points of the calculation due to the initial values set at each starting position. As illustrated in Figure 17c, the trajectory across the $(4\times 4)$ square degree shell is similar to that of the trajectory for Test 1; however, this test produced more data points.

Finally, a long space orbit propagation was carried out for all three tests. For this, a time span of 360 min was used. Figure 18 shows the long-span predicted trajectory of the satellite, as apparent from Test 3 for a radius of 100 km above the lunar surface as an example. This long-span trajectory is shown in red, while the initial small shell trajectory is in blue. The shape of this conic section is representative of typical trajectories found with two-body applications: [56,57,58].

This application of FEM shows promising results and demonstrates a potential path to significant advancements in computational speedup and precision for differential equations in astrodynamics.

5. Conclusions

This research introduced a new approach to modeling the Moon’s surface gravity anomaly, improving approximation methods for Laplace’s differential equations related to gravitational forces. The orthogonal FEM was used to approximate a gravitational acceleration model, utilizing data from NASA’s Lunar Prospector mission up to 100 spherical harmonic degrees. The FEM model was then evaluated against the traditional SH solution in terms of accuracy and computational cost. The findings demonstrated the efficiency of orthogonal approximation techniques, including Kronecker-factorable regression matrices for higher dimensions using algebraic identities on one-dimensional orthogonal least-squares operators. The use of Chebyshev polynomials enabled efficient selection of orthogonal basis functions, supporting an effective sampling mechanism. Two radial sampling techniques are also evaluated, maintaining a tolerance of ${10}^{-9}\mathrm{m}{\mathrm{s}}^{-2}$ on the resulting disturbance acceleration and saving the corresponding radially adaptive coefficient values spanning from the lunar surface to 200 km above. For instance, with a polynomial order N = 18, 361 orthogonal polynomial terms are deduced on the lunar surface. As a result, an average surface gravity acceleration value of $1.624\pm 7.626\times {10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$ was found for a polynomial order of 18 and at 200 km above the surface of the Moon, an average value of $0.5017\pm 1.36595\times {10}^{-9}$$\mathrm{m}{\mathrm{s}}^{-2}$. A significant highlight was the computational speedup achieved with the FEM approach. The results indicate that the FEM model offers a five-order-of-magnitude reduction in computational time compared to the original SH method. Consequently, there is a growing interest in applying these techniques to precise long-term orbit propagation for spacecraft and satellites in lunar missions. Three numerical orbital tests were conducted for a satellite orbiting the Moon. In all three tests, the coefficient values produced from the original FEM model were used to approximate the satellite’s trajectory as it goes across a (4 × 4) square shell 100 km above the surface. These accelerations showed small errors [${10}^{-6}–{10}^{-10}$ $\mathrm{m}{\mathrm{s}}^{-2}$] across all three tests with a four-order-of-magnitude speedup in computation time. A long-span trajectory was also achieved using FEM for one of the tests, displaying promising results in the interpolation of orbital mechanics representing the two-body problem. These findings are crucial for precise orbit propagation in various astrodynamics simulations. This research lays the groundwork for future advancements in orthogonal approximation theory and orbit propagation. Potential directions include code parallelization for computational efficiency, FEM-based machine learning benchmarks, n-body orbital propagation for satellites, and FEM modeling of the entire lunar surface for improved precision in gravitational disturbance. This approach could also be applied to model the surface gravity of other astronomical bodies like Mars, Europa, and comets. Overall, the method shows great promise for advancing astronomy and astrodynamics.