# Low-Thrust Nonlinear Orbit Control for Very Low Lunar Orbits
## 2 Orbit Dynamics
Then, the gravitational acceleration can be computed in the LH-frame,  
\[\bar{\mathbf{g}}=\
abla U,\text{ where }\
abla=\hat{\mathbf{r}}\frac{\partial}{ \partial r}+\frac{\hat{\mathbf{E}}}{r\cos\phi}\frac{\partial}{\partial\lambda_{g }}+\frac{\hat{N}}{r}\frac{\partial}{\partial\phi}. \tag{5}\]  
Finally, by means of Equation (2), it is straightforward to obtain the components of the perturbing acceleration in the LVLH-frame. The strongest gravitational perturbation is associated with the coefficient \(J_{2}\). It is useful to express the perturbing acceleration due to the zonal harmonic \(J_{2}\) in the LVLH-frame,  
\[\bar{\mathbf{a}}_{J_{2}}=\frac{3\mu_{m}}{r^{4}}R_{m}^{2}J_{2}\Bigg{[} \begin{matrix} \frac{3\sin^{2}\left(\omega+\theta\right)\sin^{2}i-1}{2}\\ -\sin^{2}i\sin\left(\omega+\theta\right)\cos\left(\omega+\theta\right)\\ -\sin i\cos i\sin\left(\omega+\theta\right)\end{matrix} \Bigg{]}^{T}\!\!\!\! \begin{matrix}\hat{\mathbf{r}}\\ \hat{\mathbf{\theta}}\\ \hat{\mathbf{h}}\end{matrix} =\Big{[}a_{r}^{J_{2}}\quad a_{\theta}^{J_{2}}\quad a_ {h}^{J_{2}}\Bigg{[} \begin{matrix}\hat{\mathbf{r}}\\ \hat{\mathbf{\theta}}\\ \hat{\mathbf{h}}\end{matrix} \Bigg{]}, \tag{6}\]  
where \(\mu_{m}\) is the Moon gravitational parameter. The \(J_{2}\) harmonic is responsible for the time variation of the orbital elements, as predicted by the well-known Gauss equations [28],  
Figure 1: MCI and LVLH frames. \[\dot{p}=2\sqrt{\frac{p}{\mu_{m}}}ra_{\theta}^{l_{2}}, \tag{7}\] \[\dot{e}=\sqrt{\frac{p}{\mu_{m}}}a_{r}^{l_{2}}\sin\theta+\sqrt{\frac {p}{\mu_{m}}}a_{\theta}^{l_{2}}\frac{e+e\cos^{2}\theta+2\cos\theta}{1+e\cos \theta},\] (8) \[\dot{i}=ra_{\theta}^{l_{2}}\frac{\cos\left(\omega+\theta\right)}{ h},\] (9) \[\dot{\Omega}=ra_{h}^{l_{2}}\frac{\sin\left(\omega+\theta\right)}{ h\sin i},\] (10) \[\dot{\omega}=-ra_{h}^{l_{2}}\frac{\sin\left(\omega+\theta\right) \cos i}{h\sin i}-\frac{a_{r}^{l_{2}}}{e}\cos\theta\sqrt{\frac{p}{\mu_{m}}}+a_ {\theta}^{l_{2}}\sqrt{\frac{p}{\mu_{m}}}\sin\theta\frac{e\cos\theta+2}{e(1+e \cos\theta)},\] (11) \[\dot{\theta}=\sqrt{\frac{\mu_{m}}{p^{3}}}(1+e\cos\theta)^{2}+ \frac{a_{r}^{l_{2}}}{e}\cos\theta\sqrt{\frac{p}{\mu_{m}}}-a_{\theta}^{l_{2}} \sqrt{\frac{p}{\mu_{m}}}\sin\theta\frac{e\cos\theta+2}{e(1+e\cos\theta)}, \tag{12}\]  
where \(p\) is the semilatus rectum and \(e\) the eccentricity. However, it is possible to apply averaging techniques in order to filter short period oscillations. This procedure yields the average time derivatives of the orbital elements, which equal zero for the semimajor axis, the eccentricy, and the inclination. Instead, the RAAN, the argument of pericenter and the true anomaly exhibit a secular variation. Therefore, the proposed strategy is to employ the control effort to maintain a constant semimajor axis, eccentricity and inclination, while complying with the variation of the RAAN, argument of pericenter and true anomaly, in order to achieve propellant saving. The orbit dynamics can be described by means of modified equinoctial elements \(x_{i}\), \(i=1,\ldots,6\),  
\[x_{1}=a(1-e^{2}), \tag{13}\] \[x_{2}=e\cos\left(\Omega+\omega\right),\] (14) \[x_{3}=e\sin\left(\Omega+\omega\right),\] (15) \[\dot{p}=2\sqrt{\frac{p}{\mu_{m}}}ra_{\theta}^{l_{2}}, \tag{7}\]0 (16) \[\dot{p}=2\sqrt{\frac{p}{\mu_{m}}}ra_{\theta}^{l_{2}}, \tag{7}\]1 (17) \[\dot{p}=2\sqrt{\frac{p}{\mu_{m}}}ra_{\theta}^{l_{2}}, \tag{7}\]2  
where \(a\) is the semi-major axis. This set of elements is well-suited for numerical integration, as it does not exhibit singularities for every keplerian trajectories, with the exception of equatorial retrograde orbits. Given the vector of non Keplerian accelerations \(\mathbf{a}= \begin{bmatrix}a_{r}&a_{\theta}&a_{h}\end{bmatrix} ^{T}\) expressed in the LVLH frame and \(\mathbf{x}\coloneqq \begin{bmatrix}x_{1}&x_{2}&x_{3}&x_{4}&x_{5}\end{bmatrix} ^{T}\), the dynamical equations can be written as  
\[\dot{\mathbf{x}}=\mathbf{G}\left(\mathbf{x},x_{6}\right)a, \tag{19}\]where \(\eta=1+x_{2}\cos x_{6}+x_{3}\sin x_{6}\) and  
\[\mathbf{G}=\sqrt{\frac{x_{1}}{\mu_{m}}}