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

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.

Thank you for posting this code. I am trying to understand the units on the variables. What are the units in the gravitational constant, G, so that its value is 3? Also, what are the units on the masses, initial positions, and velocities?

ReplyDeleteHi Roy, anybody answer you to this post? I have the same question...

DeleteThanks!

I'm pretty sure they're just random. As its so difficult to get stable orbits, people will just play around with the constants to see what works.

ReplyDelete