function gausfit, raw,mask,debug=debug
;NAME: 	gausfit
;PURPOSE:  
;	An interpolating function implementing a gaussian surface fit.
;	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.  
;EXPLANATIONS:
;	Fit a 2D gaussain function to the good data. By pluging 
;	in the coordinates of the centered pixel as variables to the 
;	fitted 2D gaussian equation, the interpolated value of the 
;	centered pixel is obtained.
;	   2D gaussian equation: 
;	z = (1/2*pi*thx)*exp(-.5*( ((x-mx)/thx)^2+((y-my)/thy)^2 ))- c
;	where mx,my are the means, thx,thy are the widths 
;       of the gaussian, and c is a positive number.
;       Rearranging the equation, we obtain the following.
;	 ln(z+c) = (-.5/thx^2)x^2 + (-mx/thx^2)x +
;		   (-.5/thy^2)y^2 + (-my/thy^2)y +
;		   (-.5((mx/thx)^2 +(my/thy)^2)+ln(1/2*pi*thx))
;	   Let x = a vector of the x coordinates of the good pixels,
;	       y = a vector of the y coordinates of the good pixels,
;	       z = a vector of the image data corresponding to (x,y), 
;	and coefs = a vector containing the coeficients [q r s t u] where
;	q = (-.5/thx^2) r = (-mx/thx^2) s = (-.5/thy^2) t = (-my/thy^2)
;       v =  (-.5((mx/thx)^2 +(my/thy)^2)+ln(1/2*pi*thx)).  		
;       We obtain the following linear system.
;	     ln(z+c) = [x^2 x y^2 y 1]*coefs
;	Call matrix [x^2 x y^2 y 1], matrix A with m rows and n colums.  
;	m must be greater than 5 in order for the system to be
;       solvable.  Since m corresponds to the number of good
;	pixels in the image, m may be greater then n, indicating that 
;       the matrix may have too many constraints. We use the method of 
;	least square to find the best fit coefs.
;               coefs = inverse(transpose(A)*A)*(ln(z+c))
;	   Let (xo,yo) be the coordinates of the centered pixel.
; 	        result = exp(qxo^2 + rxo + syo^2 + tyo + u) - c
;CALLING SEQUENCE:
;	Result = gausfit(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:
;	debug	Set to 1 to show the coeficients of the 2D gaussian equation.
;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
;	g = gausfit(raw,mask)
;FUNCTIONS USED:     coord.pro
;MODIFCATION HISTORY:
;	WRITTEN: Siree Vatanavigkit. June 26,1999
if not(keyword_defined(debug)) then debug = 0
;determine position of the centered bad pixel
info = size(raw)
xo = info[1]/2
yo = info[2]/2
;determine points for interpolation
points = coord(mask)
sp = size(points)
np = sp[2]	; num of fitting points

;restrictions of gaussian surface fit
if (np ge 5) then fit = 1 else fit = 0

if (fit) then begin
	z = raw(points(0,0:np-1),points(1,0:np-1))  ;values in z-direction
	;shift z to make the lowest point in ln(z) > zero
	if ( min(z) le 0 ) then c = abs(min(z))+1 else c = 0
        z = alog(z + c)
	A = dblarr(5,np)+1
	A[0,0:np-1] = (points[0,0:np-1])^2
	A[1,0:np-1] =  points[0,0:np-1]
	A[2,0:np-1] = (points[1,0:np-1])^2
	A[3,0:np-1] =  points[1,0:np-1]
	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)
	if(debug) then print, "gauscoefs",coefs
	val = coefs[0]*(xo^2)+coefs[1]*xo+coefs[2]*(yo^2)+coefs[3]*yo+coefs[4] 
        val = exp(val)-c
 	str = {image:val,fixed:1}
endif else str = {image:raw[xo,yo],fixed:0}
return, str
end