function Par = KarimakiAlg(XY,pole)

%--------------------------------------------------------------------------
%  
%    Karimaki circle fit
%   V. Karimaki,
%   "Effective circle fitting for particle trajectories",
%    Nucl. Instr. Meth. Phys. Res. A., Vol. 305, pp. 187--191, (1991)
%
%     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 version is modified from the original one; it uses algebraic parameters
%        No corective steps
%
%--------------------------------------------------------------------------

n = size(XY,1);      % number of data points
XY = XY - repmat(pole,n,1);   % transfer coordinates to the pole
Z = XY(:,1).*XY(:,1) + XY(:,2).*XY(:,2);

ZXY = [Z XY];
means = mean(ZXY);
ZXY0 = ZXY - repmat(means,n,1);
M = ZXY0' * ZXY0;
C = [det(M(1:2,1:2)) det(M(1:2,1:2:3)); det(M(1:2,1:2:3)) det(M(1:2:3,1:2:3))];
[E,D] = eig(C);
[Dsort,ID] = sort(diag(D));
if (Dsort(1)<0) 
    disp('Error in Karimaki: the smallest e-value is negative...')
end
BC = E(:,ID(1))';
A = -BC*M(2:3,1)/M(1,1);
Dp = -A*means(1) - BC*means(2:3)';

Par = zeros(1,3);
Par(1:2) = -BC/2/A;
Par(3) = sqrt(norm(BC)^2-4*A*Dp)/2/abs(A);
Par(1:2) = Par(1:2) + pole;

end   %  KarimakiAlg