Two semi-analytic solutions for a perturbed two-body 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, non-singular 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.
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 two-body problem. However, the motion of the satellite in real environment is also dictated by the gravitational and non-gravitational 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 semi-analytical 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 semi-analytic 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 two-body 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 two-body 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. (1-6).
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
zonal non-spherical harmonic gravitational perturbation are analyzed. For Earth gravity perturbations, the disturbing function is merely the negative of the disturbing part of geo-potential (Blitzer 1975). The disturbing function can be described as in Eq. (7).
eeccentricity of orbit
iinclination of orbital plan
ωargument of perigee
However, the disturbing function,
is not a function of COE. The partial derivatives of R in Eqs. (1-6) should survive only when
a, e, i,
ω, M, Ω). The
can be substituted with the functions of
a, e, i, υ, ω
as in Eqs. (8) and (9).
μ gravitational constant multiplied by
the mass of the Earth
rradius of satellite
Φlatitude of the sub-satellite point
R2radius of the Earth
P2Legendre polynomial degree 2
υ true anomaly
With the perturbation by the
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. (1-6). 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. (2-6), small eccentricities or small inclinations can make the equation unstable. The eccentricity and the
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,
ω, and Ω, others are functions of
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
differ slightly by authors. Arsenault et al. (1970) introduced the name ‘equinoctial elements’ for this type of non-singular elements. In this study, a set of equinoctial elements introduced by Walker et al. (1985) shown in Eqs. (10-15) were used for the simulation.
The transformation from the equinoctial elements to the COE is given in Eqs. (16-21).
With the calculation of Lagrange bracket with the EOE and the differentiation of Eqs. (10-15) with respect to time, the derivatives of EOE can be found in terms of COE and their derivatives. Eqs. (1-6) can be replaced by the new set of LPE, Eqs. (22-27).
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. (32-37).
The description of the orbit evolution by the
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)
The acceleration due to the
zonal harmonic non-spherical gravitational perturbation also can be applied to the force vector as in Eq. (40). Point mass two-body gravitational acceleration and
disturbing acceleration term both can be accumulated to form a total acceleration for the satellite.
oblateness coefficient of the Earth
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 semi-analytic type. And the perturbation functions and two-body equations do not have any singular point except
= 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 semi-analytical 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
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
represents an orbit with very small eccentricity and small
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
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 two-body 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). h0: 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,ω, M0,Ω
should not be changed as time advances in both solutions from analytic method and DNI method in two-body 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
. In case of DNI method, the integration time step must be chosen with careful consideration for classical Rung-Kutta 4th order method to minimize calculation load, round-off 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
, 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 semi-major 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 semi-major axis grows almost linearly as time advances for all four cases as shown in
. The shorter period the orbit has (i.e. the orbit is closer to the
The error evolution of semi-major axis from the state vector direct numerical integration method with integration step size: 60 seconds compared to the two-body analytical solution for seven days.
The error evolution of semi-major axis from the state vector direct numerical integration method with step size: 10 seconds compared to the two-body analytical solution for seven days.
The error evolution of eccentricity from the state vector direct numerical integration method compared to the two-body analytical solution for seven days.
The error evolution of inclination from the state vector direct numerical integration method compared to the two-body analytical solution for seven days.
The error evolution of mean anomaly at epoch from the state vector direct numerical integration method compared to the two-body analytical solution for seven days
The error evolution of the argument of perigee from the state vector direct numerical integration method compared to the two-body 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 two-body 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
. The sinusoidal variation of each case shows same periodicity of orbital revolution. Over all tendencies are similar to those in semi-major 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
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
. As shown in
, 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 two-body problem do not exceed 10
ratio in case of
ratio in case of M
, ω, 10
ratio in case of a, 10
ratio in case of
, Ω for 7 days. For the comparison of LPE with two different sets of orbital elements and state variable DNI, we tested the DNI with two-body exact solution first. In conclusion, state variable with
perturbation DNI can be used as a reference to compare two general perturbation method using COE and EOE with the confidence level shown in
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
perturbation was integrated by a classical Runge-Kutta method with 10 second integration time step. Though the characteristics of the two-body DNI is shown by comparing it with the analytic solution in the previous section, the two-body problem with
perturbation cannot be expected to act in same manner. However, the
, 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
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
. 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
). The saw tooth shaped trace of mean anomaly is mainly due to the mean motion (
The difference evolution in semi-major axes of two semi-analytic solutions from the DNI result is shown in
. The gaps at the epoch (t = 0) in
were produced by the off-set 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
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
. 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
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 semi-analytical solutions and direct numerical integration (DNI) method with the J2 perturbation for 13000 seconds in case II. From left upper corner (a) semi-major axis (b) eccentricity (c) inclination (d) right ascension of ascending node (e) mean anomaly and (f) argument of perigee. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
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 semi-major axes of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
With the seven day duration of time, only few details of variation trend can survive.
,you can find the various types of error evolutions of semi-major 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 semi-major axis as shown in
a, b, and e. Also, you can see similar variation in right ascension of ascending node as in
a-c, and e. In case of argument of perigee, same kind of periodic variations can be seen in
b and c. But it cannot be confirmed as a long period perturbation with-
The differences in right ascension of ascending nodes of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
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 semi-major axis and right ascension of ascending node as shown in
a and e,
a and e. Also, you can see secular like variation in semi-major axis as in
d. The divergence trend can be seen in the error evolutions in semi-major axis as in
b and c, in right ascension of ascending node as in
b-d, in argument of perigee as in
In summary, the description of a perturbed motion
The differences in arguments of perigees of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
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 (semi-major axis), shapes (eccentricity), and orbital plane position (inclination and right ascension of ascending node).
In this paper, we selected only three elements, semi-major 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 semi-major axes of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in one day. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
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, (semi-major 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 semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in one day. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
In this study, we compared the three methods to describe the motion of a satellite with the
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 semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in one day. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
The differences in semi-major axes of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in seven days. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
The differences in right ascensions of ascending nodes of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in seven days. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
The differences in arguments of perigees of two semi-analytical solutions from direct numerical integration (DNI) result with the J2 perturbation for all cases in seven days. Red line by COE LPE-DNI blue line by EOE LPE-DNI.
A. Short periodic perturbation on orbital elements by
Initial condition in Cartesian coordinate.
Initial condition in Cartesian coordinate.
B. Secular perturbation on orbital elements by the
Orbit determination using analytic partial derivatives of perturbed motion
DOI : 10.2514/3.5597
Handbook of orbital perturbation
University of Arizona Press
Methods of celestial mechanics
The comparison of numerical integration methods for the KASIOPEA Part II
Bull Korean Space Sci Soc
The 1-centimeter Orbit: Jason-1 precision orbit determination using GPS SLR DORIS and altimeter data
An introduction to celestial mechanics
Temporal variations of the Earth's gravitational field from satellite laser ranging to Lageos
DOI : 10.1029/93GL00169
Celestial mechanics: a computational guide for the practitioner
Fundamentals of astrodynamics and applications
El Segundo CA
A set modified equinoctial orbit elements
DOI : 10.1007/BF01227493