Analytical solutions of the relative orbital motion in unperturbed and in J₂ - perturbed elliptic orbits

Solving and modeling the relative motion problem between satellites or space crafts is of great importance in the field of formation flying, rendezvous and disturbed satellite systems which in turn play a significant rule in space missions. Analytical expressions for propagating relative motion under the influence of J₂ have been developed by several authors. Tschauner-Hempel [1] derived equations of relative motion valid for vehicles in nearby eccentric orbits. Melton [2] presents a time-explicit solution for relative motion for elliptic orbits using a perturbation method . Valaddi et al. [3] also presents a time-explicit perturbation solution to the relative motion in elliptic orbits by retaining quadratic nonlinearities of the gravitational field. Alternatively, Gim and Alfriend [4] provide a state transition matrix for the solution of the linearized relative motion problem using curvilinear coordinates. Garrison et al. [5] consider a geometric method for deriving the state transition matrix, utilizing small differences in orbital elements between two satellites. Also Srinivas R. Vadali [6] uses the geometric method but under the influence of J₂-perturbation. Schaub [7] presents analytical expressions for propagating the linearized orbital element differences using the true anomaly as the independent variable. Time-explicit solutions have been developed by Sabol et al. [8] using eccentricity expansions. This approach is valid for small eccentricities.

In this research article, the true anomaly of the Chief is used as the independent variable, rather than time in the modelling of the relative equations of motion of a deputy with respect to a chief satellite, such that both of them rotate under the gravitational field of the Earth in elliptic orbits, considering the effect of the Earth oblateness J₂ gravity perturbation into our calculations. Then we use Laplace transformation to solve the equations of motion. Finally we will illustrate the effect of the perturbation part using numerical example.

In order to describe the relative motion of the deputy satellite with respect to the chief satellite, two Cartesian coordinate systems are used as shown in Fig. 1. The first is the Earth centered inertial (ECI) coordinate frame (X, Y, Z) which is spanned by the unit vectors (I, J, K). In this system the X-axis points toward the vernal equinox, the Z-axis points toward the North-pole of the Earth, and the Y-axis is perpendicular to the other two and completes the right-handed coordinate system. The plane determined by I and J coincides with the equatorial plane. The chief orbit is defined in the ECI frame using orbital elements (r_C, θ, i, ω, Ω), which include geocentric distance r, argument of latitude θ, inclination i, argument of the peri-apse ω, and the right ascension of the ascending node Ω.

Fig. 1

The second frame is the Local-Vertical-Local-Horizontal (LVLH) frame (ξ, η, ζ) which is attached to the chief satellite S₀ and spanned by the unit vectors e_r,e_θ,e_h) $(\underline{e}_{r},\underline{e}_{\theta},\underline{e}_{h})$ where e_r $\underline{e}_{r}$ lies in the chief’s radial direction from the Earth center, e_h $\underline{e}_{h}$ is parallel to the orbit angular momentum vector, and e_θ $\underline{e}_{\theta}$ completes the right-handed orthogonal triad. Thus,

e_r=r_c|r_c|,e_h=h_c|h_c|,e_θ=e_h∧e_r$$ \begin{equation} \label{basic} \underline{e}_{r} \, =\, \frac{\underline{r}_{c}}{\left|\underline{r}_{c}\right|} ,\, \, \, \, \underline{e}_{h } =\frac{\underline{h}_{c} }{\left|\underline{h}_{c} \right|} \, \, ,\, \, \, \, \, \, \underline{e}_{\theta } =\underline{e}_{h} \wedge \underline{e}_{r} \end{equation} $$(1)

Such that that h_C=r_C∧v_C $\underline{h}_{C} =\underline{r}_{C} \wedge\underline{v}_{C}$ , is the chief’s angular momentum vector per unit mass in its orbit around the central body (Earth), and

rC=p1+ecosf,p=a(1−e2)$$ \begin{equation} r_{C} =\frac{p}{1+e\, \cos f} \, \, \, ,\, \, \, p=a\left(1-e^{2} \right) \end{equation} $$(2)

where p, a and e are the semi-latus rectum, semi-major axis and the eccentricity of the chief orbit, respectively.

The position vector of the deputy satellite relative to the chief satellite can be expressed by Cartesian coordinate vector σ_ $\mathrm{\sigma}$ :

σ_=[ξ η ζ]T$$ \begin{equation} \underline{\sigma}=\left[\xi \, \, \, \, \eta \, \, \, \, \zeta \right]^{T} \end{equation} $$(3)

Where ξ, η and ζ denote the components of the position vector along the radial, transverse, and the normal to the orbital plane of the chief directions respectively. The position vectors of the chief and deputy with respect to the Earth center are given by r→C and r→D $\vec{r}_{C} \, and\, \vec{r}_{D}$ respectively. Assuming a central gravity field, the chief frame rotates with an angular velocity ω_=[0 0 θ˙]T $\underline{\omega}=\left[0\, \, \, 0\, \, \, \dot{\theta}\right]^{T}$ , where θ=ω+f, is the argument of latitude, ω is the argument of periapsis and f is the true anomaly. Under the assumption of a central gravity field, ω is constant, and θ˙=f˙ $\dot{\theta}\, =\, \dot{f}$ .

From Figure 1, the deputy satellite position vector with respect to the Earth can be written as:

r→D=r→C+σ→=(rC+ξ)er+ηeθ+ζeh$$ \begin{equation} \vec{r}_{D} =\, \, \vec{r}_{C} +\vec{\sigma}=(r_{C} +\xi )\, e_{r} +\eta \, e_{\theta} +\zeta \, e_{h} \end{equation} $$(4)

The equations of motion of the chief and deputy satellites with respect to the ECI-frame can be written respectively as:

r→¨C=−μrC3r→C+a→C(r→C)$$ \begin{equation} \ddot{\vec{r}}_{C} \, =-\frac{\mu}{r_{C}^{3}} \, \vec{r}_{\, C} +\vec{a}_{C} \left(\, \vec{r}_{C} \right) \end{equation} $$(5)r→¨D=−μrD3r→D+a→D (r→D)$$ \begin{equation} \ddot{\vec{r}}_{D} =-\frac{\mu}{r_{D}^{3}} \vec{r}_{\, D} +\vec{a}_{D} \left(\vec{r}_{\, D} \right) \end{equation} $$(6)

where a→C and a→D $\vec{a}_{\, C} \, \, and\, \, \vec{a}_{D}$ are the perturbing accelerations acting on the chief and deputy satellites due to the Earth oblateness.

From equation 4, we can write

r→¨D=r→¨C+σ→¨ECI$$ \begin{equation} \ddot{\vec{r}}_{D} =\ddot{\vec{r}}_{C} +\ddot{\vec{\sigma}}_{ECI} \end{equation} $$(7)

Now we need to transfer the acceleration σ→¨ $\ddot{\vec{\sigma}}$ in the LVLH - frame to the acceleration in the ECI-frame (inertial frame). The transformation is well known, for example see Wiesel’s book [9].

σ→¨ECI=σ→¨LVLH+2ω→∧σ→˙LVLH+ω→˙∧σ→LVLH+ω→∧(ω→∧σ→LVLH)$$ \begin{equation} \ddot{\vec{\sigma}}_{ECI} =\ddot{\vec{\sigma}}_{LVLH} +2\, \vec{\omega}\wedge \dot{\vec{\sigma}}_{LVLH} +\dot{\vec{\omega}}\wedge \vec{\sigma}_{LVLH} +\vec{\omega}\wedge \left(\vec{\omega}\wedge \vec{\sigma}_{LVLH} \right) \end{equation} $$(8)

Where ω→ $\vec{\omega}$ is the angular velocity of the leader satellite in the LVLH - frame with respect to the ECI -frame, where

ω→=[00f˙]Tω→˙=[00f¨]Tr→C=[rC00]T$$ \begin{equation} \vec{\omega}=\left[\begin{array}{ccc} {0} & {0} & {\dot{f}} \end{array}\right]^{T} \,\,\,\, \dot{\vec{\omega}}=\left[\begin{array}{ccc} {0} & {0} & {\ddot{f}} \end{array}\right]^{T} \,\,\,\, \vec{r}_{C} =\left[\begin{array}{ccc} {r_{C}} & {0} & {0} \end{array}\right]^{T} \end{equation} $$(9)

Using equations 4, 5, 6, 7, 8, to 9, we can write

−μrD3r→D+a→D(r→D)=a→C(r→C)+(ξ¨−2η˙f˙−ηf¨−ξf˙2−μrC2)e→r+(η¨+2ξ˙f˙+ξf¨−ηf˙2)e→θ+ζ¨e→h$$ \begin{equation} \begin{array}{rl} -\frac{\mu}{r_{D}^{3}} \vec{r}_{\, D} +\vec{a}_{D} \left(\vec{r}_{\, D} \right)=&\vec{a}_{C} \left(\vec{r}_{C} \right)+\left(\ddot{\xi}-2 \dot{\eta}\dot{f}-\eta \ddot{f}-\xi \dot{f}^{2} -\frac{\mu}{r_{C}^{2}} \right)\vec{e}_{r}\\ &+\left(\ddot{\eta}+2\dot{\xi}\, \dot{f}+\xi \ddot{f}-\eta \, \dot{f}^{2} \right)\, \vec{e}_{\theta} +\ddot{\zeta} \, \vec{e}_{h} \end{array} \end{equation} $$(10)

In order to simplify the previous equation we will eliminate f¨ $\ddot{f}$ using the equation of motion per unit mass of the chief in its Keplerian motion around the Earth in the form

r→¨C=(r¨C−rCf˙2)e→r+(rCf¨+2r˙Cf˙)e→θ=−μrC2e→r$$ \ddot{\vec{r}}_{C} =(\ddot{r}_{C} -r_{C} \dot{f}^{2} )\vec{e}_{r} +(r_{C} \ddot{f}+2\dot{r}_{C} \dot{f})\vec{e}_{\theta} =-\frac{\mu}{r_{C}^{2}} \vec{e}_{r} $$

Equating the components on both sides, we obtain the following relations:

r¨C=rCf˙2−μrC2,f¨=−2r˙CrCf˙$$ \begin{equation} \ddot{r}_{C} =r_{C} \dot{f}^{2} -\frac{\mu}{r_{C}^{2}} \, \, \, \, \, \, \, \, \, ,\, \, \, \, \, \, \, \, \, \ddot{f}=-2\frac{\dot{r}_{C}}{r_{C}} \dot{f} \end{equation} $$(11)

Using equations 4 and 11, we can say that

a→D(r→D)−a→C(r→C)=(ξ¨−2η˙f˙+2r˙CrCηf˙−ξf˙2−μrC2+μrD3(rC+ξ))e→r+(η¨+2ξ˙f˙−2r˙CrCξf˙−ηf˙2+μrD3η)e→θ+(ζ¨+μrD3ζ)e→h$$ \begin{equation} \begin{array}{rl} \vec{a}_{D} \left(\vec{r}_{\, D} \right)-\vec{a}_{C} \left(\, \vec{r}_{C} \right)=&\left(\ddot{\xi}-2\, \dot{\eta}\dot{f}+2\frac{\dot{r}_{C}}{r_{C}} \eta \dot{f}-\xi \, \dot{f}^{2} -\frac{\mu}{r_{C}^{2}} +\frac{\mu}{r_{D}^{3}} \left(r_{C} +\xi \right)\right)\vec{e}_{r}\\ &+ \left(\ddot{\eta}+2\dot{\xi}\, \dot{f}-2\frac{\dot{r}_{C}}{r_{C}} \xi \dot{f}-\eta \, \dot{f}^{2} +\frac{\mu}{r_{D}^{3}} \eta \right) \vec{e}_{\theta} +\left(\ddot{\zeta}+\frac{\mu}{r_{D}^{3}} \zeta \right)\, \vec{e}_{h} \end{array} \end{equation} $$(12)

Or we can write 12 in its components form as the following

ξ¨−2(η˙−r˙CrCη)f˙−ξf˙2−μrC2+μrD3(rC+ξ)=(aD)r−(aC)rη¨+2(ξ˙−r˙CrCξ)f˙−ηf˙2+μrD3η=(aD)θ−(aC)θζ¨+μrD3ζ=(aD)h−(aC)h$$ \begin{equation} \begin{array}{l} {\ddot{\xi}-2(\dot{\eta}-\frac{\dot{r}_{C}}{r_{C}} \eta )\dot{f}-\xi \, \dot{f}^{2} -\frac{\mu}{r_{C}^{2}} +\frac{\mu}{r_{D}^{3}} (r_{C} +\xi )=\left(a_{D} \right)_{r} -\left(a_{C} \right)_{r}} \\ {\ddot{\eta}+2(\dot{\xi}-\frac{\dot{r}_{C}}{r_{C}} \xi )\dot{f}-\eta \dot{f}^{2} +\frac{\mu}{r_{D}^{3}} \eta =\left(a_{D} \right)_{\theta} -\left(a_{C} \right)_{\theta}} \\ {\ddot{\zeta}+\frac{\mu}{r_{D}^{3}} \zeta =\left(a_{D} \right)_{h} -\left(a_{C} \right)_{h}} \end{array} \end{equation} $$(13)

In order to write these equations with dimensionless coordinates, we introduce

[ξ η ζ]T=rC[xyz]T$$ \begin{equation} \left[\xi \, \, \, \, \eta \, \, \, \, \zeta \right]^{T} =r_{C} \left[x\, \, \, \, y\, \, \, \, z\right]^{T} \end{equation} $$(22)

And by these relations, we have

ξ ′ = r′_Cx+r_Cx′ and ξ ″ = r″_Cx + 2r′_Cx′+r_Cx″, with similar relations for η and ζ

After many simplifications including the perturbed part we can write 21 as

[x″y″z″]T=[2y′+31+ecosfx−2x′−z]+k(1+ecosf)[brbθbh]T$$ \begin{equation} \left[\begin{array}{ccc} {x''} & {y''} & {z''} \end{array}\right]^{T} =\left[\begin{array}{c} {2y'+\frac{3}{1+e\cos f} x} \\ {-2x'} \\ {-z} \end{array}\right]+k\left(1+e\, \cos f\right)\left[\begin{array}{ccc} {b_{r}} & {b_{\theta}} & {b_{h}} \end{array}\right]^{T} \end{equation} $$(23)

Where

k=−3J2Re22a2(1−e2)2[brbθbh]T=[1−3sin2isin2fsin2isin2fsin2isinf]$$ \begin{equation} k=\frac{-3\, J_{2} R_{e}^{2}}{2\, a^{2} \left(1-e^{2} \right)^{2}} \, \left[\begin{array}{ccc} {b_{r}} & {b_{\theta}} & {b_{h}} \end{array}\right]^{T} =\, \left[\begin{array}{c} {1-3\sin ^{2} i\, \, \sin ^{2} f} \\ {\sin ^{2} i\, \, \sin 2f\,} \\ {\sin 2i\, \, \sin f} \end{array}\right] \end{equation} $$(24)

Which is the final perturbed form of the system of linear relative equations of motion of a deputy with respect to a chief satellite, using the chief’s true anomaly as the independent variable.

Starting from the second equation of 23, which is

y″=−2x′+kA(1+ecosf)sin2f, where A=sin2i$$ y''=-2x'+k\, A\left(1+e\, \cos f\right)\, \, \, \sin 2f\, ,\, \, where\, A=\sin ^{2} i $$

Integrating with respect to the true anomaly , , we obtain:

y′=−2x+kA(sin2f−2e3cos3f)+d, where d=y′0+2x0+23ekA$$ \begin{equation} y'=-2x+k\, A\, \, (\sin ^{2} f-\frac{2e}{3} \cos ^{3} f)+d\, ,\, where\, \, \, d=y'_{0} +2x_{0} +\frac{2\,}{3} e\, k\, A \end{equation} $$(25)

Substitution into the first equation in 23 gives

(1+ecosf)x″=−(1+4ecosf)x+k(1+ecosf)2[1−3sin2isin2f]+2ksin2i(1+ecosf)(sin2f−2e3cos3f)+2d(1+ecosf)$$ \begin{equation} \begin{array}{rl} (1+e \cos f) x''=&-(1+4 e \cos f)x +k(1+e \cos f)^{2} \left[1-3 \sin ^{2} i \sin ^{2} f \right] \\ &+2k \sin ^{2} i (1+e\cos f)(\sin ^{2} f-\frac{2e}{3} \cos ^{3} f)+2d (1+e\cos f) \\ \end{array} \end{equation} $$(26)

Which can be written as

(1+ecosf)x″=−(1+4ecosf)x+k12(1+ecosf)(12−6A+3e(4−7A)cosf+6Acos2f+5eAcos3f)+2d(1+ecosf)$$ \begin{equation} \begin{array}{rl} (1+e\cos f)x''=&-(1+4 e\cos f)x +\frac{k}{12} \left(1+e\cos f\right)\left(12\right. -6A+3e\left(4-7A\right)\cos f \\ &+6A\cos 2f+5 e\, A\left. \cos 3f\right)+2 d\, (1+e\cos f) \end{array} \end{equation} $$(27)

Applying Laplace transformation on 27, such that l{g(f)} = G(s), where g(f), is any function of f We can construct table 1., taking into account the well known relation cosf=eif+e−if2 $\cos f=\frac{e^{if} +e^{-if}}{2}$

Table 1

Laplace transformation on the relative equation of motion of the Deputy in x-direction

And for the perturbed terms that contain k, we have:

ℒ{k12(1+ecosf)(12−6A+3e(4−7A)cosf+6Acos2f+5eAcos3f)}=k24(12(2−A)+3e2(4−7A)s+48e(1−A)s1+s2+4(e2(3−4A)+3A)s4+s2+16eAs9+s2+5e2As16+s2)$$ \begin{equation} \begin{array}{rl} {\rm \mathscr{L}}\left\{\frac{k}{12} \left(1+e\cos f\right)\left(12\right. -6A+3e\left(4-7A\right)\cos f+6A\cos 2f+5\, e\, A\left. \cos 3f\right)\right\} \\ =\frac{k}{24} \left(\frac{12\left(2-A\right)+3e^{2} \left(4-7A\right)}{s} +\frac{48e\left(1-A\right)s}{1+s^{2}} +\frac{4\left(e^{2} \left(3-4A\right)+3A\right)s}{4+s^{2}} +\frac{16eAs}{9+s^{2}} +\frac{5e^{2} As}{16+s^{2}} \right) \end{array} \end{equation} $$(28)

After some simplifications, we can get

X(s)=−[(s−i)2+42(s−i)(s+i)]eX(s−i)−[(s+i)2+42(s−i)(s+i)]eX(s+i)+2ds(s2+1)+(1+e)x0ss2+1+(1+e)x′0s2+1+k24(1+s2)[12(2−A)+3e2(4−7A)s+48e(1−A)s1+s2+4(e2(3−4A)+3A)s4+s2+16eAs9+s2+5e2As16+s2]$$ \begin{equation} \begin{array}{c} X(s)=-\left[\frac{(s-i)^{2} +4}{2(s-i)(s+i)} \right]e\, X(s-i)-\left[\frac{(s+i)^{2} +4}{2(s-i)(s+i)}\right]e\,X(s+i)+\frac{2d}{s(s^{2} +1)}\\ +\frac{(1+e)\, x_{0} \, s}{s^{2} +1} +\frac{(1+e)x'_{0}}{s^{2} +1} +\frac{k}{24\left(1+s^{2} \right)} \left[\frac{12\left(2-A\right)+3e^{2} \left(4-7A\right)}{s}\right.\\ +\frac{48e\left(1-A\right)s}{1+s^{2}} +\frac{4\left(e^{2} \left(3-4A\right)+3A\right)s}{4+s^{2}} +\frac{16eAs}{9+s^{2}} +\left.\frac{5e^{2} As}{16+s^{2}} \right] \end{array} \end{equation} $$(29)

And now applying the inverse Laplace transformation, recalling that

ℒ−1{G(s−a)}=eafg(f)$$ {\rm \mathscr{L}}^{-1} \{G(s-a)\} =e^{af} g(f) $$

We can get after simplifying

2e∫0fx(τ)dτ+(cscf+ecotf)x=2dcscf+[(1+e)x0−2d]cotf+[(1+e)x′0+edf]+k72sinf[9(4(2−A)+(4−7A)e2)+6(A(8+e+8e2)−4(3+e2))cosf+4((4A−3)e2−3A)cos2f−Ae(6cos3f+ecos4f)+72(1−A)efsinf]$$ \begin{equation} \begin{array}{rl} 2\, e\int _{0}^{f}x(\tau )\, d\tau +(\csc f+e\cot f)\, x=2d\csc f+[(1+e)x_{0} -2d]\cot f+[(1+e)x'_{0} +e\, d\, f] \\ +\frac{k}{72\sin f} \left[9\left(4\left(2-A\right)+\left(4-7A\right)e^{2} \right)\right. +6\left(A\left(8+e+8e^{2} \right)-4\left(3+e^{2} \right)\right)\cos f \\ +4\left(\left(4A-3\right)e^{2} -3A\right)\cos 2f-Ae\left(6\cos 3f+e\cos 4f\right)+\left. 72\left(1-A\right)e\, f\sin f\right] \end{array} \end{equation} $$(30)

By differentiating this equation with respect to f, we can get

x′+[2e−cscfcotf−ecsc2fcscf+ecotf]x=ed−2dcscfcotf−[(1+e)x0−2d]csc2fcscf+ecotf+ksinf12(1+ecosf)[(4e2+A(4−5e2))cosf+4(3+e2−2A(1+e2))1+cosf+e(12−12A+4Acos2f+Aecos3f)]$$ \begin{equation} \begin{array}{rl} x'+\left[\frac{2e-\csc f\cot f-e\, \csc ^{2} f}{\csc f+e\cot f} \right]x=&\frac{e\, d-2d\, \, \csc f\cot f-[(1+e)\, x_{0} -2\, d]\csc ^{2} f}{\csc f+e\cot f} \\ &+\frac{k\, \sin f}{12\left(1+e\cos f\right)} \left[\left(4e^{2} +A\left(4-5e^{2} \right)\right)\cos f\right.\\ &+\frac{4\left(3+e^{2} -2A\left(1+e^{2} \right)\right)}{1+\cos f}+\left. e\left(12-12A+4A\cos 2f+Ae\cos 3f\right)\right] \end{array} \end{equation} $$(31)

Which is in the common form of the linear first order differential equation. To solve such type, we need to get the integrating factor as following

exp{∫[2e−cscfcotf−ecsc2fcscf+ecotf]df}=1sinf(1+ecosf)$$ \begin{equation} \exp \left\{{\int {\left[ {\frac{{2e - \csc f\cot f - e\,\,{{\csc}^2}f}}{{\csc f + e\cot f}}} \right]df}} \right\} = \frac{1}{{\sin f(1 + e\cos f)}}\ \end{equation} $$(32)

Thus, the solution of 31 can be written as

x(f)sinf(1+ecosf)=∫ed sin2f−2dcosf−[(1+e)x0−2d]sin2f(1+ecosf)2df+∫k12(1+ecosf)2[(4e2+A(4−5e2))cosf+4(3+e2−2A(1+e2))1+cosf+e(12−12A+4Acos2f+Aecos3f)]df+c1$$ \begin{equation} \begin{array}{rl} \frac{x(f)}{\sin f(1+e\cos f)} =&\int \frac{e\, d\sin ^{2} f-2d\cos f-[(1+e)x_{0} -2d]}{\sin ^{2} f(1+e\cos f)^{2}} df +\int \frac{k}{12\left(1+e\cos f\right)^{2}} \left[\left(4e^{2} +A\left(4-5e^{2} \right)\right)\cos f\right.\\ &+\frac{4\left(3+e^{2} -2A\left(1+e^{2} \right)\right)}{1+\cos f} +\left.e\left(12-12A+4A\cos 2f+Ae\cos 3f\right)\right] df\, +c_{1} \end{array} \end{equation} $$(33)

In order to accomplish the integral, we use the eccentric anomaly E, rather than the true anomaly,f, defined by E=2tan−1[1−e1+etanf2] $E=2\tan ^{-1}\left[\sqrt{\frac{1-e}{1+e}} \tan \frac{f}{2} \right]$ , then we can obtain x(f) as:

x(f)=xo2(1+e)(1+cosf)(1+ecosf)+e2(1−e)2[(d+de−exo)(1+e)+k(1+e−2Ae)3]sin2f+sinf(1+ecosf){3e(exo−d−de)2(1−e)5/2(1+e)3/2E+(4d−(1+e)xo)2(1−e)2tanf2+13k[Asinf+3e(e(2A−1)−1)(1−e)21−e2E+(3+e2−2A(1+e2))(1−e)2tanf2]+c1}$$ \begin{equation} \begin{array}{rl} x(f)=&\frac{x_{o}}{2\left(1+e\right)} \left(1+\cos f\right)\left(1+e\cos f\right)+\frac{e^{2}}{\left(1-e\right)^{2}} \left[\frac{\left(d+d\, e-ex_{o} \right)}{\left(1+e\right)} +\frac{k\, \, \left(1+e-2Ae\right)}{3} \right]\sin ^{2} f \\ &+\sin f\left(1+e\cos f\right)\left\{\frac{3e\left(e\, x_{o} -d-de\right)\,}{2\left(1-e\right)^{5/2} \left(1+e\right)^{3/2}} E\right. +\frac{\left(4d-\left(1+e\right)x_{o} \right)}{2\left(1-e\right)^{2}} \tan \frac{f}{2} \\ &+\frac{1}{3} k\left[A\sin f+\frac{3e\left(e\left(2A-1\right)-1\right)}{\left(1-e\right)^{2} \sqrt{1-e^{2}}} E+\frac{\left(3+e^{2} -2A\left(1+e^{2} \right)\right)}{\left(1-e\right)^{2}} \right. \left. \left. \tan \frac{f}{2} \right]+\, \, c_{1} \right\} \end{array} \end{equation} $$(34)

To get the value of c₁, we can find dxdf|f=0⇒c1=x′01+e $\left. \frac{dx}{df} \right|_{f=0} \, \, \, \Rightarrow \, \, c_{1} =\, \, \frac{x'_{0}}{1+e}$ , then

x(f)=(1+ecosf)[xo(1+cosf)2(1+e)+kA sin2f3]+e2sin2f(1−e)2[d−exo1+e+k(1+e−2Ae)3]+sinf(1+ecosf){eE(1−e)5/2(1+e)1/2[3exo2(1+e)−32d+k(e(2A−1)−1)]+12(1−e)2tanf2[(4d−(1+e)xo)+2k3(3+e2−2A(1+e2))]+x′o1+e}$$ \begin{equation} \begin{array}{rl} x(f)=&\left(1+e\cos f\right)\left[\frac{x_{o} \left(1+\cos f\right)}{2\left(1+e\right)} +\frac{kA\, \sin ^{2} f}{3} \right]+\frac{e^{2} \sin ^{2} f}{\left(1-e\right)^{2}} \left[d-\frac{e\, x_{o}}{1+e} +\frac{k\, \left(1+e-2Ae\right)}{3} \right] \\ &+\sin f\left(1+e\cos f\right)\left\{\frac{e\, E}{\left(1-e\right)^{5/2} \left(1+e\right)^{1/2}} \left[\frac{3\, e\, x_{o} \,}{2\left(1+e\right)} -\frac{3}{2} d+k\left(e\left(2A-1\right)-1\right)\right]\right. \\ &+\frac{1}{2\left(1-e\right)^{2}} \tan \frac{f}{2} \left[\left(4d-\left(1+e\right)x_{o} \right)+\frac{2k}{3} \left(3+e^{2} -2A\left(1+e^{2} \right)\right)\right]\left. +\, \, \frac{x'_{o}}{1+e} \right\} \end{array} \end{equation} $$(35)

Now to obtain the solution for y, we substitute the obtained solution of x(f) into equation 25, and integrating with respect to f, we obtain

⇒y=∫[−2x+kA(sin2f−2e3cos3f)+d]df+c2$$ \begin{equation} \Rightarrow y=\int \left[-2x+k\, A\, \, (\sin ^{2} f-\frac{2e}{3} \cos ^{3} f)+d\right] \, \, df+c_{2} \end{equation} $$(36)

To evaluate ∫ E sin f (1+e cos f)df, we again use the eccentric anomaly

I=∫Esinf(1+ecosf)df=(1−e2)2∫EsinE(1−ecosE)3dE$$ \begin{equation} I=\int E\sin f(1+e\cos f)df =(1-e^{2} )^{2} \int \frac{E\, \sin E}{(1-e\cos E)^{3}} dE \end{equation} $$(37)

Integrating 37 by parts, we get

I=12e[1−e2(f+esinf)−E(1+ecosf)2)$$ \begin{equation} I=\frac{1}{2e\,}\,\left[ \sqrt{1-{{e}^{2}}}\,\left( f+e\,\,\sin f \right)-E{{\left( 1+e\cos f \right)}^{2}}\, \right) \end{equation} $$(38)

The constant of integration c₂ in 36, can be calculated from y|f=0⇒c2=y0−2x′01+e $\left. y\right|_{f=0} \, \, \Rightarrow \, \, \, c_{2} =y_{0} -\frac{2\, x'_{0}}{1+e}$ The solution for y(f) can now be expressed as:

y(f)=x′01+e(2cosf−e sin2f)−x04(1+e)(2f(2+e)+4(1+e)sinf+esin2f)−kA18(6f+3esinf−3sin2f−esin3f)−e22(1−e)2(d−ex01+e+k(1+e−2Ae)3)(2f−sin2f)−14(1−e)2(4d−(1+e)x0+2k3(3+e2−2A(1+e2)))(2f(2−e)−4(1−e)sinf−esin2f)−1(1−e)51+e(3ex02(1+e)−3d2+k(e(2A−1)−1))(1−e2(f+esinf)−E(1+ecosf)2)−136Ak(9(−2f+sin2f)+2e(9sinf+sin3f))+df+y0−2x′01+e$$ \begin{equation} \begin{array}{rl} y(f)=&\frac{x'_{0}}{1+e} \left(2\cos f-e\sin ^{2} f\right)-\frac{x_{0}}{4\left(1+e\right)} \left(2f(2+e)+4(1+e)\sin f+e\sin 2f\right) \\ &-\frac{kA}{18} \left(6f+3e\sin f-3\sin 2f-e\sin 3f\right)-\frac{e^{2}}{2\left(1-e\right)^{2}} \left(d-\frac{e\, x_{0}}{1+e} +\frac{k\left(1+e-2Ae\right)}{3} \right)\left(2f-\sin 2f\right) \\ &-\frac{1}{4\left(1-e\right)^{2}} \left(4d-\left(1+e\right)x_{0} +\frac{2k}{3} \left(3+e^{2} -2A\left(1+e^{2} \right)\right)\right)\left(2f\left(2-e\right)-4\left(1-e\right)\sin f-e\sin 2f\right) \\ &-\frac{1}{\sqrt{\left(1-e\right)^{5}} \sqrt{1+e}} \left(\frac{3\, e\, x_{0}}{2\left(1+e\right)} -\frac{3d}{2} +k\left(e\left(2A-1\right)-1\right)\right)\left(\sqrt{1-e^{2}} \left(f+e\sin f\right)-E\left(1+e\cos f\right)^{2} \right) \\ &-\frac{1}{36} Ak\left(9\left(-2f+\sin 2f\right)+2e\left(9\sin f+\sin 3f\right)\right)+d\, f+y_{0} -\frac{2\, x'_{0}}{1+e} \end{array} \end{equation} $$(39)

Finally for the third equation in 23,

z″+z=k(1+ecosf)sin2isinf$$ [z''+z=k\left(1+e\, \cos f\right)\, \sin 2i\, \, \sin f $$

Applying Laplace transformtion such that ℒ{z(f)}=Z(s) ${\rm \mathscr {L}}\{z(f)\} =Z(s)$ , we can get

Z(s)=ksin(2i)[1(s2+1)2+e(s2+1)(s2+4)]+z0ss2+1+z′0s2+1$$ \begin{equation} Z(s)=k\sin \left(2i\right)\, \left[\frac{1}{\left(s^{2} +1\right)^{2}} +\frac{e}{\left(s^{2} +1\right)\left(s^{2} +4\right)} \right]+\frac{z_{0} \, s}{s^{2} +1} +\frac{z'_{0}}{s^{2} +1} \end{equation} $$(40)

The inverse Laplace transform of equation 40, gives the solution of z(f) as:

z=z′0sinf+z0cosf+k6(sin2i)(3sinf−3fcosf+2esinf−esin2f)$$ \begin{equation} z=z'_{0} \sin f+z_{0} \cos f+\frac{k}{6} \left(\sin 2i\right)\left(3\sin f-3f\cos f+2e\sin f-e\sin 2f\right) \end{equation} $$(41)

Equations 35, 39 and 41 are the solution of the equations of motion of the deputy relative to the chief including perturbation due to oblateness of the Earth.

In this section, we will present a simulation of the relative equations of motion of the deputy with respect to the chief satellite in in both the unperturbed and the perturbed cases and we will show the difference between them. Where the unperturbed solution can be written as:

xunp(f)=xo2(1+e)(1+cosf)(1+ecosf)+e2(1−e)2(d−exo1+e)sin2f+sinf(1+ecosf){3e(exo−d−de)2(1−e)5/2(1+e)3/2E+(4d−(1+e)xo)2(1−e)2tanf2+x′o1+e}$$ \begin{equation} \begin{array}{rl} x_{unp} (f)=&\frac{x_{o}}{2\left(1+e\right)} \left(1+\cos f\right)\left(1+e\cos f\right)+\frac{e^{2}}{\left(1-e\right)^{2}} \left(d-\frac{e\, x_{o}}{1+e} \right)\sin ^{2} f \\ &+\sin f\left(1+e\cos f\right)\left\{\frac{3\, e\, \left(e\, x_{o} -d-d\, e\right)}{2\left(1-e\right)^{5/2} \left(1+e\right)^{3/2}} \right. E+\frac{\left(4d-\left(1+e\right)x_{o} \right)}{2\left(1-e\right)^{2}} \tan \frac{f}{2} \left. +\, \, \frac{x'_{o}}{1+e} \right\} \end{array} \end{equation} $$(42)yunp(f)=x′01+e(2cosf−e sin2f)−x04(1+e)(2f(2+e)+4(1+e)sinf+esin2f)−(d+de−ex0)e22(1−e)2(2f−sin2f)−(4d−(1+e)x0)4(1−e)2(2f(2−e)−4(1−e)sinf−esin2f)−3(ex0−d−de)2(1−e)5/21+e(1−e2(f+esinf)−E(1+ecosf)2)+df+y0−2x′01+e$$ \eqalign{ y_{unp} (f)= & \frac{x'_{0} }{1+e} \left(2\cos f-e\sin ^{2} f\right)-\frac{x_{0} }{4\left(1+e\right)} \left(2f(2+e)+4(1+e)\sin f+e\sin 2f\right) \\ & -\frac{\left(d+d\, e-e\, x_{0} \right)e^{2} }{2\left(1-e\right)^{2} } \left(2f-\sin 2f\right)-\frac{\left(4d-\left(1+e\right)x_{0} \right)}{4\left(1-e\right)^{2} } \left(2f\left(2-e\right)-4\left(1-e\right)\sin f-e\sin 2f\right) \\ & -\frac{3\left(e\, x_{0} -d-d\, e\right)}{2\left(1-e\right)^{{5/2} } \sqrt{1+e} } \left(\sqrt{1-e^{2} } \left(f+e\sin f\right)-E\left(1+e\cos f\right)^{2} \right)+d\, f+y_{0} -\frac{2\, x'_{0} }{1+e}} $$(43)zunp=z′0sinf+z0cosf$$ \begin{equation} {{z}_{unp}}={{{z}'}_{0}}\sin f+{{z}_{0}}\cos f \end{equation} $$(44)

According to the following initial conditions [11] e = 0.1, x₀ = 0.1, x′₀, y₀ =0, y′0=−21110 $ y'_{0} =\frac{-21}{110}$ , z₀ = 0.08 andz′₀ = 0, a = 7000km we can get the relations between the chief true anomaly with the x, y and z components of the deputy relative position in both perturbed and unperturbed cases as shown in figures 2, 3 and 4 respectively. Also in figure 5 these three components are parametrized with the chief true anomaly.

Fig. 2

Fig. 3

Fig. 4

Fig. 5

An explicit solution including the perturbation due to the oblateness of the Earth of the relative equations of motion of a deputy or follower satellite relative to a chief or leader satellite is expressed in terms of the eccentricity of the chief orbit and it’s true anomaly as the independent variable. We have the following notes: