# some subroutines
# Compute the coefficient of restitution for a pair of particles bouncing off each other
# as a function of the velocity difference v between the two objects
# and r the distance between them (at impact)
# units: v [m/s] and r [m]
# returns the coefficient of restitution eps
# which is probably defined as the ratio of post and pre relative velocity magnitudes
def cor_bridges(r, v):
v32 = 0.01 # m/s velocity giving eps = 0.32, to increase damping lower this number
eps = 0.32*pow(abs(v)/v32 + 1e-6,-0.234) # high velocity encounters have lower coefficents of restitution
if eps>1.: # keep at or below 1
eps=1.
if eps<0.01:
eps=0.01 # always bounce a little bit
return eps
# eps is 0.32 at v = 0.01 m/s, eps is lower than this at higher velocities
# You could modify this routine, or create an extra similar routine if you want to adjust the dissipation rate
# return a random variable that is drawn from a probability distribution
# that is a power-law with probability p propto v^slope
# and lies within [min_v,max_v]
def powerlaw(slope, min_v, max_v):
y = np.random.uniform()
pow_max = pow(max_v, slope+1.)
pow_min = pow(min_v, slope+1.)
return pow((pow_max-pow_min)*y + pow_min, 1./(slope+1.))
# make a pretty plot of simulation particle positions in x,y plotting circles with their radii
# This is the same as Hanno Rein's function in the Saturn example, but here x and y are not flipped
def plotParticles(sim):
fig = plt.figure(figsize=(8,8))
ax = plt.subplot(111,aspect='equal')
ax.set_xlabel("radial coordinate [m]")
ax.set_ylabel("azimuthal coordinate [m]")
boxsize = sim.boxsize.x
ax.set_ylim(-boxsize/2.,boxsize/2.)
ax.set_xlim(-boxsize/2.,boxsize/2.)
for i, p in enumerate(sim.particles): # loop over particles
circ = patches.Circle((p.x, p.y), p.r, facecolor='darkgray', edgecolor='black')
ax.add_patch(circ) # plot a circle
# notice that p.x and p.y were flipped in the plotting routine by Hanno
# plot some arrows
ax.arrow( boxsize*0.5*0.9,0,0,-boxsize/4,width=2) # to show direction of motion
ax.arrow(-boxsize*0.5*0.9,0,0, boxsize/4,width=2)
# create an initial distribution of particles
# the particles are added to the rebound simulation sim
# surface density [kg m^-2] and particle density [kg m^-3] are the total mass surface density and
# the particle mass density
# particles are added until the expected total mass in the box is reached
# the radius distribution of the particles is a powerlaw of slope = slope
# min and max radii for this distribution are min_v and max_v in m
# The mass distribution has powerlaw slope of slope/3?
# I modified this routine so each particle is asigned a unique hash number
# we need to do this if we want to id particles as the tree routine reaaranges the indices
# To adjust the mass/radius distributions you would vary the max or min radii or slope
def mkparticles(surface_density,particle_density,slope,min_v,max_v,sim):
boxsize = sim.boxsize.x
total_mass = 0.
i=0
while (total_mass < surface_density*(boxsize**2)):
radius = powerlaw(slope=slope, min_v=min_v, max_v=max_v) # [m]
mass = particle_density*4./3.*np.pi*(radius**3) # [kg m^-3]
x = np.random.uniform(low=-boxsize/2., high=boxsize/2.) # a random number in [-boxsize/2,boxside/2]
sim.add(
hash = i, m=mass, r=radius,
x=x,
y=np.random.uniform(low=-boxsize/2., high=boxsize/2.),
z=np.random.normal(),
vx = 0.,
vy = -3./2.*x*sim.ri_sei.OMEGA,
vz = 0.)
total_mass += mass
i+= 1
# you can see shear direction here from how vy depends upon x
# particles are begun at zero epicyclic amplitude
# the shear motion is down on the right (+x) and up on the left (-x)
# the (unique) hash number is set equal to the index of the loop
# the z value is a normal distribution with dispersion of 1 m
print('Number of particles= ',sim.N);
# compute residual (remainder) of x modulo m but in range [-m/2,m/2]
# this routine is useful for keep quantities within the periodic shearing box
def resid(x,m):
z = (x+ m/2)%m - m/2
return z
