Wednesday, March 18, 2015


Here's a python script "" that I wrote based on the "bouncing balls" demo at:

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("")
# ..."ipython --pylab" should do the following, according to
# 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
#axis([-10, 10, -10, 10])
# Scale the axes:
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)
    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!

1 comment:

  1. The custom essay writing service is now an important task for college students. They have to write many essays in their academic years. But some lengthier essay task will make students in a lazy mood write their future essays.