;+
; NAME:
;	lsqQuadraticFit
; PURPOSE:
;	Fits points to a quadratic polynomial
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
	FUNCTION lsqQuadraticFit, x, y, chisq=chisq
; INPUTS:
;	x		array[n,m] or array[2,n,m]; type: any
;				m sequences of n x-coordinates
;				if m=1 then the last dimension may be absent
;				if y not specified then 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 sequences of n y-coordinates matching the x-array
; OUTPUTS:
;	R		array[3] or array[3,m]; type: double
;				the 3 fitted constants in
;					y = R[0]+R[1]*x+R[2]*x^2
; OPTIONAL OUTPUT PARAMETERS:
;	chisq=chisq
;			scalar, or array[m]; type: double
;				residual chisq.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	SubArray, SyncArgs
; PROCEDURE:
; >	If both x and y is called then it is allowed to specify only
;	one set of x-values, x[n] and multiple sets of y-values, y[n,m]
;	In that case each y-series is combined with x.
; >	The same can be accomplished with the IDL function poly_fit.
;	polyfit uses a matrix inversion; here we use the analytic solution.
; MODIFICATION HISTORY:
;	JAN-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

CASE IsType(y,/defined) OF
0: BEGIN
	xx = SubArray(x, element=0)
	yy = SubArray(x, element=1)
END
1: BEGIN
	xx = x
	yy = y
	SyncArgs, xx, yy
END
ENDCASE

xx  = double(xx)
yy  = double(yy)

sz  = size(xx)
N   = sz[1]
M   = sz[sz[0]+2]/sz[1]
one = replicate(1,N)

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

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

xscale = xmean		; Only array structure matters
yscale = ymean

FOR i=0,M-1 DO BEGIN
	xscale[i] = max(abs(xx[*,i]))
	yscale[i] = max(abs(yy[*,i]))
ENDFOR


xx /= one#xscale
yy /= one#yscale
x2 = xx*xx

Sx4 = total(x2*x2,1)/N
Sx3 = total(x2*xx,1)/N
Sx2 = total(x2	 ,1)/N
Sx1 = total(xx	 ,1)/N

Sy1 = total(yy	 ,1)/N
Sxy = total(xx*yy,1)/N
Sx2y = total(x2*yy,1)/N

D1203 = Sx1*Sx2-Sx3

D2314 = Sx2*Sx3-Sx1*Sx4
D0422 = Sx4-Sx2*Sx2

D1322 = Sx1*Sx3-Sx2*Sx2
D0211 = Sx2-Sx1*Sx1

D2433 = Sx2*Sx4-Sx3*Sx3

; The three expressions for the denominator are all the same

D = D2433+Sx1*D2314+Sx2*D1322
;D = Sx1*D2314+Sx2*D0422+Sx3*D1203
;D = Sx2*D1322+Sx3*D1203+Sx4*D0211

B = (Sy1*D2314+Sxy*D0422+Sx2y*D1203)/D
C = (Sy1*D1322+Sxy*D1203+Sx2y*D0211)/D

A = (Sy1*D2433+Sxy*D2314+Sx2y*D1322)/D
;A = Sy1-B*Sx1-C*Sx2

chisq = one#A+(one#B)*xx+(one#C)*x2-yy
chisq = sqrt( (chisq*chisq)/N )

xmean /= xscale
ymean /= yscale

AA = yscale*( ymean+xmean*(xmean*C-B)+A )

yscale /= xscale

BB = yscale*(B-2*xmean*C)
CC = yscale*(C/xscale)

RETURN, transpose( [[AA],[BB],[CC]] )  &  END