(*by Nasser M. Abbasi, version 8/29/13*)

Manipulate[None,
 Item[Grid[{
    {
     Dynamic[Row[{
        gTick;
        forcingFrequencyInRadian = forcingFrequency*2*Pi;
        
        {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
          timeConstant} = 
         convertToStandardValues[k0, mass, forcingFrequencyInRadian, 
          forceCritical, zeta];
        
        {transient, steadyState} = 
         oneDegreeOfFreedomSolution[u0, v0, k0, zeta, f0, 
          forcingFrequencyInRadian, naturalFrequency, t];
        
        magnificationFactor = 
         calculateMagnificationFactor[zeta, f0, 
          forcingFrequencyInRadian, naturalFrequency];
        dynamicMagnificationPlot = 
         makeGenericDynamicMagnificationFactorPlot[{0.1, .7, 1, zeta},
           r];
        phasePlot = 
         makeGenericPhasePlot[{0.01, 0.1, 0.7, zeta}, 
          r];(*list of zeta values to use*)
        
        phaseLagPlot = If[phaseDiagramType == "standard",
          
          makePhaseDifferencePlot[zeta, f0, forcingFrequencyInRadian, 
           naturalFrequency]
          ,
          
          makeArgandDiagram[zeta, f0, forcingFrequencyInRadian, 
           naturalFrequency, tscale]
          ];
        
        }]
      ]
     },
    {
     Grid[{
       {Text[Style["differential equation", Bold, 11]]},
       {
        Grid[{
          {Text[
            TraditionalForm@
             Style[HoldForm[
                u''[t] + 2 \[Zeta] Subscript[\[Omega], n]  u'[t] + \!\(
\*SubsuperscriptBox[\(\[Omega]\), \(n\), \(2\)] \(u[t]\)\)] == 
               HoldForm[\[CapitalDigamma]/m] Sin[
                 HoldForm[\[CurlyPi] t]], Medium]]
           }
          }, Spacings -> {.7, .5}, Alignment -> Center, Frame -> True,
          FrameStyle -> Directive[Thickness[.005], Gray]]
        },
       {Text[Style["system parameters", Bold, 11]]},
       {Grid[{
          {
           
           Row[{Style["damping", 11], Spacer[3], 
             Text@Style[\[Zeta], 11]}] ,
           
           Manipulator[
            Dynamic[
             zeta, {zeta = #; forceCritical = False; 
               gTick += del} &], {0, 2, .01}, ImageSize -> Tiny, 
            ContinuousAction -> True
            ],
           Style[Dynamic@padIt2[zeta, {3, 2}], 11],
           
           Button[Style["1", 9], {zeta = 1; forceCritical = True; 
             gTick += del}, Appearance -> "Palette", 
            Background -> LightBlue, ImageSize -> {25, 18}], 
           SpanFromLeft
           }
          ,
          {
           
           Row[{Style["stiffness", 11], Spacer[3], 
             Text@Style["k", Italic]}] ,
           
           Manipulator[
            Dynamic[k0, {k0 = #; gTick += del} &], {1, 10, 0.1}, 
            ImageSize -> Tiny, ContinuousAction -> True],
           Style[Dynamic@padIt2[k0, {3, 1}], 11], SpanFromLeft
           }
          ,
          {
           
           Row[{Style["mass", 11], Spacer[3], 
             Text@TraditionalForm@Style[m]}],
           
           Manipulator[
            Dynamic[mass, {mass = #; gTick += del} &], {1, 10, 0.1}, 
            ImageSize -> Tiny, ContinuousAction -> True], 
           Style[Dynamic@padIt2[mass, {3, 1}], 11]
           }
          ,
          {
           Row[{Style["load", 11], Spacer[3], \[CapitalDigamma]}],
           
           Manipulator[
            Dynamic[f0, {f0 = #; gTick += del} &], {0, 10, 0.1}, 
            ImageSize -> Tiny, ContinuousAction -> True], 
           Style[Dynamic@padIt2[f0, {3, 1}], 11]
           },
          {
           Row[{Style["load", 11], Spacer[3], \[CurlyPi]}],
           
           Manipulator[
            Dynamic[
             forcingFrequency, {forcingFrequency = #; 
               forcingFrequencyInRadian = forcingFrequency*2*Pi; 
               r = forcingFrequencyInRadian/ Sqrt[k0/mass]; 
               gTick += del} &], {0, 1, 0.01}, ImageSize -> Tiny, 
            ContinuousAction -> True], 
           Style[Dynamic@padIt2[N@forcingFrequency, {4, 3}], 11],
           
           Button[Style["\!\(\*SubscriptBox[\(\[Omega]\), \(n\)]\)", 
             9], {forcingFrequencyInRadian = Sqrt[k0/mass]; 
             forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
             r = 1; gTick += del}, Appearance -> "Palette", 
            Background -> LightBlue, ImageSize -> {25, 18}], 
           SpanFromLeft
           }
          }, Alignment -> Left, Spacings -> {.3, .2}, Frame -> True, 
         FrameStyle -> Directive[Thickness[.005], Gray], 
         Alignment -> Center
         ]
        },
       {Text[Style["initial conditions", Bold, 11]]},
       {Grid[{
          {
           
           Row[{Style["initial", 11], Spacer[3], 
             Text@TraditionalForm@Style[u[t], 11]}],
           
           Manipulator[
            Dynamic[u0, {u0 = #; gTick += del} &], {-2, 2, 0.1}, 
            ImageSize -> Tiny, ContinuousAction -> True],
           Style[Dynamic@padIt1[u0, {2, 1}], 11],
           
           Button[Style["0", 9], {u0 = 0; gTick += del}, 
            Appearance -> "Palette", Background -> LightBlue, 
            ImageSize -> {25, 18}], SpanFromLeft
           }
          ,
          {
           
           Row[{Style["initial", 11], Spacer[3], 
             Text@TraditionalForm@Style[u'[t], 11]}],
           
           Manipulator[
            Dynamic[v0, {v0 = #; gTick += del} &], {-2, 2, 0.1}, 
            ImageSize -> Tiny, ContinuousAction -> True],
           Style[Dynamic@padIt1[v0, {2, 1}], 11],
           
           Button[Style["0", 9], {v0 = 0; gTick += del}, 
            Appearance -> "Palette", Background -> LightBlue, 
            ImageSize -> {25, 18}], SpanFromLeft
           }
          }, Alignment -> Left, Spacings -> {.3, .2}, Frame -> True, 
         FrameStyle -> Directive[Thickness[.005], Gray], 
         Alignment -> Center
         ]
        },
       {Text[Style["model information", Bold, 11]]},
       {
        Grid[{
          {Text@
            Style[Row[{"frequency ratio", Spacer[3], \[CurlyPi], "/", 
               Subscript[\[Omega], Style[n, Italic]]}], 11],
           Style[Dynamic@padIt2[N@r, {3, 2}], 11], Spacer[39]},
          
          {Text@
            Style[Row[{"natural frequency", Spacer[3], 
               Subscript[\[Omega], Style[n, Italic]]}], 11],
           
           Style[Dynamic@padIt2[N@naturalFrequency/(2.0*Pi), {4, 3}], 
            11], Style["Hz", 11]},
          
          {Text@
            Style[Row[{"damped frequency", Spacer[3], 
               Subscript[\[Omega], Style[d, Italic]]}], 11],
           
           Style[Dynamic@padIt2[N@dampedFrequency/(2.0*Pi), {4, 3}], 
            11], Style["Hz", 11]},
          
          {Text@
            Style[Row[{"natural period", Spacer[3], 2 \[Pi], "/", 
               Subscript[\[Omega], Style[n, Italic]]}], 11],
           
           Style[Dynamic@padIt2[(2. Pi)/naturalFrequency, {5, 3}], 
            11], Style["sec", 11]},
          
          {Text@
            Style[Row[{"damped period", Spacer[3], 2 \[Pi], "/", 
               Subscript[\[Omega], Style[d, Italic]]}], 11],
           
           Style[Dynamic[
             If[tdd === Infinity, Infinity, padIt2[N@tdd, {5, 3}]]], 
            11], Style["sec", 11]},
          
          {Text@
            Style[Row[{"damping coefficient", Spacer[3], 
               Style[c, Italic]}], 11],
           Style[Dynamic@padIt2[N@zeta*2*Sqrt[mass*k0], {5, 3}], 11], 
           ""},
          
          {Text@
            Style[Row[{"magnification factor", Spacer[3], \[Beta]}], 
             11] ,
           Style[Dynamic[If[magnificationFactor === Infinity,
              magnificationFactor, 
              padIt2[N@magnificationFactor, {7, 3}]]], 11], ""},
          
          {Text@
            Style[Row[{"static displacement", 
               Spacer[3], \[CapitalDigamma], "/", Style[k, Italic]}], 
             11],
           Style[Dynamic[padIt2[N[f0/k0], {7, 3}]], 11], ""},
          
          {Text@Style[Row[{"time constant", Spacer[3], \[Tau]}], 11],
           
           Style[Dynamic[
             If[timeConstant === Infinity, Infinity, 
              padIt2[N@timeConstant, {7, 3}]]], 11], Style["sec", 11]}
          
          }, Spacings -> {.7, .4}, Frame -> {All, All}, 
         FrameStyle -> Directive[Thickness[.5], Gray], 
         Alignment -> Left
         ]
        },
       {Grid[{
          {Text[Style["test cases", Bold, 11]],
           PopupMenu[Dynamic[testCase, {testCase = #;
               Which[testCase == 1,
                k0 = 2.; mass = 4.77; forcingFrequencyInRadian = 0.6; 
                forceCritical = False; zeta = 0.;
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
                u0 = 0.; v0 = 0.; f0 = 1.;
                tscale = 500; responsePlotType = "full solution"; 
                gTick += del
                ,
                testCase == 2,
                k0 = 1.; mass = k0/0.4^2; 
                forcingFrequencyInRadian = 0.4; forceCritical = False;
                 zeta = 0.;
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0; v0 = 0.0; f0 = 1.; 
                forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
                tscale = 300; responsePlotType = "full solution"; 
                gTick += del
                ,
                testCase == 3,
                k0 = 1; mass = 1; naturalFrequency = Sqrt[k0/mass]; 
                zeta = 2*(0.01)/(2*Sqrt[mass*k0]); 
                cc = zeta*(2 mass naturalFrequency); 
                forcingFrequencyInRadian = 
                 Sqrt[k0/mass - cc^2/(2 k0 mass)]; 
                forceCritical = False;
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0.; v0 = 0.; f0 = 1; 
                forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
                tscale = 90; responsePlotType = "full solution"; 
                gTick += del
                ,
                testCase == 32,
                k0 = 1.; mass = 1; naturalFrequency = Sqrt[k0/mass]; 
                zeta = 2*(0.1)/(2*Sqrt[mass*k0]); 
                cc = zeta*(2 mass naturalFrequency); 
                forcingFrequencyInRadian = 
                 Sqrt[k0/mass - cc^2/(2 k0 mass)]; 
                forceCritical = False;
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0.; v0 = 0.; f0 = 1; tscale = 90; 
                responsePlotType = "full solution"; 
                forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
                gTick += del
                ,
                testCase == 4,
                k0 = 7.8; mass = 6.07; zeta = 4.18/(2*Sqrt[mass*k0]);
                forcingFrequencyInRadian = 0 ; forceCritical = False;
                
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0.; v0 = 0.; f0 = 1.; forcingFrequency = 0; 
                tscale = 30; responsePlotType = "full solution"; 
                gTick += del,
                
                testCase == 5,
                k0 = 7.8; mass = 2.4; forcingFrequencyInRadian = 0; 
                forceCritical = False; zeta = 4.18/(2*Sqrt[mass*k0]);
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0.; v0 = 0.; f0 = 1.; forcingFrequency = 0 ; 
                tscale = 10; responsePlotType = "full solution"; 
                gTick += del,
                
                testCase == 6,
                k0 = 7.8; mass = 6.07; forcingFrequencyInRadian = 0; 
                forceCritical = False; zeta = 10/(2*Sqrt[mass*k0]);
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 0.; v0 = 0.; f0 = 1.; forcingFrequency = 0 ; 
                tscale = 30; responsePlotType = "full solution"; 
                gTick += del,
                
                testCase == 7,
                k0 = 1; mass = 9.96; 
                forcingFrequencyInRadian = 0.317 ; 
                forceCritical = False; zeta = 5.68/(2*Sqrt[mass*k0]);
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 1.; v0 = 1.; f0 = 1.; 
                forcingFrequency = forcingFrequencyInRadian/(2*Pi); 
                tscale = 200; responsePlotType = "full solution"; 
                gTick += del,
                
                testCase == 8,
                k0 = 1; mass = 10; forcingFrequencyInRadian = 0 ; 
                forceCritical = False; zeta = .7/(2*Sqrt[mass*k0]);
                {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
                  timeConstant} = 
                 convertToStandardValues[k0, mass, 
                  forcingFrequencyInRadian, forceCritical, zeta];
                
                u0 = 1; v0 = 0; forcingFrequency = 0; r = 0; 
                tscale = 80; f0 = 0; phaseDiagramType = "argand"; 
                responsePlotType = "full solution"; gTick += del
                ];
               gTick += del} &],
            { 
             1 -> Style["beating phenomenon", 10],
             2 -> Style["resonance no damping", 10],
             3 -> Style["resonance underdamped (1)", 10],
             32 -> Style["resonance underdamped (2)", 10],
             4 -> Style["step response underdamped", 10],
             5 -> Style["step response critical", 10],
             6 -> Style["step response overdamped", 10],
             7 -> Style[Row[{"response phase lag by 90", Degree}], 10],
             8 -> Style["free underdamped response", 10]
             }, ImageSize -> All]
           }
          }, Alignment -> Left,
         Spacings -> {.1, .4}
         ]
        }
       }]
     }
    }
   ], ControlPlacement -> Left]
 , (* ABOVE IS THE LEFT PANEL*)
 Item[
  Panel[
   Grid[{
     {
      Grid[{
        {Dynamic[
          gTick;
          Show[
           dynamicMagnificationPlot,
           If[f0 == 0 || r == 0, Sequence @@ {},
            Graphics[{Red,
              PointSize[.04],
              
              Point[{r, 
                If[magnificationFactor > 5, 5, magnificationFactor]}]
              }]
            ]
           ]]}
        }](*1,1*)
      ,
      Grid[{
        {Dynamic[
          gTick;
          Show[
           phasePlot,
           If[f0 == 0 || r == 0, Sequence @@ {},
            
            Graphics[{Red, PointSize[.04], 
              Point[{r, 
                If[r == 1, -90, 180./Pi ArcTan[1 - r^2, -2 zeta r^2]]}]
              }]
            ]
           ]
          ]
         }
        }](*1,2*)
      },
     {
      Grid[{
        {
         PopupMenu[
          Dynamic[responsePlotType, {responsePlotType = #; 
             gTick += del} &],
          { 
           "transient" -> Style["transient", 10],
           "steady state" -> Style["steady state", 10],
           
           "transient+steady state (separate)" -> 
            Style["transient+steady state (separate)", 10],
           
           "full solution" -> 
            Style["transient+steady state (combined)", 10],
           "load with response" -> Style["load with response", 10]
           }]
         },
        {
         Dynamic[
          gTick;
          
          
          makeResponsePlot[responsePlotType, f0, k0, r, v0, u0, zeta, 
           naturalFrequency, tscale, transient, 
           steadyState, {206, 206}, t]
          ]
         },
        {
         Grid[{
           {Style["time", 10],
            
            Manipulator[
             Dynamic[tscale, {tscale = #; gTick += del} &], {0.1, 500,
               0.1}, ImageSize -> Tiny], 
            Style[Dynamic@padIt2[N[tscale], {4, 1}], 11], Spacer[2], 
            Style["(sec)", 10]
            }
           }, Spacings -> {.2, 0}, Alignment -> Left]
         }
        }, Spacings -> {0, 0}, Alignment -> Center, Frame -> None
       ](*2,1*)
      ,
      Grid[{
        {Dynamic[gTick; phaseLagPlot]},
        {
         RadioButtonBar[
          Dynamic[phaseDiagramType, {phaseDiagramType = #} &],
          {"argand" -> Style["Argand", 10], 
           "standard" -> Style["standard", 10]
           }, Appearance -> "Horizontal"]
         }
        }, Spacings -> {0, 0}
       ](*2,2*)
      }
     }, Spacings -> {0, 0},
    Frame -> All, FrameStyle -> Directive[Thickness[.005], LightGray]
    ], Background -> White, Alignment -> Center, FrameMargins -> 0
   ], ControlPlacement -> Right],
 
 {{gTick, 0}, None},
 {{del, $MachineEpsilon}, None},
 {{cc, 0.7}, None},
 {{u0, 1.}, None},
 {{v0, 1.}, None},
 {{f0, 1.}, None},
 {{forcingFrequencyInRadian, 0.1}, None},
 {{forcingFrequency, .1/(2*Pi)}, None},
 {{k0, 1.}, None},
 {{mass, 10.}, None},
 {{tscale, 200.}, None},
 {{forceCritical, False}, None},
 {{testCase, 1}, None},
 {{specialPlot, 1}, None},
 {{responsePlotType, "full solution"}, None},
 {{phaseDiagramType, "argand"}, None},
 
 {{zeta, 0.111}, None},
 {{naturalFrequency, 0.316}, None},
 {{dampedFrequency, 0.314}, None},
 {{r, .1}, None},(*frequency ratio*)
 {{tdd, (2 Pi)/0.314}, None},
 {{timeConstant, 3.162}, None},
 
 {{magnificationFactor, 0}, None},
 {dynamicMagnificationPlot, None},
 {phasePlot, None},
 {phaseLagPlot, None},
 {{setIC, True}, None},
 {transient, None},
 {steadyState, None},
 {sol, None},
 
 SynchronousUpdating -> True,
 AppearanceElements -> "ManipulateMenu",
 ControlPlacement -> Left,
 Alignment -> Center,
 ImageMargins -> 0,
 FrameMargins -> 1,
 ContentSize -> {0},
 SynchronousInitialization -> True,
 ContinuousAction -> True,
 Alignment -> Center,
 Paneled -> True,
 Frame -> False,
 TrackedSymbols :> {gTick},
 AutorunSequencing -> Automatic,
 
 Initialization :>
  {
   (*definitions used for parameter checking*)
   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[Chop[v],
     f, NumberSigns -> {"-", "+"}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   padIt2[v_?numeric, f_List] := AccountingForm[Chop[v],
     f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   padIt2[v_?numeric, f_Integer] := AccountingForm[Chop[v],
     f, NumberSigns -> {"", ""}, NumberPadding -> {"0", "0"}, 
     SignPadding -> True];
   (*--------------------------------------------*)
   makePhaseDifferencePlot[zeta_?numericPositive,
     f0_?numericPositive, forcingFrequency_?numericPositive, 
     naturalFrequency_?numericStrictPositive] := 
    Module[{r, phaseAngle, z, data},
     
     If[f0 == 0, r = 0, r = forcingFrequency/naturalFrequency];
     phaseAngle = If[r == 1, Pi/2, ArcTan[1 - r^2, 2 zeta r]];
     
     data = 
      Table[{{z, If[f0 == 0, 0, Sin[z - phaseAngle]]}, {z, 
         If[f0 == 0, 0, Sin[z]]}}, {z, 0, 2 Pi, Pi/10}];
     
     ListLinePlot[{data[[All, 1]], data[[All, 2]]},
      PlotLabel -> Style[Grid[{
          {"displacement phase relative to load"},
          {"red color represents load"},
          {Row[{\[Theta], " = ", 
             padIt1[N[-(180/Pi) phaseAngle], {4, 1}], Degree}]
           }},
         Spacings -> {0, .2}], 11],
      PlotStyle -> {{Dashed, Black}, Red}, PlotRange -> All, 
      ImagePadding -> {{10, 12}, {10, 10}}, GridLines -> Automatic, 
      GridLinesStyle -> Directive[Thickness[.001], LightGray], 
      AspectRatio -> 0.9, 
      Ticks -> {{0, Pi/2, Pi, Pi + 1/2 Pi, 2 Pi}, None}
      ]
     ];
   (*--------------------------------------------*)
   makeArgandDiagram[
     zeta_?numericPositive,
     f0_?numericPositive,
     forcingFrequency_?numericPositive,
     naturalFrequency_?numericStrictPositive,
     maxTime_?numericStrictPositive] := 
    Module[{phaseAngle, r, tipOfForce},
     
     If[f0 == 0, r = 0, r = forcingFrequency/naturalFrequency];
     phaseAngle = If[r == 1, -(Pi/2), ArcTan[1 - r^2, -2 zeta r]];
     tipOfForce = {Cos[forcingFrequency  maxTime - Pi/2], 
       Sin[forcingFrequency  maxTime - Pi/2]};
     ListLinePlot[{{{0, 0}}, {{0, 0}}},
      PlotLabel -> Style[
        
        Grid[{
          {"displacement phase relative to load"},
          {"red color represents load"},
          {Row[{\[Theta], " = ", 
             padIt1[N[180/Pi (phaseAngle)], {4, 1}], Degree}]
           }
          },
         Spacings -> {0, .2}
         ], 11],
      
      AspectRatio -> .9,
      ImagePadding -> {{10, 12}, {10, 10}},
      GridLines -> {Range[-1.2, 1.2, .2], Range[-1.2, 1.2, .2]}, 
      GridLinesStyle -> Directive[LightGray, Thickness[.001]],
      PlotRange -> {{-1.2, 1.2}, {-1.2, 1.2}},
      Ticks -> None,
      Axes -> True,
      Epilog -> {
        {Red, Thick, Arrowheads[.1], Arrow[{{0, 0}, tipOfForce}]},
        {Blue, Arrowheads[.1], 
         Arrow[{{0, 
            0}, {Cos[forcingFrequency  maxTime + phaseAngle - Pi/2], 
            Sin[forcingFrequency  maxTime + phaseAngle - Pi/2]}}]},
        {Black, Circle[{0, 0}, 1]},
        If[phaseAngle < -Pi/16,
         circularArrow[
          Circle[{0, 0}, 
           0.3, {forcingFrequency  maxTime - Pi/2, 
            forcingFrequency  maxTime + phaseAngle - Pi/2}]],
         Sequence @@ {}]
        }
      ]
     ];
   (*--------------------------------------------*)
   (*Function below by belisarius from stackoverflow used in the \
above function to draw an arced arrow*)
   circularArrow[s_Circle] := s /. Circle[a_, r_, {start_, end_}] :> (
       {
          {Blue, s},
          {Blue, Arrow[{# - r/10^6 {Sin@end, -Cos@end}, #}]}
          } &[a + r {Cos@end, Sin@end}]
       );
   (*--------------------------------------------*)
   convertToStandardValues[k0_?numericStrictPositive,
     mass_?numericStrictPositive,
     forcingFrequency_?numericPositive,
     isCriticalDamping_?bool,
     zetaOld_?numeric] := Module[{naturalFrequency,
      r, cc, dampedFrequency, tdd, timeConstant, zeta = zetaOld},
     
     naturalFrequency = Sqrt[k0/mass];
     r = forcingFrequency/naturalFrequency;
     cc = zeta*(2 mass naturalFrequency);
     
     If[isCriticalDamping || Abs[zeta - 1] <= $MachineEpsilon,
      zeta = 1;
      dampedFrequency = 0;
      tdd = Infinity;
      timeConstant = 1/naturalFrequency;
      cc = 2 mass naturalFrequency
      ,
      If[zeta - 1 > $MachineEpsilon,
       (*overdamped*)
       dampedFrequency = 0;
       tdd = Infinity;
       timeConstant = 1.0/(naturalFrequency (zeta - Sqrt[zeta^2 - 1]));
       cc = zeta (2 mass naturalFrequency)
       ,
       (*must be underdamped or zero*)
       dampedFrequency = naturalFrequency Sqrt[1.0 - zeta^2];
       tdd = (2 \[Pi])/dampedFrequency;
       timeConstant = 1/naturalFrequency;
       cc = zeta (2 mass naturalFrequency)]
      ];
     
     {naturalFrequency, r, cc, zeta, dampedFrequency, tdd, 
      timeConstant}
     ];
   
   (*--------------------------------------------*)
   makeGenericPhasePlot[zeta_List, actualR_] := 
    Module[{i, r, data, legend, opt = {Right, Top}},
     r = If[actualR < 2, 2, actualR];
     
     data = Table[{i,
          
          If[i == 1 || # == 0, -90, 
           180./Pi ArcTan[1 - i^2, -2 # i^2]]}, {i, 0, r, .01}
         ] & /@ zeta;
     
     legend = 
      Row[{\[Xi], " = ", Style[Dynamic@padIt2[#, {3, 2}], 11]}] & /@ 
       zeta;
     
     ListLinePlot[data,
      PlotRange -> All,
      PlotLegends -> Placed[LineLegend[legend,
         LegendMargins -> 1,
         LegendLayout -> Function[{x},
           Grid[x, Spacings -> {0.1, 0}, Alignment -> Left]]], opt
        ],
      GridLines -> Automatic,
      GridLinesStyle -> Directive[LightGray, Thickness[.001]],
      Frame -> True,
      FrameLabel -> {{\[Theta], 
         None}, {Row[{"frequency ratio", Spacer[5], \[CurlyPi], 
           " / ", \[Omega]}], 
         Style["phase angle vs. frequency ratio", 12]}},
      ImageMargins -> 0,
      ImagePadding -> {{20, 12}, {35, 25}},
      AspectRatio -> 1.05,
      RotateLabel -> False,
      FrameTicksStyle -> 8,
      AxesStyle -> {Dashed, Gray},
      FrameTicks -> {{{-180, -135, -90, -45, 0}, None},
        {{0, 0.5, 1, 1.5, 2, 2.5, 3}, None}},
      PerformanceGoal -> "Speed"
      ]
     ];
   
   (*--------------------------------------------*)
   makeGenericDynamicMagnificationFactorPlot[zeta_List, actualR_] := 
    Module[{r, i,
      data, legend, opt = {Right, Top}, z},
     
     r = If[actualR < 2, 2, actualR];
     
     data = 
      Table[{i, 
          If[(# == 0 && i == 1), 
           5, (z = 1/Sqrt[(1 - i^2)^2 + (2 # i)^2]; 
            If[z > 5, 5, z])]}, {i, 0, r, .01}] & /@ zeta;
     
     legend = 
      Row[{\[Xi], " = ", Style[Dynamic@padIt2[#, {3, 2}], 11]}] & /@ 
       zeta;
     
     ListLinePlot[data,
      PlotRange -> All,
      PlotLegends -> Placed[LineLegend[legend,
         LegendMargins -> 1,
         LegendLayout -> Function[{x}, Grid[x, Spacings -> {0.1, 0},
            Alignment -> Left]]], opt],
      GridLines -> Automatic,
      GridLinesStyle -> Directive[LightGray, Thickness[.001]],
      Frame -> True,
      FrameLabel -> {{\[Beta], 
         None}, {Row[{"frequency ratio", Spacer[5], \[CurlyPi], 
           " / ", \[Omega]}], 
         Style["dynamic magnification factor", 12]}},
      ImagePadding -> {{25, 12}, {35, 25}},
      ImageMargins -> 0,
      AspectRatio -> 0.88,
      ImageSize -> {206, 206},
      RotateLabel -> False,
      PerformanceGoal -> "Speed"]
     ];
   (*--------------------------------------------*)
   calculateMagnificationFactor[zeta_?numericPositive,
     f_?numericPositive,
     forcingFrequency_?numericPositive,
     naturalFrequency_?numericStrictPositive] := 
    Module[{r = forcingFrequency/naturalFrequency},
     Which[f == 0, 0,
      zeta == 0, If[r == 1, Infinity, 1/(1 - r^2)],
      True, 1/Sqrt[(1 - r^2)^2 + (2 zeta r)^2]
      ]
     ];
   
   (*--------------------------------------------*)
   makeResponsePlot[responsePlotType_String, f0_, k_?numericPositive, 
     r_?numericPositive, v0_, u0_, zeta_, naturalFrequency_, tscale_, 
     transient_, steadyState_, imageSize_List, t_] := 
    Module[{sol, plotLegend, plotStyle, plot, envelop, 
      dampedFrequency, logarithmicDecrement, uAtEndOfCycle, 
      epilog = {}},
     
     Which[responsePlotType == "transient",
      sol = transient;
      plotStyle = Red,
      
      responsePlotType == "steady state",
      sol = steadyState;
      plotStyle = Blue,
      
      responsePlotType == "transient+steady state (separate)",
      sol = {transient, steadyState};
      plotStyle = {Red, Blue},
      
      responsePlotType == "full solution",
      sol = transient + steadyState;
      plotStyle = Blue,
      
      responsePlotType == "load with response",
      
      sol = {transient + steadyState,
        If[r == 0,
         f0/k,
         f0/k Sin[naturalFrequency*r t] (*must be harmonic loading*)
         ]
        };
      
      plotStyle = {Blue, Red}
      
      ];
     (*make envelope*)
     If[f0 == 0 && (u0 != 0 || v0 != 0) && 
       responsePlotType != "steady state" && zeta > 0 && zeta < 1,
      (
       plotLegend = None;
       dampedFrequency = naturalFrequency*Sqrt[1 - zeta^2];
       uAtEndOfCycle = (transient + steadyState) /. 
         t -> (2 Pi/dampedFrequency);
       logarithmicDecrement = Log[u0/uAtEndOfCycle];
       epilog = 
        Text[Style[
          Row[{"logrithmic decrement", Spacer[1], 
            padIt2[logarithmicDecrement, {4, 2}]}], 11],
         Scaled[{0.5, 0.9}], {0, 0}];
       envelop = 
        Sqrt[u0^2 + ((v0 + zeta*naturalFrequency*u0)/
             dampedFrequency)^2];
       plot = Plot[{envelop Exp[-zeta naturalFrequency t],
          -envelop Exp[-zeta naturalFrequency tt]}, {t, 0, tscale},
         PlotStyle -> {{Magenta, Dashed}, {Magenta, Dashed}},
         PlotRange -> All,
         PerformanceGoal -> "Speed"
         ]
       ),
      (
       plot = {};
       If[f0 > 0,
        (
         If[responsePlotType == "transient+steady state (separate)",
          (
           plotLegend = Placed[LineLegend[
              Style[#, 10] & /@ {"transient", "steady state"},
              LegendMargins -> 1,
              
              LegendLayout -> 
               Function[{x}, 
                Grid[x, Spacings -> {0.1, 0}, 
                 Alignment -> Left]]], {Center, Top}]
           ),
          (
           If[responsePlotType == "load with response",
            (
             
             plotLegend = 
              Placed[LineLegend[Style[#, 10] & /@ {"response", "load"},
                LegendMargins -> 1, 
                LegendLayout -> 
                 Function[{x}, 
                  Grid[x, Spacings -> {0.1, 0}, 
                   Alignment -> Left]]], {Center, Top}]
             ),
            (
             plotLegend = None
             )
            ]
           )
          ]
         ),
        (
         plotLegend = None
         )
        ]
       )
      ];
     
     Show[
      Plot[Tooltip[Chop@Evaluate@sol, 
        Text[Style[TraditionalForm[Chop@sol], 14]]], {t, 0, tscale},
       PlotRange -> {{0, tscale}, All},
       ImagePadding -> {{30, 15}, {10, 20}},
       ImageMargins -> {{0, 0}, {0, 0}},
       FrameLabel -> {
         {None, None},
         {None, 
          Row[{"system response ", Style["u", Italic], "(", 
            Style["t", Italic], ")", " vs. time"}]
          }},
       Frame -> True,
       LabelStyle -> 12,
       RotateLabel -> False,
       GridLines -> Automatic,
       GridLinesStyle -> Directive[Thickness[.001], LightGray],
       AxesOrigin -> {0, 0}, AxesStyle -> {Dashed, Gray},
       FrameTicksStyle -> 8,
       AspectRatio -> 1.05,
       ImageSize -> imageSize,
       PlotStyle -> plotStyle,
       PlotLegends -> plotLegend,
       PerformanceGoal -> "Speed",
       Exclusions -> None,
       Epilog -> epilog,
       Evaluated -> True
       ],
      plot
      ]
     ];
   
   (*--------------------------------------------*)
   oneDegreeOfFreedomSolution[u0_?numeric, v0_?numeric, 
     k_?numericStrictPositive, zeta_?numericPositive, 
     f0_?numericPositive, forcingFrquency_?numericPositive, 
     w_?numericStrictPositive, t_] := 
    Module[{wd, z1, z2, a, b, r, phaseAngle, p1, p2, steadyState, 
      transient},
     
     r = forcingFrquency/w;
     
     If[f0 == 0,
      steadyState = 0;
      transient = Which[
        zeta == 0, u0 Cos[w t] + v0/w Cos[w t],
        
        zeta < 1,
        wd = w Sqrt[1 - zeta^2];
        Exp[-zeta w t ] (u0 Cos[wd t] + (v0 + u0 zeta w)/wd Sin[wd t]),
        
        zeta == 1,
        (u0 + (v0 + u0 w) t) E^(-w t),
        
        True,
        p1 = -w zeta + w Sqrt[zeta^2 - 1];
        p2 = -w zeta - w Sqrt[zeta^2 - 1];
        a = (v0 - u0 p2)/(2 w Sqrt[zeta^2 - 1]);
        b = (-v0 + u0 p1)/(2 w Sqrt[zeta^2 - 1]);
        Exp[p1 t] + b Exp[p2 t]
        ],
      (*F is not zero now*)
      {steadyState, transient} = Which[
        zeta == 0,(*undamped*)
        If[forcingFrquency == 0,(*constant force*)
         {f0/k, (u0 - f0/k) Cos[w t] + v0/w Sin[w t]}
         ,
         If[Abs[r - 1] <= 0.01,(*resonance*)
          {-(f0/k) (forcingFrquency t)/2 Cos[forcingFrquency t],
           
           u0 Cos[forcingFrquency t] + 
            v0/forcingFrquency Sin[forcingFrquency t]
           },
          {(f0/k) 1/(1 - r^2) Sin[forcingFrquency t],
           u0 Cos[w t] + (v0/w - ((f0/k) r/(1 - r^2))) Sin[w t]}]
         ],
        
        zeta < 1,
        phaseAngle = If[r == 1, Pi/2, ArcTan[1 - r^2, 2 zeta r]];
        z1 = Sqrt[(1 - r^2)^2 + (2 zeta r)^2];
        wd = w Sqrt[1 - zeta^2];
        a = If[forcingFrquency == 0,
          u0 - (f0/k),
          u0 + (f0/k)/z1 Sin[phaseAngle]
          ];
        b = If[forcingFrquency == 0, (v0 + a zeta w)/wd,
          
          v0/wd + (u0 zeta w)/
            wd + (f0/k)/(wd z1) (zeta w Sin[phaseAngle] - 
              forcingFrquency Cos[phaseAngle])
          ];
        
        If[forcingFrquency == 0,
         {(f0/k), Exp[-zeta w t]*(a Cos[wd t] + b Sin[wd t])}
         ,
         {(f0/k)/z1 Sin[forcingFrquency t - phaseAngle], 
          Exp[-zeta w t] (a Cos[wd t] + b Sin[wd t])}],
        
        (*critical*)
        zeta == 1,
        phaseAngle = If[r == 1, Pi/2, ArcTan[1 - r^2, 2 r]];
        z1 = Sqrt[(1 - r^2)^2 + (2 r)^2];
        a = Which[forcingFrquency == 0,
          u0 - (f0/k),
          True,
          u0 + (f0/k)/z1 Sin[phaseAngle]
          ];
        b = Which[forcingFrquency == 0,
          v0 + u0 w - (f0/k) w,
          True,
          
          v0 + u0 w + (f0/k)/
             z1 (w Sin[phaseAngle] - forcingFrquency Cos[phaseAngle])
          ];
        
        If[forcingFrquency == 0,
         {(f0/k), (a + b t) Exp[-w t]}
         ,
         {(f0/k)/z1 Sin[
            forcingFrquency t - phaseAngle], (a + b t) Exp[-w t]}
         ],
        
        (*overdamped*)True,
        phaseAngle = If[r == 1, Pi/2, ArcTan[1 - r^2, 2 zeta r]];
        z1 = (f0/k)/Sqrt[(1 - r^2)^2 + (2 zeta r)^2];
        z2 = Sqrt[zeta^2 - 1];
        a = (v0 + u0 w zeta + u0 w z2 + 
            
            z2 (w (zeta + z2) Sin[phaseAngle] - 
               forcingFrquency Cos[phaseAngle]))/(2 w z2);
        b = -((v0 + u0 w zeta - u0 w z2 + 
              z2 (w (zeta - z2) Sin[phaseAngle] - 
                 forcingFrquency Cos[phaseAngle]))/(2 w z2));
        
        If[forcingFrquency == 0,
         p1 = -w zeta + w Sqrt[zeta^2 - 1];
         p2 = -w zeta - w Sqrt[zeta^2 - 1];
         b = ((f0/k) p1 - u0 p1 + v0)/(p2 - p1); a = u0 - (f0/k) - b;
         {
          (f0/k),
          a Exp[p1 t] + b Exp[p2 t]
          }
         ,
         {z1 Sin[forcingFrquency t - phaseAngle],
          a Exp[(-zeta + z2) w t] + b Exp[(-zeta - z2) w t]
          }]
        ]
      ];
     {transient, steadyState}
     ];
   }
 ]