# Calculating orbital elements from Cartesian Vectors

I'm having serious troubles with getting my calculations to work. I'm calculating the orbital elements from a position vector, and a velocity vector. (I also know the distance and mass of the bodies around)

The equation for angular momentum is as follows: h=r*v

The equation for eccentricity is: e=((v^2−μ/|r|)r-(r⋅|v|)v/µ

I plug in the position and velocity to get h, and I plug them into the other equation to get the eccentricity. These SHOULD be constants, and not change unless another force is put on them... yet they are not constant. In fact they change as the orbit goes around. In the diagram I have just left a trail to show the orbital path. I should be able to calculate the required information... however they rely on the constants, which are not at all constant. I reference the rigidbody's position and velocity in them, and I've entered the equations in correctly. My code for applying force between planets is as follows:

void CalculateForce(Transform body) {
float m1 = rbody.mass;
float m2 = body.GetComponent<Rigidbody> ().mass;
float distance = Vector3.Distance (transform.position, body.position);
float f = (m1 * m2) / Mathf.Pow (distance, 2);
Vector3 direction = rbody.position - body.GetComponent<Rigidbody> ().position;
}


The calculations are coming back incorrect. The eccentricity varies between 0.9 and 1.1. Anything over a 1 is not actually in orbit, and anything less is in orbit. However this remains in orbit as it travels around.

Can anybody help? (I've already posted to the space stack exchange and had no luck)

Update:

Code used to calculate eccentricity (1001f is the combined weight of both objects):

    Vector3 pos = GetComponent <Rigidbody2D> ().position;
Vector3 vel = GetComponent <Rigidbody2D> ().velocity;
Vector3 eccentricity = (Mathf.Pow (vel.magnitude, 2.0f) - (1001f / pos.magnitude)) * pos;
eccentricity -= Vector3.Dot (pos, vel) * vel;
eccentricity /= 1001f;

• Sounds like inaccuracies are piling up due to the physics integration. In the past, when I've needed to simulate orbits accurately, I've driven them straight off of Kepler's equations (ie. each frame I plug in the current time parameter, and solve for the position of the body at that moment). This is a little tricky since it's non-linear (calculating the anomalies is iterative), but ensures the object stays on the correct orbit forever. Would something like this be an option in your case? – DMGregory Jul 9 '15 at 23:48
• It could well be the cause. I've started a fresh project, and it's looking to give me a constant now. I'll look into putting the Kepler's equations in, if it doesn't work. – Jacob Morris Jul 9 '15 at 23:58
• @JacobM You might try changing float f = (m1 * m2) / Mathf.Pow (distance, 2); to float f = (m1 * m2) / Mathf.Pow (distance, 2.0f); Sometimes error can come from doing operations between integers and floats.(had the hardest time trying to code block in comment.. guess its not possible) – Stephan Jul 10 '15 at 15:41
• I've tried changing the 2 to 2.0f at all occurrences, but it appears to have no difference. I've updated my post with the code used to calculate eccentricity, if it helps. – Jacob Morris Jul 10 '15 at 16:47
• @Stephan Mathf.Pow() only takes floats, so a 2 is implicitly cast to a float. – DMGregory Jul 10 '15 at 17:32