function [xc,yc,R] = cirfit(x,y)
%  x=[x1,x2,x3,...,xn], y=[y1,y2,y3,...,yn]
n=length(x);
xx=x.*x;
yy=y.*y;
xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy)...
    sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -0.5*a(1);
yc = -0.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));
end