Gravitate

Here's a python script "gravitate.py" that I wrote based on the "bouncing balls" demo at:
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 overwriting
clf()
#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!