Calculating 2D Orbital Paths in Newtonian Gravity Simulation

Question

I am currently working on a project (for fun) for simulating Newtonian gravity and want to be able to visualize the future path of an object orbiting around a single attractor. At any given simulation step I know

1. The object's current position relative to the center of the attractor as a 2D Vector
2. The object's current velocity as a 2D Vector
3. The mass of the attractor
4. The mass of the object

Using these parameters I would like to be able to predict up to N steps ahead in the simulation so that I can draw a line between each point along the orbit.

The Position and Velocity in my simulation are measured in m and m/s while the mass is unitless. The Gravitational Constants is set to 1 and the masses for each object is 1000, for the attractor, and 10 for the object.

The simulation itself produces the result that I want, now I just want to be able to predict the future path.

During my own research, I have found that I need to use the Keplerian Elements in some way. The various examples, questions etc... I have found on stackexchange and elsewhere do not provide sufficient explanation for me to be able to work it into my simulation or they are specific for 3-dimensional geometry and provide incorrect outputs when I attempt to calculate them (NaN and Infinities).

For example, in a few places I have read that the semi-major axis and the total orbital energy should remain constant and from these I can get the eccentricity and various other properties, but when calculating them at any given step using the above-mentioned properties I have access to, they vary wildly or produce numbers that are so ridiculously high or low that they are essentially useless. I resolved the last problem and have a semi-major axis value that makes sense in the context of my simulation, but again it varies from step to step.

As an example:

Initial conditions

Attractor

• Mass = 1000
• Position = 0,0

Object

• Mass = 10
• Position = 5, 0
• Velocity = 0, 5

This produces an elliptical orbit. The semi-major axis should remain constant as no other forces are acting on the object, however, the value of the semi-major axis varies from 2-4 meters.

Answer 1

I have figured it out and to any future visitors to this question I am (hopefully) going to make your life much, much easier.

If you are like me, fairly lay in terms of astrodynamics then the sheer number of equations, their very technical explanations, and usage is a little daunting. With a little bit of effort, I now know how to iterate equations to be able to calculate future orbital positions.

What we will need are the velocity and positional coordinates of the object you want to predict as well as the mass of the orbital body and the mass of the object itself.

These need to be in m and m/s for no other reason than for this implementation they are assumed to be m and m/s but the mass can be unitless (at first).

Terms:

Attractor - The nearest massive body that is attracting the Object.

Object - The orbital body that is being attracted to the Attractor.

Step 1: Calculate the True Mass of your objects

True mass is defined as the KG mass of your Attractor and Object. To find this value you need to convert your unitless mass, using the gravitational constant, into kilograms.

If you don't know how to do this, first calculate the ratio of unitless mass to the gravitational constant in your simulation and use the value of that ratio to find true mass using the real gravitational constant.

(I cannot post images or else I would include the equations)

After this point when I refer to the gravitational constant, I will be using our real universe's value.

Step 2: Calculate μ

μ is a ratio of masses and the gravitational constant. If the mass of the attractor is large enough μ = GM, where G is the gravitational constant and M is the mass of the attractor. For my own implementation I found that μ can be calculated as:

μ = G ( M * M * M) / ((M+m)*(M+m))

Where:

M is the mass of the Attractor

m is the mass of the Object

G is the gravitational constant.

It is very very important that this value is correct. If what you are rendering looks weird, it is probably because this value was not calculated correctly and it would be the first thing I check. In my initial attempts, this was off by an order of magnitude, it should have been around 100 but it was outputting as 1000 because the mass unit conversion value I calculated was 10x larger than it should have been.

Step 3: Calculate the Orbital Elements

First up is the Semi-Major Axis. This is the radius from the furthest point in the orbit to the center (Aprox). It is represented with the symbol a and calculated as:

a = -(μ * len(pos)) / (len(pos) * len(vel) * len(vel) - (2 * μ))

Man I wish I could post images so these equations looked better, sorry!

Next up is the period of the orbit. This is used later so that we can step through the orbit incrementally. It is represented by the symbol T and its formula is:

T = 2 * π * sqrt((a * a * a) / μ)

Next, we need angular momentum. Now, I know that the post says 2D but in my own code, I am using Vector3's to store my information in case I decide to make it 3D later. This is important because normally to get the angular momentum you would take the cross product of the velocity and position. Since we are in 2D there is no true cross product for vectors. Instead, you can calculate the momentum as:

p = pos.x * vel.y - pos.y * vel.x

Now onto orbital energy, represented by the symbol ε. We use this to find the eccentricity of the orbit.

ε = (len(vel) * len(vel)) / 2 - (μ / len(pos))

After that we can get eccentricity, represented with the symbol e:

e = sqrt(1 + ((2 * ε * p * p) / (μ * μ)))

The second to last thing we need to calculate is the eccentricity vector, represented with ev:

ev = ((((len(vel) * len(vel)) / μ) - (1 / len(pos))) * pos) - ((dot(pos, vel) / μ) * vel)

I know, there are a lot of brackets but I want to make it clear what goes where. Blame StackOverflow for not letting me post equation images.

The final value is the argument of periapsis. Because we are in 2D we need to use this formula to calculate it:

w = atan2(ev.y, ev.x)

we will use this to orient our orbital path correctly.

Step 4: Keplers Equation

Now is where we get to predict the orbital positions using all the information we calculated above.

First, we need n. n is the radians / second that the Object is traveling around the Attractor. It is fairly straightforward to calculate, it is 2π/T where T is the period. This value is used when we calculate the mean anomaly in the predictive step.

Now we get into the loop. You can set the number of steps you want to integrate over to whatever you want, the higher the number the closer each point along the orbit will be. In my case, I have a step of 100, so in my for loop declaration, I will loop 100 times before exiting. In addition, you will want to define a max tries value because we need to do an integration (don't worry, not as scary as it sounds) to be able to get the correct Eccentric Anomaly.

Now that we are in the loop, we need to calculate the current mean anomaly:

ma = n * (T * (i / steps))

Where:

T is the period of the orbit in seconds

i is the current iterative step

n is speed from earlier (radians/second)

Now we do the integration to solve Kepler's Equation. This can be done using Newton's method. Basically, we take Kepler's Equation and its derivative and iterate over them until the delta value between the two of them is arbitrarily low.

This is why we need max tries, it is possible that this while loop iterates forever and never converges on a value in which case we need to bail out so we don't get stuck.

//Eccentric Anomaly
var ea = ma;
while(tries < maxTries)
{
    tries++;
    //Keplers equation
    var dx = ea - e * sin(ea) - ma;
    //Derivative of Keplers equation
    var dy = 1 - e * cos(ea);
    var dt = dx / dy;
    ea -= dt;
    if(abs(dt) < 10^-8) break;
}

See, the integration wasn't so hard. This loop is an implementation of Newton's method and will oscillate toward to correct value. It should be noted that his method is quick and dirty and produces okay results. But this is the only method I understand so this is why I used it.

Step 5: Calculating Orbital Position

Next up is where we calculate the orbital position. This is simple and I stole it directly from this post.

p = a * (cos(ea) - e)

q = a * sin(ea) * sqrt(1 - (e * e))

You can plot these two values as X and Y and you will get an orbit that looks correct in some cases, but we need to rotate it using the argument of periapsis we calculated earlier. Again, stolen from this post.

x = cos(w) * p - sin(w) * q

y = sin(w) * p - cos(w) * q

And that's it. You're done. You have the X and Y coordinate at the current step and you can plot it on the screen however you like. As the for loop continues to iterate, you will get the full orbital path from start to finish.

I hope this essay helps someone out there!

All the equations I used can be found on Wikipedia just by going to this page and navigating to the corresponding element.