Some Aspects of Rigid Body Motion

Off[General::spell1]
Off[Solve::ratnz]

rotatedBox[q_] :=  {
		{Red,Thickness[0.01],Line[{{0,0,0},
				 q.{10,0,0}}]},
  	       {Green, Thickness[0.01],Line[{{0,0,0},
				 q.{0,10,0}}]},
		{Blue, Thickness[0.01],Line[{{0,0,0},
                                     q.{0,0,10}}]},

		{Yellow,Polygon[q.#&  /@ {{2,6,-8},{2,6,8},{2,-6,8},{2,-6,-8}}
	                		]},
		{Yellow,Polygon[q.#&  /@ {{-2,-6,-8},{-2,-6,8},{-2,6,8},{-2,6,-8}}
	                		]},
		{Yellow,Polygon[q.#&  /@ {{2,6,-8},{-2,6,-8},{-2,6,8},{2,6,8}}
	                		]},
		{Yellow,Polygon[q.#&  /@ {{2,-6,-8},{2,-6,8},{-2,-6,8},{-2,-6,-8}}
	                		]},
		{Yellow,Polygon[q.#&  /@ {{2,6,-8},{2,-6,-8},{-2,-6,-8},{-2,6,-8}}
	                		]},
		{Yellow,Polygon[q.#&  /@ {{2,6,8},{-2,6,8},{-2,-6,8},{2,-6,8}}
	                		]}
                                              };

pointList= Flatten[Table[{a,b,c},{a,-2,2,.2},{b,-6,6,.6},{c,-8,8,.8}],2];
nPoints=Length[pointList] ;
massList=Table[1,{nPoints}];
stdBasis={{1,0,0},{0,1,0},{0,0,1}};

IT=Chop[Table[(∑_(k=1)^nPoints)massList[[k]](((Norm[pointList[[k]]])^2)KroneckerDelta[i,j]-(pointList[[k]].stdBasis[[i]])( pointList[[k]].stdBasis[[j]])),{i,1,3},{j,1,3}]];

{eigs,P}=Eigensystem[IT];
I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];

Note : To use this notebook, click on Evaluation>Evaluate Notebook.  Answer yes to any questions.  Then, as you go through the notebook, click on buttons as you encounter them.

Rigid Body Motion with a Point Fixed

A solid body ℬ moves in space so that one point in the body (taken to be the origin) remains fixed. 

Consider an initial position of ℬ and a final position of ℬ:  there is a mapping f  that takes each point in ℬ from its initial to its final position. 

Theorem (Euler, 1775)  The mapping f is accomplished by a rotation about some axis through the origin.

To rephrase:  There is an orthogonal matrix A with determinant 1 such that f(x)=A x for all x.

Kinematics of Rigid Body Motion

A rigid body begins moving at t=0, still assuming the origin is left fixed by the motion.  

The initial position of the body we call the reference configuration. Let r(t) be the position at time t of a particular point in the body.  

Initially, this point is located at r(0) and at time t>0 it is at r(t).  According to Euler, there is a rotation matrix A(t) such that

		r(t)=A(t) r(0).				(1)
						
We have
		A(0)=I            (I is the 3×3 identity matrix).
						
Differentiate (1) to get

		r'(t)=A'(t) r(0)
           =A'(t)(A^(-1))(t) r(t)
           =A'(t)(A^T)(t) r(t).
						
We define
		B(t)=A'(t)(A^T)(t),
						
and conclude that the position vector r(t) of a particle in the rigid body satisfies the system of differential equations

		r'(t)=B(t) r(t).

Recall:  matrix P is antisymmetric  if (P^T)=-P. 
					
Lemma:  For each t, B(t) is an antisymmetric matrix   ( a matrix P is antisymmetric  if (P^T)=-P).

Proof:  We know that A(t) is orthogonal for each t.  Therefore we have A(t) (A^T)(t)=I for all t.  Differentiating this equation gives

		A'(t)(A^T)(t)+A(t) ((A^T))'(t)=0.
Therefore,

		B(t)=-A(t) ((A^T))'(t).
											
Then (omitting the t to simplify notation)

		(B^T)=(((-A((A^T))'))^T)
       =-(((((A(A'))^T)))^T)
       =-(((A'(A^T)))^( ))
       =-B.
	
Since B(t) is antisymmetric, we have functions (Subscript[ω, 1])(t),(Subscript[ω, 2])(t), and (Subscript[ω, 3])(t), such that

		B(t)=(({{0, -(Subscript[ω, 3])(t), (Subscript[ω, 2])(t)}, {(Subscript[ω, 3])(t), 0, -(Subscript[ω, 1])(t)}, {-(Subscript[ω, 2])(t), (Subscript[ω, 1])(t), 0}})).	
								
If the position vector of a particle in the rigid body is

		r(t)=x(t) i+y(t)j+z(t) k,
									
and if we define the vector ω(t) by

		ω(t)=(Subscript[ω, 1])(t) i+(Subscript[ω, 2])(t) j + (Subscript[ω, 3])(t) k,
						
then we have
										    		     						                            			r'(t)=B(t) r(t)
           =(({{0, -(Subscript[ω, 3])(t), (Subscript[ω, 2])(t)}, {(Subscript[ω, 3])(t), 0, -(Subscript[ω, 1])(t)}, {-(Subscript[ω, 2])(t), (Subscript[ω, 1])(t), 0}}))(({{x(t)}, {y(t)}, {z(t)}}))
           =(({{(Subscript[ω, 2])(t)z(t)-(Subscript[ω, 3])(t)y(t)}, {(Subscript[ω, 3])(t)x(t)-(Subscript[ω, 1 ])(t)z(t)}, {(Subscript[ω, 1])(t)y(t)-(Subscript[ω, 2]) (t)x(t)}}))
           =ω(t) ×r(t).
						
This is a fundamental result: If a particle moves in space so that its position vector at time t is obtained from the position vector at time 0 by a rotation matrix, that is

		r(t)=A(t) r(0),
						
then there is a vector function of time ω(t), called the angular velocity vector, such that the particle's velocity vector r'(t) satisfies the equation

		r'(t)=ω(t) ×r(t).				
						
The vector ω(t) determines an instantaneous axis of rotation : those vectors that are at rest at time t are precisely those whose position at time t lies along the line through the origin determined by ω(t).

Dynamical Equations for Systems of Particles

View our rigid body as consisting of masses (Subscript[m, 1]), (Subscript[m, 2]),...,(Subscript[m, N]) with position vectors (Subscript[r, 1])(t), ...,(Subscript[r, N])(t).  Let the total force on (Subscript[m, i]) be

		(Subscript[F, i])=(Subscript[F, i]^ext)+(∑_(j≠i))(Subscript[F, i j])			
					
where (Subscript[F, i]^ext) is an external force, exerted from outside the system, and (Subscript[F, i j]) is an internal force exerted on mass (Subscript[m, i]) due to the presence of (Subscript[m, j]).  We define the total external force

		(F^ext)=(∑_(i=1)^N)(Subscript[F, i]^ext)
			
and we assume 

		(F^ext)=0.

The angular momentum of the i-th particle  about the origin is 

		(Subscript[L, i])=(Subscript[m, i]) (Subscript[r, i])(t)×(Subscript[r, i])'(t).				
					
Define the quantities:

		L=(∑_(i=1)^N)(Subscript[L, i])     	(total angular momentum)	
	 	N=(∑_(i=1)^N)(Subscript[r, i])(t)×(Subscript[F, i]^ext)    (total torque)		
	 		
We can then derive the dynamical law
					
		((d L)/(d t))=N
.

The Inertia Tensor of a Rigid Body

The rigid body that we will consider is the rectangular solid shown below.

Button["Rectangular Solid",CellPrint[Show[Graphics3D[rotatedBox[IdentityMatrix[3]]],BoxRatios->{1,1,1},AxesLabel->{Style["X",14,"Bold"],Style["Y",14,"Bold"],Style["Z",14,"Bold"]},Axes->True,Ticks->None,PlotRange->{{-10,10},{-10,10},{-10,10}},PlotLabel->Style["Rectangular Solid Body",FontSize->18,FontWeight->"Bold",FontColor->RGBColor[.6,0,0],FontFamily->"Times New Roman"]]
]]

A rigid body is a system of particles satisfying the additional condition that for each mass we have

		r'(t)=ω(t) ×r(t).				
					
The total angular momentum is then

	L=(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])(t)×(Subscript[r, i])'(t)=(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])(t)×(ω(t) ×(Subscript[r, i])(t)).				

Given any list of masses (Subscript[m, 1]), (Subscript[m, 2]),...,(Subscript[m, N]) and any set of N vectors (Subscript[r, 1]), ...,(Subscript[r, N]), we can define a linear transformation

		I:(ℝ^3)->(ℝ^3)
					
by

		I(v)=(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])×(v×(Subscript[r, i])).
					
This is the inertia tensor for the given masses and position vectors.  If the position vectors are changing with time, so will this inertia tensor be a function of time. Therefore, for a moving rigid body, we have the inertia tensor

		(Subscript[I, t])(v)=(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])(t)×(v×(Subscript[r, i])(t)),
					
and the following relationship holds:

		L(t)=(Subscript[I, t])(ω(t))
					
or, as it is often written,

		L=I·ω.
		
So the total angular momentum of a rigid body is expressible in terms of the inertia tensor and the angular velocity vector.

The kinetic energy E of the body is  similarly expressible:

		E=(1/2)ω(t)·(Subscript[I, t])(ω(t)).
		
Because:

		(1/2)ω(t)·(Subscript[I, t])(ω(t))=(1/2)ω(t)·(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])(t)×(ω(t) ×(Subscript[r, i])(t))
                              =(1/2)ω(t)·(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])(t)×(Subscript[r, i])'(t)
                              =(1/2)(∑_(i=1)^N)(Subscript[m, i]) ω(t)·(Subscript[r, i])(t)×(Subscript[r, i])'(t)
                              =(1/2)(∑_(i=1)^N)(Subscript[m, i]) ω(t)×(Subscript[r, i])(t)·(Subscript[r, i])'(t)
                              =(1/2)(∑_(i=1)^N)(Subscript[m, i]) (Subscript[r, i])'(t)·(Subscript[r, i])'(t)
                              =E.

Torque-Free Motion with One Point Fixed

If a rigid-body moves so that one point is fixed and there is no external torque, then

		((dL)/dt)=0,
so that
			
		L=(Subscript[L, 0])=constant   (Conservation of angular momentum)
					
Then

		(d/dt)(I·ω)=0.

Equations of Motion in Body-Coordinates

A body-coordinate system is a right-handed Cartesian coordinate system that moves with the body.

 Let {(Subscript[u, 1])(t),(Subscript[u, 2])(t), (Subscript[u, 3])(t)}  be the unit vectors for a body-coordinate system. 

Components of the inertia tensor are constant when expressed in body-coordinates and the matrix is symmetric and positive definite. 
			
Choose another set of body-axes so that ((Subscript[I, i j])) is diagonal  

		((Subscript[I, i j]))=(({{(Subscript[I, 1]), 0, 0}, {0, (Subscript[I, 2]), 0}, {0, 0, (Subscript[I, 3])}}))    ((Subscript[I, 1])>(Subscript[I, 2])>(Subscript[I, 3]) >0).
				
In terms of components in body coordinates, we have

		 0=(d/dt)((Subscript[I, t])·ω(t))=(d/dt)((∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])(t)(Subscript[u, i])(t))
		  =(∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])'(t)(Subscript[u, i])(t)+(∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])(t)(Subscript[u, i])'(t)
                               =(∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])'(t)(Subscript[u, i])(t)+(∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])(t)ω(t)×(Subscript[u, i])(t).	
				
Now the second sum can be written as	

		ω(t)×(∑_(i=1)^3)(Subscript[I, i ])(Subscript[ω, i])(t)(Subscript[u, i])(t)=((Subscript[ω, 1])(t)(Subscript[u, 1])(t)+(Subscript[ω, 2])(t)(Subscript[u, 2])(t)+(Subscript[ω, 3])(t)(Subscript[u, 3])(t))×((Subscript[I, 1])(Subscript[ω, 1])(t)(Subscript[u, 1])(t)+(Subscript[I, 2])(Subscript[ω, 2])(t)(Subscript[u, 2])(t)+(Subscript[I, 3])(Subscript[ω, 3])(t)(Subscript[u, 3])(t))
                                                      =((Subscript[I, 3])-(Subscript[I, 2]))(Subscript[ω, 2])(Subscript[ω, 3])(Subscript[u, 1])(t)+((Subscript[I, 1])-(Subscript[I, 3]))(Subscript[ω, 1])(Subscript[ω, 3])(Subscript[u, 2])(t)+((Subscript[I, 2])-(Subscript[I, 1]))(Subscript[ω, 1])(Subscript[ω, 2])(Subscript[u, 3])(t).
							
Setting to zero each component of the above equation, we get Euler’s equations of torque-free rigid-body motion:

		(Subscript[I, 1])(Subscript[ω, 1])'+((Subscript[I, 3])-(Subscript[I, 2]))(Subscript[ω, 2])(Subscript[ω, 3])=0   
		(Subscript[I, 2])(Subscript[ω, 2])'+((Subscript[I, 1])-(Subscript[I, 3]))(Subscript[ω, 1])(Subscript[ω, 3])=0         (Euler’s equations)                               
		(Subscript[I, 3])(Subscript[ω, 3])'+((Subscript[I, 2])-(Subscript[I, 1]))(Subscript[ω, 1])(Subscript[ω, 2])=0.                                               

Euler’s equations comprise a system of three ODEs for the three unknowns (Subscript[ω, 1])(t), (Subscript[ω, 2])(t), and (Subscript[ω, 3])(t).  

These functions are components of the angular velocity vector ω(t) with respect to the body-fixed axes.

Equations of Motion in Space Coordinates

To describe the motion in space, we need to know the rotation matrix A(t) that rotates the body from its reference position to its position at time t.  This can be done by solving the system of differential equations consisting of the three Euler equations and the equations

		(Subscript[u, i])'(t)=ω(t)×(Subscript[u, i])(t)             i=1,2,3
		
If we use initial-conditions specifying that (Subscript[u, i])(0)=(Subscript[e, i]) (where (Subscript[e, 1]),(Subscript[e, 2]),(Subscript[e, 3]) are the unit vectors of a spatially-fixed reference coordinate system), then we have

		A(t)=[(Subscript[u, 1])(t) (Subscript[u, 2])(t) (Subscript[u, 3])(t)].

Stationary Rotations

Eigenvectors of ((Subscript[I, i j])) give solutions constant in body-coordinates.  So an eigenvector for the eigenvalue (Subscript[I, 1]) is a constant solution 

		((Subscript[ω, 1])(t), (Subscript[ω, 2])(t),(Subscript[ω, 3])(t))=((Subscript[C, 1]),0,0).

A Stationary Rotation in Space

Button["A Stationary Rotation",CellPrint[
ω0={(Sqrt[((2 )/I1)]),0,0};   
tMax=5000;
 Clear[u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3,u,w,eqList];
eqList={u11'[t]==-w2[t] u13[t]+w3[t] u12[t],u21'[t]==-w2[t] u23[t]+w3[t] u22[t],u31'[t]==-w2[t] u33[t]+w3[t] u32[t],u12'[t]==w1[t] u13[t] -w3[t] u11[t],u22'[t]==w1[t] u23[t] -w3[t] u21[t],u32'[t]==w1[t] u33[t] -w3[t] u31[t],u13'[t]==w2[t] u11[t]-w1[t] u12[t],u23'[t]==w2[t] u21[t]-w1[t] u22[t],u33'[t]==w2[t] u31[t]-w1[t] u32[t],I1 w1'[t]==(I2-I3)w2[t] w3[t],I2 w2'[t]==(I3-I1)w1[t] w3[t],I3 w3'[t]==(I1-I2)w2[t] w1[t],u11[0]==1,u12[0]==0,u13[0]==0,u21[0]==0,u22[0]==1,u23[0]==0,u31[0]==0,u32[0]==0,u33[0]==1,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]};
{u11[t_],u21[t_],u31[t_],u12[t_],u22[t_],u32[t_],u13[t_],u23[t_],u33[t_],w1[t_],w2[t_],w3[t_]}={u11[t],u21[t],u31[t],u12[t],u22[t],u32[t],u13[t],u23[t],u33[t],w1[t],w2[t],w3[t]}/.NDSolve[eqList,{u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3},{t,0,tMax}][[1]];

u[t_]:={{u11[t],u12[t],u13[t]},{u21[t],u22[t],u23[t]},{u31[t],u32[t],u33[t]}};
w[t_]:={w1[t],w2[t],w3[t]};


v[t_]:=10w[t]/Norm[w[t]] ;
angMom[t_]:=Line[{{0,0,0},u[t].IT.w[t]}];
angVel[t_]:=Line[{{0,0,0},u[t].v[t]}];
Manipulate[Show[Graphics3D[{Blue,Thickness[.02],angVel[t]}],Graphics3D[{Red,Thickness[.02],angMom[t]}],Graphics3D[rotatedBox[u[t]]],BoxRatios->{1,1,1},PlotRange->{{-10,10},{-10,10},{-10,10 }},AxesLabel->{Style["X",14,"Bold"],Style["Y",14,"Bold"],Style["Z",14,"Bold"]},Axes->True,Ticks->None,ViewPoint->{1.300, -2.400, 2.000}],{t,0,tMax,Appearance->"Open"}]
]]

The Energy and Momentum Ellipsoids

Any solution of the Euler equations satisfies conservation of kinetic energy and angular momentum.  That is, the quantities

		(Subscript[I, 1])((Subscript[ω, 1])^2)+(Subscript[I, 2])((Subscript[ω, 2])^2)+(Subscript[I, 3])((Subscript[ω, 3])^2)
		
		((Subscript[I, 1])^2)((Subscript[ω, 1])^2)+((Subscript[I, 2])^2)((Subscript[ω, 2])^2)+((Subscript[I, 3])^2)((Subscript[ω, 3])^2)
				
are constant along any solution.  

Pick an initial point satisfying

		(Subscript[ℰ, T]): (Subscript[I, 1])(x^2)+(Subscript[I, 2])(y^2)+(Subscript[I, 3])(z^2)=2(Subscript[E, 0])     (KE (Subscript[E, 0]))
				
Compute the momentum quantity for the chosen point and then our solution must lie in the intersection of the ellipsoid (Subscript[ℰ, T])with the momentum ellipsoid

		(Subscript[ℰ, M]): ((Subscript[I, 1])^2)((Subscript[ω, 1])^2)+((Subscript[I, 2])^2)((Subscript[ω, 2])^2)+((Subscript[I, 3])^2)((Subscript[ω, 3])^2)=((Subscript[L, 0])^2)

Energy and Momentum Ellipsoids

I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];
ϵ=0.001`;
ω0={ϵ,(Sqrt[((2-I1 (ϵ^2))/I2)]),0};
b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];

ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100] ;ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100] ;
Show[ellipEnergy,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]
Show[ellipMomentum,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]
Show[ellipEnergy,ellipMomentum,PlotLabel->Style["Intersecting Ellipsoids",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]

Solution Initially near the Long Axis

Energy and Momentum Ellipsoids

Button["Intersecting Ellipsoids",CellPrint[
ϵ=0.0008`;
ω0={0,ϵ,(Sqrt[((2-I2 (ϵ^2))/I3)])};
b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];

ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->50]; 

ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->50] ;

Show[ellipEnergy,ellipMomentum,PlotLabel->Style["Intersecting Ellipsoids",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]
]]

Motion in the Body Reference Frame

Button["Motion in the Body Frame",CellPrint[
Clear[w1,w2,w3];
ϵ=0.0008`;
ω0={0,ϵ,(Sqrt[((2-I2 (ϵ^2))/I3)])};
b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];

tMax=5000;
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->50]; 

ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->50] ;

{w1[t_],w2[t_],w3[t_]}={w1[t],w2[t],w3[t]}/.NDSolve[{eigs[[1]] w1'[t]+(eigs[[3]]-eigs[[2]]) w2[t] w3[t]==0,eigs[[2]] w2'[t]+(eigs[[1]]-eigs[[3]]) w1[t] w3[t]==0,eigs[[3]] w3'[t]+(eigs[[2]]-eigs[[1]]) w2[t] w1[t]==0,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]},{w1,w2,w3},{t,0,tMax}][[1]];
Manipulate[Show[ellipEnergy,ellipMomentum,ParametricPlot3D[1.02{ w1[t], w2[t],w3[t] },{t,0,tMax},BoxRatios->{1,1,1},PlotStyle->Directive[Thickness[.004],Yellow]],Graphics3D[{Yellow,PointSize[0.02],Point[{1.05` w1[τ],1.05` w2[τ],1.05` w3[τ]}]}],PlotLabel->Style["Body Trajectory on the Intersection",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],AxesLabel->{"X","Y","Z"},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}}],{τ,0,tMax,100,Appearance->"Open"}]
]]

Motion in the Spatial Reference Frame

Button["Motion in Space",CellPrint[
Clear[u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3,u,w,eqList];
eqList={u11'[t]==-w2[t] u13[t]+w3[t] u12[t],u21'[t]==-w2[t] u23[t]+w3[t] u22[t],u31'[t]==-w2[t] u33[t]+w3[t] u32[t],u12'[t]==w1[t] u13[t] -w3[t] u11[t],u22'[t]==w1[t] u23[t] -w3[t] u21[t],u32'[t]==w1[t] u33[t] -w3[t] u31[t],u13'[t]==w2[t] u11[t]-w1[t] u12[t],u23'[t]==w2[t] u21[t]-w1[t] u22[t],u33'[t]==w2[t] u31[t]-w1[t] u32[t],I1 w1'[t]==(I2-I3)w2[t] w3[t],I2 w2'[t]==(I3-I1)w1[t] w3[t],I3 w3'[t]==(I1-I2)w2[t] w1[t],u11[0]==1,u12[0]==0,u13[0]==0,u21[0]==0,u22[0]==1,u23[0]==0,u31[0]==0,u32[0]==0,u33[0]==1,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]};
{u11[t_],u21[t_],u31[t_],u12[t_],u22[t_],u32[t_],u13[t_],u23[t_],u33[t_],w1[t_],w2[t_],w3[t_]}={u11[t],u21[t],u31[t],u12[t],u22[t],u32[t],u13[t],u23[t],u33[t],w1[t],w2[t],w3[t]}/.NDSolve[eqList,{u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3},{t,0,tMax}][[1]];

u[t_]:={{u11[t],u12[t],u13[t]},{u21[t],u22[t],u23[t]},{u31[t],u32[t],u33[t]}};
w[t_]:={w1[t],w2[t],w3[t]};


v[t_]:=10w[t]/Norm[w[t]] ;
angMom[t_]:=Line[{{0,0,0},u[t].IT.w[t]}];
angVel[t_]:=Line[{{0,0,0},u[t].v[t]}];
Manipulate[
Show[Graphics3D[{Yellow,Thickness[.02],angVel[τ]}],Graphics3D[{Orange,Thickness[.02],angMom[τ]}],Graphics3D[rotatedBox[u[τ]]],BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"},PlotRange->{{-12,12},{-12,12},{-12,12}}],{τ,0,tMax,Appearance->"Open"}]
]]

Solution Initially near the Short Axis

Energy and Momentum Ellipsoids

Button["Intersecting Ellipsoids",CellPrint[
I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];
ϵ=0.001`;
ω0={(Sqrt[((2-I2 (ϵ^2))/I1)]),ϵ,0};b1=(Sqrt[(2/I1)]);

b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,PlotPoints->100,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]];
ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,PlotPoints->100,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]];
Show[ellipMomentum,ellipEnergy,PlotLabel->Style["Intersecting Ellipsoids",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]
]]

Motion in the Body Reference Frame

Button["Motion in the Body Frame",CellPrint[
Clear[w1,w2,w3];
I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];
ϵ=0.001`;
ω0={(Sqrt[((2-I2 (ϵ^2))/I1)]),ϵ,0};b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];

tMax=5000;
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,PlotPoints->100,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]];
ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotPoints->100,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]];
{w1[t_],w2[t_],w3[t_]}={w1[t],w2[t],w3[t]}/.NDSolve[{eigs[[1]] w1'[t]+(eigs[[3]]-eigs[[2]]) w2[t] w3[t]==0,eigs[[2]] w2'[t]+(eigs[[1]]-eigs[[3]]) w1[t] w3[t]==0,eigs[[3]] w3'[t]+(eigs[[2]]-eigs[[1]]) w2[t] w1[t]==0,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]},{w1,w2,w3},{t,0,tMax}][[1]];
Manipulate[Show[ellipEnergy,ellipMomentum,ParametricPlot3D[{1.01` w1[t],1.01` w2[t],1.01` w3[t]},{t,0,tMax},BoxRatios->{1,1,1},PlotStyle->Directive[Thickness[.003],Yellow]],Graphics3D[{Yellow,PointSize[0.02],Point[{1.05` w1[τ],1.05` w2[τ],1.05` w3[τ]}]}],AxesLabel->{"X","Y","Z"},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},PlotLabel->Style["Body trajectory on the Intersection",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]],{τ,.1,tMax,Appearance->"Open"}]
]]

Motion in the Spatial Reference Frame

Button["Motion in Space",CellPrint[
Clear[u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3,u,w,eqList];
eqList={u11'[t]==-w2[t] u13[t]+w3[t] u12[t],u21'[t]==-w2[t] u23[t]+w3[t] u22[t],u31'[t]==-w2[t] u33[t]+w3[t] u32[t],u12'[t]==w1[t] u13[t] -w3[t] u11[t],u22'[t]==w1[t] u23[t] -w3[t] u21[t],u32'[t]==w1[t] u33[t] -w3[t] u31[t],u13'[t]==w2[t] u11[t]-w1[t] u12[t],u23'[t]==w2[t] u21[t]-w1[t] u22[t],u33'[t]==w2[t] u31[t]-w1[t] u32[t],I1 w1'[t]==(I2-I3)w2[t] w3[t],I2 w2'[t]==(I3-I1)w1[t] w3[t],I3 w3'[t]==(I1-I2)w2[t] w1[t],u11[0]==1,u12[0]==0,u13[0]==0,u21[0]==0,u22[0]==1,u23[0]==0,u31[0]==0,u32[0]==0,u33[0]==1,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]};
{u11[t_],u21[t_],u31[t_],u12[t_],u22[t_],u32[t_],u13[t_],u23[t_],u33[t_],w1[t_],w2[t_],w3[t_]}={u11[t],u21[t],u31[t],u12[t],u22[t],u32[t],u13[t],u23[t],u33[t],w1[t],w2[t],w3[t]}/.NDSolve[eqList,{u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3},{t,0,tMax}][[1]];

u[t_]:={{u11[t],u12[t],u13[t]},{u21[t],u22[t],u23[t]},{u31[t],u32[t],u33[t]}};
w[t_]:={w1[t],w2[t],w3[t]};


v[t_]:=10w[t]/Norm[w[t]] ;
angMom[t_]:=Line[{{0,0,0},u[t].IT.w[t]}];
angVel[t_]:=Line[{{0,0,0},u[t].v[t]}];
Manipulate[Show[Graphics3D[{Yellow,Thickness[.02],angVel[τ]}],Graphics3D[{Orange,Thickness[.02],angMom[τ]}],Graphics3D[rotatedBox[u[τ]]],BoxRatios->{1,1,1},PlotRange->{{-10,10},{-10,10},{-10,10 }}],{τ,.1,tMax,Appearance->"Open"}]
]]

Solution Starting Near the Intermediate Axis

Energy and Momentum Ellipsoids

Button["Intersecting Ellipsoids",CellPrint[
I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];
ϵ=0.0002`;
ω0={ϵ,(Sqrt[((2-I1 (ϵ^2))/I2)]),0};b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);

rad=1.1` Max[b1,b2,b3,β1,β2,β3];
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100];
ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100];
Show[ellipMomentum,ellipEnergy,PlotLabel->Style["Intersecting Ellipsoids",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]]
]]

Motion in the Body Reference Frame

Button["Motion in the Body Frame",CellPrint[
Clear[w1,w2,w3,α,β];
I1=eigs[[1]];
I2=eigs[[2]];
I3=eigs[[3]];
ϵ=0.0002`;
ω0={ϵ,(Sqrt[((2-I1 (ϵ^2))/I2)]),0};
b1=(Sqrt[(2/I1)]);
b2=(Sqrt[(2/I2)]);
b3=(Sqrt[(2/I3)]);
L0=Norm[eigs ω0];
β1=(L0/I1);
β2=(L0/I2);
β3=(L0/I3);
rad=1.1` Max[b1,b2,b3,β1,β2,β3];
tMax=10000;

sol=Solve[{I1 (α^2)+I3 (β^2)==2,(I1^2) (α^2)+(I3^2) (β^2)==(L0^2)},{α,β}][[4]];

{α,β}={α,β}/.sol;
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100];
ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->100];
{w1[t_],w2[t_],w3[t_]}={w1[t],w2[t],w3[t]}/.NDSolve[{eigs[[1]] w1'[t]+(eigs[[3]]-eigs[[2]]) w2[t] w3[t]==0,eigs[[2]] w2'[t]+(eigs[[1]]-eigs[[3]]) w1[t] w3[t]==0,eigs[[3]] w3'[t]+(eigs[[2]]-eigs[[1]]) w2[t] w1[t]==0,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]},{w1,w2,w3},{t,0,tMax}][[1]];
Manipulate[Show[ellipEnergy,ellipMomentum,ParametricPlot3D[{1.01 w1[t],1.01 w2[t],1.01 w3[t]},{t,0,tMax},BoxRatios->{1,1,1},PlotStyle->Directive[Thickness[.002],Yellow]],Graphics3D[{Yellow,PointSize[0.02],Point[{1.04` w1[τ],1.04` w2[τ],1.04` w3[τ]}]}],PlotLabel->Style["Body trajectory on the Intersection",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],AxesLabel->{"X","Y","Z"},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},ViewPoint->{2.962`,-1.222`,1.104`}],{τ,0,tMax,Appearance->"Open"}]
]]

Motion in the Spatial Reference Frame

Button["Motion in Space",CellPrint[
Clear[u11, u21, u31, u12, u22, u32, u13, u23, u33, w1, w2, w3, u, w, eqList];
eqList = {u11'[t] == -w2[t] u13[t] + w3[t] u12[t], u21'[t] == -w2[t] u23[t] + w3[t] u22[t], u31'[t] == -w2[t] u33[t] + w3[t] u32[t], u12'[t] == w1[t] u13[t] - w3[t] u11[t], u22'[t] == w1[t] u23[t] - w3[t] u21[t], u32'[t] == w1[t] u33[t] - w3[t] u31[t], u13'[t] == w2[t] u11[t] - w1[t] u12[t], u23'[t] == w2[t] u21[t] - w1[t] u22[t], u33'[t] == w2[t] u31[t] - w1[t] u32[t], I1 w1'[t] == (I2 - I3) w2[t] w3[t], I2 w2'[t] == (I3 - I1) w1[t] w3[t], I3 w3'[t] == (I1 - I2) w2[t] w1[t], u11[0] == 1, u12[0] == 0, u13[0] == 0, u21[0] == 0, u22[0] == 1, u23[0] == 0, u31[0] == 0, u32[0] == 0, u33[0] == 1, w1[0] == ω0[[1]], w2[0] == ω0[[2]], w3[0] == ω0[[3]]};
{u11[t_], u21[t_], u31[t_], u12[t_], u22[t_], u32[t_], u13[t_], u23[t_], u33[t_], w1[t_], w2[t_], w3[t_]} = {u11[t], u21[t], u31[t], u12[t], u22[t], u32[t], u13[t], u23[t], u33[t], w1[t], w2[t], w3[t]} /. NDSolve[eqList, {u11, u21, u31, u12, u22, u32, u13, u23, u33, w1, w2, w3}, {t, 0, tMax}][[1]];

u[t_] := {{u11[t], u12[t], u13[t]}, {u21[t], u22[t], u23[t]}, {u31[t], u32[t], u33[t]}};
w[t_] := {w1[t], w2[t], w3[t]};


v[t_] := 10 w[t]/Norm[w[t]] ;
angMom[t_] := Line[{{0, 0, 0}, u[t].IT.w[t]}];
angVel[t_] := Line[{{0, 0, 0}, u[t].v[t]}];
Manipulate[Show[Graphics3D[{Yellow, Thickness[.02], angVel[τ]}], Graphics3D[{Orange, Thickness[.02], angMom[τ]}], Graphics3D[rotatedBox[u[τ]]], BoxRatios -> {1, 1, 1}, PlotRange -> {{-10, 10}, {-10, 10}, {-10, 10 }}, AxesLabel -> {"X", "Y", "Z"}, ViewPoint -> {3.106, -1.044, 0.845}] ,{τ,0,tMax,Appearance->"Open"}]
]]

A Heteroclinic Solution

If we choose the initial value ω(0) on the intermediate body axis, then the intersection of the ellipsoids will contain that point.  There are then four solutions of the Euler equations whose trajectories join ω(0) to -ω(0).  If we choose ω(0) to be on one of these intersection curves and very near the intermediate axis we get a solution the approaches one critical point as t->-∞  and the one on the opposite side of the ellipsoids as t->∞ .   Such a solution is called heteroclinic.

Energy and Momentum Ellipsoids

Button["Intersecting Ellipsoids",CellPrint[
ω0={0,(Sqrt[(2/I2)]),0};b1=(Sqrt[(2/I1)]);b2=(Sqrt[(2/I2)]); b3=(Sqrt[(2/I3)]);L0=Norm[eigs ω0];β1=(L0/I1); β2=(L0/I2); β3=(L0/I3); rad=1.1 Max[b1,b2,b3,β1,β2,β3];
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->200];

ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->200];

Show[ellipMomentum,ellipEnergy,PlotLabel->Style["Intersecting Ellipsoids",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],ViewPoint->{0.785`,-3.123`,1.04`}]
]]

Motion in the Body Reference Frame

Take (Subscript[ω, 0]) on the intersection; Then solution approaches stationary rotation about the intermediate axis.

Button["Motion in the Body Frame",CellPrint[
Clear[w1,w2,w3,α,β];

ϵ=0.0001;
γ=(Sqrt[((2-ϵ)/I2)]);
ω0;
{α,β}={α,β}/.Solve[{I1 (α^2)+I2 (γ^2)+I3 (β^2)==2,(I1^2) (α^2)+(I2^2) (γ^2)+(I3^2) (β^2)==(L0^2)},{α,β}][[4]];
tMax=10000;
ω0={α,γ,β};
ellipEnergy=ParametricPlot3D[{b1 Cos[θ] Sin[ϕ],b2 Sin[θ] Sin[ϕ],b3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Blue,Mesh->None,PlotLabel->Style["Energy Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->200];ellipMomentum=ParametricPlot3D[{β1 Cos[θ] Sin[ϕ],β2 Sin[θ] Sin[ϕ],β3 Cos[ϕ]},{θ,0,2 π},{ϕ,0,π},BoxRatios->{1,1,1},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},Ticks->None,AxesLabel->{"X","Y","Z"},PlotStyle->Red,Mesh->None,PlotLabel->Style["Momentum Ellipsoid",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"],PlotPoints->200];{w1[t_],w2[t_],w3[t_]}={w1[t],w2[t],w3[t]}/.NDSolve[{eigs[[1]] w1'[t]+(eigs[[3]]-eigs[[2]]) w2[t] w3[t]==0,eigs[[2]] w2'[t]+(eigs[[1]]-eigs[[3]]) w1[t] w3[t]==0,eigs[[3]] w3'[t]+(eigs[[2]]-eigs[[1]]) w2[t] w1[t]==0,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]},{w1,w2,w3},{t,0,tMax}][[1]];
Manipulate[Show[ellipEnergy,ellipMomentum,ParametricPlot3D[{1.01` w1[t],1.01` w2[t],1.01` w3[t]},{t,0,tMax},BoxRatios->{1,1,1},PlotStyle->Directive[Thickness[.002],Yellow]],Graphics3D[{Yellow,PointSize[0.02`],Point[{1.04 w1[τ],1.04w2[τ],1.04w3[τ]}]}],AxesLabel->{"X","Y","Z"},PlotRange->{{-rad,rad},{-rad,rad},{-rad,rad}},ViewPoint->{3.052`,-1.026`,1.04`},PlotLabel->Style["Body trajectory on intersection",FontSize->18,FontWeight->"Bold",FontColor->Red,FontFamily->"Times New Roman"]],{τ,0,tMax,.1,Appearance->"Open"}]
]]

Motion in the Spatial Reference Frame

Button["Motion in Space",CellPrint[
Clear[u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3,u,w,eqList];
eqList={u11'[t]==-w2[t] u13[t]+w3[t] u12[t],u21'[t]==-w2[t] u23[t]+w3[t] u22[t],u31'[t]==-w2[t] u33[t]+w3[t] u32[t],u12'[t]==w1[t] u13[t] -w3[t] u11[t],u22'[t]==w1[t] u23[t] -w3[t] u21[t],u32'[t]==w1[t] u33[t] -w3[t] u31[t],u13'[t]==w2[t] u11[t]-w1[t] u12[t],u23'[t]==w2[t] u21[t]-w1[t] u22[t],u33'[t]==w2[t] u31[t]-w1[t] u32[t],I1 w1'[t]==(I2-I3)w2[t] w3[t],I2 w2'[t]==(I3-I1)w1[t] w3[t],I3 w3'[t]==(I1-I2)w2[t] w1[t],u11[0]==1,u12[0]==0,u13[0]==0,u21[0]==0,u22[0]==1,u23[0]==0,u31[0]==0,u32[0]==0,u33[0]==1,w1[0]==ω0[[1]],w2[0]==ω0[[2]],w3[0]==ω0[[3]]};
{u11[t_],u21[t_],u31[t_],u12[t_],u22[t_],u32[t_],u13[t_],u23[t_],u33[t_],w1[t_],w2[t_],w3[t_]}={u11[t],u21[t],u31[t],u12[t],u22[t],u32[t],u13[t],u23[t],u33[t],w1[t],w2[t],w3[t]}/.NDSolve[eqList,{u11,u21,u31,u12,u22,u32,u13,u23,u33,w1,w2,w3},{t,0,tMax}][[1]];

u[t_]:={{u11[t],u12[t],u13[t]},{u21[t],u22[t],u23[t]},{u31[t],u32[t],u33[t]}};
w[t_]:={w1[t],w2[t],w3[t]};

v[t_]:=10w[t]/Norm[w[t]] ;
angMom[t_]:=Line[{{0,0,0},u[t].IT.w[t]}];
angVel[t_]:=Line[{{0,0,0},u[t].v[t]}];
Manipulate[Show[Graphics3D[{Yellow,Thickness[.02],angVel[τ]}],Graphics3D[{Orange,Thickness[.02],angMom[τ]}],Graphics3D[rotatedBox[u[τ]]],BoxRatios->{1,1,1},AxesLabel->{"X","Y","Z"},PlotRange->{{-10,10},{-10,10},{-10,10 }},ViewPoint->{3.130, -0.839, 0.975}],{τ,0,tMax,Appearance->"Open"}]
]]