# I am having trouble calculating my classic orbital elements and am at a loss on where to look

I am using Unity engine and C#, and am trying to work out a way to plot future 2D orbital trajectories. I currently have a two-body setup. The planet has a mass of 1000000, the satellite a mass of 1. In the image above you can see a pink line indicating the orbital path (PAST, NOT FUTURE) created by a trail renderer, which indicates the math is wrong since the eccentricity is not anywhere between 0 and 1.

I have also noticed that the AP and PE are incorrect since the PE is showing up as -AP. I know my math is wrong somewhere but have no idea where.

If someone much smarter than I is able to figure it out/fix this, it would be great if they could also point me in the direction of how I would go about using this to plot a future trajectory in real time.

void JustKeplerThings(Vector3 curPos, Vector3 shipCurVel, GameObject orbiting) {
//Gravitational Constant
float G = 6.67408f * Mathf.Pow(10, -11);

//Mass of ship
float m1 = thisShip.mass;

//Mass of the larger body
float m2 = -orbiting.GetComponent<PointEffector2D>().forceMagnitude;
float M = m2;

//Standard Gravitational Parameter
float µ = G * (m1 + m2);
//float µ = G * M;

//Relative Position Vector
Vector3 r = FindRelativePosition(transform, orbitee.transform.position);

//Relative Velocity Vector
Vector3 v = shipCurVel - new Vector3(orbitee.GetComponent<Rigidbody2D>().velocity.x, orbitee.GetComponent<Rigidbody2D>().velocity.y, 0);

//Specific Angular Momentum
Vector3 h = Vector3.Cross(r, v);

//Eccentricity Vector
Vector3 evec = (Vector3.Cross(v, h) / µ) - (r / Vector3.Magnitude(r));

//Eccentricity
float e = Vector3.Magnitude(evec);

//Mechanical Energy
float E = (Mathf.Pow(Vector3.Magnitude(v), 2) / 2) - (µ / Vector3.Magnitude(r));

//Semi-Major Axis
float a = µ / (2 * E);
if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
a = Mathf.Infinity;

//Semi-Latus Rectum Check
float p = a * (1 - Mathf.Pow(e, 2));
if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
p = Mathf.Pow(Vector3.Magnitude(h), 2) / µ;

//Inclination (2D)
float i = 0;

//Longitude of Ascending Node (2D)
float Ω = 0;

//Argument of Periapsis
float ω = Mathf.Atan2(evec.y, evec.x);
//ω = (2 * Mathf.PI) - ω;

//True Anomaly
float t = Mathf.Acos((Vector3.Dot(evec, r)) / (e * Vector3.Magnitude(r)));
if (Vector3.Dot(r, v) < 0)
t = (2 * Mathf.PI) - t;

//Apoapsis
float ap = a * (1 + e);

//Periapsis
float pe = a * (1 - e);

//Orbital Period
float T = (2 * Mathf.PI) * Mathf.Sqrt(Mathf.Pow(a, 3) / µ);
}


-

public static Vector3 FindRelativePosition(Transform origin, Vector3 position) {
Vector3 distance = position - origin.position;
Vector3 relativePosition = Vector3.zero;
relativePosition.x = Vector3.Dot(distance, origin.right.normalized);
relativePosition.y = Vector3.Dot(distance, origin.up.normalized);
relativePosition.z = Vector3.Dot(distance, origin.forward.normalized);
return relativePosition;
}

• I don't know the language you're using so I am only guessing. You have specific angular momentum as r X v? Maybe it should be (2 * r X v)/(period of orbit). – HopDavid Jun 11 '16 at 13:47
• I had to move things around in order to even calculate things that were dependent on eachother, but this does not work as T (orbital period) is returned as NaN. – Wafer Jun 11 '16 at 16:53
• Please could you state the language you have used. – Puffin Jun 11 '16 at 17:03
• Yea sorry, I edited my question but forgot to add it to the comment as well, C# in Unity. – Wafer Jun 11 '16 at 17:04
• First you need to conquer the NaNs-- they poison any following computation. After each calculation, print the result to the console with a line like Debug.Log( "x = "+x ); -- attack the first obviously wrong result you get, over and over, in turn, until everything looks right. There are more sophisticated methods of debugging but this is the one I still use most. – Russell Borogove Jun 11 '16 at 21:09

The computation you're doing is for a 1 million kg body with a satellite 150 meters from its center moving past it at 83 meters a second. 1 million kg is almost nothing in gravitational terms; Earth masses 6e24 kg. So the satellite is moving faster than escape velocity. This means orbital energy is by convention positive, so semi-major axis "appears" negative, and there's no solution possible for orbital period, because it ain't coming back.

Unity's physics simulation is using a force (maybe expressed in newtons if you assume Unity distance units are meters -- but maybe not, because it's Unity), not a mass; you're treating it as a mass, which doesn't make sense.

Since Unity caps 2D movement speeds in its physics system, and doesn't document the basis distance at which PointEffector2D force is equal to the value specified, it'll be problematic to get your computations to line up with the output of Unity's simulation, so I recommend writing your own gravitational simulation -- it's fairly straightforward. You'll have to be conscious of the scale of your simulation, i.e. decide how many meters are represented by one Unity distance unit and stick to it.

Regarding plotting the future path of a body, you have a couple of options.

One, once you've got your own gravity sim working, is to simply remember the last N positions of the sim, and put the orbiting body at the location of the N-step-old results, and draw future-path lines up to the latest position. You'll have to run the simulation N times in advance before displaying the body, of course.

Option two, if you only have one or two massive bodies in play, is to tackle the equations here or here to translate between your orbital state vectors and future-time coordinates. With more than two bodies, you can't solve it this way. (You can have one massive body and any number of infinitesimal ones that don't affect each other gravitationally, of course.)

• So I should keep my equation calculating a as -a = µ / (2 * E) but fix the scaling on my values? I'm going to need to write my own gravity then, as the mass of the planet is currently the pointeffector2D's max. – Wafer Jun 12 '16 at 2:40
• Unity caps 2d movement speeds in its physics system, and doesn't document the basis distance at which PointEffector2D force is equal to the value specified. You might be able to figure it out by trial and error, but I'd write my own gravity if I were you. You'll have to be conscious of the scale of your simulation, i.e. decide how many meters are represented by one Unity distance unit and stick to it. – Russell Borogove Jun 12 '16 at 2:47
• I've went ahead and tried to simulate my own gravity, I used the ISS orbiting earth and scaled things down by 100. But for some reason I haven't worked out yet, it's effected by frame-rate and time-scale. I also needed to make massive changes to the Gravitational Constant manually changing it until the numbers worked. Here is what I came out with Earth Radius: 63.71, Earth Mass: 5.972E+22, Ship Mass: 4194.55, Ship Velocity: 76.5, Ship Altitude: 4, Gravitational Constant: 1.59508e-21. – Wafer Jun 15 '16 at 20:57
• Are you using Time.deltaTime to get the elapsed time per frame? (Or the fixed-update equivalent?) You need to take elapsed time into account, and decide how much real time you want it to take for your simulated ISS to make an orbit to establish a time scale. Your v is in 100m/s units but your linear dimensions are in 100km units. – Russell Borogove Jun 16 '16 at 3:04