Manipulate[
(*version 9/18/2023*)
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]]+2zeta'[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*)

(*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;

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

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]},

{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];

currentThetaD=(theta'/.sol)[delT];
currentZetaD=(zeta'/.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->{{"θ",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->{{"ψ",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"},
{padIt2[cg,{3,2}],
padIt2[omegaVector[[1]],{5,2}],
padIt2[omegaVector[[2]],{5,2}],
padIt2[omegaVector[[3]],{5,2}],
padIt1[PE,{9,2}],
padIt2[KE,{9,2}],
padIt1[PE+KE,{9,2}]
},
{"spin (hz)","(Subscript[θ, 1])","(Subscript[θ, 2])","ψ'(hz)","θ(t)","n"},
{padIt1[phiW/(2*Pi),{5,2}],
padIt2[theta0,{4,2}],
padIt2[theta2*180/Pi,{4,2}],
padIt1[zetaW/(2*Pi),{4,2}],
padIt1[currentTheta*180/Pi,{4,2}],
padIt1[n,{4,2}]},
{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->{50,40}],
Button[Text@Style["step",12],{state="STEP";tick=Not[tick]},ImageSize->{50,40}],
Button[Text@Style["stop",12],{state="STOP";tick=Not[tick]},ImageSize->{50,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;
currentThetaD=0;
currentZetaD=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->{50,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;
currentThetaD=0;
currentZetaD=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;
currentThetaD=0;
currentZetaD=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;

currentZetaD=0;
currentPhiD=phiD0 2*Pi;
currentThetaD=0;

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;

currentZetaD=0;
currentPhiD=phiD0 2*Pi;
currentThetaD=0;

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;

currentZetaD=0;
currentPhiD=phiD0 2*Pi;
currentThetaD=0;

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;

currentZetaD=0;
currentPhiD=phiD0 2*Pi;
currentThetaD=0;

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],Dynamic[padIt1[op,{2,1}]]},
{"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},
{{currentThetaD,0},None},
{{currentZetaD,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->False,
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]&);
(*--------------------------------------------*)
padIt1[v_?numeric,f_List]:=AccountingForm[v,
f,NumberSigns->{"-","+"},NumberPadding->{"0","0"},SignPadding->True];
(*--------------------------------------------*)
padIt1[v_?numeric,f_Integer]:=AccountingForm[Chop[v],
f,NumberSigns->{"-","+"},NumberPadding->{"0","0"},SignPadding->True];
(*--------------------------------------------*)
padIt2[v_?numeric,f_List]:=AccountingForm[v,
f,NumberSigns->{"",""},NumberPadding->{"0","0"},SignPadding->True];
(*--------------------------------------------*)
padIt2[v_?numeric,f_Integer]:=AccountingForm[Chop[v],
f,NumberSigns->{"",""},NumberPadding->{"0","0"},SignPadding->True]
(*--------------------------------------------*)
)
(*reference: page 238, applied mechanics, vol 2, dynamics, by Housner and Hudon*)
]