{VERSION 1 0 "X11/Motif" "1.0"}{GLOBALS 3 3}{FONT 0 "-adobe-helve
tica-bold-r-normal--14-*" "helvetica" "Helvetica-Bold" 8 12 0 "He
lvetica-Bold" 12}{FONT 1 "-adobe-times-medium-r-normal--14-*" "ti
mes" "Times-Roman" 4 12 64 "Times-Roman" 12}{FONT 2 "-adobe-couri
er-medium-r-normal--14-*" "courier" "Courier" 4 12 192 "Courier" 
12}{SCP_R 1 0 23{COM_R 2 0{TEXT 1 2608 "\"The Polhode Rolls witho
ut Slipping on the Herpolhode Lying in the Invariable Plane\"\015
\015The motion of a rigid body upon which no external forces (or \+
torques) act is a standard problem in\015intermediate classical m
echanics. Because the center of mass is unaccelerated, the proble
m may be\015treated in a reference frame in which it is and remai
ns at the origin, so the fundamental degrees of freedom  are thos
e describing the rotation of the body from a fixed initial positi
on, represented by an orthogonal\015matrix A. There are only thre
e independent degrees of freedom here, and the velocities can be \+
\015taken as the angular velocity vector $\\omega$. But the confi
guration space is the group manifold of the\015rotation group, a \+
three dimensional space but not a flat one, and the motion can be
 quite complex. \015\015In this simulation, the angular velocity \+
vector, expressed in body fixed coordinates, $\\omega_\{body\}$, \+
which is called Wbl in the program, has its motion determined by \+
Euler's equations, below. The time evolution of the matrix A is d
etermined by $dA_\{im\}/dt = \\epsilon_\{ijk\} Wbl_k A_\{jm\}$. V
ectors in the lab (inertial) frame are related to those in the bo
dy fixed frame by $V_\{body\} = A V_\{lab\}$.\015\015The problem \+
is somewhat simplified by the conservation of two quantities, the
 kinetic energy \015$T = \\omega\\cdot  I \\cdot \\omega / 2$, an
d the angular momentum $L = I \\cdot \\omega$, where I is the mom
ent of inertia tensor. The angular momentum is a vector with comp
onents conserved only in the inertial lab frame. Because of the e
nergy conservation, $\\omega$ at any instant is confined to an el
lipsoid given by the quadratic constraint. This ellipsoid is fixe
d in the body coordinate system (where the moment of inertia tens
or is fixed), and is called the inertia ellipsiod. On the other h
and, as $T = L \\omega/2$ and T and L are conserved in laboratory
 coordinates, one component of $\\omega$ is fixed in laboratory c
oordinates, confining $\\omega_\{lab\}$ to the invariable plane. \+
The locus of points mapped on the inerta ellipsoid by $\\omega_\{
body\}(t)$ is called the polhode, and the locus of points mapped \+
by $\\omega_\{lab\}(t)$ on the invariable plane is called the her
polhode. As the momentary position of $\\omega$ is on the axis of
 rotation at the moment, and points on this axis are momentarily \+
at rest, the polhode rolls without slipping on the herpolhode.\01
5\015This program generates an object with three independant prin
cipal moments of inertia, and shows the\015body, its inertia elli
psoid, and the polhode and herpolhode. It integrates the equation
s of motion, and can show the configuration at various times.\015
\015\015"}}{COM_R 3 0{TEXT 1 87 " \\hskip 2in                    \+
                                           Procedure:\015\015\015
"}}{COM_R 4 0{TEXT 1 644 "First, we load some libraries and defin
e parameters. We choose an initial orientation of the body and \0
15an Angular momentum (chosen in the -z direction). eps and ff ar
e fudge factors to place the polhode \015slightly outside the ine
rtia ellipsoid and the herpolhode slightly above the invariable p
lane, so that they\015can be seen. Deltat is the integration step
 size, Bx is the scale by which the object is enlarged to plot\01
5on axes whose units are angular velocities. The body is a unifor
m rectangular solid of dimensions\0152 sxb2  by  2 syb2 by 2 szb2
. \015\015Most important:  Adjust the movie parameters so you don
't overtax Maple. These are on next line."}}{INP_R 5 0 "> "{TEXT 
0 844 "startnum:=890; stopnum:=900; bynum:=10;\015Deltat:=.01;\01
2topnum:=1000;\012with(linalg):\012with(plots):\012initdeg:=89;\0
11# rotate z axis towards y by this angle, in degrees\012inittilt
:=evalf(Pi * initdeg/180);\012eps:=.001;\012ff:=0.999;\012Lr:=vec
tor(3,[0,0,-2]);\012A:=matrix(3,3,[1,0,0,0,cos(inittilt),sin(init
tilt),\\\012   0,-sin(inittilt),cos(inittilt)]);\012Bx:=4;\012dis
p:=vector(3,[-10,10,0]);\012Lb:=evalm(A &* Lr);\012sxb2:=evalf((3
/19)^(1/2));\012syb2:=evalf((17/19)^(1/2));\012szb2:=1;\012I11:=(
syb2^2+szb2^2)/3; I22:=(sxb2^2+szb2^2)/3; I33:=(sxb2^2+syb2^2)/3;
\012Wb:=vector(3,[Lb[1]/I11,Lb[2]/I22,Lb[3]/I33]);\012T:=dotprod(
Lb,Wb)/2;\012rt2TbI11:=evalf(ff*(2*T/I11)^(1/2));\012rt2TbI22:=ev
alf(ff*(2*T/I22)^(1/2));\012rt2TbI33:=evalf(ff*(2*T/I33)^(1/2));\
012bsph:=(theta,phi) -> [rt2TbI11 *sin(theta)*cos(phi),\\\012   r
t2TbI22*sin(theta)*sin(phi),\\\012   rt2TbI33*cos(theta)];\012Wr:
=evalm(transpose(A) &* Wb);\012"}}{COM_R 6 0{TEXT 1 120 "The abov
e has set up initial values. Next we give Euler's equations dewi.
 The number of time steps\015evaluated is topnum. "}}{INP_R 7 0 "
> "{TEXT 0 159 "dew1:=  diff(Wb1(t),t)= (I22-I33)*Wb2(t)*Wb3(t)/I
11;\012dew2:=  diff(Wb2(t),t)= (I33-I11)*Wb3(t)*Wb1(t)/I22;\012de
w3:=  diff(Wb3(t),t)= (I11-I22)*Wb1(t)*Wb2(t)/I33;\012"}}{COM_R 8
 0{TEXT 1 234 "Now we solve the differential equations (numerical
ly) for omega(t) in the lab system. and place a \015table of the \+
results in Wbl, with dimensions Wbl[topnum+1][3]. Then we integra
te the equations for\015A, using an improved Euler's method."}}
{INP_R 9 0 "> "{TEXT 0 163 "Wsol:= dsolve(\{dew1,dew2,dew3, \\\01
2\011Wb1(0)=Wb[1],Wb2(0)=Wb[2],Wb3(0)=Wb[3]\}, \\\012\011[Wb1(t),
Wb2(t),Wb3(t)],type=numeric, \\\012\011value=array([seq(Deltat*t2
,t2=0..topnum)])):\012"}}{INP_R 10 0 "> "{TEXT 0 499 "Wbl:=conver
t(submatrix(evalm(Wsol[2,1]),1..topnum+1,2..4),listlist):\012At[1
]:=evalm(A);\012\012for n from 1 to topnum do\012   told:=n*Delta
t:\012   tnew:=told + Deltat;\012   F:=evalm(matrix(3,3,[0,Wbl[n]
[3],-Wbl[n][2],-Wbl[n][3],0,Wbl[n][1], \\\012\011     Wbl[n][2], \+
-Wbl[n][1],0]) &* evalm(At[n])):\012   G:=evalm(matrix(3,3,[0,Wbl
[n+1][3],-Wbl[n+1][2],-Wbl[n+1][3],0,Wbl[n+1][1],\\\012\011     W
bl[n+1][2], -Wbl[n+1][1],0]) &* ( evalm(At[n] +evalm(Deltat*(F)))
)):\012   At[n+1]:=evalm(At[n] + evalm((Deltat/2)*evalm(F+G)));\0
12   od:\012"}}{COM_R 11 0{TEXT 1 112 "Now we generate the omega \+
table in the lab coordinates, called Wrl, and the transpose (inve
rse) of\015A, called ATt"}}{INP_R 12 0 "> "{TEXT 0 160 "Wrl:=vect
or(topnum+1);\012for n from 1 to topnum+1 do\012\011k:='k':\012\0
11ATt[n]:=transpose(evalm(At[n]));\012\011Wrl[n]:=[seq(sum(ATt[n]
[i,k]*Wbl[n][k],k=1..3),i=1..3)];\012    od:\012"}}{COM_R 13 0
{TEXT 1 130 "next group generates a herpolhode plot. remove the p
lotsetup command if you want it interactively\015instead of as a \+
postscript file."}}{INP_R 14 0 "> "{TEXT 0 144 "plotsetup(ps,plot
options=`color`,plotoutput=`herpy`.initdeg.`.ps`);\012plot(map(x-
>[op(1,x),op(2,x)],Wrl),scaling=CONSTRAINED,title=`Herpolhode`);\
012"}}{COM_R 15 0{TEXT 1 182 "next group generates a polhode plot
. The map of the ellipsoid into the plane is somewhat \015arbitra
ry - I map the north pole into an sllipse the shape of the ellips
oid's xy projection."}}{INP_R 16 0 "> "{TEXT 0 318 "rt1:=evalf((I
11 /(2*T))^(1/2));\012rt2:=evalf((I22 /(2*T))^(1/2));\012r2by1:=r
t2/rt1;\012rt3:=evalf((I33 /(2*T))^(1/2));\012phpj:=vector(topnum
+1);\012for n from 1 to topnum+1 do\012  phi := arctan( Wbl[n][2]
 * rt2, Wbl[n][1] * rt1 );\012  theta := arccos( Wbl[n][3] * rt3 \+
);\012  phpj[n]:=[theta * cos(phi), r2by1 * theta * sin(phi)];\01
2  od:\012"}}{COM_R 17 0{TEXT 1 90 "Plot the projected polhode. I
f you want it interactively, once again, remove the plotsetup"}}
{INP_R 18 0 "> "{TEXT 0 133 "plotsetup(ps,plotoptions=`color`,plo
toutput=`plhdy`.initdeg.`.ps`);\012plot(phpj,scaling=CONSTRAINED,
title=`Polhode mapped to plane`);\012\012"}}{COM_R 19 0{TEXT 1 
611 "We make  room for the polhode in lab space. \015It can take \+
a lot of time. You need to check where you want\015to start (star
tnum), where you want to stop (stopnum), and check that bynum, th
e number of time intervals between pictures drawn, is such that y
ou don't generate more pictures than your computer will be happy \+
with. Possible approaches to the pictures are to display them ind
ividually, as a movie from within maple, or saved as individual p
ostscript images for processing outside of maple.  If the movie i
s shown\015interactively, some changes need to be made, and the l
imitations on performance more noticable"}}{INP_R 20 0 "> "{TEXT 
0 1807 "Polhode:=vector(topnum+1);\012\012_OPT:=style=PATCH,scali
ng=CONSTRAINED,projection=1,axes=FRAME,\\\012labels=[`x`,`y`,`z`]
,orientation=[44,73],view=[-18..10,-10..15,-7.8..8];\012\012for n
 from startnum+1 by bynum to stopnum+1 do \012   lsph:=(theta,phi
) -> evalm(ATt[n] &* evalm(bsph(theta,phi)));\012   P1:=plot3d(\{
convert(lsph(u,v),list),\\\012      [6*u-3*Pi,3*v-3*Pi,Wr[3]-eps]
\},u=0..Pi,v=0..2*Pi,_OPT):\012#  translate body so as to not be \+
hidden\012   bside1:=(u,v)->evalm(disp+evalm(ATt[n] &* [Bx*sxb2*u
,Bx*syb2*v,-Bx*szb2]));\012   bside2:=(u,v)->evalm(disp+evalm(ATt
[n] &* [Bx*sxb2*u,Bx*syb2*v,Bx*szb2]));\012   bside3:=(u,v)->eval
m(disp+evalm(ATt[n] &* [Bx*sxb2*u,-Bx*syb2,Bx*szb2*v]));\012   bs
ide4:=(u,v)->evalm(disp+evalm(ATt[n] &* [Bx*sxb2*u,Bx*syb2,Bx*szb
2*v]));\012   bside5:=(u,v)->evalm(disp+evalm(ATt[n] &* [Bx*sxb2,
Bx*syb2*u,Bx*szb2*v]));\012   bside6:=(u,v)->evalm(disp+evalm(ATt
[n] &* [-Bx*sxb2,Bx*syb2*u,Bx*szb2*v]));\012   P4:=plot3d(bside1(
u,v),u=-1..1,v=-1..1,color=ORANGE,_OPT,grid=[2,2]);\012   P5:=plo
t3d(bside2(u,v),u=-1..1,v=-1..1,color=YELLOW,_OPT,grid=[2,2]);\01
2   P6:=plot3d(bside3(u,v),u=-1..1,v=-1..1,color=AQUAMARINE,_OPT,
grid=[2,2]);\012   P7:=plot3d(bside4(u,v),u=-1..1,v=-1..1,color=B
LUE,_OPT,grid=[2,2]);\012   P8:=plot3d(bside5(u,v),u=-1..1,v=-1..
1,color=PINK,_OPT,grid=[2,2]);\012   P9:=plot3d(bside6(u,v),u=-1.
.1,v=-1..1,color=RED,_OPT,grid=[2,2]);\012   P2:=spacecurve([seq(
Wrl[k],k=1..n)],color=BLACK,thickness=3,_OPT):\012   for k from 1
 to n do\012      Polhode[k]:=convert(evalm(ATt[n] &* Wbl[k]),lis
t); od:\012   P3:=spacecurve([seq(Polhode[k],k=1..n)],color=BLACK
,thickness=3,_OPT):\012#  for images within maple use       PP[n]
:=display3d(\{P.(1..9)\},_OPT):\012#  set up to record the pictur
es in postscript\012   plotsetup(ps,plotoptions=`color`,plotoutpu
t=`img`.initdeg.`y`.n.`.ps`);\012   print(display3d(\{P.(1..9)\},
_OPT));\012   print(n);\012od:\012"}}{COM_R 21 0{TEXT 1 73 "if yo
u want maple to show you a movie, and you generated PP[n] above, \+
use"}}{INP_R 22 0 "> "{TEXT 0 76 "# display3d([seq(PP[k1*bynum +1
],k1=0..(topnum/bynum))],insequence=true);\012\012\012"}}{COM_R 
23 0{TEXT 1 527 "To get a really good movie, it is best to genera
te many postscript files. These can be combined \015into an MPEG \+
using the \\verb\"mpeg_encode\" program with a filter that uses g
hostscript to generate\015ppmraw files, then pbm facilities to ro
tate and crop the images. A sequence of 251 frames was\015done th
at way, with the mpeg taking under an hour of Sparcstation 5 time
. Maple has trouble \015generating all the images at once, howeve
r, and it is best to do about 30 frames per session in\015Maple, \+
on my Sparcstation 5, which has 32MB ram. "}}{SEP_R 24 0}}{END}