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

Abstract

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 [3–8]. 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 [13–18]. An analytical solution gives the theoretical integrals and shows the physical effects with very clear spectral properties [19–24]. This may give a direct insight into the physical phenomenon of the disturbance [25–27]. 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𝐹^→_𝐷=−12𝐶^𝐷𝐴𝑚𝜌‖‖‖^→𝑉‖‖‖𝑉^→,(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^→𝑣⋅^→𝜎×^→𝑟+^→𝜎×^→𝑟².(2.2)

If (𝑖,̂𝑗,̂𝑘̂) are the base vectors of the rectangular geocentric equatorial frame associated with the spherical coordinates (𝑟,𝜆,𝛿), then 𝜎=𝜎𝑘̂^→ with 𝜎=7.27205217×10⁻⁵rad/s, and𝑉²=𝑣²−2𝜎ℎcos𝐼+𝜎²𝑟²cos²𝛿,(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.4)

where𝑘_𝑅=1−2𝜎ℎcos𝐼𝑣²+𝜎²𝑟²cos²𝛿𝑣².(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−𝑒²𝑒sin𝑓,_𝑡=𝜇𝑎1−𝑒²(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𝑆_𝐷=−12𝐶^𝐷𝐴𝑚𝜌‖‖‖^→𝑉‖‖‖𝜇𝑇𝑎1−𝑒²𝑒sin𝑓,_𝐷=−12𝐶^𝐷𝐴𝑚𝜌‖‖‖^→𝑉‖‖‖⎡⎢⎢⎣𝜇𝑎1−𝑒²(1+𝑒cos𝑓)−𝑟𝜔^𝑎𝑊cos𝐼_𝐷⎤⎥⎥⎦,=−12𝐶^𝐷𝐴𝑚𝜌‖‖‖^→𝑉‖‖‖𝑟𝜔_𝑎sin𝐼cos(𝑓+𝜔),(3.6)

where‖‖‖^→𝑉‖‖‖=⎡⎢⎢⎣𝜇𝑎1−𝑒²𝑒²+𝑟₂+2𝑒cos𝑓+1𝜔^2𝑎cos²𝐼+sin²𝐼cos²(𝑓+𝜔)−2𝑟𝜔_𝑎cos𝐼𝜇𝑎1−𝑒²(1+𝑒cos𝑓)⎤⎥⎥⎦^1/2.(3.7)

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

Equations (3.6), can be written in the form𝑆_𝐷=−𝑏𝑛1−cos𝐸𝑎𝑒𝜌‖‖‖^→𝑇𝑉‖‖‖sin𝐸,_𝐷=−𝑏𝑛1−𝑒²1−cos𝐸𝑎𝜌‖‖‖^→𝑉‖‖‖1−𝑑(1−𝑒𝑑cos𝐸)²𝑊_𝐷1−𝑒²,=−𝑏𝜔_𝑎𝑎𝜌‖‖‖𝑉^→‖‖‖(1−𝑒cos𝐸)sin𝐼cos(𝑓+𝜔),(3.10)

where‖‖‖^→𝑉‖‖‖=𝜇𝑎1+𝑒cos𝐸1−𝑒cos𝐸1−𝑑11+−𝑒𝑒ccooss𝐸𝐸+2𝜔𝑛^2𝑎211−+𝑒𝑒ccooss𝐸𝐸²𝑑=𝜔_𝑎+1−𝑒²cos²𝐸sin²𝐼cos²(𝑓+𝜔)+𝑒²cos²𝐼sin²𝐼+⋯,𝑛√1−𝑒²cos𝐼,(3.11)

where 𝑛=√𝜇/𝑎³ 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𝐸𝑏=11−𝑒cos𝐸1−𝑑11−+𝑒𝑒ccooss𝐸𝐸,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𝑎√1−𝑒²𝑒̇=√𝑛𝑟𝑇,1−𝑒²sin𝑓𝑛𝑎𝑆+√1−𝑒²𝑛𝑎𝑒𝑎1−𝑒²𝑟−𝑎𝑟𝑇,_⋅𝐼=𝑟cos(𝑓+𝜔)𝑛𝑎²√1−𝑒²Ω̇=𝑟sin(𝑓+𝜔𝑊),𝜔̇=−𝑛𝑎Ω²√̇co1s−𝐼𝑒−²s√in𝐼𝑊,1−𝑒²𝑀̇=𝑛−2𝑛𝑎𝑒𝑆cos𝑓−𝑇1+𝑎1𝑟−𝑒²sin𝑓,𝑛𝑎²𝑆𝑟−𝜔̇√1−𝑒²−Ω̇√1−𝑒²cos𝐼.(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 𝐼=180^o 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𝑃_𝑛(sin𝛿)=12^𝑛𝑁_𝑗=0(−1)^𝑗(2𝑛−2)!(sin𝛿)^𝑛−2𝑗𝑗!(𝑛−𝑗)!(𝑛−2𝑗)!,(𝑛=0,1,2,3,…),(5.2) where 𝑁=⎧⎪⎨⎪⎩2𝑛𝑛−1𝑛even,2𝑛odd.(5.3)

After some lengthy straightforward algebraic calculation we can evaluate the force components in the radial 𝑆_𝐺, transverse 𝑇_𝐺, and normal 𝑊_𝐺 components as 𝑆_𝐺=−34𝐽²𝜇𝑅^2⊕𝑟⁴+𝐽₃2−3sin²𝐼(1−cos2(𝑓+𝜔))𝑅^3⊕𝜇2𝑟⁵+𝐽₄sin𝐼15sin²𝐼−12sin(𝑓+𝜔)−5sin²𝐼sin3(𝑓+𝜔)5𝑅^4⊕𝜇64𝑟⁶×2c4os−21(𝑓20+si𝜔n²)𝐼++3510s5ins₄in⁴𝐼+206sin²𝐼−7sin⁴𝐼𝑇𝐼cos4(𝑓+𝜔)_𝐺,=−3𝐽₂𝜇𝑅^2⊕𝑟⁴sin𝐼cos𝐼sin𝐹_1,1cos𝐹_1,1−𝐽₃3𝑅^3⊕𝜇4𝑟⁵−𝐽₄sin𝐼cos𝐹_1,1−2+5sin²𝐼−sin²𝐼cos𝐹_2,25𝑅^4⊕𝜇8𝑟⁶𝑊_𝐺sin²𝐼cos𝐹_1,1−12+21sin²𝐼sin𝐹_1,1−7sin²𝐼sin𝐹_3,3,=−3𝐽₂𝜇𝑅^2⊕𝑟⁴−𝐽₃sin𝐼cos𝐼sin(𝑓+𝜔)3𝑅^3⊕𝜇4𝑟⁵−𝐽₄cos𝐼−2+5sin²𝐼−sin²𝐼cos2(𝑓+𝜔)5𝑅^4⊕𝜇8𝑟⁶sin𝐼cos𝐼−12+21sin²𝐼sin(𝑓+𝜔)−7sin²𝐼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 𝑂(𝑒²), 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𝜋𝑎²𝐻_𝑜1−2𝜔^𝑎𝑛cos𝐼+321+𝜔𝑛^𝑎cos𝐼𝑒²𝐼^𝑜+2𝑒𝐼₁+1232+𝜔𝑛^𝑎cos𝐼𝐼²Δ𝑎_𝐺𝑒²,=𝐽₂₅_2𝑖=1_𝑗=−2𝛼^{′𝑖𝑗}cos𝐹_𝑖,𝑗+𝐽₃₇₃_𝑘=1_𝑙=−3𝛽^{𝑘𝑙}sin𝐹_𝑘,𝑙+𝐽₄₉₄_𝑚=1_𝑛=−4Δ𝑒_𝐷𝛾^{𝑚𝑛}cos𝐹_𝑚,𝑛,=−4𝜋𝑎𝑏𝐻₀2𝑒1+52𝜔𝑛^𝑎cos𝐼𝐼⁰+1−2𝜔^𝑎𝑛cos𝐼−𝑒²58−2𝜔𝑛^𝑎cos𝐼𝐼¹+𝑒21+32𝜔𝑛^𝑎cos𝐼𝐼²+𝑒²Δ𝑒_𝐺8𝐼³,=𝐽₂₅_2𝑖=1_𝑗=−2𝛿_𝑖𝑗cos𝐹_𝑖,𝑗+𝐽₃𝑋¹𝑓+⁷₃_𝑚=1_𝑛=−3𝜎_𝑚𝑛sin𝐹_𝑚,𝑛+𝐽⁴𝑋²𝑓+⁹₄_𝑘=1_𝑙=−4𝑄_𝑘𝑙cos𝐹_𝑘,𝑙Δ𝐼_𝐷,=−𝜋𝑎𝑏𝜔_𝑎𝐻_𝑜sin𝐼𝑛+𝑒1−𝜔𝑛^𝑎cos𝐼₂249−3𝜔𝑛^𝑎cos𝐼+187cos2𝜔+558𝜔𝑛^𝑎cos𝐼cos2𝜔𝐼^𝑜+𝑒−1+3𝜔^𝑎𝑛cos𝐼+215𝜔𝑛^𝑎cos𝐼−9cos2𝜔𝐼¹+1−𝜔^𝑎+𝑒𝑛cos𝐼₂25−17𝑛𝜔^𝑎cos𝐼+1−21𝑛𝜔^𝑎cos𝐼cos2𝜔𝐼²+𝑒2𝜔𝑛^𝑎cos𝐼−1cos2𝜔𝐼³+𝑒²Δ𝐼_𝐺169−11𝑛𝜔^𝑎cos𝐼cos2𝜔𝐼⁴,=𝐽₂₃_𝑖=1𝑍_𝑖cos𝐹_𝑖,2+𝐽₃𝑋³𝑓+⁵₃_𝑘=1_𝑙=−1𝐶_𝑘𝑙sin𝐹_𝑘,𝑙+𝐽⁴𝑋⁴𝑓+⁷_4𝑝=1_𝑞=−2𝑆_𝑝𝑞cos𝐹_𝑝,𝑞ΔΩ_𝐷,=−𝑎𝑏𝜔_𝑎𝐻₀𝜔_𝑎𝑛+𝑒34𝑒1−𝜔𝑛^𝑎cos𝐼+𝑒8²41−64𝑛𝜔^𝑎cos𝐼sin2𝜔₂21−𝜔𝑛^𝑎cos𝐼cos2𝜔𝐼^𝑜+𝑒29−7𝜔𝑛^𝑎cos𝐼−𝑒2²6−11𝑛𝜔^𝑎cos𝐼sin2𝜔𝐼¹+1−𝜔^𝑎+𝑒𝑛cos𝐼+𝑒1−𝜔𝑛^𝑎cos𝐼₂49−33𝑛𝜔^𝑎cos𝐼sin2𝜔+𝑒2²1−𝜔𝑛^𝑎cos𝐼cos2𝜔𝐼²+𝑒23𝜔𝑛^𝑎cos𝐼−1−𝑒4²4−9𝜔𝑛^𝑎cos𝐼×sin2𝜔𝐼³+𝑒82−2𝜔𝑛^𝑎cos𝐼+5−13𝑛𝜔^𝑎cos𝐼sin2𝜔𝐼⁴+𝑒²ΔΩ_𝐺4𝜔𝑛^𝑎cos𝐼sin2𝜔𝐼⁵,=𝐽₂𝑋⁸𝑓+³_2𝑏=1_𝑐=0𝜌_𝑏𝑐sin𝐹_𝑏,𝑐+𝐽³𝑋⁹𝑓+⁵_3𝛼=1_𝛽=−1𝑞_𝛼𝛽cos𝐹_𝛼,𝛽+𝐽₄𝑋¹⁰𝑓+⁷_4𝑖=1_𝑗=−2𝑝_𝑖𝑗sin𝐹_𝑖,𝑗Δ𝜔_𝐷,Δ𝜔_𝐺=ΔΩ_𝐷cos𝐼,=𝐽₂𝑋⁵𝑓+⁵_2𝑝=1_𝑞=−2𝑈_𝑝𝑞sin𝐹_𝑝,𝑞+𝐽³𝑋⁶𝑓+⁷_3𝑘=1_𝑙=−3𝑉_𝑘𝑙cos𝐹_𝑘,𝑙+𝐽₄𝑋⁷𝑓+⁹₄_𝑚=1_𝑛=−4𝑍_𝑚𝑛sin𝐹_𝑚,𝑛Δ𝑀=−2𝜋𝑛,,(6.2)

where𝐼₀(𝑥)=1𝜋^𝜋0cos𝑛𝐸exp(𝑥cos𝐸)𝑑𝐸(6.3)

is the first kind modified Bessel function of order 𝑛, and by using the values of 𝐼₀(𝑐) and 𝐼₁(𝑐) one finds that𝐼₂=𝐼₀−2𝑐𝐼¹𝐼₃,=−4𝑐𝐼⁰+1+8𝑐²𝐼¹𝐼₄,=1+24𝑐²𝐼⁰−1−6𝑐²𝐼¹𝐼₅,=−12𝑐1+1𝑐6²𝐼⁰+1+72𝑐²+384𝑐⁴𝐼¹,(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 𝐽₄).

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𝐼<90^o) 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×10⁻⁶. The initial orbital elements (semimajor axis, eccentricity, and the inclination) of the satellite are 𝑎₀=6989.2057km, 𝑒₀=0.04367712, 𝐼₀=44.67198^o𝜔₀=239.3378^okm, Ω₀=174.1602^o, and 𝑀₀=25.63974^o. And one orbital revolution was elapsed in 1.58835208 hours. The required physical and dynamical constants are𝐽₃𝑅_Å=6378.135km,𝜇=398600.8km³/sec⁻²,𝐽₂=1.08263×10⁻³,=−2.5𝜌3648×10⁻⁶,𝐽₄=−1.6233×10⁻⁶,𝜔_𝑎=7.27205217×10⁻⁶rad/sec,₀(300km)𝜌₀=0.0221kgkm⁻³,𝜌₀(250km)=0.0681kgkm⁻³,(200km)=0.271kgkm⁻³,𝐶_𝐷=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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

(a)

(b)

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.

Acknowledgments

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).