;;; Wrapper function for MPFIT
function gauss2d_wrap, p, x=x, y=y, err=err, data=data, loud=loud
fit = gauss2d(x, y, p)
return, reform((data-fit)/err, n_elements(data))
end

;;; Calculate the full width at half max of a 1-D Gaussian
function gaussfit2d_fwhm, x, fin
f = fin-min(fin)
mx = max(f)
w = where(f gt mx/2., nw)
dx = median(x - shift(x,1))
fwhm = float(nw)*dx
return, fwhm
end

;;; Main fitting program
function gaussfit2d, data, par, x, y
; How big is the image?
sz = size(data, /dim)
x = indgen(sz[0])
y = indgen(sz[1])

;;; Parameter Guesses
mx = max(data, imx)
mn = min(data, imn)
cr = array_indices(data, imx) ;Column and Row of peak
;;; The FWHM of a Gaussian is 2*sqrt(2*alog(2)) * sigma 
;;; or 2.355 * sigma
sigx = gaussfit2d_fwhm(x, data[*, cr[1]])/2.355
sigy = gaussfit2d_fwhm(y, data[cr[0], *])/2.355

pguess = [mn, mx-mn, sigx, sigy, cr[0], cr[1]]
;;; Function Arguments structure
fa = {x:x, y:y, err:sqrt(data), data:data}
;;; Levenberg-Marquardt Non-Linear Least-Squares Minimization
;;; I just love saying that!
par = mpfit('gauss2d_wrap', pguess, funct=fa, quiet=1-keyword_set(loud))
fit = gauss2d(x, y, par)
return, fit
end