# Calculating a two-dimensional orbital path with infinite granularity (non-Euler integration)

For a game I am making, I am trying to calculate the position of an orbiting object around one or more bodies. I have successfully implemented this gravity simulation by calculating the force, then the acceleration, and then the velocity of the object, but the problem with this implementation is that the discretization of time (because the gravity is calculated during the update of each frame) causes instability of the orbit.

Therefore, I am looking for (a) function(s) of time that tells me the position of the orbiting object, therefore avoiding the granularity of (what I have come to understand as) the Euler integration. Any thoughts? The stuff I have read on here has been too complex to understand, so try to keep it at as low a level as you possibly can.

Edit: A parallel to what I am asking for is something similar to what is done in Kerbal Space Program, where the orbit is calculated ahead of time, except my game is in two dimensions.

• Great question. I wrote a seemingly identical program to yours and ran into the same problem. It's especially problematic when your planets come too close to each other. They shoot away from each other at ungodly velocities because of the discretized time imbalances that lead to asymmetric forces. This wouldn't happen in nature. There's got to be a simple way of countering this problem. I remember thinking about incorporating "conservation of energy" into my program because unimaginably high velocities violate conservation of energy. – Inquisitive Feb 9 '15 at 23:26
• Yeah, I sort of solved the asymmetry problem by keeping the orbiting object reasonably far away from the center of mass (by having the body being orbited around have a large enough diameter and letting the orbiting object collide with the body). It's not graceful, but it is close enough. – ben Feb 9 '15 at 23:38
• I put this question in the Game Development section of Stack Exchange: gamedev.stackexchange.com/q/93781/60448. I am not sure if it is a no-no to double post, but the person who responded left a pretty good (but unfortunate) answer. – ben Feb 9 '15 at 23:41
• My program was a 3-body 2D program. I tested it using the masses, distances, and actual initial conditions of the sun, earth, and moon. When I ran it, it worked! I was so pleased! Later on, I tested some disaster scenarios and, when my bodies got too close to each other, they shot off at warp speed in opposite directions. It was hilarious! – Inquisitive Feb 9 '15 at 23:46
• I take it the orbit has a non-zero eccentricity? If the eccentricity is $0$, though, you can just find the tangential velocity and use $\omega = \frac{v}{r}$ to find the change in the angle over a period of time $t$. – HDE 226868 Feb 10 '15 at 1:25

For two bodies this is relatively easy as the equations of motion describe a conic (an ellipse for a closed orbit, a hyperbola for an "open" orbit). You can use the vis viva equation to get the parameters of the orbit (semi major axis etc) from the given initial conditions, and the rest follows.

For an ellipse, you can also express the position as a function of time as

$$v_x = a \sin \omega t\\ v_y = b \cos \omega t$$

where your challenge is to find $a$, $b$ and $\omega$. For this you would like to know the velocity and distance at the furthest point, but you can in principle solve this for any place along the orbit.

For example, for a satellite in elliptical orbit above a massive planet (mass $M$), for given initial velocity $v$ in the radial direction, at a distance $r$ from the center of the planet, we can know the following (see http://en.wikipedia.org/wiki/Vis-viva_equation):

$$v^2 = GM\left(\frac{2}{r}-\frac{1}{a}\right)$$

In this expression, $a$ is the semimajor axis. We can calculate it from the given initial conditions,

$$\frac{1}{a} = \frac{2}{r}-\frac{v^2}{2GM}$$

Now we need $b$; again, following the wikipedia page quoted above, and for the conditions given (v perpendicular to radial vector at distance r, so we know angular momentum of the orbit $L=mvr$), we can write

$$b = vr\sqrt{\frac{a}{GM}}$$

With $a$ and $b$ calculated, we just need to find $\omega$ which follows again from the initial conditions, since we know

$$\omega = \frac{v}{r}$$

at the initial position. And with $a$, $b$ and $\omega$ known you can substitute into the above equations.

Incidentally it is more practical, when you have more than one gravitational object, to do the integration "properly". I recommend using a fourth-order Runge-Kutta integration, which is quite stable for this kind of problem. You can see http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf for some ideas on how to do this - it's a bit hard to follow. I actually wrote a piece of code that did this in Visual Basic a little while ago - it was actually taking into account both gravity and atmospheric drag. I reproduce it here without warranty...

Option Explicit

Const gm = 9.8 * 6300000# * 6300000# ' approximate value of GM on earth surface
' the # signifies a double precision "integer"
Const C_d = 0.5 * 0.03 * 0.03 * 0.2  ' drag coefficient for particular object
Const R_e = 6300000#                 ' approximate radius of earth
Const mass = 25                      ' mass of object
Const k1 = 12000#                    ' constants used to fit density of atmosphere
Const k2 = 22000#                    ' again, approximate!

Function rho(h)
' return density as a function of height
If h < 0 Then rho = 1.2 Else rho = 1.225 * Exp(-(h / k1 + (h / k2) ^ (3 / 2)))
End Function

Sub accel(x, y, vx, vy, ByRef ax, ByRef ay)
' compute acceleration of object at location x,y with velocity vx, vy
' return values in ax, ay
Dim r As Double, v2 As Double, v As Double, r3 As Double, h As Double
Dim density
r = Sqr(x * x + y * y) ' sqrt() in most languages... this is VBA
v2 = vx * vx + vy * vy
v = Sqr(v2)
r3 = 1# / r ^ 3
density = rho(r - R_e)
' don't need the second term if you don't care about drag!
ax = -gm * x * r3 - vx ^ 3 * C_d * density / (v * mass)
ay = -gm * y * r3 - vy ^ 3 * C_d * density / (v * mass)
End Sub

Function rk(ByRef x, ByRef y, ByRef vx, ByRef vy, dt)
' implement one Runge-Kutta fourth order stop
' this requires calculating the acceleration at verious locations
' for a single "step" in the algorithm
Dim ax As Double, ay As Double
Dim kx1, kx2, kx3, kx4, ky1, ky2, ky3, ky4, lx1, lx2, lx3, lx4, ly1, ly2, ly3, ly4, dt2
' get acceleration at initial point
accel x, y, vx, vy, ax, ay
' half time step:
dt2 = dt / 2
kx1 = dt2 * ax
ky1 = dt2 * ay
lx1 = dt2 * vx
ly1 = dt2 * vy
' get acceleration at new location
accel x + lx1, y + ly1, vx + kx1, vy + ky1, ax, ay
kx2 = dt2 * ax
ky2 = dt2 * ay
lx2 = dt2 * (vx + kx1)
ly2 = dt2 * (vy + ky1)
' get acceleration at half way point:
accel x + lx2, y + ly2, vx + kx2, vy + ky2, ax, ay
kx3 = dt * ax
ky3 = dt * ay
lx3 = dt * (vx + kx2)
ly3 = dt * (vy + ky2)
' compute acceleration for third combination of position and velocity:
accel x + lx3, y + ly3, vx + kx3, vy + ky3, ax, ay
kx4 = dt2 * ax
ky4 = dt2 * ay
lx4 = dt2 * (vx + kx3)
ly4 = dt2 * (vy + ky3)
' compute best approximation to new x, y, vx, vy
x = x + (lx1 + 2# * lx2 + lx3 + lx4) / 3#
y = y + (ly1 + 2# * ly2 + ly3 + ly4) / 3#
vx = vx + (kx1 + 2# * kx2 + kx3 + kx4) / 3#
vy = vy + (ky1 + 2# * ky2 + ky3 + ky4) / 3#

End Function


Calling the RK function multiple times you can now track the orbit, and you will find it is quite stable. You can easily add multiple massive objects (just update the acceleration formula) and it will still work (although orbits will become chaotic).

• Why not use a simple leap frog scheme instead of the RK4? It's symplectic and thus very stable and also very easy to implement. – Simeon Carstens Feb 10 '15 at 10:23
• @SimeonCarstens that would probably work but I happened to have this code lying around and it know it works well... – Floris Feb 10 '15 at 10:40