Hockey Stick in Harmonic Motion

Hey everyone - paste this code into a Visual Python window to see a hockey stick oscillating in harmonic motion.


from __future__ import division # makes sure is decimal and not integer
from visual import *
from visual.graph import * # graphing capability
scene.title = 'Simple Harmonic Motion with "Hockey Stick"'
pendulum = frame()
foot = box(frame = pendulum, pos = (.0885,.009,0), length=.15, height=.018, width=.006, color = (.5, .5, .5), materials = materials.marble)
stick = box(frame = pendulum, pos = (.0135,.518,0), length=.027, height=1, width = .006, color=color.orange, material=materials.wood)
#distance from top of hole to top end of stick
dhole = .001
#starting position
startAngle = 10
theta = math.radians(startAngle)
pendulum.rotate(angle = theta, axis = (0,0,1), origin = (0.5*stick.length, stick.height + foot.height - dhole ,0))
angVel = 0
angAcc = 0
# masses
foot.mass = .178
stick.mass = .129
# moments
xcmLength = (foot.mass*(0.5*foot.length + 0.5*stick.length) + stick.mass*(0.5*stick.length))/(stick.mass + foot.mass)
ycmLength = (foot.mass*(0.5*foot.height) + stick.mass*(0.5*stick.height + foot.height)) / (foot.mass + stick.mass)
stick.I = ((1/12)*stick.mass*stick.height*stick.height + stick.mass*(0.5*stick.height - dhole)*(0.5*stick.height - dhole))
foot.I = ((1/3)*foot.mass*foot.length*foot.length + foot.mass*(stick.height - dhole + 0.5*foot.height)*(stick.height - dhole + 0.5*foot.height))
I = stick.I + foot.I
#set up angle graph
graph1 = gdisplay(x=150, y=600, width=400, height=300,
title='Angle vs. Time', xtitle='time (s)', ytitle='angle (rad)',
xmax=10., xmin=0., ymax=1.1*theta, ymin=-1.1*theta,
foreground=color.black, background=color.white)
acurve = gcurve(gdisplay = graph1, color = color.black)
#time
time = 0
dt = 0.01
while True:
rate(100)
time += dt
gravTorque = -1 * math.pow((xcmLength - 0.5*stick.length)*(xcmLength - 0.5*stick.length) + (stick.height + foot.height - dhole - ycmLength)*(stick.height + foot.height - dhole - ycmLength), 0.5) * (foot.mass + stick.mass) * 9.8 * sin(theta)
angAcc = gravTorque / I
angVel += angAcc * dt
theta += angVel * dt
pendulum.rotate(angle = angVel * dt, axis = (0,0,1), origin = (0.5*stick.length, stick.height + foot.height - dhole ,0))
acurve.plot(pos = (time, theta))