{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "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.undef(\"nbextensions/vpython_libraries/glow.2.2.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([\"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": {
      "application/javascript": [
       "require([\"nbextensions/vpython_libraries/glow.2.2.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": {
      "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",
    "# Gyroscope sitting on a pedestal\n",
    "\n",
    "# The analysis is in terms of Lagrangian mechanics.\n",
    "# The Lagrangian variables are polar angle theta,\n",
    "# azimuthal angle phi, and spin angle psi.\n",
    "\n",
    "scene.width = 800\n",
    "scene.height = 600\n",
    "scene.range = 1.2\n",
    "scene.title = \"A precessing, nutating gyroscope\"\n",
    "\n",
    "Lshaft = 1 # length of gyroscope shaft\n",
    "r = Lshaft/2 # distance from support to center of mass\n",
    "Rshaft = 0.03 # radius of gyroscope shaft\n",
    "M = 1 # mass of gyroscope (massless shaft)\n",
    "Rrotor = 0.4 # radius of gyroscope rotor\n",
    "Drotor = 0.1 # thickness of gyroscope rotor\n",
    "I3 = 0.5*M*Rrotor**2 # moment of inertia of gyroscope about its own axis\n",
    "I1 = M*r**2 + .5*I3 # moment of inertia about a line through the support, perpendicular to the axis\n",
    "hpedestal = Lshaft # height of pedestal\n",
    "wpedestal = 0.1 # width of pedestal\n",
    "tbase = 0.05 # thickness of base\n",
    "wbase = 3*wpedestal # width of base\n",
    "g = 9.8\n",
    "Fgrav = vector(0,-M*g,0)\n",
    "top = vector(0,0,0) # top of pedestal\n",
    "\n",
    "shaft = cylinder(length=Lshaft, radius=Rshaft, color=color.orange)\n",
    "rotor = cylinder(pos=vector(Lshaft/2-Drotor/2,0,0), axis=vector(Drotor, 0, 0),\n",
    "                 radius=Rrotor, color=color.gray(0.9))\n",
    "base = sphere(color=shaft.color, radius=Rshaft)\n",
    "end = sphere(pos=vector(Lshaft,0,0), color=shaft.color, radius=Rshaft)\n",
    "gyro = compound([shaft, rotor, base, end])\n",
    "gyro_center = gyro.pos\n",
    "gyro.texture = textures.metal\n",
    "tip = sphere(pos=shaft.axis, radius=shaft.radius/2,  make_trail=True, retain=250)\n",
    "tip.trail_color = color.green\n",
    "tip.trail_radius = 0.15*Rshaft\n",
    "\n",
    "pedestal = box(pos=top-vector(0,hpedestal/2+shaft.radius/2,0), \n",
    "               height=hpedestal-shaft.radius, length=wpedestal, width=wpedestal, texture=textures.wood)\n",
    "pedestal_base = box(pos=top-vector(0,hpedestal+tbase/2,0), \n",
    "                    height=tbase, length=wbase, width=wbase, texture=textures.wood)\n",
    "\n",
    "theta = 0\n",
    "thetadot = 0\n",
    "psi = 0\n",
    "psidot = 0\n",
    "phi = 0\n",
    "phidot = 0\n",
    "pureprecession = False\n",
    "\n",
    "def reset():\n",
    "    global theta, thetadot, psi, psidot, phi, phidot\n",
    "    theta = 0.3*pi # initial polar angle of shaft (from vertical)\n",
    "    thetadot = 0 # initial rate of change of polar angle\n",
    "    psi = 0 # initial spin angle\n",
    "    psidot = 30 # initial rate of change of spin angle (spin ang. velocity)\n",
    "    phi = -pi/2 # initial azimuthal angle\n",
    "    phidot = 0 # initial rate of change of azimuthal angle\n",
    "    if pureprecession: # Set to True if you want pure precession, without nutation\n",
    "        a = (1-I3/I1)*sin(theta)*cos(theta)\n",
    "        b = -(I3/I1)*psidot*sin(theta)\n",
    "        c = M*g*r*sin(theta)/I1\n",
    "        phidot = (-b+sqrt(b**2-4*a*c))/(2*a)\n",
    "    gyro.axis = gyro.length*vector(sin(theta)*sin(phi),cos(theta),sin(theta)*cos(phi))\n",
    "    A = norm(gyro.axis)\n",
    "    gyro.pos = 0.5*Lshaft*A\n",
    "    tip.pos = Lshaft*A\n",
    "    tip.clear_trail()\n",
    "\n",
    "reset()\n",
    "#scene.waitfor('textures') # not yet implemented\n",
    "\n",
    "dt = 0.0001\n",
    "t = 0\n",
    "Nsteps = 20 # number of calculational steps between graphics updates\n",
    "\n",
    "while True:\n",
    "    rate(200)\n",
    "    for step in range(Nsteps): # multiple calculation steps for accuracy\n",
    "        # Calculate accelerations of the Lagrangian coordinates:\n",
    "        atheta = sin(theta)*cos(theta)*phidot**2+( M*g*r*sin(theta)-I3*(psidot+phidot*cos(theta))*phidot*sin(theta))/I1\n",
    "        aphi = (I3/I1)*(psidot+phidot*cos(theta))*thetadot/sin(theta)-2*cos(theta)*thetadot*phidot/sin(theta)\n",
    "        apsi = phidot*thetadot*sin(theta)-aphi*cos(theta)\n",
    "        # Update velocities of the Lagrangian coordinates:\n",
    "        thetadot += atheta*dt\n",
    "        phidot += aphi*dt\n",
    "        psidot += apsi*dt\n",
    "        # Update Lagrangian coordinates:\n",
    "        theta += thetadot*dt\n",
    "        phi += phidot*dt\n",
    "        psi += psidot*dt\n",
    "\n",
    "    gyro.axis = gyro.length*vector(sin(theta)*sin(phi),cos(theta),sin(theta)*cos(phi))\n",
    "    # Display approximate rotation of rotor and shaft:\n",
    "    gyro.rotate(angle=psidot*dt*Nsteps)\n",
    "    A = norm(gyro.axis)\n",
    "    gyro.pos = 0.5*Lshaft*A\n",
    "    tip.pos = Lshaft*A\n",
    "    t = t+dt*Nsteps\n",
    "    \n",
    "\n"
   ]
  },
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}