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

%  corrections

u = 1+rho*d;
Sa = M(2,4)*sp-M(3,4)*cp;
Sb = M(2,4)*cp+M(3,4)*sp;
cc = cp*cp;  ss = sp*sp;  cs = cp*sp;
Sg = (ss-cc)*M(2,3)+cs*(M(2,2)-M(3,3));
Sd = sp*M(1,2)-cp*M(1,3);
Sh = cp*M(1,2)+sp*M(1,3);
Se = ss*M(2,2)+cc*M(3,3)-2*cs*M(2,3);
Sf = cc*M(2,2)+ss*M(3,3)+2*cs*M(2,3);

W = zeros(3,3);
W(1,1) = M(1,1)/4-d*(Sd-d*(Se+M(1,1)/2-d*(Sa-d*n/4)));
W(2,2) = u*u*Sf;
W(3,3) = rho*(rho*Se-2*u*Sa)+u*u*n;
W(1,2) = -u*(Sh/2-d*(Sg-d*Sb/2));  W(2,1) = W(1,2);
W(1,3) = rho*(d*Se-Sd/2)+u*M(1,4)/2-d*((2*u+rho*d)*Sa-u*d*n)/2;  W(3,1) = W(1,3);
W(2,3) = u*(rho*Sg-u*Sb);  W(3,2) = W(2,3);

V = inv(W);

K = 2*u*Se - rho*Sd - d*(1+u)*Sa;

Drho = -(d*V(1,1) + rho*V(1,3))*K/2;
Dphi = -(d*V(1,2) + rho*V(2,3))*K/2;
Dd   = -(d*V(1,3) + rho*V(3,3))*K/2;

rho = rho + Drho;
phi = phi + Dphi;
d   = d   + Dd;

%  end of corrections

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   %  KarimakiOrigCorr