newtonian gravity - 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.

Answer

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).