diff --git a/site-packages/visual/examples/doublependulum.py b/site-packages/visual/examples/doublependulum.py
index d3c5d29..ef58c17 100644
--- a/site-packages/visual/examples/doublependulum.py
+++ b/site-packages/visual/examples/doublependulum.py
@@ -1,69 +1,76 @@
-from __future__ import print_function, division
-from visual import *
-#from visual.graph import *
+#!/usr/bin/env python
+#coding:utf-8
+"""Double pendulum
+
+The analysis is in terms of Lagrangian mechanics.
+The Lagrangian variables are angle of upper bar and angle of lower bar,
+measured from the vertical.
 
-# Double pendulum
+Bruce Sherwood
 
-# The analysis is in terms of Lagrangian mechanics.
-# The Lagrangian variables are angle of upper bar and angle of lower bar,
-# measured from the vertical.
+Corrections to the Lagrangian calculations by Owen Long, UC. Riverside
+
+"""
+from __future__ import print_function, division
+import visual as vp
 
-# Bruce Sherwood
-#
-# Corrections to the Lagrangian calculations by Owen Long, UC. Riverside
-#
+print(__doc__)
 
-scene.title = 'Double Pendulum'
-scene.height = scene.width = 800
+vp.scene.title = 'Double Pendulum'
+vp.scene.height = vp.scene.width = 800
 
-g = 9.8
+G = 9.8
 M1 = 2.0
 M2 = 1.0
 L1 = 0.5 # physical length of upper assembly; distance between axles
 L2 = 1.0 # physical length of lower bar
 I1 = M1*L1**2/12 # moment of inertia of upper assembly
 I2 = M2*L2**2/12 # moment of inertia of lower bar
-d = 0.05 # thickness of each bar
-gap = 2*d # distance between two parts of upper, U-shaped assembly
-L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle
-L2display = L2+d/2 # show lower bar a bit longer than physical, to overlap axle
-
-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
-top = vector(0,0,0) # top of inner bar of U-shaped upper assembly
-scene.center = top-vector(0,(L1+L2)/2.,0)
-
-theta1 = 2.1 # initial upper angle (from vertical)
-theta2 = 2.4 # initial lower angle (from vertical)
-theta1dot = 0 # initial rate of change of theta1
-theta2dot = 0 # initial rate of change of theta2
-
-pedestal = box(pos=top-vector(0,hpedestal/2.,offset),
-                 height=1.1*hpedestal, length=wpedestal, width=wpedestal,
-                 color=(0.4,0.4,0.5))
-base = box(pos=top-vector(0,hpedestal+tbase/2.,offset),
-                 height=tbase, length=wbase, width=wbase,
-                 color=pedestal.color)
-axle = cylinder(pos=top-vector(0,0,gap/2.-d/4.), axis=(0,0,-offset), radius=d/4., color=color.yellow)
-
-frame1 = frame(pos=top)
-bar1 = box(frame=frame1, pos=(L1display/2.-d/2.,0,-(gap+d)/2.), size=(L1display,d,d), color=color.red)
-bar1b = box(frame=frame1, pos=(L1display/2.-d/2.,0,(gap+d)/2.), size=(L1display,d,d), color=color.red)
-axle1 = cylinder(frame=frame1, pos=(L1,0,-(gap+d)/2.), axis=(0,0,gap+d),
-                 radius=axle.radius, color=axle.color)
-frame1.axis = (0,-1,0)
-frame2 = frame(pos=frame1.axis*L1)
-bar2 = box(frame=frame2, pos=(L2display/2.-d/2.,0,0), size=(L2display,d,d), color=color.green)
-frame2.axis = (0,-1,0)
-frame1.rotate(axis=(0,0,1), angle=theta1)
-frame2.rotate(axis=(0,0,1), angle=theta2)
-frame2.pos = top+frame1.axis*L1
-
-dt = 0.00005 # energy constancy check fails for dt > 0.0003
-t = 0.
+D = 0.05 # thickness of each bar
+GAP = 2*D # distance between two parts of upper, U-shaped assembly
+L1_DISPLAY = L1+D # show upper assembly a bit longer than physical, to overlap axle
+L2_DISPLAY = L2+D/2 # show lower bar a bit longer than physical, to overlap axle
+
+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
+TOP = vp.vector(0, 0, 0) # top of inner bar of U-shaped upper assembly
+vp.scene.center = TOP- vp.vector(0, (L1+L2)/2., 0)
+
+THETA1 = 2.1 # initial upper angle (from vertical)
+THETA2 = 2.4 # initial lower angle (from vertical)
+THETA1DOT = 0 # initial rate of change of theta1
+THETA2DOT = 0 # initial rate of change of theta2
+
+PEDISTAL = vp.box(pos=TOP-vp.vector(0, HPEDESTAL/2., OFFSET),
+                  height=1.1*HPEDESTAL, length=WPEDESTAL, width=WPEDESTAL,
+                  color=(0.4, 0.4, 0.5))
+BASE = vp.box(pos=TOP-vp.vector(0, HPEDESTAL+TBASE/2., OFFSET),
+              height=TBASE, length=WBASE, width=WBASE,
+              color=PEDISTAL.color)
+AXLE = vp.cylinder(pos=TOP-vp.vector(0, 0, GAP/2.-D/4.), axis=(0, 0, -OFFSET),
+                   radius=D/4., color=vp.color.yellow)
+
+FRAME1 = vp.frame(pos=TOP)
+BAR1 = vp.box(frame=FRAME1, pos=(L1_DISPLAY/2.-D/2., 0, -(GAP+D)/2.),
+              size=(L1_DISPLAY, D, D), color=vp.color.red)
+BAR1B = vp.box(frame=FRAME1, pos=(L1_DISPLAY/2.-D/2., 0, (GAP+D)/2.),
+               size=(L1_DISPLAY, D, D), color=vp.color.red)
+AXLE1 = vp.cylinder(frame=FRAME1, pos=(L1, 0, -(GAP+D)/2.), axis=(0, 0, GAP+D),
+                    radius=AXLE.radius, color=AXLE.color)
+FRAME1.axis = (0, -1, 0)
+FRAME2 = vp.frame(pos=FRAME1.axis*L1)
+BAR2 = vp.box(frame=FRAME2, pos=(L2_DISPLAY/2.-D/2., 0, 0),
+              size=(L2_DISPLAY, D, D), color=vp.color.green)
+FRAME2.axis = (0, -1, 0)
+FRAME1.rotate(axis=(0, 0, 1), angle=THETA1)
+FRAME2.rotate(axis=(0, 0, 1), angle=THETA2)
+FRAME2.pos = TOP+FRAME1.axis*L1
+
+DT = 0.00005 # energy constancy check fails for dt > 0.0003
+T = 0.
 
 C11 = (0.25*M1+M2)*L1**2+I1
 C22 = 0.25*M2*L2**2+I2
@@ -75,34 +82,34 @@
 ##gE = gcurve(color=color.red)
 
 while True:
-    rate(1/dt)
+    vp.rate(1/DT)
     # Calculate accelerations of the Lagrangian coordinates:
-    C12 = C21 = 0.5*M2*L1*L2*cos(theta1-theta2)
-    Cdet = C11*C22-C12*C21
-    a = .5*M2*L1*L2*sin(theta1-theta2)
-    A = -(.5*M1+M2)*g*L1*sin(theta1)-a*theta2dot**2
-    B = -.5*M2*g*L2*sin(theta2)+a*theta1dot**2
-    atheta1 = (C22*A-C12*B)/Cdet
-    atheta2 = (-C21*A+C11*B)/Cdet
+    C12 = C21 = 0.5*M2*L1*L2* vp.cos(THETA1-THETA2)
+    CDET = C11*C22-C12*C21
+    A0 = .5*M2*L1*L2* vp.sin(THETA1-THETA2)
+    A = -(.5*M1+M2)*G*L1* vp.sin(THETA1)-A0*THETA2DOT**2
+    B = -.5*M2*G*L2* vp.sin(THETA2)+A0*THETA1DOT**2
+    ATHETA1 = (C22*A-C12*B)/CDET
+    ATHETA2 = (-C21*A+C11*B)/CDET
     # Update velocities of the Lagrangian coordinates:
-    theta1dot += atheta1*dt
-    theta2dot += atheta2*dt
+    THETA1DOT += ATHETA1*DT
+    THETA2DOT += ATHETA2*DT
     # Update Lagrangian coordinates:
-    dtheta1 = theta1dot*dt
-    dtheta2 = theta2dot*dt
-    theta1 += dtheta1
-    theta2 += dtheta2
+    DTHETA1 = THETA1DOT*DT
+    DTHETA2 = THETA2DOT*DT
+    THETA1 += DTHETA1
+    THETA2 += DTHETA2
 
-    frame1.rotate(axis=(0,0,1), angle=dtheta1)
-    frame2.pos = top+frame1.axis*L1
-    frame2.rotate(axis=(0,0,1), angle=dtheta2)
-    t = t+dt
+    FRAME1.rotate(axis=(0, 0, 1), angle=DTHETA1)
+    FRAME2.pos = TOP+FRAME1.axis*L1
+    FRAME2.rotate(axis=(0, 0, 1), angle=DTHETA2)
+    T = T+DT
 
 ## Energy check:
 ##    K = .5*((.25*M1+M2)*L1**2+I1)*theta1dot**2+.5*(.25*M2*L2**2+I2)*theta2dot**2+\
 ##        .5*M2*L1*L2*cos(theta1-theta2)*theta1dot*theta2dot
 ##    U = -(.5*M1+M2)*g*L1*cos(theta1)-.5*M2*g*L2*cos(theta2)
-##    gK.plot(pos=(t,K))
-##    gU.plot(pos=(t,U))
-##    gE.plot(pos=(t,K+U))
+##    gK.plot(pos=(t, K))
+##    gU.plot(pos=(t, U))
+##    gE.plot(pos=(t, K+U))