State Vectors to Orbital Elements
Racist Celestial Mechanics, Lecture 2
My purpose with this page is to show how to determine the elements of an orbit when its heliocentric position and velocity are known for the same instant of time.
—Jerry Abbott
Revised: 9 August 2003
State vector is a combination of the components of the position and velocity vectors of an orbit, referred to the same coordinate system and to the same instant of time. I will show distances in both astronomical units and meters, since different equations require one or the other. I will express the speeds in meters per second wherever possible.
The default coordinate system will be classical heliocentric ecliptic coordinates, which has Earth's orbit as its XY plane, with its X axis positive in the direction of the Vernal Equinox and its Z axis positive in the direction of Earth's angular momentum vector.
My earlier paper on this subject can be found here.
We begin with the following state vector.
[ x , y , z , VX , VY , VZ ]
We first find the heliocentric distance and the sun-relative speed from
r = { x2 + y2 + z2 }1/2
V = { VX2 + VY2 + VZ2 }1/2
We appreciate MKS value of the gravitational parameter,
GM = 1.32712440018E+20 m3 sec-2
We solve for the semimajor axis of the orbit by putting MKS values for (r) and for (V) into the Vis Viva equation:
a = { 2/r - V2/GM }-1
|
Alert! Different systems of units should not mix in the same equation, just as different races of people should not mix in the same country. When using the MKS value of GM, the distances must be converted from astronomical units to meters. One AU equals 1.49597870691E+11 meters. You can re-convert the resulting distance back into AU once the calculation is finished. |
| If a>0 then... | The orbit is an ellipse. |
| If a<0 then... | The orbit is a hyperbola. |
Although I dislike the idea of representing any scalar distance by a negative number, and generally reverse the sign to make the semimajor axis of a hyperbola positive, in this lecture I will maintain it as a negative scalar. Otherwise, it would be necessary to partition too much of the following math into two separate cases.
The angular momentum per unit mass, h, is equal to the "cross product" of the heliocentric radius vector and the velocity vector at the apsidal endpoint: r x V. Angular momentum is a conserved quantity, and the components of h are constant for the entire transfer orbit.
hX = y VZ - z VY
hY = z VX - x VZ
hZ = x VY - y VX
As you can see, the cross product is a particular method for finding differences in ordered products of the components of two vectors. In fact, the method corresponds to finding the determinant of a three-dimensional matrix. I show the procedure, just above, in explicit algebra because I promised that I would not skip any steps or use any high-order notation.
h2 = hX2 + hY2 + hZ2
We obtain two intermediate quantities,
p = h2 / GM
q = x VX + y VY + z VZ
By the way, the quantity (q) is often called the "dot product" of the position vector and the velocity vector, written as r . V, because it is the sum of the products of the components (grouped by dimension) of the two vectors.
From here, we calculate the eccentricity of the orbit.
e = { 1 - p / a }1/2
If the orbit is an ellipse (i.e., if the semimajor axis found above is positive), there is a second method for finding the eccentricity:
eX = 1 - r / a
eY = q / ( a GM )1/2
e = { eX2 + eY2 }1/2
The values for the eccentricity, e, found by the two methods should agree with each other for elliptical orbits. For hyperbolic orbits (with negative values for the semimajor axis), this second method does not apply.
The inclination, i, is found from
i = ArcCos { hZ / h }
The longitude of the ascending node, W, of the transfer orbit can be calculated as follows:
| W = arctan2(hX , -hY) |
|
Definition of the function: arctan2 function arctan2( sin g , cos g )begin if cos g = 0.0 thenend The arctan function is the usual, principal value arctangent function, like the one on your calculator. The arctan2 function is the two-argument version of the arctangent, and it is superior because it corrects the output angle for the quadrant indicated by the input sine and cosine of the angle g (gamma). |
The next element we seek is the argument of perihelion, w. On the way, we will find the true anomaly, q.
| qX = h2 / (r GM ) - 1 | |
| qY = h q / (r GM ) | |
| q = arctan2( qY , qX ) | cos(w+q) = (x cos W + y sin W) / r |
| If i=0 or i=p radians, then... | sin(w+q) = (y cos W - x sin W) / r |
| If i is neither 0 nor p radians, then... | sin(w+q) = z / (r sin i) |
| w = arctan2[ sin(w+q) , cos(w+q) ] - q | |
Where q is the true anomaly of the object's current position in its orbit. The true anomaly is zero at perihelion, but it will be or p radians for the aphelion. Correct w to the interval [0, 2p), if necessary.
Elliptical Orbits.
If the orbit is an ellipse, we find the eccentric anomaly, u, from
u = arctan2( eY , eX )
And the mean anomaly, M, is found from Kepler's equation:
M = u - e sin u
The mean anomaly (along with a specified time, called an "epoch") is sometimes used instead of a time of perihelion passage. We've completed the determination of the orbit's elements as far as we can, unless we assign a specific time to the original state vector and then calculate the time of perihelion passage, T, as follows:
P = (365.256898326 days) a3/2
Enter (a) in astronomical units, and get the orbital period, P, in days. (Note that this symbol is a capital P, different than the smaller p used above as an intermediate quantity.)
m = 2p / P
The mean motion, m, will be in radians per day.
T = t - M / m
where (T) is the time of perihelion passage and (t) is the epoch of the mean anomaly (and also of the original state vector), both expressed in Julian Date or in some other specified running count of days.
Hyperbolic Orbits.
If the orbit is a hyperbola, we find the eccentric anomaly from
u' = ArcCosh { ( e + cos q ) / ( 1 + e cos q ) }
In case it is inconvenient for someone to calculate the ArcCosh (inverse hyperbolic cosine) function, the equivalent expression in natural logarithms is
Q = ( e + cos q ) / ( 1 + e cos q )
u' = ln { Q + ( Q2 - 1 )1/2 }
The sign of the hyperbolic eccentric anomaly is the same as the sign of the sine of the true anomaly:
If sin q < 0 then u = -u'
If sin q = 0 then u = 0
If sin q > 0 then u = u'
One does not correct the eccentric anomaly of a hyperbolic orbit to the interval [0,2p); if it comes out negative, it must be left that way. Furthermore, while it is desirable to keep the angles in radians throughout these calculations, it is necessary to do so when using Kepler's equation in either the elliptic or the hyperbolic forms.
The mean anomaly of is found from
M = e sinh u - u
The mean motion is
m = (86400 sec/day) { GM / a3 }1/2
One should use MKS values for (a) and GM to get the result in the customary radians per day. If the factor (86400 sec/day) had not been present, the result would have been in radians per second. Then the time of perihelion passage is found from
T = t - M / m
where (T) is the time of perihelion passage and (t) is the epoch of the mean anomaly (and also of the original state vector), both expressed in Julian Date or in some other specified running count of days.
Worked Example #1
We begin with the following state vector.
| x | +1.000212261 AU |
| y | -0.098871817 AU |
| z | +0.000000037 AU |
| VX | -17921.9 m/sec |
| VY | +27790.4 m/sec |
| VZ | +129.6 m/sec |
We find the heliocentric distance and the sun-relative velocity.
r = 1.005087162 AU = 1.503588993E+11 meters
V = 33068.4 m/sec
We find the semimajor axis with the Vis Viva equation.
a = +1.975599349E+11 meters = +1.320606597 AU
Since (a) is positive, the orbit is an ellipse. We next find the angular momentum vector and its magnitude:
hX = -1.917069146E+12 m2 / sec
hY = -1.939209853E+13 m2 / sec
hZ = +3.893184055E+15 m2 / sec
h = 3.893232823E+15 m2 / sec
We determine the intermediate quantities.
p = 1.142113114E+11 meters = 0.763455462 AU
q = -3.092695342E+15 m2 / sec
eX = +0.238920081
eY = -0.603992973
We determine the eccentricity by both methods.
| e = 0.649530843 | (by the first method) |
| e = 0.649530843 | (by the second method) |
If we didn't know which of these results was closer to the true value, we might consider their average to be the correct value:
e = 0.649530843
Notice that e is in the interval [0,1) and is therefore the orbit is an ellipse in agreement with the sign of the semimajor axis resulting from the Vis Viva equation. If it had been otherwise, something would be wrong with our calculations! Also, if the first and second methods had given very different results for the eccentricity, we would know that something was wrong somewhere.
We find the inclination to the ecliptic.
i = 0.005005277 radians
We find the longitude of the ascending node.
W = 6.184647216 radians
We find the true anomaly
qX = -0.240408702
qY = -0.603401999
q = 4.333243586 radians
We find the argument of the perihelion.
cos (w + q) = 1.000000000
sin (w + q) = 0.000007355
w = 1.949949076 radians
We find the eccentric anomaly.
u = 5.089068535 radians
We find the mean anomaly.
M = 5.693061509 radians
The period and mean motion of the orbit are
P = 1.517632591 years = 554.3175392 days
m = 0.011334993 radians/day
We find the time since the last perihelion passage.
dt = 502.255 days
Worked Example #2
We begin with the following state vector.
| x | +0.603293460 AU |
| y | -2.093152513 AU |
| z | -0.010132850 AU |
| VX | +17432.1 m/sec |
| VY | +69547.6 m/sec |
| VZ | +355.1 m/sec |
We find the heliocentric distance and the sun-relative velocity.
r = 2.178383143 AU = 3.258814797E+11 meters
V = 71699.9 m/sec
We find the semimajor axis with the Vis Viva equation.
a = -3.067509859E+10 meters = -0.205050369 AU
Since (a) is negative, the orbit is a hyperbola. We next find the angular momentum vector and its magnitude:
hX = -5.768951470E+12 m2 / sec
hY = -5.847277550E+13 m2 / sec
hZ = +1.173530313E+16 m2 / sec
h = 1.173545022E+16 m2 / sec
We determine the intermediate quantities.
p = 1.0377383748E+12 meters = 6.936852577 AU
q = -2.020478714E+16 m2 / sec
We determine the eccentricity
e = 5.901694093
Notice that e>1 and is therefore the orbit is a hyperbola in agreement with the sign of the semimajor axis resulting from the Vis Viva equation. If it had been otherwise, something would be wrong with our calculations!
We find the inclination to the ecliptic.
i = 0.005006788 radians
We find the longitude of the ascending node.
W = 6.184843098 radians
We find the true anomaly
qX = +2.184404268
qY = -5.482551519
q = 5.091539802 radians
We find the argument of the perihelion.
cos (w + q) = +0.369949675
sin (w + q) = -0.929051795
w = 6.282989337 radians
This value for the argument of the perihelion is almost equal to 2p radians, which if obtained would have been corrected to zero. In other words, this orbit has its perihelion and its ascending node at almost the same place.
We find the eccentric anomaly.
u' = 1.299193115
u = -1.299193115
The eccentric anomaly, u, is opposite in sign to (u') because sin q is negative.
We find the mean anomaly.
M = -8.714758278 radians
The mean motion of the orbit is
m = 0.185263818 radians/day
We find the since the perihelion passage.
dt = -47.040 days
The minus sign means that the current position on the orbit is prior to the perihelion by 47.040 days.