(*By Nasser M. Abbasi. 6/28/2014*)
Manipulate[
 (*by Nasser M. Abbasi, 6/28/148*)
 tick;
 Module[{mass = m*h*Pi r^2, eq1, eq2, eq3, r2, angle, t, to, Ix, Iy, Iz, thetax, thetay, thetaz, sol, g, debug = False, lengthOfw},
  Iz = (1/2) mass r^2;
  Iy = (1/4)  mass r^2 + 1/12 mass h^2;
  Ix = Iy;
  
  If[state == "RESET",
   currentThetax = 0;
   currentThetay = 0;
   currentThetaz = 0;
   Which[
    spin == "X-axes", currentThetaDotx = initialSpinRate*2*Pi/60; currentThetaDoty = 0; currentThetaDotz = 0,
    spin == "Y-axes", currentThetaDoty = initialSpinRate*2*Pi/60; currentThetaDotx = 0; currentThetaDotz = 0,
    spin == "Z-axes", currentThetaDotz = initialSpinRate*2*Pi/60; currentThetaDoty = 0; currentThetaDotx = 0
    ]
   ];
  
  (*Euler equations for 3D, zero torque*)
  eq1 = Ix *thetax''[t] + (Iz - Iy)*thetay'[t] thetaz'[t] == 0;
  eq2 = Iy *thetay''[t] + (Ix - Iz)*thetaz'[t] thetax'[t] == 0;
  eq3 = Iz *thetaz''[t] + (Iy - Ix)*thetay'[t] thetax'[t] == 0;
  
  sol = First@NDSolve[{eq1, eq2, eq3,
      thetax[0] == currentThetax,
      thetay[0] == currentThetay,
      thetaz[0] == currentThetaz,
      thetax'[0] == currentThetaDotx,
      thetay'[0] == currentThetaDoty,
      thetaz'[0] == currentThetaDotz}, {thetax, thetay, thetaz, thetax', thetay', thetaz'}, {t, 0, delT}];
  thetax = thetax /. sol;
  thetay = thetay /. sol;
  thetaz = thetaz /. sol;
  lengthOfw = Norm[{currentThetaDotx, currentThetaDoty, currentThetaDotz}];
  If[lengthOfw <= $MachineEpsilon, lengthOfw = 1];
  
  r2 = Which[
    spin == "X-axes",
    angle = (Dot[{1, 0, 0}, {currentThetaDotx, currentThetaDoty, currentThetaDotz}])/lengthOfw;
    to = {angle, 0, 0};
    r2 = lengthOfw*Sin[ArcCos[ angle]],
    
    spin == "Y-axes",
    angle = (Dot[{0, 1, 0}, {currentThetaDotx, currentThetaDoty, currentThetaDotz}])/lengthOfw;
    to = {0, angle, 0};
    r2 = lengthOfw*Sin[ArcCos[ angle]],
    
    spin == "Z-axes",
    angle = (Dot[{0, 0, 1}, {currentThetaDotx, currentThetaDoty, currentThetaDotz}])/lengthOfw;
    to = {0, 0, angle};
    r2 = lengthOfw*Sin[ArcCos[ angle]]
    ];
  
  g = Grid[{
     {Grid[{
        {"time (sec)", "Ix", "Iy", "Iz", "body cone radius"},
        {padIt2[currentTime, {5, 2}],
         padIt2[Ix, {4, 3}],
         padIt2[Iy, {4, 3}],
         padIt2[Iz, {4, 3}],
         padIt2[r, {3, 2}]
         }
        }, Alignment -> Center, Frame -> All, Spacings -> {.5, .7}]
      }
     ,
     {Grid[{
        {"\[Theta]x(deg)", "\[Theta]y(deg)", "\[Theta]z(deg)", "\[Theta]x'(rpm)", "\[Theta]'y(rpm)", "\[Theta]'z(rpm)"},
        {
         padIt2[N@Mod[currentThetax*180/Pi, 360], {4, 1}],
         padIt2[N@Mod[currentThetay*180/Pi, 360], {4, 1}],
         padIt2[N@Mod[currentThetaz*180/Pi, 360], {4, 1}],
         padIt1[N[currentThetaDotx/(2*Pi)] *60, {4, 2}],(*RPM*)
         padIt1[N[currentThetaDoty/(2*Pi)]*60, {4, 2}],
         padIt1[N[currentThetaDotz/(2*Pi) ]*60, {4, 2}]
         }
        }, Alignment -> Center, Frame -> All, Spacings -> {.5, .7}
       ]
      },
     {Framed@
       Graphics3D[
        Rotate[GraphicsGroup[Rotate[GraphicsGroup[Rotate[GraphicsGroup
              [{
               
               If[showAxes,
                {
                 {Red, Arrowheads[Medium], Arrow[{{0, 0, 0}, {.7, 0, 0}}]},
                 
                 Inset[Graphics[
                   Text[Style["x", Red, 14]]
                   ], {0.75, 0, 0}],
                 
                 
                 {Red, Arrowheads[Medium], Arrow[{{0, 0, 0}, {0, .7, 0}}]},
                 
                 Inset[Graphics[
                   Text[Style["y", Red, 14]]
                   ], {0, 0.75, 0}],
                 
                 {Red, Arrowheads[Medium], Arrow[{{0, 0, 0}, {0, 0, .7}}]},
                 
                 Inset[Graphics[
                   Text[Style["z", Red, 14]]
                   ], {0, 0, 0.75}]
                 }
                ],
               
               {Opacity[op], Cylinder[{{0, 0, -h/2}, {0, 0, h/2}}, r]},
               Sphere[{0, 0, 0}, .02],
               
               If[showW,
                {Blue, Arrow[{{0, 0, 0}, {currentThetaDotx, currentThetaDoty, currentThetaDotz}/lengthOfw}]}
                ],
               
               
               If[showCone,
                {EdgeForm[Red], Opacity[.1], Cone[{to, {0, 0, 0}}, r2]}
                ]
               
               }], thetax[0], {1, 0, 0}
             ]
            ], thetay[0], {0, 1, 0}]
          ], thetaz[0], {0, 0, 1}
         ],
        Axes -> False, AxesLabel -> {"x", "y", "z"}, PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}, {-1.1, 1.1}},
        SphericalRegion -> True, Boxed -> False, ImagePadding -> 2, ImageSize -> 350
        ]
      }
     }
    ];
  
  currentThetax = thetax[delT];
  currentThetay = thetay[delT];
  currentThetaz = thetaz[delT];
  
  currentThetaDotx = (thetax' /. sol)[delT];
  currentThetaDoty = (thetay' /. sol)[delT];
  currentThetaDotz = (thetaz' /. sol)[delT];
  
  If[debug, Print["currentThetax=", currentThetax]];
  If[debug, Print["currentThetay=", currentThetay]];
  If[debug, Print["currentThetaz=", currentThetaz]];
  
  Which[state == "RUN",
   currentTime += delT;
   If[currentTime >= 1000, currentTime = 0];
   tick = Not[tick]
   ];
  
  If[Abs[currentThetaDotx] > 10 || Abs[currentThetaDoty] > 10 || Abs[currentThetaDotz] > 10,
   state = "STOP"
   ];
  
  g
  ],
 
 Grid[{
   {
    Grid[{
      {"radius", Manipulator[Dynamic[r, {r = #, tick = Not[tick]} &], {.1, .5, .1}, ImageSize -> Small], Dynamic[padIt2[r, {2, 1}]]},
      {"height", Manipulator[Dynamic[h, {h = #, tick = Not[tick]} &], {.1, 1, .1}, ImageSize -> Small], Dynamic[padIt2[h, {2, 1}]]},
      {"density", Manipulator[Dynamic[m, {m = #, tick = Not[tick]} &], {.1, 50, .1}, ImageSize -> Small], Dynamic[padIt2[m, {3, 1}]]}
      }, Frame -> True, FrameStyle -> Gray
     ]
    },
   (*
   {
   Row[{"Initial angular positions"}]
   },
   {
   Grid[{
   {"\[Theta]x(0)",Manipulator[Dynamic[\[Theta]x,{\[Theta]x=#;currentThetax=\[Theta]x*Pi/180;tick=Not[tick]}&],{-15,15,1},
   ImageSize\[Rule]Small],Dynamic[padIt1[\[Theta]x,2]]},
   {"\[Theta]y(0)",Manipulator[Dynamic[\[Theta]y,{\[Theta]y=#;currentThetay=\[Theta]y*Pi/180;tick=Not[tick]}&],{-15,15,1},
   ImageSize\[Rule]Small],Dynamic[padIt1[\[Theta]y,2]]},
   {"\[Theta]z(0)",Manipulator[Dynamic[\[Theta]z,{\[Theta]z=#;currentThetaz=\[Theta]z*Pi/180;tick=Not[tick]}&],{-15,15,1},
   ImageSize\[Rule]Small],Dynamic[padIt1[\[Theta]z,2]]}
   },Frame\[Rule]True]
   },
   *)
   {
    Grid[{
      {Style["select spin axes", 12],
       PopupMenu[Dynamic[spin, {spin = #;
           Which[
            spin == "X-axes", currentThetaDotx = initialSpinRate*2*Pi/60; currentThetaDoty = 0; currentThetaDotz = 0; 
            currentThetax = 0; currentThetay = 0; currentThetaz = 0,
            
            spin == "Y-axes", currentThetaDoty = initialSpinRate*2*Pi/60; currentThetaDotx = 0; currentThetaDotz = 0; 
            currentThetax = 0; currentThetay = 0; currentThetaz = 0,
            
            spin == "Z-axes", currentThetaDotz = initialSpinRate*2*Pi/60; currentThetaDotx = 0; currentThetaDoty = 0; 
            currentThetax = 0; currentThetay = 0; currentThetaz = 0
            
            ];
           tick = Not[tick]} &], {"X-axes", "Y-axes", "Z-axes"},
        ImageSize -> All]
       },
      {Style["Initial spin rate (RPM)", 10], SpanFromLeft},
      {Manipulator[Dynamic[initialSpinRate, (initialSpinRate = #;
           Which[
            spin == "X-axes", currentThetaDotx = initialSpinRate*2*Pi/60; currentThetaDoty = 0; currentThetaDotz = 0,
            
            spin == "Y-axes", currentThetaDoty = initialSpinRate*2*Pi/60; currentThetaDotx = 0; currentThetaDotz = 0,
            
            spin == "Z-axes", currentThetaDotz = initialSpinRate*2*Pi/60; currentThetaDoty = 0; currentThetaDotx = 0
            ];
           tick = Not[tick]) &], {-10, 10, .1}, ImageSize -> Small
        ],
       Dynamic[padIt1[N@initialSpinRate, {4, 2}]]
       }
      }, Frame -> True, FrameStyle -> Gray
     ]
    },
   {
    Grid[{
      { Button[Text@Style["run", 12], {state = "RUN"; tick = Not[tick]}, ImageSize -> {40, 40}],
       Button[Text@Style["step", 12], {state = "STEP"; tick = Not[tick]}, ImageSize -> {40, 40}],
       Button[Text@Style["stop", 12], {state = "STOP"; tick = Not[tick]}, ImageSize -> {40, 40}],
       Button[Text@Style["reset", 12], {state = "RESET"; currentThetax = 0; currentThetay = 0; currentThetaz = 0; 
         currentThetaDotx = 0; currentThetaDoty = 0; currentThetaDotz = 0; op = 1; h = .5; r = .5; m = 1; spin = "Z-axes"; 
         initialSpinRate = 10; delT = 0.1; currentTime = 0; tick = Not[tick]}, ImageSize -> {40, 40}](*fix*)
       }
      }, Spacings -> {.2, 0}, Frame -> True, FrameStyle -> Gray
     ]
    },
   {
    Grid[{
      {
       Button[Text@Style["perturbe", Bold], {
         Which[
          spin == "X-axes", currentThetaDoty = (0.02*currentThetaDotx); currentThetaDotz = (0.02*currentThetaDotx),
          spin == "Y-axes", currentThetaDotx = (0.02*currentThetaDoty); currentThetaDotz = (0.02*currentThetaDoty),
          spin == "Z-axes", currentThetaDoty = (0.02*currentThetaDotz); currentThetaDotx = (0.02*currentThetaDotz)
          ];
         tick = Not[tick]}, ImageSize -> {80, 40}]
       }
      }]
    }
   ,
   {Grid[{
      {Row[{"opacity", Spacer[3], Manipulator[Dynamic[op, {op = #; tick = Not[tick]} &], {.1, 1, .1}, ImageSize -> Small], 
         Spacer[3], Dynamic[padIt1[op, {1, 1}]]}]},
      {Row[{"slow", Spacer[3], Manipulator[Dynamic[delT, {delT = #; tick = Not[tick]} &], {0.01, 0.2, .01}, ImageSize -> Small], 
         Spacer[3], "fast"}]},
      {Row[{"show Axes", Spacer[3], Checkbox[Dynamic[showAxes, {showAxes = #; tick = Not[tick]} &]]}]},
      {Row[{"show angular velocity direction", Spacer[2], Checkbox[Dynamic[showW, {showW = #; tick = Not[tick]} &]]}], SpanFromLeft},
      {Row[{"show body cone", Spacer[3], Checkbox[Dynamic[showCone, {showCone = #; tick = Not[tick]} &]]}]}
      }, Frame -> True, Alignment -> Left, FrameStyle -> Gray
     ]
    }
   }, Frame -> False, Alignment -> Center, FrameStyle -> Gray
  ],
 
 {{tick, False}, None},
 {{state, "RESET"}, None},
 {{currentTime, 0}, None},
 {{delT, 0.1}, None},
 {{op, 1}, None},
 {{spin, "Z-axes"}, None},
 {{initialSpinRate, 10}, None},
 {{showAxes, True}, None},
 {{showCone, True}, None},
 
 (*{{\[Theta]x,0},None},
 {{\[Theta]y,0},None},
 {{\[Theta]z,0},None},
 *)
 
 {{r, .5}, None},
 {{h, .5}, None},
 {{m, 1}, None},
 {{showW, True}, None},
 
 {{currentThetax, 0}, None},
 {{currentThetay, 0}, None},
 {{currentThetaz, 0}, None},
 
 {{currentThetaDotx, 0}, None},
 {{currentThetaDoty, 0}, None},
 {{currentThetaDotz, 0}, None},
 
 {{currentThetaDotDotx, 0}, None},
 {{currentThetaDotDoty, 0}, None},
 {{currentThetaDotDotz, 0}, None},
 TrackedSymbols :> {tick},
 
 ControlPlacement -> Left, Alignment -> Center, ImageMargins -> 0, FrameMargins -> 0,
 Initialization :>
  (
   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];
   (*--------------------------------------------*)
   
   getIx[m_, w_, h_, d_] := 1/12 m (h^2 + d^2);
   getIy[m_, w_, h_, d_] := 1/12 m (h^2 + w^2);
   getIz[m_, w_, h_, d_] := 1/12 m (w^2 + d^2)
   )
 ]