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)
