# How to programmatically calculate orbital elements using position/velocity vectors?

I would like to build some orbital mechanical software from scratch. I feel that this would be a great way to learn the steps required to calculate different Kepler orbital elements of an object, plot orbits, and predict where the object will be at some future time.

Specifically, I want to start with calculating the Keplarian elements. The inputs I would give the program would be the position and velocity vectors, along with a time. These input vectors will be relative to the center of the Earth, so I may also need to do a coordinate transfer if I want to use a specific location on the surface as a reference point.

I have seen the math for calculating Kepler orbital elements from this book, and I know lots of software has been developed throughout the years to calculate them, but I am having a hard time bridging the two. The math in the book is slightly confusing, and I think it would be easier for me to understand if I saw the steps "written out" in a programming language.

• I can use R or Python. Vectorized code from most languages I should be able to translate to one of these two. Thanks! @Chris I am going to update the question once more with a little more info on my problems with translating the math into code. – Stu Sep 11 '13 at 15:52
• Checkout orsa.sourceforge.net for their solutions/methods. Long discussion here physicsforums.com/showthread.php?t=232778. And for Python basic F=ma simulation: fiftyexamples.readthedocs.org/en/latest/gravity.html – user6972 Sep 12 '13 at 1:15
• FWIW, if your goal is computing future position, there is normally little reason to convert from position and velocity vector to Keplerian elements. Just compute and apply drag and the force of gravity at your current location and velocity and integrate forward in small steps. – RickNZ Oct 7 '18 at 23:36

Given Earth-centered, inertial (ECI) position and velocity vectors $\vec{r}$ and $\vec{v}$, you can directly solve for the classical orbital elements $(a,e,i,\Omega,\omega,\nu)$ as follows (algorithms first, followed by pseudocode at the bottom):

First, solve for the angular momentum $$\vec{h}=\vec{r} \times \vec{v}$$

then the node vector $$\hat{n}=\hat{K} \times \vec{h}$$ which will be used later.

The eccentricity vector is then $$\vec{e} = \frac{(v^2-\mu /r)\vec{r}-(\vec{r} \cdot \vec{v})\vec{v}}{\mu}$$

and $e=|\vec{e}|$.

Specific mechanical energy is $$E=\frac{v^2}{2}-\frac{\mu}{r}$$

If $e\neq 1$, then $$a = -\frac{\mu}{2E}$$ $$p=a(1-e^2)$$ Otherwise, $$p=\frac{h^2}{\mu}$$ $$a=\infty$$

Now, $$i=\cos^{-1}{\frac{h_K}{h}}$$ $$\Omega=\cos^{-1}{\frac{n_I}{n}}$$ $$\omega=\cos^{-1}{\frac{\vec{n}\cdot\vec{e}}{ne}}$$ $$\nu=\cos^{-1}{\frac{\vec{e}\cdot\vec{r}}{er}}$$

And you'll need to make the following checks: If $n_J<0$, then $\Omega=360^{\circ}-\Omega$,

If $e_K<0$, then $\omega=360^{\circ}-\omega$, and

If $\vec{r}\cdot\vec{v}<0$, then $\nu=360^{\circ}-\nu$.

Note that you'll run into problems (singularities) for certain cases: circular orbits ($e\approx 0$), and equatorial orbits ($i\approx 0$), particularly. In these cases you normally introduce a new, less troublesome variable, like mean longitude or true longitude of perigee.

h=cross(r,v)
nhat=cross([0 0 1],h)

evec = ((mag(v)^2-mu/mag(r))*r-dot(r,v)*v)/mu
e = mag(evec)

energy = mag(v)^2/2-mu/mag(r)

if abs(e-1.0)>eps
a = -mu/(2*energy)
p = a*(1-e^2)
else
p = mag(h)^2/mu
a = inf

i = acos(h(3)/mag(h))

Omega = acos(n(1)/mag(n))

if n(2)<0
Omega = 360-Omega

argp = acos(dot(n,evec)/(mag(n)*e))

if e(3)<0
argp = 360-argp

nu = acos(dot(evec,r)/(e*mag(r))

if dot(r,v)<0
nu = 360 - nu


Note: this follows from the method laid out in Fundamentals of Astrodynamics and Applications, by Vallado, 2007.

• Finally recreated in Python. Thanks for the help! It was pretty easy and now I understand the math much better. – Stu Mar 24 '14 at 14:27
• Is there any difference between nhat and n? I'm not sure what n is (assuming nhat is the node vector) but it's been used for calculating the argument of periapsis. I assume n is the same as nhat and n was used by mistake? – 9a3eedi Nov 24 '14 at 12:44
• For programming, I would prefer equivalent relations using the atan2 function, because that will perform quadrant resolutions automatically and protects you from needing to check the range before most of the uses of acos above. (Sometimes, numerical precision will cause the argument to an acos function to be very slightly autside the valid range of -1.00000 to +1.00000. When that happens, the program as shown will crash.) It is also possible for h to be zero, allowing a divide by zero error. – Matt Jessick Feb 4 '16 at 0:39
• What khat in the second line? I see in the sample code it's [0 0 1] is that the normal vector of the earth equatorial plane? – Ben Lu Sep 15 '16 at 7:15
• What is the node vector n^? What is K^? What is μ? What is r and how is it different from r⃗? What about v and v⃗? What is p? – Aaron Franke Apr 14 '19 at 7:29

The first key to figuring this out is getting your coordinate system correct. There are two commonly used coordinate systems for such things. They are the Earth Centered Earth Fixed (ECEF), and Earth Centered Interial (ECI) frames. At midnight, these two line up exactly, but they diverge in other times, based on the rotation of the Earth. ECEF works best for things on the Earth (If you are not moving, you should have 0 velocity. ECEF takes this into account, ECI velocity will have you move with the rotation of the Earth), ECI works best for things in Orbit (Orbiting objects don't care about the rotation of the Earth, at least, the physics doesn't care). Make sure the coordinate systems are correct!

Okay, so you have a position and velocity in ECI coordinates, what do you do? There is an excellent paper that describes the entire process, which I will copy the end formulas here. There is also a few good sources here, here, and here. I highly recommend reading them carefully. Uncertainty is much more difficult, so let's just assume you have a perfect knowledge of velocity and position. Specifically, the 6 classical Keplarian Elements are eccentricity (e), inclination (i), right ascension of the ascending node ($\Omega$), argument of perigee ($\omega$), semi-major axis (a), and time of perigee passage ($T_O$).

I should mention that I am primarily following the Laplace Method of Orbital Determination, there is a competing methodology known as the Gauss Method. But finally this came down to deciphering Matlab code.

Semi-Major Axis

$W_s = \frac{1}{2}*v^2s - \text{mus}./r;$

$a = -mus/2./W_s$; %semi-major axis

Eccentricity

L = [rs(2,:).*vs(3,:) - rs(3,:).*vs(2,:);...
rs(3,:).*vs(1,:) - rs(1,:).*vs(3,:);...
rs(1,:).*vs(2,:) - rs(2,:).*vs(1,:)]; %angular momentum


$p = \sum{L^2}./mus;$ %semi-latus rectum

$e = \sqrt{1 - p/a}; %eccentricity$

Inclination

$I = atan(\frac{\sqrt{L(1,:)^2 + L(2,:)^2}}{L(3,:)});$

Arguments of Pericenter

$\omega = atan2(\frac{(vs(1,:).*L(2,:) - vs(2,:).*L(1,:))./mus - rs(3,:)./r)./(e.*sin(I))}{((\sqrt{L2s}.*vs(3,:))./mus - (L(1,:).*rs(2,:) - L(2,:).*rs(1,:))./(\sqrt{L2s}.*r))./(e.*sin(I)))}$

Longitude of Ascending Node

$\Omega = atan2(-L(2,:),L(1,:));$

Time of Perigee Passing:

$T_0 = -(E - e.*sin(E))./\sqrt{mus.*a.^-3}$

• Laplace's method is one of initial orbit determination from angles measurements, and (I think) far beyond the scope of what OP is looking for. If you have ECI position and velocity, getting Keplerian elements is just a simple coordinate transformation. – user29 Sep 11 '13 at 15:45
• @Chris: I knew there had to be an easier way to do the transformation... Sigh. – PearsonArtPhoto Sep 11 '13 at 15:49
• What about the coordinate system in terms of handedness and orientation? From what I can tell, most guides are using Z-is-up. How would these formulas be implemented in a Y-is-up left-handed system, like Unity, or a Y-is-up right-handed system, like Godot? – Aaron Franke Apr 17 '19 at 8:04

OrbitalPy has a handy elements_from_state_vector function that does just that:

https://github.com/RazerM/orbital/blob/0.7.0/orbital/utilities.py#L252

You can check that the math is the same as in user29 answer.

• Hey, that's pretty cool! I'm looking forward to trying it. In the second example in the docs, titled "Create molniya orbit", can OrbitalPy implement precession? Is there a place to add a J2? (and I think Molniya should be capitalized - I think it qualifies as a proper name). – uhoh Apr 9 '16 at 0:47
• I'm not familiar with OrbitalPy myself, but from the limited analysis I did of the source code it appears it does pure two-body keplerian propagation, without any perturbation, so no ellipsoid corrections. – alexamici Apr 9 '16 at 12:15
• OK, thanks for the info - I'll give it a whirl. – uhoh Apr 10 '16 at 7:07
• Is this written with a "Z-is-up" coordinate system in mind? If I have a "Y-is-up" coordinate system, should I swap out h.z with h.y and [0, 0, 1] with [0, 1, 0]? – Aaron Franke Apr 15 '19 at 5:26
• "Argument of periapsis is the angle between eccentricity vector and its x component." Why? – Aaron Franke Apr 15 '19 at 5:46