{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "application/javascript": [ "require.undef(\"nbextensions/vpython_libraries/glow.min\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require.undef(\"nbextensions/vpython_libraries/glowcomm\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");})" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");})" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");})" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from vpython import *\n", "\n", "# Hard-sphere gas.\n", "\n", "# Bruce Sherwood\n", "\n", "win = 400\n", "\n", "Natoms = 100 # change this to have more or fewer atoms\n", "\n", "# Typical values\n", "L = 1 # container is a cube L on a side\n", "gray = color.gray(0.7) # color of edges of container\n", "mass = 4E-3/6E23 # helium mass\n", "Ratom = 0.03 # wildly exaggerated size of helium atom\n", "k = 1.4E-23 # Boltzmann constant\n", "T = 300 # around room temperature\n", "dt = 1E-5\n", "\n", "animation = canvas( width=win, height=win, align='left')\n", "\n", "d = L/2+Ratom\n", "r = 0.005\n", "boxbottom = curve(color=gray, radius=r)\n", "boxbottom.append(pos=[vector(-d,-d,-d), vector(-d,-d,d), vector(d,-d,d), vector(d,-d,-d), vector(-d,-d,-d)])\n", "boxtop = curve(color=gray, radius=r)\n", "boxtop.append([vector(-d,d,-d), vector(-d,d,d), vector(d,d,d), vector(d,d,-d), vector(-d,d,-d)])\n", "vert1 = curve(color=gray, radius=r)\n", "vert2 = curve(color=gray, radius=r)\n", "vert3 = curve(color=gray, radius=r)\n", "vert4 = curve(color=gray, radius=r)\n", "vert1.append([vector(-d,-d,-d), vector(-d,d,-d)])\n", "vert2.append([vector(-d,-d,d), vector(-d,d,d)])\n", "vert3.append([vector(d,-d,d), vector(d,d,d)])\n", "vert4.append([vector(d,-d,-d), vector(d,d,-d)])\n", "\n", "Atoms = []\n", "apos = []\n", "p = []\n", "pavg = sqrt(2*mass*1.5*k*T) # average kinetic energy p**2/(2mass) = (3/2)kT\n", "r = 0.005*L\n", " \n", "for i in range(Natoms):\n", " x = L*random()-L/2\n", " y = L*random()-L/2\n", " z = L*random()-L/2\n", " if i == 0:\n", " Atoms.append(sphere(pos=vector(x,y,z), radius=Ratom, color=color.cyan,\n", " make_trail=True, retain=100, trail_radius=0.3*Ratom))\n", " else: Atoms.append(sphere(pos=vector(x,y,z), radius=Ratom, color=gray))\n", " apos.append(vector(x,y,z))\n", " theta = pi*random()\n", " phi = 2*pi*random()\n", " px = pavg*sin(theta)*cos(phi)\n", " py = pavg*sin(theta)*sin(phi)\n", " pz = pavg*cos(theta)\n", " p.append(vector(px,py,pz))\n", "\n", "animation.title = 'A \"hard-sphere\" gas'\n", "s = ' Theoretical and averaged speed distributions (meters/sec).\\n'\n", "s += ' Initially all atoms have the same speed, but collisions\\n'\n", "s += ' change the speeds of the colliding atoms. One of the atoms is\\n'\n", "s += ' marked and leaves a trail so you can follow its path.'\n", "animation.caption = s\n", "animation.range = L\n", "\n", "deltav = 100 # binning for v histogram\n", "\n", "def barx(v):\n", " return int(v/deltav) # index into bars array\n", "\n", "nhisto = int(4000/deltav)\n", "histo = []\n", "for i in range(nhisto): histo.append(0.0)\n", "histo[barx(pavg/mass)] = Natoms\n", "\n", "graph( width=500, height=300, xmax=3000, ymax=Natoms*deltav/1000, align='right' )\n", "theory = gcurve( color=color.cyan )\n", "dv = 10\n", "for v in range(0,3001+dv,dv): # theoretical prediction\n", " 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 )\n", "\n", "accum = []\n", "for i in range(int(3000/deltav)): accum.append([deltav*(i+.5),0])\n", "vdist = gvbars(color=color.red, delta=deltav )\n", "\n", "def interchange(v1, v2): # remove from v1 bar, add to v2 bar\n", " barx1 = barx(v1)\n", " barx2 = barx(v2)\n", " if barx1 == barx2: return \n", " histo[barx1] -= 1\n", " histo[barx2] += 1\n", " \n", "def checkCollisions(apos):\n", " hitlist = []\n", " Natoms = len(apos)\n", " r2 = 2*Ratom\n", " r2 *= r2\n", " for i in range(Natoms):\n", " ai = apos[i]\n", " for j in range(i+1, Natoms) :\n", " aj = apos[j]\n", " dr = ai - aj\n", " if mag2(dr) < r2: hitlist.append([i,j])\n", " return hitlist\n", "\n", "nhisto = 0 # number of histogram snapshots to average\n", "\n", "while True:\n", " rate(150)\n", " # Accumulate and average histogram snapshots\n", " for i in range(len(accum)): accum[i][1] = (nhisto*accum[i][1] + histo[i])/(nhisto+1)\n", " if nhisto % 10 == 0:\n", " vdist.data = accum\n", " nhisto += 1\n", "\n", " # Update all positions\n", " for i in range(Natoms): Atoms[i].pos = apos[i] = apos[i] + (p[i]/mass)*dt\n", " \n", " # Check for collisions\n", " hitlist = checkCollisions(apos)\n", "\n", " # If any collisions took place, update momenta of the two atoms\n", " for ij in hitlist:\n", " i = ij[0]\n", " j = ij[1]\n", " ptot = p[i]+p[j]\n", " vi = p[i]/mass\n", " vj = p[j]/mass\n", " vrel = vj-vi\n", " a = vrel.mag2\n", " if a == 0: continue; # exactly same velocities\n", " rrel = apos[i]-apos[j]\n", " b = 2*rrel.dot(vrel)\n", " c = rrel.mag2-4*Ratom*Ratom\n", " d = b*b-4*a*c\n", " if d < 0: continue # something wrong; ignore this rare case\n", " deltat = (-b+sqrt(d))/(2*a) # t-deltat is when they made contact\n", " apos[i] = apos[i]-vi*deltat # back up to contact configuration\n", " apos[j] = apos[j]-vj*deltat\n", " mtot = 2*mass\n", " pcmi = p[i]-ptot*mass/mtot # transform momenta to cm frame\n", " pcmj = p[j]-ptot*mass/mtot\n", " rrel = norm(rrel)\n", " pcmi = pcmi-2*pcmi.dot(rrel)*rrel # bounce in cm frame\n", " pcmj = pcmj-2*pcmj.dot(rrel)*rrel\n", " p[i] = pcmi+ptot*mass/mtot # transform momenta back to lab frame\n", " p[j] = pcmj+ptot*mass/mtot\n", " apos[i] = apos[i]+(p[i]/mass)*deltat # move forward deltat in time\n", " apos[j] = apos[j]+(p[j]/mass)*deltat\n", " interchange(vi.mag, p[i].mag/mass)\n", " interchange(vj.mag, p[j].mag/mass)\n", " \n", " for i in range(Natoms):\n", " loc = apos[i]\n", " if abs(loc.x) > L/2:\n", " if loc.x < 0: p[i].x = abs(p[i].x)\n", " else: p[i].x = -abs(p[i].x)\n", " \n", " if abs(loc.y) > L/2:\n", " if loc.y < 0: p[i].y = abs(p[i].y)\n", " else: p[i].y = -abs(p[i].y)\n", " \n", " if abs(loc.z) > L/2:\n", " if loc.z < 0: p[i].z = abs(p[i].z)\n", " else: p[i].z = -abs(p[i].z)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "VPython", "language": "python", "name": "vpython" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }