HardSphereGas.ipynb

{
 "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": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/vpython_libraries/glowcomm\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");})"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");})"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");})"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "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
}