function [xc,yc,R,AvgAngSpeed] = AngVelRadRot(FileName)

clc
close all

%% parameters

% plot font size
FontSize = 14;
figure;hold on;

% Objective funtion tolerance value
ObjFunTol = 1;
MinArcLength = pi/2;

% initialize arrays, variables
SelectTimes = [];
xc = [];
yc = [];
R =[];
AvgAngSpeed = [];

%% measurement
% load the front information at all times
load(FileName,'Front_All','Time_All')

% get the centroid at all snapshots
[x_c, y_c, ~, ~] = GetAreaCentroid(Front_All);

% Estimating the circular motion initiation
[pks,locs] = findpeaks(x_c);

if isempty(pks)
  plot(x_c,y_c,'o');
  return
end

if length(pks)>1
  
  MinObjFun = ObjFunTol;
  index = [];
  for j=1:length(pks)
    [~,~,R,ObjFun] = IsCircle(x_c(locs(j):end),y_c(locs(j):end),x_c(locs(j)),y_c(locs(j)));
    % check for the actual arc length of the fitted data
    xdiff = diff(x_c(locs(j):end));
    ydiff = diff(y_c(locs(j):end));
    dist = sqrt(xdiff.^2+ydiff.^2);
    ArcLength = sum(acos(1-0.5*(dist/R).^2));
    
    if ObjFun<MinObjFun && ArcLength>MinArcLength && isreal(ArcLength)
      MinObjFun = ObjFun;
      index = j;
    end
  end
  if isempty(index)
    plot(x_c,y_c,'o');
  else
    SelectTimes = Time_All(locs(index):end);
  end
  
else
  
  [~,~,R,ObjFun] = IsCircle(x_c(locs:end),y_c(locs:end),x_c(locs),y_c(locs));
  % check for the actual arc length of the fitted data
  xdiff = diff(x_c(locs:end));
  ydiff = diff(y_c(locs:end));
  dist = sqrt(xdiff.^2+ydiff.^2);
  ArcLength = sum(acos(1-0.5*(dist/R).^2));  
  
  if ObjFun<ObjFunTol && ArcLength>MinArcLength && isreal(ArcLength)
    SelectTimes = Time_All(locs(1):end);
  else
    plot(x_c,y_c,'o');
  end
  
end

if isempty(SelectTimes)
    R = [];
    return
end

% original data indices related to the selected times
ind = zeros(1,length(SelectTimes));
% simulation dt
dt = Time_All(2)-Time_All(1);

for i=1:length(SelectTimes)
  ind(i) = find(abs(Time_All-SelectTimes(i)) < 0.5*dt, 1 );
end 

% fit the centroid data to a circle
Par = CircleFitByPratt([x_c(ind) y_c(ind)]);
xc = Par(1);yc = Par(2);R = Par(3);

% angular velocity
xdiff = diff(x_c(ind));
ydiff = diff(y_c(ind));

dist = sqrt(xdiff.^2+ydiff.^2);
AvgSpeed = sum(dist) / (max(Time_All(ind))-min(Time_All(ind)));
AvgAngSpeed = AvgSpeed/R;

%% plot all
% define the fitted circle for visualization
th = linspace(0,2*pi,360)';
xe = R*cos(th)+xc;
ye = R*sin(th)+yc;
plot(x_c(ind),y_c(ind),'o');
plot(x_c(1:ind-1),y_c(1:ind-1),'.');
plot([xe;xe(1)],[ye;ye(1)],'r-','linewidth',3);
plot(xc,yc,'rx','markersize',10)
set(gca,'FontSize',FontSize)
axis equal tight
