2. THE THREE-BODY GRAVITY SIMULATION

Hi all!


This second programs is about a famous problem that Newton proposed because he couldn't find an analytical solution, and in fact there isn't any! The only way to calculate this movement is through computer simulations using numerical methods. Actually if you download the program you will see that the planets move in a crazy way! I's very difficult to find an stable solution as you have to fix it yourself by changing some parameters. If you find any interesting solution don't doubt to post a comment.

Picture:

Red and blue planets are orbiting around the green one and all are moving upwards.

Data sheet:

Involved physics:

- Newtonian gravity

Needed libraries:
-Visual

-Math


Explaining the code:

As always we have to import the libraries needed in this program. Also set up the scene dimensions, auto-scaling and title:

from visual import *
import math
scene.width = 700
scene.height = 600
scene.title = 'BOUNCING BALL'
scene.autoscale= False
scene.fullscreen = False


After that, lets us to define some parameters that can be changed for different results. They are the mass of the 3 planets/objects, their radius according to their mass (mass is proportional to volume), the gravitational constant G, dt, initial velocities and positions :

#PARAMETERS

#Masses
m1 = 3 
m2 = 3
m3 =  3

e1 = 0.1 #Escale factor
e2 = 0.1 #Escale factor
e3 = 0.1 #Escale factor

#The radius as a funtion of mass (volume)
r1 = pow(m1, 1/3)*e1 
r2 = pow(m2, 1/3)*e2
r3 =  pow(m3, 1/3)*e3
G = 3
dt =0.01

#Initial positions
pos1 = vector(-3, 0, 0) #red ball
pos2 = vector(0, 0, 0) #green ball
pos3 =vector(3, 0, 0) #blue ball

#velocitats inicials
v1= vector(0, 0, 1.5) #red
v2 = vector(0, 1.5, 0) #green
v3 = vector(0, 0, -1.5) #blue


Done that it's time to start drawing the objects as spheres and their respective trail to see the trajectory:
#DRAWING
P1 = sphere(pos = pos1, radius = r1, color = color.red)
P2 = sphere(pos = pos2, radius = r2, color = color.green)
P3 = sphere(pos = pos3, radius = r3, color = color.blue)
P1.trail = curve(color = P1.color)
P2.trail = curve(color = P2.color)
P3.trail = curve(color = P3.color)


Now comes the movement calculations, when physics come into. As there are three objects I made a function called "grav" to calculate the acceleration due to any of the three spheres. The needed input variables for that function are the position of the body that it's attracting (p_m) and the position of the attracted (p) and it's mass (m):

#MOVEMENT
#Defining an acceleration function:
def grav(p, p_m, m):
    r = p-p_m
    r_mag = mag(r)
    r_norm = norm(r)
    a = (-G*m/(r_mag**2)) * r_norm
    return a


Once defined the functions lets start with the while loop. Basically what it does here is to calculate the position of the three balls by "differential equations":

#Computing
while True:
    rate(60)
    #P1
    a1 = grav(P1.pos, P2.pos, m2) + grav(P1.pos, P3.pos, m3)
    dv1 = a1*dt
    v1 = v1 + dv1
    dp1 = v1*dt
    #P2
    a2 = grav(P2.pos, P1.pos, m1) + grav(P2.pos, P3.pos, m3)
    dv2 = a2*dt
    v2 = v2 + dv2
    dp2 = v2*dt
    #P3
    a3 = grav(P3.pos, P2.pos, m2) + grav(P3.pos, P1.pos, m1)
    dv3 = a3*dt
    v3 = v3 + dv3
    dp3 = v3*dt


Then, Without leaving the while loop, it's drown the new position for planets and it's updated the trail position. Notice that as the equations are written as a function of the sphere's position would be an error to update the position (i.e, P1.pos = P1.pos + dp1) just after calculating it, as the other planets would ask for it and the correct position would be the after-updating one!

As all positions for this iteration are calculated lets apply the changes:
    #Updating positions
    P1.pos = P1.pos + dp1
    P1.trail.append(pos=P1.pos)
    P2.pos = P2.pos + dp2
    P2.trail.append(pos=P2.pos)
    P3.pos = P3.pos + dp3
    P3.trail.append(pos=P3.pos)
    scene.center = P2.pos

And finally an "if" condition to break the while loop if any ball gets to many close to another one:
    #Colision condition
    if  mag(P1.pos-P2.pos)<(r1+r2) or mag(P1.pos-P3.pos)<(r1+r3) or (mag(P3.pos-P2.pos)<(r3+r2)):
        print "Collision !?"
        break

So that's all, ask anything you didn't understand!
Have fun!

You can download the script from here.