function chebyshfit, raw,mask,n=n,maxn=maxn,debug=debug
;+
;NAME: 	chebyshfit
;PURPOSE:  
;	An interpolating function implementing a 2D chebyshef
;	polynomial interpolation.
;	Given an image data array(2D) and a mask(of the same size)
;	which identifies a bad pixel on the image with a 0 and a 
;	good pixel with a 1, this function attempts to correct the
;       bad pixel at the center of the grid.  
;EXPLANATION:
;	Use curve fitting to find a function that best defines the 
;	values of the good pixels in the centered row and another 
;	function that best defines the values in the centered column. 
;	Plug in the centered coordinates to these two functions to
;	obtain two interpolated values. The result is the average of 
;	those numbers.
;	In general, all functions can be written in the following terms.
;		f(x) = sum(Ti*Ci)  - .5*Co,
;	where i = 0,1,...infinity, Ti`s are Chebyshef polynomials,
;	Ci's are constants.  
;	Chebyshef Polynomials:
;	To(x) = 1, T1(x) = x, Tn+1(x) = 2*x*Tn - Tn-1 
;	Let (xo,yo) be the coordinates of the centered pixel.
;	Let    x = a vector of the x coordinates of the good pixels,
;	       z = a vector of the image data corresponding to (x,yo), 
;	       coefs = a vector containing Ci's.
;       We obtain the following linear systems.
;	    	z = [To(x)-.5, T1(x), T2(x), ..., Tn(x)]*coefs
;	Call matrix [To(x)-.5, T1(x), T2(x), ..., Tn(x)], matrix A 
;	with m rows and n colums.  m must be greater than n+1 in order 
;	for the system to be solvable.  Since m corresponds to the number 
;	of good pixels in the centered row, m may be greater than n,
;	indicating that the matrix may be overconstrained. We use 
;	the method of least squares to find the best fit coefs.
;               coefs = inverse(transpose(A)*A)*z
; 	        result1 = sum(Ti(xo)*coefs(i)) - .5*coefs(0),
;	where i = 0,1,...n
;	Result2 in the y direction is obtained in the similar fashion with
;	the change from x to y and row to column.
;       	result = (result1+result2)/2
;CALLING SEQUENCE:
;	Result = chebyshfit(raw,mask)
;INPUTS:raw	A 2-dimensional data array.
;	mask	A 2-dimensional integral array of the same size as 
;		raw data.  Allowed values are 1 and 0. 1 specifies
;		that the corresponding pixel in rawdata is good.
;		0 specifies that it is bad.
;KEYWORDS:maxn	maximum degree of Chebyshef polynomial.  If not set, maxn
;		is default to 10. 
;	debug	Set to 1 to display the coeficients of the interpolated
;		functions. 
;OUTPUT:str     A structure with the following fields:
;		image 	Value of the corrected bad pixel.
;		fixed	Set to 1 if the bad pixel is corrected, 0 otherwise.
;EXAMPLE:
;	raw = indgen(5,7)	;create image data array
;	raw(2,3) = -99
;	mask = intarr(5,7)+1	;create mask
;	mask(2,3) = 0
;	ch = chebyshfit(raw,mask)
;FUNCTIONS USED: ludc,lusol,chvarchg
;MODIFICATION HISTORY:
;	WRITTEN: Siree Vatanavigkit. June 26,1999
;	MODIFIED: Alex Ruane. November 20, 2000
;		bug fixes, chvarchg application	
;       MODIFIED: W. Landsman October 2001, Use TOTAL() instead of ROWSUM	
;-
if not(keyword_defined(debug)) then debug = 0
if not(keyword_defined(maxn)) then maxn = 10

;determine position of badpixel to be interpolated
info = size(raw)
rangex = [0,(info[1]-1)]		
rangey = [0,(info[2]-1)]		
xo = info[1]/2
yo = info[2]/2

;determine points for interpolation
row = total(mask,1)    ;Row sum
col = total(mask,2)    ;Column sum
hn   = row(yo)
vn   = col(xo)

;restrictions for fit
fit = 1
;Need at least four points to create a good polynomial(at least degree3)
if (min([hn,vn]) lt 4) then fit =0

if fit then begin
	;row
	if (hn gt maxn) then hn = maxn
	pixel = where(mask(*,yo) eq 1)*1d
	z = raw(pixel,yo)
	for k = 0, n_elements(pixel)-1 do begin
	  pixel[k] = chvarchg(pixel[k], rangex)
	endfor
	A = dblarr(hn,n_elements(pixel))
	for j = 0,n_elements(pixel)-1 do begin
		for i = 0,hn-1 do begin
	    		case i of
			0: A[i,j] = 1d
			1: A[i,j] = pixel(j)
		     else: A[i,j] = 2*pixel(j)*A(i-1,j)-A(i-2,j)
		    endcase
		endfor
		A[0,j] = 0.5d
	endfor
	AA = transpose(A)##A  ;AA is sym pos definite
	b  = transpose(A)##z
	b  = transpose(b)
	;solve for least square ie. find x that minimizes (A^TAx=A^Tz)
	ludc,AA,in
	coefs = lusol(AA,in,b)
	xcoefs = coefs		;acr addition
	if (debug) then print,"Xcoefs",coefs
	;calculate val1
	xof = chvarchg(xo, rangex)
	T=dblarr(hn)
	T[0] =	1d
	T[1] =  xof
	for i = 2,hn-1 do begin
	     T[i] = 2*xof*T(i-1)-T(i-2)
	endfor
	T[0] = 0.5d
	T = transpose(T)
	val1 = coefs##T

	;col
	if (vn gt maxn) then vn = maxn		;hn =maxn
	pixel = where(mask(xo,*) eq 1)*1d
	z = raw(xo,pixel)
	for k = 0, n_elements(pixel)-1 do begin
	  pixel[k] = chvarchg(pixel[k], rangey)
	endfor
	A = dblarr(vn,n_elements(pixel))
	for j = 0,n_elements(pixel)-1 do begin
		for i = 0,vn-1 do begin
	    		case i of
			0: A[i,j] = 1d
			1: A[i,j] = pixel(j)
		     else: A[i,j] = 2*pixel(j)*A(i-1,j)-A(i-2,j)
		    endcase
		endfor
		A[0,j] = 0.5d
	endfor
	AA = transpose(A)##A  ;AA is sym pos definite
	b  = transpose(A)##z
	b  = transpose(b)
	;solve for least square ie. find x that minimizes (A^TAx=A^Tz)
	ludc,AA,in
	coefs = lusol(AA,in,b)
	ycoefs = coefs
	If (debug) then print,"Ycoefs",coefs
	;calculate val2
	yof = chvarchg(yo, rangey)
	T=dblarr(vn)
	T[0] =	1.0d
	T[1] =  yof
	for i = 2,vn-1 do begin
	     T[i] = 2*yof*T(i-1)-T(i-2)
	endfor
	T[0] = 0.5d
	T = transpose(T)
	val2 = coefs##T
	str = {image:(val1+val2)/2.,fixed:1}
endif else str = {image:raw(xo,yo),fixed:0}
return, str
end