;+
; NAME:
;	lsqNormalFit
; PURPOSE:
;	Find least squares fit to a data series x,y
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
	FUNCTION lsqNormalFit, xin, yin , $
		plotx	= plotx , $
		oplotx	= oplotx, $
		worst	= worst , $
		xfit	= xfit	, $
		yfit	= yfit	, $
		_extra	= _extra
; INPUTS:
;	x,y	arrays to be fitted
; OPTIONAL INPUT PARAMETERS:
;	/plot	data points and fit are plotted
;	/oplot	the fit is plotted using the current user scale
;	/worst	will also plot the worst possible fit
; OUTPUTS:
;	A	fltarr(5,2)
;		A[*,0] represents the best fit; A[*,1] represents the worst fit
;
;	The equation for the straight line fit:
;		A[0,i]*(x-A[1,i])+A[2,i]*(y-A[3,i]) = 0
;
;	A[1,i] is the average of the x-array
;	A[3,i] is the average of the y-array
;
;	A[4,i] is the residual average distance of the `normalized' data points to the fit.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; PROCEDURE:
;	The x and y arrays are `normalized' by dividing the arrays by their
;	respective standard deviations. The resulting arrays are fitted by
;	minimizing the perpendicular distance to a straight line fit.
;	The result is symmetric: interchanging x and y arrays will give the same fit.
; MODIFICATION HISTORY:
;	JUN-1994, Paul Hick (UCSD/CASS)
;	AUG-1999, Paul Hick (UCSD/CASS); added title keyword
;	APR-2000, Paul Hick (UCSD/CASS); added _extra keyword
;	FEB-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu); changed the main plotting command to
;		use the input arrays directly for plotting.
;-

InitVar, plotx	, /key
InitVar, oplotx , /key
InitVar, worst	, /key

n = n_elements(xin)
IF n_elements(yin) NE N THEN message, 'both arrays must have same size

x	= xin
y	= yin

Vx  = total(x)/n
Vy  = total(y)/n

xx	= x-Vx
yy	= y-Vy
Dxx = sqrt( total(xx*xx)/n )
Dyy = sqrt( total(yy*yy)/n )

x	= x/Dxx
y	= y/Dyy

Tx  = total(x)/n
Ty  = total(y)/n

xx  = x-Tx
yy  = y-Ty
Rx  = total(xx*xx)/n
Ry  = total(yy*yy)/n
Rxy = total(xx*yy)/n

Rp  = (Rx+Ry)/2.
Rm  = (Rx-Ry)/2.

Zxy = (Rm*Rm+Rxy*Rxy) > 0
Zxy = sqrt(Zxy)
;print, rp,'  ',rm,'  ',zxy

; The equation for the fit is A*(x-Tx)+B*(y-Ty)=0

A1	= -Tx*Rxy+Ty*(Rm+Zxy)		; Worst fit
B1	=  Ty*Rxy+Tx*(Rm-Zxy)
F1	= sqrt(Rp+Zxy)				; Root of average of dist^2

A2	= -Tx*Rxy+Ty*(Rm-Zxy)		; Best fit
B2	=  Ty*Rxy+Tx*(Rm+Zxy)
F2	= sqrt( (Rp-Zxy) > 0 )		; Root of average of dist^2

d	= sort(x)					; Do not sort the arrays prior to the calculation.
x	= x[d]						; This destroys the symmetry between x and y.
y	= y[d]

;if plotx then plot, Dxx*x,Dyy*y, /ynozero, psym=2, _extra=_extra
IF plotx THEN plot, xin[d], yin[d], /ynozero, psym=2, _extra=_extra

IF plotx OR oplotx THEN BEGIN
	d  = 3.*sqrt(Rp)
	s  = sqrt(A2*A2+B2*B2)
	dx =  d*(b2/s)
	dy = -d*(a2/s)
	oplot, Dxx*(Tx+dx*[-1,1]), Dyy*(Ty+dy*[-1,1])
	IF worst AND (A1 NE 0 OR B1 NE 0) THEN BEGIN
		s  = sqrt(A1*A1+B1*B1)
		dx =  d*(b1/s)
		dy = -d*(a1/s)
		oplot, Dxx*(Tx+dx*[-1,1]), Dyy*(Ty+dy*[-1,1])
	ENDIF
	catch, ierr
	IF ierr EQ 0 THEN s = poly_fit(Dxx*x,Dyy*y,1,d)
	catch, /cancel
	;oplot, Dxx*x, d, linestyle=3
ENDIF

yfit = Vy-(A2/B2)*(Dyy/Dxx)*(x-Vx)
xfit = Vx-(B2/A2)*(Dxx/Dyy)*(y-Vy)

RETURN, [[A2/Dxx,Vx,B2/Dyy,Vy,F2],[A1/Dxx,Vx,B1/Dyy,Vy,F1]]  &  END