function [Par,Var] = GGS(XY)

%--------------------------------------------------------------------------
%  
%     Algebraic circle fit by Gander,Golub,Strebel (with constraint A^2+B^2+C^2+D^2=1)
%      W. Gander, G.H. Golub, R. Strebel, "Least squares fitting of circles and ellipses",
%      BIT, Vol. 34, pages 558-578, (1994)
%
%     Input:  XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%--------------------------------------------------------------------------

Z = XY(:,1).*XY(:,1) + XY(:,2).*XY(:,2);
ZXY1 = [Z XY ones(size(XY,1),1)];

[U,S,V]=svd(ZXY1,0);
A = V(:,4);
Par = zeros(1,3);
Par(1:2) = -(A(2:3))'/A(1)/2;
Par(3) = sqrt(A(2)*A(2)+A(3)*A(3)-4*A(1)*A(4))/abs(A(1))/2;

end   %  GGS