#!/usr/bin/python3
# ===============================================================
# From: www.youtube.com/watch?v=bD8E-hWY9Yk
# Astrophysics: N-Body Simulation in Python
# ===============================================================
#
# Gravity (or the acceleration due to gravity) is 9.81 meters
# per second squared (32 feet per second squared), on the
# surface of Earth, because of the size of Earth and the
# distance we are on its surface from its center.
#
# The acceleration due to gravity on Earth is not a constant.
# It varies from place to place and from altitude to altitude.
# The standard value at sea level is 9.80665 m/s². That's still
# not considered a constant.
#
# ===============================================================
#
# VPython math functions and vectors documentation:
# www.glowscript.org/docs/VPythonDocs/math.html
#
# ===============================================================
from vpython import *
# ---- set every thing to 1 to simplify things
# ---- (this is only a simulation demo)
G = 1 # gravitational constant
M = 1 # total mass of all stars
R = 1 # radius of universe
# ---- create N stars at random locations (between -1 to +1)
n = 0 # loop counter
N = 5 # number of stars
stars = [] # list of stars
while n < N:
# ---- random vector (rt)
rt = R*vector(2*random()-1,2*random()-1,2*random()-1)
# ---- create sphere (coordinates,radius, ...)
stars = stars + [sphere(pos=rt,radius=R/30,
make_trail=True,retain=100)]
n += 1
# ---- display stars (coordinates,radius)
##for star in stars: print(star.pos,star.radius)
# ---- add mass, velocity, and net force on each star
# ---- m = star's mass (all stars have the same mass)
# ---- p = star's momentum
# ---- F = net force on the star
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)
# ---- run the simulation
t = 0 # start time
dt = 0.01 # delta time (time step)
while t < 10: # 10 seconds
# ---- animation rate
rate(100)
# ---- initial net force vector for each star (0,0,0)
# ---- (each time through the loop recalculate net force)
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 and position
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