SetDirectory[NotebookDirectory[]];

dt=; (* s, time interval between frames*)

dataDirectory="";

lastFrame=Length[FileNames[dataDirectory<>"/*.tif"]]-1;

images=Table[Import[dataDirectory<>"/t"<>IntegerString[i,10,4]<>".tif"],{i,0,lastFrame}];

findAngle[image_]:=Module[{morphed,pts1,pts2,pts},
morphed=MorphologicalComponents[image];
pts1=N@Reverse[Mean[Position[morphed,1,2]]];
pts2=N@Reverse[Mean[Position[morphed,2,2]]];
pts=SortBy[{pts1,pts2},#[[2]]&];
ArcTan[((#[[1]])/(#[[2]]))]&@(pts[[2]]-pts[[1]])
];

θList=findAngle/@images;
timeList=Range[0,Length[θList]-1]dt;

θvsTime=Transpose[{timeList,θList}];

θvsTimePlot=ListLinePlot[θvsTime,Axes->False,Frame->True,ImageSize->600,FrameStyle->Directive[18,Black,Thick],FrameLabel->{"Time (s)","θ"},PlotStyle->Blue]

Export["Theta vs Time.png",θvsTimePlot,ImageSize->600];

Clear[θ0,b,ω0,t,ϕ];
fit=NonlinearModelFit[θvsTime,θ0 (E^(-((b ω0)/2)t))Cos[ω0 t+ϕ],{θ0,b,ω0,ϕ},t];
fit["ParameterTable"]
fitω0=fit["ParameterTableEntries"][[3,1]]; (* (s^(-1)), fit natural frequency *)
text="Oscillation Period = "<>ToString[Round[((2π)/(fitω0)),0.0001]]<>" s";
Print[text];

Export["Oscillation Period.txt",text];

comparisonPlot=Show[{θvsTimePlot,Plot[fit[t],{t,0,Max[timeList]},PlotStyle->RGBColor[1,0,0,0.5]]}]

Export["Theta vs Time with Fitting.png",comparisonPlot,ImageSize->600];

ωList=Differences[θList]/dt; (* (s^(-1)), angular velocity *)
meanθList=MovingAverage[θList,2];

ωvsθ=Transpose[{meanθList,ωList}];

ωvsθimePlot=ListLinePlot[ωvsθ,Axes->False,Frame->True,ImageSize->600,FrameStyle->Directive[18,Black,Thick],FrameLabel->{"θ (°)"},PlotStyle->Blue]

Export["Phase Diagram.png",ωvsθimePlot,ImageSize->600];