Double_pendulum.ipynb

Kernel: Python 2 (Ubuntu Linux)

from __future__ import division
from vpython import *
from math import pi, sqrt,sin,cos
# Double pendulum

# The analysis is in terms of Lagrangian mechanics.
# The Lagrangian variables are angle of upper bar, angle of lower bar,
# measured from the vertical.

# Bruce Sherwood

# Corrections to the Lagrangian calculations by Owen Long, UC. Riverside

scene.width = scene.height = 600
scene.range = 1.8
scene.title = "A double pendulum"

def display_instructions():
    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"
    s += "\nTouch screen= pinch/extend to zoom, swipe or two-finger rotate."
    scene.caption = s

# Display text below the 3D graphics:
display_instructions()

g = 9.8
M1 = 2.0
M2 = 1.0
d = 0.05 # thickness of each bar
gap = 2*d # distance between two parts of upper, U-shaped assembly
L1 = 0.5 # physical length of upper assembly; distance between axles
L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle
L2 = 1 # physical length of lower bar
L2display = L2+d/2 # show lower bar a bit longer than physical, to overlap axle
# Coefficients used in Lagrangian calculation
A = (1/4)*M1*L1**2+(1/12)*M1*L1**2+M2*L1**2
B = (1/2)*M2*L1*L2
C = g*L1*(M1/2+M2)
D = M2*L1*L2/2
E = (1/12)*M2*L2**2+(1/4)*M2*L2**2
F = g*L2*M2/2

hpedestal = 1.3*(L1+L2) # height of pedestal
wpedestal = 0.1 # width of pedestal
tbase = 0.05 # thickness of base
wbase = 8*gap # width of base
offset = 2*gap # from center of pedestal to center of U-shaped upper assembly
pedestal_top = vec(0,hpedestal/2,0) # top of inner bar of U-shaped upper assembly

theta1 = 1.3*pi/2 # initial upper angle (from vertical)
theta1dot = 0 # initial rate of change of theta1
theta2 = 0 # initial lower angle (from vertical)
theta2dot = 0 # initial rate of change of theta2

pedestal = box( pos=pedestal_top-vec(0,hpedestal/2,offset),
                size=vec(wpedestal,1.1*hpedestal,wpedestal),
                color=vec(0.4,0.4,0.5) )
base = box( pos=pedestal_top-vec(0,hpedestal+tbase/2,offset),
                 size=vec(wbase,tbase,wbase),
                 color=pedestal.color )
axle1 = cylinder( pos=pedestal_top-vec(0,0,gap/2-d/4), axis=vec(0,0,-1),
                 size=vec(offset,d/4,d/4), color=color.yellow )

bar1 = box( pos=pedestal_top+vec(L1display/2-d/2,0,-(gap+d)/2), 
                 size=vec(L1display,d,d), color=color.red )
bar1.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1.pos.z) )
bar1.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1.pos.z) )

bar1b = box( pos=pedestal_top+vec(L1display/2-d/2,0,(gap+d)/2), 
                 size=vec(L1display,d,d), color=bar1.color )
bar1b.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1b.pos.z) )
bar1b.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, bar1b.pos.z) )

pivot1 = vec(axle1.pos.x, axle1.pos.y, 0)

axle2 = cylinder( pos=pedestal_top+vec(L1,0,-(gap+d)/2), axis=vec(0,0,1), 
                size=vec(gap+d,axle1.size.y/2,axle1.size.y/2), color=axle1.color )
axle2.rotate( angle=-pi/2, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, axle2.pos.z) )
axle2.rotate( angle=theta1, axis=vec(0,0,1), origin=vec(axle1.pos.x, axle1.pos.y, axle2.pos.z) )

bar2 = box( pos=axle2.pos+vec(L2display/2-d/2,0,(gap+d)/2), 
        size=vec(L2display,d,d), color=color.green )

bar2.rotate( angle=-pi/2,  axis=vec(0,0,1), origin=vec(axle2.pos.x, axle2.pos.y, bar2.pos.z) )
bar2.rotate( angle=theta2,  axis=vec(0,0,1), origin=vec(axle2.pos.x, axle2.pos.y, bar2.pos.z) )

dt = 0.001
t = 0

while True:
    rate(1/dt) 
    # Calculate accelerations of the Lagrangian coordinates=
    atheta1 = ((E*C/B)*sin(theta1)-F*sin(theta2))/(D-E*A/B)
    atheta2 = -(A*atheta1+C*sin(theta1))/B
    # Update velocities of the Lagrangian coordinates=
    theta1dot = theta1dot+atheta1*dt
    theta2dot = theta2dot+atheta2*dt
    # Update Lagrangian coordinates=
    dtheta1 = theta1dot*dt
    dtheta2 = theta2dot*dt
    theta1 = theta1+dtheta1
    theta2 = theta2+dtheta2
    
    bar1.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )
    bar1b.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )
    pivot2 = vec(axle2.pos.x, axle2.pos.y, pivot1.z)
    axle2.rotate( angle=dtheta1, axis=vec(0,0,1), origin=pivot1 )
    bar2.rotate( angle=dtheta2, axis=vec(0,0,1), origin=pivot2 )
    pivot2 = vec(axle2.pos.x, axle2.pos.y, pivot1.z)
    bar2.pos = pivot2 + bar2.axis/2
    
    t = t+dt

MIME type unknown not supported

MIME type unknown not supported

MIME type unknown not supported

MIME type unknown not supported

MIME type unknown not supported

MIME type unknown not supported

MIME type unknown not supported

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-1-b51e67e16045> in <module>()
     91 
     92 while True:
---> 93     rate(1/dt)
     94     # Calculate accelerations of the Lagrangian coordinates=
     95     atheta1 = ((E*C/B)*sin(theta1)-F*sin(theta2))/(D-E*A/B)
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/vpython/vpython.pyc in __call__(self, N)
    209         self.rval = N
    210         if self.rval < 1: raise ValueError("rate value must be greater than or equal to 1")
--> 211         super(_RateKeeper2, self).__call__(self.rval) ## calls __call__ in rate_control.py
    212 
    213 if sys.version > '3':
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/vpython/rate_control.pyc in __call__(self, maxRate)
    205             dt = self.lastSleep + self.calls*(self.userTime + self.callTime + self.delay) + \
    206                  renders*self.renderTime + sleeps*self.interactionPeriod - _clock()
--> 207             _sleep(dt)
    208             self.lastSleep = _clock()
    209             self.calls = 0
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/vpython/rate_control.pyc in _sleep(dt)
     47         dtsleep = nticks*_tick
     48         t = _clock()
---> 49         time.sleep(dtsleep)
     50         t = _clock()-t
     51         dt -= t
KeyboardInterrupt: