Fully analytical propagator for lunar satellite orbits in closed form

TL;DR

This work delivers a fully analytical, long-term propagator for lunar satellite orbits by combining a secular Hamiltonian framework with two closed-form element transformations: osculating→mean and mean→proper. The model uses a gravity description including 12 lunar harmonics and Earth's quadrupole tide with a refined lunicentric ephemeris, achieving accuracy comparable to semi-analytical approaches across altitudes $300\mathrm{km}$–$3000\mathrm{km}$ and decades-long horizons, while avoiding intermediate numerical propagation. The authors provide an open-source Python/symbolic implementation (legmaseft2025) and demonstrate precision tests against Cartesian propagation, elucidating error sources and the impact of secular resonances (notably the resonances $\mathcal{R}_{2,1,0,1,-2}$ and $\mathcal{R}_{2,2,-2,0,0}$) on small divisors. The approach is extensible to include more gravity terms, offers non-singular variable handling, and yields a practical, analytically tractable tool for planning and studying lunar satellite missions over multi-decade timescales.

Abstract

We present a fully analytical propagator for the orbits of lunar artificial satellites in a lunar gravity and third-body model sufficiently precise for a wide range of practical applications. The gravity model includes the twelve most important lunar gravity harmonics as well as the Earth's quadrupole tidal terms with a precise representation of the Earth's lunicentric ephemeris, and it gives an accuracy comparable to the way more extended semi-analytical propagator SELENA [6] for satellite orbits at altitudes from 300 to 3000 km. Extra terms of a more complete gravity model are straightforward to include using the formulas of the presently discussed analytical theory. The theory is based on deriving an approximate analytical solution of the secular part of the equations of motion using a Hamiltonian normal form in closed form. In total, we have two types of element transformations: from osculating to mean elements (as in [6]), and from mean to proper elements. The solution of the problem in proper elements is trivial, and, through the inverses of the above transformations, it allows to recover the position and velocity of a satellite analytically at any time t given initial conditions of the osculating elements at time $t_0$ without any intermediate numerical propagation. The propagator model is valid in time spans of several decades, and for every initial condition leading to no-fall on the Moon's surface, except for identified thin zones around a set of secular resonances corresponding to commensurabilities between the satellite's secular frequencies and the secular frequencies of the lunicentric Earth's orbit. Open software python and symbolic routines implementing our propagator are provided in the repository [14]. Precision tests with respect to fully numerical orbital propagation in Cartesian coordinates are reported.

Figures (10)

• Figure 1: Maps of the maximum value of $\log_{10}(|\Delta\underline{e}(t)|)$ in the timespan $0\leq t\leq 178$ days, where $\Delta\underline{e}(t)=\underline{e}_a(t)-\underline{e}_n(t)$, with $\underline{e}_{a}(t)=$ analytically propagated, and $\underline{e}_{n}(t)=$ numerically propagated (with the Hamiltonian $\mathcal{H}_{sec}$, Eq. \ref{['hamsm_simply']}) mean eccentricity of an orbit. The value of $\log_{10}(|\Delta\underline{e}(t)|)$ is given in color scale as indicated in each figure. The map is computed over a $100\times 100$ grid of initial conditions in $\underline{e}(0),\underline{i}(0)$ within the limits indicated in each panel, with the remaining mean elements being $\underline{\ell}(0)=0$, $\underline{g}(0)=-0.4$ rad, $\underline{h}(0)= 0.7$ rad, and $\underline{a}(0)=a_p=const=R_{\leftmoon}+\delta a$, where $\delta a = (300,500,700,900,2000,3000)~$km. Orbits with error smaller than $10^{-3}$ are classified as of error $10^{-3}$, while those with error larger than $10^{-1}$ are classified as of error $10^{-1}$.
• Figure 2: Same as in Fig. \ref{['fig:errormapsecc']}, but for the maximum error $\Delta\underline{i}(t)=\underline{i}_a(t)-\underline{i}_n(t)$ in the analytical versus numerical propagation of the mean inclination of the trajectories. The error is given in radians, hence the allowance to use a logarithmic color scale.
• Figure 3: Analytical (black line) versus numerical (red line) propagation of the mean elements $(\underline{e}(t),\underline{i}(t))$ for four orbits with initial conditions $\underline{\ell}(0)=0$, $\underline{g}(0)=-0.4$ rad, $\underline{h}(0)=0.7$ rad, $\underline{e}(0)=0.1$, $\underline{i}(0)=15^\circ$ and (constant) semi-major axis corresponding to the group of panels $\underline{a}=R_{\leftmoon}+500$ km (top left), $\underline{a}=R_{\leftmoon}+700$ km (top right), $\underline{a}=R_{\leftmoon}+900$ km (bottom left), $\underline{a}=R_{\leftmoon}+2000$ km (bottom right). In each group of panels, the first row shows the evolution of the corresponding elements within the first year of propagation, the second row within the fifth year of propagation, and the third row shows the corresponding errors of the comparison of the analytical and numerical solution in $\log_{10}$ scale.
• Figure 4: The error in eccentricity $|\Delta\underline{e}(t)|$ of the analytical versus numerical propagation of the mean elements as a function of time, shown in linear scale, for the trajectories with initial mean elements $\underline{\ell}(0)=0$, $\underline{g}(0)=-0.4$ rad, $\underline{h}(0)=0.7$ rad, $\underline{e}(0)=0.1$, $\underline{i}(0)=15^\circ$, and (left) $\underline{a}=500$ km, (right) $\underline{a}=300$ km.
• Figure 5: Analytical (black line) versus numerical (red line) propagation of the mean elements $(\underline{e}(t),\underline{i}(t))$ for two orbits with groups of panels similar as in Fig. \ref{['fig:orbitsgood']} and initial conditions $\underline{\ell}(0)=0$, $\underline{g}(0)=-0.4$ rad, $\underline{h}(0)=0.7$ rad and (left) $\underline{a}=R_{\leftmoon}+300$ km, $\underline{e}(0)=0.1$, $\underline{i}(0)=15^\circ$, (right) $\underline{a}=R_{\leftmoon}+3000$ km, $\underline{e}(0)=0.35$, $\underline{i}(0)=40^\circ$.