Effects of Geopotential and Atmospheric Drag Effects on Frozen Orbits Using Nonsingular Variables

The concept of frozen orbit has been applied in space missions mainly for orbital tracking and control purposes. This type of orbit is important for orbit design because it is characterized by keeping the argument of perigee and eccentricity constant on average, so that, for a given latitude, the satellite always passes at the same altitude, benefiting the users through this regularity. Here, the system of nonlinear differential equations describing the motion is studied, and the effects of geopotential and atmospheric drag perturbations on frozen orbits are taken into account. Explicit analytical expressions for secular and long period perturbations terms are obtained for the eccentricity and the argument of perigee. The classical equations of Brouwer and Brouwer and Hori theories are used. Nonsingular variables approach is used, which allows obtaining more precise previsions for CBERS (China Brazil Earth Resources Satellite) satellites family and similar satellites (SPOT, Landsat, ERS, and IRS) orbital evolution.


Introduction
The orbital motion of an the artificial satellite free of disturbing effects, named Keplerian motion, would be ellipses with constant sizes and eccentricities, in permanent planes, where the satellites would stay indefinitely.Some of the major perturbations that affect the orbital elements of an artificial Earth satellite are as follows: nonhomogeneity of Earth's mass distribution, sun-moon gravitational attraction, solar radiation pressure (direct and albedo), atmospheric drag, forces due to Earth's tides, Poynting-Robertson drag, and Yarkovsky effect.
For orbital control purposes, it can be important that some orbital elements stay frozen, that is, with near constant values, in order to facilitate some maneuvers adjustment.Especially for maneuvers that have been carried out at INPE (National Institute for Space Research) Satellite Tracking and Control Center (TSCC), with CBERS family, it is important that eccentricity and argument of perigee stay frozen.
The notion of frozen orbits goes back to many years, the history of which is well explained in [1,2].Due to their scientific interest, frozen orbits concept has been also investigated around planetary satellites and asteroids [3,4].Many investigators have contributed research; in order to predict the long-term evolution of the eccentricity and argument of perigee without numerical integration, an analytical method based on the Delaunay variables was employed in [5]; and a simple solution has been developed for the longterm behavior of a near-circular orbit in a zonal gravity field in [6], linearizing the singly averaged variational equations of motion and eliminating a degree of freedom with an integral of motion.The constraint equation for Earth frozen orbits, based on the Lagrange's planetary equations, can be obtained for the gravity model involving higher degree zonals in [2].
Basically, the odd terms of geopotential, responsible for the long period effects, give rise to the frozen orbits.The influence of geopotential perturbations on frozen orbits has been studied by several authors, including Cutting et al. [7], where secular and long period terms of geopotential up to  3 are considered.
In the present work, such study is extended in order to verify the influence of perturbations due to geopotential up to  5 terms and atmospheric drag on frozen orbits.Nonsingular variables [8] are also introduced, turning analytical formulas suitable for very low eccentricity orbits application.Brouwer theory [9] is used to express analytically secular and long period terms due to geopotential for eccentricity and argument of perigee.The influence of atmospheric drag is also analytically expressed, using Brouwer and Hori theory [10].
In artificial satellite theory, several analytical models have been proposed to describe the atmospheric density [11].However, in general, the more realistic are the models, the more difficult is the solution of the motion equations.Both Brouwer and Brouwer and Hori theories are convenient for analytical developments and furnish reasonable orders of magnitude accuracy.
These analytical expressions enable carrying out the temporal behavior analysis of eccentricity and argument of perigee in the neighborhood of a frozen orbit, as well as calculating the magnitude of the perturbations due to geopotential and atmospheric drag.
Such development will allow more precise predictions for the orbital maintenance and evolution of CBERS family and similar satellites.

Frozen Orbits Concept
A frozen orbit is characterized by keeping (or trying to keep) the argument of perigee and the eccentricity of an orbit constant, so that, for a given latitude, the satellite always passes at the same altitude, benefiting the users through this regularity.In another way, this type of orbit maintains almost constant altitude over any point on Earth surface.The design of frozen orbits involves selecting the correct values of eccentricity and argument of perigee, for a given semimajor axis and orbital inclination (in which critical value of 63.43 ∘ is an intrinsic singularity in the artificial satellite theory [12]).
The analytical nonlinear system of equations considering geopotential perturbations up to  3 terms are as follows [7]: where a, e, i, and  are the classical orbital elements: semi major axis, eccentricity, inclination, and argument of perigee, respectively;   is the Earth equatorial radius;  = √/ 3 is the satellite mean motion;  is the Earth gravitational constant; and  2 and  3 are the second and third gravity coefficients, respectively.Equation (1) shows that perturbations effects in the eccentricity vanish for values of  equal to 90 ∘ and 270 ∘ .

Secular Perturbations and Long Period
Terms Using Nonsingular Variables Brouwer [9] and Brouwer and Hori [10] analytical developments allow analyzing the effects of the perturbations due to geopotential up to  5 order and to atmospheric drag on the neighborhood of a frozen orbit.Such equations are expressed in terms of the Keplerian orbital elements and are related to the disturbances on eccentricity and argument of perigee.However, these equations are placed so that singularity points may occur.Especially for the eccentricity values used by INPE, which applies frozen orbits theory on CBERS family, it indeed occurs [13,14].With a view to eliminate the singularity points, a variables change is needed.Brouwer and Brouwer and Hori analytical equations were rewritten using nonsingular variables approach.As consequence, eccentricity and argument of perigee no longer explicitly appear in the rewritten equations.
The nonsingular variables,  and , used by INPE TSCC and defined in terms of the eccentricity and the argument of perigee, are as follows [8]: Consequently, the eccentricity and the argument of perigee are given by In order to carry out the variables change for eliminating the singularity points, other trigonometric expressions were rewritten in terms of nonsingular variables: The analytical equations for the frozen orbits dynamics, ω and ė , are presented.They were developed according to Brouwer [9] and Brouwer and Hori [10] theories and take into account only secular and long period terms of the perturbations.Due to the requirements on the accuracy of orbit determination, it is convenient to write such equations using nonsingular variables.This is because the eccentricity values are usually very small and, for some types of orbits, are near to zero, which results in indetermination in the equations of motion, if placed as the equations of Brouwer and Brouwer and Hori theories.The rewritten equations using nonsingular variables approach are presented in the Appendix, since the formula are too large to be displayed explicitly in the text.

Results and Analysis
The solutions for ω and ė equations, which are explicitly written in the Appendix, were obtained using CBERS-1 satellite real orbital data.Perturbations due to geopotential up to  5 and atmospheric drag were analytically included and rewritten in terms of nonsingular variables.Results were plotted using a phase plane of eccentricity and argument of perigee, which is analogous to the system state evolving over time.Here, the period of time corresponds to 300 days.Such results allowed evaluating frozen orbits behavior.
In this application, the orbit was considered frozen when the argument of perigee was  = 90 ∘ .The graphics afterwards show that the orbit starts escaping from its initial position if the arguments of perigee values are much different from 90 ∘ and that it is inclined to be confined and cycle-limited when such values are closer to 90 ∘ .
First, a dynamic model including perturbations due to geopotential up to  3 only was considered.In a second step, the model took into account the effects due to geopotential up to  5 .Finally, a model including perturbations due to geopotential up to  5 terms and atmospheric drag was analyzed.The drag effects analysis occurred for two different geopotential model accuracy levels: up to  3 terms only and up to  5 terms.
With the dynamical model development (i.e., a model consisting of more disturbing effects), a greater accuracy was expected for the frozen elements.First of all, a greater prevision capacity of these elements variations was expected, when compared to the results including terms up to  3 only.
In fact, the amplitudes of variation of argument of the perigee and eccentricity decreased when geopotential terms up  5 and atmospheric drag effects were included, as it will be shown next.
Figures 1 and 2 show curves obtained for five different initial conditions.Such graphics were based on the solution of ( 2) and ( 4), both rewritten in nonsingular variables.Figure 1 shows the results for two models: geopotential up to  3 only and geopotential up to  3 included in atmospheric drag.In Figure 2, the results for other two models are presented: geopotential up to  5 and a more complete dynamic model that takes into account geopotential up to  5 and atmospheric drag.Regarding the atmospheric drag, the solar flux used was F10.7 of 250, which denotes the maximum solar activity and manifests itself on augmented atmospheric densities and, consequently, on higher atmospheric drag effects.
From Figures 1 and 2, it is noticeable that the curves of Figure 2 have less amplitude of variation for the frozen elements than the first one.Besides, if the two graphics of Figure 1 are compared amongst themselves, the inclusion of the atmospheric drag perturbations also decreases eccentricity and argument of perigee amplitudes of variation, despite being less meaningful than the one caused by the inclusion of geopotential terms up to  5 .The same behavior is detected when the atmospheric drag is added to the perturbations due to geopotential up to  5 , as Figure 2 showed.Table 1 shows the maximum amplitudes of variation of eccentricity and argument of perigee, Δ and Δ, for five initial conditions.From Table 1, it is possible to observe that the magnitudes decrease when the perturbations due to geopotential up to  5 and atmospheric drag are included in the dynamic model, although the greater decreases occur with the inclusion of perturbations due to geopotential up to  5 .
From Table 1, the amplitude of variation of the argument of perigee decreases when geopotential terms up to  5 and atmospheric drag terms are included.This behavior happens again for all initial conditions.If the argument of perigee initial values are far from the frozen orbits conditions ( = 90 ∘ ), the cycles including  5 and atmospheric drag are still complete, while the cycles including  3 only start demeaning, as can be seen in Figure 1.In other words, the inclusion of the effects due to geopotential up to  5 and atmospheric drag improves the prevision precision of the argument of perigee.With respect to the eccentricity, the reduction is subtler, but it still occurs, since one of the atmospheric drag effects is to circularize the orbit (because the drag acts reducing the semimajor axis of the orbit).
In practice, the theory using only terms factored up to  3 can persuade to wrong needs of orbital maneuver corrections.If it is supposed, for instance, that the mission requires a perigee value between 90 ∘ ± 10 ∘ , the  3 dynamic model would foresee a corrective maneuver, while the model including  5 and atmospheric drag excludes the need of a maneuver, as can be seen on the first line of Table 1.So, the conclusion is that it is imperative for the inclusion of geopotential terms up to  5 and atmospheric drag terms to improve the precision on planning of maneuvers carried out by INPE TSCC.
As shown in Figure 2 and Table 1, the more meaningful decreases on the frozen elements amplitudes of variation occurred when the geopotential was factored up to  5 .This way, Figure 3 isolates the geopotential up to  3 and up to  5 effects from the atmospheric drag ones, and only results due to geopotential will be shown.Because of the major interest on the frozen orbits conditions, only  equal to 90 ∘ and 100 ∘ will be considered from this point on.Figure 3 phase plane shows that including geopotential terms up to  5 causes a great diminution on the magnitude of eccentricity and argument of perigee amplitude of variations.
The atmospheric drag effects in eccentricity and argument of perigee over time are emphasized in Figures 4  and 5. Secular and long period components, whose outline is senoidal with exponential overlay, are noticeable.The dynamic model includes disturbing effects due to geopotential up to  5 and atmospheric drag.There was a maximum solar activity (solar flux 250), and again only  equal to 90 ∘ and 100 ∘ results were analyzed.
As the analytical development of geopotential and atmospheric drag perturbations may lead to singularities, the equations were rewritten using nonsingular variables approach.Figure 6 was built in an effort to verify whether ω and ė equations were correctly written in terms of nonsingular variables.The results obtained from analytical equations of geopotential up to  5 , that is, developed using Brouwer theory, were compared to the results obtained from the equations written in nonsingular variables.As can be verified in Figure 6, the transformation to nonsingular variables was developed correctly, since the curves plotted extremely resembles.In fact, the curves even superpose each other.And, considering that the interest rests on frozen orbits concept, only  being equal to 90 ∘ and 100 ∘ was considered in Figure 6.

Conclusions
Analytical expressions for the frozen orbits, containing secular and long period components of the orbital perturbations,   were obtained explicitly.To the dynamic model including geopotential up to  3 terms only were included in perturbations due to geopotential up to  5 and atmospheric drag.
The inclusion of such perturbations followed Brouwer and Brouwer and Hori theories.The perturbations analytical expressions were then rewritten in terms of nonsingular variables, in order to avoid singularities due to the eccentricity values used by INPE TSCC.Since the analytical equations developed for frozen orbits maintenance showed a small additional computational burden, the inclusion of terms due to perturbations of geopotential up to  5 and atmospheric drag in INPE TSCC operation software is justified.The inclusion of geopotential perturbations up to  5 was responsible for the major decrease in the frozen elements amplitudes of variation, for any initial condition.The atmospheric drag introduces secular and long period components with senoidal outline of exponential overlay amplitude in eccentricity and argument of perigee.Whether to take or not into account the atmospheric drag showed relevant impact for long periods of time.
The frozen elements smaller amplitudes of variation, when  5 terms and atmospheric drag perturbations are included, improve the precision for prediction of eccentricity and argument of perigee.This means an enhanced precision not only on maneuver's computations but also on maneuver's prediction, which contributes to a better performance in the orbital operations conducted at the INPE TSCC.
Frozen orbits feature is extremely important for the users of the images obtained by the camera onboard CBERS family type of satellite, which has a sun-synchronous orbit with an altitude of 778 km, completing about 14 revolutions per day.The local solar time at the equator crossing is always 10:30 a.m., thus providing the same solar illumination conditions for comparing images taken in different days [15].The images from different days can be compared for the same latitude and then be used to preview harvests, to detect fire on forests, and to locate underground airports, as well as others utilities.

Equations of Motion
The analytical equations for the frozen orbits dynamics, ω and ė , which were developed according to Brouwer and Brouwer and Hori theories and were rewritten using nonsingular variables approach are presented as follows:

Figure 1 :
Figure 1: Phase plane diagram considering perturbations due to geopotential up to  3 and due to  3 and atmospheric drag for 300 days.

Figure 2 :
Figure 2: Phase plane diagram considering perturbations due to geopotential up to  5 and due to  5 and atmospheric drag for 300 days.

Figure 3 :
Figure 3: 300-day phase plane diagram for perturbations due to geopotential up to  3 and  3 +  5 .

Figure 4 :
Figure 4: Eccentricity variation in terms of secular and long period components.

Figure 5 :Figure 6 :
Figure 5: Argument of perigee variation in terms of secular and long period components.

Table 1 :
Frozen elements maximum variation for each dynamic model.