Coding Challenge: The Three Body Problem

The three body problem is to solve the differential equations:

x(double dot) - 2*y(dot)*w^2 = (1-mu)/r1^3*(x1-x) + mu/r2^3*(x2-x)
y(double dot) +2*x(dot)*w - y*w^2 = (1-mu)/r1^3*(y1-y) + mu/r2^3*(y2-y)
z(double dot) = (1-mu)/r1^3*(z1-z) + mu/r2^2*(z2-z)

where "double dot" and "dot" is the second or first time derivative
r1 = ((x+mu)^2 + y^2 + z^2)^(1/2)
r2 = ((x-1+mu)^2 + y^2 + z^2)^(1/2)
mu = is the mass of larger body/Total mass (Earth-Moon system is roughly 0.02, but must be <=0.5)
w = angular frequency (but always 1)

x1-x = (current position in x) + mu
x2 - x = (current position in x) - 1 + mu
you may replace x with y or z

These equations give a solution to the three body problem in a rotating coordinate system. With the x-axis always pointing towards the next largest body and the origin being the largest body. For example, in the Earth-Moon system, Earth is the origin and as the Moon rotates around the Earth the x-axis will always point towards the moon; y-axis is in the orbital plain but perpendicular to x-axis and z-axis is perpendicular to the orbital plain.

Because it solution is in a rotating frame, to make sense of it you'll need to rotate it back to the Earth-fixed frame. To do this use a 3d rotation matrix on the solution you get from the diff eqs. with the angle being the time.

The units are such that the distance between the two largest bodies is always equal to one, which is why the angular frequency is always 1.

The solution gives the orbit of a small object, such a satellite. However it is assumed that body one is massive enough that its orbit is not perturbed massively by the other bodies, the second object is in a circular orbit, and the third body is much much smaller than the other 2. Again, the Earth, Moon, and satellite is an example of this.

***important*** solving these diff eqs requires an very accurate approach so you'll need to use an rk45 algorithm or something of the like with low error tolerance.

So finally the challenge is to animate the orbit of the Moon and the satellite in the three body system.
The Moon's orbit can be animated by a circle with the angle being time and radius being 1 since its a circular orbit.

This gives a very detailed explanation and also goes into zero velocity curves which would also be another fun challenge. Page 10 gives the equations I gave here.
[http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf](url)

here is an example from my code of what the results should look like with initial conditions of r0 = [1.13 0 0] v0 = [0 0 0] mu = 0.02. Note this represents a Lagrange point


another: r0 = [-0.04 0 0] v0 = [0 -9.74 0] mu = 0.02. Note this is roughly escape velocity for the system.