{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [ { "data": { "application/javascript": [ "require.undef(\"nbextensions/jquery-ui.custom.min\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require.undef(\"nbextensions/glow.2.1a.min\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require.undef(\"nbextensions/glow.2.1b.min\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require.undef(\"nbextensions/glowcomm\");" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": [ "require([\"nbextensions/glowcomm\"], function(){console.log(\"glowcomm 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" } ], "source": [ "from vpython import *\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 VPython programs:\\n\"\n", " s += \" Rotate the camera by dragging with the right mouse button,\\n or hold down the Ctrl key and drag.\\n\"\n", " s += \" To zoom, drag with the left+right mouse buttons,\\n or hold down the Alt/Option key and drag,\\n or use the mouse wheel.\\n\"\n", " s += \"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.001\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\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 }