6
$\begingroup$

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: enter image description here

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)

$\endgroup$
  • 2
    $\begingroup$ Can you include some more information? At least one state vector $(\mathbf{r}, \ \dot{\mathbf{r}})$ and the $GM$ you are using and a screen shot of the final "wrong" orbit that results. Maybe I can try to reproduce the problem, and maybe someone else will recognize it. You should post some code - with a statement like "I'm sure the problem lies in that I'm not using the math correctly; code-wise..." showing your code would be an obvious step if your asking for help fixing it. $\endgroup$ – uhoh Dec 8 '16 at 13:50
  • 2
    $\begingroup$ If you can provide example values, I can place them into the orbital software I use therefore we can find out which one of your functions are producing the different results. $\endgroup$ – VolkanOzcan Dec 8 '16 at 14:10
  • 1
    $\begingroup$ I've found a better way to frame my question after trying to understand what's going on...should I make a new thread or edit my existing question? $\endgroup$ – echl Dec 10 '16 at 3:43
  • $\begingroup$ @echl However, if it is distinct, so that it requires a different answer, then you can definitely post a new separately, mention and link this question and explain why it is new and different and needs a different kind of answers. $\endgroup$ – uhoh Dec 10 '16 at 10:23
  • $\begingroup$ @uhoh Okay, thank you. I will remember to do that for future questions. Sorry, I'm new to posting here. $\endgroup$ – echl Dec 10 '16 at 19:00
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Very helpful answer (and yay for python)! I wanted to give it a test-drive, but my OpenOffice couldn't seem to render the equations in the Kep to Cart doc. However uploading to google docs once allows re-downloading in .odt and .pdf formats and the equations seem to be just fine there. $\endgroup$ – uhoh Dec 9 '16 at 0:20
  • $\begingroup$ Care to try another question? $\endgroup$ – uhoh Dec 9 '16 at 0:43
  • $\begingroup$ How did you calculate EA in kep2cart? I don't see it anywhere there. $\endgroup$ – echl Dec 12 '16 at 23:21
  • 1
    $\begingroup$ EA is calculated in the 11th step of cart_2_kep, then passed as an argument to kep_2_cart $\endgroup$ – user17622 Dec 14 '16 at 12:31
  • $\begingroup$ I was not able to corroborate results here with output from JPL Horizons e.g. for Earth's orbital parameters vs. vectors. It would be really great if you could check your results against that and add an example. $\endgroup$ – imallett Oct 27 '17 at 21:44
3
$\begingroup$

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);
$\endgroup$
1
$\begingroup$

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.

enter image description here

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

enter image description here

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.

enter image description here

above: From Wikipedia.

enter image description here

above: From Wikipedia.

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

enter image description here

above: Cartesian coordinate handedness from here

enter image description here

above: Cartesian Coordinate handedness from here

$\endgroup$
  • $\begingroup$ Ah yes, you're right about that...but somehow my orbits still do go in the correct direction. Right now it just seems like the orbit is always drawn starting from the periapsis of orbit; I'm not sure if that has to do with the left-handedness of Unity's coordinate system, or the way that I am calculating the state vectors. Could the problem lie in my code where I calculate mean anomaly Mt and eccentric anomaly E? For t=0, Mt=0 -> E = 0 Therefore rc = a * (1 - e * Mathf.Cos(E)) = minimum value Therefore at time t0, the position vector is always the periapsis. $\endgroup$ – echl Dec 10 '16 at 19:23
  • $\begingroup$ This isn't an answer, it's just a very long observation that Unity does indeed use a left-handed coordinate system. $\endgroup$ – Aaron Franke Apr 17 at 8:24
  • $\begingroup$ @AaronFranke you might notice that there isn't a question either. The title is "Converting Orbital Elements to Cartesian State Vectors" unpunctuated, and the only question mark anywhere in the OP"s post is a "follow-up question." In this case the OP was having some trouble in several ways, and four answers have been posted to help them address the technical problems they were having. It's a little unconventional, and most of the time people with programming problems never get a helpful answer, or they never reply to comments and just disappear. This one was much more productive than average. $\endgroup$ – uhoh Apr 17 at 14:56
1
$\begingroup$

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:

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

And finally:

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

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

$\endgroup$
  • $\begingroup$ I can change the starting point of the orbit arbitrarily, but the problem I am having is that it needs to match the position of the ship at t=t0, which it doesn't always. I edited my post to show some modifications I made. $\endgroup$ – echl Dec 12 '16 at 22:58
  • 1
    $\begingroup$ @echl OK I think I see what you mean now. I'm working on it (your code translated to python (I am allergic to curly braces)) and I'll get back to you soon. To double check, I think you want a function that takes a physical starting point and calculates the correct $t_0$ so that when everything starts at time $t=0$ the orbit begins from that starting point. And maybe there might or might not be some issues with the nu = 2pi - nu flipping. $\endgroup$ – uhoh Dec 13 '16 at 11:43
  • 1
    $\begingroup$ That's exactly correct! I've possibly narrowed down to where I calculate E from nu (to recap, the calculation "pipeline" is: r,v->nu->E->Mt->Mt advanced by (t - t0)->E->nu->r,v). I tried two equations from [Wikipedia, in the section "From the true anomaly". Calculating E from CosE =... works for some orbits, and calculating E from TanE = ... works for some other orbits. $\endgroup$ – echl Dec 13 '16 at 23:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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