INTRODUCTION
The drag force is one of the major perturbing forces on satellites in low earth orbit where the atmospheric density is the most variable element. A Low Earth Orbit (LEO) satellite is defined as a satellite that orbits below an altitude of approximately 1000 km (Samwel, 2014). Although the air density is much lower than that near the earth’s surface, the air resistance in those layers of the atmosphere, especially at altitudes ranging from 300800 km above the earth’s surface, is still strong enough to produce drag and pull the satellites closer to the earth. Hence, disregarding the geopotential perturbing forces, the air drag force is considered to be the predominant among other perturbing forces as stated by Eshagh and Alamdari (2007).
Many efforts were made in the past to model the air drag using different mathematical solution techniques. Brouwer and Hori (1961) concerned the presence of the oblateness perturbations and the drag effect in a single solution. They used the canonical variables which have particular advantages in the treatment of this problem because of the ease with which the necessary transformation of variables can be made. They were the first to treat the joint effects of drag and earth oblateness. Thereafter, Delhaise (1991), developed an analytical solution for the motion of artificial earth satellites subjected to the combined effects of earth oblatness and atmospheric drag.
Khalil (2002) has proposed a theory of an Artificial Earth Satellite under the joint effects of earth oblateness and atmospheric drag. The Hamilton’s equations of motion are derived and the atmospheric model was taken as an oblate rotating model. Also, Saad et al. (2008) used KS regular variables to predict orbital lifetime numerically. They considered the perturbations due to the nonaspherecity of the earth and the atmospheric drag.
The impact of the solar activity variation on the satellite decay in LEO was addressed in many studies. Sehnal (1975) paid attention to the difficulty to obtain a precise analytic continuous solution due to the possible and unexpectable change of the density with the activity of the sun and due to the rapid increase of the drag constant above 400 km from the surface of the earth. Pardini et al. (2006) analyzed the orbital decay of eleven spherical satellites in the 1501500 km altitude range over a full solar activity cycle to highlight the performances of the JR71 and MSISE90 density models. They concluded that the two models were not able to correctly represent the air density at each altitude and environmental condition. They found that below 400 km, both models overestimated the atmosphere density. However, MSISE90 was proved to be the best model to compute the air density, in low solar activity conditions while JR71 was generally more accurate over the maximum of the 23rd solar cycle. Lately Nwankwo and Chakrabarti (2013) studied the effects of the adverse space weather condition induced by solar events, especially the solar flares and Coronal Mass Ejections (CMEs) on the orbital decay of low earth satellite orbits. They showed that the result strongly depends on the phase of the solar cycle. They also showed that a major CME event can cause sufficient heating and expansion of the atmosphere so that the orbital radius may go down by a few km in a single event.
The unwanted orbital changes occur due to the drag force are known to be related to atmospheric density changes which are primarily induced by solar extreme ultraviolet irradiance variations (Tobiska, 2003). This ambient atmosphere is considered to be important environment at the LEO and polar earth orbit altitudes. For the other orbits, it is too low to be of significance (Hasting and Garrett, 1996).
Therefore, in the present study, the variations of the orbital elements which define the orbital trajectory of the satellite due to the effect of the air drag force at LEO environment during maximum and minimum solar activity are evaluated. The variation of the rotation velocity of the atmosphere is considered. The Lagrange planetary equation in Guass form is used to calculate the variation of the orbital elements numerically and the NRLMSISE00 model which maps out the atmosphere from sea level to 1000 km is used to estimate the density of the atmosphere.
As a numerical application, the orbital elements of the Hubble Space Telescope (HST), CORONAS I and PRIRODA satellites, which express the Low Earth Orbit (LEO) satellites are used. The variation of the orbital elements of these satellites under the effect of the drag force for both maximum and minimum solar activity are calculated.
METHODOLOGY
Atmospheric Density model (NRLMSISE00): The NRLMSISE is an empirical, global model that models temperatures and densities of atmospheric components from the ground to the exosphere. It is developed by the U.S. Naval Research Laboratory (NRL) and based on the data of mass spectrometer and Incoherent Scatter Radar (Hedin, 1991; Picone et al., 2002). The main advantage of the usage of these datasets is that they consist of independent observations of both temperature and number densities for the atmospheric constituents. The MSIS86 model replaced Jacchia71. An extension to the model was published as MSISE 90 (Hedin, 1991), which is identical to MSIS86 for the thermosphere region, but extends down to zero height. At the end of the 1990s, NRLMSISE00 model was released (Picone et al., 2002) including additional mass spectrometer and incoherent scatter radar data, accelerometer data, as well as database of Drag Temperature Model (DTM) and Jacchia satellite orbit decay (Doornbos, 2012).
The NRLMSISE00 comprises the main drivers of the upper atmosphere: The solar extreme ultraviolet (EUV) flux and geomagnetic heating. The 10.7 cm solar radio flux (F_{10.7}) is an indicator of the overall solar activity level. It measures the radio flux at a wavelength of 10.7 cm in units of 10^{–22} W m^{–2} Hz^{–1} and is used as a proxy of the extreme ultraviolet radiation (Samwel, 2014) while the A_{p} daily geomagnetic index measures the geomagnetic component of space weather. The F_{10.7} index is obtained from solar data at the National Geophysical Data Center (NGDC) http://www.swpc.noaa.gov/ and The National Oceanic and Atmospheric Administration world data center for geomagnetism in Kyoto http://wdc.kugi.kyotou.ac.jp/.
The NRLMSISE00 has also been compared to the previous MSIS models and Jacchia70 and represented a noticeable improvement (Picone et al., 2002). In the present study, the NRLMSISE00 model for density calculations is used. The variation of the total mass density (ρ) with the true anomaly (f), during the maximum and minimum solar activity, is obtained and fitted using Fourier series of degree 4 as represented in Eq. 1. The values of the fitting coefficients (a_{0}, a_{i} and b_{i}) are outlined in Table 2:
Drag force: The aerodynamic drag force per unit mass acting on a satellite can be represented by:
where, is the velocity of the satellite relative to the atmosphere, C_{D} is a dimensionless drag coefficient, A/m is the area to mass ratio of the satellite and ρ is density of the atmosphere.
The Ballistic number (B) which is obtained from the TLE outlined in SpaceTrack website, is used (https://www.spacetrack.org/) where:
It is usually assumed that the atmosphere rotates with the same angular velocity as the earth. Hence, its velocity at a position vector But this assumption implies that at infinite distance from the earth, particles are coupled to it in the same way as those close to its surface. So, it will be assumed that the atmosphere rotates with the angular velocity:
which realizes that as:
and as:
If the satellite is at with instantaneous velocity then:
If are the base vectors of the rectangular geocentric equatorial frame associated with spherical coordinates then:
where:
and:
where, are unit vectors in the direction respectively.
In the case of Lagrange’s planetary equations, the perturbing force decomposed into radial, transverse and normal directions is needed. The unit vector normal to the orbit is collinear with the angular momentum unit vector
with:
where, μ is the earth’s gravitational parameter, a is the semimajor axis, e is the eccentricity of the satellite orbit and the transverse unit vector can be calculated from the righthanded set:
From Eq. 411 we get:
and:
The drag force in the radial, transverse and normal components (R, T, N) are:
where:
where, i, ω, f are the inclination of the orbit on the equator, argument of the perigee and the true anomaly of the spacecraft orbit, respectively.
Lagrange planetary equations with air drag force: The Lagrange planetary equations are useful only for evaluating forces that can be presented in terms of a potential. However, there is a Gaussian form of these equations that can be used to study the perturbing forces (air drag, in our case) (Brouwer and Clemence, 1961).
The evaluation of such perturbing forces analytically can become quite complex. Therefore, we used the numerical integration to calculate the variation of the orbital elements due to air drag force. The Lagrange’s planetary equations in Gauss form are:
where, R, T, N are defined by Eq. 1416, Ω is longitude of ascending node, M is mean anomaly and n is mean motion:
The above system of differential equations (Eq. 2328) numerically is solved, using the fourth order RungeKutta method to get variations of orbital elements (Δa, Δe, Δi, ΔΩ, Δω, ΔM) as functions of true anomaly, due to the impact of the drag force. Introducing the true anomaly f as follows:
RESULTS AND DISCUSSION
In the present section, the total atmospheric density using the NRLMSIS00 model during both maximum and minimum solar activity is calculated, taking three LEO satellites to express the LEO environment at different altitudes. In addition, the variations of the orbital elements due to the impact of air drag force are calculated.
Atmospheric density calculations: The NRLMSIS00 is used to calculate the total atmospheric density for the maximum and minimum solar activity for the Hubble Space Telescope (HST), CORONAS I and PRIRODA satellites, as a representative of the LEO satellites.
The NRLMSIS00 model needs a time input to reflect the current state of the atmosphere. It also accounts for the fluctuations in the atmosphere from latitude, longitude, F_{10.7} (Previous day), F_{10.7} (81 day average) and the geomagnetic index A_{p}.
The Two Line Elements (TLE) of the three satellites are obtained from SpaceTrack website https://www.spacetrack.org/ and their orbital elements are outlined in Table 1 during maximum and minimum solar activity.
The maximum level of solar activity employed for this simulation, is taken to start at the orbital epoch 2000, March 31. Thus, considering F_{10.7} (Previous day) = 205.1×10^{–22} Wm^{–2} Hz^{–1}, F_{10.7} (81 day average) = 189.1593 and A_{p} = 23.
The Minimum level of solar activity employed for this simulation, is considered to start at the orbital epoch 1997, March 31. Thus, considering F_{10.7 } (Previous day) = 73.8×10^{–22} Wm^{–2} Hz^{–1}, F_{10.7 }(81 day average) = 73.56049 and Ap = 5.
Figures 1 and 2 draw the total atmospheric density variation with the true anomaly through one orbital revolution for the three satellites during solar maximum and minimum respectively, using the NRLMSIS00 model through the SPENVIS website https://www.spenvis.oma.be/. The figures represent that the total atmospheric density is higher during the epoch of maximum solar activity than that during minimum solar activity and that it is increasing with increasing altitude.
The data obtained is fitted using Eq. 1, The values of the fitting coefficients obtained in Eq. 1 for the three satellites are outlined in Table 2.
Orbital element variations due to air drag force impact: In the present subsection, the variation of orbital elements (Δa, Δe, Δi, ΔΩ, Δω, ΔM), due to air drag force, as function of true anomaly are obtained by solving Eq. 2328 numerically for the three satellites mentioned above. The orbital element variations are calculated during both maximum and minimum solar activity.
The variation in the semi major axis and eccentricity due to the air drag force during the maximum and minimum solar activity for the three satellites, HST, CORONASI and PRIRODA are represented in Fig. 35, respectively. Figure 3 represents the orbital variation of two orbital elements, namely, the semi major axis (a) and eccentricity (e) of the HST, for one revolution, during both maximum and minimum solar activity. In addition, Fig. 4 and 5 represent the variation in the same orbital elements for the satellites CORONASI and PRIRODA, respectively. The orders of magnitude of Δi, ΔΩ, Δω and ΔM are outlined in Table 3.
Table 1:  Orbital elements of HST, CORONAS I and PRIRODA during maximum and minimum solar activity 


Fig. 1(ac): 
Total atmospheric density variation with the true anomaly through one orbital revolution for (a) HST, (b) CORONASI and (c) PRIRODA satellites during the maximum solar activity 
Table 2:  Coefficients of the fitting using Fourier series of the total mass density for the three satellite HST, CORONASI and PRIRODA 

Eshagh and Alamdari (2007), in their study, represented that the air drag force is reducing the semi major axis of the satellite CHAMP by 6 m in one revolution which is comparable to the result obtained in the present study for the satellite CORONASI at the epoch of maximum solar activity. The variation of the other orbital elements due to the air drag force obtained by Eshagh and Alamdari (2007), are nearly similar to the values obtained in the present study as represented in Table 3 and Fig. 4. They used the direct onorbit measurements of the nongravitational satellite acceleration instead of the air density models and used the highlow technique of satellitetosatellite tracking to determine the precise orbit of CHAMP which was launched in summer 2000 into a circular polar orbit at an altitude of about 454 km. As it is represented in Table 1, satellite CORONASI is orbiting the earth in a nearly polar orbit with an altitude of about 450 km during the epoch of solar maximum.
From Table 3 and Fig. 35, it is found that the differences in the variations of the orbital elements at maximum and minimum solar activity are about three orders of magnitude.

Fig. 2(ac): 
Variation of the total atmospheric density with the true anomaly through one orbital revolution for (a) HST, (b) CORONAS I and (c) PRIRODA satellites during the minimum solar activity 

Fig. 3(ad): 
Variations of the semi major axis and eccentricity during (ab) Maximum solar activity and (cd) Minimum solar activity for HST 

Fig. 4(ad): 
Variations of semi major axis and eccentricity during (ab) Maximum solar activity and (cd) Minimum solar activity for CORONAS I satellite 

Fig. 5(ad): 
Variations of semi major axis and eccentricity during (ab) Maximum solar activity and (cd) Minimum solar activity for PRIRODA satellite 
Table 3:  Order of magnitude of (Δi, ΔΩ, Δω and ΔM) due to the air drag force during maximum and minimum solar activity for the three satellites, HST, CORONAS I and PRIRODA 

This difference may be interpreted as when the sun is active, the sun adds extra energy to the atmosphere, hence the low density layers of the atmosphere at the LEO altitudes rise and are replaced by higher density layers that were previously at lower altitudes. As a result, the satellite flies through higher density layer and experiences a stronger drag force.
In addition, Fig. 35 and Table 3 clarify that the variation of the orbital elements are much larger in a case of low altitude satellites such as PRIRODA, incomparable to higher altitude satellites, such as CORONAS I and HST (Table 1). This can be interpreted that the density of the atmosphere decreases exponentially with altitude and consequently the air drag force has a larger impact on satellites orbiting in lower altitude that at higher altitudes.
CONCLUSION
In the present study, The variations of the orbital elements of the satellites in Low Earth Orbit (LEO) environment due to the effect of the air drag force during maximum and minimum solar activity are addressed. The variation of the rotation velocity of the atmosphere was considered. The Lagrange’s planetary equations in the Gauss form to analyze the time rates of classical orbital elements resulting from the total air drag force is employed. The (NRLMSISE00) model was used to estimate the density of the atmosphere. As a numerical application, HST, CORONAS I, PRIRODA were taken to express the LEO satellites.
It is found that:
• 
The variations of the orbital element during the epoch of the maximum solar activity is much larger than those during the epoch of the minimum solar activity. These differences in the variations reach three orders of magnitude 
• 
Due to the air drag force, the variations of the orbital elements for a low altitude satellite is much larger than those of higher altitude satellites 
Therefore, the present study shows that the air drag force has a larger impact on the satellites during the maximum solar activity and on lower altitude satellites.
The present study, is a preliminary investigation of the orbital element variations for one revolution and will be extended for more than one revolution for a comprehensive study for a future work.