from visual import *

# Gyroscope sitting on a pedestal

# The analysis is in terms of Lagrangian mechanics.
# The Lagrangian variables are polar angle theta,
# azimuthal angle phi, and spin angle alpha.

scene.width=800
scene.height=800
scene.title='Nutating Gyroscope'

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
Dsquare = 1.4*Drotor # thickness of square that turns with rotor
I = 0.5*M*Rrotor**2. # moment of inertia of gyroscope
hpedestal = Lshaft # height of pedestal
wpedestal = 0.1 # width of pedestal
tbase = 0.05 # thickness of base
wbase = 3.*wpedestal # width of base
g = 9.8
Fgrav = vector(0,-M*g,0)
top = vector(0,0,0) # top of pedestal

theta = pi/3. # initial polar angle of shaft (from vertical)
thetadot = 0 # initial rate of change of polar angle
alpha = 0 # initial spin angle
alphadot = 15 # initial rate of change of spin angle (spin ang. velocity)
phi = -pi/2. # initial azimuthal angle
phidot = 0 # initial rate of change of azimuthal angle
# Comment in following line to get pure precession
##phidot = (-alphadot+sqrt(alphadot**2+2*M*g*r*cos(theta)/I))/cos(theta)

pedestal = box(pos=top-vector(0,hpedestal/2.,0),
                 height=hpedestal, length=wpedestal, width=wpedestal,
                 color=(0.4,0.4,0.5))
base = box(pos=top-vector(0,hpedestal+tbase/2.,0),
                 height=tbase, length=wbase, width=wbase,
                 color=pedestal.color)

gyro=frame(axis=(sin(theta)*sin(phi),cos(theta),sin(theta)*cos(phi)))
shaft = cylinder(axis=(Lshaft,0,0),
                 radius=Rshaft, color=(0,1,0), frame=gyro)
rotor = cylinder(pos=(Lshaft/2 - Drotor/2, 0, 0),
                 axis=(Drotor, 0, 0),
                 radius=Rrotor, color=(1,0,0), frame=gyro)

trail = curve(radius=Rshaft/8., color=(1,1,0))

scene.autoscale = 0

r = Lshaft/2.
dt = 0.0001
t = 0.
Nsteps = 20 # number of calculational steps between graphics updates

while 1:
    rate(100)
    for step in range(Nsteps): # multiple calculation steps for accuracy
        # Calculate accelerations of the Lagrangian coordinates:
        atheta = (phidot**2*sin(theta)*cos(theta)
                  -2.*(alphadot+phidot*cos(theta))*phidot*sin(theta)
                  +2.*M*g*r*sin(theta)/I)
        aphi = 2.*thetadot*(alphadot-phidot*cos(theta))/sin(theta)
        aalpha = phidot*thetadot*sin(theta)-aphi*cos(theta)
        # Update velocities of the Lagrangian coordinates:
        thetadot = thetadot+atheta*dt
        phidot = phidot+aphi*dt
        alphadot = alphadot+aalpha*dt
        # Update Lagrangian coordinates:
        theta = theta+thetadot*dt
        phi = phi+phidot*dt
        alpha = alpha+alphadot*dt

    gyro.axis = vector(sin(theta)*sin(phi),cos(theta),sin(theta)*cos(phi))
    # Display approximate rotation of rotor and shaft:
    gyro.rotate(angle=alphadot*dt*Nsteps)
    trail.append(pos=gyro.pos + gyro.axis * Lshaft)
    t = t+dt*Nsteps