Double Pendulum SImulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glow.min\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/glowcomm\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require.undef(\"nbextensions/vpython_libraries/jquery-ui.custom.min\");}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glow.min\"], function(){console.log(\"GLOW LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/glowcomm\"], function(){console.log(\"GLOWCOMM LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") {require([\"nbextensions/vpython_libraries/jquery-ui.custom.min\"], function(){console.log(\"JQUERY LOADED\");});}else{element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from vpython import *\n",
    "\n",
    "# GlowScript 2.9 VPython\n",
    "# Double pendulum\n",
    "\n",
    "# The analysis is in terms of Lagrangian mechanics.\n",
    "# The Lagrangian variables are angle of upper bar, angle of lower bar,\n",
    "# measured from the vertical.\n",
    "\n",
    "# Bruce Sherwood\n",
    "\n",
    "# Corrections to the Lagrangian calculations by Owen Long, UC. Riverside\n",
    "\n",
    "scene.width = scene.height = 600\n",
    "scene.range = 1.8\n",
    "scene.title = \"A double pendulum\"\n",
    "\n",
    "def display_instructions():\n",
    "    s = \"\"\"In GlowScript programs:\n",
    "To rotate \"camera\", drag with right button or Ctrl-drag.\n",
    "To zoom, drag with middle button or Alt/Option depressed, or use scroll wheel.\n",
    "  On a two-button mouse, middle is left + right.\n",
    "To pan left/right and up/down, Shift-drag.\n",
    "Touch screen: pinch/extend to zoom, swipe or two-finger rotate.\"\"\"\n",
    "    scene.caption = s\n",
    "\n",
    "# Display text below the 3D graphics:\n",
    "display_instructions()\n",
    "\n",
    "g = 9.8\n",
    "M1 = 2.0\n",
    "M2 = 1.0\n",
    "d = 0.05 # thickness of each bar\n",
    "gap = 2*d # distance between two parts of upper, U-shaped assembly\n",
    "L1 = 0.5 # physical length of upper assembly; distance between axles\n",
    "L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle\n",
    "L2 = 1 # physical length of lower bar\n",
    "L2display = L2+d/2 # show lower bar a bit longer than physical, to overlap axle\n",
    "# Coefficients used in Lagrangian calculation\n",
    "A = (1/4)*M1*L1**2+(1/12)*M1*L1**2+M2*L1**2\n",
    "B = (1/2)*M2*L1*L2\n",
    "C = g*L1*(M1/2+M2)\n",
    "D = M2*L1*L2/2\n",
    "E = (1/12)*M2*L2**2+(1/4)*M2*L2**2\n",
    "F = g*L2*M2/2\n",
    "\n",
    "hpedestal = 1.3*(L1+L2) # height of pedestal\n",
    "wpedestal = 0.1 # width of pedestal\n",
    "tbase = 0.05 # thickness of base\n",
    "wbase = 8*gap # width of base\n",
    "offset = 2*gap # from center of pedestal to center of U-shaped upper assembly\n",
    "pedestal_top = vec(0,hpedestal/2,0) # top of inner bar of U-shaped upper assembly\n",
    "\n",
    "theta1 = 1.3*pi/2 # initial upper angle (from vertical)\n",
    "theta1dot = 0 # initial rate of change of theta1\n",
    "theta2 = 0 # initial lower angle (from vertical)\n",
    "theta2dot = 0 # initial rate of change of theta2\n",
    "\n",
    "pedestal = box( pos=pedestal_top-vec(0,hpedestal/2,offset),\n",
    "                size=vec(wpedestal,1.1*hpedestal,wpedestal),\n",
    "                color=vec(0.4,0.4,0.5) )\n",
    "base = box( pos=pedestal_top-vec(0,hpedestal+tbase/2,offset),\n",
    "                 size=vec(wbase,tbase,wbase),\n",
    "                 color=pedestal.color )\n",
    "axle1 = cylinder( pos=pedestal_top-vec(0,0,gap/2-d/4), axis=vec(0,0,-1),\n",
    "                 size=vec(offset,d/4,d/4), color=color.yellow )\n",
    "\n",
    "bar1 = box( pos=pedestal_top+vec(L1display/2-d/2,0,-(gap+d)/2), \n",
    "                 size=vec(L1display,d,d), color=color.red )\n",
    "bar1.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1.pos.z) )\n",
    "bar1.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1.pos.z) )\n",
    "\n",
    "bar1b = box( pos=pedestal_top+vec(L1display/2-d/2,0,(gap+d)/2), \n",
    "                 size=vec(L1display,d,d), color=bar1.color )\n",
    "bar1b.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1b.pos.z) )\n",
    "bar1b.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1b.pos.z) )\n",
    "\n",
    "pivot1 = vec(axle1.pos.x, axle1.pos.y, 0)\n",
    "\n",
    "axle2 = cylinder( pos=pedestal_top+vec(L1,0,-(gap+d)/2), axis=vec(0,0,1), \n",
    "                size=vec(gap+d,axle1.size.y/2,axle1.size.y/2), color=axle1.color )\n",
    "axle2.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, axle2.pos.z) )\n",
    "axle2.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, axle2.pos.z) )\n",
    "\n",
    "bar2 = box( pos=axle2.pos+vec(L2display/2-d/2,0,(gap+d)/2), \n",
    "        size=vec(L2display,d,d), color=color.green )\n",
    "\n",
    "bar2.rotate( angle=-pi/2,  axis=vec(0,0,1), origin=vec(axle2.pos.x, axle2.pos.y, bar2.pos.z) )\n",
    "bar2.rotate( angle=theta2,  axis=vec(0,0,1), origin=vec(axle2.pos.x, axle2.pos.y, bar2.pos.z) )\n",
    "\n",
    "dt = 0.0002\n",
    "t = 0\n",
    "\n",
    "while True:\n",
    "    rate(1/dt) \n",
    "    # Calculate accelerations of the Lagrangian coordinates=\n",
    "    atheta1 = ((E*C/B)*sin(theta1)-F*sin(theta2))/(D-E*A/B)\n",
    "    atheta2 = -(A*atheta1+C*sin(theta1))/B\n",
    "    # Update velocities of the Lagrangian coordinates=\n",
    "    theta1dot = theta1dot+atheta1*dt\n",
    "    theta2dot = theta2dot+atheta2*dt\n",
    "    # Update Lagrangian coordinates=\n",
    "    dtheta1 = theta1dot*dt\n",
    "    dtheta2 = theta2dot*dt\n",
    "    theta1 = theta1+dtheta1\n",
    "    theta2 = theta2+dtheta2\n",
    "    \n",
    "    bar1.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )\n",
    "    bar1b.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )\n",
    "    pivot2 = vec(axle2.pos.x, axle2.pos.y, pivot1.z)\n",
    "    axle2.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )\n",
    "    bar2.rotate( angle=dtheta2, axis=vec(0,0,1), origin=pivot2 )\n",
    "    pivot2 = vec(axle2.pos.x, axle2.pos.y, pivot1.z)\n",
    "    bar2.pos = pivot2 + bar2.axis/2\n",
    "    \n",
    "    t = t+dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}