function Par = Nievergelt(XY)

%--------------------------------------------------------------------------
%  
%     Algebraic circle fit by Nievergelt  (with constraint A^2+B^2+C^2=1)
%      Y. Nievergelt, "Hyperspheres and hyperplanes fitted seamlessly
%                      by algebraic constrained total least-squares",
%      Linear Algebra Appl., Vol. 331, pages 43-59, (2001)
%
%     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
%
%--------------------------------------------------------------------------

centroid = mean(XY);   % the centroid of the data set

X = XY(:,1) - centroid(1);  %  centering data
Y = XY(:,2) - centroid(2);  %  centering data
Z = X.*X + Y.*Y;
Zmean = mean(Z);
Z0 = Z-Zmean;
ZXY = [Z0 X Y];
[U,S,V]=svd(ZXY,0);
A = V(:,3);
A = [A ; -Zmean*A(1)];
Par = zeros(1,3);
Par(1:2) = -(A(2:3))'/A(1)/2 + centroid;
Par(3) = sqrt(A(2)*A(2)+A(3)*A(3)-4*A(1)*A(4))/abs(A(1))/2;

end   %  Nievergelt