(*
Spinning disk pendulum that swings on top of a rotating table
by Nasser M. Abbasi, version June 25, 2011*)

Manipulate[
 
 (
  {sol\[Theta], sol\[Phi], sol\[Psi]} = 
   getSolution[rSmall, hSmall, \[Rho]Small, rLarge, 
    hLarge, \[Rho]Large, 
    len, \[Theta]0, \[Theta]0Speed 2 Pi, \[Psi]0, \[Psi]0Speed 2 Pi, \
\[Phi]0, \[Phi]0Speed 2 Pi ];
  
  (*set parameters according to the test case to run*)
  
  If[Not[testCase == 0],
   Which[
    testCase == 
     1, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = maxrSmall; rLarge = maxrLarge; hSmall = maxhSmall; 
     hLarge = 0.5 maxhLarge; \[Theta]0 = 15; \[Theta]0Speed = 
      0.5 max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = 
      0.2 max\[Omega]; \[Phi]0 = 40; \[Phi]0Speed = .3 max\[Omega]; 
     viewPoint = {1.3, -2.4, .5}; angularMomentumOption = 0; 
     showI = True; zoom = 1; animRate = 0.03},
    
    testCase == 
     2, {len = 
      lPost; \[Rho]Small = \[Rho]min; \[Rho]Large = \[Rho]min; 
     rSmall = maxrSmall; rLarge = maxrLarge; hSmall = minhSmall; 
     hLarge = 0.5 maxhLarge; \[Theta]0 = 0; \[Theta]0Speed = 
      max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = max\[Omega]; \[Phi]0 = 
      0; \[Phi]0Speed = 0; viewPoint = {1.3, -2.4, .5}; 
     angularMomentumOption = 1; showI = False; zoom = 1; 
     animRate = 0.01},
    
    testCase == 
     3, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = maxrSmall; rLarge = 0.6 maxrLarge; 
     hSmall = 0.5 maxhSmall; 
     hLarge = maxhLarge; \[Theta]0 = 30; \[Theta]0Speed = 
      0.3 max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = 0; \[Phi]0 = 
      0; \[Phi]0Speed = 0; viewPoint = {2, -2, 0}; 
     angularMomentumOption = 1; showI = True; zoom = 1; 
     animRate = 0.05},
    
    testCase == 
     4, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = 0.65 maxrSmall; rLarge = 0.6 maxrLarge; 
     hSmall = 0.5 maxhSmall; 
     hLarge = maxhLarge; \[Theta]0 = 0; \[Theta]0Speed = 
      max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = 
      0.85 max\[Omega]; \[Phi]0 = 0; \[Phi]0Speed = max\[Omega]; 
     viewPoint = {Pi, Pi/2, 2}; angularMomentumOption = 0; 
     showI = False; zoom = .73; animRate = 0.02},
    
    testCase == 
     5, {len = lPost; \[Rho]Small = 10; \[Rho]Large = \[Rho]min; 
     rSmall = 0.65 maxrSmall; rLarge = 0.6 maxrLarge; 
     hSmall = minhSmall; 
     hLarge = maxhLarge; \[Theta]0 = 0; \[Theta]0Speed = 
      max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = max\[Omega]; \[Phi]0 = 
      0; \[Phi]0Speed = max\[Omega]; viewPoint = {Pi, Pi/2, 2}; 
     angularMomentumOption = 0; showI = False; zoom = 0.74; 
     animRate = 0.002},
    
    testCase == 
     6, {len = lPost; \[Rho]Small = 10; \[Rho]Large = \[Rho]max; 
     rSmall = 0.65 maxrSmall; rLarge = minrLarge; hSmall = minhSmall; 
     hLarge = minhLarge; \[Theta]0 = 133; \[Theta]0Speed = 
      0.5 max\[Omega]; \[Psi]0 = 280; \[Psi]0Speed = 
      min\[Omega]; \[Phi]0 = 135; \[Phi]0Speed = max\[Omega]; 
     viewPoint = {Pi, Pi/2, 2}; angularMomentumOption = 0; 
     showI = False; zoom = 0.76; animRate = 0.002},
    
    testCase == 
     7, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = maxrSmall; rLarge = maxrLarge; hSmall = minhSmall; 
     hLarge = 0.25 maxhLarge; \[Theta]0 = 
      0; \[Theta]0Speed = -0.3; \[Psi]0 = 0; \[Psi]0Speed = 
      0; \[Phi]0 = 0; \[Phi]0Speed = 0; viewPoint = {Pi, Pi/2, 2}; 
     angularMomentumOption = 0; showI = False; zoom = 1; 
     animRate = 0.06},
    
    testCase == 
     8, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = maxrSmall; rLarge = maxrLarge; hSmall = maxhSmall; 
     hLarge = 0.5 maxhLarge; \[Theta]0 = 15; \[Theta]0Speed = 
      0.5 max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = 
      0.2 max\[Omega]; \[Phi]0 = 40; \[Phi]0Speed = .3 max\[Omega]; 
     viewPoint = {1.3, -2.4, .5}; angularMomentumOption = 9; 
     showI = True; zoom = 1; animRate = 0.03},
    
    testCase == 
     9, {len = 
      lPost; \[Rho]Small = \[Rho]max; \[Rho]Large = \[Rho]min; 
     rSmall = .2 maxrSmall; rLarge = .8 maxrLarge; hSmall = maxhSmall;
      hLarge = 0.5 maxhLarge; \[Theta]0 = 0; \[Theta]0Speed = 
      max\[Omega]; \[Psi]0 = 0; \[Psi]0Speed = 0; \[Phi]0 = 
      83; \[Phi]0Speed = 0; viewPoint = {1.3, -2.4, .5}; 
     angularMomentumOption = 0; showI = False; zoom = 1; 
     animRate = 0.018}
    ]
   ];
  
  Dynamic[
   update[len, \[Rho]Small, \[Rho]Large, rSmall, rLarge, hSmall, 
    hLarge, currentTime, viewPoint, boxIt, angularMomentumOption, 
    showI, zoom, testCase, sol\[Theta], sol\[Phi], sol\[Psi], 
    traceThickness, isTraceOn],
   TrackedSymbols :> {currentTime, sol\[Theta], sol\[Phi], sol\[Psi], 
     viewPoint, boxIt, angularMomentumOption, showI, zoom, testCase}]
  
  ),
 
 (*------------------------*)                       
 (* L E F T    P A N E L   *)
 (*------------------------*)                       
 Item[
  Grid[{
    
    {
     Button[Style["min", 10], len = 0.1 lPost, ImageSize -> Tiny],
     Button[Style["max", 10], len = lPost, ImageSize -> Tiny],
     Control[{{len, lPost, Text[Style["len", fontSizeForControl]]}, 
       0.3 lPost, lPost, .1, ImageSize -> Tiny, 
       Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Rho]Small = \[Rho]min, 
      ImageSize -> Tiny],
     Button[Style["max", 10], \[Rho]Small = \[Rho]max, 
      ImageSize -> Tiny],
     Control[{{\[Rho]Small, \[Rho]max, 
        Text[Style["\!\(\*SubscriptBox[\(\[Rho]\), \(m\)]\)", 
          fontSizeForControl]]}, \[Rho]min, \[Rho]max, .1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Rho]Large = 1, ImageSize -> Tiny],
     Button[Style["max", 10], \[Rho]Large = 10, ImageSize -> Tiny],
     Control[{{\[Rho]Large, \[Rho]min, 
        Text[Style["\!\(\*SubscriptBox[\(\[Rho]\), \(M\)]\)", 
          fontSizeForControl]]}, \[Rho]min, \[Rho]max, .1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], rSmall = minrSmall, ImageSize -> Tiny],
     Button[Style["max", 10], rSmall = maxrSmall, ImageSize -> Tiny],
     Control[{{rSmall, Mean[{minrSmall, maxrSmall}], 
        Text[Style["r", Italic, fontSizeForControl]]}, minrSmall, 
       maxrSmall, .1, ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], rLarge = minrLarge, ImageSize -> Tiny],
     Button[Style["max", 10], rLarge = maxrLarge, ImageSize -> Tiny],
     Control[{{rLarge, maxrLarge, 
        Text[Style["R", Italic, fontSizeForControl]]}, minrLarge, 
       maxrLarge, .1, ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], hSmall = minhSmall, ImageSize -> Tiny],
     Button[Style["max", 10], hSmall = maxhSmall, ImageSize -> Tiny],
     Control[{{hSmall, maxhSmall, 
        Text[Style["h", Italic, fontSizeForControl]]}, minhSmall, 
       maxhSmall, .1, ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], hLarge = minhLarge, ImageSize -> Tiny],
     Button[Style["max", 10], hLarge = maxhLarge, ImageSize -> Tiny],
     Control[{{hLarge, 0.5 maxhLarge, 
        Text[Style["H", Italic, fontSizeForControl]]}, minhLarge, 
       maxhLarge, .1, ImageSize -> Tiny, Appearance -> "Labeled"}]
     }
    
    }, Spacings -> 0.5, Frame -> All, 
   FrameStyle -> Directive[Thickness[.001], Gray]], 
  ControlPlacement -> Left
  ],
 
 Item[Grid[{
    (*Initial conditions*)
    {
     Button[Style["min", 10], \[Theta]0 = 0, ImageSize -> Tiny],
     Button[Style["mid", 10], \[Theta]0 = 180, ImageSize -> Tiny],
     Control[{{\[Theta]0, 133, 
        Text[Style["\[Theta]  ", fontSizeForControl]]}, 0, 360, 1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Psi]0 = 0, ImageSize -> Tiny],
     Button[Style["mid", 10], \[Psi]0 = 180, ImageSize -> Tiny],
     Control[{{\[Psi]0, 186, 
        Text[Style["\[Psi]  ", fontSizeForControl]]}, 0, 360, 1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Phi]0 = 0, ImageSize -> Tiny],
     Button[Style["mid", 10], \[Phi]0 = 180, ImageSize -> Tiny],
     Control[{{\[Phi]0, 112, 
        Text[Style["\[Phi]  ", fontSizeForControl]]}, 0, 360, 1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    (*initial speeds*)
    {
     Button[Style["min", 10], \[Theta]0Speed = 0, ImageSize -> Tiny],
     Button[Style["max", 10], \[Theta]0Speed = max\[Omega], 
      ImageSize -> Tiny],
     Control[{{\[Theta]0Speed, 0.4 max\[Omega], 
        Text[Style[
          "\!\(\*OverscriptBox[\(\[Theta]\), \(\[Bullet]\)]\)  ", 
          fontSizeForControl]]}, min\[Omega], max\[Omega], .1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Psi]0Speed = 0, ImageSize -> Tiny],
     Button[Style["max", 10], \[Psi]0Speed = max\[Omega], 
      ImageSize -> Tiny],
     Control[{{\[Psi]0Speed, 0.2 max\[Omega], 
        Text[Style[
          "\!\(\*OverscriptBox[\(\[Psi]\), \(\[Bullet]\)]\)  ", 
          fontSizeForControl]]}, min\[Omega], max\[Omega], .1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     },
    
    {
     Button[Style["min", 10], \[Phi]0Speed = 0, ImageSize -> Tiny],
     Button[Style["max", 10], \[Phi]0Speed = max\[Omega], 
      ImageSize -> Tiny],
     Control[{{\[Phi]0Speed, .3 max\[Omega], 
        Text[Style[
          "\!\(\*OverscriptBox[\(\[Phi]\), \(\[Bullet]\)]\)  ", 
          fontSizeForControl]]}, min\[Omega], max\[Omega], .1, 
       ImageSize -> Tiny, Appearance -> "Labeled"}]
     }
    
    }, Spacings -> 0.5, Frame -> All, 
   FrameStyle -> Directive[Thickness[.001], Gray]],
  ControlPlacement -> Left
  ],
 
 Item[
  Grid[{
    
    {
     Control[{{animRate, 0.025, Style["step size", 12]}, .001, 
       0.1, .001, Appearance -> "Labeled", ImageSize -> Tiny}]
     },
    
    {
     Control[{{currentTime, 0, Style["run", 12]}, 0, 
       maxSimulationTime , Dynamic[animRate], ControlType -> Trigger, 
       DisplayAllSteps -> True, ImageSize -> Tiny, 
       AppearanceElements -> {"ProgressSlider", "ResetPlayButton", 
         "PauseButton", "StepLeftButton", "StepRightButton", 
         "ResetButton"}}]
     },
    
    {
     Grid[{
       {
        Control[{{viewPoint, {1.3, -2.4, .5}, 
           Style["viewpoint", 12]}, {{1.3, -2.4, .5} -> 
            Style["1", 11], {-0.57, -2.7, 0.063} -> 
            Style["2", 11], {2, -2, 0} -> Style["3", 11], {2, 0, 0} ->
             Style["4", 11], {0, 0, 2} -> 
            Style["5", 11], {1, -1, 1} -> 
            Style["6", 11], {3.96, -0.54, 0.54} -> 
            Style["7", 11], {Pi, Pi/2, 2} -> Style["8", 11]}, 
          ControlType -> SetterBar, ImageSize -> Small}]
        }
       }, Spacings -> .5, Frame -> None]
     },
    
    {
     Grid[{
       {
        
        Grid[{{
           Style["display", 11],
           Control[{
              {angularMomentumOption, 0, ""},
             {0 -> Text@Style["bob only", 11],
              1 -> Text@Style["L", Italic, 11],
              
              2 -> Text@
                Style["{\!\(\*SubscriptBox[\(L\), \
\(x\)]\),\!\(\*SubscriptBox[\(L\), \(y\)]\),\!\(\*SubscriptBox[\(L\), \
\(z\)]\)}", Italic, 10],
              
              3 -> Text@
                Row[{Style["L", Italic, 10], " + ", 
                  Style["L", Italic, 10], Style[" components", 10]}],
              
              4 -> Text@
                Row[{Style["d", 10], Style["L", Italic, 10], 
                  Style["/d", 10], Style["t", Italic, 10]}],
              
              5 -> Text@
                Row[{Style["d", 10], Style["L", Italic, 10], 
                  Style["/d", 10], Style["t", Italic, 10], 
                  Style[" components", 10]}],
              6 -> Text@Style["\[Omega]", 10],
              
              7 -> Text@
                Style["{\!\(\*SubscriptBox[\(\[Omega]\), \
\(x\)]\),\!\(\*SubscriptBox[\(\[Omega]\), \
\(y\)]\),\!\(\*SubscriptBox[\(\[Omega]\), \(z\)]\)}", 10],
              
              8 -> Text@
                Row[{Style["\[Omega] and ", 10], 
                  Style["L", Italic, 10]}],
              
              9 -> Text@
                Row[{Style["d{\[Theta],\[Psi],\[Phi]}/d", 10], 
                  Style["t", Italic, 10]}],
              
              10 -> Text@
                Row[{Style["d", 10], Style["L", Italic, 10], 
                  Style["/d", 10], Style["t", Italic, 10], 
                  Style[" and ", 10], Style["L", Italic, 10]}],
              
              11 -> Text@
                Row[{Style["d", 10], Style["L", Italic, 10], 
                  Style["/d", 10], Style["t", Italic, 10], 
                  Style[" and ", 10], Style["L", Italic, 10], 
                  Style[" and \[Omega]", 10]}]
              }, ControlType -> PopupMenu, ImageSize -> All}
            ]}}, Spacings -> -.5, Frame -> None, Alignment -> Center],
        
        Grid[{
          {
           Style["test", 11],
           Control[{ {testCase, 0, ""},
             {0 -> Style["0", Small],
              1 -> Style["1", Small],
              2 -> Style["2", Small],
              3 -> Style["3", Small],
              4 -> Style["4", Small],
              5 -> Style["5", Small],
              6 -> Style["6", Small],
              7 -> Style["7", Small],
              8 -> Style["8", Small],
              9 -> Style["9", Small]
              }, ControlType -> PopupMenu, ImageSize -> All}
            ]}}, Spacings -> -.5, Frame -> None, Alignment -> Center]
        
        }
       }, Spacings -> .5, Frame -> None]
     },
    
    {
     Grid[{
       {
        Style["time (sec)", 11],
        Dynamic[
         Style[PaddedForm[currentTime, {5, 3}, 
           NumberSigns -> {"-", ""}, NumberPadding -> {"0", "0"}, 
           SignPadding -> True], 11]]
        }}, Alignment -> Center, Spacings -> .8
      ]
     }
    
    }, Alignment -> Left, Spacings -> .5, Frame -> All, 
   FrameStyle -> Directive[Thickness[.001], Gray]],
  ControlPlacement -> Left
  ],
 
 (*---------------------------*)                       
 (* R I G H T    P A N E L    *)
 (*---------------------------*)    
 Item[Grid[{
    {Text[Style["zoom", 10]]},
    {Control[{ {zoom, 1, ""}, .73, 1, .01, 
       ControlType -> VerticalSlider, ImageSize -> Small}]},
    {""},
    {Text[Style["info", 11]]},
    {Control[{ {showI, True, ""}, {True, False}, 
       ControlType -> Checkbox, ImageSize -> Tiny}]},
    {Text[Style["box", 11]]},
    {Control[{ {boxIt, False, ""}, {True, False}, 
       ControlType -> Checkbox, ImageSize -> Tiny}]},
    {""},
    {Grid[{
       {Text[Style["trace", 11]]},
       {Control[{ {isTraceOn, False, ""}, {True, False}, 
          ControlType -> Checkbox, ImageSize -> Tiny}]},
       {Text@Style["length", 10]},
       {Control[{ {currentMaximumTraceSize, defaultTraceSize, ""}, 1, 
          maxTraceSize, 1, ControlType -> VerticalSlider, 
          ImageSize -> Tiny}]},
       {Text[Style["thickness", 10]]},
       {Control[{ {traceThickness, defaultTraceThickness, ""}, 0.001, 
          0.01, 0.001, ControlType -> VerticalSlider, 
          ImageSize -> Small}]}
       }, Frame -> True, FrameStyle -> {Thin, Gray}, Spacings -> 0]
     }
    }, Alignment -> Center, Frame -> {{1, 2} -> True}, 
   Spacings -> .2], ControlPlacement -> Right
  ],
 
 {{sol\[Theta], {}}, ControlType -> None},
 {{sol\[Phi], {}}, ControlType -> None},
 {{sol\[Psi], {}}, ControlType -> None},
 
 {{previousTestCaseNumber, 0}, ControlType -> None},
 {{maxSimulationTime , 100}, ControlType -> None},
 {{lPost, 10}, 
  ControlType -> None}, (*length of post below main table*)
 {{rPost, 0.1 lPost}, ControlType -> None}, (*radius of post*)
 {{minhLarge, 0.1 lPost}, 
  ControlType -> None},(*minumum height of table*)
 {{maxhLarge, lPost}, 
  ControlType -> None}, (*maximum height of table*)
 {{minrLarge, 11 rPost}, 
  ControlType -> None},(*minumum radius of table*)
 {{maxrLarge, 20 rPost}, 
  ControlType -> None}, (*maximum radius of table*)
 {{minrSmall , 2 rPost}, 
  ControlType -> None},(*minumum radius of bob disk*) 
 {{maxrSmall , 10 rPost}, 
  ControlType -> None},(*maximum radius of bob disk*) 
 {{minhSmall, 0.1 lPost}, 
  ControlType -> None},(*minumum height of bob disk*) 
 {{maxhSmall , 0.5 lPost}, 
  ControlType -> None},(*maximum height of bob disk*) 
 {{max\[Omega], 1}, 
  ControlType -> None},(*maximum angular velocity in hz*) 
 {{min\[Omega], -1}, 
  ControlType -> None},(*minumum angular velocity in hz*) 
 {{\[Rho]min, 1}, ControlType -> None},(*minumum density kg/m^3*) 
 {{\[Rho]max, 10}, ControlType -> None},(*maximum density kg/m^3*) 
 
 (*these below is data and variables to track the center of mass of \
the pendulum*)
 (*if trace is selected *)
 {{defaultTraceThickness, 0.006}, ControlType -> None},
 {{maxTraceSize, 1000}, 
  ControlType -> None}, (*maximum trace points to keep*)
 {{defaultTraceSize, 200}, ControlType -> None},
 
 Alignment -> Center,
 SynchronousUpdating -> True,
 SynchronousInitialization -> True,
 FrameMargins -> 1,
 ImageMargins -> 1,
 Initialization :> (
   
   traceBuffer = 
    Table[0, {maxTraceSize}]; (*where to store the trace coordinates*)

   
   previousMaxTraceSize = currentMaximumTraceSize;
   isFirstScan = True;
   currentTraceSize = 0;
   isSolutionChanged = False;
   
   fontSizeForControl = 11;
   
   (*--------------------------------------------*)
   (* helper function for formatting             *)
   (*--------------------------------------------*)
   padIt1[v_, f_List] := 
    AccountingForm[Chop[v] , f, NumberSigns -> {"-", "+"}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True];
   
   (*--------------------------------------------*)
   (* helper function for formatting             *)
   (*--------------------------------------------*)
   padIt2[v_, f_List] := 
    AccountingForm[Chop[v] , f, NumberSigns -> {"", ""}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True];
   
   (*---------------------------------------------------*)
   (* main entry to find the numerical solution               *)
   (*---------------------------------------------------*)
   getSolution[rSmall_, hSmall_, \[Rho]Small_, rLarge_, 
     hLarge_, \[Rho]Large_, 
     len_, \[Theta]0_, \[Theta]0Speed_, \[Psi]0_, \[Psi]0Speed_, \
\[Phi]0_, \[Phi]0Speed_] :=
    Module[{mSmall, mLarge, Id, Icg, Io, kinetic, v, g = 9.8, 
      lagrangian, eqs, initialConditions, 
      sol, \[Theta], \[Phi], \[Psi], t},
     
     {mSmall, mLarge, Id, Icg, Io} = 
      findMassesAndMomentsOfInertia[rSmall, hSmall, \[Rho]Small, 
       rLarge, hLarge, \[Rho]Large, len];
     
     (* find the solution using numerical solver*)
     (*Find kinetic and potential energy and then the Lagrangian*)
     kinetic = 
      1/2 Id \[Phi]'[t]^2 + 
       1/2 mSmall (   (len Sin[\[Theta][t]] \[Phi]'[
              t])^2 + (len \[Theta]'[t])^2 ) + 
       1/2 Icg[[3, 
         3]] ( \[Psi]'[t] + \[Phi]'[t] Cos[\[Theta][t]])^2 + 
       1/2 Icg[[2, 2]] (\[Phi]'[t] Sin[\[Theta][t]])^2 + 
       1/2 Icg[[1, 1]] \[Theta]'[t]^2;
     
     v = len (1 -  Cos[\[Theta][t]]) mSmall g;
     lagrangian = kinetic - v;
     
     (*write down the 3 equations of motion using the above \
Lagrangian*)
     (*no generalized forces, life is simple *)
     
     eqs = 
      Apply[D[ D[lagrangian, #1], t] - D[lagrangian, #2] == 
         0 & , {{\[Theta]'[t], \[Theta][t]}, {\[Psi]'[t], \[Psi][
          t]}, {\[Phi]'[t], \[Phi][t]}}, 1];
     
     (*solve using NDSolve with the initial conditions from the user*)

     
     initialConditions = {\[Theta][0] == \[Theta]0*Pi/180, \[Theta]'[
         0] == \[Theta]0Speed, \[Psi][0] == \[Psi]0*Pi/180,
       \[Psi]'[0] == \[Psi]0Speed, \[Phi][0] == \[Phi]0*
         Pi/180, \[Phi]'[0] == \[Phi]0Speed};
     
     sol = 
      First@NDSolve[
        Flatten@{eqs, initialConditions}, {\[Theta], \[Phi], \[Psi]},
        {t, 0, maxSimulationTime}, MaxSteps -> Infinity, 
        PrecisionGoal -> 7];
     
     isSolutionChanged = True;
     
     {\[Theta] /. sol, \[Phi] /. sol, \[Psi] /. sol}
     
     ];
   
   (*---------------------------------------------------*)
   (* called before numerically solving the system            *)
   (* to calculates masses and moments of inertia             *)
   (*---------------------------------------------------*)
   findMassesAndMomentsOfInertia[rSmall_, hSmall_, \[Rho]Small_, 
     rLarge_, hLarge_, \[Rho]Large_, len_] := 
    Module[{mSmall, mLarge, Id, Icg, Io, Icg1, Icg2, Icg3, Io1, Io2, 
      Io3},
     
     (*calculate mass of small and large wheel*)
     mSmall = (Pi rSmall^2) hSmall \[Rho]Small;
     mLarge = (Pi rLarge^2) hLarge \[Rho]Large;
     
     (* moments of inertia of table around its z-axis*)
     Id = ( mLarge rLarge^2)/2;
     
     Icg1 = 1/12 mSmall (3 rSmall^2 + hSmall^2); (*Ix*)
     Icg2 = Icg1; (*Iy*)
     Icg3 = (mSmall rSmall^2)/2; (*Iz*)
     
     (*apply parallel axis theorem to find I with reference to point \
o. Point o*)
     (*point o is where the rod of the pendulum is attached to the \
frame*)
     Io1 = Icg1 + mSmall len^2;
     Io2 = Io1;
     Io3 = Icg3;
     
     Icg = {{Icg1, 0, 0}, {0, Icg2, 0}, {0, 0, Icg3}};
     Io = {{Io1, 0, 0}, {0, Io2, 0}, {0, 0, Io3}};
     
     {mSmall, mLarge, Id, Icg, Io}
     ];
   
   (*---------------------------------------*)
   (* Generate title grid                        *)
   (*---------------------------------------*)
    generateTitle[current\[Theta]_, current\[Phi]_, current\[Psi]_, 
     current\[Theta]Der_, \[Phi]Der_, \[Psi]Der_, len_, Id_, Icg_, 
     mSmall_] := 
    Module[{currentKE, currentPE, title, totalEnergy, 
      currentKEAsPercentage, currentPEAsPercentage, 
      currentKEformattedAsPercentage, currentPEformattedAsPercentage, 
      currentKEformattedAsPercentageV1, 
      currentPEformattedAsPercentageV1, g = 9.8},
     
     currentKE = 
      1/2 Id \[Phi]Der^2 + 
       1/2 mSmall (   (len Sin[
              current\[Theta]] \[Phi]Der)^2 + (len \
current\[Theta]Der)^2 ) + 
       1/2 Icg[[3, 
         3]] ( \[Psi]Der + \[Phi]Der Cos[current\[Theta]])^2 + 
       1/2 Icg[[2, 2]] (\[Phi]Der Sin[current\[Theta]])^2 + 
       1/2 Icg[[1, 1]] current\[Theta]Der^2;
     
     currentPE = len (1 -  Cos[current\[Theta]]) mSmall g;
     totalEnergy = currentKE + currentPE;
     
     If[totalEnergy <= $MachineEpsilon,(*special case, 
      system at rest*)
      {
       currentKEAsPercentage  = 0;
       currentPEAsPercentage  = 0;
       },
      {
       currentKEAsPercentage  = currentKE/totalEnergy 100;
       currentPEAsPercentage  = currentPE/totalEnergy 100;
       }
      ];
     
     currentKEformattedAsPercentage = 
      Text@Row[{padIt2[currentKEAsPercentage, {2, 1}], " %"}];
     currentPEformattedAsPercentage = 
      Text@Row[{padIt2[currentPEAsPercentage, {2, 1}], " %"}];
     currentKEformattedAsPercentageV1 = 
      Text@Row[{padIt2[currentKEAsPercentage, {2, 1}], "%"}];
     currentPEformattedAsPercentageV1 = 
      Text@Row[{padIt2[currentPEAsPercentage, {2, 1}], "%"}];
     
     title = Text@Style[Grid[{
          {
           "",
           Text["\[Theta]"],
           Text["\[Psi]"],
           Text["\[Phi]"],
           Text@Row[{Style["P.E.", Blue], " (kJ)"}],
           Text@Row[{Style["K.E.", Red], " (kJ)"}]
           },
          
          { (*angular positions*)
           Text[Style["position (deg)", 9]],
           padIt2[Mod[current\[Theta] 180./Pi, 360], {6, 3}],
           padIt2[Mod[current\[Psi] 180./Pi, 360], {6, 3}],
           padIt2[Mod[current\[Phi] 180./Pi, 360], {6, 3}],
           padIt2[currentPE/1000, {8, 0}],
           padIt2[currentKE/1000, {8, 0}]
           },
          
          {(*angular velocities*)
           Text[Style["\[Omega] (hz)", 9]],
           padIt1[current\[Theta]Der/(2. Pi), {5, 3}],
           padIt1[\[Psi]Der/(2. Pi), {5, 3}],
           padIt1[\[Phi]Der/(2. Pi), {5, 3}],
           currentPEformattedAsPercentage,
           currentKEformattedAsPercentage
           }
          
          }, Frame -> All,
         FrameStyle -> Gray,
         Spacings -> 1,
         ItemSize -> {{All, 2 ;; -1} -> 6},
         Alignment -> Center], 11
        ];
     
     {title, currentKE , currentPE, currentKEformattedAsPercentage, 
      currentPEformattedAsPercentage, currentPEAsPercentage , 
      currentKEAsPercentage , currentKEformattedAsPercentageV1, 
      currentPEformattedAsPercentageV1}
     ];
   
   (*---------------------------------------*)
   (* calculate L and L' with reference to pt    *)
   (* which is the point where the pendulum rod  *)
   (* is attached to the hanger. Also generate   *)
   (* grid table containing formatted information*)
   (*---------------------------------------*)
   calculateAngularMomentum[pt_, ptcg_, Io_, 
     scaleAmount_, \[Theta]_, \[Phi]_, \[Phi]Der_, \[Psi]Der_, \
\[Theta]Der_, \[Theta]DerDer_, \[Psi]DerDer_ , \[Phi]DerDer_] := 
    Module[{Lf, Lx, Ly, Lz, L, norm, inertiaTableDisplay, LDot, LfDot,
       LxDot, LyDot, LzDot, 
      omegaDotVector, \[Omega], \[Omega]Vector, \[Omega]xComp, \
\[Omega]yComp, \[Omega]zComp, \[Theta]Vector, \[Psi]Vector, \
\[Phi]Vector, \[Theta]VectorAnnotation, \[Psi]VectorAnnotation, \
\[Phi]VectorAnnotation, maxVelocityComponent, currentptcg, 
      LAnnotation, LDotAnnotation, \[Omega]VectorAnnotation, normL, 
      r0, r1},
     
     (* find coordinates in inertia space of the cg *)
     r0 = RotationTransform[\[Theta], {1, 0, 0}, pt];
     r1 = RotationTransform[\[Phi], {0, 0, 1}, {0, 0, 0}];
     currentptcg = r1[r0[ptcg]];
     
     (* resolve the angular veclocity of the bob along the 3 \
principal axis*)
     \[Omega] = {\[Theta]Der, \[Phi]Der Sin[\[Theta]], \[Psi]Der + \
\[Phi]Der  Cos[\[Theta]]};
     
     (* resolve the rate of angular veclocity of the bob along the 3 \
principal axis*)
     (* by taking derivative w.r.t time of the omega vector, 
     this is the angular acceleration *)
     omegaDotVector = {
       \[Theta]DerDer,
       \[Phi]DerDer Sin[\[Theta]] + \[Phi]Der Cos[\[Theta]] \
\[Theta]Der ,
       \[Psi]DerDer + \[Phi]DerDer Cos[\[Theta]] - \[Phi]Der Sin[\
\[Theta]] \[Theta]Der};
     
     (* find the angular momentum relative to fixed point 0, 
     this is the fixed point in space*)
     (* that the bob is attached to*)
     
     L = Chop[Io.\[Omega], 10^-6];
     
     (* find the rate of angular momentum relative to fixed point 0*)

     
     LDot = Io.omegaDotVector;
     
     (*due to rotation of table, 
     to find ABSOLUTE rate of angular momentum      *)
     (*we need to use Subscript[(d/dt A), absolute] = Subscript[(d/
     dt A), xyz] + cross[\[Phi],A] formula        *)
     (*where in this case A is L, 
     and the LHS above is absolute rate of change  *)
     (*note: \[Phi]Der is used below, 
     since L is on the fixed point 0, which rotates *)
     (*by \[Phi]Der relative to the ground *)
     
     LDot = Chop[LDot + Cross[{0, 0, \[Phi]Der }, L], 10^-5];
     
     inertiaTableDisplay = Text@Grid[{{
          Grid[{
            {
             Row[{Style["I", Italic, 11], " = "}],
             
             Style[TraditionalForm[{
                {padIt2[Io[[1, 1]], {9, 1}], 0, 0},
                {0, padIt2[Io[[2, 2]], {9, 1}], 0},
                {0, 0, padIt2[Io[[3, 3]], {9, 1}]}}], 11]
             }}, Spacings -> 0, Frame -> None, Alignment -> Left],
          
          Grid[{
            
            {Style[
              "\[LeftDoubleBracketingBar]\[Omega]\
\[RightDoubleBracketingBar] = ", 11],
             Style[TraditionalForm[padIt2[Norm@\[Omega], {11, 1}]], 10]
             },
            
            {Row[{Style["\[LeftDoubleBracketingBar]", 11], 
               Style["L", Italic, 11], 
               Style["\[RightDoubleBracketingBar] = ", 11]}],
             Style[TraditionalForm[padIt2[Norm@L, {11, 1}]], 10]},
            
            {Row[{Style["\[LeftDoubleBracketingBar]", 11], 
               Style["\!\(\*FractionBox[\(dL\), \(dt\)]\)", Italic, 
                11], Style["\[RightDoubleBracketingBar] = ", 11]}],
             Style[TraditionalForm[padIt2[Norm@LDot, {11, 1}]], 10]}
            }, Spacings -> 0, Frame -> None, Alignment -> Left]
          },
         
         {
          Grid[{
            {
             Style["\[Omega] = ", 11],
             Style[TraditionalForm[padIt1[List@\[Omega], {9, 1}]], 11]
             }}, Spacings -> 0, Frame -> None, Alignment -> Left], 
          SpanFromLeft
          },
         
         {
          Grid[{
            {
             
             Row[{Style["L = I", Italic, 11], 
               Style["\[Omega] = ", 11]}],
             Style[TraditionalForm[padIt1[List@L, {10, 1}]], 11]
             }}, Spacings -> 0, Frame -> None, Alignment -> Left], 
          SpanFromLeft
          },
         
         {
          Grid[{
            {
             
             Style["\!\(\*FractionBox[\(dL\), \(dt\)]\) = ", Italic, 
              11],
             Style[TraditionalForm[padIt1[List@LDot, {10, 1}]], 11]
             }}, Spacings -> 0, Frame -> None, Alignment -> Left], 
          SpanFromLeft
          }
         
         }, Frame -> All, Spacings -> 1, Alignment -> Left, 
        FrameStyle -> Gray
        ];
     
     (* generate the vector representation of \
\[Theta]',\[Psi]',\[Phi]' from the cg of the bob *)
     maxVelocityComponent = 
      Max[Abs[{\[Theta]Der, \[Psi]Der, \[Phi]Der}]];
     If[\[Theta]Der <= $MachineEpsilon, \[Theta]Vector = {0, 0, 
        0}, \[Theta]Vector = ({\[Theta]Der, 0, 0}/
         maxVelocityComponent)   1.5 scaleAmount];
     \[Theta]VectorAnnotation = 
      If[\[Theta]Der <= $MachineEpsilon, Null, 
       Text[Style[
         "\!\(\*OverscriptBox[\(\[Theta]\), \(\[Bullet]\)]\)", Red, 
         14], ptcg + \[Theta]Vector + 
         0.005 Norm[\[Theta]Vector], {0, -1}]];
     \[Theta]Vector = {ptcg, ptcg + \[Theta]Vector};
     \[Theta]Vector = {Blue, Arrowheads[0.04], 
       Arrow[Tube[\[Theta]Vector, .2]]};
     
     If[\[Psi]Der <= $MachineEpsilon, \[Psi]Vector = {0, 0, 
        0}, \[Psi]Vector = ({0, 0, \[Psi]Der}/
         maxVelocityComponent)    1.5 scaleAmount];
     \[Psi]VectorAnnotation = 
      If[\[Psi]Der <= $MachineEpsilon, Null, 
       Text[
        Style["\!\(\*OverscriptBox[\(\[Psi]\), \(\[Bullet]\)]\)", Red,
          14], ptcg + \[Psi]Vector + 
         0.005  Norm[\[Psi]Vector], {0, -1}]];
     \[Psi]Vector = {ptcg, ptcg + \[Psi]Vector};
     \[Psi]Vector = {Blue, Arrowheads[0.04], 
       Arrow[Tube[\[Psi]Vector, .2]]};
     
     If[\[Phi]Der <= $MachineEpsilon, \[Phi]Vector = {0, 0, 
        0}, \[Phi]Vector = ({0, 0, \[Phi]Der}/
         maxVelocityComponent)   1.5 scaleAmount];
     \[Phi]VectorAnnotation = 
      If[\[Phi]Der <= $MachineEpsilon, Null, 
       Text[Style["\!\(\*OverscriptBox[\(\[Phi]\), \(\[Bullet]\)]\)", 
         Red, 14], 
        currentptcg + \[Phi]Vector + 
         0.005  Norm[\[Phi]Vector], {0, -1}]];
     
     \[Phi]Vector = {currentptcg, currentptcg + \[Phi]Vector};
     \[Phi]Vector = {Blue, Arrowheads[0.04], 
       Arrow[Tube[\[Phi]Vector, .2]]};
     
     (* generate the vector \[Omega] and its components for the \
angular velocity*)
     norm = Norm[\[Omega]];
     If[norm <= 
       2 $MachineEpsilon, \[Omega] = {0, 0, 
        0}, \[Omega] = (\[Omega]/norm) 1.5 scaleAmount];
     \[Omega]Vector = {ptcg, ptcg + \[Omega]};
     \[Omega]VectorAnnotation = 
      If[norm <= $MachineEpsilon, Null, 
       Text[Style["\[Omega]", Red, 15], 
        ptcg + \[Omega] + 0.005  Norm[\[Omega]], {0, -1}]];
     
     \[Omega]Vector = {Green, Arrowheads[0.04], 
       Arrow[Tube[\[Omega]Vector, .4]]};
     \[Omega]xComp = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{ptcg, ptcg + {\[Omega][[1]], 0, 0} }, .1]]};
     \[Omega]yComp = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{ptcg, ptcg + {0, \[Omega][[2]], 0} }, .1]]};
     \[Omega]zComp = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{ptcg, ptcg + {0, 0, \[Omega][[3]]} }, .1]]};
     
     (* generate the vector and its components for the angular \
momentum L*)
     normL = Norm[L];
     If[normL <= 2 $MachineEpsilon, L = {0, 0, 0}, 
      L = (L/normL)  scaleAmount];
     Lf = {pt, pt + L};
     
     LAnnotation = 
      If[normL <= $MachineEpsilon, Null, 
       Text[Style["L", Red, 15], pt + L + 0.005  Norm[Lf], {0, -1}]];
     
     Lf = {Red, Arrowheads[0.04], Arrow[Tube[Lf, .2]]};
     Lx = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {L[[1]], 0, 0} }, .1]]};
     Ly = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {0, L[[2]], 0} }, .1]]};
     Lz = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {0, 0, L[[3]]} }, .1]]};
     
     (* generate the vector and its components for the rate of \
angular momentum dL/dt*)
     norm = Norm[LDot];
     If[norm < 10^-6 normL, norm = 0]; (* Force norm to zero. 
     Was due to some numerical errors*)
     
     If[norm <= 2 $MachineEpsilon, LDot = {0, 0, 0}, 
      LDot = (LDot/norm)* 0.8 scaleAmount];
     LfDot = {pt, pt + LDot};
     LDotAnnotation = 
      If[norm <= $MachineEpsilon, Null, 
       Text[Style["\!\(\*OverscriptBox[\(L\), \(\[Bullet]\)]\)", 
         Black, 15], pt + LDot + 0.005  Norm[LfDot], {0, -1}]];
     
     LfDot = {Black, Arrowheads[0.04], Arrow[Tube[LfDot, .2]]};
     LxDot = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {LDot[[1]], 0, 0} }, .1]]};
     LyDot = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {0, LDot[[2]], 0} }, .1]]};
     LzDot = {Blue, Arrowheads[0.03], 
       Arrow[Tube[{pt, pt + {0, 0, LDot[[3]]} }, .1]]};
     
     {inertiaTableDisplay, Lf, Lx, Ly, Lz, LfDot, LxDot, LyDot, 
      LzDot, \[Omega]Vector, \[Omega]xComp, \[Omega]yComp, \
\[Omega]zComp, \[Theta]Vector, \[Psi]Vector, \[Phi]Vector, \
\[Theta]VectorAnnotation, \[Psi]VectorAnnotation, \
\[Phi]VectorAnnotation, currentptcg, LDotAnnotation, 
      LAnnotation, \[Omega]VectorAnnotation}
     ];
   
   (*---------------------------------------*)
   (* Manage trace buffer                        *)
   (*---------------------------------------*)
   refeshTraceBuffer[tnow_, isTraceOn_] := Module[{},
     
     If[(previousMaxTraceSize != currentMaximumTraceSize || 
        isSolutionChanged == True || Length[traceBuffer] == 0  ),
      {
       isSolutionChanged = False;
       traceBuffer = Table[0, {currentMaximumTraceSize}];
       previousMaxTraceSize = currentMaximumTraceSize;
       isFirstScan = True;
       currentTraceSize = 0
       }
      ];
     
     If[tnow <= $MachineEpsilon || 
       Not[isTraceOn], {currentTraceSize = 0; isFirstScan = True}]; 
     
     ];
   
   (*---------------------------------------*)
   (* Called by Manipulate main expression  *)
   (*---------------------------------------*)
   update[len_, \[Rho]Small_, \[Rho]Large_, rSmall_, rLarge_, hSmall_,
      hLarge_, tnow_, viewPoint_, boxIt_, angularMomentumOption_, 
     showI_, zoom_, testCase_, sol\[Theta]_, sol\[Phi]_, sol\[Psi]_, 
     traceThickness_, isTraceOn_] :=
    Module[{Id, mSmall, title, gr, g1, g2, g3, pt0, pt1, pt2, pt3, 
      pt4, pt5, pt6, pt6a, pt7, pt8, pt9, pt10, pt11, pt12, pt13, 
      pt14, z0, gextraCylinderOnTopOfHanger, frameRadius = 0.6, 
      currentKE, currentPE, peke, totalScale, 
      currentKEformattedAsPercentage, currentPEformattedAsPercentage, 
      currentPEAsPercentage, currentKEAsPercentage, a, b, ghangers1, 
      ghangers2, ghangers3, gLargeCylinder, line1, gline1 , gPost , 
      gWheel, gLine2, referencePointX, referencePointY, gXYZ, labels, 
      currentKEformattedAsPercentageV1, 
      currentPEformattedAsPercentageV1, g0, inertiaTableDisplay, 
      LfDot, LxDot, LyDot, LzDot, Lf, Lx, Ly, Lz, imageSize, opacity, 
      Io, Icg, \[Omega]Vector, \[Omega]xComp, \[Omega]yComp, \
\[Omega]zComp, \[Theta]Vector, \[Psi]Vector, \[Phi]Vector, \
\[Theta]VectorAnnotation, \[Psi]VectorAnnotation, \
\[Phi]VectorAnnotation, gextraCylinderOnTopOfHangerSphere, base, 
      currentptcg, p, LDotAnnotation, 
      LAnnotation, \[Omega]VectorAnnotation, \[Theta]Der, \[Psi]Der, \
\[Phi]Der, \[Theta]DerDer, \[Psi]DerDer, \[Phi]DerDer, 
      mLarge, \[Theta], \[Psi], \[Phi], t},
     
     refeshTraceBuffer[tnow, isTraceOn];
     
     (*this value is the largest vertical value for the overall 3D \
image. Will use as *)
     (*measuring stick for zooming action and other layout to measure \
things against*)
     
     totalScale = 3.2 lPost + hLarge + len + hSmall + rSmall;
     
     (* The masses and moments of inertia are now calculated from \
user input parameters *)
     {mSmall, mLarge, Id, Icg, Io} = 
      findMassesAndMomentsOfInertia[rSmall, hSmall, \[Rho]Small, 
       rLarge, hLarge, \[Rho]Large, len];
     
     (* Use the solution passed in, which was allread found *)
     (* Evaluate the solution are the current time*)
     
     \[Theta] = Chop@sol\[Theta][tnow];
     \[Phi] = Chop@sol\[Phi][tnow];
     \[Psi] = Chop@sol\[Psi][tnow];
     \[Theta]Der = Chop@(sol\[Theta]'[t] /. t -> tnow);
     \[Psi]Der = Chop@(sol\[Psi]'[t] /. t -> tnow);
     \[Phi]Der = Chop@(sol\[Phi]'[t] /. t -> tnow);
     \[Theta]DerDer = Chop@(sol\[Theta]''[t] /. t -> tnow);
     \[Psi]DerDer = Chop@(sol\[Psi]''[t] /. t -> tnow);
     \[Phi]DerDer = Chop@(sol\[Phi]''[t] /. t -> tnow);
     
     {title, currentKE , currentPE, currentKEformattedAsPercentage, 
       currentPEformattedAsPercentage, currentPEAsPercentage , 
       currentKEAsPercentage , currentKEformattedAsPercentageV1, 
       currentPEformattedAsPercentageV1} = 
      generateTitle[\[Theta], \[Phi], \[Psi], \[Theta]Der, \[Phi]Der, \
\[Psi]Der , len, Id, Icg, mSmall];
     
     (*set the coodinates of the main points to use to draw the 3D \
graphics*)
     (*these are the coordinates of main markers in the system as \
things look at*)
     (*rest and all initial conditions are zero*)
     
     z0 = lPost + hLarge + 2 lPost;
     pt0 = {0, 0, 0}; pt1 = {0, 0, lPost}; 
     pt2 = {0, 0, lPost + hLarge}; pt3 = {0, 0, z0 - len - hSmall}; 
     pt4 = {0, 0, z0 - len}; pt5 = {0, 0, z0}; 
     pt6 = {0.95 rLarge, 0, lPost}; 
     pt6a = {rLarge, 0, lPost + hLarge}; pt7 = {0.95 rLarge, 0, z0}; 
     pt8 = {-0.95 rLarge, 0, z0}; pt9 = {- 0.95 rLarge, 0, lPost}; 
     pt10 = {rLarge, 0, lPost}; pt11 = {rSmall, 0, z0 - len}; 
     pt12 = {rSmall, 0, z0 - len - hSmall}; 
     pt13 = {0.1 rLarge, 0, z0}; pt14 = {-0.1 rLarge, 0, z0};
     
     (*only calculate angular momentum L if needed to display*)
     If[(Not[angularMomentumOption == 0] || showI || isTraceOn),
      {inertiaTableDisplay, Lf, Lx, Ly, Lz, LfDot, LxDot, LyDot, 
         LzDot, \[Omega]Vector, \[Omega]xComp, \[Omega]yComp, \
\[Omega]zComp, \[Theta]Vector, \[Psi]Vector, \[Phi]Vector, \
\[Theta]VectorAnnotation, \[Psi]VectorAnnotation, \
\[Phi]VectorAnnotation, currentptcg, LDotAnnotation, 
         LAnnotation, \[Omega]VectorAnnotation} = 
        calculateAngularMomentum[pt5, pt4, Io, 
         2 zoom rSmall, \[Theta], \[Phi], \[Phi]Der, \[Psi]Der, \
\[Theta]Der, \[Theta]DerDer, \[Psi]DerDer , \[Phi]DerDer];
      ];
     
     
     If[isTraceOn,
      {
       If[++currentTraceSize > 
         currentMaximumTraceSize, {isFirstScan = False; 
         currentTraceSize = 1}];
       If[DEBUG, 
        Print["isTraceOn True, updated currentTraceSize now=", 
         currentTraceSize]];
       traceBuffer[[currentTraceSize]] = currentptcg;
       If[DEBUG, Print["isFirstScan=", isFirstScan]];
       }
      ];
     
     (* start making the 3D graphics                   *)
     (*make the main post which the large table sits on*)
     base = {RGBColor[.1, .8, .8], 
       Cylinder[{pt0, pt0 + {0, 0, rPost}}, 5 rPost]};
     gPost = {base, Cylinder[{pt0, pt2}, rPost]};
     
     (*make the large table*)
     gLargeCylinder = {Opacity[.8], Cylinder[{pt1, pt2}, rLarge]};
     
     (*line drawn on top of table*)
     line1 = {Thickness[0], Red, Line[{pt2, pt6a, pt10}]};
     
     (* draw the hanger to attach the pendulum on*)
     opacity = 1;
     
     ghangers1 = {Opacity[.8], Cylinder[{pt6, pt7}, frameRadius]};
     ghangers2 = {Opacity[opacity], 
       Cylinder[{pt7 + {0.05 pt7[[1]], 0, 0}, 
         pt8 - {0.05 pt8[[1]], 0, 0}}, frameRadius]};
     ghangers3 = {Opacity[.8], Cylinder[{pt8, pt9}, frameRadius]};
     
     (*make the little extra pump to show where the pendulum hangs*)
     gextraCylinderOnTopOfHanger = {Opacity[opacity], 
       Cylinder[{pt13, pt14}, 3 frameRadius]};
     
     (*make the end small balls at the connection of the frame joints*)

     
     gextraCylinderOnTopOfHangerSphere = {Red, 
       Opacity[1], {Sphere[pt7, 2 frameRadius], 
        Sphere[pt8, 2 frameRadius]}};
     
     (*draw the pendulum rod itself*)
     gline1 = Cylinder[{pt4, pt5}, .4];
     
     (*draw the pendulum bob, which is a cylinder in this case*)
     gWheel = {Yellow, Opacity[.8], Cylinder[{pt3, pt4}, rSmall]};
     
     (*red line on top of the above, 
     to make it easy to see it spining*)
     gLine2 = If[angularMomentumOption == 0,
       {Thickness[.01], Red, Line[{pt4, pt11, pt12 }]},
       If[
        angularMomentumOption == 1 || angularMomentumOption == 2 || 
         angularMomentumOption == 3 || angularMomentumOption == 4 || 
         angularMomentumOption == 9 || angularMomentumOption == 10 || 
         angularMomentumOption == 11,
         {Thickness[.01], Green, Line[{pt4, pt11, pt12 }]}, Null]
       ];
     
     (*--  Now start applying the solution to rotate items    -*)
     (*--  Use Rotate and GeometricTransformation with        -*)
     (*--  the RotationTransform based on angles found from   -*)
     (*--  the numerical solution above                       -*)
     
     (* start by rotating the pendulum bob itself on its z axis*)
     g0 = Rotate[{gLine2, gWheel}, \[Psi], {0, 0, 1}, pt3];
     
     (* check what display to show.*)
     g1 = If[Not[angularMomentumOption == 0],
       Which[
        angularMomentumOption == 
         1, {Rotate[{Lf}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         2, {Rotate[{Lx, Ly, Lz}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         3, {Rotate[{Lf, Lx, Ly, Lz}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         4, {Rotate[{LfDot}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         5, {Rotate[{LxDot, LyDot, LzDot}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         6, {Rotate[{\[Omega]Vector}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         7, {Rotate[{\[Omega]xComp, \[Omega]yComp, \[Omega]zComp}, \
\[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         8, {Rotate[{Lf, 
           LAnnotation, \[Omega]Vector, \[Omega]VectorAnnotation}, \
\[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         9, {Rotate[{\[Psi]Vector, \[Psi]VectorAnnotation}, \[Psi], \
{0, 0, 1}, pt5]},
        angularMomentumOption == 
         10, {Rotate[{LfDot, Lf, LDotAnnotation, 
           LAnnotation}, \[Psi], {0, 0, 1}, pt5]},
        angularMomentumOption == 
         11, {Rotate[{LfDot, Lf, LDotAnnotation, 
           LAnnotation, \[Omega]Vector, \[Omega]VectorAnnotation}, \
\[Psi], {0, 0, 1}, pt5]}
        ],
       Null
       ];
     
     (*transform the whole pendulum with its rod by theta*)
     g2 = 
      GeometricTransformation[{gline1, g0, g1, 
        If[angularMomentumOption == 
          9, {\[Theta]Vector, \[Theta]VectorAnnotation}]}, 
       RotationTransform[\[Theta], {1, 0, 0}, pt5]];
     
     (*now rotate everything by \[Phi], the table rotation angle*)
     g3 = 
      Rotate[{gLargeCylinder, line1, ghangers1 , ghangers2 , 
        ghangers3, gextraCylinderOnTopOfHanger, 
        gextraCylinderOnTopOfHangerSphere, g2 }, \[Phi], {0, 0, 1}, 
       pt1];
     
     (*Now check which test case is run. 
     Some test cases does not need PE/KE display*)
     If[zoom == 1 && 
       Not[testCase == 3 || testCase == 1 || testCase == 7],
      {
       referencePointX = rLarge + 3 rPost;
       referencePointY = -rLarge - 4 rPost;
       
       (*make the small XYZ coordinate frame on the side for \
reference*)
       labels = {Text[
          Style["X", 12], {referencePointX + 6 rPost, referencePointY,
            0}], Text[
          Style["Y", 12], {referencePointX, referencePointY + 6 rPost,
            0}], Text[
          Style["Z", 12], {referencePointX, referencePointY, 
           6 rPost}]};
       
       gXYZ = {Thin, Red, Line[{
           {referencePointX, referencePointY, 0},
           {referencePointX + 5 rPost, referencePointY, 0},
           {referencePointX, referencePointY, 0},
           {referencePointX, referencePointY + 5 rPost, 0},
           {referencePointX, referencePointY, 0},
           {referencePointX, referencePointY, 5 rPost}}]};
       
       (*now make PE/KE illustration on the side*)
       a = currentKEAsPercentage/100 totalScale/2;
       b = currentPEAsPercentage/100 totalScale/2;
       
       referencePointX = -rLarge - 10 rPost;
       referencePointY = -rLarge - 10 rPost;
       
       (*draw the PE and KE illustrations on the side of the main \
plot*)
       peke = {
         {Red, 
          Cuboid[{referencePointX, referencePointY, 
            0}, {-rLarge - 8 rPost, -rLarge - 8 rPost, 
            currentKEAsPercentage/100 totalScale/2 }]
          },
         
         {Blue, 
          Cuboid[{-rLarge - 7 rPost, -rLarge - 7 rPost, 
            0}, {-rLarge - 4 rPost, -rLarge - 4 rPost, 
            currentPEAsPercentage/100 totalScale/2 }]
          },
         Text[
          currentKEformattedAsPercentageV1, {referencePointX, 
           referencePointY, a + 4 rPost}, {0, 0}],
         Text[
          currentPEformattedAsPercentageV1, {-rLarge - 
            4 rPost, -rLarge - 4 rPost, b + 4 rPost}, {0, 0}]
         }
       },
      peke = Null;
      gXYZ = Null;
      labels = Null
      ];
     
     (*-- Done making the graphics parts. now make the final display --*)

     
     imageSize = If[showI, {345, 270}, {345, 480}];
     
     p = If[isTraceOn,
       ListPointPlot3D[
        If[isFirstScan, traceBuffer[[1 ;; currentTraceSize]], 
         traceBuffer[[1 ;; currentMaximumTraceSize]]], 
        PlotStyle -> {PointSize[traceThickness], Blue}]
       ];
     
     gr = 
      Graphics3D[{gPost, g3, labels, gXYZ, peke, 
        If[angularMomentumOption == 
          9, {\[Phi]VectorAnnotation, \[Phi]Vector}]},
       PlotRange -> {
         {-zoom 2 Max[rLarge, len + hSmall] , 
          zoom 1.7 Max[rLarge, len + hSmall]}, {-zoom 1.9 Max[rLarge, 
            len + hSmall] , zoom 1.4 Max[rLarge, len + hSmall]},
         {If[zoom < 1, lPost, 0], totalScale}},
       ImageSize -> imageSize,
       Axes -> False,
       Boxed -> boxIt,
       AxesOrigin -> {0, 0, 0},
       ImageMargins -> 2,
       ImagePadding -> 2,
       PlotRangePadding -> 1,
       ViewPoint -> viewPoint,
       ViewAngle -> All
       ];
     
     
     If[showI,
      Grid[{{title}, {inertiaTableDisplay}, {If[isTraceOn, 
          Show[gr, p], gr]} },
       Spacings -> 0,
       Frame -> None 
       ],
      Grid[{{If[isTraceOn, Show[gr, p], gr]}  },
       Spacings -> 0,
       Frame -> None 
       ]
      
      ]
     ]
   
   )
 ]