The dynamical system of a satellite motion perturbed by both atmospheric drag and gravitational attraction is nonlinear, non conservative in form and the integration of the system, in general, is analytically intractable. To predict the motion precisely a mathematical representation for these forces must be selected for integrating the resulting differential equations of motion. Some of the early studies and analytical difficulties for the coupled problem were addressed by de Nike [1]. Hoots [2] used the gravitational and atmospheric models as used by Lane [3] and arrived at an improved analytical solution. Well known and commonly used models [4-9]. The KS total-energy elements equations [10] is a very powerful method for numerical solution with respect to any type of perturbing forces as the equations are less sensitive to round-off and truncation errors in the numerical integration algorithm. Sharma worked with these KS element equations to compute very accurate short-periodic terms due to J₂, even for very high eccentricity orbits [11,12]. Sharma [13,14] expanded analytic solutions by series expansion method using analytical models for oblate exponential atmospheric density model, and a model of the same with the effect of diurnal bulge.

In this paper my attempt is to get analytical solutions with the KS element equations of motion for long-term motion by considering the perturbations due to the combined effect of Earth’s zonal harmonics J₂ to J₄ and atmospheric drag. The model used here is an oblate diurnally varying atmosphere with variation of the scale height depending on altitude almost similar to my work with Sharma [15]. Using series expansion method third-order terms in e, eccentricity, c, a small parameter depending on the ellipticity of the atmosphere and second order terms in μ, gradient of the scale height altitude are collected. Only one of the eight equations is solved analytically to obtain the state vector at the end of each revolution due to symmetry in KS equations. Numerical studies with test cases reveal that there is a good comparison between the analytical (ANAL) and numerically integrated (NUM) values of the position as well as velocity vectors and .

The KS element equations of motion of a satellite under the effect of perturbing potential V and additional perturbing force [10] are

(1)

(2)

(3)

(4)

, specifying attraction between two masses M and m, E, ω,t,r and k² are, respectively,the eccentric anomaly, angular frequency, physical time, radial distance and the gravitational constant.

The perturbing potential V[10] and the aerodynamic drag force [16] per unit mass acting on a satellite of mass m are respectively

(5)

R, equatorial radius, J_n ’s, dimensionless constants known as zonal harmonics. Using Equation (5),

where

As in [16], the effective area of the satellite, CD, the drag coefficient, ρ , the atmospheric density at the point of calculating atmospheric drag force and , the velocity of the satellite relative to the ambient air. If is the velocity of the satellite relative to the Earth’s centre, then where is the velocity of the air relative to the Earth’s centre is assumed to be west to east, Λ ,the rotational rate of the atmosphere about the Earth’s axis and i0, the initial inclination, is the initial perigee radius, the velocity at the initial perigee. Then the drag force per unit mass tangential to the orbit can be written as

Where and

(6)

Following [17] the density function for an oblateness atmosphere together with day-to-night density variation is

(7)

To express cosφ in terms of the true anomaly θ and then in terms of the eccentric anomaly E let and are the right ascension and declination of the sun and the day time bulge respectively, then

. We can write

(8)

(9)

(10)

The scale height H is known to increase with altitude and this variation of H will have an influence upon its motion. The value of H may be taken as =where |μ |< 0.2 and for any particular value of r_p. To sum up, expression for the density, similar to that of Swinerd Boulton [17], is

(11)

In terms of E,

(12)

(13)

(14)

The integrals available in the above theory are of the form

and 0 if either m or n is odd.

Knowing the position and velocity vectors and at the instant t = 0 , the values of and can be computed [10], (pp. 91-92), and by adopting E = 0 as the initial value of the eccentric anomaly, we obtain

In the entire test cases reported here, the values of ω , Right Ascension Node, Ω and mean anomaly, M are 60, 30 and 0 degrees respectively. The value of K², R, J₂, J₃ and J₄ utilized for numerical computations are 398600.8km³s⁻² , 6378.135 kilometer and ,−2.53648D − 06 and -1.52D-06 respectively. Jacchia (1977) atmospheric density model, which is relatively easier to use, is employed to compute the values of ρ_p0 , the density at the perigee and H, the density scale height at the end of each revolution. Arbitrarily 22 August 2002 is chosen as the initial epoch. The values of ε, ∧ and , utilized during the computations are 0.00335, 1.2 and 50.0 respectively. In this model sin i approaches maximum value 0.2042 at e = 0.003 and i= 80° while minimum value 0.0232 at e = 0.005 and i= 20° respectively at the end of 100 revolutions (Figures 1 and 2). The values of 10.7cm solar flux (F_10.7) and averaged geomagnetic index (A_P) are taken as 150 and 10, respectively, which approximately represents an average density and results in exospheric temperatures between 1000 and 1100 K for the different cases we considered (Table 1).

Figure 1: Comparison between NUM and ANAL values at perigee ht. 400 km, i=50 and e from 0.001-0.04.

Figure 2: % error between NUM and ANAL values at perigee ht. 400 km, i = 50 and e from 0.001-0.04.

Table 1: Decrease in Argument of Perigee.

We have transformed equations for and using equations (11), (12), (13), (14) and programmed in double precision arithmetic to compute the KS elements and (i = 1,2,3,4). Once and are known, they are converted into and and , which are further converted to the osculating orbital elements (Figures 3 and 4). Percentage Error = (ANAL – NUM)/NUM is calculated to check validity of the work. The algebraic computations are made with MAPLE12 mathematical software (Table 2).

Figure 3: Comparison between NUM and ANAL values at perigee ht 350 km, e = 0.01 and i from 10-80 degree.

Figure 4: Comparison between NUM and ANAL values at perigee ht.400 km, e=0.1 and i from 10-80 degree.

Table 2: Comparison of ANAL to NUM with % error at perigee ht 400 km, i = 80 after 1and 100 revolutions.

The KS element equations are integrated analytically by a series expansion method by assuming an oblate diurnal atmosphere when density scale height varies with altitude and by including the terms corresponds to Earth’s zonal harmonics J₂, J₃ and J₄. A wide range of eccentricity and inclination is considered for calculating the change in argument of perigee by present analytical theory and by numerical integration. Comparison between analytically and numerically integrated values for 1 and 100 revolutions shows that the analytically integrated values are reasonably accurate and thus highlights the usefulness of the analytical expressions. Graphical representation as well as the table presented here emphasizes the importance of developing the theory to find the decrease in argument of perigee.