N-body Simulation

What are N-body simulations?

According to physics.princeton.edu, an N-body simulation replicates the motion of particles that interact with one another through some type of physical forces. These particles can range from celestial bodies to individual atoms or molecules, but for our code, we will define the particles as physical celestial bodies.

Here is the code for Python:

3-Body Simulation:

from visual import*
 R= 6.378e6
 G=6.67e-11

M1= 6.0e24
 M2= M1
 M3= M1

r1= vector(0,0,0)
 r2= vector(-10*R,0,0)
 r3= vector(10*R,0,0)

object1= sphere(pos=(r1),
 radius=R,
 color=color.white,
 make_trail=True,
 trail_type="points",
 interval=10)

object2= sphere(pos=(r2),
 radius=R,
 color=color.red,
 make_trail=True,
 trail_type="points",
 interval=10)

object3= sphere(pos=(r3),
 radius=R,
 color=color.blue,
 make_trail=True,
 trail_type="points",
 interval=10)

v0= sqrt(2*G*M1/R)

v1= vector(0,v0,0)
 v2= v1
 v3= v1

v_characteristic= sqrt(G*M1/mag(r1-r2))
 T = mag(r1-r2)/v_characteristic

dt = 0.001*T

p_object1=M1*v1
 p_object2=M2*v2
 p_object3=M3*v3

t=0
 small= 0.1*R

while t< 50*T:
 rate(100)

#1 on 2
 r12 = mag(r1 - r2)
 if r12<small:
 r12=small

F12= G*M1*M2/(r12*r12)*norm(r1-r2)
 F21 = -F12

p_object1=p_object1+F21*dt
 p_object2=p_object2+F12*dt

v1=p_object1/M1
 v2=p_object2/M2

r1= r1+v1*dt
 r2= r2+v2*dt

object1.pos= object1.pos+v1*dt
 object2.pos= object2.pos+v2*dt

#1 on 3
 r13 = mag(r1- r3)
 if r13<small:
 r13=small

F13= G*M1*M3/(r13*r13)*norm(r3-r1)
 F31 = -F13

p_object1=p_object1+F31*dt
 p_object3=p_object3+F13*dt

v1=p_object1/M1
 v3=p_object3/M3

r1= r1+v1*dt
 r3= r3+v3*dt

object1.pos= object1.pos+v1*dt
 object3.pos= object3.pos+v3*dt

#2 on 3
 r23 = mag(r2- r3)
 if r23<small:
 r23=small

F23= G*M2*M3/(r23*r23)*norm(r3-r2)
 F32 = -F23

p_object2=p_object2+F32*dt
 p_object3=p_object3+F23*dt

v2=p_object2/M2
 v3=p_object3/M3

r2= r2+v2*dt
 r3= r3+v3*dt

object2.pos= object2.pos+v2*dt
 object3.pos= object3.pos+v3*dt

t=t+dt

4-Body Simulation:

from visual import*
 R= 6.378e6
 G=6.67e-11

M1= 6.0e24
 M2= M1
 M3= M1
 M4= M1

r1= vector(10*R,-10*R,0)
 r2= vector(20*R,0,0)
 r3= vector(0,0,0)

object1= sphere(pos=(r1),
 radius=R,
 color=color.white,
 make_trail=True,
 trail_type="points",
 interval=10)

object2= sphere(pos=(r2),
 radius=R,
 color=color.yellow,
 make_trail=True,
 trail_type="points",
 interval=10)

object3= sphere(pos=(r3),
 radius=R,
 color=color.red,
 make_trail=True,
 trail_type="points",
 interval=10)

v1= vector(0,0,0)
 v2= v1
 v3= v1

v_characteristic= sqrt(G*M1/mag(r1-r2))
 T = mag(r1-r2)/v_characteristic

dt = 0.0001*T

p_object1=M1*v1
 p_object2=M2*v2
 p_object3=M3*v3

t=0
 small= 0.5*R

while t< 50*T:
 rate(10000)

r12 = mag(r2- r1)
 r13 = mag(r3- r1)
 r23 = mag(r3- r2)

if r12<small:
 r12=small

if r13<small:
 r13=small

if r23<small:
 r23=small

F12= G*M2*M1/(r12*r12)*norm(r1-r2)
 F21 = -F12
 F13= G*M3*M1/(r13*r13)*norm(r1-r3)
 F31 = -F13
 F23= G*M3*M2/(r23*r23)*norm(r2-r3)
 F32 = -F23

F1= F21+F31
 F2= F12+F32
 F3= F23+F13

p_object1=p_object1+F1*dt
 v1=p_object1/M1
 r1= r1+v1*dt
 object1.pos= r1

p_object2=p_object2+F2*dt
 v2=p_object2/M2
 r2= r2+v2*dt
 object2.pos= r2

p_object3= p_object3+F3*dt
 v3=p_object3/M3
 r3= r3+v3*dt
 object3.pos= r3

t=t+dt

5-Body Simulation:

from visual import*
 R= 6.378e6
 G=6.67e-11

M1= 6.0e24
 M2= M1
 M3= M1
 M4= M1
 M5= M1

r1= vector(10*R,-5*R,0)
 r2= vector(15*R,0,0)
 r3= vector(0,0,0)
 r4= vector(7*R,8*R,R)
 r5= vector(8*R,10*R,0)

object1= sphere(pos=(r1),
 radius=R,
 color=color.white)

object2= sphere(pos=(r2),
 radius=2*R,
 color=color.blue)

object3= sphere(pos=(r3),
 radius=R,
 color=color.cyan)

object4= sphere(pos=(r4),
 radius=4*R,
 color=color.red)

object5= sphere(pos=(r5),
 radius=R,
 color=color.magenta)

v1= vector(0,0,0)
 v2= v1
 v3= v1
 v4= v1
 v5= v1

v_characteristic= sqrt(G*M1/mag(r1-r2))
 T = mag(r1-r2)/v_characteristic

dt = 0.0001*T

p_object1=M1*v1
 p_object2=M2*v2
 p_object3=M3*v3
 p_object4=M4*v4
 p_object5=M5*v5




t=0
 small= 0.1*R

while t< 500*T:
 rate(10000)

r12 = mag(r2- r1)
 r13 = mag(r3- r1)
 r14 = mag(r4- r1)
 r15 = mag(r5- r1)
 r23 = mag(r3- r2)
 r24 = mag(r4- r2)
 r25 = mag(r5- r2)
 r34 = mag(r4- r3)
 r35 = mag(r5- r3)
 r45 = mag(r5- r4)

if r12<small:
 r12=small

if r13<small:
 r13=small

if r14<small:
 r14=small

if r15<small:
 r15=small

if r23<small:
 r23=small

if r24<small:
 r24=small

if r25<small:
 r25=small

if r34<small:
 r34=small

if r35<small:
 r35=small

if r45<small:
 r45=small

F12= G*M2*M1/(r12*r12)*norm(r1-r2)
 F21 = -F12

F13= G*M3*M1/(r13*r13)*norm(r1-r3)
 F31 = -F13

F14= G*M1*M4/(r14*r14)*norm(r1-r4)
 F41 = -F14

F15= G*M1*M5/(r15*r15)*norm(r1-r5)
 F51 = -F15

F23= G*M2*M3/(r23*r23)*norm(r2-r3)
 F32 = -F23

F24= G*M2*M4/(r24*r24)*norm(r2-r4)
 F42 = -F24

F25= G*M2*M5/(r25*r25)*norm(r2-r5)
 F52 = -F25

F34= G*M3*M4/(r34*r34)*norm(r3-r4)
 F43 = -F34

F35= G*M3*M5/(r35*r35)*norm(r3-r5)
 F53 = -F35

F45= G*M4*M5/(r45*r45)*norm(r4-r5)
 F54 = -F45




F1= F21+F31
 F2= F12+F32
 F3= F23+F13
 F4= F14+F34
 F5= F15+F45

p_object1=p_object1+F1*dt
 v1=p_object1/M1
 r1= r1+v1*dt
 object1.pos= r1

p_object2=p_object2+F2*dt
 v2=p_object2/M2
 r2= r2+v2*dt
 object2.pos= r2

p_object3= p_object3+F3*dt
 v3=p_object3/M3
 r3= r3+v3*dt
 object3.pos= r3

p_object4= p_object4+F4*dt
 v4=p_object4/M4
 r4= r4+v4*dt
 object4.pos= r4

p_object5= p_object5+F5*dt
 v5=p_object5/M5
 r5= r5+v5*dt
 object5.pos= r5

t=t+dt