The present work aims at constructing an atlas of the balanced Earth satellite orbits with respect to the secular and long periodic effects of Earth oblateness with the harmonics of the geopotential retained up to the 4th zonal harmonic. The variations of the elements are averaged over the fast and medium angles, thus retaining only the secular and long periodic terms. The models obtained cover the values of the semi-major axis from 1.1 to 2 Earth’s radii, although this is applicable only for 1.1 to 1.3 Earth’s radii due to the radiation belts. The atlas obtained is useful for different purposes, with those having the semi-major axis in this range particularly for remote sensing and meteorology.

1. Introduction

The problem of the motion of an artificial satellite of the Earth was not given serious attention until 1957. At this time, little was known about the magnitudes of the coefficients of the tesseral and sectorial harmonics in the Earth’s gravitational potential. It was pretty well known at this time (1957–1960) that the contributions of the 3rd, 4th, and 5th zonal harmonics were of order higher than the contribution of the 2nd zonal harmonic, but the values of the coefficients C30, C40, and C50 were not very well established. No reliable information was available for the tesseral or the sectorial coefficients except that the observations of orbiting satellites indicated that these coefficients must be small, certainly no more than the first order with respect to C20.

For low Earth orbits within an altitude less than 480 km, if the satellite attitude is stabilized, or at least a mean projected area could be estimated, the perturbative effects of atmospheric drag should be included. Unfortunately, the literature is still void of even a mention of this topic of balancing this kind of very low Earth orbits. The reason may be the present increased interest in space communications and broadcasting, which still make use of the geostationary orbits that lie beyond the effects of atmospheric drag, though they still suffer the effects of drift solar radiation pressure.

With the advance of the space age, it became clear that most space applications require fixing, as strictly as possible, the areas covered by the satellite or the constellation of satellites. In turn, fixing the coverage regions requires fixed nodes and fixed apsidal lines. This in turn leads to the search for orbits satisfying these requirements. The families of orbits satisfying such conditions are called “frozen orbits” [15]. Clearly, the design of such orbits includes the effects of the perturbing influences that affect the motion of the satellite. As the present work is interested in low Earth orbits, only the effect of Earth oblateness is taken into concern. These have been extensively treated in the literature [610].

This paper is aiming at constructing an atlas of the balanced low Earth satellite orbits, which fall in the range from 600 km to 2000 km above sea surface, in the sense that the variations of the elements are averaged over the fast angle to keep only long periodic and secular variations that affect the orbit accumulatively with time. In this paper, a model is given for the averaged effects (over the mean anomaly) of Earth oblateness. Then, the Lagrange planetary equations for perturbations of the elements are investigated to get sets of orbital values at which the variations of the elements can be cancelled simultaneously.

2. Earth Potential

The actual shape of the Earth is that of an eggplant. The center of mass does not lie on the spin axis, and neither the meridian nor the latitudinal contours are circles. The net result of this irregular shape is to produce a variation in the gravitational acceleration to that predicted using a point mass distribution. This variation reaches its maximum value at latitude 45 deg and approaches zero at latitudes 0 and 90 deg.

The motion of a particle around the Earth can be visualized best by resolving it into individual motions along the meridian and the latitudinal contours. The motion around the meridian can be thought of as consisting of different periodic motions called “zonal harmonics.” Similarly, the motion along a latitudinal contour can be visualized as consisting of different periodic motions called “tesseral harmonics.” The zonal harmonics describe the deviations of a meridian from a great circle, while the tesseral harmonics describe the deviations of a latitudinal contour from a circle.

At points exterior to the Earth, the mass density is zero, thus, at external points the gravitational potential satisfieswhere V is a scalar function representing the potential. Also, the gravitational potential of the Earth must vanish as we recede to infinitely great distances. With these conditions on the above equation, the potential V at external points can be represented in the following form [11]:

This expression of the potential is called “Venti potential,” and it was adopted by the IAU (International Astronomical Union) in 1961. The terms arising in the above equation are Cnm and Snm are harmonic coefficients (they are bounded as is always the case in physical problems), R is the equatorial radius of the Earth, is the Earth’s gravitational constant, G is the universal constant of gravity, M is the mass of the Earth, and (r, α, δ) are the geocentric coordinates of the satellite (Figure 1, [12]) with α measured east of Greenwich, and represents the associated Legendre polynomials.The terms with m = 0 are called “zonal harmonics.”The terms with 0 < m < n are called “tesseral harmonics.”The terms with m = n are called “sectorial harmonics.”

The case of axial symmetry is expressed by taking m = 0, while if equatorial symmetry is assumed, we consider only even harmonics since P2n+1 (−x) = −P2n+1 (x). Also, the coefficients C21 and S21 are vanishingly small. Further if the origin is taken at the center of mass, the coefficients C10, C11, and S11 will be equal to zero.

Considering axial symmetry, with origin at the center of mass, we can writewhere Jn = −Cn0.

Taking terms up to j4, we can write V in the following form:where

It is a purely geometrical transformation to express the potential function V(r, δ), given by the above equations, as a function of the Keplerian orbital elements a, e, i, Ω, ω, and I in their usual meanings (Figure 2, [12]), where a and e are the semi-major axis and the eccentricity of the orbit, respectively, i is the inclination of the orbit to the Earth’s equatorial plane, Ω and ω describes the position of the orbit in space where Ω is the longitude of the ascending node and ω is the argument of perigee, and finally, l is the mean anomaly to describe the position of the satellite with respect to the orbit.

Then, V (a, e, i, Ω, ω, I) is in a form suitable to use in Lagrange’s planetary equations, and in canonical perturbations methods through the relation. From the spherical trigonometry of the celestial sphere, we havewhere f is the true anomaly, ω is the argument of perigee, and Si= sin i, Ci= cos i.

Substituting for δ, in Pn (sin δ), we get

Thus, we get V2, V3, and V4 as functions of the orbital elements.

We now proceed to evaluate the effects of Earth oblateness, considering the geopotential up to the zonal harmonic J4.

In the present solution, we consider only the secular and long periodic terms, averaging over the mean anomaly l.

2.1. The Disturbing Function

The disturbing function is defined as follows:where V2, V3, and V4 are the 2nd, 3rd, and 4th terms of the geopotential, namelywhere , , and is the true anomaly.

Since terms depending only on the fast variable l will not affect the orbit in an accumulating way with time, we average the perturbing function R over the mean anomaly with its period 2π. The average function is defined by

As the perturbing function R is a function of the true anomaly f not the mean anomaly l, we use the relationwhere both angles have the same period and the same end points 0 and 2π, and therefore the average function will be given by

Applying the required integrals, we get the averaged disturbing function where

2.2. Lagrange Equations for the Averaged Variations of the Elements

The Lagrange planetary equations for the variations of the elements for a disturbing potential R are as follows [13]:where n is the mean motion given by .

Substituting for the averaged disturbing function due to Earth oblateness, the Lagrange equations becomewhere the equation for l is neglected because we concentrate on balancing the orbit position not the satellite motion in the orbit.

2.3. Variations of the Elements due to Earth Oblateness

Substituting for the required derivatives in equations (16) to (20) yields


We getor


We get

Defining as in equation (28), we get equation (29):

Equations (22)–(29) give the average effects of Earth oblateness including the zonal harmonics of the geopotential up to J4 on the Keplerian elements of the satellite orbit.

3. Balanced Low Earth Satellite Orbits

In what follows, we try to find orbits that are balanced in the sense that the averaged (over the fast variable l) variations of the orbit elements are set equal to zero. In equation (23), we put it equal to zero and get a relation between the argument of perigee ω and the inclination i, while treating the eccentricity e and the semi-major axis a as parameters. This will give a range of values for ω and i at different values of e and a, which are all give balanced orbits with respect to both the eccentricity e and the inclination i. The same is done with equations (26) and (29), while putting  = 0 and  = 0.

The applicable ranges for this model of the semi-major axis a are , where the range 1.4R–2R is avoided due to the predominance of the radiation belts at these levels, to avoid the damages of the equipment that it may produce, besides its fatal effects on human life (for inhabited spacecrafts). The values for the eccentricity e are taken between 0.01 (almost circular orbit) and 0.5.

3.1. Orbits with Fixed Eccentricity and Inclination

By equating the variation of e by zero, we get

This implies

So either orwhere

Equations (32) and (33) give the family of low orbits that have both the eccentricity and the inclination fixed.

The condition for the existence of such orbits is clearly thatwhich is guaranteed by

We put it as a condition on the eccentricity i only when

3.2. Orbits with Fixed Node

The families of orbits for which dΩ/dt = 0 (i.e with fixed nodes) are obtained from

This implies

This can be arranged as a second order equation for sin(ω), giving a family of orbits with fixed argument of perigee for different values of ω as a function of the inclination i, the semi-major axis a, and the eccentricity e.

When solving for sin(ω), we get the family of orbits for which the longitude of the node is balanced,where

The condition for having the orbit if real solution exists is again that , which gives

3.3. Orbits with Fixed Perigee

For the argument of perigee to balance, we solve .We substitute from equation (26) into equation (29) then expand cos(2ω) and collect terms with respect to sin(ω). We get

Thus for , we getwhere

Equations (44)–(47) give the relation between sin(ω) and the inclination i, the semi-major axis a, and the eccentricity e, which gives the family of orbits that balance the argument of perigee ω, subject of course to the restriction that is a real value between −1 and 1.

4. Numerical Results

In this section, numerical results and graphs are obtained for the case of seasat a = 7100 km by putting ω as a function of e and i from equations (32) and (33). The curves are plotted within the possible range given by condition (36) to give curves of balanced e and i. The curves are against the eccentricity e in the range [0.01, 0.5]. The numerical values involved are J2 = 0.001082645, J3 = −0.000002546, J4 = −0.000001649, R = 6378.165 km, α = 7100 m, and μ = 398600.5 km3sec−2.

The condition (36) gives the upper and lower bounds of the function F(i) as a function of e, and since the function is increasing with e and has no critical points in the interval [0.01, 0.5], then the minimum value occurs at e= 0.01 and the maximum at e= 0.5, which gives . The graphs are plotted for different values of F(i), which corresponds to specific values of i found by solving the equation resulting from setting F(i) equal to the required values. After that we plot  = 0 and  = 0 simultaneously for the same values of F(i) at each specific i-value, to find the orbit values at which we have nearest values for  = 0 and  = 0. In the graphs of and , the relation (32) was kept to ensure that e and i are already balanced.

Five selected values of F(i) are chosen: F(i) = 0.01, 0.5, 0.10, 0.15, and 0.20. Negative values will give the same results with negative sign since F(i) is odd with respect to i. Also for each value, the solution of F(i) = x, with x equals one of the above values we get four real values of i on the form: i1, 180 − i1, i2, and −180 − i2, where i1 and i2 are near the critical inclination one of them is positive and the other is negative. Table 1 gives the values of i corresponding to the selected values of F(i).

We note that as F(i) increases, i gets away from the critical inclination, and the curve gets shorter indicating less stability of ω as expected.

The graphs are plotted for each value of F(i) first for the balanced values of e and i, then for the corresponding four values of i, four graphs are plotted for and to find the nearest values of zero variation for both elements.

Figures 312 show the possibility of balancing ω with e and i, while Ω will have a variation of order 10−6/sec or it must be balanced alone at i = 90 deg (as shown in equation (26)), according to the orbit kind. Figures 3, 5, 7, 9, and 11 show the possible values of ω(e) that balance e and i, while Figures 4, 6, 8, 10, and 12 show the curves of and at the values of i at which .

5. Conclusion

Let balanced orbits be defined as those for which the orbital elements are set equal to zero under the effect of secular and long periodic perturbations. In this work, the effect of Earth oblateness is the considered perturbing force because of dealing with low Earth orbits. The above analysis then shows that such an orbit will be balanced within a reliable tolerance only for few weeks since we are forced to accept the motion of either the node or the perigee by about 10−6 deg/sec. The reason is that under the influence of the Earth oblateness, (exactly) for , while only near the critical inclination deg. Hence, the best procedure is to design a satellite constellation for which the nodal shifts due to the perturbative effects and Earth rotation are modeled to yield continuous coverage. The perigees are either fixed or arranged to realize that the perigee (or the apogee) be overhead the coverage region (regions) in due times. This may require near commensurability with the admitted nodal periods.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was funded by the University of Jeddah, Saudi Arabia, under grant no. UJ-26-18-DR. Thus, the author therefore acknowledges with thanks the university technical and financial support.