https://python4astronomers.github.io/contest/bounce.html
This script qualitatively illustrates the gravitational motion of massless test particles in response to two uniform-density spherical potentials of radius R=1 and mass M=1. In units where G=1 as well, the acceleration outside the spheres is 1/d^2, and within the spheres is ((d/R)^3)/d^2 = d/R^3 = d.
Notice that the structures which have collapsed by the end of the script are elongated as expected.
Updated: script below runs in ipython.
# run using:
# ipython --pylab
# execfile("gravitate.py")
#
# ..."ipython --pylab" should do the following, according to http://ipython.org/ipython-doc/2/api/generated/IPython.core.magics.pylab.html?highlight=pylab
# import numpy as np
# from matplotlib import pyplot as plt
# from matplotlib import pylab, mlab
# from IPython.display import display
# from IPython.core.pylabtools import figsize, getfigs
#
figure(1) # change number to create new fig instead of overwritingclf()
#axis([-10, 10, -10, 10])
# Scale the axes:
siz=2.5
axis([-siz*10, siz*10, -siz*10, siz*10])
plt.gca().set_aspect('equal', adjustable='box')
# Define properties of the "gravitating spheres"
n = 999
# Factor 1.5 means starting positions can be off screen
pos = 1.5*siz*(20 * random_sample(n*2) - 10).reshape(n, 2)
vel = (0.003 * normal(size=n*2)).reshape(n, 2)
radius = 10.0
# radius=0.10 # uncomment to plot particle tracks
sizes = radius * random_sample(n) + radius
# Colors where each row is (Red, Green, Blue, Alpha). Each can go
# from 0 to 1. Alpha is the transparency.
colors = random_sample([n, 4])
# Draw all the circles and return an object ``circles`` that allows
# manipulation of the plotted circles.
circles = scatter(pos[:,0], pos[:,1], marker='o', s=sizes, c=colors)
# Gravitational clustering of massless particles near two mass overdensities
# located at lower left (ll) and upper right (ur) of figure area
ten=7.5 #10.0 # initial overdensity coords: -ten,-ten and +ten,+ten
hlf=1.0 #1.0 # scale factor for gravitational force inside overdensity
scl=0.1 #0.1 # scale factor for overall acceleration on test particles
for i in range(200):
ten=ten-scl*sqrt(hlf)/ten**2.0 # move the overdensities closer over time
# Renormalize positions from (-ten,+ten) to (0,2)*sqrt(hlf)
posll=(pos+ten)*sqrt(hlf)/10. ; posur=(ten-pos)*sqrt(hlf)/10. # -posll
# Compute distances in renormalized units to overdensity centres
dll = sqrt(posll[:,0]**2.+posll[:,1]**2.)
dur = sqrt(posur[:,0]**2.+posur[:,1]**2.)
# Set flags: dll1=1 for dll>1, dll0=1 for dll<1; same for ur
dll1 = 0.5*(1.0+(dll-1.0)/abs(dll-1.0)) ; dll0 = 1.0-dll1
dur1 = 0.5*(1.0+(dur-1.0)/abs(dur-1.0)) ; dur0 = 1.0-dur1
# At d<1, x acceleration \propto (x/d)*(rho*d**3.)*(1/d**2) = x
# At d>1, x acceleration \propto (x/d)*(1/d**2) = x/d**3
# ...and similarly for y acceleration.
# Use flags to set one acceleration term to zero for each sphere
xllaccel = (-posll[:,0]*(dll0+sqrt(hlf)*dll1/(dll)**3.0)).reshape(n,1)
yllaccel = (-posll[:,1]*(dll0+sqrt(hlf)*dll1/(dll)**3.0)).reshape(n,1)
xuraccel = (posur[:,0]*(dur0+sqrt(hlf)*dur1/(dur)**3.0)).reshape(n,1)
yuraccel = (posur[:,1]*(dur0+sqrt(hlf)*dur1/(dur)**3.0)).reshape(n,1)
#print(pos,xuraccel,yuraccel)
llaccel = np.append(xllaccel,yllaccel,axis=1)
uraccel = np.append(xuraccel,yuraccel,axis=1)
pos = pos + vel + scl*(llaccel+uraccel)
#bounce = abs(pos) > 10*siz # Find balls that are outside walls
#vel[bounce] = -vel[bounce] # Bounce if outside the walls
circles.set_offsets(pos) # Change the positions
draw() # comment out this line & uncomment next line & above to plot tracks
#scatter(pos[:,0], pos[:,1], marker='o', s=sizes, c=colors) # will be slow!
