#!/usr/bin/python3
# ===============================================================
# https://www.youtube.com/watch?v=bD8E-hWY9Yk
# Astrophysics: N-Body Simulation in Python
# ===============================================================
from vpython import *
# ---- set every thing to 1 to simplify things
G = 1 #
M = 1 # total mass of all stars
R = 1 #
# ---- create N random stars
# ---- (location is -1 to +1)
n = 0
N = 5
stars = []
while n < N:
rt = R*vector(2*random()-1,2*random()-1,2*random()-1)
stars = stars + [sphere(pos=rt,radius=R/30,
make_trail=True,retain=100)]
n += 1
# ---- display stars (coords,raidus,make_trail)
##for star in stars: print(star.pos)
# ---- add mass, velocity, and net force on each star
# ---- total star mass (M=1) is the same no mater how many stars
for star in stars:
star.m = M/N
star.p = star.m*vector(0.1*(2*random()-1),0,0)
star.F = vector(0,0,0)
t = 0 # start time
dt = 0.01 # delta time (time step)
while t < 10: # 10 seconds
# ---- number of calculations per second
rate(100)
# ---- initial force vector for each star
for star in stars: star.F = vector(0,0,0)
# ---- take each star and calculate how it interacts with
# ---- the other star (calculate the net force on
# ---- each star)
# ---- ------------------------------------------------------
# ---- note
# ---- the force between A,B is the same as as the force
# ---- between B,A. we calculate it twice anyway to
# ---- simplify the code and avoid a lot of if tests.
for i in range(len(stars)):
for j in range(len(stars)):
if i != j:
rji = stars[i].pos - stars[j].pos
stars[i].F = stars[i].F - \
G * stars[i].m * stars[j].m * norm(rji)/mag(rji)
# ---- update each star's momentum
for star in stars:
star.p = star.p + star.F*dt
star.pos = star.pos + star.p*dt/star.m
# ---- next time step
t = t + dt