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!