FUNCTION lsqCircleFit, x, y, chisq=chisq
;+
; NAME:
;	lsqCircleFit
; PURPOSE:
;	Least square fitting a circle to a set of points
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
;	R = lsqCircleFit(x,y)
; INPUTS:
;	x		array[n,m] or array[2,n,m]; type: any
;				m sets of n x-coordinates
;				if y is not specified than x must be of the form x[2,n,m]
;				with x[0,*,*] the x-values and x[1,*,*] the y-values. 
;	y		array[n,m]; type: any
;				m sets of n y-coordinates matching the x-array
; OUTPUTS:
;	R		array[3] or array[3,m]; type: float or double
;				x-coordinate and y-coordinate of center, followed by the radius
; OPTIONAL OUTPUT PARAMETERS:
;	chisq=chisq
;			scalar or array[m]; type: float
;				residual 'chi-square' = sqrt ( sqrt( (1/N)*Sum( ((x-x0)^2+(y-y0)^2-R^2)^2 ) ) )
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	SubArray, SyncArgs
; PROCEDURE:
; >	If both x and y is specified than it is allowed to specify only one set of x-value
;	x[n] with multiple sets of y-values y[n,m]. Each of the sets of y-values is combined
;	with the set of x-values.
; >	Standard least square fit.
;	Each row in the input arrays x and y is fitted to a circle. A separate center and radius
;	is returned for each of the rows.
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS)
;	FEB-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu), added SyncArgs call.
;-

IF n_elements(y) EQ 0 THEN BEGIN
	xx = SubArray(x, element=0)
	yy = SubArray(x, element=1)
ENDIF ELSE BEGIN
	xx = x
	yy = y
	SyncArgs, xx, yy
ENDELSE

sz  = size(xx)
N   = sz[1]
one = replicate(1,N)

xmean = total(xx,1)/N
ymean = total(yy,1)/N

xx = xx-one#xmean
yy = yy-one#ymean




scale = max(abs([xx,yy]))
xx = xx/scale
yy = yy/scale

x2 = xx*xx
y2 = yy*yy

Sx3  = total(x2*xx,1)/N
Sx2y = total(x2*yy,1)/N
Sxy2 = total(xx*y2,1)/N
Sy3  = total(yy*y2,1)/N

Sx2  = total(x2   ,1)/N
Sxy  = total(xx*yy,1)/N
Sy2  = total(y2   ,1)/N

Sx   = total(xx   ,1)/N
Sy   = total(yy   ,1)/N

Dx   = Sx2-Sx*Sx
Dy   = Sy2-Sy*Sy
Dxy  = Sxy-Sx*Sy

Sr2  = Sx2+Sy2

x0 = (Sx3+Sxy2)*Dy-(Sy3+Sx2y)*Dxy-Sr2*(Sx*Sy2-Sy*Sxy)
y0 = (Sy3+Sx2y)*Dx-(Sx3+Sxy2)*Dxy-Sr2*(Sy*Sx2-Sx*Sxy)

D  = Dx*Dy-Dxy*Dxy
x0 = 0.5*x0/D
y0 = 0.5*y0/D
R  = sqrt( x0*x0+y0*y0+Sr2-2*(x0*Sx+y0*Sy) > 0 )

; Some care must be taken in calculating the residual differences between
; data points and fit, since for a very tight fit it involves calculating the
; the differences between large values.

dx = xx-one#x0
dy = yy-one#y0
dr = sqrt( dx*dx+dy*dy )

tmp = one#R
dx = dx-tmp*(dx/dr)
dy = dy-tmp*(dy/dr)
chisq = dx*dx+dy*dy

chisq = sqrt( total( chisq*chisq, 1 ) / N )	; Estimate for difference in dR^2
chisq = scale*sqrt(chisq)					; Estimate for difference in dR

x0 = xmean+scale*x0
y0 = ymean+scale*y0
R  = R*scale

;plot, x,y
;angle = findgen(360)/!radeg
;oplot, x0+R*cos(angle), y0+R*sin(angle)

RETURN, transpose([[x0], [y0], [R]])  &  END