```
Manipulate[
(*Nasser M.ABbasi,version July 3,2014*)
tick;
Module[{t, phi, theta, zeta, L2, thetas, phis, zetas, thetaW, phiW, zetaW, graph, g, sol, eq1, eq2, eq3, Iz, I0, x, mass1, mass2,
cg, coneOff = 10, distToGroud, ic, omegaVector, bodyToInertia, yBody, zBody, L0body, L2body, cgbody, PE, KE, initialSpin, N0,
beta, theta2, n, tmp, omegaVectorBody, inertiaToBody, stripLinesSideOfRotor, stripLinesTopOfRotor, plots},

initialSpin = phiD0  (2*Pi);
mass1 = density*Pi diskr^2*diskThick;(*mass of rotor*)
mass2 = density*Pi (diskr/10)^2*L0;(*mass of bar*)
g = 9.81;
cg = ( mass2*L0/2 + mass1*(L0 + diskThick/2))/(mass1 + mass2);
I0 = 1/12 mass1 (3*diskr^2 + diskThick^2) + mass1* (L0 + diskThick/2)^2;(*parallel axes*)
I0 += 1/12 mass2 (3*(diskr/10)^2 + (L0)^2) + mass2*(L0/2)^2;
Iz = mass1*diskr^2/2;
Iz += mass2*(diskr/10)^2/2;

eq1 = (mass1 + mass2)*g*(cg) Sin[theta[t]] ==
I0 theta''[t] + Iz (phi'[t] + zeta'[t] Cos[theta[t]]) zeta'[t] Sin[theta[t]] - I0 (zeta'[t])^2 Sin[theta[t]] Cos[theta[t]];
(*changed sign*)
eq2 = 0 == I0 (zeta''[t] Sin[theta[t]] + 2 zeta'[t] theta'[t] Cos[theta[t]]) - Iz theta'[t] (phi'[t] + zeta'[t] Cos[theta[t]]);
eq3 = 0 == Iz (phi''[t] + zeta''[t] Cos[theta[t]] - zeta'[t] theta'[t] Sin[theta[t]]);

ic = {theta[0] == currentTheta, theta'[0] == currentThetaD, zeta[0] == currentZeta, zeta'[0] == currentZetaD,
phi[0] == currentPhi, phi'[0] == currentPhiD};
sol = First@NDSolve[Flatten@{eq1, eq2, eq3, ic}, {theta, zeta, phi, theta', zeta', phi'}, {t, 0, delT},
Method -> "StiffnessSwitching"];

thetas = (theta /. sol)[0];
zetas = (zeta /. sol)[0];
phis = (phi /. sol)[0];
thetaW = (theta' /. sol)[0];
zetaW = (zeta' /. sol)[0];
phiW = (phi' /. sol)[0];

n = (phiW + zetaW Cos[thetas]);
N0 = (Iz/I0) n ;
beta = 2 (mass1 + mass2)*g*cg/I0;
tmp = N0^2/(2 beta);
theta2 = ArcCos[tmp - Sqrt[1 - 2 tmp Cos[theta0 Degree] + tmp^2]];(*maximum nutation*)

(*inertiaToBody={{1,0,0},{0,Cos[thetas],Sin[thetas]},{0,-Sin[thetas],Cos[thetas]}}.{{Cos[zetas],Sin[zetas],0},{-Sin[zetas],Cos[
zetas],0},{0,0,1}};*)
(*bodyToInertia=Inverse[inertiaToBody];*)

(*this transformation for only Euler angles zeta and theta (first two) *)
bodyToInertia = {{Cos[zetas], -Cos[thetas] Sin[zetas], Sin[thetas] Sin[zetas]}, {Sin[zetas],
Cos[thetas] Cos[zetas], -Cos[zetas] Sin[thetas]}, {0, Sin[thetas], Cos[thetas]}};
L2 = L0 + diskThick;
(*
inertiaToBody={{Cos[phis]Cos[zetas]-Sin[phis]Cos[thetas] Sin[zetas],Cos[phis]Sin[zetas]+Sin[phis]Cos[thetas]Cos[zetas],Sin[
phis]Sin[thetas]}
,{-Sin[phis]Cos[zetas]-Cos[phis]Cos[thetas] Sin[zetas],-Sin[phis]Sin[zetas]+Cos[phis]Cos[thetas]Cos[zetas],Cos[phis]Sin[thetas]}
,{Sin[thetas]Sin[zetas],-Sin[thetas]Cos[zetas],Cos[thetas]}
};
bodyToInertia={{Cos[phis]Cos[zetas]-Sin[phis]Cos[thetas] Sin[zetas],-Sin[phis]Cos[zetas]-Sin[zetas]Cos[thetas]Cos[phis],Sin[
thetas]Sin[zetas]}
,{Cos[phis]Sin[zetas]+Sin[phis]Cos[thetas] Cos[zetas],-Sin[phis]Sin[zetas]+Cos[phis]Cos[thetas]Cos[zetas],-Sin[thetas]Sin[zetas]}
,{Sin[thetas]Sin[phis],Sin[thetas]Cos[phis],Cos[thetas]}};
*)

L0body = {0, 0, L0};
L2body = {0, 0, L2};
cgbody = {0, 0, cg};
distToGroud = L0 Cos[thetas] - diskr Sin[thetas];

(*case 2*)
(*omegaVector=bodyToInertia.{thetaW Cos[phis]+zetaW*Sin[thetas]Sin[phis],-thetaW Sin[phis]+zetaW Sin[thetas]Cos[phis], zetaW Cos[
thetas]+phiW};*)

(*case 1*)
(*omegaVector=bodyToInertia.{thetaW ,zetaW Sin[thetas]Cos[phis], zetaW Cos[thetas]+phiW};*)
omegaVector = {thetaW , zetaW Sin[thetas], n};
omegaVectorBody = bodyToInertia.omegaVector;
stripLinesSideOfRotor = {bodyToInertia.{diskr Cos[#], diskr Sin[#], L2}, bodyToInertia.{diskr Cos[#], diskr Sin[#], L0}} & /@
Range[0, 2 Pi, Pi/4];
stripLinesTopOfRotor = {bodyToInertia.{diskr Cos[#], diskr Sin[#], L2}, bodyToInertia.{0, 0, L2}} & /@ Range[0, 2 Pi, Pi/4];

graph = Graphics3D[
{
Rotate[GraphicsGroup[
{
{Lighting -> "Neutral", FaceForm[LightGray], Opacity[op],
Cylinder[{{bodyToInertia.L0body, bodyToInertia.L2body}}, .98 diskr]},
{Blue, Line[stripLinesSideOfRotor]},
{Blue, Line[stripLinesTopOfRotor]},

(*
If[showTexture,
First@ParametricPlot3D[bodyToInertia.{diskr Cos[theta],diskr Sin[theta],rho},{theta,-Pi,Pi},{rho,L0,L2},
PlotStyle\[Rule]Directive[Specularity[White,30],Texture[mplt]],TextureCoordinateFunction\[Rule]({#1,#3}&),
Lighting\[Rule]"Neutral",Mesh\[Rule]None,PlotRange\[Rule]All,TextureCoordinateScaling\[Rule]True]
],
*)

{Lighting -> "Neutral", FaceForm[LightGray],
Cylinder[{bodyToInertia.L0body/coneOff , bodyToInertia.L0body}, diskr/10]}, (*rod*)
Cone[{{(bodyToInertia.L0body/coneOff), {0, 0, 0}}}, diskr/10],(*botton cone*)
If[showCG, {Red, Sphere[bodyToInertia.{0, 0, cg}, 1.1 (diskr/10)]}]
}], phis, bodyToInertia.{0, 0, 1}
],

If[showTorque,
{Red, Arrow[{bodyToInertia.{0, 0, cg}, bodyToInertia.{L0, 0, cg}}]}
],

If[showH,
{Black, Arrowheads[Medium], Arrow[{{0, 0, 0}, 1.3 Norm[L2body] (omegaVectorBody/Norm[omegaVectorBody])}]}
],

If[showBodyAxes,
{
Cylinder[{bodyToInertia.{0, 0, -.2}, bodyToInertia.{0, 0, -.1 .5}}, 1.5],
{Red, Arrowheads[Small], Arrow[{{0, 0, 0}, bodyToInertia.{1, 0, 0}}]},
Inset[Graphics[
Text["x"]
], 1.1 bodyToInertia.{1, 0, 0}],

{Red, Arrowheads[Small], Arrow[{{0, 0, 0}, bodyToInertia.{0, 1, 0}}]},
Inset[Graphics[
Text["y"]
], 1.1 bodyToInertia.{0, 1, 0}],

{Red, Arrowheads[Small], Arrow[{{0, 0, 0}, bodyToInertia.{0, 0, 1}}]},
Inset[Graphics[
Text["z"]
], 1.1 bodyToInertia.{0, 0, 1}]
}],

If[showPath,
{Red, Line[tipTimeHistory[[1 ;; timeHistoryIndex]]]}
],

If[showInertiaAxes,
{
{Arrowheads[Small], Arrow[{{0, 0, 0}, {1, 0, 0}}]},
Inset[Graphics[
Text["x"]
], {1.1, 0, 0}],

{Arrowheads[Small], Arrow[{{0, 0, 0}, {0, 1, 0}}]},
Inset[Graphics[
Text["y"]
], {0, 1.1, 0}],

{Arrowheads[Small], Arrow[{{0, 0, 0}, {0, 0, 1}}]},

Inset[Graphics[
Text["z"]
], {0, 0, 1.1}]
}],

If[showW,
{Red, Arrowheads[Medium], Arrow[{{0, 0, 0}, 1.2 Norm[L2body] (omegaVectorBody/Norm[omegaVectorBody])}]}
],

If[showWC,
{
{Blue, Arrowheads[Small], Arrow[{{0, 0, 0}, 1.2 Norm[L2body] ({omegaVectorBody[[1]], 0, 0}/Norm[omegaVectorBody])}]},
{Blue, Arrowheads[Small], Arrow[{{0, 0, 0}, 1.2 Norm[L2body] ({0, omegaVectorBody[[2]], 0}/Norm[omegaVectorBody])}]},
{Blue, Arrowheads[Small], Arrow[{{0, 0, 0}, 1.2 Norm[L2body] ({0, 0, omegaVectorBody[[3]]}/Norm[omegaVectorBody])}]}
}
],

{Lighting -> "Neutral", FaceForm[LightGray], Cylinder[{{0, 0, -.2}, {0, 0, -.1 .5}}, 1.79 L0]},

If[showSphere,
{
{Opacity[.1], Sphere[{0, 0, 0}, 1.79 L0]},

{Red, FaceForm[None],
Cylinder[{{0, 0, 1.78 L0 Cos[theta0 Degree]}, {0, 0, 1.79 L0 Cos[theta0 Degree]}}, 1.79 L0 Sin[theta0 Degree]]},

{Red, FaceForm[None], Cylinder[{{0, 0, 1.78 L0 Cos[theta2]}, {0, 0, 1.79 L0 Cos[theta2]}}, 1.79 L0 Sin[theta2]]}
}
],

If[distToGroud <= 0, Text[Style["Crash!", 14, Red], {3, 3, 0}]]
},
PlotRange -> {{-zoom, zoom}, {-zoom, zoom}, {-.5, 1.9 L0}},
SphericalRegion -> True, ImagePadding -> .1, ImageMargins -> 0, If[showPlots, ImageSize -> 250, ImageSize -> 350],
ViewPoint -> viewPoint, Boxed -> True
];

Which[state == "RUN" || state == "STEP",
currentTheta = (theta /. sol)[delT];
currentZeta = (zeta /. sol)[delT];
currentPhi = (phi /. sol)[delT];

currentPhiD = (phi' /. sol)[delT];

If[currentTime > 0, timeHistoryIndex++];
thetaTimeHistory[[timeHistoryIndex]] = {currentTime, thetas*180/Pi};
zetaTimeHistory[[timeHistoryIndex]] = {currentTime, zetas*180/Pi};
tmp = (bodyToInertia.{0, 0, cg});
tipTimeHistory[[timeHistoryIndex]] = 1.8 L0 (tmp/Norm[tmp]);

If[timeHistoryIndex == 500,
timeHistoryIndex = 1;
currentTime = 0
,
currentTime = currentTime + delT
];

If[state == "RUN" && distToGroud > 0,
tick = Not[tick]
]
];

PE = (mass1 + mass2)*g*(cg) Cos[thetas];
KE = 1/2 Iz (initialSpin + zetaW Cos[thetas])^2 + I0 (thetaW^2 + zetaW^2 Sin[thetas]^2);
If[showPlots, plots =
Row[{ListLinePlot[ thetaTimeHistory[[1 ;; timeHistoryIndex]], PlotRange -> {{0, 500*delT}, {0.9 theta0, 1.1 theta2*180/Pi}} ,
Frame -> True, AspectRatio -> 0.3, FrameLabel -> {{"\[Theta]", None}, {"time (sec)", "Nutation angle (deg) vs. time"}},
ImageSize -> 250, ImagePadding -> {{45, 10}, {30, 45}}, ImageMargins -> 0],
ListLinePlot[ zetaTimeHistory[[1 ;; timeHistoryIndex]], PlotRange -> {{0, 500*delT}, {0, 360}} , Frame -> True,
AspectRatio -> 0.3, FrameLabel -> {{"\[Psi]", None}, {"time (sec)", "precession angle (deg) vs. time"}}, ImageSize -> 250,
ImagePadding -> {{45, 10}, {30, 45}}, ImageMargins -> 0]
}]
];

Grid[{
{Grid[{
{"c.g.", "w1", "w2", "w3", "P.E.", "K.E.", "total energy"},
},
{"spin (hz)", "\!\(\*SubscriptBox[\(\[Theta]\), \(1\)]\)", "\!\(\*SubscriptBox[\(\[Theta]\), \(2\)]\)", "\[Psi]'(hz)",
"\[Theta](t)", "n"},
{If[showPlots,
Grid[{
{plots},
{graph}
}],
Grid[{
{graph}
}]
], SpanFromLeft
}
}, Frame -> All, FrameStyle -> Gray]
}
}]],

Grid[{
{
Grid[{
{ Button[Text@Style["run", 12], {state = "RUN"; tick = Not[tick]}, ImageSize -> {60, 40}],
Button[Text@Style["step", 12], {state = "STEP"; tick = Not[tick]}, ImageSize -> {60, 40}],
Button[Text@Style["stop", 12], {state = "STOP"; tick = Not[tick]}, ImageSize -> {60, 40}],
Button[Text@Style["reset", 12], {state = "RESET"; delT = 0.01;
density = 1;
diskr = 1;
diskThick = .6;
L0 = 2.5;
showCG = True;
nStrips = 1;
currentTheta = 35*Pi/180;
theta0 = 35;
currentZeta = 0;
currentPhi = 0;
currentPhiD = 2 Pi*15;
phiD0 = 15;
showWC = False;
showW = False;
showInertiaAxes = True;
showBodyAxes = True;
delT = 0.035;
timeHistoryIndex = 1;
showPlots = False;
currentTime = 0;
showPath = True;
showH = False;
showTorque = False;
tick = Not[tick]}, ImageSize -> {60, 40}](*fix*)
}
}, Spacings -> {.5, 0}, Frame -> True, FrameStyle -> Gray
]
},
{
Text@Grid[{
{"Initial disk spin (hz)",
Manipulator[Dynamic[phiD0, {phiD0 = #, currentPhiD = phiD0 2*Pi; timeHistoryIndex = 1; currentTime = 0;
currentTheta = theta0*Pi/180;
currentZeta = 0;
currentPhi = 0;
tick = Not[tick]} &], {1, 50, .1}, ImageSize -> Small], Dynamic[padIt1[phiD0, {3, 1}]]},

{"Initial disk angle",
Manipulator[Dynamic[theta0, {theta0 = #, currentTheta = theta0*Pi/180; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentPhiD = phiD0 2*Pi;
tick = Not[tick]} &], {1, 45, 1}, ImageSize -> Small], Dynamic[padIt1[theta0, 2]]},

{"simulation time step",
Manipulator[Dynamic[delT, {delT = #; timeHistoryIndex = 1; currentTime = 0; tick = Not[tick]} &], {.001, 0.05, .001},
ImageSize -> Small], Dynamic[padIt1[delT, {4, 3}]]},

{"disk density", Manipulator[Dynamic[density, {density = #; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {1, 100, 1}, ImageSize -> Small], Dynamic[padIt1[density, 3]]},

{"disk radius", Manipulator[Dynamic[diskr, {diskr = #; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {.1, 1, .1}, ImageSize -> Small], Dynamic[padIt1[diskr, {2, 1}]]},

{"disk thickness", Manipulator[Dynamic[diskThick, {diskThick = #; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {.01, .6, .01}, ImageSize -> Small], Dynamic[padIt1[diskThick, {2, 1}]]},

{"bar length", Manipulator[Dynamic[L0, {L0 = #; timeHistoryIndex = 1; currentTime = 0;
currentZeta = 0;
currentPhi = 0;
currentTheta = theta0*Pi/180;

currentPhiD = phiD0 2*Pi;

tick = Not[tick]} &], {1, 2.5, .1}, ImageSize -> Small], Dynamic[padIt1[L0, {2, 1}]]},

{"rotor opacity", Manipulator[Dynamic[op, {op = #, tick = Not[tick]} &], {0, 1, .01}, ImageSize -> Small],
{"select viewpoint",
SetterBar[
Dynamic[viewPoint, {viewPoint = #, tick = Not[tick]} &], {{2, 0, 3} -> 1, {1, -2, 1} -> 2, {0, -2, 2} -> 3, {-2, -2, 0} ->
4, {2, -2, 0} -> 5, {Pi, Pi/2, 2} -> 6}
]
},
{"zoom",
Manipulator[Dynamic[zoom, {zoom = #, tick = Not[tick]} &], {1, 8, .1}, ImageSize -> Small],
""
},
{Grid[{
{"show c.g.", Checkbox[Dynamic[showCG, {showCG = #; tick = Not[tick]} &]],
"show w components", Checkbox[Dynamic[showWC, {showWC = #; tick = Not[tick]} &]]},
{"show w", Checkbox[Dynamic[showW, {showW = #; tick = Not[tick]} &]],
"show inertial axes", Checkbox[Dynamic[showInertiaAxes, {showInertiaAxes = #; tick = Not[tick]} &]]},
{"show body axes", Checkbox[Dynamic[showBodyAxes, {showBodyAxes = #; tick = Not[tick]} &]],
"show plots", Checkbox[Dynamic[showPlots, {showPlots = #; tick = Not[tick]} &]]},
{"show space path", Checkbox[Dynamic[showPath, {showPath = #; tick = Not[tick]} &]],
"show sphere", Checkbox[Dynamic[showSphere, {showSphere = #; tick = Not[tick]} &]]},
{"show torque vector", Checkbox[Dynamic[showTorque, {showTorque = #; tick = Not[tick]} &]],
"show angular momentum", Checkbox[Dynamic[showH, {showH = #; tick = Not[tick]} &]]}
}, Frame -> All, FrameStyle -> Gray], SpanFromLeft
}
}, Frame -> True, FrameStyle -> Gray, Alignment -> Left
]
}
}, Alignment -> Left],

{{showH, False}, None},
{{showTorque, False}, None},
{{showPath, True}, None},
{{showSphere, True}, None},
{{tick, False}, None},
{{showPlots, False}, None},
{{currentTime, 0}, None},
{{density, 1}, None},
{{diskr, 1}, None},
{{diskThick, .6}, None},
{{L0, 2.5}, None},
{{state, "STOP"}, None},
{{showCG, True}, None},
{{nStrips, 1}, None},
{{zoom, 4.6}, None},
{{op, 1}, None},
{{showWC, False}, None},
{{showW, False}, None},
{{showInertiaAxes, True}, None},
{{showBodyAxes, True}, None},

{{theta0, 35}, None},
{{currentTheta, 35.0*Pi/180}, None},
{{currentZeta, 0}, None},
{{currentPhi, 0}, None},
{{currentPhiD, 2 Pi*13}, None},
{{phiD0, 13}, None},

{{viewPoint, {Pi, Pi/2, 2}}, None},
{{delT, 0.035}, None},
{{thetaTimeHistory, Table[{0, 0}, {500}]}, None},
{{zetaTimeHistory, Table[{0, 0}, {500}]}, None},
{{tipTimeHistory, Table[{0, 0, 0}, {500}]}, None},
{{timeHistoryIndex, 1}, None},

SynchronousUpdating -> True,
ControlPlacement -> Left, Alignment -> Center, ImageMargins -> 0, FrameMargins -> 0,
TrackedSymbols :> {tick},
Initialization :>
(
mplt = MatrixPlot[Table[Sin[x y/100], {x, -5, 5}, {y, -5, 5}], Frame -> False, ImagePadding -> 0, PlotRangePadding -> 0];
integerStrictPositive = (IntegerQ[#] && # > 0 &);
integerPositive = (IntegerQ[#] && # >= 0 &);
numericStrictPositive = (Element[#, Reals] && # > 0 &);
numericPositive = (Element[#, Reals] && # >= 0 &);
numericStrictNegative = (Element[#, Reals] && # < 0 &);
numericNegative = (Element[#, Reals] && # <= 0 &);
bool = (Element[#, Booleans] &);
numeric = (Element[#, Reals] &);
integer = (Element[#, Integers] &);
(*--------------------------------------------*)
f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, SignPadding -> True];
(*--------------------------------------------*)