Orbital Elements to State Vectors

Racist Celestial Mechanics, Lecture 1

Revised: 9 August 2003

An orbit is the path followed by an unaccelerated object under the influence of external gravitational fields. Throughout this web page, I consider the celestial mechanics of objects orbiting the sun, which is assumed to be the only important source of gravitational attraction. As a result, our problem is the relatively simple problem of two bodies.

All orbits considered here will be confined to planes, and, furthermore, they will all be conic sections: ellipses or hyperbolas. I will not deal with parabolas just yet. The sun will be at a focus of these conic sections, not at their geometical center.

The default coordinate system used on this page is classical heliocentric ecliptic coordinates. The classical HEC system has its origin at the sun, its XY plane being the plane of Earth's orbit, the X axis being positive in the direction of the Vernal Equinox, and the Z axis being positive in the direction of the angular momentum vector of Earth's orbit.

My earlier paper on approximately this subject can be found here.

Six numbers are required to define an orbit's shape, size, orientation in space, and the specific position of the planet (or spaceship) along the path at a particular time. These six numbers are called "elements of the orbit," or simply "orbital elements."

The Semimajor Axis, a.

Traditionally, the first-listed orbital element is the semimajor axis (a).

For elliptical orbits, the semimajor axis is the distance from the center of the ellipse to either of the "long ends" of the ellipse. The "semi" in the phrase means half. An ellipse's major axis is the longest line segment that can be drawn within the ellipse. The semimajor axis is half of the major axis. (We celestial mechanics will let the philosophers figure out which half it is.)

An elliptical orbit is periodic. A planet moving in an elliptical orbit goes around the sun repeatedly along the same path, completing each revolution in the same amount of time. This "period" of the orbit can be calculated by knowing the gravitational constant, G, the sun's mass, M, and the semimajor axis, a.

P = 2p { a³ / GM }^1/2

Where (P) is in seconds, (a) is in meters, and the product of GM has a value of 1.32712440018E+20 m³ sec^-2.

However, if you input the semimajor axis in astronomical units, the conversion factor cancels out the other constants, leaving only a numerical relationship (units suppressed) of

P = (365.256898326 days) a^3/2

To convert a distance in meters to astronomical units, you must know that 1 AU equals 1.49597870691E+11 meters.

The mean motion (average angular speed) is

m = 2p radians / P

The sun-relative speed in an elliptical orbit is found from

V = { GM ( 2/r - 1/a ) }^1/2

where (r) is the current distance from the sun. MKS input values for (r), (a) and GM will give (V) in meters per second.

An elliptical orbit has two apsides: a perihelion and an aphelion. The perihelion is the point on the orbit closest to the sun. The aphelion is the point on the orbit furthest from the sun. The perihelion and aphelion are on opposite sides of the sun. The perihelion and the aphelion are also on opposite sides of the geometrical center of the ellipse, at either end of the major axis.

For hyperbolic orbits, the semimajor axis is the distance from the vertex (the point where the hyperbolic asymptotes cross) and the nearest point on the hyperbolic curve.

A hyperbolic orbit is non-periodic. A spaceship moving in a hyperbolic orbit makes only one pass by the sun before it leaves the solar system forever, unless something happens to make it slow down; e.g., it fires its retro-rockets or crashes into a planet or something of that nature.

The mean motion along a hyperbolic orbit remains mathematically useful, although it is calculated directly from the gravitational parameter, GM, and the semimajor axis.

m = 86400 sec/day { GM / a³ }^1/2

where (a) is in meters and the MKS value is used for for GM, the mean motion from the equation will be in radians per day. (Without the 86400 factor, it would have been radians per second.)

The sun-relative speed in a hyperbolic orbit is found from

V = { GM ( 2/r + 1/a ) }^1/2

A hyperbolic orbit has only one apside: a perihelion. As with elliptical orbits, the perihelion of a hyperbolic orbit is the point on the orbit that is closest to the sun. It just so happens that the perihelion is the point on a hyperbola that is closest to both the vertex and to the focus (where the sun is), however the sun and the vertex are on opposite sides of the perihelion and at unequal distances. The distance from the vertex to the perihelion is, by definition, the semimajor axis of the hyperbolic orbit.

Many conic sections have two foci, not only one. But only one focus will be the location of the sun. The other focus will be physically unimportant. You can write a science-fiction story in which strange things happen there, if you want to.

The semimajor axis is usually expressed in astronomical units (AU). However, it is sometimes necessary to convert it into meters for use in a formula.

The semimajor axis of Earth's orbit is 1.0000001 AU.

The Eccentricity, e.

The second orbital element is the eccentricity (e). The eccentricity is a pure number that determines which type of conic section an orbit will be. If e=0, the orbit is a circle. If e is between zero and one, the the orbit is an ellipse. if e=1, the orbit is a parabola. If e>1, the orbit is a hyperbola.

The distance from the sun-focus to the perihelion of an elliptical orbit is

r_P = a ( 1 - e )

The distance from the sun-focus to the aphelion of an elliptical orbit is

r_A = a ( 1 + e )

The general equation for the distance of a point on an elliptical orbit is

r = a ( 1 - e² ) / ( 1 + e cos q )

The distance from the sun-focus to the perihelion of a hyperbolic orbit is equal to

r_P = a ( e - 1 )

The general equation for the distance of a point on a hyperbolic orbit is

r = a ( e² - 1 ) / ( 1 + e cos q )

The angle q is called the true anomaly. It is the angle, subtended at the sun-focus, between a ray extended through the perihelion and another ray extended through the current position of the planet or spaceship.

Although the equations for heliocentric distance, r, are similar with respect to elliptic and hyperbolic orbits, the fact that the eccentricities are on either side of one makes a big difference in how the curves graph.

Ellipse eccentricities being always less than one means that the (e cos q) term in the denominator never gets negatively substantial enough to cause the whole denominator to become zero. The equation of the orbit thus never "blows up" from a divide-by-zero condition.

On the other hand, hyperbolic eccentricities are larger than one, which means that the (e cos q) term can cause precisely that kind of trouble, when cos q = -1/e. The whole denominator becomes zero, and consequently the heliocentric distance becomes infinite. The two values for q that causes this to occur correspond to lines parallel to each of the asymptotes of the hyperbola.

Since objects moving at finite speed will never reach the other side of infinity, hyperbolic orbits are thus non-periodic. Even a nigger could probably figure that out.

The eccentricity of Earth's orbit is 0.016708.

The Inclination, i.

The semimajor axis and eccentricity of an orbit define its size and shape. The next three elements of the orbit are angles that define the orbit's physical orientation in space. The first of these angles is the inclination to the ecliptic (i).

The inclination is an angle between two planes, that of the orbit and that of the ecliptic, and it can be measured at any point along the line of intersection. This line of intersection is called "the line of nodes" in celestial mechanics. The inclination can take any value from zero to p radians.

The ecliptic plane is classically taken to be the plane of Earth's orbit. Both the sun and the Earth are always in the ecliptic plane, so defined. The sun, furthermore, is taken as the origin of the heliocentric ecliptic coordinate (HEC) system. The "line of nodes" for any heliocentric orbit passes through the sun.

On this page, as well as other pages that I write on this subject, the inclination of Earth's orbit is defined as zero.

The Longitude of the Ascending Node, W.

The longitude of the ascending node (W) is the angle, subtended at the sun, measured in the ecliptic plane in the direction of increasing heliocentric longitude, from the Vernal Equinox to the point where the planet (or spaceship) crosses the ecliptic plane going from ecliptic south to ecliptic north.

Heliocentric longitude works a lot like the more familiar terrestrial longitude does. The universe is imagined to be compressed on to the surface of a sphere of unit radius having the sun at its center. On this imaginary sphere, lines are drawn to represent the "equator" (which is the ecliptic, the plane of Earth's orbit), as well as parallels of latitude and meridians of longitude.

The Vernal Equinox is a direction in space heading somewhere off in the direction of Pisces. There aren't any bright stars in that direction, unfortunately, so I have no fixed marker to use in reference to it. The Vernal Equinox is along the line of intersection between Earth's orbit (i.e., the ecliptic) and the plane in which Earth's equator lies. The season of spring officially begins at the moment the sun, as seen from Earth, moves into position along the Vernal Equinox.

The "north" side of the ecliptic plane is the side from which, looking back down on the solar system, you'd see Earth (and the other major planets) moving in a counterclockwise direction. When an object crosses the ecliptic heading "north," we picture it as "rising above" the ecliptic plane, and so this point on the line of nodes is called the ascending node.

The longitude of the ascending node for Earth's orbit is defined as zero.

The Argument of the Perihelion, w.

The argument of the perihelion (w) is the angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the ascending node to the perihelion.

The argument of the perihelion of Earth's orbit is 1.797879352 radians.

Time of Perihelion Passage, T.

The time of perihelion passage is, obviously, the moment at which the planet happens to be at the perihelion of its own orbit. Celestial mechanics like to represent T as decimal days (e.g., 4.75 days instead of 4 days 18 hours) in order to avoid the necessity of converting years, months, hours, minutes, and seconds into decimal day form. The classical day-numbering system is called Julian Date. The first day of the Julian Date system began in ancient times because some medieval idiot thought that the universe was "created" on that day.

Just for reference, midnight beginning the day of 1 January 2000 at Greenwich, England, corresponds to Julian Date 2451544.500, exactly.

Earth was at the perihelion of its orbit on JD 2451547.5, which corresponds on the calendar to midnight beginning the day of 4 January 2000.

Given the Elements of an Elliptical Orbit
Solve for a HEC State Vector.

Finding the Mean Anomaly.

We first determine the number of days the planet has been moving since it was at a (preferably recent) perihelion.

dt = t - T

where (t) is the moment of interest. The period, P, of the orbit is

P = (365.256898326 days) a^3/2

We next calculate the mean motion, m of the planet along its orbit. The mean motion is the angular speed at which the planet would move constantly if its orbit were a circle having radius equal to the semimajor axis of the actual, elliptical orbit.

m = 2p / P

We determine the amount by which the mean anomaly has changed since the planet was at this particular perihelion.

DM = m dt

The mean anomaly is the angle, subtended at the sun, measured in the plane of the orbit in the direction of motion, from the perihelion to where the planet would have been if it moved constantly around the sun at its average angular speed.

The current mean anomaly is

M = M_perihelion + DM

Since the mean anomaly is defined such that it's value is zero at perihelion, M=DM. However, we correct the mean anomaly to the interval [0,2p) by adding or subtracting the appropriate multiple of 2p radians.

Now that you've seen it broken into a series of absurdly easy steps, I'll combine the equations for ease of use.

M = 2p ( t - T ) / (365.256898326 days a^3/2)

Input days for (t) and (T), and use astronomical units for (a). Correct (M) to the interval [0,2p), if necessary.

Finding the Elliptical Eccentric Anomaly.

The equation relating the mean anomaly to the eccentric anomaly is called "Kepler's equation" because Johannes Kepler is the one who discovered it.

M = u - e sin u

Where (u) is the eccentric anomaly, which must be in radians. Notice that if we already knew the value of (u) and (e) and wanted to get the value of (M), the problem would be simple. However, we wish to go in the opposite direction: given (M) and (e), find (u). This is a much more difficult problem, and it cannot be solved with algebra alone.

Because the equation we must solve is transcendental in the variable we're trying to solve for, we implement some sort of differential calculus method for converging to the actual value of the eccentric anomaly. The simplest method of this type is Newton's method.

However, Danby's Method is superior because it is less prone to skittering, which is a failure to converge that sometimes occurs with Newton's Method for high eccentricity orbits.

The converged value for the eccentric anomaly (the final u_i+1) is the one assigned to the quantity (u) hereafter.

Finding the Canonical Position Vector.

The canonical position vector is the heliocentric position vector (from the sun to the object) referred to the coordinate system having the sun at the origin, with the XY plane being that of the orbit, with the X axis running positively through the perihelion, and with the Z axis being positive in the direction of the planet's angular momentum.

I present the canonical position vector as the triple-primed vector:

Rotating the Position to Heliocentric Ecliptic Coordinates (HEC).

We firstly adjust the components of the canonical (triple-primed) position vector by a rotation of the coordinate system around the Z''' negatively by the argument of the perihelion, w.

We secondly adjust the resulting vector by a rotation of the double-primed coordinate system around the X'' axis negatively by the inclination, i.

We thirdly adjust the resulting vector by a rotation of the single-primed coordinate system around the Z' axis negatively by the longitude of the ascending node, W.

The position vector [ x , y , z ] is the heliocentric position of the spaceship (or planet) at the moment of interest (t).

Finding the Canonical Velocity Vector.

The canonical velocity vector is the velocity of the object referred to the coordinate system that I defined when I was finding the canonical position vector (q.v.).

We obtain the true anomaly of the object's current position from the components of the triple-primed position vector.

q = arctan2( y ''' , x ''' )

The arctan2 function is the two-argument version of the arctangent, giving an angle corrected for the quadrant indicated by the input sine (the first argument) and cosine (the second argument).

I present the canonical velocity vector as the triple-primed vector:

Rotating the Velocity to Heliocentric Ecliptic Coordinates (HEC).

The canonical velocity vector is adjusted by the same coordinate rotations used (above) to adjust the position to HEC.

We secondly adjust the resulting vector by a rotation of the double-primed coordinate system around the X'' axis negatively by the inclination, i.

We thirdly adjust the resulting vector by a rotation of the single-primed coordinate system around the Z' axis negatively by the longitude of the ascending node, W.

The velocity vector [ V_X , V_Y , V_Z ] is the sun-relative velocity of the spaceship (or planet) at the moment of interest (t). A check on the magnitude of this velocity is available from the Vis Viva equation:

V_X² + V_Y² + V_Z² = GM { 2 / ( x² + y² + z² )^1/2 - 1/a }

And, again, ensure that the units within the same equation are consistently to the same scale.

The state vector of the object in an elliptical orbit of elements [ a , e , i, W , w , T ] at the moment of interest (t) is

[ x , y , z , V_X , V_Y , V_Z ]

Worked Example for an Elliptical Orbit

We are given these elements and time of interest:

We find the mean anomaly.

M = 5.693069656 radians

We use the mean anomaly as our initial approximation to the eccentric anomaly. Then,

The eccentric anomaly, u, is 5.089077456 radians.

We get the canonical position vector.

We rotate the triple-primed coordinate system around its Z axis negatively by the argument of the perihelion, getting

We rotate the double-primed coordinate system around its X axis negatively by the inclination, getting

The Z component of the single-primed position vector is only about five-and-a-half kilometers above the ecliptic. Hmm. We may have caught the object near one of its nodes! We next rotate the single-primed coordinate system around its Z axis negatively by the longitude of the ascending node, getting the position in heliocentric celestial coordinates.

The true anomaly is

q = 4.333250151 radians

The canonical velocity vector is

We rotate the coordinates for the sun-relative velocity the same way that we did for the heliocentric position.

The state vector for this object at this point in its orbit is

Given the Elements of a Hyperbolic Orbit
Solve for a HEC State Vector.

Finding the Mean Anomaly.

We first determine the number of days the spaceship has been moving since it was at it's perihelion.

dt = t - T

where (t) is the moment of interest. If the spaceship has not yet reached perihelion, then dt will be negative.

There is, of course, no period associated with a hyperbolic orbit. But we can determine an equivalent to the mean motion.

m = { GM / a³ }^1/2

Remember to convert (a) from astronomical units to meters before solving for the mean motion or for the mean anomaly. One astronomical unit is 1.49597870691E+11 meters. The value of the gravitational parameter GM is 1.32712440018E+20 m³ sec^-2.

Likewise, the mean anomaly has no special geometric meaning for a hyperbolic orbit, but it nonetheless remains mathematically convenient as an intermediate quantity. The mean anomaly is zero at the perihelion, negative prior to perihelion, and positive after perihelion.

M = m dt

I'll combine the equations for the mean anomaly, skipping an explicit reference to the mean motion.

M = { GM / a³ }^1/2 ( t - T )

In the equation, just above, use meters for (a) and input seconds for (t) and (T). Alternatively,

M = { 0.01720209895 AU^3/2 sec^-1 } ( t - T ) / a^3/2

Where, in this equation, (a) is entered as astronomical units and the times are entered in days. Important: We do not correct the mean anomaly of hyperbolic orbits to the interval [0,2p). If it comes out negative, we leave it that way.

Finding the Hyperbolic Eccentric Anomaly.

Kepler's equation for hyperbolic orbits is

M = e sinh u - u

Where (u) is the hyperbolic eccentric anomaly, which must be in radians. As it was with the elliptical case, the equation is transcendental in the variable that we are trying to find, and we must use a differential calculus method for solving it.

We will use Danby's Method for hyperbolic orbits.

You could use u₀ = M as the initial guess in Danby's Method, and it would probably still converge.

The converged value for the eccentric anomaly (the final u_i+1) is the one assigned to the quantity (u) hereafter. We do not correct the hyperbolic eccentric anomaly to the interval [0,2p). It should have the same sign that the mean anomaly has.

Finding the Canonical Position Vector.

The canonical position vector is the heliocentric position vector (from the sun to the object) referred to the coordinate system having the sun at the origin, with the XY plane being that of the orbit, with the X axis running positively through the perihelion, and with the Z axis being positive in the direction of the planet's angular momentum.

We determine the true anomaly, q, from the eccentricity and the eccentric anomaly.

The heliocentric distance is

r = a (e cosh u - 1)

I present the canonical position vector as the triple-primed vector:

Rotating the Position to Heliocentric Ecliptic Coordinates (HEC).

The coordinate rotations to move the position vector from canonical coordinates [ x''' , y''' , z''' ] to HEC [ x , y , z ] is done for hyperbolic orbits in exactly the same way it is for elliptical orbits. See the section above where I show how to do it.

Finding the Canonical Velocity Vector.

The canonical velocity vector is the velocity of the object referred to the coordinate system that I defined when I was finding the canonical position vector (q.v.).

I present the canonical velocity vector as the triple-primed vector:

Rotating the Velocity to Heliocentric Ecliptic Coordinates (HEC).

The coordinate rotations to move the velocity vector from canonical coordinates [ V_X''' , V_Y''' , V_Z''' ] to HEC [ V_X , V_Y , V_X ] is done for hyperbolic orbits in exactly the same way it is for elliptical orbits. See the section above where I show how to do it.

A check on the magnitude of the velocity is available from the Vis Viva equation:

V_X² + V_Y² + V_Z² = GM { 2 / ( x² + y² + z² )^1/2 + 1/a }

Notice that there is a sign change in the Vis Viva equation for hyperbolic orbits, as compared with the Vis Viva equation for elliptic orbits. That's because I considered the semimajor axis to be positive. And, again, ensure that the units within the same equation are consistently to the same scale by converting the distances to meters.

The state vector of the object in a hyperbolic orbit of elements [ a , e , i, W , w , T ] at the moment of interest (t) is

[ x , y , z , V_X , V_Y , V_Z ]

Worked Example for a Hyperbolic Orbit

We are given these elements and time of interest:

We find the mean anomaly.

M = -8.714915420

We use the mean anomaly as our initial approximation to the eccentric anomaly. Then,

The eccentric anomaly, u, is -1.299202502.

We determine the true anomaly.

q ' = +1.191649715 radians

q = +5.091535592 radians

The heliocentric distance is

r = 2.178398513 AU

We get the canonical position vector.

We rotate the triple-primed coordinate system around its Z axis negatively by the argument of the perihelion. However, since the argument of the perihelion is zero, the vector components do not change. (When w = 0, the perihelion is located exactly on the ecliptic at the ascending node.)

We rotate the double-primed coordinate system around its X axis negatively by the inclination, getting

We next rotate the single-primed coordinate system around its Z axis negatively by the longitude of the ascending node, getting the position in heliocentric ecliptic coordinates:

The canonical velocity vector is

We rotate the coordinates for the sun-relative velocity the same way that we did for the heliocentric position. Again, note that since w=0, the first rotation doesn't change the vector components.

The state vector for this object at this point in its orbit is

Racist Celestial Mechanics
Jerry's Aryan Battle Page