function Par = RiemannSWFL(XY)

%--------------------------------------------------------------------------
%  
%     Riemann circle fit, the SWFL version
%   A. Strandlie, J. Wroldsen, R. Fruhwirth, and B. Lillekjendlie,
%   "Particle tracks fitted on the Riemann sphere",
%   Computer Physics Commun., Vol. 131, pp. 95-108,  (2000)
%
%     Input:  XY(n,2) is array of coordinates of n points <XY(i,1),XY(i,2)>
%
%     Output: Par = [a b R] is the fitting circle:
%                           center (a,b) and radius R
%
%     This is the original version
%
%--------------------------------------------------------------------------

n = size(XY,1);      % number of data points
centroid = mean(XY);
X = XY(:,1) - centroid(1);
Y = XY(:,2) - centroid(2);
Z = X.*X + Y.*Y;
factor = 2*sqrt(sum(Z)/n);
X = X/factor;
Y = Y/factor;
Z = Z/factor/factor;
XP = X./(Z+1);
YP = Y./(Z+1);
ZP = Z./(Z+1);
W = (Z+1).*(Z+1);
R = [XP YP ZP];
Rbar = W'*R/sum(W);
Rcen = R - repmat(Rbar,n,1);
RcenW = Rcen; 
for k=1:3
    RcenW(:,k) = RcenW(:,k).*W; 
end
Scatter = RcenW'*Rcen;
[E,D] = eig(Scatter);
[Dsort,ID] = sort(diag(D));
if (Dsort(1)<0) 
    disp('Error in SWFL: the smallest e-value is negative...')
end
V = E(:,ID(1))';
c = -Rbar*V';
Par = -V(1:2)/2/(c+V(3));
Par = [Par , sqrt(Par*Par'-c/(c+V(3)))];
Par = Par*factor;
Par(1:2) = Par(1:2) + centroid;

end   %  RiemannSWFL