# Solve the pendulum using numerical approximation
# Copyright (c) 2008 David M. Harrison
# The next line is an internal revision control id:
# $Date: 2008/05/10 10:33:45 $, $Revision: 1.1 $
# Import the visual library
from visual import *
# The initial angle in radians.
theta = pi/2.0
# The initial angular velocity
omega = 0
# Set g and the length of the pendulum
g = 9.80
L = 1.00
# These four lines control the size of the window of
# the animation and the scale. The details of these lines
# are not important for our purposes.
scene.autoscale = 0
scene.height = 600
scene.width = 600
scene.range = vector(2.0,2.0,2.0)
#
# Now we build the pendulum which we will animate.
#
# The support for the pendulum
support = cylinder( pos = (0, 0, -0.5), axis = (0,0,1), radius = 0.02)
# The "frame" construct groups two or more objects into a single one.
# Here we group the cylinder and the sphere into a single object
# which is the pendulum.
pendulum = frame()
cylinder(frame=pendulum, pos=(0,0,0), radius=0.01, length=1, color=color.cyan)
sphere(frame=pendulum, pos=(1,0,0), radius=0.1, color=color.red)
# Position the pendulum.
pendulum.pos = (0,0,0)
# Rotate the pendulum about the z axis. Note that VPython measures
# angles with respect to the x (horizontal) axis. We are measuring
# angles with respect to the vertical (-y axis) so we subtract
# pi/2.0 radians from the angle.
pendulum.rotate(axis = (0,0,1), angle = theta - pi/2.0)
# The time
t = 0.
# Below we will want to store the old value of the time.
# Set it the "impossible" value of -1 initially.
t_old = -1.
# The time step
dt = 0.0005
# The value "1" is equivalent to true. So this causes the while
# loop to run forever.
while 1:
# Set the rate of the animation in frames per second
rate(1/dt)
# The angular acceleration, i.e. the second derivative of the
# angle with respect to time.
alpha = -(g/L) * sin(theta)
# The new value of the angular velocity
omega = omega + alpha * dt
# The change in the angle of the pendulum
d_theta = omega * dt
# A rough and ready way to estimate the period of the oscillation.
# It the angle is positive and adding d_theta to it will make
# it negative, then it is going through the vertical
# from right to left.
if(theta > 0 and theta + d_theta < 0) :
# If t_old is > 0, then this is not the first cycle of
# the oscillation. The difference between t and t_old
# is the period within the resolution of the time step dt
# and rounding errors. Print the period.
if(t_old > 0):
print "Estimated Period =", t - t_old, "s"
# Store the current value of the time in t_old
t_old = t
# Rotate the pendulum about the z axis by the change in the angle
pendulum.rotate(axis = (0,0,1), angle = d_theta)
# Update the value of the angle
theta = theta + d_theta
# Update the time
t = t + dt
