go to next section ("Physics")
I used this skeleton for my front-end project, which uses JavaScript for the calculations and ReactJS (with functional components and hooks) to render the results. The most important piece of this simulation is the clock, inspired by this and implemented as follows:
useEffect(() => {
let interval;
if (running) interval = setInterval(() => setTime(time + dt/1000), dt);
if (!running && time !== 0) clearInterval(interval);
return () => clearInterval(interval);
}, [running, time, dt]);
... with the running condition controlled by this button:
<button onClick={() => setRunning(!running)}>{running ? "Stop" : "Start"}</button>
The most basic user-control for this simulation is the value of the requisite "time-step" dt, for which the following pair of considerations must always be balanced:
- Too large will lead to choppy animations and inaccurate calculations.
- Too small will lead to a backup of the message queue which'll manifest itself by the clock's running slowly.
To make the animation less choppy (even for larger values of dt), I use css transitions (with transition-timing-function: linear) to interpolate the object's orientation between successive timesteps in order to allow for large values of Δt without sacrificing animation smoothness. Below is the relative code from the Body component
<div className="body"
style={{ transitionDuration: `${!running ? 0 : dt / 1000}s` }}
>
... and from index.css:
.body {transition-timing-function: linear; transition-property: all; }
In the absence of air resistance (not present in this simulation) and aforementionned errors associated with the nonzero size of the time-step dt, the kinetic energy (monitored at the bottom of the page) should be constant. Hence monitoring the constancy of this energy is the best way for the user to determine if dt is sufficiently small.
return to previous section ("Introduction")
go to next section ("Linear algebra")
Just as the position of an object is specified by the values of three numbers (the x-, y-, and z- coordinates of the center of mass), the angular orientation of a rigid body is specified by three numbers: the Euler angles, commonly designated by the Greek letters φ, θ, and ψ. In the absence of any forces, the translational motion of an object is quite simple: the values of x, y, and z each change at a constant rates. These (independent) rates are the Cartesian components of the object's velocity. In a similar manner the rotation of a rigid body involves the rates of change of the Euler angles. Unlike for translational motion, however, these rates of change are not constant, even when there is no torque. Rigid-body rotation is much more complicated than rotation about a fixed axis, examples of the latter being the motion of a globe or of a hinged door.
Of fundamental importance in rotations is the inertia matrix. This 3 x 3 matrix is usually expressed in terms of three quantities: the moment of inertia for rotation about each of the body's three principal axes. The three principal axes are mutually perpendicular and - for symmetric objects - point along the object's axes of symmetry. In my app the object is a rectangular cuboid, an example of such a symmetric structure. For these shapes the moment of inertia equal m(a2 + b2)/12, in which m is the mass, and a and b are the dimensions of the cuboid face which is perpendicular to the particular axis. However in my app the user specifies the moments of inertia, not the face sizes, so I need to invert that formula in order to obtain cuboid dimensions d which I can use for rendering the cuboid. To invert the formulas I set m/12 = 1 (to simplify the arithmetic) and arrive at the result that the length of any cuboid axis equals the square root of the difference between half the sum of the (three) moments of inertia minus the moment for rotation about that particular axis. Below are the lines of code from my calcD helper function in the App component.
let sumMom = moms[0] + moms[1] + moms[2];
let newD = moms.map(mom => Math.max(0.000001, Math.sqrt((sumMom / 2 - mom))));
Because singularities occur if any cuboid dimension equals zero (as would be the case for a flat shape cut from sheet metal of infinitessimal thickness), I manually prevent this thru the use of Math.max. Note that one physical constraint which must be satisfied by the moments of inertia is that no single moment of inertia may exceed the sum of the other two. The line of code below (from the input handler handlerMom) makes - and sets in state - a determination of whether or not that constraint has been violated for the moments of inertia moms
setAreLegalMoms(newMoms.reduce((legal, mom, i, moms) => (legal && mom <= (moms[(i+1)%3] + moms[(i+2)%3])), true));
If such a violation occurs, the user is warned accordingly:
{areLegalMoms ? null :
<div className="message">
No single moment of inertia should exceed the sum of the other two.
</div>
}
The formula (found on p. 7 of the linked reference) for the rate of change of the Euler angles is a nonlinear function of the angles themselves. The code for this function is below, in which ths is the three-component array with the instantaneous values of the radian measure of the Euler angles, moms is the three-component array of the moments of inertia, and Lz is the z-component of the angular momentum. (The other two components of the angular momentum are zero.)
const Fs = ths => {
let cs = [];
let ss = [];
for (const th of ths) {
cs.push(Math.cos(th));
ss.push(Math.sin(th));
};
let Fs = []
Fs[0] = Lz * (cs[2] * cs[2] / moms[1] + ss[2] * ss[2] / moms[0]);
Fs[1] = Lz * (1 / moms[0] - 1 / moms[1]) * ss[1] * ss[2] * cs[2];
Fs[2] = Lz * (1 / moms[2] - cs[2] * cs[2] / moms[1] - ss[2] * ss[2] / moms[0]) * cs[1];
// many unrelated lines removed for the sake of brevity
return Fs;
}
This constitutes a system of three nonlinear first-order differential equations, which I solve using a 4th-order Runge-Kutta method.
useEffect(() => {
if (!running) return;
let Fs1 = Fs(ths);
let Fs2 = nextFs(Fs1, 2);
let Fs3 = nextFs(Fs2, 2);
let Fs4 = nextFs(Fs3, 1);
setThs([...ths].map((th, i) => th + (Fs1[i] + Fs4[i] + 2 * (Fs2[i] + Fs3[i])) * dt/ 1000 / 6));
}, [time, running]);
Because this method is 4th-order, the algorithm's time complexity is proportional to the inverse of the fourth root of the simulation's desired accuracy for ths. Put differently, the accuracy of its calculation of ths is proportional to the fourth power of the user-specified timestep dt.
return to previous section ("Physics")
go to next section ("3-d Rendering")
The instantaneous value of the cuboid's rotation is characterized by a 3 x 3 "total" rotation matrix, which itself is the product of three simpler rotation matrices], each parametrized by one of the three Euler angles. Below are the simpler rotation matrices, first for rotation about the z-axis (used twice) and second for rotation about the x-axis (used once).:
const zRot = th => {
let [c, s] = [Math.cos(th), Math.sin(th)];
return [[c, s, 0], [-s, c, 0], [0, 0, 1]];
}
const xRot = th => {
let [c, s] = [Math.cos(th), Math.sin(th)];
return [[1, 0, 0], [0, c, s], [0, -s, c]];
}
For subsequent steps are needed four helper functions for operations of linear algebra: scalar multiplication (also known as "dot product") of two vectors, multiplication of one matrix and one vector, transposition of a matrix, and multiplication of two matrices. Below are my codes for these, each accomplished with one or two list comprehensions.
const dotproduct = (vec1, vec2) => vec1.reduce((dot, comp, i) => dot + comp * vec2[i], 0);
const mult1 = (mat, vec) => mat.map(row => dotproduct(row, vec));
const transpose = mat => mat[0].map((blah, i) => mat.map(row => row[i]));
const mult2 = (mat1, mat2) => mat1.map(x => transpose(mat2).map(y => dotproduct(x, y)));
Below is the calculation of the total rotation matrix, the return from which I will henceforth simply call the "rotation matrix" and symbolize by mat.
const invRot=ths=> mult2(mult2(zRot(-ths[0]),xRot(-ths[1])), zRot(-ths[2]));
Any rotation matrix such as mat may be characterized by two things: an axis (of rotation) and an angle, and I must obtain these two things to render the rotated cuboid using the css function rotate3d(...). The absolute value of the angle can be obtained easily from the trace of mat, as seen in the following lines from the helper function rotate.
let trace = mat[0][0] + mat[1][1] + mat[2][2];
let angle = Math.acos((trace - 1) / 2);
It is less straightforward to calculate the rotation axis, which can be any non-null vector that is parallel to the eigenvector of mat whose eigenvalue equals 1. To find this eigenvector I use the ml-matrix package.
import { EigenvalueDecomposition, Matrix } from "ml-matrix";
The following line is another from the helper function rotate.
let vectors = new EigenvalueDecomposition(new Matrix(mat)).eigenvectorMatrix.transpose().data;
The vectors array contains the three eigenvectors of mat, but I only need the eigenvector whose eigenvalue equals 1, ie the vector vec with the property that mat·vec = vec. This is accomplished by the following lines in rotate.
let dVectors = vectors.map(vector => mult1(mat, vector).map((comp, i) => comp - vector[i]));
let mags = dVectors.map(dVector => dVector.reduce((mag, comp) => mag + comp * comp, 0));
let min = mags.reduce((min, mag, i) => mag < min[1] ? [i, mag] : min, [-1, Infinity]);
let axisVec = vectors[min[0]];
Recall that the earlier formula used for the rotation angle only provides its absolute value. In order to determine whether to multiply its absolute value by -1 I need to (1) define a vector (Vec) which is perpendicular to the rotation axis, (2) rotate that vector by multiplying it by the rotation matrix mat, (3) evaluate the cross-product of those two vectors, which'll then point either parallel to or antiparallel to the rotation axis, and (4) determine the dot product of that cross-product with the rotation-axis vector. The dot product will be either +/-1, ie the factor by which I'll need to multiply angle.
let vec = vectors[(min[0] + 1) % 3];
let rVec = mult1(mat, vec);
// cross product of rVec with vec
let rVecCrossVec = [rVec[1] * vec[2] - rVec[2] * vec[1],
rVec[2] * vec[0] - rVec[0] * vec[2],
rVec[0] * vec[1] - rVec[1] * vec[0]];
angle *= Math.sign(dotproduct(axisVec, rVecCrossVec));
return to previous section ("Linear algebra")
go to next section ("Examples")
The coordinate system used for this project has x pointing to the right and y pointing down. In order to maintain a right-handed coordinate system, I choose the z-axis to point into the screen. Hence a positive value for the angular momentum means that the object will rotate clockwise. In order to render the rotation within this coordinate system I followed this example. The App component contains the following lines of code to render the rotating cuboid
<div className="container" style={{
height:`${npx}px`,
width:`${npx}px`,
perspective:`${perspective}px`
}}>
<Dot x={npx/2} y={npx/2} d={10} />
<Line npx={npx} r={omfLat * npx / 2} angle={omfAng} dt={dt} time={time} />
<Body npx={npx} angleVec={angleVec} d={d} dt={dt} mids={mids0} degeneracies={degeneracies} running={running} />
</div>
... here is from index.css: .body {position: absolute; transform-style: preserve-3d; }
The extent of the cuboid's [graphical perspective](https://en.wikipedia.org/wiki/Perspective_(graphical)#:~:text=Perspective%20works%20by%20representing%20the,seen%20directly%20onto%20the%20windowpane.) is achieved straightforwardly via the parent div's <tt>perspective</tt> property, which takes a value set by the user. This value is manifest by how far away parallel lines seem to converge together (at a "vanishing point").
Note above that the parent div is positionned relatively so that absolute positionning can be used by each of its three child components: <tt>Dot</tt>, <tt>Line</tt>, and <tt>Body</tt>. The <tt>Dot</tt> component simply depicts the cuboid's [center of mass](https://en.wikipedia.org/wiki/Center_of_mass), and I'll describe the other two components next.
The <tt>Line</tt> component depicts the cuboid's instantaneous angular-velocity vector. (Of course it only represents the the *x*- and *y*-components of the vector.) The code below is taken from <tt>Fs</tt>, shown earlier to calculate the derivatives of the Euler angles. This helper function is also convenient for calculating the body's angular velocity (typically symbolized by the Greek letter ω). The letters <tt>Omf</tt> stand for "omega in the fixed frame". The user is in a fixed frame of reference, and I want to calculate ω in the same frame.
let newOmfs = []; newOmfs[0] = Fs[2] * ss[1] * ss[0] + Fs[1] * cs[0]; newOmfs[1] =-Fs[2] * ss[1] * cs[0] + Fs[1] * ss[0]; newOmfs[2] = Fs[2] * cs[1] + Fs[0]; setOmfs(newOmfs); let newOmf = Math.sqrt(newOmfs.reduce((om2, om) => om2 + om * om, 0)); setOmf(newOmf); setomfLat(Math.sqrt(newOmfs[0] * newOmfs[0] + newOmfs[1] * newOmfs[1]) / newOmf); let newOmfAng = Math.atan2(omfs[1], omfs[0]);
The <tt>Lat</tt> suffix signifies the magnitude of ω's lateral components, which is used (along with the angular direction of the vector, calculated trigonometrically above) for rendering the angular velocity as a line segment emanating from the cuboid's center of mass.
Below is <tt>App</tt>'s call to the <tt>Line</tt> component.
<Line npx={npx} r={omfLat * npx / 2} angle={omfAng} dt={dt} time={time} />
Below is the relevant code in the <tt>Line</tt> component, which uses props and the CSS functions <tt>rotate</tt> and <tt>translateX</tt> to compute the instantaneous value of the CSS <tt>transform</tt> property when rendering the lateral components of ω as a rotated line segment.
const Line = ({ npx, r, angle, dashed, dt }) => {
return (
<div className="line" style={{
width:${r}px,
left: ${npx / 2 - r / 2}px,
transform: rotate(${angle * 180 / Math.PI}deg) translateX(${r / 2}px),
}}/>
)
}
The <tt>Body</tt> component contains - first and foremost - a parent div whose opening tag is below.
The three pairs of the cuboid's parallel faces are each initially perpendicular to the z- (red), x- (green), and y-axes (blue) respectively of the fixed frame. The values of the width and height properties reflect the pixel-dimensions of the particular face of the cuboid. The values of the left and top properties are then set in order to center the particular face properly within the parent div. The value of the transform property is then set by the appropriate translate# and rotate# function in order to rotate and translate the face into its appropriate start position and orientation. Subsequent rotation is controlled by the transform property of the parent div within the Body component, utilizing the rotate3d function shown earlier.
<div className="side"
style={{
width:`${2 * d[0]}px`, height:`${2 * d[1]}px`,
left: `${-d[0]}px`, top: `${-d[1]}px`,
transform: `translateZ(${d[2]}px)`, background: "rgba(100,30,30,0.8)",
}}
/>
<div className="side"
style={{
width:`${2*d[0]}px`, height:`${2*d[1]}px`,
left: `${-d[0]}px`, top: `${-d[1]}px`,
transform: `translateZ(${-d[2]}px)`, background: "rgba(100,30,30,0.8)",
}}
/>
<div className="side"
style={{
width:`${2 * d[2]}px`, height:`${2*d[1]}px`,
left: `${-d[2]}px`, top: `${-d[1]}px`,
transform: `rotateY(90deg) translateZ(${d[0]}px)`, background: "rgba(40,100,40,0.8)",
}}
/>
<div className="side"
style={{
width:`${2*d[2]}px`, height:`${2*d[1]}px`,
left: `${-d[2]}px`, top: `${-d[1]}px`,
transform: `rotateY(-90deg) translateZ(${d[0]}px)`, background: "rgba(40,100,40,0.8)",
}}
/>
<div className="side"
style={{
width:`${2*d[0]}px`, height:`${2*d[2]}px`,
left: `${-d[0]}px`, top: `${-d[2]}px`,
transform: `rotateX(90deg) translateZ(${d[1]}px)`, background: "rgba(50,50,100,0.8)",
}}
/>
<div className="side"
style={{
width:`${2*d[0]}px`, height:`${2*d[2]}px`,
left: `${-d[0]}px`, top: `${-d[2]}px`,
transform: `rotateX(-90deg) translateZ(${d[1]}px)`, background: "rgba(50,50,100,0.8)",
}}
/>
The six remaining children depict the cuboid's principal axes in an analogous manner.
{degeneracies[0] ? null
:<>
<div className="line" style={{
left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
<div className="line" style={{
transform: `rotateX(90deg)`, left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
</>
}
{degeneracies[1] ? null
:<>
<div className="line" style={{
transform: `rotateZ(90deg)`, left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
<div className="line" style={{
transform: `rotateZ(90deg) rotateX(90deg)`, left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
</>
}
{degeneracies[2] ? null
:<>
<div className="line" style={{
transform: `rotateY(90deg)`, left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
<div className="line" style={{
transform: `rotateY(90deg) rotateX(90deg)`, left: `${-dmin}px`, width: `${2 * dmin}px`
}} />
</>
}
While it may seem sufficient to render each axis with a single div with a single 1-px border, I discovered that the subsequent "strip" would disappear when viewed "edge-on". Accordingly I added a second div rotated 90 degrees from the first one, so that both divs would never be edge-on at the same time.
An axis will not render if it is "degenerate" with another, by which I mean if the two axes have identical moments of inertia. In this case it is not appropriate to speak in terms of separate principal axes, so it's best not to render them at all. The code below determines whether or not a particular principal axis is degenerate with another.
let newDegeneracies = newMoms.map((momI, i) => {
return newMoms.reduce((degenerate, momJ, j) => {
return degenerate || (momJ === momI && i !== j);
}, false);
})
return to previous section ("3-d rendering")
Several examples are worth examining because their qualitative results are often described in theoretical discussions of this system. Each example differs in terms of its shape (isotropic, axisymmetric, or asymmetric), the particular axis which is chosen to be initially close (in direction) to the z-axis, and the initial value of θ.
The body undergoes pure rotation with θ = 0 (always), and ω is constant and always pointing parallel to the angular momentum (ie, into the page, which makes the lateral component of ω invisible).
The body's motion is qualitatively idential to the previous example.
The body undergoes periodic "precession" around the z-axis. θ never changes, so the precession is "stable".
The body undergoes unstable precession, because θ diverges.
The body undergoes precession which is either stable or unstable, considering the particular axis which is near the z-axis, as described in the next two bullets. See this for a thorough discussion, including a fascinating video from a zero-g environment.
θ changes but does not diverge, so the precession is stable.
The body's motion is chaotic, described better as "tumbling" rather than "precession".