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.

  • 2
    $\begingroup$ 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. $\endgroup$
    – Stu
    Commented Sep 11, 2013 at 15:52
  • $\begingroup$ 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 $\endgroup$
    – user6972
    Commented Sep 12, 2013 at 1:15
  • $\begingroup$ 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. $\endgroup$
    – RickNZ
    Commented Oct 7, 2018 at 23:36

3 Answers 3


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.

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)
   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.

  • 2
    $\begingroup$ 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? $\endgroup$
    – 9a3eedi
    Commented Nov 24, 2014 at 12:44
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Feb 4, 2016 at 0:39
  • 2
    $\begingroup$ 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? $\endgroup$
    – Ben Lu
    Commented Sep 15, 2016 at 7:15
  • 2
    $\begingroup$ 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? $\endgroup$ Commented Apr 14, 2019 at 7:29
  • 2
    $\begingroup$ What is h_K and n_I? $\endgroup$ Commented Apr 15, 2019 at 3:03

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


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$


$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}$

  • 3
    $\begingroup$ 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. $\endgroup$
    – user29
    Commented Sep 11, 2013 at 15:45
  • $\begingroup$ @Chris: I knew there had to be an easier way to do the transformation... Sigh. $\endgroup$
    – PearsonArtPhoto
    Commented Sep 11, 2013 at 15:49
  • $\begingroup$ 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? $\endgroup$ Commented Apr 17, 2019 at 8:04

OrbitalPy has a handy elements_from_state_vector function that does just that:


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

  • $\begingroup$ 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). $\endgroup$
    – uhoh
    Commented Apr 9, 2016 at 0:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – alexamici
    Commented Apr 9, 2016 at 12:15
  • $\begingroup$ 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]? $\endgroup$ Commented Apr 15, 2019 at 5:26
  • $\begingroup$ "Argument of periapsis is the angle between eccentricity vector and its x component." Why? $\endgroup$ Commented Apr 15, 2019 at 5:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.