Converting Orbital Elements to Cartesian State Vectors

Question

Updated 2016-12-12

I'm writing an orbital mechanics program in Unity. I convert an orbiting object's position and velocity into orbital elements, then converting the orbital elements back into cartesian position vectors so that I can plot the entire orbit. I followed the equations on these two links:

Cartesian State Vectors to Keplerian Elements

Keplerian Elements to Cartesian State Vectors

The equations I am following to calculate state vectors assume the reference plane is X-Y. In Unity the reference plane is X-Z, with Y being the "up direction". So what I did was subtract 90 degrees from the inclination such that an orbit on the X-Z plane would be zero inclination, which ended up being the cause of my issues. So I deleted that modification, and now most of my orbits display correctly, for orbits that begin at the periapsis. But orbits beginning at the apoapsis are reversed. See the images below:

Here is a properly drawn orbit. I have a little spaceship orbiting, and it starts on the -X axis. I've drawn an incomplete orbit so that I can see the starting point of the orbit, which you can see begins on the -X axis.

Here I changed the starting state of the ship so that it starts at apoapsis. As you can see the starting point of the orbit is 180 degrees off, lying on the +X axis instead of the -X axis. In fact, the orbit is always drawn starting from the periapsis. What I'm trying to do is have the first point of the orbit be the ship's position at the point in time when I call the ConvertToCartesian function

Here's my code:

Vector3 ConvertToCartesian(float e, float a, float i, float O, float w, float t, float t0){
    float mu = G * M;
    float Mt; //Mean anomaly at time t
    float T = 2 * Mathf.PI * Mathf.Sqrt(a * a * a / mu); //Orbital period
    if (t == t0) {
        t = t0;
        Mt = 0; 
    } else {
        float deltaT = (T/60) * (t - t0); //divide time increments into 1/60th of the orbital period
        Mt = deltaT * Mathf.Sqrt (mu / Mathf.Pow (a, 3));
    }
    //Calculate eccentric anomaly using Newton's method
    float E = Mt;
    float F = E - e * Mathf.Sin (E) - Mt;
    int j = 0, maxIter = 30;
    float delta = 0.000001f;
    while (Mathf.Abs(F) > delta && j < maxIter) {
        E = E - F / (1 - e * Mathf.Cos (E));
        F = E - e * Mathf.Sin (E) - Mt;
        j++;
    }
    float nu = 2 * Mathf.Atan2 (Mathf.Sqrt (1 + e) * Mathf.Sin (E / 2), Mathf.Sqrt (1 - e) * Mathf.Cos (E / 2)); //True anomaly

    float rc = a * (1 - e * Mathf.Cos(E)); //Distance to central body

    Vector3 o = new Vector3(rc * Mathf.Cos(nu), rc * Mathf.Sin(nu), 0); 
                
    Vector3 odot = new Vector3 (Mathf.Sin (E), Mathf.Sqrt (1 - e * e) * Mathf.Cos (E), 0); //Velocity vector in the orbital frame
    odot = (Mathf.Sqrt (mu * a) / rc) * odot;

    float rx, ry, rz;
    rx = o.x; ry = o.y; rz = o.z;

    rx = ( o.x * (Mathf.Cos (w) * Mathf.Cos (O) - Mathf.Sin (w) * Mathf.Cos (i) * Mathf.Sin (O)) -
            o.y * (Mathf.Sin (w) * Mathf.Cos (O) + Mathf.Cos (w) * Mathf.Cos (i) * Mathf.Sin (O)));
    ry = (o.x * (Mathf.Cos (w) * Mathf.Sin (O) + Mathf.Sin (w) * Mathf.Cos (i) * Mathf.Cos (O)) +
        o.y * (Mathf.Cos (w) * Mathf.Cos (i) * Mathf.Cos (O) - Mathf.Sin (w) * Mathf.Sin (O)));
    rz = (o.x * (Mathf.Sin (w) * Mathf.Sin (i)) + o.y * (Mathf.Cos (w) * Mathf.Sin (i)));

    Vector3 r = new Vector3(rx, ry, rz); //Position vector

    return r / objectScale;

}


/* Convert a body's cartesian state vectors (position, r, and velocity, v) into Kepler elements
 */
void ConvertToKeplerElements(Vector3 r, Vector3 v){
    float mu = G * M;
    Vector3 h = Vector3.Cross (r, v); //Orbital momentum
    Vector3 n = Vector3.Cross (Vector3.forward, h); 
    eVector = ((v.magnitude * v.magnitude - mu / r.magnitude) * r - Vector3.Dot(r, v) * v) / mu; //Eccentricity
    float emag = eVector.magnitude;
    float E = v.magnitude * v.magnitude / 2 - mu / r.magnitude; //Specific mechanical energy
    if (eVector.magnitude != 1) {
        a = -mu / (2 * E); //Semi major axis
        float p = a * (1 - emag * emag);
    } else {
        //a = infinity
        a = 0;
        float p = h.magnitude * h.magnitude / mu;
    }

    inc = Mathf.Acos (h.z / h.magnitude); //Inclination

    //Longitude of ascending node (LAN), angle of body to node vector
    if (inc == 0 || inc == Mathf.PI) { 
        O = 0; //For an equatorial orbit (i = 0), LAN is undefined. By convention, set to 0. 
    } else {
        O = Mathf.Acos (n.x / n.magnitude); 
    }
    if (n.y < 0) {
        O = 2 * Mathf.PI - O;
    }

    //Argument of periapsis
    if (eVector.magnitude == 0) {
        w = 0; //For a circular orbit, by convention place w at ascending node (w = 0)
    } else {
        if (inc == 0 || inc == Mathf.PI) {
            w = Mathf.Atan2 (eVector.y, eVector.x);
        } else {
            w = Mathf.Acos (Vector3.Dot (n, eVector) / (n.magnitude * eVector.magnitude)); 
        }
    }
    if (eVector.z < 0) {
        w = 2 * Mathf.PI - w;
    }

    //True anomaly
    if (eVector.magnitude == 0) {
        if (inc == 0 || inc == Mathf.PI) {
            nu = Mathf.Acos (r.z/ r.magnitude); //For a circular non-inclined orbit, nu is angle of body from the x axis
            //(True Longitude)
        } else {
            //For a circular inclined orbit, nu is the angle between ascending node position vectors, (Argument of Latitude)
            nu = Mathf.Acos (Vector3.Dot (n, r) / (n.magnitude * r.magnitude));
        }
    } else {
        nu = Mathf.Acos (Vector3.Dot (eVector, r) / (eVector.magnitude * r.magnitude)); //True anomaly is angle between eccentricity and position vectors
        if (Vector3.Dot (r, v) < 0) {
            nu = 2 * Mathf.PI - nu;
        }
    }
}

To draw the orbit line I just call the function ConvertToCartesian in a loop with increasing values of t to get an array of position vectors.

Values: M = 3.3e16, GM ~= 2.2e6, r = (-50000, 0,0) For a circular orbit, v = (0,0, sqrt(GM/50000))

Follow up question: since the equations assume X-Y reference plane, how can I change it so that it conforms to the X-Z reference plane I'm using in my program?

Update

I updated the method ConvertToCartesian, such that it takes the True Anomaly (nu) calculated from the ConvertToKeplerElements function as an input. Then it calculates Eccentric Anomaly (E), then calculates Mean Anomaly (Mt) from that. When ConvertToCartesian is called in subsequent loops, Mt is updated with the relevant time step. Then I calculate E again from the updated Mt, then calculate nu from the updated E. This solves most cases. Problem is, sometimes nu is calculated 180 degrees off from where it should be, making the starting point of the displayed orbit 180 degrees off.

In my ConvertToKeplerElements function, according to the math, nu = 2pi - nu if r.v < 0. This happens, for example, when I set the initial r, v vectors such that the ship starts neither at periapsis or apoapsis but somewhere between. I know that this is correct, but when I calculate nu->E->nu, that 180 degree flip is lost in the calculation.

Here's the relevant parts of the updated method ConvertToCartesian.

        float mu = G * M;
    float E = Mathf.Acos((e + Mathf.Cos(nu)) / (1 + e * Mathf.Cos(nu)));
    float Mt = E - e * Mathf.Sin (E);

    //Mean anomaly at time t
    float T = 2 * Mathf.PI * Mathf.Sqrt(a * a * a / mu); //Orbital period
    float deltaT = (T/60) * (t - t0); //divide time increments into 1/60th of the orbital period
    Mt = Mt + deltaT * Mathf.Sqrt (mu / Mathf.Pow (a, 3));

    //Calculate eccentric anomaly using Newton's method
    float F = E - e * Mathf.Sin (E) - Mt;
    int j = 0, maxIter = 30;
    float delta = 0.000001f;
    while (Mathf.Abs(F) > delta && j < maxIter) {
        E = E - F / (1 - e * Mathf.Cos (E));
        F = E - e * Mathf.Sin (E) - Mt;
        j++;
    }
        
    //True anomaly
    nu = 2 * Mathf.Atan2 (Mathf.Sqrt (1 + e) * Mathf.Sin (E / 2), Mathf.Sqrt (1 - e) * Mathf.Cos (E / 2)); 

Here's the bad result I'm getting:

The true anomaly at the start of the orbit should be 1.5pi (~4.71), but when calculating it from E it is 0.5pi(~1.57)

Answer 1

You can find better instructions here:

Cartesian to Keplerian (PDF)

Keplerian to Cartesian (DOC)

I codded all steps given in these documents to check if I could find your error but it seems to work. I will give my code below (in Python) and perhaps if you compare it to yours you can find your error.

Running the code I got the follow output:

[ -9.31322575e-10  -1.09686156e-20   4.12053763e-05]
[  2.80478539e-13   9.09494702e-13  -8.42143002e-19]

Considering that $r$ is in meters and $\dot{r}$ is in meters per second these errors are quite acceptable.

OBS: If you want to transform $r$ and $\dot{r}$ to Earth Centered and Fixed there is an additional transformation described in the document!

from astropy.constants import G, M_earth, R_earth
from astropy import units as u
import numpy as np

mu = G.value*M_earth.value
Re = R_earth.value

#Test vectors
r_test = np.array([Re + 600.0*1000, 0, 50])
v_test = np.array([0, 6.5 * 1000, 0])
t = 0


def cart_2_kep(r_vec,v_vec):
    #1
    h_bar = np.cross(r_vec,v_vec)
    h = np.linalg.norm(h_bar)
    #2
    r = np.linalg.norm(r_vec)
    v = np.linalg.norm(v_vec)
    #3
    E = 0.5*(v**2) - mu/r
    #4
    a = -mu/(2*E)
    #5
    e = np.sqrt(1 - (h**2)/(a*mu))
    #6
    i = np.arccos(h_bar[2]/h)
    #7
    omega_LAN = np.arctan2(h_bar[0],-h_bar[1])
    #8
    #beware of division by zero here
    lat = np.arctan2(np.divide(r_vec[2],(np.sin(i))),\
    (r_vec[0]*np.cos(omega_LAN) + r_vec[1]*np.sin(omega_LAN)))
    #9
    p = a*(1-e**2)
    nu = np.arctan2(np.sqrt(p/mu) * np.dot(r_vec,v_vec), p-r)
    #10
    omega_AP = lat - nu
    #11
    EA = 2*np.arctan(np.sqrt((1-e)/(1+e)) * np.tan(nu/2))
    #12
    n = np.sqrt(mu/(a**3))
    T = t - (1/n)*(EA - e*np.sin(EA))

    return a,e,i,omega_AP,omega_LAN,T, EA

def kep_2_cart(a,e,i,omega_AP,omega_LAN,T, EA):

    #1
    n = np.sqrt(mu/(a**3))
    M = n*(t - T)
    #2
    MA = EA - e*np.sin(EA)
    #3
    #
    # ERROR WAS HERE
   #nu = 2*np.arctan(np.sqrt((1-e)/(1+e)) * np.tan(EA/2))
    nu = 2*np.arctan(np.sqrt((1+e)/(1-e)) * np.tan(EA/2))
    #4
    r = a*(1 - e*np.cos(EA))
    #5
    h = np.sqrt(mu*a * (1 - e**2))
    #6
    Om = omega_LAN
    w =  omega_AP

    X = r*(np.cos(Om)*np.cos(w+nu) - np.sin(Om)*np.sin(w+nu)*np.cos(i))
    Y = r*(np.sin(Om)*np.cos(w+nu) + np.cos(Om)*np.sin(w+nu)*np.cos(i))
    Z = r*(np.sin(i)*np.sin(w+nu))

    #7
    p = a*(1-e**2)

    V_X = (X*h*e/(r*p))*np.sin(nu) - (h/r)*(np.cos(Om)*np.sin(w+nu) + \
    np.sin(Om)*np.cos(w+nu)*np.cos(i))
    V_Y = (Y*h*e/(r*p))*np.sin(nu) - (h/r)*(np.sin(Om)*np.sin(w+nu) - \
    np.cos(Om)*np.cos(w+nu)*np.cos(i))
    V_Z = (Z*h*e/(r*p))*np.sin(nu) + (h/r)*(np.cos(w+nu)*np.sin(i))

    return [X,Y,Z],[V_X,V_Y,V_Z]

a,e,i,omega_AP,omega_LAN,T, EA = cart_2_kep(r_test,v_test)
r_test2, v_test2 = kep_2_cart(a,e,i,omega_AP,omega_LAN,T, EA)


print(r_test2 - r_test)
print(v_test2 - v_test)

Edit: I was able to plot these against JPL's Horizon data and with the small edit they are indeed accurate.

Edit2: There was an error in the implementation of #9 in cart_2_kep. It would never output all possible angles. The alternate formula using atan2 fixes this.

Edit3: Equation of V_z corrected, from a "-" to a "+"

Answer 2

I'm not sure - maybe I'm missing something, but your image seems to show left-handed cartesian coordinates. This could cause a complete orbit to look right (especially if there are no arrows to show the direction of motion), but an incomplete orbit would indeed look wrong because it would highlight the direction of motion (a plot of the velocity vector part of the state vector could also show this!)

Here is your image - if I used the right-hand rule, my thumb seems to point in your negative z direction.

This is what I have in Beldner, and I'll bet AutoCad and other 3D engineering programs look similar:

above: Screenshot from the 3D window of Blender.

Here are random images of the right hand Cartesian Coordinates. They look just like the Blender coordinates.

above: From Wikipedia.

above: From Wikipedia.

And these shows the difference between right-handed and left-handed:

above: Cartesian coordinate handedness from here

above: Cartesian Coordinate handedness from here

Answer 3

In both your original code and your update, your root approximation code is not taking advantage of the fact that E is between 0 and 360 degrees (0 and 2pi).

Changing the calculation of E to this may suffice if this the remaining problem:

  E = (E - F / (1 - e * Mathf.Cos(E))) % (Mathf.PI * 2);

Answer 4

First of all kudos for showing your work and the results illustrating your question so clearly!

Do the two suggestions at the bottom shed any light on it? Lighty edited from Wikipedia the Keplerian orbital elements are:

---

The main two elements that define the shape and size of the ellipse:

1. Eccentricity (e)—shape of the ellipse...
2. Semimajor axis (a)—the sum of the periapsis and apoapsis distances divided by two.

Two elements define the orientation of the orbital plane in which the ellipse is embedded:

3. Inclination (i)—vertical tilt of the ellipse with respect to the reference plane
4. Longitude of the ascending node (ω or Ω)—horizontally orients the ascending node with respect to the reference frame's vernal point

And finally:

5. Argument of periapsis (ω) the orientation of the ellipse in the orbital plane, measured from the ascending node to periapsis

6. True anomaly (ν, θ, or f) at epoch ($\mathbf{M_0}$ defines the position at a specific time (also known as the "epoch").

---

check your units:

float deltaT = (T/60) * (t - t0); //divide time increments into 1/60th of the orbital period
Mt = deltaT * Mathf.Sqrt (mu / Mathf.Pow (a, 3));

If 60 is unitelss, then deltaT has units of time squared. That gives Mt units of time. Then that means:

float E = Mt;
float F = E - e * Mathf.Sin (E) - Mt;

gives E units of time, and things that go inside trig functions like Sin() should be unitless.

note; the use of 60 to convert seconds to minutes for example, which would be called seconds per minute is unitless.

try different true anomalies ( t0's )

The true anomaly (#6) is called t0 in your program - have you tried just passing different values? If you pass different values between 0 and the period T, you should get the orbit to start at different places. $t0 = 0.5T$ should start as apoapsis.