(*
Rigid Body Pendulum on a Flywheel
by Nasser M. Abbasi, version June 8, 2011*)

Manipulate[
 Quiet[process[L, m, massOfFlyWheel, \[Theta]0, \[Phi]0, der\[Theta]0,
    der\[Phi]0, viewPoint, currentTime], InterpolatingFunction::dmval],
 {{L, 7.5, "rod length (m)"}, 2, 7.5, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny},
 {{m, 3.5, "mass of rod (kg)"}, 1, 100, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny},
 {{massOfFlyWheel, 3.5, "mass of flywheel (kg)"}, 1, 500, .1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 Delimiter,
 {{\[Theta]0, 45, "\[Theta](0) (deg)"}, -90, 90, 1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{\[Phi]0, 90, "\[Phi](0) (deg)"}, 0, 360, 1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{der\[Theta]0, 24, "\[Theta]'(0) (deg/sec)"}, 0, 90, 1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{der\[Phi]0, 50, "\[Phi]'(0) (deg/sec)"}, 0, 90, 1, 
  Appearance -> "Labeled", ImageSize -> Tiny},
 {{animRate, 1, "animate rate"}, .1, 2, .1, Appearance -> "Labeled", 
  ImageSize -> Tiny},
 Delimiter,
 {{currentTime, 0, "start simulation"}, 0, 100, .001, 
  ControlType -> Trigger, AnimationRate -> Dynamic[animRate], 
  ImageSize -> Tiny},
 Delimiter,
 {{viewPoint, {Pi, Pi/2, 2}, 
   "select viewpoint"}, {{1.3, -2.4, 2} -> 1, {1, -2, 1} -> 
    2, {0, -2, 2} -> 3, {0, -2, -2} -> 4, {-2, -2, 0} -> 
    5, {2, -2, 0} -> 6, {Pi, Pi/2, 2} -> 7}, ControlType -> SetterBar,
   ImageSize -> Small},
 
 SynchronousUpdating -> True,
 SynchronousInitialization -> True,
 AutorunSequencing -> {9},
 
 Initialization :> 
  (
   g = 9.8;
   h1 = .4;
   r1 = 6*h1;
   h0 = 3*h1; r0 = r1/3;
   h2 = 4*h1;
   r2 = h2/2;
   r3 = r2/4;
   
   process[L_, m_, massOfFlyWheel_, \[Theta]0_, \[Phi]0_, 
     der\[Theta]0_, der\[Phi]0_, viewPoint_, currentTime_] := 
    Module[{eq1, eq2, Id, sol, sol\[Theta], 
      sol\[Phi], \[Theta], \[Phi], t, lines, i, title},
     Id = massOfFlyWheel*r1^2/2;
     
     (*The equations of motions of for the 2 generalized coordinates \
have beeen*)
     (*derived using Lagrangian energy method and will now be \
numerically solved*)
     eq1 = 
      Id \[Phi]''[
         t] + (1/3 m L^2) (\[Phi]''[t] Sin[\[Theta][t]]^2 + 
          2 \[Phi]'[t] Sin[\[Theta][t]] Cos[\[Theta][t]] \[Theta]'[t]);
     eq2 = 
      1/3 m L^2 \[Theta]''[t] - 
       1/3 m L^2 \[Phi]'[t]^2 Sin[\[Theta][t]] Cos[\[Theta][t]] + 
       m g L/2 Sin[\[Theta][t]];
     
     sol = 
      First@Quiet@
        NDSolve[{eq1 == 0, 
          eq2 == 0, \[Theta][0] == \[Theta]0 Degree, \[Theta]'[0] == 
           der\[Theta]0 Degree, \[Phi][0] == \[Phi]0 Degree, \[Phi]'[
            0] == der\[Phi]0 Degree}, {\[Theta], \[Phi]}, {t, 0, 100}];
     
     sol\[Theta] = \[Theta] /. sol;
     sol\[Phi] = \[Phi] /. sol;
     
     With[{current\[Theta] = sol\[Theta][currentTime], 
       current\[Phi] = sol\[Phi][currentTime]},
      
      (*make few lines to draw around the pendulum rod in order to \
make it *)
      (*more clear that it is spining on its axes*)
      lines = Table[Line[{
          {r3 Cos[current\[Phi] + i*2 Pi/5], 
           r3 Sin[current\[Phi] + i*2 Pi/5] , -(h1/2 + h2 - r2/2) + 
            r3 Sin[current\[Theta]] Cos[current\[Phi] + i*2 Pi/5]  },
          {r3 Cos[current\[Phi] + i*2 Pi/5] + L Sin[current\[Theta]], 
           r3 Sin[current\[Phi] + i*2 Pi/5], - L Cos[
              current\[Theta]] - (h1/2 + h2 - r2/2) + 
            r3 Sin[current\[Theta]] Cos[current\[Phi] + i*2 Pi/5] }
          }], {i, 1, 5}];
      
      title = Text@Style[Grid[{
           {"time (sec)", "\[Theta] (deg)", "\[Phi] (deg)"},
           
           {Row[{
              
              Text@PaddedForm[currentTime, {4, 2}, 
                NumberSigns -> {"-", ""}, NumberPadding -> {"0", "0"},
                 SignPadding -> True]
              }],
            
            Row[{
              
              Text@PaddedForm[current\[Theta] 180./Pi , {5, 2}, 
                NumberSigns -> {"-", "+"}, 
                NumberPadding -> {"0", "0"}, SignPadding -> True]
              }],
            
            Row[{
              
              Text@PaddedForm[Mod[current\[Phi] 180./Pi, 360], {5, 2},
                 NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"},
                 SignPadding -> True]
              }]
            }}, Frame -> All, Spacings -> 1, ItemSize -> 8]
         , 12];
      
      (*build the 3D graphics, one by one from top to bottom *)
      Grid[{
        {title},
        {
         Graphics3D[
          {
           Cylinder[{{0, 0, h1/2}, {0, 0, h1/2 + h0}}, r0],
           Cylinder[{{0, 0, h1/2}, {0, 0, -h1/2}}, r1],
            Cuboid[{-r2, -r2/2, -h2}, {r2, r2/2, 0}] ,
            Cylinder[{{0, r2, -h2}, {0, -r2, -h2}}, r2] , 
           
           {  EdgeForm[{Thick, Blue}], FaceForm[{Pink, Opacity[1]}], 
            Cylinder[{{0, 
               0, -(h1/2 + h2 - r2/2)}, {L Sin[current\[Theta]], 
               0, -(h1/2 + h2 - r2/2 + L Cos[current\[Theta]])}}, 
             r3] },
           
           (*line around cylinder*)
           {Thickness[.005], Blue, lines},
           
           {Thickness[.02], Red, 
            Line[{{0, 0, h1/2}, {r1 Cos[current\[Phi]], 
               r1 Sin[current\[Phi]], h1/2}}]},
           {Thickness[.02], Red, 
            Line[{{r1 Cos[current\[Phi]], r1 Sin[current\[Phi]], 
               h1/2}, {r1 Cos[current\[Phi]], 
               r1 Sin[current\[Phi]], -h1/2}}]}
           },
          
          PlotRange -> {{-5, 5}, {-3, 3}, {-8, 1}}, 
          ImageSize -> {330, 330}, Axes -> False, Boxed -> False, 
          AxesOrigin -> {0, 0, 0}, ImageMargins -> 2, 
          ImagePadding -> 2, PlotRangePadding -> 1,
          ViewPoint -> viewPoint, ViewCenter -> {0, 0, 0}
          ]}
        }, Alignment -> Center, Spacings -> 0, 
       Frame -> {False, -1 -> True}]
      ]
     ]
   )
 ]