; NAME :	planefit
; PURPOSE :
;	Fit the data by a plane and give the fitted data.
;	The result can be later used to subtract the linear trend.
; CATEGORY :
; 	image processing
; CALLING SEQUENCE :
;	fittedplane = planefit(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 linear function,
;	        fittedplane = coeff(0)*x + coeff(1)*y + coeff(2)
; OPTIONAL OUTPUTS :
;	none
; COMMON BLOCKS : 	none
; SIDE EFFECTS :	none
; RESTRICTIONS :	none
; PROCEDURE :	linear fit by least square procedure
;
; MODIFICATION HISTORY :
;	T.Sakurai, 1993
;
;------------------------------------------------------------------------
;
function planefit, array, where=where, coeff=coeff
;
; assume a linear function z(x,y)
; z = a*x + b*y +c
;
; the least square solution for (a,b,c) is obtained by solving
;
; [x x] [x y] [x 1]   a    [x z]
; [y x] [y y] [y 1] * b  = [y z]
; [1 x] [1 y] [1 1]   c    [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
	one = fltarr(nx,ny)+1
	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))
	x = x-xmean
	y = y-ymean

	x1 = total(x(where))
	y1 = total(y(where))
	z1 = total(z(where))
	xx = total(x(where)^2)
	yy = total(y(where)^2)
	xy = total(x(where)*y(where))
	xz = total(x(where)*z(where))
	yz = total(y(where)*z(where))
	ones = total(one(where))

	mat = [ [ xx, xy, x1], [xy, yy, y1], [x1, y1, ones] ]
	bb  = [xz, yz, z1]
	ludcmp, mat, index, dd
	lubksb, mat, index, bb

	a=bb(0)
	b=bb(1)
	cc=bb(2)

	z = a*x + b*y + cc
; 		mean values are subtracted from x and y
	c = a*xmean + b*ymean +cc
	coeff=[a,b,c]

	return, z
	end