

In [1]:

    
from __future__ import division, print_function
from ivisual import *
from math import *



In [11]:

    
#rotation options
rotation = 0 #barbell orientation doesn't change (mounted on frictionless axle, no torque acts on it)
#rotation = 1 #barbell rotates at same rate as rod
#rotation = 2 #barbell rotates a lot

#direction of rotation of barbell about its CM
direction = 1 #counterclockwise
#direction = -1 #clockwise

#create scene
scene=canvas(title="Barbell Angular Momentum")

#set background color
#scene.background = color.white

#mass of particles on barbell
rmass = 0.2

#length of the barbell
L = 2

#rod from origin to CM of barbell
rod=cylinder(pos=(0,0,0), axis=(4,0,0), radius=0.03, color=(1,.9,0))

#create a frame
barbell = frame()

#particle 1
s1=sphere(frame = barbell, pos=(0,L/2.,0), radius=rmass, color=(1,0,0))
s1.mass=0.01

#particle 2
s2=sphere(frame = barbell, pos=(0,-L/2.,0), radius=rmass, color=(0,0,1))
s2.mass=0.01

#barbell rod
rd=cylinder(frame = barbell, pos=s1.pos, axis=(s2.pos-s1.pos),
            color=(1,1,1), radius=0.04)

#barbell moment of inertia about CM
barbell.Icm = 2*s1.mass*(L/2)**2
barbell.pos = rod.pos+rod.axis

#moment of inertia of CM of barbell about origin
barbell.Iorig = 2*s1.mass*(mag(rod.axis)**2)

#angular velocity of the barbell about its CM
omegaCM = vector(0,0,pi)

#angular velocity of the CM of the barbell about the origin
omega = vector(0,0,pi/5)

t = 0.0
dt = 0.01
scene.range=5
Lscale = 2.0/0.1

LT=arrow(pos=rod.pos, axis=(0,0,0), color=color.cyan,
         shaftwidth = 0.2)
LR=arrow(pos=barbell.pos, axis=(0,0,0), color=color.green,
         shaftwidth = 0.2)

theta=0
phi=0
dphi=mag(omega)*dt
dtheta=mag(omegaCM)*dt

T=2*pi/mag(omega)

while t<2*T:
    rate(150)
    theta = theta + mag(omegaCM)*dt
    phi = phi + mag(omega)*dt
    rod.axis=4*vector(cos(phi), sin(phi),0)
    barbell.pos = rod.pos+rod.axis
    Ltrans = barbell.Iorig*omega
    Lrot = vector(0,0,0)
    if rotation == 1:
        barbell.rotate(angle=dphi, axis=omegaCM, origin=(barbell.pos))
        Lrot = barbell.Icm*omega 
    if rotation == 2:
        barbell.rotate(angle=direction*dtheta, axis=omegaCM, origin=(barbell.pos))
        Lrot = direction*barbell.Icm*omegaCM
    LT.axis = Ltrans*Lscale
    LR.pos = barbell.pos
    LR.axis = Lrot*Lscale
    t = t+dt