# Lunar Satellite Constellations in Frozen Low Orbits
## 2 Dynamical Model  
When numerically propagating low lunar orbits, the major perturbation stems from the complex lunar gravitational potential. It is very important to properly truncate it in order to avoid excessive computations while retaining high propagation accuracy. The most systematic approach to truncating the gravitational potential shares the philosophy of the perturbation theory: the magnitude of truncated terms relative to the central (Newtonian) part of the lunar gravitational acceleration should be bounded by some predefined small parameter \(\varepsilon\). In order to estimate the magnitude of any single harmonic, we devised the following procedure: the circumlunar spherical shell from 50 km to 2 lunar radii (\(\approx\)3500 km) is uniformly sampled at 100 values of altitude. Specifically, a sphere corresponding to each of these values is sampled at 9802 different points uniformly distributed in both latitude (100 values from \(-90^{\circ}\) to \(90^{\circ}\)) and longitude (100 values from \(-180^{\circ}\) to \(180^{\circ}\)). For a given harmonic, one can calculate the maximum perturbing acceleration it generates over 9802 points of equal altitude. The lowest degree (the coefficients of the lunar spherical harmonic model GRGM1200A used in this study were retrieved from [[https://pgda.gsfc.nasa.gov/products/50](https://pgda.gsfc.nasa.gov/products/50)]([https://pgda.gsfc.nasa.gov/products/50](https://pgda.gsfc.nasa.gov/products/50)), accessed on 4 November 2024) for which the spherical harmonics of any order produce the acceleration with a relative magnitude less than the threshold \(\varepsilon=10^{-7}\), is portrayed as a function of orbital altitude in Figure 1. The non-smooth curve is well approximated by the simple expression:  
\[N=\left\lceil\left(25/h\right)^{0.8}\right\rceil\]  
with altitude \(h\) measured in thousands km. The half square brackets denote the ceiling function. The approximation is accurate for low and medium lunar orbits up to 3000 km. The same sampling points of the circumlunar space were used to estimate the magnitude of gravitational perturbations due to the Earth and the Sun, as well as the solar radiation pressure effect. The corresponding dependencies are shown as functions of altitude in Figure 2. The cumulative perturbing effect due to lunar harmonics is found to dominate over all the other perturbations at lunar altitudes up to 2000 km. As will be explained later, two types of low lunar orbits are considered in this research: with an altitude of about 260 km and 520 km. The dynamical model adopted for the propagation of such orbits includes the \(N\times N\) lunar gravity field (\(N\) is selected according to the above given approximate formula: \(N=23\) for 260 km LLOs, whereas \(N=39\) for 520 km LLOs), the third-body gravitational perturbations due to the Earth and the Sun, and solar radiation pressure. The positions of celestial bodies were retrieved from JPL's DE430 ephemeris model [10]. The area-to-mass ratio of 0.02 m\({}^{2}\)/kg was assumed. To validate sufficiently high accuracy of the truncated dynamical model, we conduct a numerical experiment: two sets of 60 and 12 circular LLOs with an altitude of 260 km or, respectively, 520 km are numerically propagated for one year starting from 1 January 2022. Figure 1: Degree and order of the lunar gravitational potential expansion required to meet the truncated perturbing acceleration threshold \(\epsilon=10^{-7}\). Figure 2: Relative magnitude of perturbations as a function of the orbital altitude. All the orbits share the same inclination (\(84^{\circ}\)) and differ in only the right ascension of the ascending node (RAAN), whose values are sampled uniformly over the \([0,360]\) deg range.