1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
from visual import *
from visual.text import *
# Gyroscope hanging from a spring
# Bruce Sherwood
print """
Click to pause, click to proceed.
Click on rotor to see angular momentum and impulse arrows.
Click again on rotor to hide these arrows.
"""
arrowsvisible = 0 # angular momentum and angular impulse arrows
scene.title = 'Gyroscope'
scene.visible = 0
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
Dsquare = 1.4*Drotor # thickness of square that turns with rotor
I = 0.5*M*Rrotor**2. # moment of inertia of gyroscope
omega = 40 # angular velocity of rotor along axis
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=(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=(1,0.7,0.2))
gyro = frame(pos=top+spring.axis) # gyro.pos at end of spring
gyro.axis = vector(1,0,0)
shaft = cylinder(pos=gyro.pos, axis=Lshaft*gyro.axis,
radius=Rshaft, color=(0.85,0.85,0.85))
rotor = cylinder(pos=0.5*gyro.axis*(Lshaft-Drotor),
axis=gyro.axis*Drotor,
radius=Rrotor, color=(0.5,0.5,0.5))
stripe1 = curve(frame=gyro, color=color.black,
pos=[rotor.pos+1.03*rotor.axis+vector(0,Rrotor,0),
rotor.pos+1.03*rotor.axis-vector(0,Rrotor,0)])
stripe1 = curve(frame=gyro, color=color.black,
pos=[rotor.pos-0.03*rotor.axis+vector(0,Rrotor,0),
rotor.pos-0.03*rotor.axis-vector(0,Rrotor,0)])
gyro.rotate(axis=(0,1,0), angle=pi)
cm = gyro.pos+0.5*Lshaft*gyro.axis # center of mass of shaft
Lrot = I*omega*gyro.axis
p = pprecess
dt = 0.01
t = 0.
Lrotarrow = arrow(length=0, shaftwidth=Rshaft,
color=color.red, visible=arrowsvisible)
Lrotscale = 0.2
rotimpulsearrow = arrow(length=0, shaftwidth=Lrotarrow.shaftwidth,
color=color.cyan, visible=arrowsvisible)
rotimpulsescale = 5.
Lrotlabel = text(string='L', height=0.06, depth=0.25, visible=arrowsvisible,
justify='center', color=Lrotarrow.color)
Lrotimpulselabel = text(string='DL', justify='center', visible=arrowsvisible,
height=0.06, depth=0.25, color=rotimpulsearrow.color)
while 1:
rate(50)
if scene.mouse.clicked: # pause the animation with a mouse click
m = scene.mouse.getclick()
if m.pick is rotor:
arrowsvisible = not arrowsvisible
Lrotarrow.visible = arrowsvisible
rotimpulsearrow.visible = arrowsvisible
for obj in Lrotlabel.objects:
obj.visible = arrowsvisible
for obj in Lrotimpulselabel.objects:
obj.visible = arrowsvisible
else:
while 1:
if scene.mouse.clicked:
scene.mouse.getclick()
break
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.) # spin easier to see if slower than actual omega
shaft.pos = gyro.pos
shaft.axis = Lshaft*gyro.axis
rotor.pos = gyro.pos+0.5*gyro.axis*(Lshaft-Drotor)
rotor.axis = gyro.axis*Drotor
# Update arrows representing angular momentum and angular impulse
Lrotarrow.pos = gyro.pos+0.5*Lshaft*gyro.axis+vector(0,2.*Rshaft,0)
Lrotarrow.axis = Lrot*Lrotscale
rotimpulsearrow.pos = Lrotarrow.pos+Lrotarrow.axis
rotimpulsearrow.axis = torque*dt*rotimpulsescale
Lrotlabel.frame.pos = Lrotarrow.pos+Lrotarrow.axis/2.+vector(0,Rshaft,0)
Lrotimpulselabel.frame.pos = rotimpulsearrow.pos+rotimpulsearrow.axis-vector(0,3.*Rshaft,0)
if t == 0.: # make sure everything is set up before first visible display
scene.visible = 1
scene.autoscale = 0
t = t+dt
|
