Tumbling Video (experiment)

My friend Van of Warren Design Visions told me to post my tumbling video on youtube. Let's hope this works. If the embedded video doesn't show up, try

http://www.youtube.com/watch?v=Zne5cVXzPtA

A Famous Instability

Take a hardcover book that you don't mind dropping on the floor and put a couple of rubber bands around it to keep it closed. Observe that there are three axes of symmetry of the book: longest, shortest, and intermediate. Toss the book in the air so that it spins about the shortest axis. That's the axis that passes from the front of the book to the back (unless you have picked a very strangely shaped book). Try to catch it, too, but notice that the book starts off spinning around this axis, and continues spinning around this axis. This is a stable spinning motion. Now toss the book in the air so that it spins around the longest axis, most probably the axis from the top of the book to the bottom. Stable. Now, one more time, spin about the middle axis. Whoops! It won't stay spinning about that axis, but tumbles all over the place. UNSTABLE about the middle axis.

Somewhere in my vast library I have a theorem about this, which I'll try to resurrect for another blog. But this time I just want to point out that I've reproduced this behavior in my little Mathematica laboratory of rigid-body motion. I'm not sure how to publish a Flash movie on this blog, or, in fact, anywhere else on the web, but I do have one. Here's a snapshot:

I've got three versions of this movie, one showing stable spinning about the x (red) axis, another showing stable spinning about the y (green) axis, and the third showing UNstable spinning about the z (blue) axis. All three of them use the same 4th-order Runge-Kutta integrator of the rigid-body equations of motion for free spinning, which I coded up by hand in Mathematica. I'll be explaining this in future installments.

Fun with Shapes

rcoord:={Random[],Random[],Random[]}

eloof[f_]:=TranslateShape[
rotw[(2rcoord-1)Pi,AffineShape[wcube,rcoord]]
, f rcoord];

Depicting Rotations via Euler Angles

Here are the promised OBJECT-rotation operators. This actually took me much longer to explain than to create.

Using Mathematica once again, Consider a vector w representing angular velocity by its direction and magnitude. The direction specifies the axis of rotation and the magnitude specifies a positive twist rate in radians per unit time about this axis via the right-hand rule: thumb pointing along the axis, fingers curling in the sense of positive twist. To animate rotations, multiply w by dt to get dw, a small but finite actual rotation. Use "RotateShape," which documentation says "takes Euler angles," to depict such rotations.

Given dw's components, then, what arguments to supply to "RotateShape?" There are lots of different Euler-angle sets. This one turns out to be zxz frame-wise:

xaxis=Line[{{0,0,0},{1,0,0}}];
yaxis=Line[{{0,0,0},{0,1,0}}];
zaxis=Line[{{0,0,0},{0,0,1}}];

Look at the x axis after rotating about +z counterclockwise by a small angle, say p/6. If RotateShape rotates the FRAME, expect a small negative y contribution. If RotateShape rotates the OBJECT, expect a small positive y contribution:


Look at the y axis after a rotation about +x by a small angle. Expect small negative z contribution:

Picture an all-positive w, in the first octant of E3. To rotate the image of an object about such a w via "RotateShape," first get the frame's y axis directly under w by a negative twist about z, then do a negative twist about new x to get z aligned with w, then twist by negative mag[w], to twist the OBJECT rather than the frame, then undo the first two twists.

phi[w] is the angle w makes with the y axis, measured CLOCKwise from y to w in the xy plane. The first Euler angle should be a negative z twist by this angle.

theta[w] is the angle w makes with the z axis, measured counterclockwise from z to w in their common plane. The second Euler angle should be a negative x twist to align the new z axis with w.

Here's the rotation op:

And an animator for it:

And some evidence that it's working right

where

s3p=Show[Graphics3D[#
,FaceGrids->All
,Axes->True
,AxesLabel->{"x","y","z"}
,Lighting->False]
,ImageSize->Automatic]&;

This is good enough to support animation of rigid bodies in motion.

Models for Orientation Gadgets

Here is some Mathematica code, without explanation, for improved orientation-visualization gadgetry. I'll do the display operators later, but I'm rushing out the door:

vcube=
TranslateShape[
{FaceForm[RGBColor[0,0,1],RGBColor[1,1,0]]
,Polygon[{{0,0,0},{1,0,0}(*xy*)
,{1,1,0},{0,1,0}}]
,Polygon[{{0,0,1},{1,0,1}
,{1,1,1},{0,1,1}}]
,FaceForm[RGBColor[0,1,0],RGBColor[1,0,1]]
,Polygon[Reverse[{{0,1,0},{1,1,0}
,{1,1,1},{0,1,1}}]](*zx*)
,Polygon[Reverse[{{0,0,0},{1,0,0}
,{1,0,1},{0,0,1}}]]
,FaceForm[RGBColor[1,0,0],RGBColor[0,1,1]]
,Polygon[{{0,0,0},{0,1,0}
,{0,1,1},{0,0,1}}]
,Polygon[{{1,0,0},{1,1,0},{1,1,1}
,{1,0,1}}]
},{-.5,-.5,-.5}];

x2pols={FaceForm[RGBColor[0,0,0],RGBColor[.5,.5,.5]]
,Polygon[
Reverse[{{1,1,0},{2.5,6,0},{3,6,0},{3,4+1/3,0},{2,1,0}}/12.]]
,Polygon[
Reverse[{{4,1,0},{3,4+1/3,0},{3,6,0},{3.5,6,0},{5,1,0}}/12.]]
,Polygon[
Reverse[{{3,6,0},{2.5,6,0},{1,11,0},{2,11,0},{3,7+2/3,0}}/12.]]
,Polygon[
Reverse[{{3,6,0},{3,7+2/3,0},{4,11,0},{5,11,0},{3.5,6,0}}/12.]]
};

y2pols={FaceForm[RGBColor[0,0,0],RGBColor[.5,.5,.5]]
,Polygon[
Reverse[{{2.5,1,0},{2.5,6,0},{3.5,6,0},{3.5,1,0}}/12.]]
,Polygon[
Reverse[{{3,6,0},{2.5,6,0},{1,11,0},{2,11,0},{3,7+2/3,0}}/12.]]
,Polygon[
Reverse[{{3,6,0},{3,7+2/3,0},{4,11,0},{5,11,0},{3.5,6,0}}/12.]]
};

z2pols={FaceForm[RGBColor[0,0,0],RGBColor[.5,.5,.5]]
,Polygon[
Reverse[{{1,1,0},{1.3,2,0},{5,2,0},{5,1,0}}/12.]]
,Polygon[
Reverse[{{1.3,2,0},{3.7,10,0},{4.7,10,0},{2.3,2,0}}/12.]]
,Polygon[
Reverse[{{1,11,0},{5,11,0},{4.7,10,0},{1,10,0}}/12.]]
};

combp=Join[#1,TranslateShape[#2,{.5,0,0}]]&;

xy2p=combp[x2pols,y2pols];
yz2p=combp[y2pols,z2pols];
zx2p=combp[z2pols,x2pols];

face2labels=
TranslateShape[
{TranslateShape[RotateShape[zx2p,Pi,Pi/2,0], {1,0,0}]
,TranslateShape[RotateShape[zx2p,Pi,Pi/2,0],{1,1,0}]
,RotateShape[yz2p,0,-Pi/2,-Pi/2]
,TranslateShape[RotateShape[yz2p,0,-Pi/2,-Pi/2],{1,0,0}]
,xy2p
,TranslateShape[xy2p,{0,0,1}]},{-.5,-.5,-.5}];

wcube=Chop[Join[AffineShape[vcube,{0.99,0.99,0.99}],face2labels]];

thing3={Thickness[0.010]
,RGBColor[1,0,0]
,Line[{{0,0,0},{1,0,0}}]
,RGBColor[0,1,0]
,Line[{{0,0,0},{0,1,0}}]
,RGBColor[0,0,1]
,Line[{{0,0,0},{0,0,1}}]
,RGBColor[0,1,1]
,Line[{{0,0,0},{-1,0,0}}]
,RGBColor[1,0,1]
,Line[{{0,0,0},{0,-1,0}}]
,RGBColor[1,1,0]
,Line[{{0,0,0},{0,0,-1}}]
};

xcube=Join[AffineShape[thing3,{2,2,2}],wcube];

Visualizing Orientation in E3

I'm working toward a little laboratory for rigid-body simulation. For that, we need some more-or-less generic rigid bodies and a way of visualizing their positions and orientations in Euclidean 3-space (E3). E3 is good for now; we might get to curved spacetime later.

To get started, let's use Mathematica. Right out of the box, investigating 3D plotting from the help file, do this:

rcoord := {Random[], Random[], Random[]}
blarg[f_] := Module[{loc = f rcoord, siz = rcoord},
    Cuboid[loc,
      Table[loc[[i]] + siz[[i]], {i, 3}]]]
Show[Graphics3D[Table[blarg[3], {20}]], FaceGrids -> All, 

Axes -> True, AxesLabel -> {"x", "y", "z"},
    ViewPoint -> {1.3, -2.4, 2.}];

Great. We have a way of spraying cubes around E3, and we can see their positions. Now, we need a way to track their orientations. The help file says "RotateShape" is the way to go. I found out by experimentation that "Cuboid" doesn't rotate like other graphics, so I made my own.

I wanted a color coding for the orientation of each of the six planes of a cube, so I came up with this. Let the unit vectors of the cube-fixed reference frame be x, y, z; in that order. Assign to their positive directions the primary colors, Red, Green, and Blue (RGB), in that order. Assign to their negative directions the complementary colors, Cyan, Magenta, and Yellow, respectively (CMY). The following illustrates:

thing2 = {Thickness[.02]
      , RGBColor[1, 0, 0]
      , Line[{{0, 0, 0}, {1, 0, 0}}]
      , RGBColor[0, 1, 0]
      , Line[{{0, 0, 0}, {0, 1, 0}}]
      , RGBColor[0, 0, 1]
      , Line[{{0, 0, 0}, {0, 0, 1}}]
      , RGBColor[0, 1, 1]
      , Line[{{0, 0, 0}, {-1, 0, 0}}]
      , RGBColor[1, 0, 1]
      , Line[{{0, 0, 0}, {0, -1, 0}}]
      , RGBColor[1, 1, 0]
      , Line[{{0, 0, 0}, {0, 0, -1}}]
      };

Embrace the right-hand rule and the wedge product for 2-forms, that is, for oriented planes, defining positive orientation for planes. Because of the lucky accidents that in E3 there are three basis vectors AND three planes, x^y, y^z, and z^x, plus three primary colors AND three secondary colors, we can unify color coding for vectors and planes:

x ~ y^z ~ Red   ~ {1,0,0},  -x ~ z^y ~ Cyan    ~ {0,1,1}
y ~ z^x ~ Green ~ {0,1,0},  -y ~ x^z ~ Magenta ~ {1,0,1}
z ~ x^y ~ Blue  ~ {0,0,1},  -z ~ y^x ~ Yellow  ~ {1,1,0}

Just to nail this down, hold up your right hand, thumb toward your face, and imagine the wrist bones on your pinky side touching the origin of the y^z plane, so that the plane's y basis vector points to the right and the plane's z basis vector points upward.

Let your fingers curl or sweep from the y direction to the z direction. Then, your thumb represents the x direction. If your thumb points at your face, then you should see Red, and if your thumb points away from your face, then you should see Cyan. That little procedure accounts for the first line above. Repeat for the other two lines, then convince yourself that you understand the following pictures. Their construction is a little too lengthy to reproduce in a blog entry.

The final image below shows the body rotated with respect to the "world frame" and sets the stage for visualizing Newtonian rigid-body dynamical simulation. Again, make sure you can visually "parse" this image. The Green plane represents the body-fixed positive y ~ z^x axis and plane pointing at you. The Cyan plane represents the body-fixed negative x ~ y^z axis-plane pointing at you, or the positive ones pointing away from you. The Yellow plane represents the body-fixed positive z ~ x^y axis-plane pointing away from you, or, likewise, the negative ones pointing toward you.

Euclid Tougher than Minkowski

We've done a few Geometry Expressions demos showing off Minkowski geometry and how it solves problems in Special Relativity. How about Euclidean geometry? The usual 2D rotation diagram is easy enough to set up. Let point B have coordinates x,y. Let the red frame be rotated by angle theta with respect to the blue frame and ask "What are the coordinates of B in the red frame?" We all remember the answer: x' = x cos theta + y sin theta, y' = - x sin theta + y cos theta. What does GX say?


JF is obviously (BD = x) cos theta, as is IC = (AC = x) cos theta, just by looking at the picture. Ditto HF = (BC = y) sin theta, and the y' coordinate is just as easy to see. Mathematica's simplification machinery is better than GX's, and it's a simple matter to copy-as and paste back to Mathematica,


I'm getting quite accustomed to this procedure. It's amusing as well as useful.

GX file for you to download is at http://home.comcast.net/~brianbec/Euclidean3.gx

Velocity Addition for Free

Ever go through the traditional derivation of the relativistic formula for addition of velocities? Painful, because there are a bunch of things you have to hold in your head for the duration of a calculation. Modern tools make this easy, almost obvious. First, recall this sketch:

in which Geometry Expressions calculates the Lorentz transformation for free. Now, suppose the point B really represents the endpoint of a particle path, the beginning-point being the origin (0,0). In the blue frame -- our lab frame -- the (average) velocity of the particle over this path is x/t, which we can write as u. What's the velocity of the same particle path, this time measured in the red frame, moving with velocity v w.r.t. our blue frame? Why, let's let Mathematica compute it for us. Use the "Edit / Copy As / Mathematica" menu item in Geometry Expressions on the yellow expressions in the sketch above, and just paste into Mathematica:

Beautiful simplicity itself, eh? Once again, I've been lazy with absolute-value expressions, just sidestepping them. This means the signs -- but only the signs -- may turn out wrong and we would have to use common sense to correct them.