function Par = KarimakiOrig(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 is the original version, without corrective step
%
%--------------------------------------------------------------------------

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);

ZXY1 = [Z XY ones(n,1)];
M = ZXY1' * ZXY1;
Cxx = det(M(2:2:4,2:2:4));
Cyy = det(M(3:4,3:4));
Cxy = det(M(2:2:4,3:4));
Czx = det(M(1:3:4,2:2:4));
Czy = det(M(1:3:4,3:4));
Czz = det(M(1:3:4,1:3:4));
frac = 2*(Czz*Cxy - Czx*Czy)/(Czz*(Cxx-Cyy)-Czx*Czx+Czy*Czy);
phi = atan(frac)/2;
if (Czz*Cxy - Czx*Czy)*sin(2*phi)<0
   phi = phi+pi/2;
end
sp = sin(phi);  cp = cos(phi);
kappa = (Czx*sp - Czy*cp)/Czz;
mu = (-M(1,4)*kappa + M(2,4)*sp - M(3,4)*cp)/n;
rho = 2*kappa/sqrt(1-4*kappa*mu);
d = 2*mu/(1+sqrt(1-4*kappa*mu));

Par = zeros(1,3);
Par(1) = (d+1/rho)*sin(phi);
Par(2) =-(d+1/rho)*cos(phi);
Par(1:2) = Par(1:2) + pole;
Par(3) = 1/abs(rho);

end   %  KarimakiOrig