# Exploiting Lunar Navigation Constellation for GNC Enhancement in Landing Missions
## 3 GNC Algorithms
### Navigation
\[\mathbf{f}_{S/C}(\mathbf{x}_{S/C},\bar{\mathbf{a}}_{IMU}) =\mathbf{f}_{S/C}(\mathbf{x}_{S/C})\ _{\ \mathbb{D}}+\mathbf{f}_{S/C}(\mathbf{x}_{S/C})_{\ \ \mathbb{D}}+\bar{\mathbf{a}}_{IMU} \tag{17}\] \[\mathbf{h}(\mathbf{x},\bar{\mathbf{x}}_{s,i}) =\left[\rho_{1},\ \rho_{1},\ \rho_{2},\ \rho_{2},\ \ldots,\ \rho_{n},\ \rho_{n},\ \zeta\right]^{\top} \tag{18}\]  
Equation (17) shows the strategy of dynamic model replacement to leverage the accelerometer within the navigation filter, encompassing both Solar Radiation Pressure (SRP) and the thrusters' acceleration contributions. It is also possible to see that only gravitational terms of Moon and Earth point mass are introduced. The Moon irregular mass distribution is then simply modelled by adding the oblateness, \(J_{2}\), contribution, in order to reduce the computational effort. The spacecraft state is estimated in the Lunar Mean Equator at J2000 (LME2000) frame, whose origin is the Moon's centre of mass, and its X-Y plane is defined by the Moon equator at J2000. In particular, the +Z axis is parallel to the Moon's rotation axis, pointing toward the north side of the invariant plane. The +X axis is defined by the intersection of the Moon's equator at the J2000 with the Earth Mean Equator of J2000, and the +Y axis completes the right-handed triad. Moreover, despite the fact that the broadcast ephemeris of traditional GNSS constellations are provided in the Earth-Centred-Earth-Fixed (ECEF) frame [32], in this scenario it is further assumed that the retrieved states of the LNS satellites are also directly expressed in the LME2000 frame, as a proper definition of an official Moon-Centred-Moon-Fixed (MCMF) frame is still not available. Additionally, it is worth underlining that the measurement function of Equation (18) integrates the estimated clock bias and drift to refine the pseudorange and range-rate a priori predictions, guided by Equations (9) and (10). The propagation model for these parameters is expressed by \(\mathbf{f}_{c}\), as defined in Equation (19). \[\mathbf{f}_{c}= \begin{cases}b_{c}&=d_{c}\\ d_{c}&=0\end{cases}  \tag{19}\]Combining Equations (17) and (19) results in the comprehensive dynamical model for the filter, denoted as \(\mathbf{f}=[\mathbf{f}_{S/C}^{\top},\ \mathbf{f}_{C}^{\top}]^{\top}\). Initialisation of the filter entails setting the initial conditions for the state, \(\hat{\mathbf{x}}_{0}\), and the state covariance matrix, \(\mathbf{P}_{0}\), at the beginning of the simulation. For \(\mathbf{P}_{0}\), we considered the values of \(\sigma_{r,0}=1\,\mathrm{km}\) and \(\sigma_{v,0}=1\,\mathrm{m}\,\mathrm{s}^{-1}\), which also serve as the standard deviations for \(\hat{\mathbf{x}}_{0}\). \(\mathbf{Q}\) and \(\mathbf{R}\) represent the process and measurement noise covariance matrices. If the former can be regarded as a tuning parameter, the latter is dictated by the accuracy of the sensors. The evaluation of the innovation vector defined by \(\mathbf{\hat{y}}_{k+1}-h(\hat{\mathbf{x}}_{k+1}^{-},\hat{\mathbf{x}}_{s,j,k+1})\) is influenced by two distinct errors: the pseudorange and range-rate errors, resulting from the physics of signal transmission and RF terminals, along with the errors arising in the ephemerides of the servicers, an outcome of their individual navigation budgets. Given the complex and more dynamically challenging scenario with respect to the natural motion studied in [14], the filter update frequency has been set to \(10\,\mathrm{Hz}\).