

In [ ]:

    
from vpython import *

# Gyroscope hanging from a spring
# Bruce Sherwood

scene.title = 'Gyroscope'
scene.background = color.gray(0.9)

top = vector(0,1,0) # where top of spring is held
ks = 100       # spring stiffness
Lspring = 1    # unstretched length of spring
Rspring = 0.03 # radius of spring
Dspring = 0.03 # thickness of spring wire
Lshaft = 1     # length of gyroscope shaft
Rshaft = 0.03  # radius of gyroscope shaft
M = 1          # mass of gyroscope (massless shaft)
Rrotor = 0.4   # radius of gyroscope rotor
Drotor = 0.1   # thickness of gyroscope rotor
omega = 40     # angular speed of rotor around axis
Dsquare = 1.4*Drotor # thickness of square that turns with rotor
I = 0.5*M*Rrotor**2  # moment of inertia of gyroscope
g = 9.8
Fgrav = vector(0,-M*g,0)
precession = M*g*(Lshaft/2)/(I*abs(omega))  # exact precession angular velocity
phi = atan(precession**2*(Lshaft/2)/g)      # approximate angle of spring to vertical
s = M*g/(ks*cos(phi))                       # approximate stretch of spring
# Refine estimate of angle of spring to vertical:
phi = 1/( ((I*abs(omega))/(M*Lshaft/2))**2/(g*Lshaft/2)-(Lspring+s)/(Lshaft/2) )
# Refine again:
s = M*g/(ks*cos(phi))
phi = 1/( ((I*abs(omega))/(M*Lshaft/2))**2/(g*Lshaft/2)-(Lspring+s)/(Lshaft/2) )
pprecess = vector(0,-1,M*precession*(Lshaft/2+(Lspring+s)*sin(phi)))
# Momentum required for completely smooth precession:
##pprecess = vector(0,0,M*precession*(Lshaft/2+(Lspring+s)*sin(phi)))
if omega < 0:
    pprecess = -pprecess

support = box(pos=top+vector(0,0.01,0), size=vector(0.2,0.02,0.2), color=color.green)
spring = helix(pos=top, axis=vector(-(Lspring+s)*sin(phi),
                                    -(Lspring+s)*cos(phi),0), coils=10, radius=Rspring, 
                                     thickness=Dspring, color=vector(1,0.7,0.2))
gyropos = top+spring.axis 
gyroaxis = vector(1,0,0)
shaft = cylinder(pos=vector(0,0,0), axis=Lshaft*gyroaxis, 
                 radius=Rshaft, color=color.gray(0.8))
rotor = cylinder(pos=0.5*gyroaxis*(Lshaft-Drotor),
                 axis=gyroaxis*Drotor, radius=Rrotor, color=color.gray(0.8))

stripe1 = box(color=color.black, size=vector(.01,2*Rrotor,.01),
                pos=rotor.pos+1.02*rotor.axis)
stripe2 = box(color=color.black, size=vector(.01,2*Rrotor,.01),
                pos=rotor.pos-0.02*rotor.axis)
gyro = compound([shaft,rotor,stripe1,stripe2], 
                pos=gyropos, axis=gyroaxis)
gyro.texture = textures.metal
gyro.rotate(axis=vector(0,1,0), angle=pi)

# center of mass of shaft
cm = gyro.pos+0.5*Lshaft*gyro.axis 
Lrot = I*omega*gyro.axis
p = pprecess
dt = 0.01
t = 0

while True:
    rate(50) 
    Fspring = -ks*norm(spring.axis)*(mag(spring.axis)-Lspring)
    # Calculate torque about center of mass:
    torque = cross(-0.5*Lshaft*gyro.axis,Fspring)
    Lrot = Lrot+torque*dt
    p = p+(Fgrav+Fspring)*dt
    cm = cm+(p/M)*dt

    # Update positions of shaft, rotor, spring, stripes
    if omega > 0:
        gyro.axis = norm(Lrot)
    else:
        gyro.axis = -norm(Lrot)
    gyro.pos = cm-0.5*Lshaft*gyro.axis # shaft rotated, adjust connection to spring
    spring.axis = gyro.pos - top
    gyro.rotate(angle=omega*dt/4, axis=gyro.axis) # spin easier to see if slower than actual omega