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