# Autonomous and Earth-Independent Orbit Determination for a Lunar Navigation Satellite System
## 4 Lunar Orbital Filter
### Orbital Propagation  
Orbital dynamics with third-body perturbations have been studied since the Apollo program [58], but considerations of a permanent orbiting service infrastructure are more recent. The frozen orbit, a mathematical outcome of three-body dynamics where gravity perturbations are minimized in certain orbital configurations, was proposed and explored more recently in the 2000s for the Moon by [26]. The ELFO is one such configuration, where LunaNet [4], Moonlight [5], and the Japanese LCNS [6] have all utilised this orbital 'anomaly' within their system architecture. The three-body problem is notorious for not being solvable algebraically. However, techniques have been developed to provide a numerical estimate. The Runge-Kutta four-step method is applied in this instance for calculating the propagation from time step \(k\), \(\mathbf{x}_{k}\), to time step \(k+1\), \(\mathbf{x}_{k+1}\)[59], which is reasonable for accuracies targeted in this work. However, a linear approximation is required for most filter-based solutions,  
\[\mathbf{x}_{k}=\mathbf{\phi}(\mathbf{x}_{k-1})\approx\Phi\mathbf{x}_{k}, \tag{23}\]  
where \(\mathbf{f}\) is a vector-valued function to calculate the propagation, and \(\Phi\) is the linearised matrix form of \(\mathbf{\phi}\), the propagation matrix. In this work, the state vector is composed of \(\mathbf{x}=\left\{\mathbf{r}^{\tau}\quad\mathbf{v}^{\tau}\right\}^{\tau}\). The acceleration of the spacecraft is caused mostly by gravitation forces from the Earth, Sun, and Moon, as well as solar-derived forces, such as radiation pressure and wind. Typically, for a small spacecraft orbiting the Moon, the resulting acceleration by the solar radiation pressure is in the order of 1 \(\times\) 10\({}^{-8}\) m/s\({}^{2}\)[60]. Solar wind is even less, causing accelerations at 0.5 AU from the Sun to 1 \(\times\) 10\({}^{-11}\) m/s\({}^{2}\)[61]. At this stage, solar-derived forces are not treated, as this work focuses on a covariance analysis to predict expected levels of accuracy with a proposed set of measurements, and not on dynamic model performance. Both the gravitational field of the Earth and Moon cannot be treated as point masses. Spherical harmonic models have been developed to a high level of fidelity, where the Earth Gravitational Model (EGM) provides a 2160-degree harmonic for the Earth [62] and the GRGM900C provides a 900-degree harmonic for the Moon [63]. To reduce processing time, with an assumed negligible loss in performance (on the order of 1 m per day [28]), the EGM and LP160C [64] models have been reduced to two-degree and four-degree harmonics for the Earth and Moon, respectively. These models are taken from Montenbruck for Earth-based orbital GNSS [65] and applied to the lunar domain. A point mass assumption is established for the Sun, using the Newtonian gravitational acceleration equation,  
\[\mathbf{g}_{s}=\frac{-Gm_{s}\mathbf{r}_{s}}{|\mathbf{r}_{s}|^{3}}, \tag{24}\]  
where \(G\) is the gravitational constant, \(m_{s}\) is the mass of the Sun, and \(\mathbf{r}_{s}\) is the vector pointing from the spacecraft to the Sun. The relation between \(\mathbf{r}\) and \(\mathbf{r}_{s}\) is \(\mathbf{r}_{s}=\mathbf{p}_{s}-\mathbf{r}\). The gravitational accelerations of the Earth, Moon, and Sun, \(g_{e}\), \(g_{m}\) and \(g_{s}\), respectively, can then be added together to calculate the acceleration. The time derivatives for the state vector \(\mathbf{x}\) are written as  
\[\dot{\mathbf{x}}=\left\{ \begin{matrix}\dot{\mathbf{r}}\\ \dot{\mathbf{v}}\end{matrix} \right\}=\left\{ \begin{matrix}\mathbf{v}\\ g_{e}+g_{m}+g_{s}\end{matrix} \right\}. \tag{25}\]  
Each is calculated as a function of the state \(\mathbf{x}\), \(\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})\). As mentioned, for filtering techniques, a linear solution to the non-linear propagation function \(\mathbf{\phi}\) must be first derived. This may be calculated through a linearised form of the dynamic model \(\mathbf{f}\), where the linearised form is the dynamic matrix \(F\).