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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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