Analytic model for the long-term evolution of circular Earth satellite orbits including lunar node regression

TL;DR

This paper develops an analytic model to predict the long-term evolution of circular Earth satellite orbits, incorporating lunar node regression and gravitational perturbations, improving understanding of orbital dynamics over extended periods.

Contribution

It introduces an approximate analytic model that accounts for lunar orbital inclination and node regression, extending previous solutions for orbit evolution under Earth's $J_2$ and lunar effects.

Findings

Model agrees well with numerical results outside resonant cases.

Provides explicit analytic expressions for orbit evolution.

Highlights limitations near resonance conditions.

Abstract

This paper aims to obtain an analytic approximation to the evolution of circular orbits governed by the Earth's $J_2$ and the luni-solar gravitational perturbations. Assuming that the lunar orbital plane coincides with the ecliptic plane, Allan and Cook (1964) derived an analytic solution to the orbital plane evolution of circular orbits. Using their result as an intermediate solution, we establish an approximate analytic model with lunar orbital inclination and its node regression be taken into account. Finally, an approximate analytic expression is derived, which is accurate compared to the numerical results except for the resonant cases when the period of the reference orbit approximately equals the integer multiples (especially 1 or 2 times) of lunar node regression period.

Figures

Equations

$$\dot{\mathbf{R}}=\mathbf{R}\times \mathbf{MR}$$

$$\mathbf{M}=\sum _{j=0}^2{\omega }_j{\mathbf{Z}}_j{\mathbf{Z}}_j^T$$

$$\{\begin{array}{rl}& {\omega }_0=\frac{3}{2}\frac{n{J}_2{R}_e^2}{a^2}\\ & {\omega }_1=\frac{3}{4}\frac{{\mu }_S}{n{a}_S^3(1-e_S{)}^{3/2}}\\ & {\omega }_2=\frac{3}{4}\frac{{\mu }_M}{n{a}_M^3(1-e_M{)}^{3/2}}\end{array}$$

$$\{\begin{array}{rl}& {\mathbf{Z}}_0=(0,\sin \epsilon ,\cos \epsilon {)}^T\\ & {\mathbf{Z}}_1=(0,0,1{)}^T\\ & {\mathbf{Z}}_2=(\sin i_M\sin {\Omega }_M,-\sin i_M\cos {\Omega }_M,\cos i_M{)}^T\end{array}$$

$$\mathbf{M}={\mathbf{M}}_0+{\mathbf{M}}_1+{\mathbf{M}}_2$$

$${\mathbf{M}}_0=\left(\begin{array}{ccc}0& 0& 0\\ 0& {\omega }_0{\sin }^2\epsilon & {\omega }_0\sin \epsilon \cos \epsilon \\ 0& {\omega }_0\sin \epsilon \cos \epsilon & {\omega }_0{\cos }^2\epsilon +{\omega }_1^{\ast }\end{array}\right)$$

$${\mathbf{M}}_1=\frac{1}{2}{\omega }_2\sin 2i_M\left(\begin{array}{ccc}0& 0& \sin {\Omega }_M\\ 0& 0& -\cos {\Omega }_M\\ \sin {\Omega }_M& -\cos {\Omega }_M& 0\end{array}\right)$$

$${\mathbf{M}}_2=\frac{1}{2}{\omega }_2{\sin }^2{i}_M\left(\begin{array}{ccc}1-\cos 2{\Omega }_M& -\sin 2{\Omega }_M& 0\\ -\sin 2{\Omega }_M& 1+\cos 2{\Omega }_M& 0\\ 0& 0& 0\end{array}\right)$$

$${\lambda }_{2,3}=\frac{1}{2}({\omega }_0+{\omega }_1^{\ast })\left(1\mp \sqrt{1-4\frac{{\omega }_0{\omega }_1^{\ast }}{({\omega }_0+{\omega }_1^{\ast }{)}^2}{\sin }^2\epsilon }\right)$$

$$\{\begin{array}{rl}& {\mathbf{\eta }}_1=C_1{\left(\begin{array}{ccc}1& 0& 0\end{array}\right)}^T\\ & {\mathbf{\eta }}_2=C_2{\left(\begin{array}{ccc}0& \cos \theta & -\sin \theta \end{array}\right)}^T\\ & {\mathbf{\eta }}_3=C_3{\left(\begin{array}{ccc}0& \sin \theta & \cos \theta \end{array}\right)}^T\end{array}$$

$$\tan \theta =\frac{{\omega }_0\sin \epsilon \cos \epsilon }{{\lambda }_3-{\omega }_0{\sin }^2\epsilon }$$

$$\{\begin{array}{rl}& \dot{x}=({\lambda }_3-{\lambda }_2)yz\\ & \dot{y}=({\lambda }_1-{\lambda }_3)xz\\ & \dot{z}=({\lambda }_2-{\lambda }_1)xy\end{array}$$

$$\{\begin{array}{rl}& \mathbf{R}\cdot \mathbf{R}=x^2+y^2+z^2=1\\ & \mathbf{R}\cdot {\mathbf{M}}_0\mathbf{R}={\lambda }_1{x}^2+{\lambda }_2{y}^2+{\lambda }_3{z}^2\triangleq \lambda \end{array}$$

$$\{\begin{array}{rl}& x=\mathrm{sign}(x)\sqrt{\frac{{\lambda }_3-\lambda }{{\lambda }_3}}\mathrm{dn}u\\ & y=\sqrt{\frac{\lambda }{{\lambda }_2}}\mathrm{sn}u\\ & z=\sqrt{\frac{\lambda }{{\lambda }_3}}\mathrm{cn}u\end{array}$$

$$\kappa ={\kappa }_1=\sqrt{\frac{\lambda ({\lambda }_3-{\lambda }_2)}{{\lambda }_2({\lambda }_3-\lambda )}}\in [0,1]$$

$$u=u_0-\mathrm{sign}(x)\sqrt{{\lambda }_2({\lambda }_3-\lambda )}t$$

$$\{\begin{array}{rl}& x=\sqrt{\frac{{\lambda }_3-\lambda }{{\lambda }_3}}\mathrm{cn}u\\ & y=\sqrt{\frac{{\lambda }_3-\lambda }{{\lambda }_3-{\lambda }_2}}\mathrm{sn}u\\ & z=\mathrm{sign}(z)\sqrt{\frac{\lambda }{{\lambda }_3}}\mathrm{dn}u\end{array}$$

$$\kappa ={\kappa }_2=\sqrt{\frac{{\lambda }_2({\lambda }_3-\lambda )}{\lambda ({\lambda }_3-{\lambda }_2)}}\in [0,1]$$

$$u=u_0-\mathrm{sign}(z)\sqrt{\lambda ({\lambda }_3-{\lambda }_2)}t$$

$$\left\{\begin{aligned} \dot{x}\approx&(\lambda_{3}-\lambda_{2})yz+\epsilon\Big{[}xy\sin\Omega_{M}\cos\theta\\ &+xz\sin\Omega_{M}\sin\theta-y^{2}\cos\Omega_{M}\cos 2\theta\\ &-2yz\cos\Omega_{M}\sin 2\theta+z^{2}\cos\Omega_{M}\cos 2\theta\Big{]}\\ \dot{y}\approx&-\lambda_{3}xz-\epsilon\Big{[}x^{2}\sin\Omega_{M}\cos\theta\\ &+yz\sin\Omega_{M}\sin\theta-xy\cos\Omega_{M}\cos 2\theta\\ &-z^{2}\sin\Omega_{M}\cos\theta-xz\cos\Omega_{M}\sin 2\theta\Big{]}\\ \dot{z}\approx&\lambda_{2}xy-\epsilon\Big{[}x^{2}\sin\Omega_{M}\sin\theta\\ &-y^{2}\sin\Omega_{M}\sin\theta-xy\cos\Omega_{M}\sin 2\theta\\ &+yz\sin\Omega_{M}\cos\theta+xz\cos\Omega_{M}\cos 2\theta\Big{]}\end{aligned}\right.$$

$$\{\begin{array}{rl}& x=\Lambda (\kappa )\kappa \mathrm{cn}u\\ & y=\sqrt{\frac{{\lambda }_3}{{\lambda }_3-{\lambda }_2}}\Lambda (\kappa )\kappa \mathrm{sn}u\end{array}$$

$$w=\frac{\pi u}{2K(\kappa )}$$

$$\Lambda (\kappa )$$

$$\in \left(\sqrt{\frac{{\lambda }_3-{\lambda }_2}{{\lambda }_3}},\sqrt{\frac{{\lambda }_3-{\lambda }_2}{{\lambda }_2}}\right)$$

$$\left(\begin{array}{c}\dot{\kappa }\\ \\ \dot{u}\end{array}\right)$$

$$=\frac{1}{\det \left(\frac{\partial (x,y)}{\partial (\kappa ,u)}\right)}\left(\begin{array}{cc}\frac{\partial y}{\partial u}& -\frac{\partial x}{\partial u}\\ & \\ -\frac{\partial y}{\partial \kappa }& \frac{\partial x}{\partial \kappa }\end{array}\right)\left(\begin{array}{c}\dot{x}\\ \\ \dot{y}\end{array}\right)$$

$$\det\left(\frac{\partial(x,y)}{\partial(\kappa,u)}\right)=\sqrt{\frac{\lambda_{3}}{\lambda_{3}-\lambda_{2}}}\Big{[}1+\Lambda^{2}(\kappa)\kappa^{2}\Big{]}\Lambda^{2}(\kappa)\kappa\operatorname{dn}u$$

$$\dot{\kappa }=$$

$$+A_{0,1}(\kappa )\sin {\Omega }_M+A_{2,-1}(\kappa )\sin (2w-{\Omega }_M)$$

$$\displaystyle+A_{2,1}(\kappa)\sin(2w+\Omega_{M})+\mathcal{O}(\kappa^{2})\bigg{]}+\mathcal{O}(\epsilon^{2})$$

$$\dot{u}=-\mathrm{sign}(z)\sqrt{{\lambda }_2{\lambda }_3}\Lambda (\kappa )+\mathcal{O}(ϵ/\kappa )$$

$$\dot{w}=$$

Full text

00footnotetext: Purple Mountain Observatory, Chinese Academy Sciences, Nanjing, 21000800footnotetext: Key Laboratory of Space Object and Debris Observation, PMO, CAS, Nanjing, 210008

Keywords Orbital plane evolution; Milankovitch element; Lunar node regression; Analytic approximation

1 Introduction

The problem of the orbital plane evolution of a satellite around an oblate planet has been studied since two hundred years ago when Laplace investigated the motion of Iapetus, the third largest natural satellite of Saturn. Sixty years ago, the launch of Sputnik-1 boosted similar research in the field of artificial satellite theory. Allan and Cook (1964) presented their analytic results on the long period motion of the plane of distant circular orbits due to the Earth’s oblateness and the luni-solar gravitational perturbations. The dynamical system was doubly averaged over the mean motions of both the satellite and the disturbing bodies, and further simplified by the assumption that the moon is orbiting the Earth in the ecliptic plane. Then the classical Laplace plane was defined as the local equilibrium of the simplified dynamical system. In addition, according to the analytic results, the orbital plane of a circular orbit was precessing around the Laplace plane, with a period of several decades depending on the semi-major axis and the initial direction of the orbital pole.

A number of theories and applications were presented after the classical work of Allan and Cook (1964). Satellite dynamics on the Laplace surface of general Sun-Planet-Satellite systems were discussed in detail by Tremaine et al. (2009). The Laplace plane at the Geosynchronous Earth Orbit (GEO) altitude (i.e. $a\approx{42164}{\operatorname{\,km}}$) is approximately $7.4\operatorname{{}^{\circ}}$ inclined to the equator, which implies that the inclination of conventional GEO satellites ($a\approx{42164}{\operatorname{\,km}}$, $e\approx 0$ and $i\approx{0}\operatorname{{}^{\circ}}$) varies in the range of $[{0}{\operatorname{{}^{\circ}}},{15}{\operatorname{{}^{\circ}}}]$ with a period of about 53 years (Valk et al., 2009; Ulivieri et al., 2013; Zhao et al., 2015). The orbital plane cannot be frozen in the classical Laplace plane, and the periodicity of orbital plane evolution cannot reveal real motion, when the analysis is extended to Medium Earth Orbit (MEO), due to the influence of Moon’s real motion, especially the lunar orbital regression along the ecliptic plane (Ulivieri et al., 2013; Circi et al., 2016).

Allan and Cook (1964) warned that circular motion may be unstable for some inclinations due to the Kozai-Lidov effect or luni-solar resonance in general (Breiter, 2001). Recently, a series of papers showed that the eccentricity of the orbits for Global Navigation Satellite System (GNSS) may increase exponentially (Chao and Gick, 2004; Bordovitsyna et al., 2012). Through theoretical calculations and numerical verifications for the Chinese BeiDou Inclined Geosynchronous Earth Orbit (IGSO), we found that the eccentricity of the near-circular IGSO satellite with an initial inclination not larger than $32\operatorname{{}^{\circ}}$ can always remain small in the long-term evolution, and the initial Right Ascension of the Ascending Node (RAAN) can affect the eccentricity too (Zhao et al., 2015). Daquin et al. (2016b, a) and Gkolias et al. (2016) improved the investigation of eccentricity excitation by introducing dynamical indicators (e.g. Fast Lyapunov Indicator, or FLI) to describe the chaoticity of the dynamic systems for the highly inclined orbits. Their results verified that the instability is highly correlated with the orbital plane and demonstrated the complexity of the problem.

Therefore, an investigation of the orbital plane variation of circular orbits with the consideration of lunar node regression is crucial in understanding the orbital dynamics. This paper aims to obtain an analytic approximation to the evolution of circular orbits governed by the Earth’s oblateness and the luni-solar gravitational perturbations. The main body of this paper is organised as follows:

Section 2 briefly introduces the doubly averaged dynamic model in terms of vectors and tensors. Similar to Circi et al. (2016), the force model consists of the Earth’s $J_{2}$ perturbation and the luni-solar gravitational perturbations. The lunar node regression is also taken into consideration. Then the force model is divided into an autonomous part and a non-autonomous part (seen as the perturbation), where the former can be solved analytically by the method presented in Allan and Cook (1964). Treating their solution as an intermediate orbit, we introduce a variable transformation and establish the equation of motion using the method of constant variation. Section 3 presents the analytic solution of the doubly averaged model. A general expression describing the long-period motion of the orbital plane is derived from the equation of motion, which is suitable for non-resonant circular orbits with small or medium inclination, except for those close to the classical Laplace plane. The analytic approach for the quasi-frozen (i.e. near Laplace plane) orbits is accomplished based on the linearisation of the original equation. Section 4 applies massive numerical simulations to verify the analytic model.

2 Doubly averaged model

The long-term evolution of a near-circular orbit can be described by its orbital angular momentum, and the equation of motion takes the form of (Allan and Cook, 1964; Circi et al., 2016)

[TABLE]

where $\bm{R}$ is a unit vector along the orbital angular momentum, and $\bm{M}$ is a two dimensional symmetric tensor

[TABLE]

with $\bm{Z}_{0}$, $\bm{Z}_{1}$ and $\bm{Z}_{2}$ unit vectors along the Earth pole, the ecliptic pole and the lunar orbital angular momentum, respectively. $\bm{Z}_{j}\bm{Z}_{j}^{T}$ are symmetric dyads (refer to the appendices of Rosengren and Scheeres (2014)). The coefficients $\omega_{j}$’s take the following form

[TABLE]

where $R_{e}$ is the Earth’s equator radius; $\mu_{S}$ and $\mu_{M}$ are the gravitational constants of the sun and the moon, respectively; $a_{S}$, $e_{S}$, $a_{M}$ and $e_{M}$ are the Keplerian orbital elements of the sun and the moon in the Geocentric Ecliptic Reference System (GERS, $O\text{\,-\,}X_{1}Y_{1}Z_{1}$); $a$ is the orbital semi-major axis of the satellite and $n$ is the satellite’s mean motion. Obviously, $\omega_{j}$ are one variable functions of semi-major axis, $a$. In addition, in the doubly averaged dynamical system, $a$ is a first integral (or constant of motion), then $\omega_{j}$ can be seen as constants.

Given some reference system, the vectors $\bm{Z}_{j}$ can be expressed by $3\times 1$ matrices and thus the symmetric tensor $\bm{M}$ can be expressed by $3\times 3$ matrix. In GERS,

[TABLE]

where $\varepsilon$ is the obliquity of the ecliptic plane with respect to the equatorial plane, $i_{M}$ and $\Omega_{M}$ are the inclination and the RAAN of the lunar orbit in GERS. As $\sin i_{M}\approx 0.0894$ can be seen as small parameter, the symmetric matrix $\bm{M}$ can be expressed as follows:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

with $\omega_{1}^{*}=\omega_{1}+\omega_{2}\cos^{2}i_{M}$. The symmetric matrix $\bm{M}$ is composed of an autonomous part $\bm{M}_{0}$ and a time-dependent part $\bm{M}_{1}+\bm{M}_{2}$. The time-dependent part can be seen as perturbations with magnitudes dependent on $\sin i_{M}$.

2.1 Intermediate solution

The intermediate dynamical system is obtained by excluding the time-dependent tensor from Eq. (1). Similarly to Allan and Cook (1964), the intermediate solution is obtained through several steps listed below:

1. (1)

Eigenvalue and eigenvector analysis of $\bm{M}_{0}$

The Eigenvalues of $\bm{M}_{0}$ are $0=\lambda_{1}<\lambda_{2}<\lambda_{3}$

[TABLE]

and the corresponding eigenvectors in GERS are

[TABLE]

where $C_{j}\neq 0$ for $j=1,2,3$ are arbitrary constants, and $\theta$ is a function of semi-major axis determined by

[TABLE]

Fig. 1 illustrates the eigenvalues $\lambda_{2,3}$ and angle $\theta$ as functions of semimajor axis for readers’ convenience of reading off their magnitudes.

2. (2)

Determination of the intermediate reference system

Let $C_{1,2,3}=1$, then $\bm{\eta}_{1,2,3}$ form a right-hand Geocentric Intermediate Reference System (GIRS, $O\text{\,-\,}xyz$), in which $Ox$-axis points to the equinox, and $Oz$-axis points to the north with an angle $\theta$ from the ecliptic pole ($OZ_{1}$). Then the $O\text{\,-\,}xy$ plane is the local Laplace plane with an inclination $\alpha=\varepsilon-\theta$ to the equator. 3. (3)

Solve the equation of motion in the intermediate reference system analytically

In the $O\text{\,-\,}xyz$ reference system, the equation of motion is simplified to

[TABLE]

where $(x,y,z)^{T}$ is the coordinate of $\bm{R}$ in GIRS, and note that $\lambda_{1}=0$. There are two first integrals of the intermediate dynamical system:

[TABLE]

Finally, the solution is expressed analytically in two cases, depending on $\lambda$:

1. (a)

If $\lambda\in[0,\lambda_{2}]$, then

[TABLE]

where $\operatorname{dn}$, $\operatorname{sn}$ and $\operatorname{cn}$ are Jacobi elliptic functions with $\kappa$ the elliptic modulus defined as follows:

[TABLE]

and $u$ a linear function in time

[TABLE]

where $u_{0}$ is the initial value of $u$. The function $\operatorname{sign}(x)$ equals to 1 for positive $x$, $-1$ for negative $x$, and zero for $x=0$, and it is determined by the initial value and does not change along the solution trajectory. 2. (b)

If $\lambda\in[\lambda_{2},\lambda_{3}]$, then

[TABLE]

and the corresponding $(\kappa,u)$ are

[TABLE]

[TABLE]

2.2 Geometric view of the intermediate solution

As showed in Fig. 2, the pole of the orbital plane precesses around the axis $Ox$ or $Oz$ for case (a) and (b) respectively, and they are separated by the separatrices (dash-dot line, corresponds to $\lambda=\lambda_{2}$) connected to the $Oy$-axis. In addition, common orbits with small and medium inclination (with respect to the equator) are of case (b).

2.3 Equation of perturbed motion

The equation of motion considering the time-dependent tensor can be expressed in GIRS as follows:

[TABLE]

where $\epsilon=\omega_{2}\sin i_{M}\cos i_{M}$, and the second order tensor $\bm{M}_{2}$ is ignored.

The phase space of system (20) is two dimensional considering that the first integral $\bm{R}\cdot\bm{R}=1$ still holds. Then the equation of the perturbed motion is established based on the theory of constant variation. Specifically, $\kappa$ is a first integral (or constant of motion) of the unperturbed dynamical system, and it determines the integral curve (i.e. solid circles with an arrow in Fig. 2), while variable $u$ (linear function in time) determines the location on the curve at a specific time. Taking into account of the perturbation, $\kappa$ turns to a slowly changing variable, and $u$ becomes quasi-linear with time. According to the variable transformation $(x,y,z)\to(\lambda,u)\to(\kappa,u)$, the equation of the perturbed motion is derived easily, as outlined below.

Concentrating on case (b), the explicit expression of the transformation $(x,y)\to(\kappa,u)$ is

[TABLE]

and the phase angle $w$ is defined by

[TABLE]

where $K$ is the complete elliptic integral, and $\Lambda$ is the function of $\kappa$ ($\kappa\in(0,1)$) in the following form:

[TABLE]

The equation after transformation (21) is defined by

[TABLE]

where

[TABLE]

is the Jacobi determinant, and the factor $\kappa$ leads to vanishing divisor at $\kappa=0$ for the expression of $\dot{u}$ (see below).

After a series of cumbersome derivations and introducing the Fourier series of elliptic functions, Eq. (24) reduces to Eqs. (26) and (27):

[TABLE]

and

[TABLE]

For convenience, we replace $u$ with the phase angle $w$ as an independent variable, then

[TABLE]

Note that the neglected terms $\mathcal{O}(\kappa^{2})$ in Eq. (26) are not all the second order terms of $\kappa$, but only those from the Fourier series of Jacobi’s elliptic functions, which are of smaller amplitudes and higher frequencies than the leading terms. In addition, the changing rate of angle $u$ or $w$ is dominated by the zeroth order term of $\epsilon$ (zeroth order secular term) and the neglected terms in $\mathcal{O}(\epsilon/\kappa)$ are first order periodic terms (with the period of 18.6 years). Due to the singularity at $\kappa=0$, the above simplifications are valid for $\kappa$ that is not so small according to our numerical tests. And analytic solution for extremely small $\kappa$ is developed in another form in order to avoid the small divisor problem.

By introducing the truncated power series of elliptic integral $K$, the amplitudes $A_{i,j}$’s in Eq. (26) can be approximated expressed as follows:

[TABLE]

We do not further simplify the amplitudes by expanding $\Lambda$ and $\sqrt{1-\kappa^{2}\Lambda^{2}}$, because the ratio of eigenvalues, $\frac{\lambda_{2}}{\lambda_{3}}$, is small, to the extent that the expansion would cause severe loss of accuracy.

2.4 Model validation

The validation of the doubly averaged model in studying the long-term evolution of distant circular orbits has been confirmed in a number of papers (Allan and Cook, 1964; Rosengren et al., 2014; Ulivieri et al., 2013). We verify the model by comparing the time series of inclination due to Eq. (1) and the results from the STELA software (CNES, 2015). The Earth’s zonal harmonics up to 4th degree ($J_{2}$, $J_{3}$ and $J_{4}$ terms) and the luni-solar gravity (based on the simplified Meeus and Brown theory, a highly accurate analytic model) are taken into consideration in the STELA simulation. Fig. 3 illustrates three of our verification examples, the initial (at epoch J2000.0) mean Keplerian elements are $a={26554}{\operatorname{\,km}}$, $e=0$, $i={15}{\operatorname{{}^{\circ}}},{55}{\operatorname{{}^{\circ}}},{75}{\operatorname{{}^{\circ}}}$, and $\Omega={0}{\operatorname{{}^{\circ}}}$.

3 Analytic approach

3.1 General solution for the perturbed motion

Observing the Eqs. (26) and (28), it is obvious that $\kappa$ varies quasi-periodically with the amplitude of the order of $\mathcal{O}(\epsilon)$. Then it is reasonable to expand the right-hand side of both the equations in the vicinity of some mean value of $\kappa$ (denoted by $\bar{\kappa}$ and to be determined later) and retaining the leading terms only:

[TABLE]

and

[TABLE]

Denote the variables’ initial values by $w_{0}$ and $\kappa_{0}$, the approximate solutions are obtained as

[TABLE]

and

[TABLE]

where $\bar{B}_{i,j}=B_{i,j}(\bar{\kappa})=A_{i,j}(\bar{\kappa})/(i\dot{w}(\bar{\kappa})+j\dot{\Omega}_{M})$ are constants. The mean $\kappa$ is obtained by removing all the periodic terms from the osculating value at $t=t_{0}$, i.e.

[TABLE]

with $\Omega_{M}^{0}=\Omega_{M}(t_{0})$ the initial value. In Eq. (34), $\bar{B}_{i,j}$’s can be replaced by $B_{i,j}(\kappa_{0})$ to the first order.

Finally, the analytic model for the perturbed motion in terms of variables $(\kappa,w)$ is expressed by their approximate solutions (Eqs. (32) and (33)) along with Eq. (34).

Remark 1**.**

The assumption on the variation range of $\kappa$ does not hold when resonance happens, i.e. when $i\dot{w}+j\dot{\Omega}_{M}\approx 0$. As a result, the approximate general analytic solution fails for resonant cases.

3.2 Evolution near the classical Laplace plane

As presented in section 2.3, the equation of the perturbed motion may not be appropriate for the orbits in the neighborhood of classical Laplace plane, i.e. when $\kappa\approx 0$. Traditionally, the circular orbits on the Laplace plane are frozen (i.e. the orbital plane has no secular change). As a result, here we call the orbits near the classical Laplace plane the quasi-frozen orbits.

The analytic approach on the quasi-frozen orbits stands on linearizing Eq. (20) at $\kappa=0$. Note again that in this paper, we consider case (b) only. Then it reduces to

[TABLE]

Let $p=\sin\Omega_{M}$ and $q=\cos\Omega_{M}$, neglect higher-order terms:

[TABLE]

The eigenvalues of the coefficient matrix are \big{\{}\pm\operatorname{i}\dot{\Omega}_{M},\pm\operatorname{i}\nu_{2}\big{\}}, where $\nu_{2}=\sqrt{\lambda_{3}(\lambda_{3}-\lambda_{2})}$. If $\dot{\Omega}_{M}^{2}\neq\nu_{2}^{2}$, the solution takes the following form:

[TABLE]

The coefficients in Eq. (37) satisfy

[TABLE]

where $(x^{0},y^{0},z^{0})$ and $\Omega_{M}^{0}$ stand for the initial values. Then

[TABLE]

Finally, the orbital plane evolution of the quasi-frozen orbits for case (b) is approximately described by the analytic expression of Eqs. (37), (39), and

[TABLE]

4 Numerical tests on the analytic model

The analytic expression is verified by comparing with numerical solution of the doubly averaged model, i.e. Eq. (1).

Both the numerical and analytic solutions are converted back to Keplerian elements for the clarity of physical meaning. Then the time series of deviation in the inclination and RAAN (noted by $\Delta i(t)$ and $\Delta\Omega(t)$ respectively) of the analytic results from the numerical ones is obtained. At last, we calculate the root mean square (RMS) of $\Delta X$ by

[TABLE]

with $X$ either the inclination or RAAN, $N=5000$ in our simulation, and $t_{k}=k\times 0.1$ year.

For a given semi-major axis, the comparison sample set is generated as follows:

• •

Initial inclination: $i\in({0}{\operatorname{{}^{\circ}}},{75}{\operatorname{{}^{\circ}}})$ with the step size being ${1}{\operatorname{{}^{\circ}}}$

• •

Initial RAAN: $\Omega\in({0}{\operatorname{{}^{\circ}}},{360}{\operatorname{{}^{\circ}}})$ with the step size being ${10}{\operatorname{{}^{\circ}}}$

• •

The initial epoch is fixed at J2000.0, so $\Omega_{M}$ is not sampled to cover $({0}{\operatorname{{}^{\circ}}},{360}{\operatorname{{}^{\circ}}})$.

Then the Keplerian elements for the circular orbits of each sample are converted to vector expressions in GCRS (Geocentric Celestial Reference System), GERS, and GIRS.

The numerical integration of Eq. (1) is applied in GERS, using the classical Runge-Kutta-Dormand-Prince 8(9) method. The time range of the integration is 500 years.

A two dimensional scatter is obtained for each semi-major axis (see Figs. 4 and 5).

The general analytic solution have comparatively large error when $\kappa\approx 0$. This is rather obvious for the cases $a=26554\operatorname{\,km}$ and $a=30000\operatorname{\,km}$. For small $\kappa$, we introduce the quasi-frozen analytic expression, i.e. Eqs. (37) and (39). The mean errors in inclination and RAAN for both analytic solutions are illustrated in Figs. 6 and 7, which show a better accuracy of quasi-frozen solution compared to the general analytic approach for small $\kappa$.

To note that the y-axis in Figs. 4 and 5 are limited to no larger than ${3}{\operatorname{{}^{\circ}}}$ and ${30}{\operatorname{{}^{\circ}}}$ respectively, then some points are removed from the first three plots (specifically, points at $\kappa\approx 0.14$ for $a=20000{\operatorname{\,km}}$ , at $\kappa\approx 0.07$ and $\kappa\approx 0.2$ for $a=26554{\operatorname{\,km}}$ , at $\kappa\approx 0.15$ for $a=30000{\operatorname{\,km}}$), which corresponds to the resonances as mentioned in Remark 1.

An extended simulation is applied to cover the principle medium and high Earth orbits. The samples are generated in the three dimensional initial value space:

• •

$a\in({20000}{\operatorname{\,km}},{45000}{\operatorname{\,km}})$ with the step size of ${500}{\operatorname{\,km}}$

• •

$i\in({0}{\operatorname{{}^{\circ}}},{75}{\operatorname{{}^{\circ}}})$ with the step size of ${1}{\operatorname{{}^{\circ}}}$

• •

$\Omega\in({0}{\operatorname{{}^{\circ}}},{360}{\operatorname{{}^{\circ}}})$ with the step size of ${10}{\operatorname{{}^{\circ}}}$

• •

The initial epoch is fixed at J2000.0

Fig. 8 illustrates the histogram of errors in inclination (above), and its relation with the basic period by scatter plot (bottom). The errors in RAAN are showed in Fig. 9. The histograms show that the percentage of samples with \big{\|}\Delta i(t)\big{\|}_{2}<{0.6}{\operatorname{{}^{\circ}}} and \big{\|}\Delta\Omega\|_{2}<{25}{\operatorname{{}^{\circ}}} is about ${90}{\%}$ among our sample set. The bottom plots reveal a strong correlation between the failing analytic model and the resonances at $T_{w}\approx nT_{M}$ for $n=1,2$, where $T_{w}=2\pi/\dot{w}$ (depends on $a$ and $\kappa$) is the basic period, and $T_{M}\approx{18.6}$ year is the regression period of the lunar node.

5 Conclusion and discussion

We have established an approximate analytic theory of the orbital plane evolution of circular Earth satellite orbits with high altitude. The force model consists of the Earth’s oblateness perturbation, and the luni-solar gravitational perturbations, with lunar orbital inclination to the ecliptic plane and the node regression taken into consideration. A doubly averaged model is established using the Milankovitch elements. The approximate analytic expression for the long period orbital plane evolution is derived using the method of the variation of constants, on the basis of the classical analytic solution, in which the lunar node regression was omitted. Appropriate amendments have been applied to the orbits close to the classical Laplace plane. According to our analysis and verification, the analytic approach is valid for most medium and high Earth orbits with small eccentricity, except for the resonant cases.

Acknowledgements This research was supported by the Youth Program of the National Natural Science Foundation of China (Grant No. 11603078) and the Key Program of the National Natural Science Foundation of China (Grant No. 11533010). The authors would like to thank two anonymous reviewers for the valuable comments that help to substantially improve the manuscript. The authors would also like to thank Dr. Fangzhou Jiang for his help in improving the language. The first author would like to give special thanks to his bride, Q. Qian, for her endless support.