(* by Nasser M. Abbasi, version: April 26, 2013 *)
Manipulate[
 (
  (*If Q, the quality factor is zero, 
  then we assume no damping will be used*)
  (*This way there is no need to add an extra dialog asking if \
damping is to be used or not*)
  isDamped = If[Abs@q <= $MachineEpsilon, False, True];
  
  (* start of the finite state machine, 
  where all the logic of the program in included. 5 states, 
  8 events *)
  (* events are set using the second argument of dynamics right at \
the control level below *)
  (* this is a way to simulate event driven UI with callbacks *)
  
  Which[
   (*-------  STATE  RESET ---------------*)
   state == "reset", 
   (
    currentTime = 0;
    tick = 0;
    {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, x1,
       x2, x3, currentRunMaxSpeed} = 
     initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
      isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
      spectrumScale, currentTime, dt, maxRunTime];
    
    If[bifurcationPlot === {}, 
     bifurcationPlot = 
      makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
       x10, x20, isDamped, plotImageSize, 
       Unevaluated@bifurcationProgress]
     ];
    
    state = "pause"
    ),
   (*-------  STATE  PAUSE ---------------*)
   state == "pause",
   (
    Which[
     lastEvent == "no_event", state = "pause",
     lastEvent == "pause_button", (lastEvent = "no_event";  
      state = "pause"),
     lastEvent == "reset_button", (lastEvent = "no_event"; 
      state = "reset"; tick += DEL),
     lastEvent == 
      "initial_conditions_changed", (lastEvent = "no_event"; 
      state = "reset"; tick += DEL),
     lastEvent == "time_series_choice", (lastEvent = "no_event";
      bottomPlot = If[plotChoice == 1,
        makeTimeSeriesPlot[phaseData, plotImageSize, nSteps, 
         currentTime],
        makePowerSpectrumPlot[phaseData, plotImageSize, nSteps, omega,
          spectrumScale, currentTime]
        ]),
     
     lastEvent == "step_button",
     (
      lastEvent = "no_event";
      {currentTime, nSteps, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3} = 
       makeOneStep[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, nSteps, dt, phasePlot, bottomPlot,
         track, maxRunTime, x1, x2, x3]
      ),
     
     lastEvent == "run_button", (lastEvent = "no_event"; 
      state = "running"; tick += DEL),
     lastEvent == "fast_button", (lastEvent = "no_event"; 
      state = "fast mode"; tick += DEL),
     lastEvent == "update_bifurcation_button",
     (lastEvent = "no_event";
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress]
      ),
     
     lastEvent == "test_case_changed",
     (
      lastEvent = "no_event";
      {torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals, 
        x10, isDamped, x20} = setParameters[testCase];
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress];
      
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime]
      )
     ]
    ),
   (*-------  STATE  RUNNING ---------------*)
   state == "running",
   (
    If[poincareMap === {},
     (
      currentTime = 0;
      
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime]
      )
     ];
    
    If[bifurcationPlot === 0,
     (
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress]
      )
     ];
    
    Which[
     lastEvent == "initial_conditions_changed",
     (
      lastEvent = "no_event";
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime];
      
      tick += DEL
      ),
     
     lastEvent == "no_event" || lastEvent == "time_series_choice",
     (
      If[lastEvent == "time_series_choice", lastEvent = "no_event"];
      
      If[currentTime + dt <=  maxRunTime,
       (
        {currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
           x1, x2, x3} = 
         makeOneStep[q, torqueAmplitude, omega, phase, x10, x20, 
          isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
          spectrumScale, currentTime, nSteps, dt, phasePlot, 
          bottomPlot, track, maxRunTime, x1, x2, x3];
        
        tick += DEL
        ),
       (
        state = "pause"
        )
       ]
      ),
     
     lastEvent == "step_button",
     (
      lastEvent = "no_event";
      If[currentTime + dt <=  maxRunTime,
       (
        {currentTime, nSteps, phasePlot, bottomPlot, track, phaseData,
           x1, x2, x3} = 
         makeOneStep[q, torqueAmplitude, omega, phase, x10, x20, 
          isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
          spectrumScale, currentTime, nSteps, dt, phasePlot, 
          bottomPlot, track, maxRunTime, x1, x2, x3];
        
        state = "pause"
        )
       ]
      ),
     lastEvent == "pause_button", (lastEvent = "no_event"; 
      state = "pause"),
     lastEvent == "fast_button", (lastEvent = "no_event"; 
      state = "fast mode" ; tick += DEL),
     lastEvent == "reset_button", (lastEvent = "no_event"; 
      state = "reset"; tick += DEL),
     
     lastEvent == "update_bifurcation_button",
     (
      lastEvent = "no_event";
       bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress];
      
      tick += DEL
      ),
     
     lastEvent == "test_case_changed",
     (
      lastEvent = "no_event";
      {torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals, 
        x10, isDamped, x20} = setParameters[testCase];
      
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress];
      
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime];
      
      tick += DEL
      )
     ]
    ),
   (*-------  STATE  FAST_MODE ---------------*)
   state == "fast mode",
   (
    If[poincareMap === {},
     (
      currentTime = 0;
      
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime]
      )
     ];
    
    If[bifurcationPlot === {},
     bifurcationPlot = 
      makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
       x10, x20, isDamped, plotImageSize, 
       Unevaluated@bifurcationProgress]
     ];
    
    Which[
     lastEvent == "no_event" || 
      lastEvent == "initial_conditions_changed" || 
      lastEvent == "time_series_choice" || lastEvent == "fast_button",
     (
      If[Not[lastEvent == "no_event"], lastEvent = "no_event"];
      
      currentTime = maxRunTime;
      {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
        x1, x2, x3, currentRunMaxSpeed} = 
       initializeSystem[q, torqueAmplitude, omega, phase, x10, x20, 
        isDamped, plotImageSize, axesLabel, plotChoice, phaseData, 
        spectrumScale, currentTime, dt, maxRunTime]
      ),
     lastEvent == "pause_button", (lastEvent = "no_event"; 
      state = "pause"),
     lastEvent == "reset_button", (lastEvent = "no_event"; 
      state = "reset"; tick += DEL),
     lastEvent == "run_button", (lastEvent = "no_event"),
     lastEvent == "step_button", (lastEvent = "no_event"),
     
     lastEvent == "update_bifurcation_button",
     (lastEvent = "no_event";
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress]
      ),
     
     lastEvent == "test_case_changed",
     (
      lastEvent = "no_event";
      {torqueAmplitude, q, omega, phase, aStart, aLen, aIntervals, 
        x10, isDamped, x20} = setParameters[testCase];
      
      bifurcationPlot = 
       makeBifurcationPlot[aStart, aLen, aIntervals, q, omega, phase, 
        x10, x20, isDamped, plotImageSize, 
        Unevaluated@bifurcationProgress];
      
      tick += DEL
      )
     ]
    )
   ];
  
  generateFinalDisplay[poincareMap, phasePlot, bifurcationPlot, track,
    currentTime, plotImageSize, x1, x2, x3, bottomPlot, 
   currentRunMaxSpeed, torqueAmplitude]
  
  ),
 (*--- START OF CONTROLS ---*)
 Grid[{
   {Grid[{(* BLOCK 1 *)
      
      {Grid[{
         {Button[
           Text[Style["play", 14]], (lastEvent = "run_button"; 
            tick += DEL), ImageSize -> {75, 30}],
          
          Button[Text[
            Style["pause", 14]], (lastEvent = "pause_button"; 
            tick += DEL), ImageSize -> {75, 30}],
          
          Button[Text[Style["step", 14]], (lastEvent = "step_button"; 
            tick += DEL), ImageSize -> {75, 30}]}
         }, Spacings -> {.2, .1}, Alignment -> Left]
       },
      
       {Grid[{
         {Button[
           Text[Style["reset", 14]], (lastEvent = "reset_button"; 
            tick += DEL), ImageSize -> {75, 30}],
          
          Button[Text[Style["fast", 14]], (lastEvent = "fast_button"; 
            tick += DEL), ImageSize -> {75, 30}],
          
          
          Framed[Dynamic@
            Graphics[Text@Style[state, 12, Italic], 
             ImageSize -> {75, 30}], 
           FrameStyle -> Directive[Thickness[.005], Gray], 
           FrameMargins -> -1, BaselinePosition -> Baseline, 
           ImageMargins -> 0, 
           ContentPadding -> True]} (*the state of the system*)
         
         }, Spacings -> {.2, .1}, Alignment -> Left]}
      
      }, Spacings -> {.2, .1}, Alignment -> Center
     ]
    },
   
   {Grid[{(* BLOCK 2 *)
      {Text@Style["maximum duration", 11], Spacer[4],
       Manipulator[
        Dynamic[maxRunTime, (maxRunTime = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {1, 
         500, 1}, ImageSize -> Tiny, ContinuousAction -> False], 
       Spacer[4], Dynamic@Text@Style[padIt3[maxRunTime, {3, 0}], 11]
       },
      {Text@
        Style[TraditionalForm[HoldForm[\[CapitalDelta]\[Tau]]], 11], 
       Spacer[4], 
       Manipulator[
        Dynamic[dt, (dt = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {0.01,
          1, 0.01}, ImageSize -> Tiny, ContinuousAction -> False], 
       Spacer[4], Dynamic@Text@Style[padIt2[dt, {3, 3}], 11], 
       SpanFromLeft}
      }, Frame -> None, Spacings -> {0, 0.1}, Alignment -> Center
     ]
    },
   
   {Grid[{(* BLOCK 3 *)
      {Row[{Style["\[Theta]''(\[Tau]) +  \[Theta]'(\[Tau])/", "TR", 
          12], Style["Q", Italic, "TR", 12], 
         Style["+ sin \[Theta](\[Tau])) \[Equal] ", "TR", 12] , 
         Style["a", Italic, "TR", 12], 
         Style[" cos \[Phi](\[Tau])", "TR", 12]}]}
      }]
    },
   
   {Grid[{(* BLOCK 4 *)
      {Style["Q", Italic, "TR"], Spacer[2], 
       Manipulator[
        Dynamic[q, (q = #; {lastEvent = "initial_conditions_changed", 
            tick += DEL}; #) &], {0, 4, 0.01}, ImageSize -> Tiny, 
        ContinuousAction -> False], Spacer[2], 
       Dynamic@Text@Style[padIt2[q, {3, 3}], 12], "", ""},
      { Text@Style["\[Omega] = d\[Phi]/d\[Tau]", 12], Spacer[2], 
       Manipulator[
        Dynamic[omega, (omega = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {-4, 
         4, 0.01}, ImageSize -> Tiny, ContinuousAction -> False], 
       Spacer[2], Dynamic@Text@Style[padIt1[omega, {3, 3}], 11], "", 
       ""},
      {Text@Style["a", Italic, 12], Spacer[2], 
       Manipulator[
        Dynamic[torqueAmplitude, (torqueAmplitude = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {0, 2,
          0.01}, ImageSize -> Tiny, ContinuousAction -> False], 
       Spacer[2], 
       Dynamic@Text@Style[padIt2[torqueAmplitude, {3, 3}], 11], "", ""}
      }, Frame -> None, Spacings -> {0.2, 0.6}, Alignment -> Center]
    },
   
   {Grid[{(* BLOCK 5 *)
      {Item[Text@Style["initial conditions", 12], 
        Alignment -> Center], SpanFromLeft},
      { Text@Style["\[Theta] (0)", 12], Spacer[2], 
       Manipulator[
        Dynamic[x10, (x10 = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {-N@
           Pi, N@Pi, 0.01}, ImageSize -> Tiny, 
        ContinuousAction -> False], Spacer[2], 
       Dynamic@Text@Style[padIt1[x10, {6, 4}], 11], Spacer[1], 
       Text@Style["radian", 10]},
      {Text@Style["d\[Theta]/d\[Tau] (0)", 12], Spacer[2], 
       Manipulator[
        Dynamic[x20, (x20 = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {-Pi, 
         Pi, 0.01}, ImageSize -> Tiny, ContinuousAction -> False], 
       Spacer[2], Dynamic@Text@Style[padIt1[x20, {6, 4}], 11], "", ""},
      {Text@Style["\[Phi] (0)", 12], Spacer[2], 
       Manipulator[
        Dynamic[phase, (phase = #; {lastEvent = 
             "initial_conditions_changed", tick += DEL}; #) &], {-N@
           Pi, N@Pi, 0.01}, ImageSize -> Tiny, 
        ContinuousAction -> False], Spacer[2], 
       Dynamic@Text@Style[padIt1[phase, {5, 3}], 11], Spacer[1], 
       Text@Style["radian", 11]}
      }, Spacings -> {.4, 0.2}, Alignment -> Left
     ]},
   
   {Grid[{(* BLOCK 6 *)
      {Item[
        Text@Row[{Style["select range of parameter ", 12], 
           TraditionalForm[HoldForm[a]]}], Alignment -> Center], 
       SpanFromLeft}, {Item[Text@Style["from", 11], 
        Alignment -> Left], 
       Control[{{aStart, 0.9, ""}, 0, 2, 0.01, ImageSize -> Tiny, 
         ImagePadding -> 0}], 
       Dynamic@Text@Style[padIt2[aStart, {3, 3}], 11]},
      {Item[Text@Style["length", 11], Alignment -> Left], 
       Control[{{aLen, 0.6, ""}, 0.01, 2, 0.01, ImageSize -> Tiny, 
         ImagePadding -> 0}], 
       Dynamic@Text@Style[padIt2[aLen, {3, 3}], 11]},
      {Item[Text@Style["intervals", 11], Alignment -> Left], 
       Control[{{aIntervals, 20, ""}, 5, 300, 1, ImageSize -> Tiny, 
         ImagePadding -> 0}], 
       Dynamic@Text@Style[padIt3[Rationalize@aIntervals, {4, 0}], 11]},
      {Button[
        Text[Style["generate", 11]], {lastEvent = 
          "update_bifurcation_button"; tick += DEL; 
         Method -> "Queued"}, ImageSize -> {65, 25}, 
        BaselinePosition -> Bottom], 
       Dynamic@ProgressIndicator[bifurcationProgress, {0, 1}, 
         ImageSize -> {105, 23}, ImageMargins -> 0, 
         BaselinePosition -> Bottom], 
       Dynamic@Graphics[
         Text@Row[{ 
            Style[padIt2[bifurcationProgress*100, {4, 2}] , 11], 
            Style["%", 11]}], ImageSize -> {40, 20}, 
         BaselinePosition -> Bottom]}
      }, Frame -> None, Spacings -> {.4, 0.1}, Alignment -> Center]
    },
   
   {Grid[{(* BLOCK 8 *)
      {Text@Style["preconfigured test cases", 10]}, {PopupMenu[
        Dynamic[
         testCase, (testCase = #; {lastEvent = "test_case_changed" , 
            tick += DEL}; #) &],
        { 
         "double periodic #1" -> 
          Style["double periodic, on edge of chaotic", 10],
          
         "periodic, A=0.9,Q=2,omega=2/3" -> 
          Style["periodic, A=0.9,Q=2,omega=2/3", 10],
         "periodic, A=1.35,Q=2,omega=2/3" -> 
          Style["periodic, A=0.9,Q=2,omega=2/3", 10],
         "period doubling, A=1.12,Q=2,omega=2/3" -> 
          Style["period doubling, A=1.12,Q=2,omega=2/3", 10],
         "period doubling, A=1.07,Q=2,omega=2/3" -> 
          Style["period doubling, A=1.07,Q=2,omega=2/3", 10],
         "period doubling, A=1.45,Q=2,omega=2/3" -> 
          Style["period doubling, A=1.45,Q=2,omega=2/3", 10],
         "period quadrupling, A=1.47,Q=2,omega=2/3" -> 
          Style["period quadrupling, A=1.47,Q=2,omega=2/3", 10],
         "chaotic, A=1.15,Q=2,omega=2/3" -> 
          Style["chaotic, A=1.15,Q=2,omega=2/3", 10],
         "chaotic, A=1.5,Q=2,omega=2/3" -> 
          Style["chaotic, A=1.5,Q=2,omega=2/3", 10],
         "chaotic, A=1.5,Q=2,omega=2/3,phase=3.14" -> 
          Style["chaotic, A=1.5,Q=2,omega=2/3,phase=3.14", 10],
         "chaotic, Bifurcation #1" -> 
          Style["chaotic, Bifurcation #1", 10],
          
         "chaotic, Bifurcation #2" -> 
          Style["chaotic, Bifurcation #2", 10],
         "chaotic, Bifurcation #3" -> 
          Style["chaotic, Bifurcation #3", 10]}, ImageSize -> All, 
        ContinuousAction -> False]
       }
      }
     ]
    },
   
   {Grid[{(* BLOCK 9 *)
      {Text@Style["plot type     ", 10], Spacer[5], 
       RadioButtonBar[
        Dynamic[plotChoice, (plotChoice = #; {lastEvent = 
             "time_series_choice", tick += DEL}; #) &], {1 -> 
          Text@Style["time series", 11], 
         2 -> Text@Style["power spectrum", 11]}]}, {Text@
        Style["spectrum scale", 10], Spacer[5], 
       RadioButtonBar[
        Dynamic[spectrumScale, (spectrumScale = #; {lastEvent = 
             "time_series_choice", tick += DEL}; #) &], {1 -> 
          Text@Style["standard   ", 11], 2 -> Text@Style["log", 11]}, 
        Enabled -> Dynamic[plotChoice == 2]]}
      }, Frame -> None, Spacings -> {0.1, 0.2}, Alignment -> Left]
    }
   }, Frame -> All, FrameStyle -> Directive[Thickness[.005], Gray], 
  Spacings -> {0.2, 0.6}, Alignment -> Center
  ],
 
 {{DEL, $MachineEpsilon}, None},
 {{currentRunMaxSpeed, 0}, None},
 {{tick, 0}, None},
 {{currentTime, 300}, None},
 {{plotChoice, 1}, None},
 {{spectrumScale, 1}, None},
 {{bifurcationProgress, 0}, None},
 {{omega, 0.6667}, None},
 {{q, 2}, None},
 {{torqueAmplitude, 1.5}, None},
 {{phase, 0}, None},
 {{x10, 0.6184}, None},
 {{x20, 0}, None},
 {{dt, 0.15}, None},
 {{x1, 0}, None},
 {{x2, 0}, None},
 {{x3, 0}, None},
 {{maxRunTime, 300}, None},
 {{isDamped, True}, None},
 {{testCase, "periodic, A=0.9,Q=2,omega=2/3"}, None},
 {{state, "fast mode"}, None},
 {{lastEvent, "no_event"}, None},
 {{track, {Blue, PointSize[0.03], Point[ {0, 0}]}}, None},
 {{phaseData , Table[{0, 0, 0}, {500*101}]}, 
  None}, (*DO NOT CHANGE, DEPENDS ON FIXED SETTINGS *)
 {{poincareMap, {}}, None},
 {{phasePlot, {}}, None},
 {{bottomPlot, {}}, None},
 {{nSteps, 0}, None},
 {{plotImageSize, {340, 340}}, None},
 {{axesLabel, {Text@Style["\[Theta]", 9], 
    Text@Style[
      Row[{" ", Style["d", Italic], "\[Theta]/", Style["d", Italic], 
        "\[Tau]"}], 9]}}, None},
 {{bifurcationPlot, {}}, None},
 
 TrackedSymbols -> {tick}, (*has to be -> and NOT :> else won't work *)

 
 SynchronousUpdating -> False,
 ContinuousAction -> False,
 SynchronousInitialization -> False,
 ControlPlacement -> Left,
 Alignment -> Center, ImageMargins -> 1, FrameMargins -> 1,
 Initialization :> 
  (
   (*----------------------------------------------------------------------*)
   padIt1[v_?(NumericQ[#] && Im[#] == 0 &), f_List] := 
    AccountingForm[Chop[N@v] , f, NumberSigns -> {"-", "+"}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True];
   
   (*-----------------------------------------------------------------------*)
   padIt2[v_?(NumericQ[#] && Im[#] == 0 &), f_List] := 
    AccountingForm[Chop[N@v] , f, NumberSigns -> {"", ""}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True];
   
   (*--------------------------------------------------------------------*)
   padIt3[v_?(NumericQ[#] && Im[#] == 0 &), f_List] := 
    AccountingForm[v , f, NumberSigns -> {"", ""}, 
     NumberPadding -> {"0", "0"}, SignPadding -> True, 
     NumberPoint -> ""];
   
   (*------------------------------------------------------------------------*)
   getBobPosition[len_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &), 
     angle_?(NumericQ[#] && Im[#] == 0 &)] := Module[{bob},
     bob = {len Sin[angle], -len Cos[angle]};
     { {Red, Thick, Style[Line[{{0, 0}, bob}], Antialiasing -> True]},
      {Blue, Opacity[1], PointSize[0.06], 
       Style[Point[bob], Antialiasing -> True]}  }
     ];
   
   (*-----------------------------------------------------------------------*)
   (* keep angle from -180...+180 to make phase and poincare plots easier to make*)
   normalize[angle_?(NumericQ[#] && Im[#] == 0 &)] := 
    angle - 2 Pi Floor[(angle + Pi)/(2 Pi)];
   
   (*---------------------------------------------------------------------------------*)
   solve[q_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),  (*quality factor*)
     a_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),  (*driver amplitude*)
     omega_?(NumericQ[#] && 
         Im[#] == 0 &),           (*driver frequency*) 
     phase_?(NumericQ[#] && Im[#] == 0 &),           (*driver phase*)

     
     from_?(NumericQ[#] && 
         Im[#] == 
          0 &),            (*time to start integrating the ode*)
     to_?(NumericQ[#] && Im[#] == 0 &),              (*end time *)
     x10_?(NumericQ[#] && 
         Im[#] == 
          0 &),             (*initial condition for angle at time=0*)

     
     x20_?(NumericQ[#] && 
         Im[#] == 
          0 &),             (*initial condition for speed at time=0*)

     
     isDamped_?(Element[#, Booleans] &),           (*is Q used?*)
     accuracyGoal_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
     precisionGoal_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &)] := 
    Module[ {ic, eq1, eq2, eq3, x1, x2, x3, t, sol, ndsolveOptions},
     ndsolveOptions = {MaxSteps -> Infinity, 
       Method -> {"StiffnessSwitching", 
         Method -> {"ExplicitRungeKutta", Automatic}}, 
       AccuracyGoal -> accuracyGoal, PrecisionGoal -> precisionGoal};
     
     eq1 = x1'[t] == x2[t];
     If[isDamped,
      eq2 = x2'[t] == -1/q x2[t] - Sin[x1[t]] + a Cos[x3[t]],
      eq2 = x2'[t] == -Sin[x1[t]] + a Cos[x3[t]]
      ];
     
     eq3 = x3'[t] == omega;
     ic = {x1[0] == x10, x2[0] == x20, x3[0] == phase};
     sol = 
      First@NDSolve[ 
        Flatten@{eq1, eq2, eq3, ic}, {x1, x2, x3}, {t, from, to}, 
        Sequence@ndsolveOptions];
     {x1 /. sol, x2 /. sol, x3 /. sol}
     ];
   
   (*---------------------------------------------------------------------------------*)
   makePoincreSection[
     q_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),   (*quality factor*)
     torqueAmplitude_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),  (*driver amplitude*)
     omega_?(NumericQ[#] && 
         Im[#] == 0 &),                        (*driver frequency*) 
     phase_?(NumericQ[#] && 
         Im[#] == 0 &),                        (*driver phase*)
     x10_?(NumericQ[#] && 
         Im[#] == 0 &),    (*initial condition for angle at time=0*)
     x20_?(NumericQ[#] && 
         Im[#] == 0 &),    (*initial condition for speed at time=0*)
     isDamped_?(Element[#, Booleans] &), (*is Q used?*)
     plotImageSize_List, axesLabel_List] := 
    Module[{x1, x2, x3, i, data, accuracyGoal = 4, precisionGoal = 4, 
      pointSize = 0.01, driverPeriod, initialDriverPeriod = 150},
     
     If[Abs@omega <= $MachineEpsilon || 
       torqueAmplitude <= $MachineEpsilon,
      (
       pointSize = 0.05;
       data = {{x10, x20}}
       ),
      (
       driverPeriod = 2. Pi/Abs@omega;
       
       {x1, x2, x3} = 
        solve[q, torqueAmplitude, omega, phase, 
         initialDriverPeriod*driverPeriod, 1000*driverPeriod, x10, 
         x20, isDamped, accuracyGoal, precisionGoal];
       
       data = 
        Table[{normalize[x1[i]], x2[i]}, {i, 
          initialDriverPeriod*driverPeriod, 1000*driverPeriod, 
          driverPeriod}]
       )
      ];
     
     ListPlot[data, AspectRatio -> 1, 
      Ticks -> {{-Pi, {-Pi/2, "-\[Pi]/2"}, 0, {Pi/2, "\[Pi]/2"}, Pi}, 
        Automatic}, 
      PlotStyle -> {Blue, Directive[PointSize[pointSize]]}, 
      PerformanceGoal -> "Quality", ImageSize -> plotImageSize/2, 
      PlotLabel -> Text@Style["Poincaré section", 12, Bold], 
      PerformanceGoal -> "Quality", TicksStyle -> Directive[7], 
      AxesLabel -> {Text@Style["\[Theta]", 9], 
        Text@Style[
          Row[{" sampled ", Style["d", Italic], "\[Theta]/", 
            Style["d", Italic], "\[Tau]"}], 9]}
      ]
     ];
   
   
   makePhasePlot[data_List,(*phase data*)
     plotImageSize_List, axesLabel_List, 
     nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &)] := Module[{},
     
     ListPlot[data[[1 ;; nStep, 2 ;; 3]], Joined -> False, 
      PlotStyle -> {Red, Directive[PointSize[0.003]]}, 
      PerformanceGoal -> "Quality", 
      PlotLabel -> Text@Style["phase portrait", 12, Bold], 
      AxesLabel -> axesLabel, ImageSize -> plotImageSize/2, 
      AspectRatio -> 1, 
      PlotRange -> {{-1.1 Pi, 1.1 Pi}, {-1.2 Pi, 1.2 Pi}}, 
      Ticks -> {{-Pi, {-Pi/2, "-\[Pi]/2"}, 0, {Pi/2, "\[Pi]/2"}, Pi}, 
        Automatic}, PlotStyle -> {Blue, Directive[PointSize[0.003]]}, 
      TicksStyle -> Directive[8],
      PerformanceGoal -> "Quality"
      ]
     ];
   
   (*---------------------------------------------------------------------------------*)
   generateFinalDisplay[poincareMap_Graphics,
     phasePlot_Graphics,
     bifurcationPlot_Graphics,
     track_List,
     currentTime_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),
     plotImageSize_List, x1_, x2_, x3_,
     bottomPlot_Graphics,
     currentRunMaxSpeed_?(NumericQ[#] && Im[#] == 0 &), 
     torqueAmplitude_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &) (*driver amplitude*)] := 
    Module[{g, options, arrow, torqueArrow, speed = x2[currentTime], 
      angle = x1[currentTime], angleOffset, torqueAngleOffset, 
      gTorque},
     
     angleOffset = (Abs@speed/currentRunMaxSpeed)*0.35*Pi;
     If[angleOffset > $MachineEpsilon,
      arrow = If[speed < 0,
        arcArrow[{0, 0}, 1, angle - Pi/2 - angleOffset, angle - Pi/2 ,
          True], arcArrow[{0, 0}, 1, 
         angle - Pi/2, (angle + angleOffset) - Pi/2 , False]
        ]
      ];
     
     If[torqueAmplitude <= $MachineEpsilon,
      gTorque = {},
      (
       torqueAngleOffset = Abs@Cos[x3[currentTime]]*Pi;
       If[torqueAngleOffset > 0,
        (
         torqueArrow = If[Cos[x3[currentTime]] < 0,
           arcArrow[{0, 0}, 0.25, Pi - torqueAngleOffset , Pi, True], 
           arcArrow[{0, 0}, 0.25, -Pi, -Pi + torqueAngleOffset , False]
           ];
         gTorque = {Opacity[1], Thickness[0.010], Arrowheads[.08], 
           torqueArrow}
         )
        ];
       )
      ];
     
     g = Graphics[{
        {Dotted, RGBColor[65/255, 105/255, 225/255], Thickness[0.01], 
         Circle[{0, 0}, 1]},
        {Red, PointSize[0.03], Point[{0, 0}]},
        getBobPosition[1, angle],
        If[
         angleOffset > $MachineEpsilon, {Thick, 
          RGBColor[205/255, 92/255, 92/255], Arrowheads[.07], 
          arrow}, {}],
        gTorque
        }];
     
     options = {ImageMargins -> 1, Axes -> False, 
       PlotRange -> {{-1.2, 1.2 }, {-1.2 , 1.2}}, AspectRatio -> 1, 
       Background -> White, TicksStyle -> Directive[8], 
       ImagePadding -> 4, ImageSize -> plotImageSize/2, 
       AxesStyle -> Directive[{Blue, Thickness[0.01]}]};
     
     Grid[{
       {poincareMap, Show[phasePlot, Graphics[track]]},
       {bifurcationPlot, Show[g, Sequence@options]},
       {bottomPlot, SpanFromLeft}
       }, Frame -> All, 
      FrameStyle -> Directive[Thickness[.005], Gray], 
      Spacings -> {0.5, 0.4}
      ]
     
     ];
   (*---------------------------------------------------------------------------------*)
   makeBifurcationPlot[
     aStart_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &), (*driver amplitude starting value*)
     aLen_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] > 0 &),     (*driver amplitude range of values length*)
     aIntervals_?(IntegerQ[#] && Im[#] == 0 && 
         Re[#] > 0 &),(*how many intervals to do *)
     q_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),      (*quality factor*)
     omega_?(NumericQ[#] && 
         Im[#] == 0 &),              (*driver frequency*) 
     phase_?(NumericQ[#] && 
         Im[#] == 0 &),              (*driver phase*)
     x10_?(NumericQ[#] && 
         Im[#] == 
          0 &),                (*initial condition for angle at time=
     0*)x20_?(NumericQ[#] && 
         Im[#] == 
          0 &),                (*initial condition for speed at time=
     0*)
     isDamped_?(Element[#, Booleans] &),
     plotImageSize_List,
     bifurcationProgress_ (*by reference*)] := 
    Module[{driverPeriod, timeValues, valuesOfDriverAmplitude, x1, x2,
       x3, i, j, p, data, initialDriverPeriod = 150, 
      finalDriverPeriod = 250, plotOptions, total, count, 
      currentTorqueAmplitude, pointSize = 0.004},
     
     plotOptions = {Joined -> False, 
       PlotStyle -> {Red, Directive[PointSize[pointSize]]}, 
       PerformanceGoal -> "Quality", 
       PlotLabel -> Text@Style["bifurcation map", 12, Bold], 
       ImageSize -> plotImageSize/2, AspectRatio -> 0.85, 
       PlotRange -> {{aStart, aStart + aLen}, All}, 
       PlotStyle -> {Blue, Directive[PointSize[0.003]]}, 
       TicksStyle -> Directive[8], 
       AxesLabel -> {Text@Style["a", Italic, 9], 
         Text@Style[
           Row[{" sampled ", Style["d", Italic], "\[Theta]/", 
             Style["d", Italic], "\[Tau]"}], 9]}, 
       PerformanceGoal -> "Quality"};
     
     If[Abs@omega <= $MachineEpsilon,
      (
       data = {{0, 0}};
       p = ListPlot[data, Sequence@plotOptions]
       )
      ,
      (
       driverPeriod = 2. Pi/Abs@omega;
       valuesOfDriverAmplitude = 
        Range[aStart, aStart + aLen, aLen/aIntervals]; 
       timeValues = 
        Range[initialDriverPeriod*driverPeriod, 
         finalDriverPeriod*driverPeriod, driverPeriod ];
       
       data = 
        Table[0, {Length[valuesOfDriverAmplitude]}, {Length[
           timeValues ]}];
       total = Length[timeValues]*Length[valuesOfDriverAmplitude];
       bifurcationProgress = 0.0;
       count = 0.0;
       Do[
        (
         {x1, x2, x3} = 
          solve[q, valuesOfDriverAmplitude [[i]], omega, phase, 
           initialDriverPeriod*driverPeriod, (finalDriverPeriod)*
            driverPeriod, x10, x20, isDamped, 4, 4];
         
         currentTorqueAmplitude = valuesOfDriverAmplitude[[i]];
         Do[
          (
           
           data[[i, j]] = {currentTorqueAmplitude, 
             x2[timeValues[[j]]  ]};
           count += 1;
           bifurcationProgress = count/total
           ), {j, 1, Length[timeValues]}
          ]
         ), {i, 1, Length[valuesOfDriverAmplitude]}
        ];
       p = ListPlot[Flatten[data, 1], Sequence@plotOptions]
       )
      ];
     
     bifurcationProgress = 0;
     p
     ];
   
   (*---------------------------------------------------------------------------------*)
   makeTimeSeriesPlot[data_List, plotImageSize_List, 
     nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &), 
     currentTime_] := Module[{},
     
     ListPlot[
      Transpose[{data[[1 ;; nStep, 1]], data[[1 ;; nStep, 3]]}], 
      Joined -> True, GridLines -> Automatic, 
      PlotStyle -> {Red, Directive[PointSize[0.001]]}, 
      PerformanceGoal -> "Quality", 
      ImageSize -> {370, .44*plotImageSize}, AspectRatio -> .30, 
      PlotStyle -> {Blue, Directive[PointSize[0.003]]}, 
      AxesOrigin -> {0, 0}, ImagePadding -> {{30, 2}, {15, 35}}, 
      PlotRange -> {{0, Automatic}, All}, 
      PerformanceGoal -> "Quality", Frame -> True, 
      FrameTicksStyle -> Directive[7], Axes -> None, 
      FrameLabel -> {{None, None}, {None, 
         Text@Row[{Style["time series", 11, Bold], Spacer[5], 
            Style[Row[{Style["d", Italic], "\[Theta]/", 
               Style["d", Italic], 
               "\[Tau] vs. \[Tau] (dimensionless)"}], 9], Spacer[5], 
            Style["time: ", 10], 
            Text@Style[padIt2[currentTime, {5, 2}], 11]}]}}
      ]
     
     ];
   (*---------------------------------------------------------------------------------*)
   makePowerSpectrumPlot[data_List, plotImageSize_List,
     nStep_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
     omega_?(NumericQ[#] && Im[#] == 0 &),
     spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &), 
     currentTime_] := Module[{list, len},
     
     list = 
      Fourier[data[[1 ;; nStep, 3]], FourierParameters -> {1, -1}];
     list = If[Abs[omega] > 2*$MachineEpsilon,
       (Abs@list)^2/(nStep*Abs@omega),
       (Abs@list)^2/nStep
       ];
     
     len = Length[list];
     
     If[len > 1, list = list[[1 ;; Floor[Length[list]/2] ]]];
     
     ListPlot[If[spectrumScale == 1, list, Log[10, list]], 
      Joined -> True, GridLines -> Automatic, 
      PlotStyle -> {Red, Directive[PointSize[0.001]]}, 
      PerformanceGoal -> "Quality", 
      ImageSize -> {370, 0.44*plotImageSize}, AspectRatio -> .30, 
      PlotStyle -> {Blue, Directive[PointSize[0.003]]}, 
      AxesOrigin -> {0, 0}, PlotRange -> {{0, Automatic}, All}, 
      ImagePadding -> {{30, 2}, {15, 35}}, 
      PerformanceGoal -> "Quality", Frame -> True, 
      FrameTicksStyle -> Directive[8], Axes -> None, 
      RotateLabel -> False, 
      FrameLabel -> {{None, None}, {None, 
         Text@Row[{Style["power spectrum of angular velocity", 11, 
             Bold], Spacer[5], Style[" time: ", 10], 
            Text@Style[padIt2[currentTime, {5, 2}], 11]}]}   }
      ]
     ];
   (*---------------------------------------------------------------------------------*)
   setParameters[testCase_String] := 
    Module[{a, omega, phase, aStart, aLen, aIntervals, x10, isDamped, 
      x20, q},
     Which[
      testCase == "periodic, A=0.9,Q=2,omega=2/3", (a = 0.9; q = 2; 
       omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6; 
       aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
           
      testCase == "periodic, A=1.35,Q=2,omega=2/3", (a = 1.35; q = 2; 
       omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.7; 
       aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
      
      testCase == "period doubling, A=1.12,Q=2,omega=2/3", (a = 1.12; 
       q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.7; 
       aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
      
      testCase == "period doubling, A=1.07,Q=2,omega=2/3", (a = 1.07; 
       q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6; 
       aIntervals = 100; x10 = 45*Pi/180; isDamped = True; x20 = 0),
      
      testCase == "period doubling, A=1.45,Q=2,omega=2/3", (a = 1.45; 
       q = 2; omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6; 
       aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
      
      testCase == 
       "period quadrupling, A=1.47,Q=2,omega=2/3", (a = 1.47; q = 2; 
       omega = 2/3; phase = 0; aStart = 0.9; aLen = 0.6; 
       aIntervals = 100; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
      
      testCase == "chaotic, A=1.15,Q=2,omega=2/3",
      (a = 1.15; q = 2; omega = 2/3; phase = 0; aStart = 0.9; 
       aLen = 0.7; aIntervals = 50; x10 = 45*Pi/180.; isDamped = True;
        x20 = 0),
      
      testCase == "chaotic, A=1.5,Q=2,omega=2/3",
      (a = 1.5; q = 2; omega = 2/3; phase = 0; aStart = 0.9; 
       aLen = 0.7; aIntervals = 50; x10 = 45*Pi/180.; isDamped = True;
        x20 = 0),
      
      testCase == "chaotic, A=1.5,Q=2,omega=2/3,phase=3.14", (a = 1.5;
        q = 2; omega = 2/3; phase = 3.14; aStart = 0.9; aLen = 0.7; 
       aIntervals = 50; x10 = 45*Pi/180.; isDamped = True; x20 = 0),
      
      testCase === "chaotic, Bifurcation #1",
      (a = 1.5; q = 2; omega = 2/3; phase = 0; aStart = 0.9; 
       aLen = 0.2; aIntervals = 70; x10 = 45*Pi/180.; isDamped = True;
        x20 = 0),
      
      testCase === "chaotic, Bifurcation #2",
      (a = 1.15; q = 4; omega = 2/3; phase = 0; aStart = 0.9; 
       aLen = 0.8; aIntervals = 80; x10 = 45*Pi/180.; isDamped = True;
        x20 = 0),
      
      testCase === "chaotic, Bifurcation #3",
      (a = 1.5; q = 2; omega = -1.33; phase = 1.8; aStart = 0.9; 
       aLen = 0.6; aIntervals = 150; x10 = 50*Pi/180.; 
       isDamped = True; x20 = 1.642),
           
      testCase === "double periodic #1", (a = 1.08; q = 2; 
       omega = N[2/3]; phase = 0; aStart = 1.050; aLen = 0.08; 
       aIntervals = 300; x10 = 0; isDamped = True; x20 = 0)
      ];
     
     {a, q, omega, phase, aStart, aLen, aIntervals, x10, isDamped, x20}
     ];
   
   (*---------------------------------------------------------------------------------*)
   initializeSystem[
     q_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),(*quality factor*)
     torqueAmplitude_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &), (*driver amplitude*)
     omega_?(NumericQ[#] && 
         Im[#] == 0 &),                       (*driver frequency*) 
     phase_?(NumericQ[#] && 
         Im[#] == 0 &),                       (*driver phase*)
     x10_?(NumericQ[#] && 
         Im[#] == 0 &),   (*initial condition for angle at time=0*)
     x20_?(NumericQ[#] && 
         Im[#] == 0 &),   (*initial condition for speed at time=0*)
     isDamped_?(Element[#, Booleans] &),
     plotImageSize_List,
     axesLabel_List,
     plotChoice_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
     pphaseData_,
     spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
     currentTime_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),
     dt_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
     maxRunTime_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &)] := 
    Module[{phaseData = pphaseData, nSteps, poincareMap, phasePlot, 
      bottomPlot, track, x1, x2, x3, i, range, maxSpeed},
     
     (*solve first to find maximum speed over the whole run, 
     needed for plotting relative arrow*)
     {x1, x2, x3} = 
      solve[q, torqueAmplitude, omega, phase, 0, maxRunTime, x10, x20,
        isDamped, 4, 4];
     maxSpeed = Max@Abs@Table[x2[i], {i, 0, maxRunTime, dt}];
     
     {x1, x2, x3} = 
      solve[q, torqueAmplitude, omega, phase, 0, currentTime + 1, x10,
        x20, isDamped, 4, 4];
     nSteps = Length[Range[0, currentTime, dt]];
     phaseData[[1 ;; nSteps]] = 
      Table[{i, normalize[x1[i]], x2[i]}, {i, 0, currentTime, dt}];
     poincareMap = 
      makePoincreSection[q, torqueAmplitude, omega, phase, x10, x20, 
       isDamped, plotImageSize, axesLabel]; 
     phasePlot = 
      makePhasePlot[phaseData, plotImageSize, axesLabel, nSteps];
     
     bottomPlot = If[plotChoice == 1,
       makeTimeSeriesPlot[phaseData, plotImageSize, nSteps, 
        currentTime], 
       makePowerSpectrumPlot[phaseData, plotImageSize, nSteps, omega, 
        spectrumScale, currentTime]
       ];
     
     track = {Blue, PointSize[0.03], 
       Point[ { normalize@x1[currentTime], x2[currentTime]}]};
     {nSteps, poincareMap, phasePlot, bottomPlot, track, phaseData, 
      x1, x2, x3, maxSpeed} 
     ];
   
   (*---------------------------------------------------------------------------------*)
   makeOneStep[
     q_?(NumericQ[#] && Im[#] == 0 && Re[#] >= 0 &),(*quality factor*)
     torqueAmplitude_?(NumericQ[#] && Im[#] == 0 && 
         Re[#] >= 0 &),(*driver amplitude*)
     omega_?(NumericQ[#] && 
         Im[#] == 0 &),                      (*driver frequency*) 
     phase_?(NumericQ[#] && 
         Im[#] == 0 &),                       (*driver phase*)
     x10_?(NumericQ[#] && 
         Im[#] == 0 &), (*initial condition for angle at time=0*)
     x20_?(NumericQ[#] && 
         Im[#] == 0 &), (*initial condition for speed at time=0*)
     isDamped_?(Element[#, Booleans] &),
     plotImageSize_List,
     axesLabel_List,
     plotChoice_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
     pphaseData_,
     spectrumScale_?(IntegerQ[#] && Im[#] == 0 && Re[#] > 0 &),
     ccurrentTime_,
     nnSteps_,
     dt_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
     pphasePlot_,
     bbottomPlot_,
     ttrack_,
     maxRunTime_?(NumericQ[#] && Im[#] == 0 && Re[#] > 0 &),
     xx1_,
     xx2_,
     xx3_ ] := 
    Module[{x1 = xx1, x2 = xx2, x3 = xx3, track = ttrack, 
      bottomPlot = bbottomPlot, phasePlot = pphasePlot, 
      nSteps = nnSteps, currentTime = ccurrentTime, 
      phaseData = pphaseData},
     
     If[currentTime + dt <= maxRunTime,
      (
       nSteps += 1;
       currentTime += dt;
       {x1, x2, x3} = 
        solve[q, torqueAmplitude, omega, phase, currentTime, 
         currentTime + 1, x10, x20, isDamped, 7, 7]; 
       phaseData[[nSteps]] = {currentTime, normalize@x1[currentTime], 
         x2[currentTime]}; 
       phasePlot = 
        makePhasePlot[phaseData, plotImageSize, axesLabel, nSteps];
       
       If[plotChoice == 1,
        bottomPlot = 
         makeTimeSeriesPlot[phaseData, plotImageSize, nSteps, 
          currentTime], 
        bottomPlot = 
         makePowerSpectrumPlot[phaseData, plotImageSize, nSteps, 
          omega, spectrumScale, currentTime]
        ];
       track = {Blue, PointSize[0.03], 
         Point[ { normalize@x1[currentTime], x2[currentTime] }  ] };
       )
      ];
     
     {currentTime, nSteps, phasePlot, bottomPlot, track, phaseData, 
      x1, x2, x3} 
     ];
   
   (*---------------------------------------------------------------------------------*)
   arcArrow[centre_, radius_, phi1_, phi2_, flip_, resolution_: 10] :=
     Module[{i, data}, 
     data = Table[
       centre + radius {Cos[i], Sin[i]}, {i, phi1, 
        phi2, (phi2 - phi1)/Ceiling[(phi2 - phi1)/(Pi/resolution)]}
       ];
     If[flip, data = Reverse[data]];
     Arrow@data
     ]
   
   )
 ]