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];
   
   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 -> {{"\[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"},
       {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)", "\!\(\*SubscriptBox[\(\[Theta]\), \(1\)]\)", "\!\(\*SubscriptBox[\(\[Theta]\), \(2\)]\)", "\[Psi]'(hz)", 
        "\[Theta](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 -> {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;
         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 -> {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;
            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 -> 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] &);
   (*--------------------------------------------*)
   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*)
 ]