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[Δτ]],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["θ''(τ) +  θ'(τ)/", "TR",12], Style["Q",Italic,"TR",12], Style["+ sin θ(τ)) == ","TR",12] ,Style["a",Italic,"TR",12],Style[" cos ϕ(τ)", "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["ω = dϕ/dτ",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["θ (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θ/dτ (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["ϕ (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["θ",9],Text@Style[Row[{" ",Style["d",Italic],"θ/",Style["d",Italic],"τ"}],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,"-π/2"},0,{Pi/2,"π/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["θ",9],Text@Style[Row[{" sampled ",Style["d",Italic],"θ/",Style["d",Italic],"τ"}],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.1Pi,1.1Pi},{-1.2 Pi,1.2 Pi}},Ticks->{{-Pi,{-Pi/2,"-π/2"},0,{Pi/2,"π/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],"θ/",Style["d",Italic],"τ"}],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],"θ/",Style["d",Italic],"τ vs. τ (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
]

)
]

Theory background: θ''+c θ'+g/L sin θ=A sin((Subscript[ω, f]) t) is a second-order nonlinear differential equation that models the motion of a simple pendulum subject to a forcing function A sin(ω t), where c is the damping coefficient, L is the length of the pendulum, and θ is its angular position.
For ease of analysis, this differential equation is first converted to a nondimensional form by defining t=τ/(Subscript[ω, n]), (Subscript[ω, f])=ω (Subscript[ω, n]), Q=(Subscript[ω, n])/c, and A=a (Subscript[ω, n]^2), where (Subscript[ω, n]) is the pendulum’s undamped natural frequency given by (Subscript[ω, n])=(Sqrt[g/L]) and g is the gravitational constant. 
Applying this transformation results in dτ/dt=(Subscript[ω, n]), and using dθ/dt=(dθ/dτ)(dτ/dt)  results in θ''=(Subscript[ω, n]^2)((d^2)θ/d(τ^2)) , θ'=(Subscript[ω, n])dθ/dτ, and A sin((Subscript[ω, f]) t)=a (Subscript[ω, n]^2)cos(ω τ). Substituting these back into the original differential equation and dividing throughout by (Subscript[ω, n]^2) gives θ''+θ'/Q+sin θ=a sin(ω τ), where the derivatives are now with respect to τ. Q is called the quality factor. This is the nonlinear differential equation that is used in this Demonstration, not the original one.
The driver frequency ω is given by ω=dϕ/dτ, where ϕ is the phase angle. This system is now converted to matrix form X'=A X, which contains three first-order differential equations as follows. Let (Subscript[x, 1])=θ, (Subscript[x, 2])=dθ/dτ, and (Subscript[x, 3])=ϕ. Taking derivatives gives (Subscript[x, 1])'=(Subscript[x, 2]), (Subscript[x, 2])'=-(Subscript[x, 2])/Q-sin((Subscript[x, 1]))+a sin((Subscript[x, 3])), and (Subscript[x, 3])'=ω. The system in state can thus be written (({{(Subscript[x, 1])'}, {(Subscript[x, 2])'}, {(Subscript[x, 3])'}}))==(({{(Subscript[x, 2])}, {a cos((Subscript[x, 3]))-sin((Subscript[x, 1]))-((Subscript[x, 2])/Q)}, {ω}})).
The Mathematica built-in function NDSolve[] is now used for this system. The initial conditions are (Subscript[x, 1])(0)=θ(0), (Subscript[x, 2])(0)=θ'(0), and (Subscript[x, 3])(0)=ϕ(0). Therefore, the three dynamic variables (or state-space variables) in this model are θ(τ), dθ/dτ, and ϕ(τ).
This describes how to use the Demonstration: The top buttons are used to run and control the simulation. You can vary the simulation time up to a maximum of 500 units. Use the “pause” button to pause the simulation at any time. Use the “step” button to advance the time by one unit. Use the “reset” button to bring the system back to initial time. Use the button labeled “fast” to advance the system to the end of the simulation period immediately without going through the animation steps. Use the slider labeled “Δτ” to adjust the size of the time step. The smaller the step size, the more accurate and smoother the motion will appear, but it will take a longer time to run. 
The parameters of the system Q, a, and ω are varied using the sliders. Making Q zero will remove damping from the system, while making a zero will remove the forcing function. This effectively makes the system run as a standard simple pendulum with no damping and no forcing function.
There are five main plots on the display. The Poincaré map is always plotted to reflect the state at the end of the simulation period. The y axis represents the values of the angular velocity of the pendulum sampled at time intervals equal in magnitude to the period of the forcing function. The phase plot is updated as time is advanced and at each step, and similarly for the time series plot and the power spectrum plot. 
Switch between the time series plot and the power spectrum plot using the option selection shown in the left part of the display. The phase portrait was drawn such that the angle of the pendulum is kept between -π and π to make it easier to draw the plot.
The display of the actual physical pendulum is also shown with a small arrow attached to the bob. The length of the arrow indicates the relative magnitude of the current angular velocity. The inner arrow represents the relative strength of the external torque. The directions of the arrows represent the sign of the respective values. An anticlockwise arrow direction indicates a positive sign.
The bifurcation plot is only updated when you click the button labeled “generate”. This lets you choose your parameters before viewing the final result, as it takes time to generate the plot. The parameters used to select the bifurcation plot let you select the range of values of the parameter a (the torque amplitude). The more intervals you ask for, the clearer the plot will be, but it takes longer to generate as the system is numerically solved for each interval. Sampling of values starts after the time has reached 150 periods of the driver frequency, in order to allow transient effects to attenuate sufficiently. The solution is then run for an additional 100 periods to collect the data for the plot. This is done at each interval. Therefore, the y axis for the bifurcation is similar to the y axis of the Poincaré map; however, the x axis represents a and not θ as in the phase portrait and Poincaré map.
The current time for the simulation is displayed on the right side of the title of the bottom plot. When the system reaches the end of the simulation, the “reset” button can be used bring it back to the initial time. The simulation will not advance beyond 500 time units. 
The Demonstration comes with a number of test cases preconfigured to show different aspect of the behavior of the system. These test cases are selected from the popup menu at the top right. Once a test case is selected, the user interface parameters are updated automatically, and the final simulation result is shown based on the use of these new parameters.
Possible issues:
1. The maximum simulation duration of 500 was fixed in the program for efficient implementation. Hence this value should not be increased without adjusting the corresponding internal variable in the software that depends on this value.
2. When the current time is at the maximum duration time, clicking the “run” or “step” button has no effect. The “reset” button can be used to bring the system to the initial state in order to run a new simulation.
References
[1] R. Fitzpatrick. “Introduction to Computational Physics.” (2006) http://farside.ph.utexas.edu/teaching/329/lectures/node47.html.
[2] R. Gilmore and C. Letellier, The Symmetry of Chaos, Oxford: Oxford University Press, 2007.
[3] C. Greene. “Driven Pendulum.” (2000) http://www.phy.davidson.edu/StuHome/chgreene/chaos/Pendulum/pendulum_content_frame.htm.
[4] G. Leo. “Chaotic Dynamics of a Damped, Driven Pendulum.” (2007) http://homepage.ntlworld.com/gerry.leo/Theory.html.