; NAME :	quadfit
; PURPOSE :
;	Fit the data by a quadratic function and give the fitted data.
;	The result can be later used to subtract the quadratic trend.
; CATEGORY :
; 	image processing
; CALLING SEQUENCE :
;	fitted = quadfit(array, where=where, coeff=coeff)
; INPUTS :
;	array = 2-d array
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
;	where  = index array that points the locations where the data
;	         are valid
;	coeff = coefficient of fitted quadratic function,
;	        fitted = coeff(0)*x*x + coeff(1)*y*y + coeff(2)*x*y
;                    + coeff(3)*x + coeff(4)*y +coeff(5)
; OPTIONAL OUTPUTS :
;	none
; COMMON BLOCKS : 	none
; SIDE EFFECTS :	none
; RESTRICTIONS :	none
; PROCEDURE :	quadratic by least square procedure
;
; MODIFICATION HISTORY :
;	T.Sakurai, 1993
;
;------------------------------------------------------------------------
;
function quadfit, array, where=where, coeff=coeff
;
;
; assume a quadratic function z(x,y)
; z = a*x**2 + b*y**2 + c*x*y + d*x + e*y +f
;
; the least square solution for (a,b,c,d,e,f) is obtained by solving
;
; [x^2 x^2] [x^2 y^2] [x^2 xy] [x^2 x] [x^2 y] [x^2 1]   a    [x^2 z]
; [y^2 x^2] [y^2 y^2] [y^2 xy] [y^2 x] [y^2 y] [y^2 1]   b    [y^2 z]
; [xy  x^2] [xy  y^2] [xy  xy] [xy  x] [xy  y] [xy  1]*  c  = [xy  z]
; [x   x^2] [x   y^2] [x   xy] [x   x] [x   y] [x   1]   d    [x   z]
; [y   x^2] [y   y^2] [y   xy] [y   x] [y   y] [y   1]   e    [y   z]
; [1   x^2] [1   x^2] [1   xy] [1   x] [1   y] [1   1]   f    [1   z]
;
; where [ ] means summation
;
	nsize = size(array)
	nx = fix(nsize(1))
	ny = fix(nsize(2))
	if (n_elements(nsize) eq 4 ) then nz = fix(nsize(3)) else nz=1
	if (not keyword_set(where)) then where=lindgen(nx,ny)

	z = array
	x = z
	y = z
	for j=0, ny-1 do x(*,j)=findgen(nx)
	for i=0, nx-1 do y(i,*)=findgen(ny)
	xmean = total(x)/(float(nx)*float(ny))
	ymean = total(y)/(float(nx)*float(ny))
	xd = x-xmean
	yd = y-ymean

	xy  = fltarr(5,5)
	xyz = fltarr(3,3)
	mat = fltarr(6,6)
	bb  = fltarr(6)

	for j=0,4 do for i=0,4-j do  $
		xy (i,j) = total( xd(where)^i * yd(where)^j )
	for j=0,2 do for i=0,2-j do  $
		xyz(i,j) = total (xd(where)^i * yd(where)^j * z(where))
	
	mapx = [ 2, 0, 1, 1, 0, 0 ]
	mapy = [ 0, 2, 1, 0, 1, 0 ]

	for j=0,5 do for i=0,5 do $
		mat(i,j)  = xy( mapx(i)+mapx(j), mapy(i)+mapy(j) )
	for j=0,5 do $
		bb(j) = xyz( mapx(j), mapy(j) )

	ludcmp, mat, index, dd
	lubksb, mat, index, bb

	a=bb(0)
	b=bb(1)
	c=bb(2)
	dd=bb(3)
	ed=bb(4)
	fd=bb(5)

	z = a*xd^2 + b*yd^2 + c*xd*yd + dd*xd + ed*yd +fd
	d = dd - 2*a*xmean -c*ymean
	e = ed - 2*b*xmean -c*xmean
	f = fd + a*xmean^2 + b*ymean^2 + c*xmean*ymean - dd*xmean - ed*ymean
;	z = a*x^2 + b*y^2 + c*x*y + d*x + e*y +f
	coeff=[a,b,c,d,e,f]
	return, z
	end