Quantitative and Qualitative Results from Orbital Plane Precession of Medium-high Near-circular Orbits

According to the transformation relation between the geocentric equatorial inertial coordination system O − XYZ and the coordinate system O − 123 composed by the principal axes in Figure 1 and the theoretical analysis in Allan & Cook (1964), Zhao et al. (2015) merely established the linkage between the orbital elements and the mechanism of classical Laplace precession of the orbital plane by using Equation (3) in Zhao et al. (2015) and Equation (26a) in Allan & Cook (1964) and further derived the analytical formulae of the quantitative variation ranges of orbital inclination for near-circular geosynchronous orbits in theory. Equation (3) in Zhao et al. (2015) and Equation (26a) in Allan & Cook (1964), mentioned here, and Equations (22), (26b), (27a), and (27b) in Allan & Cook (1964), used below, are presented in the Appendix for reference.

However, as stated in Section 2, there are two types of possible mechanisms of orbital plane precession for medium-high near-circular orbits (not just the classical Laplace precession); hence, the analytical formulae of the quantitative variation ranges of orbital inclination obtained in Zhao et al. (2015) are theoretically incomplete.

For this reason, we examine these two types of possible precession mechanisms (the classical and polar Laplace precessions) of the orbital plane for medium-high near-circular orbits as a whole in this section, and we revise the analytical formulae of the quantitative variation ranges of orbital inclination derived in Zhao et al. (2015) by using a similar methodology. First, using Equation (3) in Zhao et al. (2015) and Equation (22) in Allan & Cook (1964), we establish the linkage between the orbital elements and the nontrivial integral λ₀ of the long-term evolution of the orbital plane, from which the long-term evolution of the orbital plane governed by the classical Laplace precession mechanism or the polar Laplace precession mechanism, depending on the initial orbital values, can be determined directly. Then, we introduce two parameters, ${i}_{\ast }={\left.i\right|}_{{\rm{\Omega }}={0}^{\circ }}$ and ${i}_{\ast }^{{\rm{{\prime} }}}={\left.i\right|}_{{\rm{\Omega }}={90}^{\circ }}$ (i.e., i_* and ${i}_{* }^{{\prime} }$ are the orbital inclinations when the longitude of the ascending node Ω = 0° and 90°, respectively), and present their specific expressions by using Equations (26a) and (26b) in Allan & Cook (1964). Ultimately, based on the classical and polar Laplace precession mechanisms of the orbital plane, we derive the complete analytical formulae of the quantitative variation ranges of orbital inclination for medium-high near-circular orbits as follows:

$$\begin{eqnarray}{\rm{\Delta }}i=\left\{\begin{array}{ll}\left[{i}_{* }-2\alpha ,{i}_{* }\right] & \left({i}_{* }\geqslant 2\alpha ,{\lambda }_{2}\lt {\lambda }_{0}\leqslant {\lambda }_{3}\right),\\ \left[2\alpha -{i}_{* },{i}_{* }\right] & \left(\alpha \leqslant {i}_{* }\leqslant 2\alpha ,{\lambda }_{2}\lt {\lambda }_{0}\leqslant {\lambda }_{3}\right),\\ \left[{i}_{* },2\alpha -{i}_{* }\right] & \left({i}_{* }\leqslant \alpha ,{\lambda }_{2}\lt {\lambda }_{0}\leqslant {\lambda }_{3}\right),\\ \left[180^\circ -{i}_{* }^{{\prime} },{i}_{* }^{{\prime} }\right] & \left({i}_{* }^{{\prime} }\geqslant 90^\circ ,{\lambda }_{1}\leqslant {\lambda }_{0}\lt {\lambda }_{2}\right),\\ \left[{i}_{* }^{{\prime} },180^\circ -{i}_{* }^{{\prime} }\right] & \left({i}_{* }^{{\prime} }\leqslant 90^\circ ,{\lambda }_{1}\leqslant {\lambda }_{0}\lt {\lambda }_{2}\right),\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{c}\begin{array}{c}\left\{\begin{array}{c}{\lambda }_{0}={\lambda }_{2}{\left(\sin \alpha \cos {i}_{0}-\cos \alpha \cos {{\rm{\Omega }}}_{0}\sin {i}_{0}\right)}^{2}+{\lambda }_{3}{\left(\cos \alpha \cos {i}_{0}+\sin \alpha \cos {{\rm{\Omega }}}_{0}\sin {i}_{0}\right)}^{2},\\ {\sin }^{2}\left({i}_{\ast }-\alpha \right)={\sin }^{2}\left({\left.i\right|}_{{\rm{\Omega }}={0}^{\circ }}-\alpha \right)=\displaystyle \frac{{\lambda }_{3}-{\lambda }_{0}}{{\lambda }_{3}-{\lambda }_{2}}\,\,\,\,\left({\lambda }_{2}\lt {\lambda }_{0}\leqslant {\lambda }_{3}\right),\\ {\cos }^{2}{i}_{\ast }^{{\rm{{\prime} }}}={\cos }^{2}{\left.i\right|}_{{\rm{\Omega }}={90}^{\circ }}=\displaystyle \frac{{\lambda }_{0}}{{\lambda }_{2}{\sin }^{2}\alpha +{\lambda }_{3}{\cos }^{2}\alpha }\,\,\,\,\left({\lambda }_{1}\leqslant {\lambda }_{0}\lt {\lambda }_{2}\right).\end{array}\right.\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 2 }$$

It should be specifically pointed out that the above complete analytical formulae of the quantitative variation ranges of orbital inclination for medium-high orbits are derived on the precondition of near-circular orbit. In other words, Equations (1) and (2) are valid only if the initial near-circular orbital characteristic is maintained in the long-term evolution of medium-high orbits. The above amendments of the analytical formulae of the quantitative variation ranges of orbital inclination in Zhao et al. (2015) not only have particular regard for the polar Laplace precession mechanism of the orbital plane but also present adaptable replacements of the original expressions, Equations (2) and (4) in Zhao et al. (2015), of the quantities i_* and λ₀.

For near-circular geosynchronous orbits with a nominal altitude H = 35,786 km, initial inclination i₀ = 55°, and different initial longitude of the ascending node Ω₀ in the BeiDou navigation satellite system (China Satellite Navigation Office 2013), it is seen that the orbital plane precession mechanism is the classical Laplace precession from Figure 2. Likewise, the orbital plane precession mechanism is also the classical Laplace precession for the near-circular medium orbits with a nominal altitude H = 21,528 km, initial inclination i₀ = 55°, and different initial longitude of the ascending node Ω₀ in the BeiDou navigation satellite system (China Satellite Navigation Office 2013). From Equations (1) and (2), we obtain a series of specific theoretical results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ in the long-term evolution of these orbits for different initial longitudes of the ascending node Ω₀ at intervals of 1° in the range of $\left[0^\circ ,360^\circ \right]$, upon the precondition that the initial near-circular orbital characteristic is maintained; see the blue symbols in Figures 3(a) and (b). It is observed that the quantitative variation amplitudes of orbital inclination i_max − i_min in the long-term orbital evolution are all 2α for given initial inclination i₀ = 55° and different initial longitude of the ascending node Ω₀; when Ω₀ = 0° and 360°, i_min and i_max have the minimum values, while when Ω₀ = 180°, i_min and i_max have the maximum values from Figures 3(a) and (b). These are the marked features of the classical Laplace precession of the orbital plane for medium-high near-circular orbits revealed through the analytical formulae of Equations (1) and (2).

Zoom In Zoom Out Reset image size

However, from Figure 2, it is seen that the classical and polar Laplace precession mechanisms both exist in the long-term evolution of the orbital plane for near-circular geosynchronous orbits with an altitude H = 35,786 km and initial inclination i₀ = 83° when the initial longitude of the ascending node Ω₀ takes different values. The same phenomenon can also be observed for near-circular medium orbits with an altitude H = 21,528 km and initial inclination i₀ = 87° when the initial longitude of the ascending node Ω₀ takes different values. Likewise, from Equations (1) and (2), we obtain a series of specific theoretical results of quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ in the long-term evolution of these orbits for different initial longitudes of the ascending node Ω₀ at intervals of 1° in the range of $\left[0^\circ ,360^\circ \right]$, upon the precondition that the initial near-circular orbital characteristic is maintained; see the blue and red symbols in Figures 4(a) and 5(a), respectively. The blue symbols correspond to the theoretical results from the classical Laplace precession of the orbital plane, while the red symbols correspond to the theoretical results from the polar Laplace precession of the orbital plane. It is observed that the quantitative variation amplitudes of orbital inclination i_max − i_min determined by the polar Laplace precession of the orbital plane are not a fixed value any longer as that determined by the classical Laplace precession of the orbital plane. It is a marked feature of the polar Laplace precession of the orbital plane for medium-high near-circular orbits revealed through the analytical formulae of Equations (1) and (2).

Zoom In Zoom Out Reset image size

Using the expression of λ₀ in Equation (2) above and Equations (27a) and (27b) in Allan & Cook (1964), the periods of orbital plane evolution (or the variation periods of orbital inclination) for medium-high near-circular orbits can also be calculated. For medium-high near-circular orbits with initial values as set in Figures 3(a) and (b), 4(a), and 5(a), we present a series of specific theoretical results of the periods of orbital plane evolution P upon the precondition that the initial near-circular orbital characteristic is maintained; see Figures 3(c) and (d), 4(b), and 5(b), respectively. The blue symbols in the figures correspond to the theoretical results from the classical Laplace precession of the orbital plane, while the red symbols correspond to the polar Laplace precession.

From Figures 3(c) and (d), it is seen that for the medium-high near-circular orbits with the given initial inclination i₀ and different initial longitude of the ascending node Ω₀, the period of orbital plane evolution P reaches the minimum when Ω₀ = 0° and 360° and the maximum when Ω₀ = 180°. It is another marked feature determined by the classical Laplace precession of the orbital plane for medium-high near-circular orbits.

To verify the above theoretical calculation results of quantitative variation ranges of orbital inclination based on the simplified and approximate model and to obtain information on the long-term evolution of orbital eccentricity, the first-order semianalytical calculations with simplified models are performed. Meanwhile, parallel calculations based on the simplified models with and without the practical motion of the lunar pole are conducted to reveal its influence on the long-term evolution of orbital inclination and eccentricity.

The second-order zonal harmonic J₂ of the Earth's gravitational field and the lunisolar perturbation terms up to the second-order Legendre polynomial ${P}_{2}(\cos \psi )$ are still taken into account in the simplified models. The first-order solutions take the complete first-order perturbation terms and the second-order secular terms of the perturbation factors of the simplified models. Thus, the first- and second-order secular terms, the first-order long-period terms, and the first-order short-period terms are involved in the first-order solutions (Liu 2000; Liu & Tang 2015).

In the first-order semianalytical calculations, the nonsingular elements applied to arbitrary eccentricity (0 ≤ e < 1), ${\boldsymbol{\sigma }}\,={\left(a,i,{\rm{\Omega }},\xi =e\sin \omega ,\eta =e\cos \omega ,\lambda =M+\omega \right)}^{T}$, are adopted as the basic variables, where a, e, i, Ω, ω, and M are the classical Keplerian elements. The considered perturbation models are averaged over the fast variable λ or M up to the terms of ${J}_{2}^{2}$, and then the averaged perturbation models are numerically integrated by the Runge–Kutta–Fehlberg 7(8) method to acquire the first- and second-order secular terms and the first-order long-period terms. Meantime, the first-order short-period terms at the corresponding times are calculated by using the analytical formulae in Liu et al. (2005), and then the first-order solutions are ultimately obtained. The minimum and maximum inclinations i_min and i_max and the maximum eccentricity e_max in the long-term orbital evolution are directly obtained through successive comparisons in the process of the first-order semianalytical calculations. In addition, the JPL DE430 precise ephemeris³ is directly invoked to get the required lunisolar position coordinates in the first-order semianalytical calculations.

For the medium-high near-circular orbits with altitudes $H=35,786\,$ and $21,528\,{\rm{k}}{\rm{m}}$, given initial inclinations i₀, and different initial longitudes of the ascending node Ω₀ discussed in Section 3.1, the first-order semianalytical calculation results of the quantitative variation ranges of inclination $\left[{i}_{\min },{i}_{\max }\right]$ and maximum eccentricity e_max, based on the simplified models with and without the practical motion of the lunar pole in the long-term orbital evolution, are presented in Figures 3(a), (b), (e), and (f); 4(c) and (d); and 5(c) and (d), respectively, by using different colored symbols. The cyan symbols in the figures correspond to the results without the practical motion of the lunar pole, while the magenta symbols correspond to the results with the practical motion of the lunar pole. The corresponding initial orbital values in these first-order semianalytical calculations are set as follows: the initial epoch is 00:00:00 UTC on 2019 May 5, the initial semimajor axis a₀ = 42,164.1363 or 27,906.1363 km, the initial eccentricity e₀ = 0.0001, the initial longitude of the ascending node Ω₀ is set at intervals of 1° in the range of $\left[0^\circ ,360^\circ \right]$, the initial argument of perigee ω₀ = 0°, and the initial mean anomaly M₀ = 0°. To make the time spans of the first-order semianalytical calculations cover a variation period of orbital inclination, we take the theoretical values of the periods of orbital plane evolution P as references and set the time spans of the first-order semianalytical calculations as 120 yr for Figures 3(a) and (e), 40 yr for Figures 3(b) and (f), and 500 yr for Figures 5(c) and (d). Further considering the time span that the JPL DE430 precise ephemeris covers, we set the time spans of the first-order semianalytical calculations in Figures 4(c) and (d) as 600 yr. In addition, if the orbital perigee altitude is below 80 km in the first-order semianalytical calculation, we assume that orbital decay occurs, and then the calculation will be aborted.

Zoom In Zoom Out Reset image size

From Figures 3(a) and (e), it can be seen that for certain initial longitudes of the ascending node Ω₀, the maximum eccentricities e_max in the long-term evolution up to 120 yr for the initial near-circular geosynchronous orbits with a nominal altitude H = 35,786 km and initial inclination i₀ = 55° in the BeiDou navigation satellite system can even exceed 0.5 due to lunisolar resonances (Cook 1962; Hughes 1980, 1981; Zhao et al. 2015), but the theoretical results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ accord well with the results of the first-order semianalytical calculations. Meanwhile, the first-order semianalytical calculation results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ and the maximum eccentricities e_max with and without the practical motion of the lunar pole are consistent on the whole.

From Figures 3(b) and (f), it can be seen that the long-term evolution of eccentricities up to 40 yr for the initial near-circular medium orbits with a nominal altitude $H={\rm{21,528}}\,\mathrm{km}$ and initial inclination i₀ = 55° in the BeiDou navigation satellite system can maintain their initial near-circular characteristic; the theoretical results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ accord well with the results of the first-order semianalytical calculations without the practical motion of the lunar pole but roughly with the results with the practical motion of the lunar pole, and there are relatively obvious differences. Meanwhile, the first-order semianalytical calculation results of the maximum eccentricities e_max with and without the practical motion of the lunar pole are consistent on the whole.

From Figures 4(c) and (d), it can be seen that all of the maximum eccentricities e_max in the long-term evolution up to 600 yr for the initial near-circular orbits with a geosynchronous altitude H = 35,786 km, initial inclination i₀ = 83°, and different initial longitude of the ascending node Ω₀ can exceed 0.6 due to lunisolar resonances (Cook 1962; Hughes 1980, 1981; Zhao et al. 2015), and the theoretical results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$, upon the precondition that the initial near-circular orbital characteristic is maintained, have obvious differences from the first-order semianalytical calculations. However, to a certain extent, the first-order semianalytical calculation results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ reflect that both the classical and polar Laplace precession mechanisms of the orbital plane exist for the initial near-circular orbits with a geosynchronous altitude H = 35,786 km and initial inclination i₀ = 83°. The first-order semianalytical calculation results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ and the maximum eccentricities e_max with and without the practical motion of the lunar pole are also consistent on the whole.

From Figures 5(c) and (d), it can be seen that the long-term evolution of eccentricities up to 500 yr for the initial near-circular medium orbits with an altitude H = 21,528 km, initial inclination i₀ = 87°, and different initial longitude of the ascending node Ω₀ can maintain their initial near-circular characteristic. The theoretical results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ accord well with the results of the first-order semianalytical calculations with and without the practical motion of the lunar pole on the whole, but there are obvious differences, in particular in the case of the classical Laplace precession of the orbital plane. The first-order semianalytical calculation results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ with and without the practical motion of the lunar pole are consistent on the whole, and there are few differences for the maximum eccentricities e_max.

To better compare the theoretical results and the first-order semianalytical calculation results of the quantitative variation ranges of orbital inclination $\left[{i}_{\min },{i}_{\max }\right]$ for the long-term evolution of initial near-circular medium orbits with an altitude H = 21,528 km, initial inclination i₀ = 87°, and different initial longitude of the ascending node Ω₀, we replot these results in Figure 6. From Figure 6, it can be seen that in the case of the polar Laplace precession of the orbital plane, the theoretical results are well consistent with the first-order semianalytical calculation results. Then, it is verified that the polar Laplace precession mechanism of the orbital plane does exist through the first-order semianalytical calculations. The obvious differences in the case of the classical Laplace precession of the orbital plane lie in the inconsistent results of i_min derived by the theoretical analysis and the first-order semianalytical calculations.

Zoom In Zoom Out Reset image size

We further present the specific results of the long-term evolution of inclination and eccentricity for the initial near-circular orbits with altitudes H = 35,786 and 21,528 km, initial inclination i₀ = 55°, and different initial longitude of the ascending node Ω₀ = 0°, 60°, 120°, 180°, 240°, and 300° in Figures 7 and 8, respectively, by the first-order semianalytical calculations based on the simplified models with and without the practical motion of the lunar pole. Similarly, when the initial longitude of the ascending node Ω₀ = 0°, 60°, 120°, 180°, 240°, and 300°, the specific results for the initial near-circular orbits with an altitude H = 35,786 km and initial inclination i₀ = 83° and those with an altitude H = 21,528 km and initial inclination i₀ = 87° are presented in Figures 9 and 10, respectively. The cyan curves in Figures 7–10 correspond to the results without the practical motion of the lunar pole, while the magenta curves correspond to the results with the practical motion of the lunar pole. The influence of the practical motion of the lunar pole on the long-term evolution of orbital inclination and eccentricity can be further observed from Figures 7–10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The periods of orbital plane evolution observed from the long-term evolution of orbital inclination in Figures 7, 8 and 10 are in accord with the corresponding theoretical results in Figures 3(c) and (d) and 5(b) on the whole. This verifies the theoretical calculation of the period of orbital plane evolution by the first-order semianalytical calculations based on the simplified models to a certain extent.

In addition, it seems that the curves of the long-term evolution of inclination and eccentricity are broken for the cases Ω₀ = 0°, 60°, 180°, and 300° in Figure 9 derived by the first-order semianalytical calculations based on the simplified models without the practical motion of the lunar pole, as well as for the case Ω₀ = 0° in Figure 9 derived by the first-order semianalytical calculation based on the simplified models with the practical motion of the lunar pole. The specific reason for this phenomenon is as follows. Due to lunisolar resonances (Cook 1962; Hughes 1980, 1981; Zhao et al. 2015), the orbital eccentricities e can evolve to high enough, and the related orbital perigee altitudes can evolve to low enough (below 80 km). As stated above, when the orbital perigee altitude is below 80 km in the first-order semianalytical calculations, orbital decay is supposed to occur, and then the calculations are aborted, which makes these curves seem to be broken.

It should also be pointed out that there are model error accumulations in the first-order semianalytical calculations through numerical integrations of the averaged perturbation models, which cause the differences between the results with and without the practical motion of the lunar pole. Especially, from Figure 9, it is seen that there is a good coincidence between the curves of the long-term evolution of inclination and eccentricity, derived by the first-order semianalytical calculations with and without the practical motion of the lunar pole, within the first several decades from the initial epoch in the long-term calculations. Then, with the increase of the model error accumulations, the relatively obvious differences can be observed. This is the reason that orbital decays occur for the cases Ω₀ = 60°, 180°, and 300° in Figure 9 in the first-order semianalytical calculations without the practical motion of the lunar pole, while orbital decays do not occur for these cases in the first-order semianalytical calculations with the practical motion of the lunar pole.

To compare the first-order semianalytical calculation results based on the simplified models and further verify the theoretical results based on the simplified and approximate model, we perform numerical calculations based on the precise models for the long-term evolution of initial near-circular orbits with altitudes H = 35,786 and 21,528 km; initial inclinations i₀ = 55°, 83°, and 87°; and different initial longitude of the ascending node Ω₀ = 0°, 60°, 120°, 180°, 240°, and 300°, respectively, as in Section 3.3. The corresponding long-term evolution of inclinations and eccentricities for these orbits obtained by numerical calculations is presented in Figures 7–10 with black curves.

The Earth's gravitational field up to 20 × 20 order and degree, lunisolar perturbations, solar radiation pressure, and atmospheric drag perturbation (when the orbital altitude is less than 1000 km) are taken into account in the precise models. The JGM3 model is used to get the coefficients of the Earth's gravitational field. The JPL DE430 precise ephemeris is directly invoked to get the required lunisolar position coordinates, and the U.S. Standard Atmosphere 1976 model is used to get the atmospheric mass densities in the numerical calculations. The solar radiation pressure coefficient is set as 1.2, and the A/m value of the object in orbits is set as $0.01\,{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}$. The initial orbital values of the numerical calculations are the same as the first-order semianalytical calculations. Likewise, if the orbital perigee altitude is below 80 km in the numerical calculation, we assume that orbital decay occurs, and then the calculation will be aborted.

From Figures 7–10, it is shown that although the long-term evolution of eccentricities obtained by numerical calculations with the precise models has relatively obvious differences from the first-order semianalytical calculations with the simplified models, the dynamical tendencies whether the initial near-circular characteristic can be maintained in the long-term evolution, revealed from the first-order semianalytical calculations and numerical calculations, are consistent. Meanwhile, when the initial near-circular orbital characteristic is maintained, most of the long-term evolution of inclinations obtained by numerical calculations accord with the first-order semianalytical calculations on the whole. This suggests the reliability of the first-order semianalytical calculations and further verifies the theoretical results based on the simplified and approximate model. In other words, the complete analytical formulae of the quantitative variation ranges of orbital inclination in Equations (1) and (2), upon the precondition that the initial near-circular orbital characteristic is maintained, are verified to be correct.

The orbital decay also occurs in the numerical calculation for the case Ω₀ = 0° in Figure 9, mainly due to lunisolar resonances (Cook 1962; Hughes 1980, 1981; Zhao et al. 2015), and then the numerical calculation is aborted, which makes the curves of the long-term evolution of inclination and eccentricity seem to be broken. Meanwhile, the atmospheric drag perturbation is taken into account in the numerical calculation when the orbital altitude is below 1000 km, which accelerates the orbital decay. The differences between the results derived by the first-order semianalytical and numerical calculations mainly arise from the model errors and their accumulations.

Similar to the qualitative characteristic formula for the variation of orbital inclination-longitude of the ascending node of geostationary orbits in Zhao et al. (2016), we use Equation (3) in Zhao et al. (2015), Equations (26a) and (26b) in Allan & Cook (1964), and the above parameters ${i}_{\ast }={\left.i\right|}_{{\rm{\Omega }}={0}^{\circ }}$ and ${i}_{\ast }^{{\rm{{\prime} }}}={\left.i\right|}_{{\rm{\Omega }}={90}^{\circ }}$ defined in Equation (2) and derive the generalized qualitative characteristic formulae of the variation of orbital inclination-longitude of the ascending node for medium-high near-circular orbits

$$\begin{eqnarray}&&\begin{array}{r}\left\{\begin{array}{l}\cos {\rm{\Omega }}={\left[\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\left\{\pm \left[\displaystyle \frac{{\cos }^{2}\left({i}_{\ast }-\alpha \right)}{{\sin }^{2}\alpha {\sin }^{2}i}+\displaystyle \frac{{\lambda }_{2}}{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha }\right.\right.\\ {\left.-\displaystyle \frac{{\lambda }_{2}{\cot }^{2}\alpha {\cot }^{2}i}{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\left.\,-\,{\left[\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\cot \alpha \cot i\right\}\\ =\,{\left[\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\left\{\pm {\left[\displaystyle \frac{{\lambda }_{0}-{\lambda }_{2}{\cos }^{2}i}{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha {\sin }^{2}i}-\displaystyle \frac{{\lambda }_{2}{\cot }^{2}\alpha {\cot }^{2}i}{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\right.\\ \left.-\,{\left[\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\cot \alpha \cot i\right\}\,\,({\lambda }_{2}\lt {\lambda }_{0}\leqslant {\lambda }_{3}),\end{array}\right.\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{c}\left\{\begin{array}{l}\cos {\rm{\Omega }}=\pm {\left[\displaystyle \frac{\left({\lambda }_{2}{\tan }^{2}\alpha +{\lambda }_{3}\right){\cos }^{2}i{{\rm{{\prime} }}}_{\ast }}{{\lambda }_{3}{\cos }^{2}i}-\displaystyle \frac{{\lambda }_{2}{\left(1+{\tan }^{2}\alpha \right)}^{2}}{{\lambda }_{3}{\tan }^{2}\alpha +{\lambda }_{2}}\right]}^{1/2}{\left[{\left({\tan }^{2}\alpha +\displaystyle \frac{{\lambda }_{2}}{{\lambda }_{3}}\right)}^{1/2}\tan i\right]}^{-1}-\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right)\tan \alpha }{\left({\lambda }_{3}{\tan }^{2}\alpha +{\lambda }_{2}\right)\tan i}\\ =\,\pm {\left[\displaystyle \frac{{\lambda }_{0}}{{\lambda }_{3}{\cos }^{2}\alpha {\cos }^{2}i}-\displaystyle \frac{{\lambda }_{2}{\left(1+{\tan }^{2}\alpha \right)}^{2}}{{\lambda }_{3}{\tan }^{2}\alpha +{\lambda }_{2}}\right]}^{1/2}{\left[{\left({\tan }^{2}\alpha +\displaystyle \frac{{\lambda }_{2}}{{\lambda }_{3}}\right)}^{1/2}\tan i\right]}^{-1}-\displaystyle \frac{\left({\lambda }_{3}-{\lambda }_{2}\right)\tan \alpha }{\left({\lambda }_{3}{\tan }^{2}\alpha +{\lambda }_{2}\right)\tan i}\,\,({\lambda }_{1}\leqslant {\lambda }_{0}\lt {\lambda }_{2}).\end{array}\right.\end{array}\end{eqnarray} \tag{ 4 }$$

Equation (3) applies to the classical Laplace precession of the orbital plane. For the mathematical symbol "±" in Equation (3), the plus sign corresponds to the case ${\left[\tfrac{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }{({\lambda }_{3}-{\lambda }_{2}){\sin }^{2}\alpha }\right]}^{1/2}\cos {\rm{\Omega }}\,+\,{\left[\tfrac{({\lambda }_{3}-{\lambda }_{2}){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\cot \alpha \cot i\geqslant 0$, while the minus sign corresponds to the case ${\left[\tfrac{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }{({\lambda }_{3}-{\lambda }_{2}){\sin }^{2}\alpha }\right]}^{1/2}\cos {\rm{\Omega }}\,+\,{\left[\tfrac{({\lambda }_{3}-{\lambda }_{2}){\sin }^{2}\alpha }{{\lambda }_{3}{\sin }^{2}\alpha +{\lambda }_{2}{\cos }^{2}\alpha }\right]}^{1/2}\cot \alpha \cot i\lt 0$. Equation (4) applies to the polar Laplace precession of the orbital plane. For the mathematical symbol "±" in Equation (4), the plus sign corresponds to the case $\tfrac{\left({\lambda }_{3}-{\lambda }_{2}\right)\tan \alpha }{{\left({\lambda }_{3}^{2}{\tan }^{2}\alpha +{\lambda }_{2}{\lambda }_{3}\right)}^{1/2}}\,+{\left({\tan }^{2}\alpha +\tfrac{{\lambda }_{2}}{{\lambda }_{3}}\right)}^{1/2}\tan i\cos {\rm{\Omega }}\geqslant 0$, while the minus sign corresponds to the case $\tfrac{\left({\lambda }_{3}-{\lambda }_{2}\right)\tan \alpha }{{\left({\lambda }_{3}^{2}{\tan }^{2}\alpha +{\lambda }_{2}{\lambda }_{3}\right)}^{1/2}}+{\left({\tan }^{2}\alpha +\tfrac{{\lambda }_{2}}{{\lambda }_{3}}\right)}^{1/2}\tan i\cos {\rm{\Omega }}\lt 0$. As stated above, the quantities λ₀, λ₁, λ₂, λ₃, and α are constants determined by the initial state or values of the given near-circular orbit, if the assumption that the lunar orbital plane coincides with the ecliptic is considered as in Allan & Cook (1964). Then, from the initial orbital values, the qualitative characteristic of the variation of orbital inclination-longitude of the ascending node for an arbitrary medium-high near-circular orbit can be presented by Equations (3) and (4). This shows the generality of the qualitative characteristic formulae of Equations (3) and (4).

The qualitative characteristics of the variation of orbital inclination-longitude of the ascending node, derived by the theoretical formulae of Equations (3) and (4) for these initial near-circular orbits discussed in Figures 7–10, are presented in Figures 11–14 with the blue and red curves, respectively. The blue curves correspond to the classical Laplace precession of the orbital plane, while the red curves correspond to the polar Laplace precession of the orbital plane. It is worth pointing out that the qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node governed by the above theoretical formulae are in pairs. In general cases, according to the initial orbital values, the specific and single variation curve of orbital inclination-longitude of the ascending node for the long-term evolution of the initial near-circular orbits can be determined from the paired curves. Using the initial orbital values, we directly present the specific and single qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node for the long-term evolution of the corresponding orbits in Figures 11, 12 and 14, as well as the cases Ω₀ = 120° and 240° in Figure 13 for the polar Laplace precession of the orbital plane. For the cases Ω₀ = 0°, 60°, 180°, and 300° in Figure 13 in the state of classical Laplace precession of the orbital plane, we present the paired qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node. Moreover, the paired qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node for the case Ω₀ = 180° are interwoven. When the paired qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node are close to each other or interwoven, predictably, some special dynamical phenomena may occur.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figures 11–14, it can be seen that the qualitative characteristics of the variation of orbital inclination-longitude of the ascending node in cases of classical and polar Laplace precessions of the orbital plane are markedly different. For the classical Laplace precession of the orbital plane, the long-term evolution of the orbital plane for these orbits is in a state of circulation; i.e., the variation range of the longitude of the ascending node is $\left[0^\circ ,360^\circ \right]$. For the polar Laplace precession of the orbital plane, the long-term evolution of the orbital plane for these orbits is in a state of libration; i.e., the longitude of the ascending node varies merely in a certain range, and the transition between prograde and retrograde motion for these near-circular orbits occurs. In addition, it should be noted that for geostationary orbits in the mechanism of the classical Laplace precession of the orbital plane, the variation range of the longitude of the ascending node is $\left[-90^\circ ,90^\circ \right]$ (Zhao et al. 2016), which is different from the near-circular orbits with high inclinations discussed in this paper. The features of the long-term evolution of the orbital plane in a state of libration and the transition between prograde and retrograde motion for the polar Laplace precession of the orbital plane might offer a reference for the orbit design to meet the particular mission requirement.

To verify the theoretical results of the qualitative characteristics of the variation of orbital inclination-longitude of the ascending node given by the analytical formulae of Equations (3) and (4), the results from the first-order semianalytical calculations based on the simplified models with and without the practical motion of the lunar pole and from the numerical calculations based on the precise models are also presented in Figures 11–14 with the cyan, magenta, and black curves respectively.

From Figures 11 and 12, it can be seen that the theoretical results of the qualitative characteristics of the variation of orbital inclination-longitude of the ascending node given by the analytical formula of Equation (3) from the classical Laplace precession of the orbital plane are in accord with the results of the first-order semianalytical and numerical calculations on the whole.

Although the initial near-circular orbital characteristic is no longer maintained in the long-term evolution as seen in Figure 9, from Figure 13, the theoretical results of the qualitative characteristics of the variation of orbital inclination-longitude of the ascending node for these orbits with the initial longitude of the ascending node Ω₀ = 120° and 240°, given by the analytical formula of Equation (4) from the polar Laplace precession of the orbital plane, are in accord with the results of the first-order semianalytical and numerical calculations on the whole. For the cases Ω₀ = 0°, 60°, 180°, and 300° in Figure 13 in the state of the classical Laplace precession of the orbital plane, the overall trend of the one certain characteristic curve for these orbits given by the analytical formula of Equation (3) is consistent with the results of the first-order semianalytical and numerical calculations in the initial stage of long-term evolution. The results of the first-order semianalytical calculations are also in accord with the results of numerical calculations in the initial stage of long-term evolution on the whole, and then the obvious differences are gradually revealed with time. In particular, for the case Ω₀ = 60° in Figure 13, the variation of orbital inclination-longitude of the ascending node in the long-term evolution given by numerical calculation presents the transition between the paired qualitative characteristic curves determined by the analytical formula of Equation (3). It shows that special dynamical phenomena might occur when the paired qualitative characteristic curves of the variation of orbital inclination-longitude of the ascending node determined by the theoretial formula are close to each other.

From Figure 14, it can be seen that the theoretical results for these orbits with the initial longitude of the ascending node Ω₀ = 60°, 120°, 240°, and 300° from the polar Laplace precession of the orbital plane are also consistent with the results of the first-order semianalytical and numerical calculations on the whole. The theoretical results for these orbits with the initial longitude of the ascending node Ω₀ = 0° and 180° from the classical Laplace precession of the orbital plane are well consistent with the results of the first-order semianalytical calculations. Moreover, in the initial stage of long-term evolution of the orbital plane, the theoretical results for the cases Ω₀ = 0° and 180° are well consistent with the results of the numerical calculations, and they are also well consistent between the results of the first-order semianalytical and numerical calculations; but the relatively obvious differences are gradually revealed with the increase of the timescale of long-term orbital evolution due to model errors.

Zoom In Zoom Out Reset image size

To a certain extent, the generalized qualitative characteristic formulae of Equations (3) and (4) of the variation of orbital inclination-longitude of the ascending node for medium-high near-circular orbits, derived on the basis of the simplified and approximate model, are verified by the first-order semianalytical calculations based on the simplified models and the numerical calculations based on the precise models.