The two-body problem consists of two objects under gravitational attraction. The attraction force is the newtonian force G*M1*M2/distance^2, where G is some constant, M is mass and distance is the distance between the two objects.
In the two-body problem, the path of each particle can be solved for explicitly. We can get circles or ellipses, like the earth's path around the sun, or we can get escape trajectories where the two particles pass by and never collide with each other again, taking the form of hyperbolas and parabolas.
The three-body problem is more difficult though. Three-bodies can cause chaos, and you can't solve it explicitly for hyperbolas/parabolas. The only way to solve the general case is to simulate the system over time.
It's hard to get to grips with the chaos of the three body problem, but this program tries to depict it.
It starts out with three bodies with the same mass/radius and sends them off in different velocities. From this, you can calculate the future path of the particles. This program shows the path of three objects over the full simulation, so instead of viewing three circles (the particles current positions), you see three lines (the particles paths given the starting conditions). As the initial velocity varies over time, the paths vary, and you get chaotic looking results.
When two particles collide, they're combined into the same particle, so this explains when paths combine into one line. It also explains why the problem is chaotic - if two particles are going to collide, a small nudge one way or the other could lead to a whole alternate timeline where they don't collide.
`grav
`grav.dba
`======================
type circ
radius as float
x as float //position
y as float
xv as float //velocity
yv as float
xf as float //force (mass*acceleration)
yf as float
mass as float
col as dword
removeme as integer
endtype
type screendot
x as integer
y as integer
col as dword
endtype
global w as integer
global h as integer
global numcirc as integer
global numdots as integer
global springK as float
global imagenum as integer
imagenum=1
w=screen width()
h=screen height()
numcirc=-1
numdots=-1
dim c(numcirc) as circ
dim d(numdots) as screendot
dim removequeue(-1) as integer
height as float
height=-50
cleardots()
initialize(height)
calculateLocus(200,.05,4)
sync on
sync rate 0
do
clearCircles()
cleardots()
initialize(height)
calculateLocus(4800,.1,4)
deletememimage(1)
makememimage(1,500,500)
drawall()
drawmemimage(1)
sync
height=height+.1
if height>50 then end
remstart
for m=0 to 40
handleCollision()
integrate(.01)
next m
updateDots()remend
loop
function getImageName()
if imagenum<10
ret$="000"+str$(imagenum)
exitfunction ret$
endif
if imagenum<100
ret$="00"+str$(imagenum)
exitfunction ret$
endif
if imagenum<1000
ret$="0"+str$(imagenum)
exitfunction ret$
endif
ret$=str$(imagenum)
endfunction ret$
function makeMemImage(memnum as integer,xs as integer, ys as integer)
pixels=xs*ys
bytes=12+4*pixels
make memblock memnum,bytes
write memblock dword memnum,0,xs
write memblock dword memnum,4,ys
write memblock dword memnum,8,32
endfunction
function deleteMemImage(memnum as integer)
if memblock exist(memnum) then delete memblock memnum
endfunction
function memDot(memnum as integer,xs as integer, x as integer, y as integer, color as dword)
if x>=memblock dword(memnum,0) or y>=memblock dword(memnum,4) or x<0 or y<0 then exitfunction
write memblock dword memnum,(x+xs*y)*4+12,color||(255>>24)
endfunction
function drawMemImage(memnum as integer)
if image exist(1) then delete image 1
make image from memblock 1,memnum
paste image 1,0,0
save image getImageName()+".bmp",1
inc imagenum
endfunction
function clearDots()
undim d()
numdots=-1
dim d(numdots) as screendot
endfunction
function clearCircles()
undim c()
numcirc=-1
dim c(numcirc) as circ
endfunction
function zeroMomentum()
xm as float
ym as float
m as float
xm=0
ym=0
m=0
for n=0 to numcirc
inc xm,c(n).xv*c(n).mass
inc ym,c(n).yv*c(n).mass
inc m,c(n).mass
next n
xm=xm/m
ym=ym/m
for n=0 to numcirc
dec c(n).xv,xm
dec c(n).yv,ym
next
endfunction
function initialize(height as float)
springK=1
//setCircle(0,w/2-200,h/2-height,0,.4,30,1,rgb(255,0,0))
//setCircle(1,w/2,h/2,0,-.4,30,1,rgb(0,255,0))
//setCircle(2,w/2+200,h/2,.4,-.6,30,1,rgb(0,255,0))
addCircle(w/2-200,h/2,.1,cos(height*180/50.0),30,1,rgb(255,0,0))
addCircle(w/2,h/2,0,-.9,30,1,rgb(0,255,0))
addCircle(w/2-100,h/2-100,.4,.6,30,1,rgb(0,0,255))
zeroMomentum()
endfunction
function commitRemovals()
for n=0 to array count(c())
if c(n).removeme=1 then removeCircle(n)
next n
endfunction
function removeCircle(num as integer)
array delete element c(),num
dec numcirc
endfunction
function addCircle(xpos as float, ypos as float, xvel as float, yvel as float, radius as float, mass as float, color as dword)
array insert at bottom c()
inc numcirc
n=array count(c())
setCircle(n,xpos,ypos,xvel,yvel,radius,mass,color)
endfunction
function calculateLocus(time as float, timestep as float,timeperdot as float)
curtime as float
counter as float
curtime=0
counter=0
while curtime<time
curtime=curtime+timestep
counter=counter+timestep
commitRemovals()
handleGravity()
//handleCollision()
integrate(timestep)
zeroMomentum()
if counter>=timeperdot
counter=0
updateDots()
endif
endwhile
endfunction
function updateDots()
for n=0 to numcirc
addDot(c(n).x,c(n).y,c(n).col)
next n
endfunction
function addDot(x as integer, y as integer, color as dword)
inc numdots
array insert at bottom d()
d(numdots).x=x
d(numdots).y=y
d(numdots).col=color
endfunction
function handleGravity()
for n=0 to numcirc
for m=0 to numcirc
if n<>m and c(n).removeme=0 and c(m).removeme=0
diffx as float
diffy as float
dist as float
overlap as float
forcemag as float
diffx=c(m).x-c(n).x //difference in position
diffy=c(m).y-c(n).y
dist=sqrt(diffx*diffx+diffy*diffy)
if dist<.1*(c(m).radius+c(n).radius) //then combine the two particles
c(m).removeme=1 //queue for removal
//conserve momentum
c(n).xv=(c(n).xv*c(n).mass+c(m).xv*c(m).mass)/(c(n).mass+c(m).mass)
c(n).yv=(c(n).yv*c(n).mass+c(m).yv*c(m).mass)/(c(n).mass+c(m).mass)
//reposition to center of mass
c(n).x=(c(n).x*c(n).mass+c(m).x*c(m).mass)/(c(n).mass+c(m).mass)
c(n).y=(c(n).y*c(n).mass+c(m).y*c(m).mass)/(c(n).mass+c(m).mass)
//add forces already calculated
c(n).xf=c(n).xf+c(m).xf
c(n).yf=c(n).yf+c(m).yf
//aand last but not least, add the masses
c(n).mass=c(n).mass+c(m).mass
else
diffx=diffx/dist
diffy=diffy/dist
forcemag=100.0/(dist*dist)
diffx=diffx*forcemag
diffy=diffy*forcemag
c(n).xf=c(n).xf+diffx
c(n).yf=c(n).yf+diffy
endif
endif
next m
next n
endfunction
function handleCollision()
for n=0 to numcirc
for m=0 to numcirc
if n<>m
diffx as float
diffy as float
dist as float
overlap as float
forcemag as float
diffx=c(m).x-c(n).x //difference in position
diffy=c(m).y-c(n).y
dist=sqrt(diffx*diffx+diffy*diffy)
overlap=c(n).radius+c(m).radius-dist
if overlap>0 //if the sum of the radii is greater than the distance
forcemag=-springK*overlap //simple spring repulsion force
diffx=diffx/dist //normalize the direction vector. direction from n to m
diffy=diffy/dist
diffx=diffx*forcemag
diffy=diffy*forcemag
c(n).xf=c(n).xf+diffx //the other object will loop through and change c(m) by the same amount.
c(n).yf=c(n).yf+diffy //it will be recalculated but oh well.
endif
endif
next m
next n
endfunction
function setCircle(n as integer,xpos as float, ypos as float, xvel as float, yvel as float, radius as float, mass as float, color as dword)
c(n).x=xpos
c(n).y=ypos
c(n).xv=xvel
c(n).yv=yvel
c(n).xf=0
c(n).yf=0
c(n).radius=radius
c(n).mass=mass
c(n).col=color
c(n).removeme=0
endfunction
function drawall()
//for n=0 to numcirc
//ink c(n).col,0
//circle c(n).x,c(n).y,c(n).radius
//next n
for n=0 to numdots
//ink d(n).col,0
//dot d(n).x*.5,d(n).y*.5
memdot(1,500,d(n).x*.5,d(n).y*.5+100,d(n).col)
next n
ink rgb(255,255,255),0
endfunction
function integrate(timestep as float)
for n=0 to numcirc
//assume constant force being applied.
//P=P0+V*time+.5*accel*time^2
c(n).x=c(n).x+(c(n).xv+.5*c(n).xf/c(n).mass*timestep)*timestep
c(n).y=c(n).y+(c(n).yv+.5*c(n).yf/c(n).mass*timestep)*timestep
c(n).xv=c(n).xv+c(n).xf/c(n).mass*timestep
c(n).yv=c(n).yv+c(n).yf/c(n).mass*timestep
c(n).xf=0 //and reset the force for updating next time.
c(n).yf=0
next n
endfunction
