from vpython import *
# Hard-sphere gas.
# Bruce Sherwood
win = 400
Natoms = 100 # change this to have more or fewer atoms
# Typical values
L = 1 # container is a cube L on a side
gray = color.gray(0.7) # color of edges of container
mass = 4E-3/6E23 # helium mass
Ratom = 0.03 # wildly exaggerated size of helium atom
k = 1.4E-23 # Boltzmann constant
T = 300 # around room temperature
dt = 1E-5
animation = canvas( width=win, height=win, align='left')
d = L/2+Ratom
r = 0.005
boxbottom = curve(color=gray, radius=r)
boxbottom.append(pos=[vector(-d,-d,-d), vector(-d,-d,d), vector(d,-d,d), vector(d,-d,-d), vector(-d,-d,-d)])
boxtop = curve(color=gray, radius=r)
boxtop.append([vector(-d,d,-d), vector(-d,d,d), vector(d,d,d), vector(d,d,-d), vector(-d,d,-d)])
vert1 = curve(color=gray, radius=r)
vert2 = curve(color=gray, radius=r)
vert3 = curve(color=gray, radius=r)
vert4 = curve(color=gray, radius=r)
vert1.append([vector(-d,-d,-d), vector(-d,d,-d)])
vert2.append([vector(-d,-d,d), vector(-d,d,d)])
vert3.append([vector(d,-d,d), vector(d,d,d)])
vert4.append([vector(d,-d,-d), vector(d,d,-d)])
Atoms = []
apos = []
p = []
pavg = sqrt(2*mass*1.5*k*T) # average kinetic energy p**2/(2mass) = (3/2)kT
r = 0.005*L
for i in range(Natoms):
x = L*random()-L/2
y = L*random()-L/2
z = L*random()-L/2
if i == 0:
Atoms.append(sphere(pos=vector(x,y,z), radius=Ratom, color=color.cyan,
make_trail=True, retain=100, trail_radius=0.3*Ratom))
else: Atoms.append(sphere(pos=vector(x,y,z), radius=Ratom, color=gray))
apos.append(vector(x,y,z))
theta = pi*random()
phi = 2*pi*random()
px = pavg*sin(theta)*cos(phi)
py = pavg*sin(theta)*sin(phi)
pz = pavg*cos(theta)
p.append(vector(px,py,pz))
animation.title = 'A "hard-sphere" gas'
s = ' Theoretical and averaged speed distributions (meters/sec).\n'
s += ' Initially all atoms have the same speed, but collisions\n'
s += ' change the speeds of the colliding atoms. One of the atoms is\n'
s += ' marked and leaves a trail so you can follow its path.'
animation.caption = s
animation.range = L
deltav = 100 # binning for v histogram
def barx(v):
return int(v/deltav) # index into bars array
nhisto = int(4000/deltav)
histo = []
for i in range(nhisto): histo.append(0.0)
histo[barx(pavg/mass)] = Natoms
graph( width=500, height=300, xmax=3000, ymax=Natoms*deltav/1000, align='right' )
theory = gcurve( color=color.cyan )
dv = 10
for v in range(0,3001+dv,dv): # theoretical prediction
theory.plot( v, (deltav/dv)*Natoms*4*pi*((mass/(2*pi*k*T))**1.5) *exp(-0.5*mass*(v**2)/(k*T))*(v**2)*dv )
accum = []
for i in range(int(3000/deltav)): accum.append([deltav*(i+.5),0])
vdist = gvbars(color=color.red, delta=deltav )
def interchange(v1, v2): # remove from v1 bar, add to v2 bar
barx1 = barx(v1)
barx2 = barx(v2)
if barx1 == barx2: return
histo[barx1] -= 1
histo[barx2] += 1
def checkCollisions(apos):
hitlist = []
Natoms = len(apos)
r2 = 2*Ratom
r2 *= r2
for i in range(Natoms):
ai = apos[i]
for j in range(i+1, Natoms) :
aj = apos[j]
dr = ai - aj
if mag2(dr) < r2: hitlist.append([i,j])
return hitlist
nhisto = 0 # number of histogram snapshots to average
while True:
rate(150)
# Accumulate and average histogram snapshots
for i in range(len(accum)): accum[i][1] = (nhisto*accum[i][1] + histo[i])/(nhisto+1)
if nhisto % 10 == 0:
vdist.data = accum
nhisto += 1
# Update all positions
for i in range(Natoms): Atoms[i].pos = apos[i] = apos[i] + (p[i]/mass)*dt
# Check for collisions
hitlist = checkCollisions(apos)
# If any collisions took place, update momenta of the two atoms
for ij in hitlist:
i = ij[0]
j = ij[1]
ptot = p[i]+p[j]
vi = p[i]/mass
vj = p[j]/mass
vrel = vj-vi
a = vrel.mag2
if a == 0: continue; # exactly same velocities
rrel = apos[i]-apos[j]
b = 2*rrel.dot(vrel)
c = rrel.mag2-4*Ratom*Ratom
d = b*b-4*a*c
if d < 0: continue # something wrong; ignore this rare case
deltat = (-b+sqrt(d))/(2*a) # t-deltat is when they made contact
apos[i] = apos[i]-vi*deltat # back up to contact configuration
apos[j] = apos[j]-vj*deltat
mtot = 2*mass
pcmi = p[i]-ptot*mass/mtot # transform momenta to cm frame
pcmj = p[j]-ptot*mass/mtot
rrel = norm(rrel)
pcmi = pcmi-2*pcmi.dot(rrel)*rrel # bounce in cm frame
pcmj = pcmj-2*pcmj.dot(rrel)*rrel
p[i] = pcmi+ptot*mass/mtot # transform momenta back to lab frame
p[j] = pcmj+ptot*mass/mtot
apos[i] = apos[i]+(p[i]/mass)*deltat # move forward deltat in time
apos[j] = apos[j]+(p[j]/mass)*deltat
interchange(vi.mag, p[i].mag/mass)
interchange(vj.mag, p[j].mag/mass)
for i in range(Natoms):
loc = apos[i]
if abs(loc.x) > L/2:
if loc.x < 0: p[i].x = abs(p[i].x)
else: p[i].x = -abs(p[i].x)
if abs(loc.y) > L/2:
if loc.y < 0: p[i].y = abs(p[i].y)
else: p[i].y = -abs(p[i].y)
if abs(loc.z) > L/2:
if loc.z < 0: p[i].z = abs(p[i].z)
else: p[i].z = -abs(p[i].z)