Two semianalytic solutions for a perturbed twobody problem known as Lagrange planetary equations (LPE) were compared to a numerical integration of the equation of motion with same perturbation force. To avoid the critical conditions inherited from the configuration of LPE, nonsingular orbital elements (EOE) had been introduced. In this study, two types of orbital elements, classical Keplerian orbital elements (COE) and EOE were used for the solution of the LPE. The effectiveness of EOE and the discrepancy between EOE and COE were investigated by using several near critical conditions.The near one revolution, one day, and seven days evolutions of each orbital element described in LPE with COE and EOE were analyzed by comparing it with the directly converted orbital elements from the numerically integrated state vector in Cartesian coordinate. As a result, LPE with EOE has an advantage in long term calculation over LPE with COE in case of relatively small eccentricity.
1. INTRODUCTION
The motion of the satellite around the Earth can be expressed as a simple Newtonian equation in the Earth centered coordinates. The two bodies rotating around each other without any external forces can be described by the exact mathematical solution. We generally call this simple analytically solvable problem a twobody problem. However, the motion of the satellite in real environment is also dictated by the gravitational and nongravitational forces due to the irregular shape of the Earth, other celestial bodies, the Earth’s atmosphere, Solar wind, and planetary radiation, etc. These forces other than the point mass gravitational force are called perturbation forces. Even before Classical Dynamics and Mathematics were well defined, some of giants in the science left a big footprint in describing the satellite (planetary)motion with the perturbation (Broucke & Cefola 1972, Taff 1985).
Several methods to calculate the evolution of the satellite’s orbit governed by the perturbation were discovered since the Principia of Newton. The Lagrange and Gauss planetary equations (LPE and GPE) are two most recognized methods to describe the change rates of orbital elements with the partials of the disturbing forces. With these theories, Kozai (1959) and Brouwer & Clemence(1961) derived the semianalytical solutions for the short period, long period, and secular perturbations with the variety of mathematical techniques.
With the recent development of the modern computer,many of precision orbit propagation tasks have been processed with direct numerical integration (DNI) of the Newtonian equation of motion as the convention of the Cowell’s orbit integration method. The DNI method has achieved one centimeter level of precision orbit propagation (Luthcke et al. 2003). One day arc to even over ten years arc of satellite’s orbit can be propagated by the DNI method to be estimated for the analysis of the perturbations (Nerem et al. 1993). However, the degradation of orbit propagation due to the integration and numerical errors is inevitable as the integration time span grows. Also, the computation time increases as the requested integration time span is getting inflated. To avoid the several inconvenient factors mentioned previously, the analytic or semianalytic orbit propagators are still investigated and applied in many Celestial Mechanics and Orbital Mechanics problem.
As the development of a domestic satellite orbit propagator and estimator Korea Astronomy and Space Science Institute Orbit propagator and Estimator (KASIOPEA) is right on track (Jo 2008), the characteristic figures of the equations of motion, invariants, and primary variables to be considered must be investigated.
In this study, two types of LPE orbital solutions with classical Keplerian orbital elements (COE) and EOE were reviewed. Before the effectiveness of EOE and the discrepancy between EOE and COE were investigated by using several near critical conditions, a result from DNI had been tested for a reference by the twobody analytic solution. Two types of solution and integration were performed for three different lengths of time. The 13,000 seconds, one day, and seven days evolutions of each orbital element described in LPE were analyzed by comparing it with the directly converted orbital elements from the numerically integrated state vector in Cartesian coordinate.
2. LAGRANGE PLANETARY EQUATION WITH CLASSICAL ORBITAL ELEMENTS
The general theory for finding the rates of change of the osculating elements has been known as the LPE of motion, or simply the Lagrangian variation of parameters (Vallado 2004). In twobody problem it is confirmed that the position and velocity components (in Cartesian coordinates) at any instant permit the determination of a unique set of six COE. And these COE do not change with time. However, most of problems we encountered in the satellite motion have the relatively larger acceleration from the primary body and perturbing accelerations from the others. The equation of the motion of a satellite may be put in a form of Cartesian rectangular coordinates: (
ẍ,
ÿ,z) It is a system of differential equations of the sixth order. Also, these equations of the rectangular coordinates of the elliptic motion can be expressed in terms of the time and six constants (integrals) of integration. As a result, obtained six equations, each of the first order, are exactly equivalent to the original three equations, each of the second order (Brouwer & Clemence 1961). With the six COE from these equations, the equations of the variation of parameters with respect to the small conservative disturbing forces can be obtained as Eqs. (16).
where,

asemimajor axis

eeccentricity of orbit

iinclination of orbital plan

ωargument of perigee

Mmean anomaly

Ωascending node

nmean motion

Rdisturbing function
Using this LPE, we can investigate the perturbation characteristics of specific orbital element with the specific perturbing force. In our case, the time dependent changes of the orbital elements due to the
J_{2}
zonal nonspherical harmonic gravitational perturbation are analyzed. For Earth gravity perturbations, the disturbing function is merely the negative of the disturbing part of geopotential (Blitzer 1975). The disturbing function can be described as in Eq. (7).
where,

μ gravitational constant multiplied by

the mass of the Earth

= 3.986004418×1014(m3/sec2)

(Vallado 2004)

rradius of satellite

Φlatitude of the subsatellite point

R2radius of the Earth

P2Legendre polynomial degree 2

=2/1(3sin2Φ1)
However, the disturbing function,
R
is not a function of COE. The partial derivatives of R in Eqs. (16) should survive only when
R
=
R
(
a, e, i,
ω, M, Ω). The
r
and
sin φ
can be substituted with the functions of
a, e, i, υ, ω
as in Eqs. (8) and (9).
where,
υ true anomaly
With the perturbation by the
J_{2}
term, only the first order secular and the short period disturbing function are survived (Kozai 1959). After several variable substitutions, averaging, and first order approximation, Kozai (1959) and Brouwer & Clemence (1961) successfully solved the variation of equation, Eqs. (16). The solutions for short period and secular terms derived by Kozai (1959) are presented in Appendix A.
3. LAGRANGE PLANETARY EQUATION WITH EQUINOCTIAL ORBITAL ELEMENTS54
With the COE in LPE, there can be some inconvenience with specific types of orbits. As you can see in Eqs. (26), small eccentricities or small inclinations can make the equation unstable. The eccentricity and the
sin i
in the denominators of the equations may become significant trouble with very small value. The Many efforts have been made to avoid this inconvenience by selecting new set of orbital elements. Some elements are functions of the classical elements,
e,
ω, and Ω, others are functions of
i
and Ω. Note that the elements f and g in Eqs. (11) and (12) have been discussed in Moulton (1914), Brouwer & Clemence (1961), Cohen & Hubbard (1962), and Walker et al. (1985) but that
h
and
k
differ slightly by authors. Arsenault et al. (1970) introduced the name ‘equinoctial elements’ for this type of nonsingular elements. In this study, a set of equinoctial elements introduced by Walker et al. (1985) shown in Eqs. (1015) were used for the simulation.
The transformation from the equinoctial elements to the COE is given in Eqs. (1621).
With the calculation of Lagrange bracket with the EOE and the differentiation of Eqs. (1015) with respect to time, the derivatives of EOE can be found in terms of COE and their derivatives. Eqs. (16) can be replaced by the new set of LPE, Eqs. (2227).
where
The disturbing function, Eq. (7) can be rewritten with the substitution of Eqs. (30) and (31).
New set of derivatives of the disturbing functions in terms of EOE is shown in Eqs. (3237).
where
The description of the orbit evolution by the
J_{2}
perturbation with EOE is completed.
4. DIRECT NUMERICAL INTEGRATION OF THE EQUATION OF MOTION IN CARTESIAN COORDINATE
A state vector and a force vector in Cartesian coordinate can be defined as in Eqs. (39) and (40)
where
The acceleration due to the
J_{2}
zonal harmonic nonspherical gravitational perturbation also can be applied to the force vector as in Eq. (40). Point mass twobody gravitational acceleration and
J_{2}
disturbing acceleration term both can be accumulated to form a total acceleration for the satellite.
where,
J_{2}
oblateness coefficient of the Earth
R_{e}
radius of the Earth
The equation of motion expressed in the Cartesian coordinate has an advantage in the numerical integration for its simplicity over analytic or semianalytic type. And the perturbation functions and twobody equations do not have any singular point except
r
= 0 (unusual condition with the finite size of the Earth). The straight forward formulation helps the design and the manipulation of external forces.
The transformation from the state vector to the COE can be achieved without any trouble in most elliptic case. The full procedure for this transformation can be found at Vallado (2004) and other references. Unlike semianalytical solutions, this DNI method may introduce an error as time increases: truncation error. This type of error cannot be avoidable as the computer has limited digit of significant number and the numerical integration method inherits the approximation errors.
5. NUMERICAL PROCEDURES
In this study, the comparison of the results from the each solving method of the equation of satellite motion with
J_{2}
perturbation is the primary objective. It is rather not possible to find any analytical exact solution to compare the each orbit. Remember the LPE type solution is a first order approximation by the variation of parameter.
First, four specific orbits chosen to represent particular orbit might introduce some inconvenience. The case I in the
Table 1
represents an orbit with very small eccentricity and small
sin i
terms; unusual equatorial low Earth orbit (LEO) satellite. The cases II and III show the typical condition of LEO satellite. The case IV displays an example of geodetic satellite. The COE initial conditions in the
Table 1
can be transformed into Cartesian coordinates as in Table A1 in Appendix.
Second, the accuracy of the numerical integration of the equation of the motion with twobody force only in Cartesian coordinate was compared with the analytical
Initial condition (Keplerian elements).h0: altitude at perigee;Earth Radius: 6,378.137 km (Vallado 2004).
Initial condition (Keplerian elements). h_{0}: altitude at perigee; Earth Radius: 6,378.137 km (Vallado 2004).
solution: mean motion on the orbital plane. To validate the DNI method, all cases of initial condition were tested. The value of
a, e, i,ω, M_{0},Ω
should not be changed as time advances in both solutions from analytic method and DNI method in twobody problem without any perturbation. However, in case of DNI method previously mentioned several errors will be accumulated on the state vector, as time advances or rather as we perform every step of numerical integration.
The difference (error) evolutions of orbital elements between the DNI method and the analytical solution are shown in
Figs. 1

7
. In case of DNI method, the integration time step must be chosen with careful consideration for classical RungKutta 4th order method to minimize calculation load, roundoff error, and truncation error. However, it is not a simple task to find an optimal integration step size for general cases.
In this study, fixed step sizes: 60 and 10 seconds were tested for all cases. As shown in
Figs. 1
and
2
, the accuracy was improved in 1 to 1,000 ratio as the time step was reduced in 6 to 1 ratio. With the integration step size of 60 seconds, error in the semimajor axis grows to about 200 meters for the case IV in 7 days. With the integration step size of 10 seconds, same orbit has only 0.025 meter error after 7 days. An optimal integration step size can be found with further test or an integration step size control can be applied to this study. However, 10 second step size was chosen for its obvious simplicity in decimal number with relatively acceptable accuracy.
The error in semimajor axis grows almost linearly as time advances for all four cases as shown in
Fig. 2
. The shorter period the orbit has (i.e. the orbit is closer to the
The error evolution of semimajor axis from the state vector direct numerical integration method with integration step size: 60 seconds compared to the twobody analytical solution for seven days.
The error evolution of semimajor axis from the state vector direct numerical integration method with step size: 10 seconds compared to the twobody analytical solution for seven days.
The error evolution of eccentricity from the state vector direct numerical integration method compared to the twobody analytical solution for seven days.
The error evolution of inclination from the state vector direct numerical integration method compared to the twobody analytical solution for seven days.
The error evolution of mean anomaly at epoch from the state vector direct numerical integration method compared to the twobody analytical solution for seven days
The error evolution of the argument of perigee from the state vector direct numerical integration method compared to the twobody analytical solution for seven days.
The error evolution of the right ascension of ascending node from the state vector direct numerical integration method compared to the twobody analytical solution for seven days.
Earth), the steeper slope the error shows. Since we used a fixed size integration step for all numerical integrations, the orbit with larger orbital period requires larger number of integration step to complete one revolution. It has benefits of smoother and dense derivatives, and shorter interval relative to the whole arc of orbit. Unlike other cases, the case IV shows the sinusoidal behavior with same period of orbital revolution
More complicated figures are shown in the error evolution in eccentricity for all cases as you can see in
Fig. 3
. The sinusoidal variation of each case shows same periodicity of orbital revolution. Over all tendencies are similar to those in semimajor axis case. However, the eccentricity looks more susceptible to the errors induced by the numerical integration.
The error evolutions of inclination and the right ascension of ascending node are shown in
Figs. 4
and
7
respectively. In these two cases, it is hard to find any periodicity or obvious tendency except gradual degradation. The magnitude of the errors in both cases lies close to the lower bound of significant number of digit of the computer used in this study compare to the original values of both orbital elements.
The fast variable in the LPE, mean anomaly shows typical error evolution figures in
Fig. 5
. As shown in
Fig. 6
, the effect of error on the argument of perigee is stronger as shallower the ellipse is.
In summary, the deviations of the DNI method from the exact solution for twobody problem do not exceed 10
^{6}
ratio in case of
e
, 10
^{7}
ratio in case of M
_{0}
, ω, 10
^{9}
ratio in case of a, 10
^{13}
ratio in case of
i
, Ω for 7 days. For the comparison of LPE with two different sets of orbital elements and state variable DNI, we tested the DNI with twobody exact solution first. In conclusion, state variable with
J_{2}
perturbation DNI can be used as a reference to compare two general perturbation method using COE and EOE with the confidence level shown in
Figs. 2

7
.
Third, the general perturbation solutions by using COE and EOE shown in Appendix A and B are calculated and the set of state variable with
J_{2}
perturbation was integrated by a classical RungeKutta method with 10 second integration time step. Though the characteristics of the twobody DNI is shown by comparing it with the analytic solution in the previous section, the twobody problem with
J_{2}
perturbation cannot be expected to act in same manner. However, the
J_{2}
, a dominant perturbation term by mass displacement from the perfect sphere is still 1:1,000 ratio to point mass two?body gravitational force.
6. COMPARISON RESULTS AND DISCUSSION
The time evolution of each orbital element (whether it is integral or variable) by three different approach for solution was investigated with the initial conditions suggested in the previous section. The initial condition of orbits to be investigated was chosen to check the system instability caused by the near singular condition in the denominators of the equations.
The first comparison test with 13,000 seconds duration, which represents one revolution of the longest period of test orbit (case IV), is expected to show major characteristics of short period perturbation as
J_{2}
term only produce short period and secular perturbation (Kozai 1959, Brouwer & Clemence 1961). The evolution of 6 COE and 6 EOE compared to the DNI results for the case II are presented in
Fig. 8
. Most of elements agree to those from other methods but eccentricity. The calculation of increment in eccentricity with COE LPE leads to unstable behavior in case of relatively small eccentricity (~10
^{3}
). The saw tooth shaped trace of mean anomaly is mainly due to the mean motion (
Fig. 8
).
The difference evolution in semimajor axes of two semianalytic solutions from the DNI result is shown in
Fig. 9
. The gaps at the epoch (t = 0) in
Fig. 9
were produced by the offset due to an addition of first solution from the LPE with COE. Even application of mean value to adjust osculating orbital elements to mean orbital elements (Not mean orbit) cannot get rid of these gaps completely. The arbitrary initial orbit elements as shown in
Table 1
are used for orbit simulation in this study, however, we obtain an osculating orbital element from the calculation of the Lagrange planetary equation’s solution with perturbing forces. The result from LPE with perturbing forces is always an osculating orbit. In other words, using a priori orbit elements, perturbed orbit must not be admitted to use at next epoch. This can be adopted both COE and EOE case. In the DNI case, a numerical integration of dynamics including perturbations produces an osculating orbit at any epoch naturally.
As we have mentioned, LPE uses input elements as a reference orbit and calculates each element’s variation caused by perturbations. In this study, orbit simulation of each case has been set a same initial condition, so same orbit elements can be used as an initial value for input at each case, COE, EOE, and DNI. The difference evolution of right ascension of ascending node for all cases is shown in
Fig. 10
. The inclination and the ascending node decide the spatial position of an orbital plane in an inertial frame. The trends of the evolutions of the inclination and the right ascension of ascending node are very similar naturally. As shown in
Fig. 11
a, the amplitude of error variation in argument of perigee is relatively larger than any other angle related orbital elements. However, the error of mean anomaly in the case I by COE LPE shows unreasonable behavior compared to the results from other method. For the analysis of this type of error in the case I,
The evolutions of orbital elements from two semianalytical solutions and direct numerical integration (DNI) method with the J_{2} perturbation for 13000 seconds in case II. From left upper corner (a) semimajor axis (b) eccentricity (c) inclination (d) right ascension of ascending node (e) mean anomaly and (f) argument of perigee. Red line by COE LPEDNI blue line by EOE LPEDNI.
the initial condition of the case V was selected. The initial condition of both cases has identical components except eccentricity. The eccentricity of the case V is exactly ten times of that of the case I. As the eccentricity in the case I was 0.0005, both terms in the right hand side of Eq. (2) induced some inconveniency. On the contrary, the one in the case V with the value of 0.005 did not introduce same magnitude of troubles. More rigorous analysis on the evolution of the eccentricity will be performed on a separate study.
The second comparison test was performed in same manner with one day duration to check the consistencies.Also, we can check the slight hint of long term error behavior of each orbital element from different solutions.
The differences in semimajor axes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases. Red line by COE LPEDNI blue line by EOE LPEDNI.
With the seven day duration of time, only few details of variation trend can survive.
In
Fig. 12
,you can find the various types of error evolutions of semimajor axis. With only J2 perturbation, long period perturbation behavior was not expected. The orbital trajectory from each method (LPE with COE, LPE with EOE, and DNI) did not show obvious long periodic perturbed motion in seven days. However, the cases I, II, and V show periodic variation in errors semimajor axis as shown in
Figs. 12
a, b, and e. Also, you can see similar variation in right ascension of ascending node as in
Figs 13
ac, and e. In case of argument of perigee, same kind of periodic variations can be seen in
Figs. 14
b and c. But it cannot be confirmed as a long period perturbation with
The differences in right ascension of ascending nodes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases. Red line by COE LPEDNI blue line by EOE LPEDNI.
out further analysis.
Then the third comparison test was performed in same manner with seven day duration to check the long term variation in each element. With the seven day duration of time, only few details of variation trend can survive. However, the cases I and V show periodic variation in error evolutions in semimajor axis and right ascension of ascending node as shown in
Figs. 15
a and e,
16
a and e. Also, you can see secular like variation in semimajor axis as in
Fig.15
d. The divergence trend can be seen in the error evolutions in semimajor axis as in
Figs. 15
b and c, in right ascension of ascending node as in
Figs. 16
bd, in argument of perigee as in
Figs. 17
bd.
In summary, the description of a perturbed motion
The differences in arguments of perigees of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases. Red line by COE LPEDNI blue line by EOE LPEDNI.
of satellite by LPE with COE can be erratic on several orbital elements in specific cases like small eccentricity or small inclination. As you can see in most figures, EOE LPE agrees relatively better with DNI method than COE LPE does in most elements. It is consistent with different sizes (semimajor axis), shapes (eccentricity), and orbital plane position (inclination and right ascension of ascending node).
In this paper, we selected only three elements, semimajor axis, right ascension of ascending node, and argument of perigee to present the error evolution by time. However, the position of a satellite on an orbit can be defined by mean anomaly. This non geometric angle plays an important role as a running variable instead of time
The differences in semimajor axes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in one day. Red line by COE LPEDNI blue line by EOE LPEDNI.
in LPE with COE. In this study, the trends of variation in mean anomaly for each case shows similar pattern as in argument of perigee. In the same manner, (semimajor axis, eccentricity), (right ascension of ascending node, inclination), and (argument of perigee, mean anomaly) can be loosely couple by its similarity in the variations.
The theories and equations we used in this study are consistent with Kozai (1959), Brouwer & Clemence (1961) and Vallado (2004). And the result from EOE LPE and DNI agreed each other. Not until further study and analysis are performed, the long periodic variations in the solutions of COE LPE cannot be confirmed.
The differences in right ascensions of ascending nodes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in one day. Red line by COE LPEDNI blue line by EOE LPEDNI.
7. CONCLUSIONS
In this study, we compared the three methods to describe the motion of a satellite with the
J_{2}
perturbation up to seven days. The time evolutions of orbital elements are presented to check the system stability with the several near singular conditions of the equations. As a result, Lagrange planetary equation with equinoctial orbital element agreed better with direct numerical integration method than Lagrange planetary equation with classical Kepler orbital elements in our test cases.
The differences in arguments of perigees of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in one day. Red line by COE LPEDNI blue line by EOE LPEDNI.
The differences in semimajor axes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in seven days. Red line by COE LPEDNI blue line by EOE LPEDNI.
The differences in right ascensions of ascending nodes of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in seven days. Red line by COE LPEDNI blue line by EOE LPEDNI.
The differences in arguments of perigees of two semianalytical solutions from direct numerical integration (DNI) result with the J_{2} perturbation for all cases in seven days. Red line by COE LPEDNI blue line by EOE LPEDNI.
 APPENDIX
A. Short periodic perturbation on orbital elements by
J_{2}
Initial condition in Cartesian coordinate.
Initial condition in Cartesian coordinate.
B. Secular perturbation on orbital elements by the
J_{2}
Arsenault JL
,
Ford KC
,
Koskela PE
1970
Orbit determination using analytic partial derivatives of perturbed motion
AIAAJ
8
4 
12
Blitzer L
1975
Handbook of orbital perturbation
University of Arizona Press
Tempe
28 
32
Broucke RA
,
Cefola PJ
1972
On the equinoctial orbit elements
CeMec
5
303 
310
Brouwer D
,
Clemence GM
1961
Methods of celestial mechanics
Academic Press
New York
273 
307
Cohen CJ
,
Hubbard EC
1962
A nonsingular set of orbit elements
AJ
67
10 
Jo JH
2008
The comparison of numerical integration methods for the KASIOPEA Part II
Bull Korean Space Sci Soc
17
26 
27
Kozai Y
1959
The motion of a close earth satellite
AJ
64
367 
Luthcke SB
,
Zelensky NP
,
Rowlands DD
,
Lemoine FG
,
Williams TA
2003
The 1centimeter Orbit: Jason1 precision orbit determination using GPS SLR DORIS and altimeter data
Mar Geodes
26
399 
421
Moulton FR
1914
An introduction to celestial mechanics
MacMillan Company
New York
387 
425
Nerem RS
,
et al.
,
Chao BF
,
et al.
,
Au AY
,
et al.
,
Chan JC
,
et al.
,
Klosko SM
,
et al.
1993
Temporal variations of the Earth's gravitational field from satellite laser ranging to Lageos
GeoRL
et al.
20
595 
598
Taff LG
1985
Celestial mechanics: a computational guide for the practitioner
Wiley
New York
288 
363
Vallado DA
2004
Fundamentals of astrodynamics and applications
2nd ed
Microcosm Press
El Segundo CA
491 
674
Walker MJH
,
Owens J
,
Ireland B
1985
A set modified equinoctial orbit elements
CeMec
36
409 
419