Orbit Propagation

Orbit Propagation

Rodney Long

Rodney Long has worked as an engineer and programmer for more than 15 years in industry and government, including work on the NASA Space Shuttle and U.S. Air Force Global Positioning System projects. He may be reached at 19003 Swan Drive, Gaithersburg, MD 20879.

The prediction problem in orbital mechanics is this — given the position and velocity of an orbiting object at some initial time t0, find the position and velocity at some second (usually later) time t. (A time-indexed table of such positions and velocities is called an ephemeris.)

This article discusses the prediction problem of orbital mechanics. It provides sample C code implementations for some of the classical and modern computations that solve that problem. You don't need a course in orbital mechanics to use this code, but you do need some grasp of the concepts and terminology. I have provided the essentials here and recommend Fundamentals of Astrodynamics [1] for further study.

Orbital Predictors

For orbit prediction, exact expressions are available for the solutions of objects moving in two-body or Keplerian motion, which involves modeling the motion of an object being accelerated only by the gravitational force from a second object. In two-body modeling, the two objects are treated as point-masses, which means the entire mass of each object is considered to be concentrated at one mathematical point. This is the simplest model used for orbital motion.

Analytic orbit predictors evaluate these exact expressions and are called orbit propagators. Orbit propagators require relatively few numerical calculations and execute quickly.

When you want to model the orbit motion precisely, you must use a numerical method, because orbiting objects are subject not only to the point-mass gravitational force, which produces their Keplerian motion, but also to a range of other forces that perturb or modify this motion.

Which perturbing forces are significant to model depends on the particular orbiting object and its orbit. For example, for low-altitude Earth satellites, drag forces due to the earth's atmosphere are highly significant. Other forces commonly modeled in space-tracking systems include point-mass effects of the sun and moon, perturbations due to the scattered concentrations of mass within the earth, thrust effects (if the orbiting object has engines), and solar pressure due to the bombardment of the object from particles originating in the sun.

Because no exact mathematical expressions exist to model this complex, perturbed motion, you must use numerical approximations. These techniques are orbit integrators, as distinguished from the analytic orbit propagators. Orbit integrators require many numerical calculations. Also, depending on the forces being modeled, the integrator may need to repeatedly access databases of empirical data, containing such information as ephemeris for the sun and moon, atmosphere models, etc. For these reasons, integrators have computationally-intensive, slower execution.

Integrators do the actual tracking of orbiting objects. Perturbations change the orbital motion too significantly to allow for practical use of a propagator for long-term tracking. Nevertheless, propagators are a valuable method of getting a quick first approximation of an object's path. They provide short-duration predictions for an object not subject to severe perturbing forces, such as a small, unpowered, high-altitude earth satellite. An orbit propagator may also perform a consistency check on an integrator. By setting all the non-Keplerian forces in the integrator to zero, you should get the positions and velocities output from the integrator to match those output by the propagator.

Almost all major space-tracking systems incorporate propagators as part of their systems. For example, the original space-tracking program for the space shuttle used a program called the Analytic Ehpemeris Generator (AEG), which solved the equations of motion by a variation of parameters method. The Goddard Trajectory Determination System (GTDS), used for general-purpose satellite tracking at Goddard Space Flight Center in Greenbelt, MD, contains a propagator called TWOBDY which uses a universal variables analytic algorithm. You will find a C version of TWOBDY at the end of this article.

Two-Body Orbital Mechanics

Under Keplerian motion, the fundamental equation governing the motion of this object is:

RDD = - (mu/r<SUP>3</SUP>)*R

R = (x,y,z) is the position of the object,
RDD = R¨ is the acceleration of the object,
r = the magnitude of R,
mu = g(mr + mo),
g = the universal gravitational constant,
mc = mass of the "central body" (the Earth), and
mo = mass of the orbiting object.

(Uppercase denotes vectors in the above equation.)

In the context of one very small object orbiting a much larger one, the term mu is usually approximated by g*(mc), since mc >> mo.

Objects under Keplerian motion move in paths that are conic sections; this follows from the previous fundamental equation. The central body lies at one focus of the conic section. The point in the orbit when the object most closely approaches the earth is called periapsis or perigee. The point in the orbit when the object is farthest from the earth is called apoapsis or apogee. From the fundamental equation, you may derive an equation that displays the conic nature of the orbits, namely:

r = p/(1 + e * cos(nu))

r is defined as above,
p = the semi-latus rectum of the conic section which the object travels on,
nu = the angle, measured in the orbit plane between (1) a line from focus to periapsis, and (2) a line from focus to the object's current position, and
e = the eccentricity of the conic section.

Nu is called the "true anomaly" of the object. If you know nu and you know where periapsis is, you know the object's position. See Figure 2.

The value of e determines the type of conic section that the orbit describes:

• e < 1 implies that the orbit is an ellipse.
• e = 1 is a parabolic orbit.
• e > 1 implies that the object travels on one branch of a hyperbola.

For elliptical orbits, the object's state (position and velocity) will repeat periodically. For the parabolic and hyperbolic cases, the position will obviously change indefinitely as the object continues to move away from the central body. But what happens to the velocity? The answer is, for the parabola, the velocity approaches zero as the distance from the central body increases. For the hyperbola, however, it is possible for the object to have a finite non-zero velocity regardless of how far it has traveled. Whether the hyperbolic object's velocity approaches zero or some constant non-zero value depends on how much initial velocity it is given when launched into its orbit.

Earth satellites follow elliptical orbits; many of them are circular to a high degree. Some comets follow parabolic paths, but parabolic orbits are otherwise rare in nature [1]. Earth-launched interplanetary missions follow hyperbolic paths.

Keplerian Elements

• a = the semi-major axis of the orbit
• e = the eccentricity of the orbit
• i = the inclination of the orbit plane with respect to the X-Y plane
• lan = the longitude of the ascending node, or the angle, measured in the X-Y plane, between the positive X-axis and the point at which the object crosses the X-Y plane when ascending in the positive Z direction
• w = the argument of periapsis, or the angle, measured in the orbit plane, between the X-Y plane and periapsis
• M = the mean anomaly, defined below

If you know these five elements, you know the two-body orbit completely, but you don't know where your object is in the orbit. The element M provides this information. M determines the orbiting object's location, relative to periapsis.

Knowing the six Cartesian coordinates for position and velocity is equivalent to knowing the six Keplerian elements. Any good orbital mechanics text will explain the conversion of one system to the other (see Reference 1). I have written two C functions, c2ke and k2ce, to perform the conversions in both directions. The C code in Listing 1 shows the calling sequences and sample executions.

Solving Orbit Prediction For Ellipses

M = (mu/a<SUP>3</SUP>)<SUP>1/2</SUP> * [t - tp]

mu and a are defined as before,
t = time corresponding to current position of the object, and
tp = time that the object went through its periapsis point.

If you know the first five Keplerian elements and tp, you can compute M for any desired time t.

The eccentric anomaly E is defined (for elliptical orbits) by circumscribing a circle around the ellipse, then dropping a perpendicular from a point p on the circle so it passes through the current position of the object on the ellipse (see Figure 4) . The angle between the major axis and a line drawn from the center of the circle to p is called the eccentric anomaly. (In Figure 4, F is the central body, O is the orbiting object, and C is the center of the cirsumscribed circle.)

The equation of Kepler that relates the mean anomaly M to the eccentric anomaly E is:

M = E - e*sin(E)

The relationship between E and the true anomaly nu is given by:

cos(E) = (e + cos(nu)) 1 + e*cos(nu))

cos(nu) : (e - cos(E)) e*cos(E) - 1)

Unlike the previous equations, Kepler's equation does not allow you to solve explicitly for E or sin (E). Instead, you must make a first estimate, and then iteratively improve your previous estimate, usually by Newton's method. Listing 3 is a C function sk that solves Kepler's

equation for E. The inputs are the mean anomaly M, the starting estimate E0 for the solution to Kepler's equation, the orbit eccentricity e, and two parameters, mi and t, which determine when the iteration will terminate. If two successive estimates for the eccentric anomaly differ in absolute value by less than t, the procedure stops. Otherwise, it continues for a maximum of mi iterations. Listing 3 also provides a sample execution of sk.

Now, remember that:

M = (mu/a<SUP>3</SUP>)<SUP>1/2</SUP> * [t - tp]

1. Compute the mean anomaly M.

2. Solve Kepler's equation for E.

3. Use the equations for E and nu to compute nu.

4. Use the Keplerian orbit parameters plus nu to compute position and velocity.

(To find the position and velocity at t given the eccentric anomaly E0 at any time t0, with t0 not necessarily the time of periapsis, is only slightly more difficult.)

In spite of being a time-tested method, the classic approach has at least two drawbacks, already implicit in what I have described earlier. The first difficulty is numerical. The classical method requires the iterative computation of the quantity (E - e*sin(E)), until this quantity matches a given M to within a satisfactory tolerance. Beware of subtracting nearly equal quantities. In the case where E is small, sin(E) is approximately E (remember that sin(x) = x - x³/3! + x⁵/5! - ...). Morever, if the orbit is near-parabolic, e is close to 1. Under these conditions, E is very close to e*sin(E), and the subtraction may result in a signficant loss of precision.

The second objection to the classical prediction approach is that it is not general. For elliptical orbits, you solve Kepler's equation, but for parabolas and hyperbolas you must use similar equations involving the "parabolic (hyperbolic) eccentric anomalies" D and F.

For at least 25 years, an approach that gets around both the numerical difficulty and the lack of generality has been available.

The Universal Variables Solution

The universal variables algorithm unifies the separate approaches to orbit prediction for elliptical, parabolic, and hyperbolic orbits by providing one procedure that iteratively solves for a "generalized eccentric anomaly" which plays the role of D, E, and F. Moreover, it does not suffer from the numerical instabilities of the classic method.

This generalized anomaly, denoted psi, is defined implicitly by:

    dt/d(psi) = r

t = the time variable, and
r = orbit radius at t

Through some lengthy mathematical manipulations (shown in the NASA report), you can formulate a generalized Kepler equation in terms of psi:

    t = t0 + r0 * s1 + sig0 * s2 + mu * s3

t = time value corresponding to psi,
t0 = time value corresponding to psi0,
r0 = orbit radius corresponding to psi0,

s1,s2,s3 = functions of psi which may be evaluated by series expansions,
sig0 = the dot product of the position and velocity at time t0, and
mu is as previously defined.

Now suppose you know the position and velocity (r0, v0) at time t0, and you want to know the position and velocity (r, v) at time t. Let tau be defined by tau = t - t0. Then define dtau by:

   dtau(psi) = r0 * s1(psi) + sig0 * s2(psi) +
               mu * s3(psi) - tau

Then, you will be able to compute r and v at time t, as functions of this psi. The TWOBDY algorithm does just this, using the Newton method for adjusting the estimates. When looking for the Newton step in the code, you should know that the derivative of dtau with respect to psi is given by:

   d(dtau)/d(psi) = r0 * s0 + sig0 * s1 + mu * s2

If the Newton step does not bring psi to within acceptable bounds on an iteration, TWOBDY tries several other adjustments to the current psi. If the resulting psi is acceptable, Newton iteration continues. If all of the alternative adjustments fail, iteration is halted, and the assumption is that psi may not be improved further. Reaching a psi such that dtau(psi) = 0 also halts the iteration.

The abstract and advanced concepts used in TWOBDY, and the many intermediate calculations, may obscure the basically simple logic of the modified Newton algorithm. For this reason, I have provided a brief pseudo-code summary in Figure 5.

The C version, twobdy, of the TWOBDY program follows the public domain FORTRAN version fairly closely, although C structures have been added for input and output, and some general program structure has been incorporated. Note that psi is output into the "input" structure, which facilitates using a good initial estimate of psi when making iterative calls to twobdy.

I have submitted to the C Users' Group Library expanded versions of the codes in this article. The volume includes both Keplerian-to-Cartesian and Cartesian-to-Keplerian routines, as well as a version of twobdy that outputs partial derivatives in addition to position, velocity, and psi.

Bibliography

[2] Atkinson, K., An Introduction to Numerical Analysis, John Wiley and Sons, 1978.

[3] Goodyear, W.H., "A General Method for the Computation of Cartesian Coordinates and Partial Derivatives of the Two-Body Problem," NASA Contractor Report, NASA CR-522, September 1966.

Listing 4

Sidebar: Analytic vs. Numerical