Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 489849, 18 pages
Research Article

Analytical and Semi-Analytical Treatment of the Satellite Motion in a Resisting Medium

1Department of Applied Mathematics, Faculty of Applied Science, Taibah University, Al-Madina-Al Munawwarah, Saudi Arabia
2Department of Mathematics, Faculty of Science, Tanta University, Tanta, Egypt
3Department of Mathematics, Faculty of Science, Taibah University, Al-Madina-Al Munwwarah, Saudi Arabia
4Department of Astronomy, Faculty of Science, Cairo University, Cairo 12613, Egypt

Received 11 August 2011; Revised 8 October 2011; Accepted 17 October 2011

Academic Editor: Pablo González-Vera

Copyright © 2012 S. E. Abd El-Bar and F. A. Abd El-Salam. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


The orbital dynamics of an artificial satellite in the Earth's atmosphere is considered. An analytic first-order atmospheric drag theory is developed using Lagrange's planetary equations. The short periodic perturbations due to the geopotential of all orbital elements are evaluated. And to construct a second-order analytical theory, the equations of motion become very complicated to be integrated analytically; thus we are forced to integrate them numerically using the method of Runge-Kutta of fourth order. The validity of the theory is checked on the already decayed Indian satellite ROHINI where its data are available.

1. Introduction

The problem of the orbital motion of celestial bodies in a resisting medium can be traced back to Newton's Principle [1] and Laplace's analysis of a possibly finite velocity of gravitation [2]. However, with the advent of the space age there arose the need of precisely predicting the orbital behavior of the newly launched bodies.

When the artificial satellites move in Low Earth Orbits (in brief LEO) their orbital motion is perturbed by the resisting force. Therefore, an adequate analytical description of the drag force is required to predict precisely the orbital perturbations due to atmospheric drag. These perturbations due to the atmospheric drag are more complicated and difficult to deal with than orbits perturbed by gravity because of the nonconservative property of the drag force. The adequate allowance for atmospheric drag turned out to be a problem for unforeseen complexity [38]. In order to be able to solve the equations of motion including the drag force analytically, one is forced to do many approximations. On the other hand, from the computational point of view, having an analytical solution at hand, we can jump from the initial conditions to the new state, which makes the computation extremely fast.

The order of magnitude of the atmospheric drag acting on a satellite depends on the altitude of the satellite. The orbital behavior of an LEO satellite could be strongly influenced by atmospheric drag [4, 9, 10], and therefore affects the quality of the remote sensing of the satellite [11, 12]. Traditionally, the problem of orbit disturbance by atmospheric drag is solved by numerical integration [1318]. An analytical solution gives the theoretical integrals and shows the physical effects with very clear spectral properties [1924]. This may give a direct insight into the physical phenomenon of the disturbance [2527]. Bezdeka and Vokrouhlicky [28] developed a semianalytic theory of motion for close-Earth spherical satellites including drag and gravitational perturbations. Xu et al. [29] derived an analytical solution of a satellite orbit disturbed by atmospheric drag. They first transformed the disturbance force vector and rotated it to the orbital frame so that it can be used in the simplified Gaussian equations of satellite motion. The disturbances are separated into three parts: shortperiodic terms with triangular functions of the mean anomaly𝑀, longperiodic terms with triangular functions of the argument of perigee and inclination (𝜔,𝐼), and secular terms [nonperiodic functions of semimajor axis and eccentricity (𝑎,𝑒)].

The main effects of the force are to bring out large secular decreases in the elements 𝑎, the semimajor axis, and 𝑒, the eccentricity, that cause the satellite plunges toward Earth. This perturbing influence is one of the most important perturbing forces in the altitude regime from 150 to 600 km which is called the drag force. The drag force,𝐹𝐷, per unit mass of the satellite can be represented by𝐹𝐷1=2𝐶𝐷𝐴𝑚𝜌𝑉𝑉,(1.1) where 𝐶𝐷 is a dimensionless drag coefficient, empirically determined to have a value close to 2; it can be taken as 2.2 with an error (standard deviation) which should not exceed 5% unless the assumption about the accommodation coefficient is greatly in error. But in fact the drag coefficient depends on the mean free path of the atmospheric molecules with respect to the linear dimensions of the satellite; therefore, the orbital theory should be applied with caution for large satellites in their last revolutions [4]. After that Sehnal [30] derived a theoretical expression for the dimensionless drag coefficient (see (3.1) there). 𝑉 is the velocity of the satellite relative to the atmosphere (called the ambient velocity), 𝐴 represents the effective cross-sectional area of the satellite, which is to be found by averaging all possible projected areas of the satellite onto a plane perpendicular to 𝑉, 𝑚 is the mass of the satellite, and 𝜌 is the density of the atmosphere. The density 𝜌 depends on position and time in a very complicated manner. Its variation with height depends on the form of the atmosphere, since it decreases rapidly with slight growing in the altitude; the effect of the flattening of the atmosphere is rather significant.

2. Rotation of the Atmosphere

This effect is so small, so a simple approximation is justified. The differential rotation of the atmosphere complicates the numerical calculation and it gives no sense in this very small effect. So we will assume that the atmosphere rotates with the same angular velocity of the Earth 𝜎, then𝑉=𝑣𝜎×𝑟,(2.1) where 𝑣 is the orbital velocity of the satellite, and 𝑟 its radius vector. Thus𝑉2=𝑣22𝑣𝜎×𝑟+𝜎×𝑟2.(2.2)

If (̂̂̂𝑖,𝑗,𝑘) are the base vectors of the rectangular geocentric equatorial frame associated with the spherical coordinates (𝑟,𝜆,𝛿), then ̂𝑘𝜎=𝜎 with 𝜎=7.27205217×105rad/s, and𝑉2=𝑣22𝜎cos𝐼+𝜎2𝑟2cos2𝛿,(2.3) where 𝐼 is the orbital inclination,𝛿 is the satellite's declination, and is the integral of areas. Now the ambient velocity can be written in the form:𝑉2=𝑣2𝑘𝑅,(2.4)

where𝑘𝑅=12𝜎cos𝐼𝑣2+𝜎2𝑟2cos2𝛿𝑣2.(2.5)The third term in 𝑘𝑅 may usually be neglected, being less than about 1/250, and the second term which is of order of 1/15 may be evaluated at perigee, where the density is much greater than elsewhere on the orbit; thus we may effectively regard 𝑘𝑅 as a constant.

3. The Atmospheric Drag Force Model

We have assumed that the atmospheric density is approximated as an exponential function which is given by𝜌=𝜌𝑝𝑜𝛽𝑟exp𝑝𝑜𝑟,(3.1) where 𝜌𝑝𝑜 is the atmospheric density at the initial perigee point 𝑟𝑝𝑜, 𝛽 is taken as a constant, 1/𝛽 is about 50 kilometers for most satellite orbits, and 𝑟 is given by relation: 𝑟=𝑎(1𝑒cos𝐸),(3.2) where 𝑎, 𝑒 are the orbital semimajor axis and eccentricity and 𝐸 is the eccentric anomaly of the satellite. Now𝜌=𝜌𝑝𝑜𝛽𝑎exp𝑜𝑎𝑥𝑜+𝛽𝑥cos𝐸,(3.3)

where 𝑥=𝑎𝑒,𝑥𝑜=𝑎𝑜𝑒𝑜,𝑉 is the velocity of satellite relative to ambient atmosphere and is generally through𝑉=𝑣𝜔𝑎×𝑟,(3.4) where 𝑣 is the velocity of satellite with respect to the Earth's center, 𝜔𝑎 is the west-to-east angular velocity of the atmosphere and 𝑟is the position vector of the satellite.

The radial 𝑉𝑠, transverse 𝑉𝑡, and orthogonal 𝑉𝑤 components of the velocity are given by𝑉𝑠=𝜇𝑎1𝑒2𝑉𝑒sin𝑓,𝑡=𝜇𝑎1𝑒2(1+𝑒cos𝑓)𝑟𝜔𝑎𝑉cos𝐼,𝑤=𝜔𝑎𝑟sin𝐼cos(𝑓+𝜔),(3.5)

where 𝜇=𝐺(𝑀+𝑚)𝐺𝑀, the Earth's gravitational parameter, is the product of the universal constant of gravity 𝐺 and the total mass of the Earth 𝑀, the mass of the satellite 𝑚, and 𝑓,𝜔 are the true anomaly and the argument of the perigee of the satellite. Then, the radial 𝑆𝐷, transverse 𝑇𝐷, and orthogonal 𝑊𝐷 components of the drag force are given by𝑆𝐷1=2𝐶𝐷𝐴𝑚𝜌𝑉𝜇𝑎1𝑒2𝑇𝑒sin𝑓,𝐷1=2𝐶𝐷𝐴𝑚𝜌𝑉𝜇𝑎1𝑒2(1+𝑒cos𝑓)𝑟𝜔𝑎,𝑊cos𝐼𝐷1=2𝐶𝐷𝐴𝑚𝜌𝑉𝑟𝜔𝑎sin𝐼cos(𝑓+𝜔),(3.6)


Using the relationscos𝑓=cos𝐸𝑒,1𝑒cos𝐸(3.8)sin𝑓=1𝑒2sin𝐸.1𝑒cos𝐸(3.9)

Equations (3.6), can be written in the form𝑆𝐷𝑛=𝑏1cos𝐸𝑎𝑒𝜌𝑉𝑇sin𝐸,𝐷𝑛=𝑏1𝑒21cos𝐸𝑎𝜌𝑉1𝑑(1𝑒𝑑cos𝐸)21𝑒2,𝑊𝐷=𝑏𝜔𝑎𝑎𝜌𝑉(1𝑒cos𝐸)sin𝐼cos(𝑓+𝜔),(3.10)


where 𝑛=𝜇/𝑎3 is the mean motion.

For a close Earth’s satellite the ratio 𝜔𝑎/𝑛 is about 1/15; the third term in the expansion is thus only about 1/500 of the leading term and will be ignored. For satellites that are not close, the ratio will be larger but the atmospheric density that will multiply this expression will then be so small that the product will still be unimportant. we thus take𝑉=𝜇𝑎1+𝑒cos𝐸1𝑒cos𝐸1𝑑1+𝑒cos𝐸,11𝑒cos𝐸𝑏=2𝐶𝐷𝐴𝑚.(3.12)

4. Lagrange's Planetary Equations

The orbit is uniquely determined by 6 elements, namely, the orbital elements. The Lagrange planetary equations are a set of 6 differential equations which measure the rate change of these orbital elements 𝑎,𝑒,𝑖,𝜔,Ω and any time element, for example, the mean anomaly 𝑀̇𝑎=2𝑒sin𝑓𝑛1𝑒2𝑆+2𝑎1𝑒2𝑛𝑟𝑇,̇𝑒=1𝑒2sin𝑓𝑛𝑎𝑆+1𝑒2𝑎𝑛𝑎𝑒1𝑒2𝑟𝑟𝑎𝑇,𝐼=𝑟cos(𝑓+𝜔)𝑛𝑎21𝑒2̇𝑊,Ω=𝑟sin(𝑓+𝜔)𝑛𝑎21𝑒2̇sin𝐼𝑊,̇𝜔=Ωcos𝐼1𝑒2𝑟𝑛𝑎𝑒𝑆cos𝑓𝑇1+𝑎1𝑒2,̇2sin𝑓𝑀=𝑛𝑛𝑎2𝑆𝑟̇𝜔1𝑒2̇Ω1𝑒2cos𝐼.(4.1)

There are many forms of the Lagrange planetary equations. In many cases, some of orbital elements are undefined, and hence, we have to replace this set by a new set of elements valid to define the orbit, and the planetary equations will take a new form in terms of the new element. For the equatorial orbits 𝐼=0 or 𝐼=180o or even when 𝐼, 𝑒 have small values we have near equatorial near circular orbits and the periapsis is undefined and hence the time of pericentric passage,𝑇 and 𝜔 are undefined or at least meaningless and should be replaced by other elements as the longitude of periapsis, 𝜛, defined as 𝜛=𝜔+Ω.

5. The Earth's Gravity Force Model

The gravitational field of the Earth is quite irregular. However, it can be represented by an infinite series of spherical harmonics, choosing the principal axes with origin of coordinates at centre of mass in addition assuming also that the gravitational field of the Earth is axially symmetric; then the form adopted is, [31]𝑈𝜇=𝑟1𝑛=2𝑅𝑟𝑛𝐽𝑛𝑃𝑛(sin𝛿),(5.1) where 𝑅 is the equatorial radius of the Earth, 𝜇 is the Earth’s gravitational parameter, (𝑟,𝛿) are the appeared geocentric coordinates of the satellite, 𝐽𝑛are zonal harmonic coefficients, and 𝑃𝑛(sin𝛿) the Legendre Polynomials are given by𝑃𝑛1(sin𝛿)=2𝑛𝑁𝑗=0(1)𝑗(2𝑛2)!(sin𝛿)𝑛2𝑗𝑗!(𝑛𝑗)!(𝑛2𝑗)!,(𝑛=0,1,2,3,),(5.2) where 𝑛𝑁=2𝑛even,𝑛12𝑛odd.(5.3)

After some lengthy straightforward algebraic calculation we can evaluate the force components in the radial 𝑆𝐺, transverse 𝑇𝐺, and normal 𝑊𝐺 components as 𝑆𝐺3=4𝐽2𝜇𝑅2𝑟423sin2𝐼(1cos2(𝑓+𝜔))+𝐽3𝑅3𝜇2𝑟5sin𝐼15sin2𝐼12sin(𝑓+𝜔)5sin2𝐼sin3(𝑓+𝜔)+𝐽45𝑅4𝜇64𝑟624120sin2𝐼+105sin4𝐼+206sin2𝐼7sin4𝐼×cos2(𝑓+𝜔)+35sin4,𝑇𝐼cos4(𝑓+𝜔)𝐺=3𝐽2𝜇𝑅2𝑟4sin𝐼cos𝐼sin𝐹1,1cos𝐹1,1𝐽33𝑅3𝜇4𝑟5sin𝐼cos𝐹1,12+5sin2𝐼sin2𝐼cos𝐹2,2𝐽45𝑅4𝜇8𝑟6sin2𝐼cos𝐹1,112+21sin2𝐼sin𝐹1,17sin2𝐼sin𝐹3,3,𝑊𝐺=3𝐽2𝜇𝑅2𝑟4sin𝐼cos𝐼sin(𝑓+𝜔)𝐽33𝑅3𝜇4𝑟5cos𝐼2+5sin2𝐼sin2𝐼cos2(𝑓+𝜔)𝐽45𝑅4𝜇8𝑟6sin𝐼cos𝐼12+21sin2𝐼sin(𝑓+𝜔)7sin2.𝐼sin3(𝑓+𝜔)(5.4)

6. The Analytical Solution

In this section, we proceed to calculate the first-order approximation of the orbital elements caused by the drag and the gravity forces during revolution by integrating Lagrange's planetary equations (4.1) from 𝐸=𝜋 to 𝐸=𝜋, and retain only terms to order 𝑂(𝑒2), after substitution of {𝑆=𝑆𝐷+𝑆𝐺,𝑇=𝑇𝐷+𝑇𝐺,𝑊=𝑊𝐷+𝑊𝐺} into the Lagrange planetary equations we have the perturbations in the elements asΔ𝜎=Δ𝜎𝐷+Δ𝜎𝐺,(6.1) where 𝜎(𝑎,𝑒,𝐼,𝜔,Ω,𝑀) denotes any of the orbital elements, Δ𝜎𝐷 the perturbations due to the drag force, and Δ𝜎𝐺 the perturbations due to the gravity force is given byΔ𝑎𝐷=4𝜋𝑎2𝐻𝑜𝜔12𝑎𝑛+3cos𝐼2𝜔1+𝑎𝑛𝑒cos𝐼2𝐼𝑜+2𝑒𝐼1+1232+𝜔𝑎𝑛𝐼cos𝐼2𝑒2,Δ𝑎𝐺=𝐽252𝑖=1𝑗=2𝛼𝑖𝑗cos𝐹𝑖,𝑗+𝐽373𝑘=1𝑙=3𝛽𝑘𝑙sin𝐹𝑘,𝑙+𝐽494𝑚=1𝑛=4𝛾𝑚𝑛cos𝐹𝑚,𝑛,Δ𝑒𝐷=4𝜋𝑎𝑏𝐻0𝑒21+5𝜔𝑎𝐼2𝑛cos𝐼0+12𝜔𝑎𝑛cos𝐼𝑒2582𝜔𝑎𝑛𝐼cos𝐼1+𝑒21+3𝜔𝑎𝐼2𝑛cos𝐼2+𝑒28𝐼3,Δ𝑒𝐺=𝐽252𝑖=1𝑗=2𝛿𝑖𝑗cos𝐹𝑖,𝑗+𝐽3𝑋1𝑓+73𝑚=1𝑛=3𝜎𝑚𝑛sin𝐹𝑚,𝑛+𝐽4𝑋2𝑓+94𝑘=1𝑙=4𝑄𝑘𝑙cos𝐹𝑘,𝑙,Δ𝐼𝐷=𝜋𝑎𝑏𝜔𝑎𝐻𝑜sin𝐼𝑛𝜔1𝑎𝑛+𝑒cos𝐼22943𝜔𝑎𝑛cos𝐼+178+cos2𝜔55𝜔𝑎𝐼8𝑛cos𝐼cos2𝜔𝑜+𝑒1+3𝜔𝑎𝑛1cos𝐼+25𝜔𝑎𝑛𝐼cos𝐼9cos2𝜔1+𝜔1𝑎𝑛+𝑒cos𝐼22517𝜔𝑎𝑛+cos𝐼121𝜔𝑎𝑛𝐼cos𝐼cos2𝜔2+𝑒2𝜔𝑎𝑛cos𝐼1cos2𝜔𝐼3+𝑒216911𝜔𝑎𝑛cos𝐼cos2𝜔𝐼4,Δ𝐼𝐺=𝐽23𝑖=1𝑍𝑖cos𝐹𝑖,2+𝐽3𝑋3𝑓+53𝑘=1𝑙=1𝐶𝑘𝑙sin𝐹𝑘,𝑙+𝐽4𝑋4𝑓+74𝑝=1𝑞=2𝑆𝑝𝑞cos𝐹𝑝,𝑞,ΔΩ𝐷=𝑎𝑏𝜔𝑎𝐻0𝜔𝑎𝑛3𝑒4𝜔1𝑎𝑛+𝑒cos𝐼284164𝜔𝑎𝑛+𝑒cos𝐼sin2𝜔22𝜔1𝑎𝑛𝐼cos𝐼cos2𝜔𝑜+𝑒297𝜔𝑎𝑛𝑒cos𝐼22611𝜔𝑎𝑛cos𝐼sin2𝜔𝐼1+𝜔1𝑎𝑛𝜔cos𝐼+𝑒1𝑎𝑛+𝑒cos𝐼24933𝜔𝑎𝑛𝑒cos𝐼sin2𝜔+22𝜔1𝑎𝑛𝐼cos𝐼cos2𝜔2+𝑒23𝜔𝑎𝑛𝑒cos𝐼12449𝜔𝑎𝑛cos𝐼×sin2𝜔𝐼3+𝑒822𝜔𝑎𝑛+cos𝐼513𝜔𝑎𝑛cos𝐼sin2𝜔𝐼4+𝑒24𝜔𝑎𝑛cos𝐼sin2𝜔𝐼5,ΔΩ𝐺=𝐽2𝑋8𝑓+32𝑏=1𝑐=0𝜌𝑏𝑐sin𝐹𝑏,𝑐+𝐽3𝑋9𝑓+53𝛼=1𝛽=1𝑞𝛼𝛽cos𝐹𝛼,𝛽+𝐽4𝑋10𝑓+74𝑖=1𝑗=2𝑝𝑖𝑗sin𝐹𝑖,𝑗,Δ𝜔𝐷=ΔΩ𝐷cos𝐼,Δ𝜔𝐺=𝐽2𝑋5𝑓+52𝑝=1𝑞=2𝑈𝑝𝑞sin𝐹𝑝,𝑞+𝐽3𝑋6𝑓+73𝑘=1𝑙=3𝑉𝑘𝑙cos𝐹𝑘,𝑙+𝐽4𝑋7𝑓+94𝑚=1𝑛=4𝑍𝑚𝑛sin𝐹𝑚,𝑛,Δ𝑀=2𝜋𝑛,(6.2)


is the first kind modified Bessel function of order 𝑛, and by using the values of 𝐼0(𝑐) and 𝐼1(𝑐) one finds that𝐼2=𝐼02𝑐𝐼1,𝐼34=𝑐𝐼0+81+𝑐2𝐼1,𝐼4=1+24𝑐2𝐼061𝑐2𝐼1,𝐼5=12𝑐1+16𝑐2𝐼0+1+72𝑐2+384𝑐4𝐼1,(6.4)

where 𝐻𝑜 is the scale height. It is given by𝐻𝑜=𝜌𝑝𝑜𝛽𝑎exe𝑜𝑎𝑥𝑜.(6.5)

Finally the nonvanishing coefficients included in the gravity perturbations are previously computed by the authors [32].

7. The Numerical Solution

It is well known that the numerical solution methods provide accurate ephemeris of a satellite with respect to a force of a complex mathematical formulation. So in this section we will solve numerically the equations of motion after substituting the drag force components 𝑆𝐷, 𝑇𝐷, 𝑊𝐷, given by (3.10) together with the Earth's gravity force components 𝑆𝐺, 𝑇𝐺, 𝑊𝐺, given by (5.4) (we will use Earth's oblate model and retaining terms up to fourth zonal harmonics 𝐽4).

There are many methods to get the numerical solutions of any first-order differential equations. We will use one of the best methods of the numerical integration schemes, namely, Runge-Kutta fourth-order method.

If the Earth was assumed spherical, had no atmosphere, and the satellite is isolated from the other bodies in the solar system, the orbit of a satellite would be an ellipse of constant size and shape in a plane whose direction remains fixed relative to the stars. The orbital elements 𝑎,𝑒,𝑖,𝜔,Ω would remain unchanged. This simple thought rather insipid situation is upset by the effects of the considered perturbing forces, the Earth's oblateness, and the atmospheric drag.

The Earth's oblateness leads to two major gravitational perturbations on the satellite orbits, the first is the rotation of the orbital plane about the Earth's axis in the direction opposite to the satellite's motion so that (if𝐼<90o) the angle steadily decreases while 𝐼 remains constant. The second effect is the rotation of the major axis in the orbital plane so that the argument of perigee 𝜔 increases. The chief effects of the atmospheric drag on a satellite orbit are fortunately quite different from those of the gravitational field. Since the atmospheric density decreases rapidly as height above the Earth increases, a satellite in an orbit of appreciable eccentricity is affected mostly by drag within a small section of the orbit where it is closest to the Earth. To first-order approximation, the effect of atmospheric drag is to retard the satellite as it passes perigee, with the result that it does not swing out so far from the Earth at the subsequent apogee. The apogee height is reduced while the perigee height remains almost constant. The orbit contracts and becomes more nearly circular; the orbital elements (𝑎,𝑒) decrease steadily [4].

Now, a computer program of fourth-order Runge-Kutta method is applied for the solution of this system using the initial values and a fixed step size. Also, we will choose one of the numerical examples to show the decay in two orbital elements, namely, (𝑎,𝑒). Also we will give a numerical example to the analytical solution given by (6.1) in order to compare the results of the numerical and analytical solutions.

8. Numerical Example

In this paragraph, we give a numerical example concerning the artificial satellite ROHINI (1980 62A) [33]. The data are as follows.

The satellite mass and its cross-sectional areas are, respectively, (𝑚)=35.443kg, (𝐴)=0.319019×106. The initial orbital elements (semimajor axis, eccentricity, and the inclination) of the satellite are 𝑎0=6989.2057km, 𝑒0=0.04367712, 𝐼0=44.67198o𝜔0=239.3378okm, Ω0=174.1602o, and 𝑀0=25.63974o. And one orbital revolution was elapsed in 1.58835208 hours. The required physical and dynamical constants are𝑅Å=6378.135km,𝜇=398600.8km3/sec2,𝐽2=1.08263×103,𝐽3=2.53648×106,𝐽4=1.6233×106,𝜔𝑎=7.27205217×106𝜌rad/sec,0(300km)=0.0221kgkm3,𝜌0(250km)=0.0681kgkm3,𝜌0(200km)=0.271kgkm3,𝐶𝐷=2.2.(8.1)

The results of the numerical example of the analytical solution as well as the numerical solution are represented on Figures 1(a), 2(a), 3(a), 4(a), 5(a), 6(a), 7(a), 8(a), and 9(a) to show the decay of the semimajor axis, and Figures 1(b)– 9(b) to show the change in the eccentricity. To read the figures properly we used the letters B, C, and D to denote to the numerical example of the analytical solutions while we used the letters E, F, and G to denote for the numerical solutions. The description of these letters is as follows. B: with drag and rotation,C: with drag without rotation,D: without drag,E: with drag and rotation,F: with drag without rotation, andG: without drag.

Figure 1: (a) For one hundred revolutions at 300 km height above the Earth's surface. (b) For one hundred revolutions at 300 km height above the Earth's surface.
Figure 2: (a) For five hundred revolutions at 300 km height above the Earth's surface. (b) For five hundred revolutions at 300 km height above the Earth's surface.
Figure 3: (a) For one thousand revolutions at 300 km height above the Earth's surface. (b) For one thousand revolutions at 300 km height above the Earth's surface.
Figure 4: (a) For one hundred revolutions at 250 km height above the Earth's surface. (b) For one hundred revolutions at 250 km height above the Earth's surface.
Figure 5: (a) For five hundred revolutions at 250 km height above the Earth's surface. (b) For five hundred revolutions at 250 km height above the Earth's surface.
Figure 6: (a) For one thousand revolutions at 250 km height above the Earth's surface. (b) For one thousand revolutions at 250 km height above the Earth's surface.
Figure 7: (a) For one hundred revolutions at 200 km height above the Earth's surface. (b) For one hundred revolutions at 200 km height above the Earth's surface.
Figure 8: (a) For five hundred revolutions at 200 km height above the Earth's surface. (b) For five hundred revolutions at 200 km height above the Earth's surface.
Figure 9: (a) For five thousand revolutions at 200 km height above the Earth's surface. (b) For five thousand revolutions at 200 km height above the Earth's surface.

All figures show the decrease in the semimajor axis and in the eccentricity for the rotating and nonrotating atmosphere at different altitudes.

9. Analysis of the Solutions

The following dynamical effects are very clear. The decay of the semimajor axis is shown on Figures 1(a) to  9(a). We depicted the analytical and numerical solutions on one figure at different altitudes. Similar situation is done for the eccentricity on Figures 1(b)9(b).(1)The decrease in the semimajor and in the eccentricity is significant (much bigger) at the high altitudes at the beginning of the motion but it oscillates afterwards, that is, the slope of the rate of change is sometimes high, sometimes moderate, and sometimes low. Therefore, we can conclude that the perturbation on the long range is nonlinear (see Figures 10(a) and 10(b)). (2)The effect of the atmospheric rotation is so small.(a) The atmospheric rotation appeared clearly on the figures of numerical solution because we have used the well-known fourth-order Runge-Kutta method.(b) The atmospheric rotation disappeared on the figures of the analytical solution since the approximation made in the analytical solution is much bigger than this effect.(3)Disregarding the drag force makes the rate of change of the semimajor axis and the rate of change of the eccentricity very little as seen clearly from the slopes of the figures of the solutions with drag and those without drag. This justifies the fact that at Low Earth Orbits (in brief LEO) the dominant secular perturbation is due to the drag force. (4)The difference between the solutions with drag and those without drag is sometimes big, sometimes low, and sometimes moderate. This signifies that the perturbation is nonlinear.

Figure 10: (a) Long range perturbation in the semimajor axis. (b) Long range perturbation in the eccentricity.

10. Conclusion

In this paper, we studied the effect of the atmospheric drag and the asphericity of the Earth on the orbital dynamics of an artificial satellite in a low Earth orbit. We approximated the density of the atmosphere to an exponential function including the rotation of the atmosphere. We obtained a first-order analytic atmospheric drag perturbations. We evaluated the short periodic perturbations of all orbital elements due to the geopotential. To obtain a more precise treatment, we integrated the equations of motion numerically using the Runge-Kutta method of fourth order. The validity of the theory is checked using the already decayed Indian satellite ROHINI where its data are available. A good agreement between the observed line elements of the satellite and our predicted solutions encourage us in future to apply our theory on any similar future missions.


The authors are deeply indebted and thankful to the deanship of the scientific research and his helpful and distinct team of employees at Taibah University, Al-Madinah Al-Munawwarah, Saudia Arabia. This research work was supported by a Grant no. (617/1431).


  1. I. Newton, Philosophiae Naturalis Principia Mathematica, Book 2, Section 4, University of California Press, Berkeley, Calif, USA, 1934, English translation by F. Cajori, Newton.s Principia.
  2. P. S. Laplace, Traité de Mécanique Céleste, Tom 4, Par 2 Courcier, 1805.
  3. D. Brouwer, “Review of celestial mechanics,” Annual Review of Astronomy and Astrophysics, vol. 1, pp. 219–234, 1963. View at Google Scholar
  4. D. G. King-Hele, Theory of Satellite Orbits in an Atmosphere, Butterworth, London, UK, 1964.
  5. D. G. King-Hele, A Tapestry of Orbits, Cambridge University Press, Cambridge, UK, 1992.
  6. D. Brouwer and G. M. Clemence, Methods of Celestial Mechanics, Academic Press, New York, NY, USA, 1961.
  7. V. A. Chobotov, Ed., Orbital Mechanics, AIAA, Washington, DC, USA, 1991.
  8. D. Boccaletti and G. Pucacco, Theory of Orbits. Vol. 1, Springer, Berlin, Germany, 1996. View at Zentralblatt MATH
  9. P. D. Van Kamp, Principles of Astronomy, Freemann & Co., San Francisco, Calif, USA, 1967.
  10. R. H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA, Reston, Va, usa, 1999. View at Zentralblatt MATH
  11. S. Herrick, Astrodynamics. Vol. 2, Van Nostrand Reinhold, London, uk, 1972. View at Zentralblatt MATH
  12. D. A. Vallado, Fundamentals of Astrodynamics and Applications, Springer, Berlin, Germany, 3rd edition, 2007. View at Zentralblatt MATH
  13. J. Meeus, Astronomical Algorithms, Willmann-Bell, Richmond, Va, USA, 1991.
  14. O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods and Applications, Springer, Heidelberg, Germany, 2000.
  15. G. Seeber, Satellite Geodesy: Foundations, Methods, and Applications, Walter de Gruyter, Berlin, Germany, 2003.
  16. G. Xu, GPS—Theory, Algorithms and Applications, Springer, Heidelberg, Germany, 2nd edition, 2007.
  17. G. Xu, Orbits, Springer, Heidelberg, Germany, 2008.
  18. G. Xu, T. H. Xu, T. K. Yeh, and W. Chen, “Analytic solution of satellite orbit disturbed by atmospheric drag,” Monthly Notices of the Royal Astronomical Society, vol. 410, no. 1, pp. 654–662, 2010. View at Google Scholar
  19. C. Cui, Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wissenschaften, Reihe C, 357, 1990.
  20. M. Schneider and C. F. Cui, Deutsche Geodätische Kommission bei der Bayerischen Akademie der Wissenschaften, Reihe A: Theoretische Geodäsie, 121, 2005.
  21. S. M. Kudryavtsev, “Long-term harmonic development of lunar ephemeris,” Astronomy & Astrophysics, vol. 471, pp. 1069–1075, 2007. View at Google Scholar
  22. J. Licandro, A. Alvarez-Candal, J. De León, N. Pinilla-Alonso, D. Lazzaro, and H. Campins, “Spectral properties of asteroids in cometary orbits,” Astronomy and Astrophysics, vol. 481, no. 3, pp. 861–877, 2008. View at Publisher · View at Google Scholar
  23. A. Pál, “An analytical solution for Kepler's problem,” Monthly Notices of the Royal Astronomical Society, vol. 396, no. 3, pp. 1737–1742, 2009. View at Publisher · View at Google Scholar
  24. D. Lynden-Bell, “Analytic orbits in any central potential,” Monthly Notices of the Royal Astronomical Society, vol. 402, no. 3, pp. 1937–1941, 2010. View at Publisher · View at Google Scholar
  25. W. Torge, Geodesy, Walter de Gruyter, Berlin, Germany, 1991.
  26. J. Desmars, S. Arlot, J.-E. Arlot, V. Lainey, and A. Vienne, “Estimating the accuracy of satellite ephemerides using the bootstrap method,” Astronomy and Astrophysics, vol. 499, no. 1, pp. 321–330, 2009. View at Publisher · View at Google Scholar
  27. J. R. Touma, S. Tremaine, and M. V. Kazandjian, “Gauss's method for secular dynamics, softened,” Monthly Notices of the Royal Astronomical Society, vol. 394, no. 2, pp. 1085–1108, 2009. View at Google Scholar
  28. A. Bezdeka and D. Vokrouhlicky, “Semianalytic theory of motion for close-Earth spherical satellites includingdragand gravitational perturbations,” Planetary and Space Science, vol. 52, pp. 1233–1249, 2004. View at Google Scholar
  29. G. Xu, X. Tianhe, W. Chen, and T.-K. Yeh, “Analytical solution of a satellite orbit disturbed by atmospheric drag,” Monthly Notices of the Royal Astronomical Society, vol. 410, no. 1, pp. 654–662, 2011. View at Publisher · View at Google Scholar
  30. L. Sehnal, “Dynamics and Astrometry of natural and artificial celestial bodies,” in Proceedings of the Conference on Astrometry and Celestial Mechanics, pp. 325–332, Poznan, Poland, 1994.
  31. P. M. Fitzpatrick, Principles Of Celestial Mechancics, Academic Press, London, UK, 1970.
  32. F.A. Abd El- Salam and S. Abd El-Bar, “Computation of the different errors in the ballistic missiles range,” Journal of the Applied Mathematics, vol. 2011, Article ID 349737, 16 pages, 2011. View at Publisher · View at Google Scholar
  33. R. K. Sharma and L. Mani, “Study of RS-1 orbital decay with KS differential equations,” Indian Journal of Pure and Applied Mathematics, vol. 6, pp. 833–842, 1985. View at Google Scholar