WOLFRAM|DEMONSTRATIONS PROJECT

Free Rotation of an Asymmetric Top

a=4;b=2;c=1;​​q1=0.01;q2=3;q3=1;​​m=1000;​​T=5;​​kk=0.0005;
I1=m(b*b+c*c)/12;​​I2=m(c*c+a*a)/12;​​I3=m(a*a+b*b)/12;
EE=0.5(I1*q1*q1+I2*q2*q2+I3*q3*q3);​​MM=Sqrt[I1*I1*q1*q1+I2*I2*q2*q2+I3*I3*q3*q3];
range=kkMax[Sqrt[2EE*I1],Sqrt[2EE*I2],N@Sqrt[2EE*I3]];
sgnQ1=Sign[q1];​​sgnQ3=Sign[q3];
Ω[t_]:={Ω1[t],Ω2[t],Ω3[t]};​​AngularMomentum[t_]:=rotationMatrix[t].{I1Ω1[t],I2Ω2[t],I3Ω3[t]};
If2EE*I2<
2
MM
,​​Block{ch,ss,k2,AA,Δt},​​ss=q2
I2(I3-I2)
2EEI3-MM*MM
;​​k2=
(I2-I1)(2EEI3-MM*MM)
(I3-I2)(MM*MM-2EEI1)
;​​ch=
I1I2I3
(I3-I2)(MM*MM-2EEI1)
;​​AA=-NIntegrate
1
(1-s*s)(1-k2*s*s)
,{s,0,ss};​​Δt=AAch;​​Period=4chNIntegrate
1
(1-s*s)(1-k2*s*s)
,{s,0,1};​​Ω1[t_]=sgnQ1
2EEI3-MM*MM
I1(I3-I1)
JacobiCN[(t-Δt)/ch,k2];​​Ω2[t_]=
2EEI3-MM*MM
I2(I3-I2)
JacobiSN[(t-Δt)/ch,k2];​​Ω3[t_]=sgnQ3
MM*MM-2EEI1
I3(I3-I1)
JacobiDN[(t-Δt)/ch,k2]​​;,​​​​Block{ch,ss,k2,AA,Δt},​​ss=q2
I2(I1-I2)
2EEI1-MM*MM
;​​k2=
(I2-I3)(2EEI1-MM*MM)
(I1-I2)(MM*MM-2EEI3)
;​​ch=
I1I2I3
(I1-I2)(MM*MM-2EEI3)
;​​AA=-NIntegrate
1
(1-s*s)(1-k2*s*s)
,{s,0,ss};​​Δt=AAch;​​Period=4chNIntegrate
1
(1-s*s)(1-k2*s*s)
,{s,0,1};​​Ω3[t_]=sgnQ3
2EEI1-MM*MM
I3(I1-I3)
JacobiCN[(t-Δt)/ch,k2];​​Ω2[t_]=
2EEI1-MM*MM
I2(I1-I2)
JacobiSN[(t-Δt)/ch,k2];Ω1[t_]=sgnQ1
MM*MM-2EEI3
I1(I1-I3)
JacobiDN[(t-Δt)/ch,k2]​​;;
rotationMatrix[t_]=
a11[t]
a12[t]
a13[t]
a21[t]
a22[t]
a23[t]
a31[t]
a32[t]
a33[t]
;
solution=NDSolveFlatten@Thread/@Thread
∂
t
rotationMatrix[t]rotationMatrix[t].
0
-Ω3[t]
Ω2[t]
Ω3[t]
0
-Ω1[t]
-Ω2[t]
Ω1[t]
0
,Flatten@(Thread/@Thread[rotationMatrix[0]IdentityMatrix[3]]),{a11,a12,a13,a21,a22,a23,a31,a32,a33}​​,{t,0,T},MaxStepsInfinity;
aM=First[{a11,a12,a13,a21,a22,a23,a31,a32,a33}/.solution];
rotationMatrix[t_]:={{aM[[1]][t],aM[[2]][t],aM[[3]][t]},{aM[[4]][t],aM[[5]][t],aM[[6]][t]},{aM[[7]][t],aM[[8]][t],aM[[9]][t]}};
points1=Table[{t,Ω1[t]},{t,0,T,0.1}];​​points2=Table[{t,Ω2[t]},{t,0,T,0.1}];​​points3=Table[{t,Ω3[t]},{t,0,T,0.1}];​​ΩI1=Interpolation[points1];​​ΩI2=Interpolation[points2];​​ΩI3=Interpolation[points3];
Brick[a_,b_,c_]:={Green,GraphicsComplex[{{-a,-b,-c},{a,-b,-c},{a,b,-c},{-a,b,-c},{-a,-b,c},{a,-b,c},{a,b,c},{-a,b,c}},​​Polygon[{{1,4,3,2},{3,4,8,7},{1,5,8,4},{2,3,7,6},{5,6,7,8},{1,2,6,5}}],VertexColors{Blue,Blue,Blue,Blue,Green,Green,Green,Green}]};
Arrow3D[{x0_,y0_,z0_},{x1_,y1_,z1_},​​{HeadSizeH_,HeadSizeW_,ShaftSize_,color_}]:=​​​​Module[{Vertex1,Vertex2,Vertex3,Vertex4,C,D1,D2,D3,D4,k1,k2,ch,p},​​​​k1=HeadSizeW;​​k2=HeadSizeH*Norm[{x1-x0,y1-y0,z1-z0}];​​p=1;​​ch=0;​​​​If[z1z0,If[y1y0,If[x1x0,p=2,ch=1],ch=2],ch=3];​​​​C=k1{x1-x0,y1-y0,z1-z0}+{x0,y0,z0};​​​​D1={0,0,0};​​D2={0,0,0};​​D3={0,0,0};​​D4={0,0,0};​​​​If[ch3,​​D1[[1]]=1;​​D1[[2]]=0;​​D1[[3]]=((-x0+x1)(1-x0-k1(-x0+x1))+(-y0+y1)(-y0-k1(-y0+y1)))/(z0-z1)+z0+k1(-z0+z1);​​];​​​​If[ch2,​​D1[[1]]=1;​​D1[[3]]=0;​​D1[[2]]=((-x0+x1)(1-x0-k1(-x0+x1))+(-z0+z1)(-z0-k1(-z0+z1)))/(y0-y1)+y0+k1(-y0+y1);​​];​​​​If[ch1,​​D1[[2]]=1;​​D1[[3]]=0;​​D1[[1]]=((-y0+y1)(1-y0-k1(-y0+y1))+(-z0+z1)(-z0-k1(-z0+z1)))/(x0-x1)+x0+k1(-x0+x1);​​];​​​​D1=D1-C;​​D2=Cross[{x1-x0,y1-y0,z1-z0},D1];​​D3=k2Normalize[D1];​​D4=k2Normalize[D2];​​D1=-D3;​​D2=-D4;​​D3=D3+C;​​D4=D4+C;​​D1=D1+C;​​D2=D2+C;​​Vertex1=D3;​​Vertex2=D4;​​Vertex3=D1;​​Vertex4=D2;​​​​{​​{Thickness[ShaftSize],color,​​Line[{{x0,y0,z0},{x1,y1,z1}}],​​GraphicsComplex[{{x1,y1,z1},Vertex1,Vertex2,Vertex3,Vertex4},Polygon[{{1,2,3},{1,3,4},{1,4,5},{1,5,2}}]]},​​​​{color,Point[{x0,y0,z0}]}​​}[[p]]​​​​];