(*by Nasser M. Abbasi. 12000.org  8/23/2023 *)
Manipulate[
Row[{Dynamic[Refresh[Which[state=="RESET",
(
currentTime = 0;
tick = 0;

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase];

state = "PAUSE"
),

state=="PAUSE",
(
Which[lastEvent=="duration_changed"||lastEvent=="delt_changed",
(
lastEvent="no_event";

If[showPhase,
phasePortraitPlot=makePhasePortrait[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,c,maxRunTime,dt]
]
),

lastEvent=="show_phase",
(
lastEvent="no_event";
phasePortraitPlot=makePhasePortrait[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,c,maxRunTime,dt]
),

lastEvent=="run_button",(lastEvent="no_event";state="RUNNING"; tick+=DEL),

lastEvent=="initial_conditions_changed",
(
lastEvent="no_event";
currentTime=0;

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase]
),

lastEvent=="step_button",
(
lastEvent="no_event";

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase];

If[currentTime+dt>maxRunTime,
(
state="RESET";
tick += DEL
),
(
currentTime+=dt
)
]
),
lastEvent=="reset_button", (lastEvent="no_event";state="RESET";tick += DEL),
lastEvent=="pause_button",lastEvent="no_event",
lastEvent=="mouseDown",lastEvent="no_event",
lastEvent=="mouseUp",
(
lastEvent="no_event";
currentTime=0
)
]
),

state=="RUNNING",
(
Which[lastEvent=="show_phase",
(
lastEvent="no_event";
phasePortraitPlot=makePhasePortrait[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,c,maxRunTime,dt];
tick += DEL
),

lastEvent=="duration_changed"||lastEvent=="delt_changed",
(
lastEvent="no_event";
If[showPhase,
phasePortraitPlot=makePhasePortrait[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,c,maxRunTime,dt]
];
tick += DEL
),

lastEvent=="no_event"||lastEvent=="run_button",
(
If[lastEvent=="run_button",lastEvent="no_event"];

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase];

If[currentTime+dt>maxRunTime,
(
currentTime=0
),
(
currentTime+=dt
)
];

tick += DEL
),

lastEvent=="pause_button",
(
lastEvent="no_event";
state = "PAUSE"
),

lastEvent=="mouseDown",(lastEvent="no_event"),

lastEvent=="mouseUp",
(
lastEvent="no_event";
currentTime=0;
tick += DEL
),

lastEvent=="reset_button",(lastEvent="no_event";state="RESET";tick += DEL),

lastEvent=="step_button", 
(
lastEvent="no_event";

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase];

If[currentTime+dt>maxRunTime,currentTime=0,currentTime+=dt];
state="PAUSE"
),

lastEvent=="initial_conditions_changed",
(
lastEvent="no_event";
currentTime=0;
{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}=update[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,c,maxRunTime,dt,showPhase];

tick += DEL
),

lastEvent=="pause_button",(lastEvent="no_event"; state="PAUSE")
]
)
];"",TrackedSymbols->{tick}]],

EventHandler[
Dynamic[If[showPhase,
Refresh[Grid[{{Framed[Show[phasePortraitPlot,Graphics[{Blue,PointSize[0.02],Point[{\[Theta]1[currentTime],\[Theta]1Speed[currentTime]}]}]],FrameStyle->Directive[Thickness[.005],Gray]]},{Framed[Graphics[
{
{RGBColor[{188,143,143}/255],Circle[{0,0},L1]},
{Circle[{0,0},0.05]},
getCoordinates[nBobs,bob1,bob2,bob3,m1,m2,m3]
},
graphicsOptions2D[6.2,{370,312}]
],FrameStyle->Directive[Thickness[.005],Gray]
]
}},Alignment->Center],TrackedSymbols->{tick,bob1,bob2,bob3}]
,
Refresh[Framed[Graphics[
{
{RGBColor[{188,143,143}/255],Circle[{0,0},L1]},
{Circle[{0,0},0.05]},

getCoordinates[nBobs,bob1,bob2,bob3,m1,m2,m3]
},
graphicsOptions2D[6.2,{370,497}]
],FrameStyle->Directive[Thickness[.005],Gray]],TrackedSymbols->{tick,bob1,bob2,bob3}]
]
]
,
{
"MouseDown":>
(lastEvent="mouseDown";
{\[Theta]1InitMouse,\[Theta]2InitMouse,\[Theta]3InitMouse}={\[Theta]1Init,\[Theta]2Init,\[Theta]3Init}
),

"MouseDragged":>
(
p1=MousePosition["Graphics"];

{bob1,bob2,bob3,\[Theta]1InitMouse,\[Theta]2InitMouse,\[Theta]3InitMouse}=obtainScreenPositions[p1,bob1,bob2,bob3,\[Theta]1InitMouse,\[Theta]2InitMouse,\[Theta]3InitMouse,L1,L2,L3,nBobs];
),
"MouseUp":>
(
currentTime=0;
lastEvent="mouseup";
p1=MousePosition["Graphics"];
{\[Theta]1Init,\[Theta]2Init,\[Theta]3Init}={\[Theta]1InitMouse,\[Theta]2InitMouse,\[Theta]3InitMouse};
tick +=DEL
)
}
]
}],
(*----------------- controls --------------------------------*)
Item[
Grid[{

{Grid[{(* BLOCK 1 *)
{
Button[Text[Style["play",14]],(lastEvent="run_button";tick+=DEL),ImageSize->{74,30}],
Button[Text[Style["pause",14]],(lastEvent="pause_button";tick+=DEL),ImageSize->{74,30}]
},
{
Button[Text[Style["step",14]],(lastEvent="step_button"; tick+=DEL),ImageSize->{74,30}], Button[Text[Style["reset",14]],(lastEvent="reset_button";tick+=DEL),ImageSize->{74,30}]
}
},Spacings->{.8,.3},Alignment->Center
]
},


{Grid[{
{Text@Style["duration",12],

Manipulator[Dynamic[maxRunTime,(maxRunTime=#;{lastEvent="duration_changed",tick+=DEL};#)&],{1,100,1},ImageSize->Tiny,ContinuousAction->False],
Dynamic@Text@Style[padIt2[maxRunTime,{5,0}],10]},

{Style[Row[{"\[CapitalDelta]",Style["t", Italic]}],"TR",12],
Manipulator[Dynamic[dt,(dt=#;{lastEvent="delt_changed",tick+=DEL};#)&],{0.01,0.1,0.01},ImageSize->Tiny,ContinuousAction->False],Dynamic@Text@Style[padIt2[dt,{3,3}],10]
}
}]
},

{Grid[{
{Text@Style["number of bobs",12],
RadioButtonBar[Dynamic[nBobs,(nBobs=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1->Text@Style["1",10],2->Text@Style["2",10],3->Text@Style["3",10]}]
}}
,Frame->None,Spacings->{0.2,0.4},FrameStyle->Directive[Thickness[.005],Gray]
]
},

{Grid[{(*adjust mass*)

{Button[ Text@Style["min",10],(m1=1;lastEvent="initial_conditions_changed";tick+=DEL),
ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["m",Italic],"1"],"TR",10],Manipulator[Dynamic[m1,(m1=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,20,0.01},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[m1,{5,3}],11]},

{Button[ Text@Style["min",10],(m2=1;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["m",Italic],"2"],"TR",10],Manipulator[Dynamic[m2,(m2=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,20,0.01},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[m2,{5,3}],11]},

{Button[ Text@Style["min",10],(m3=1;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["m",Italic],"3"],"TR",10],Manipulator[Dynamic[m3,(m3=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,20,0.01},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[m3,{5,3}],11]},

{Button[ Text@Style["min",10],(L1=1;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["L",Italic],"1"],"TR",10],Manipulator[Dynamic[L1,(L1=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,2,0.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[L1,{3,2}],11]},

{Button[ Text@Style["min",10],(L2=1;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["L",Italic],"2"],"TR",10],Manipulator[Dynamic[L2,(L2=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,2,0.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[L2,{3,2}],11]},

{Button[ Text@Style["min",10],(L3=1;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],Text@Style[Subscript[Style["L",Italic],"3"],"TR",10],Manipulator[Dynamic[L3,(L3=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1,2,0.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt2[L3,{3,2}],11]}

},Frame->None, Alignment->Left]
},

{Grid[{(*initial conditions*)

{Button[ Text@Style["zero",10],(\[Theta]1Init=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Subscript["\[Theta]","1"],"TR",10],Manipulator[Dynamic[\[Theta]1Init,(\[Theta]1Init=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-N@Pi,N@Pi,Pi/100.},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]1Init,{4,1}],10]},

{Button[ Text@Style["zero",10],(\[Theta]2Init=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Subscript["\[Theta]","2"],"TR",10],Manipulator[Dynamic[\[Theta]2Init,(\[Theta]2Init=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-N@Pi,N@Pi,Pi/100.},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]2Init,{4,1}],10]},

{Button[ Text@Style["zero",10],(\[Theta]3Init=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Subscript["\[Theta]","3"],"TR",10],Manipulator[Dynamic[\[Theta]3Init,(\[Theta]3Init=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-N@Pi,N@Pi,Pi/100.},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]3Init,{4,1}],10]},

{Button[ Text@Style["zero",10],(\[Theta]1SpeedInit=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Overscript[Subscript["\[Theta]","1"],"\[Bullet]"],"TR",10],Manipulator[Dynamic[\[Theta]1SpeedInit,(\[Theta]1SpeedInit=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-3,3,.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]1SpeedInit,{3,2}],10]},

{Button[ Text@Style["zero",10],(\[Theta]2SpeedInit=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Overscript[Subscript["\[Theta]","2"],"\[Bullet]"],"TR",10],Manipulator[Dynamic[\[Theta]2SpeedInit,(\[Theta]2SpeedInit=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-3,3,.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]2SpeedInit,{3,2}],10]},

{Button[ Text@Style["zero",10],(\[Theta]3SpeedInit=0;lastEvent="initial_conditions_changed";tick+=DEL),ImageSize->{40,20}, Alignment-> Bottom],Spacer[2],
Text@Style[Overscript[Subscript["\[Theta]","3"],"\[Bullet]"],"TR",10],Manipulator[Dynamic[\[Theta]3SpeedInit,(\[Theta]3SpeedInit=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{-3,3,.1},ImageSize->Tiny,ContinuousAction->False],Spacer[2],Dynamic@Text@Style[padIt1[\[Theta]3SpeedInit,{3,2}],10]}

},Frame->None,Alignment->Left]},

{Dynamic[makePeKeChart[currentPE,currentKE]]}
},Frame->All,FrameStyle->Directive[Thickness[.005],Gray],Spacings->{0.3,0.5},Alignment->Center]
,ControlPlacement->Left],

(*RIGHT PANEL*)
Item[
Grid[{

{Grid[{
{Text@Style["gravity",10]},

{PopupMenu[Dynamic[g,(g=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{1.63->Text@Style["moon",12],9.81->Text@Style["earth",11],274.68->Text@Style["sun",11]},ImageSize->Tiny]
}
},Frame->None,Alignment->Center,Spacings->{0,.6}]},

{Grid[{
{Text@Style["damping",10]},

{VerticalSlider[Dynamic[c,(c=#;{lastEvent="initial_conditions_changed",tick+=DEL};#)&],{0,3,.1},ImageSize->Tiny,ContinuousAction->False]},

{Dynamic@Text@Style[padIt2[c,{3,1}],10]}
},Frame->None,Alignment->Center,Spacings->{0,.6}]},

{Grid[{
{Text@Style["show phase",10]},
{Checkbox[Dynamic[showPhase,(showPhase=#;{lastEvent="show_phase",tick+=DEL};#)&]]}
},Frame->None,Alignment->Center,Spacings->{0,.6}]},

{Grid[{
{Text@Style["current time",10],Spacer[5]},
{Dynamic@Text[Style[padIt2[currentTime,{5,2}],10]]}
},Frame->None,Alignment->Center,Spacings->{0,.6}]}

},Frame->All,FrameStyle->Directive[Thickness[.005],Gray],Spacings->{0.3,0.5},Alignment->Center
],ControlPlacement->Right
],

{{currentPE,1},ControlType->None},
{{currentKE,0},ControlType->None},
{{showPhase,False},None},
{{phasePortraitPlot,{}},None},
{{g,9.81},None},
{{c,0},None},
{{maxRunTime,100},None},
{{dt,0.05},None},
{{bob1,{0,0}},ControlType->None},
{{bob2,{0,0}},ControlType->None},
{{bob3,{0,0}},ControlType->None},
{{\[Theta]1,0},None},
{{\[Theta]1Speed,0},None},
{{nBobs,3},None},
{{m1,20},None},
{{m2,15},None},
{{m3,20},None},
{{L1,2},None},
{{L2,2},None},
{{L3,2},None},
{{\[Theta]1InitMouse,0},None},
{{\[Theta]2InitMouse,0},None},
{{\[Theta]3InitMouse,0},None},
{{state,"RESET"},None}, (*adjust this to RESET if do not want it to start in run mode*)
{{\[Theta]1Init,0.6},None},
{{\[Theta]2Init,0.1},None},
{{\[Theta]3Init,0.1},None},
{{\[Theta]1SpeedInit,2.3},None},
{{\[Theta]2SpeedInit,0},None},
{{\[Theta]3SpeedInit,0},None},
{{p1,0},None},
{{lastEvent,"no_event"},None},
{{currentTime,0},None},
{{tick,0},None},
{{DEL,$MachineEpsilon},None},

TrackedSymbols->{None},
ContinuousAction->False,
SynchronousUpdating->False,
SynchronousInitialization->False,
ControlPlacement->Left,
Alignment->Center,
ImageMargins->0,
FrameMargins->0,

Initialization:>
(
Off[InterpolatingFunction::dmval];

(*formating functions*)
(*---------------------------------------------------------*)
padIt1[v_?(NumberQ[#]&),f_List]:=AccountingForm[Chop[N@v] ,f,NumberSigns->{"-","+"},NumberPadding->{"0","0"},SignPadding->True
];

(*---------------------------------------------------------*)
padIt2[v_?(NumberQ[#]&),f_List]:=AccountingForm[Chop[N@v] ,f,NumberSigns->{"",""},NumberPadding->{"0","0"},SignPadding->True
];

(*---------------------------------------------------------*)
graphicsOptions2D[maxExt_,imageSize_]:=Module[{frameThickness=0.001},
{
ImageSize->imageSize,ImagePadding->2,ImageMargins->0,
Axes->False,
PlotRange->{{-maxExt,maxExt },{-maxExt ,maxExt }},
AspectRatio->1,
Background->White,
TicksStyle->Small,
PlotRangePadding->None,
AxesStyle->Directive[{Blue,Thickness[frameThickness]}],
AspectRatio->1
}
];

(*---------------------------------------------------------*)
getCoordinates[nBobs_,bob1_,bob2_,bob3_,m1_,m2_,m3_]:=Module[{slope=0.03/19,c=0.01-(0.03/19)},

Which[nBobs==1,
{
{Red,Style[Line[{{0,0},bob1}],Antialiasing->True]},
{Blue,PointSize[slope*m1+c],Style[Point[bob1],Antialiasing->True]}
},
nBobs==2,
{
{Red,Style[Line[{{0,0},bob1,bob2}],Antialiasing->True]},
{Blue,PointSize[slope*m1+c],Style[Point[bob1],Antialiasing->True]},
{Blue,PointSize[slope*m2+c],Style[Point[bob2],Antialiasing->True]}
},
True,
{
{Red,Style[Line[{{0,0},bob1,bob2,bob3}],Antialiasing->True]},
{Blue,PointSize[slope*m1+c],Style[Point[bob1],Antialiasing->True]},
{Blue,PointSize[slope*m2+c],Style[Point[bob2],Antialiasing->True]},
{Blue,PointSize[slope*m3+c],Style[Point[bob3],Antialiasing->True]}
}
]
];

(*---------------------------------------------------------*)
obtainScreenPositions[p1_,bbob1_,bbob2_,bbob3_,\[Theta]\[Theta]1Init_,\[Theta]\[Theta]2Init_,\[Theta]\[Theta]3Init_,L1_,L2_,L3_,nBobs_]:=Module[{bob1=bbob1,bob2=bbob2,bob3=bbob3,\[Theta]1InitMouse=\[Theta]\[Theta]1Init,\[Theta]2InitMouse=\[Theta]\[Theta]2Init,\[Theta]3InitMouse=\[Theta]\[Theta]3Init,delx,dely,oldbob2,newbob1,newbob2,newbob3},

Which[nBobs==1,
(
bob1=positionOfNewBob[p1,L1];
\[Theta]1InitMouse=normalizedAngleFromMouseInput[ArcTan[bob1[[1]],bob1[[2]]]]
),
nBobs==2,
(
If[EuclideanDistance[p1,bob1]<EuclideanDistance[p1,bob2],
(
newbob1=positionOfNewBob[p1,L1];
delx=newbob1[[1]]-bob1[[1]];
dely=newbob1[[2]]-bob1[[2]];
bob1=newbob1;
bob2={bob2[[1]]+delx,bob2[[2]]+dely};

\[Theta]1InitMouse=normalizedAngleFromMouseInput[ArcTan[bob1[[1]],bob1[[2]]]];
\[Theta]2InitMouse=normalizedAngleFromMouseInput[ArcTan[bob2[[1]]-bob1[[1]],bob2[[2]]-bob1[[2]]]]
)
,
(
bob2=positionOfNewBob[p1,bob1,L2];
\[Theta]2InitMouse=normalizedAngleFromMouseInput[ArcTan[bob2[[1]]-bob1[[1]],bob2[[2]]-bob1[[2]]]]
)
]
),
True,
(
If[EuclideanDistance[p1,bob1]<EuclideanDistance[p1,bob2],
(
If[EuclideanDistance[p1,bob1]<EuclideanDistance[p1,bob3],
( (*bob1*)

newbob1=positionOfNewBob[p1,L1];
delx=newbob1[[1]]-bob1[[1]];
dely=newbob1[[2]]-bob1[[2]];
bob1=newbob1;
oldbob2=bob2;
bob2={bob2[[1]]+delx,bob2[[2]]+dely};
\[Theta]1InitMouse=normalizedAngleFromMouseInput[ArcTan[bob1[[1]],bob1[[2]]]];

\[Theta]2InitMouse=normalizedAngleFromMouseInput[ArcTan[bob2[[1]]-bob1[[1]],bob2[[2]]-bob1[[2]]]];

delx=bob2[[1]]-oldbob2[[1]];
dely=bob2[[2]]-oldbob2[[2]];
bob3={bob3[[1]]+delx,bob3[[2]]+dely};
\[Theta]3InitMouse=normalizedAngleFromMouseInput[ArcTan[bob3[[1]]-bob2[[1]],bob3[[2]]-bob2[[2]]]]
),
((*bob3*)
newbob3=positionOfNewBob[p1,bob2,L3];
delx=newbob3[[1]]-bob3[[1]];
dely=newbob3[[2]]-bob3[[2]];
bob3={bob3[[1]]+delx,bob3[[2]]+dely};
\[Theta]3InitMouse=normalizedAngleFromMouseInput[ArcTan[bob3[[1]]-bob2[[1]],bob3[[2]]-bob2[[2]]]]
)
]
)
,
(
If[EuclideanDistance[p1,bob2]<EuclideanDistance[p1,bob3],
( (*bob2*)

newbob2=positionOfNewBob[p1,bob1,L2];
delx=newbob2[[1]]-bob2[[1]];
dely=newbob2[[2]]-bob2[[2]];
bob2={bob2[[1]]+delx,bob2[[2]]+dely};
\[Theta]2InitMouse=normalizedAngleFromMouseInput[ArcTan[bob2[[1]]-bob1[[1]],bob2[[2]]-bob1[[2]]]];

bob3={bob3[[1]]+delx,bob3[[2]]+dely};
\[Theta]3InitMouse=normalizedAngleFromMouseInput[ArcTan[bob3[[1]]-bob2[[1]],bob3[[2]]-bob2[[2]]]]
),
((*bob3*)
newbob3=positionOfNewBob[p1,bob2,L3];
delx=newbob3[[1]]-bob3[[1]];
dely=newbob3[[2]]-bob3[[2]];
bob3={bob3[[1]]+delx,bob3[[2]]+dely};
\[Theta]3InitMouse=normalizedAngleFromMouseInput[ArcTan[bob3[[1]]-bob2[[1]],bob3[[2]]-bob2[[2]]]]
)
]
)]
)
];

{bob1,bob2,bob3,\[Theta]1InitMouse,\[Theta]2InitMouse,\[Theta]3InitMouse}
];

(* keep angle in range 0..Pi, -Pi...0 *)
(*---------------------------------------------------------*)
normalizedAngleFromMouseInput[\[Theta]_]:=Module[{},

Which[\[Theta]>=0 && \[Theta]<= Pi/2,\[Theta]+Pi/2,
\[Theta]>Pi/2 && \[Theta]<= Pi,-(Pi/2+(Pi-\[Theta])),
True,Pi/2+\[Theta]
]
];

(*Helper function to determine which bob to move on the screen*)
(*when using eventHandler*)
(*---------------------------------------------------------*)
positionOfNewBob[p1_,L1_]:=Module[{newBob,y,x,eq,pt},

eq=y-p1[[2]]==-p1[[2]]/-p1[[1]] (x-p1[[1]]);
pt=NSolve[{x^2+y^2==L1^2,eq},{x,y}];

If[EuclideanDistance[p1,{x/.pt[[1]],y/.pt[[1]]}]<
EuclideanDistance[p1,{x/.pt[[2]],y/.pt[[2]]}],
(
newBob={x/.pt[[1]],y/.pt[[1]]}
),
(
newBob={x/.pt[[2]],y/.pt[[2]]}
)
];

newBob
];

(*---------------------------------------------------------*)
positionOfNewBob[p1_,bob_,len_]:=Module[{newBob,y,x,eq,pt},
eq=y-p1[[2]]==(bob[[2]]-p1[[2]])/(bob[[1]]-p1[[1]]) (x-p1[[1]]);

pt=NSolve[{(x-bob[[1]])^2+(y-bob[[2]])^2==len^2,eq},{x,y}];

If[EuclideanDistance[p1,{x/.pt[[1]],y/.pt[[1]]}]<
EuclideanDistance[p1,{x/.pt[[2]],y/.pt[[2]]}],
(
newBob={x/.pt[[1]],y/.pt[[1]]}
),
(
newBob={x/.pt[[2]],y/.pt[[2]]}
)
];

newBob
];

(*---------------------------------------------------------*)
solve[nBobs_,\[Theta]1Init_,\[Theta]2Init_,\[Theta]3Init_,\[Theta]1SpeedInit_,\[Theta]2SpeedInit_,\[Theta]3SpeedInit_,m1_,m2_,m3_,L1_,L2_,L3_,g_,from_,to_,c_,accuracyGoal_]:=Module[{numericalSolution,ic,t,ndsolveOptions,x1,x2,x3,x1der,x2der,x3der},

ndsolveOptions = {MaxSteps->Infinity,Method->{"StiffnessSwitching",Method->{"ExplicitRungeKutta",Automatic}},AccuracyGoal->accuracyGoal};

Which[
nBobs==1,
(
ic={x1[0]==\[Theta]1Init ,x1'[0]==\[Theta]1SpeedInit};
numericalSolution=
First@NDSolve[
Flatten[{eqOneBob[L1,m1,g,t,x1,c]==0,ic}],{x1,x1'},{t,from,to},Sequence@ndsolveOptions];

x1=x1/.numericalSolution;
x1der=x1'/.numericalSolution
),
nBobs==2,
(
ic={x1[0]==\[Theta]1Init ,x1'[0]==\[Theta]1SpeedInit,x2[0]==\[Theta]2Init ,x2'[0]==\[Theta]2SpeedInit};
numericalSolution=
First@NDSolve[
Flatten[{
eqOne2Bob[L1,L2,m1,m2,g,t,x1,x2,c]==0,
eqTwo2Bob[L1,L2,m2,g,t,x1,x2,c]==0,ic}],{x1,x1',x2,x2'},{t,from,to},Sequence@ndsolveOptions];
{x1,x2,x1der,x2der}={x1/.numericalSolution,x2/.numericalSolution,x1'/.numericalSolution,x2'/.numericalSolution}
),
True,
(
ic={x1[0]==\[Theta]1Init ,x1'[0]==\[Theta]1SpeedInit,x2[0]==\[Theta]2Init ,x2'[0]==\[Theta]2SpeedInit,x3[0]==\[Theta]3Init ,x3'[0]==\[Theta]3SpeedInit};

numericalSolution=First@NDSolve[
Flatten[{
eqOne3Bob[L1,L2,L3,m1,m2,m3,g,t,x1,x2,x3,c]==0,
eqTwo3Bob[L1,L2,L3,m2,m3,g,t,x1,x2,x3,c]==0,
eqThree3Bob[L1,L2,L3,m3,g,t,x1,x2,x3,c]==0,
ic}],
{x1,x1',x2,x2',x3,x3'},
{t,from,to},
Sequence@ndsolveOptions];

{x1,x2,x3,x1der,x2der,x3der}={x1/.numericalSolution,x2/.numericalSolution,x3/.numericalSolution,x1'/.numericalSolution,x2'/.numericalSolution,x3'/.numericalSolution}
)
];

{x1,x2,x3,x1der,x2der,x3der}
];

(*---------------------------------------------------------*)
makePhasePortrait[nBobs_,\[Theta]1Init_,\[Theta]2Init_,\[Theta]3Init_,\[Theta]1SpeedInit_,\[Theta]2SpeedInit_,\[Theta]3SpeedInit_,m1_,m2_,m3_,L1_,L2_,L3_,g_,c_,maxRunTime_,dt_]:=Module[{\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed,data,imageSize={365,161},t,icPoint,finalPoint,accuracyGoal=4},


{\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed} =solve[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,0,maxRunTime,c,accuracyGoal];

data=Table[{\[Theta]1[t],\[Theta]1Speed[t]},{t,0,maxRunTime,dt}];
icPoint={{Black,Text[Style["(i.c)",12],{\[Theta]1Init,\[Theta]1SpeedInit},{0,1.5}]},{Red,PointSize[0.02],Point[{\[Theta]1Init,\[Theta]1SpeedInit}]}};
finalPoint={Black,PointSize[0.02],Point[{\[Theta]1[maxRunTime],\[Theta]1Speed[maxRunTime]}]};

Show[ListPlot[data,
Joined->True,
PlotStyle->{RGBColor[{250,128,114}/255],Directive[Thickness[0.001]]},
ImageSize->imageSize,
PlotRange->All,
ImageMargins->1,
AspectRatio->.3,
TicksStyle->Directive[7],
AxesLabel->{
Text@Grid[{{Style["\[Theta]",10]},{Style["(rad)",10]}},Alignment->Center,Spacings->{0,0}],
Text@Style[" Overscript[\[Theta], \[Bullet]] (rad/sec)",10]
}
],Graphics@icPoint,Graphics@finalPoint,PlotRange->All]

];

(*Below are the ode equestions. For different number of bobs*)
(*I derived these in a different notebook using basic Lagrangian derivation*)
(*---------------------------------------------------------*)
eqOneBob[L1_,m1_,g_,t_,x1_,c_]:=Module[{},
c (x1^\[Prime])[t]+L1 m1 (g Sin[x1[t]]+L1 (x1^\[Prime]\[Prime])[t])
];

(*---------------------------------------------------------*)
eqOne2Bob[L1_,L2_,m1_,m2_,g_,t_,x1_,x2_,c_]:=Module[{},
c x1'[t]+L2 m2 Sin[x1[t]-x2[t]] (x2^\[Prime])[t]^2+(m1+m2) (g Sin[x1[t]]+L1 (x1^\[Prime]\[Prime])[t])+L2 m2 Cos[x1[t]-x2[t]] (x2^\[Prime]\[Prime])[t]
];

(*---------------------------------------------------------*)
eqTwo2Bob[L1_,L2_,m2_,g_,t_,x1_,x2_,c_]:=Module[{},
c x2'[t]+m2 (g Sin[x2[t]]-L1 Sin[x1[t]-x2[t]] (x1^\[Prime])[t]^2+L1 Cos[x1[t]-x2[t]] (x1^\[Prime]\[Prime])[t]+L2 (x2^\[Prime]\[Prime])[t])
];

(*---------------------------------------------------------*)
eqOne3Bob[L1_,L2_,L3_,m1_,m2_,m3_,g_,t_,x1_,x2_,x3_,c_]:=Module[{},c x1'[t]+g (m1+m2+m3) Sin[x1[t]]+L2 (m2+m3) Sin[x1[t]-x2[t]] (x2^\[Prime])[t]^2+L3 m3 Sin[x1[t]-x3[t]] (x3^\[Prime])[t]^2+L1 m1 (x1^\[Prime]\[Prime])[t]+(m2+m3) (L1 (x1^\[Prime]\[Prime])[t]+L2 Cos[x1[t]-x2[t]] (x2^\[Prime]\[Prime])[t])+L3 m3 Cos[x1[t]-x3[t]] (x3^\[Prime]\[Prime])[t]
];

(*---------------------------------------------------------*)
eqTwo3Bob[L1_,L2_,L3_,m2_,m3_,g_,t_,x1_,x2_,x3_,c_]:=Module[{},
c x2'[t]-L1 (m2+m3) Sin[x1[t]-x2[t]] (x1^\[Prime])[t]^2+L3 m3 Sin[x2[t]-x3[t]] (x3^\[Prime])[t]^2+(m2+m3) (g Sin[x2[t]]+L1 Cos[x1[t]-x2[t]] (x1^\[Prime]\[Prime])[t]+L2 (x2^\[Prime]\[Prime])[t])+L3 m3 Cos[x2[t]-x3[t]] (x3^\[Prime]\[Prime])[t]
];

(*---------------------------------------------------------*)
eqThree3Bob[L1_,L2_,L3_,m3_,g_,t_,x1_,x2_,x3_,c_]:=Module[{},
c x3'[t]+m3 (g Sin[x3[t]]-L1 Sin[x1[t]-x3[t]] (x1^\[Prime])[t]^2-L2 Sin[x2[t]-x3[t]] (x2^\[Prime])[t]^2+L1 Cos[x1[t]-x3[t]] (x1^\[Prime]\[Prime])[t]+L2 Cos[x2[t]-x3[t]] (x2^\[Prime]\[Prime])[t]+L3 (x3^\[Prime]\[Prime])[t])
];

(*---------------------------------------------------------*)
updatePendulumPositions[nBobs_,\[Theta]1_,\[Theta]2_,\[Theta]3_,L1_,L2_,L3_]:=Module[{bob1=0,bob2=0,bob3=0},
(*check how many bobs we are using, and update the solution*)

Which[nBobs==1,
(
bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]}
),

nBobs==2,
(
bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]};
bob2=bob1+{L2 Sin[\[Theta]2],-L2 Cos[\[Theta]2]}
),
True,
(bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]};
bob2=bob1+{L2 Sin[\[Theta]2],-L2 Cos[\[Theta]2]};
bob3=bob2+{L3 Sin[\[Theta]3],-L3 Cos[\[Theta]3]}
)
];

{bob1,bob2,bob3}
];

(*---------------------------------------------------------*)
(*update bob positions when IC changes*)
resetUsingInitialConditions[nBobs_,\[Theta]1Init_,\[Theta]2Init_,\[Theta]3Init_,\[Theta]1SpeedInit_,\[Theta]2SpeedInit_,\[Theta]3SpeedInit_,bbob2_,bbob3_,L1_,L2_,L3_]:=Module[{\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed,bob1,bob2=bbob2,bob3=bbob3},

\[Theta]1=\[Theta]1Init*Pi/180.;
\[Theta]2=\[Theta]2Init*Pi/180.;
\[Theta]3=\[Theta]3Init*Pi/180.;
\[Theta]1Speed=\[Theta]1SpeedInit;
\[Theta]2Speed=\[Theta]2SpeedInit;
\[Theta]3Speed=\[Theta]3SpeedInit;

Which[nBobs==1,
(
bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]}
),
nBobs==2,
(
bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]};
bob2=bob1+{L2 Sin[\[Theta]2],-L2 Cos[\[Theta]2]}
),
True,
(
bob1={L1 Sin[\[Theta]1],-L1 Cos[\[Theta]1]};
bob2=bob1+{L2 Sin[\[Theta]2],-L2 Cos[\[Theta]2]};
bob3=bob2+{L3 Sin[\[Theta]3],-L3 Cos[\[Theta]3]}
)
];

{\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed,bob1,bob2,bob3}
];
(*---------------------------------------------------------*)
makePeKeChart[currentPE_,currentKE_]:=
Module[{g1,peValueAsPercentage,keValueAsPercentage,totalE},

totalE=currentPE+currentKE;

If[totalE<=$MachineEpsilon,
(
peValueAsPercentage=0;
keValueAsPercentage=0

),
(
peValueAsPercentage=Abs[currentPE]/totalE*100;
keValueAsPercentage=currentKE/totalE*100
)
];

g1=Grid[{
{
Text@Style["P.E. (J)",10],Text@Style["K.E. (J)",10]
},
{
Text@Row[{Style[padIt2[currentPE,{7,2}],10]}],
Text@Row[{Style[padIt2[currentKE,{7,2}],10]}]
},
{
Text@Row[{Style[padIt2[peValueAsPercentage,{4,2}],10],Style[" %",10]}],

Text@Row[{Style[padIt2[keValueAsPercentage,{4,2}],10],Style[" %",10]}]
}


},Frame->All,FrameStyle->Directive[Thickness[.005],Gray],Spacings->{0.4,0.2},Alignment->Center
];

Grid[{
{Grid[{
{g1},
{Text@Row[{Style["total ",11],Style[padIt2[totalE,{8,2}],11], Style[" (J)",11]}]}
},Spacings->{0,0},Alignment->Center],Graphics[
{
{Blue,Rectangle[{0,0},{60,peValueAsPercentage}]},
{Red,Rectangle[{65,0},{115,keValueAsPercentage}]},
Text[Style["PE",8],{30,1.1peValueAsPercentage},{0,-1}],
Text[Style["KE",8],{80,1.1*keValueAsPercentage},{0,-1}]
},ImageSize->{70,63},ImageMargins->0,ImagePadding->0,PlotRange->{{-20,125},{-1,140}}
]
}
},Frame->None,Spacings->{0,0},Alignment->Center
]

];

(*Calulate PE*)
(*---------------------------------------------------------*)
pe[nBobs_,m1_,mm2_,mm3_,L1_,LL2_,LL3_,g_,\[Theta]1_,\[Theta]\[Theta]2_,\[Theta]\[Theta]3_]:=Module[{m2,m3,L2,L3,\[Theta]2,\[Theta]3},

Which[
nBobs==1,(m2=0;m3=0;L2=0;L3=0;\[Theta]2=0;\[Theta]3=0),
nBobs==2,(m2=mm2;L2=LL2;m3=0;L3=0;\[Theta]2=\[Theta]\[Theta]2;\[Theta]3=0),
True,(m2=mm2;L2=LL2;\[Theta]2=\[Theta]\[Theta]2;m3=mm3;L3=LL3;\[Theta]3=\[Theta]\[Theta]3);
];

g L1 m1 (1-Cos[\[Theta]1])+g m2 (L1 (1-Cos[\[Theta]1])+L2 (1-Cos[\[Theta]2]))+g m3 (L1 (1-Cos[\[Theta]1])+L2 (1-Cos[\[Theta]2])+L3 (1-Cos[\[Theta]3]))
];

(*calculate K.E.*)
(*---------------------------------------------------------*)
ke[nBobs_,\[Theta]1_,\[Theta]\[Theta]2_,\[Theta]\[Theta]3_,\[Theta]1Speed_,\[Theta]\[Theta]2Speed_,\[Theta]\[Theta]3Speed_,m1_,mm2_,mm3_,L1_,LL2_,LL3_]:=Module[{m2,m3,L2,L3,\[Theta]2,\[Theta]3,\[Theta]2Speed,\[Theta]3Speed},

Which[

nBobs==1,
(m2=0;m3=0;L2=0;L3=0;\[Theta]2Speed=0;\[Theta]3Speed=0;\[Theta]2=0;\[Theta]3=0),

nBobs==2,
(m2=mm2;L2=LL2;\[Theta]2=\[Theta]\[Theta]2;\[Theta]2Speed=\[Theta]\[Theta]2Speed;m3=0;L3=0;\[Theta]3=0;\[Theta]3Speed=0),

True,
(m2=mm2;L2=LL2;\[Theta]2=\[Theta]\[Theta]2;\[Theta]2Speed=\[Theta]\[Theta]2Speed;m3=mm3;L3=LL3;\[Theta]3Speed=\[Theta]\[Theta]3Speed;\[Theta]3=\[Theta]\[Theta]3)
];

1/2 m1 (L1 \[Theta]1Speed)^2+1/2 m2( (L1 \[Theta]1Speed Cos[\[Theta]1]+L2 \[Theta]2Speed Cos[\[Theta]2])^2+(L1 \[Theta]1Speed Sin[\[Theta]1]+L2 \[Theta]2Speed Sin[\[Theta]2])^2)+
1/2 m3( (L1 \[Theta]1Speed Cos[\[Theta]1]+L2 \[Theta]2Speed Cos[\[Theta]2]+L3 \[Theta]3Speed Cos[\[Theta]3])^2+(L1 \[Theta]1Speed Sin[\[Theta]1]+L2 \[Theta]2Speed Sin[\[Theta]2]+L3 \[Theta]3Speed Sin[\[Theta]3])^2)
];

(*---------------------------------------------------------*)
update[nBobs_,\[Theta]1Init_,\[Theta]2Init_,\[Theta]3Init_,\[Theta]1SpeedInit_,\[Theta]2SpeedInit_,\[Theta]3SpeedInit_,m1_,m2_,m3_,L1_,L2_,L3_,g_,currentTime_,c_,maxRunTime_,dt_,showPhase_]:=Module[{currentPE,currentKE,bob1,bob2,bob3,\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed,phasePortraitPlot={},accuracyGoal=7},

{\[Theta]1,\[Theta]2,\[Theta]3,\[Theta]1Speed,\[Theta]2Speed,\[Theta]3Speed} =solve[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,currentTime,currentTime+1,c,accuracyGoal];

currentPE=pe[nBobs,m1,m2,m3,L1,L2,L3,g,\[Theta]1[currentTime], \[Theta]2[currentTime],\[Theta]3[currentTime]];
currentKE=ke[nBobs,\[Theta]1[currentTime],\[Theta]2[currentTime],\[Theta]3[currentTime],\[Theta]1Speed[currentTime],\[Theta]2Speed[currentTime],\[Theta]3Speed[currentTime],m1,m2,m3,L1,L2,L3];

If[showPhase,
phasePortraitPlot=makePhasePortrait[nBobs,\[Theta]1Init,\[Theta]2Init,\[Theta]3Init,\[Theta]1SpeedInit,\[Theta]2SpeedInit,\[Theta]3SpeedInit,m1,m2,m3,L1,L2,L3,g,c,maxRunTime,dt]
];

{bob1,bob2,bob3}=updatePendulumPositions[nBobs,\[Theta]1[currentTime],\[Theta]2[currentTime],\[Theta]3[currentTime],L1,L2,L3];

{currentPE,currentKE,phasePortraitPlot,bob1,bob2,bob3,\[Theta]1,\[Theta]1Speed}

])
]