function stgaussfit, x, y, w, a, sigmaa, Function_Name = Function_Name, $
                        itmax=itmax, iter=iter, tol=tol, chi2=chi2, $
                        noderivative=noderivative
; This is a modified version of the routines published with IDL
; to fit a Gaussian to a line.  Use of the IDL routines produced 
; incorrect fits to many lines.  For full documentation, see the
; help for the CURVEFIT routine in IDL.
       on_error,2              ;Return to caller if error

       ;Name of function to fit
       if n_elements(function_name) le 0 then function_name = "FUNCT"
       if n_elements(tol) eq 0 then tol = 1.e-3		;Convergence tolerance
       if n_elements(itmax) eq 0 then itmax = 20	;Maximum # iterations
	type = size(a)
	type = type(type(0)+1)
	double = type eq 5
	if (type ne 4) and (type ne 5) then a = float(a)  ;Make params floating

       ; If we will be estimating partial derivatives then compute machine
       ; precision
       if keyword_set(NODERIVATIVE) then begin
          res = nr_machar(DOUBLE=double)
          eps = sqrt(res.eps)
       endif

       nterms = n_elements(a)   ; # of parameters
       nfree = n_elements(y) - nterms ; Degrees of freedom
       if nfree le 0 then message, 'Curvefit - not enough data points.'
       flambda = 0.001          ;Initial lambda
       diag = lindgen(nterms)*(nterms+1) ; Subscripts of diagonal elements

;      Define the partial derivative array
       if double then pder = dblarr(n_elements(y), nterms) $
	else pder = fltarr(n_elements(y), nterms)
;
       for iter = 1, itmax do begin   ; Iteration loop

;         Evaluate alpha and beta matricies.
          if keyword_set(NODERIVATIVE) then begin
;            Evaluate function and estimate partial derivatives
             call_procedure, Function_name, x, a, yfit
             for term=0, nterms-1 do begin
                p = a       ; Copy current parameters
                ; Increment size for forward difference derivative
                inc = eps * abs(p(term))    
                if (inc eq 0.) then inc = eps
                p(term) = p(term) + inc
                call_procedure, function_name, x, p, yfit1
                pder(0,term) = (yfit1-yfit)/inc
             endfor
          endif else begin
             ; The user's procedure will return partial derivatives
             call_procedure, function_name, x, a, yfit, pder 
          endelse

          beta = (y-yfit)*w # pder
          alpha = transpose(pder) # (w # (fltarr(nterms)+1)*pder)
          chisq1 = total(w*(y-yfit)^2)/nfree ; Present chi squared.

				; If a good fit, no need to iterate
	  all_done = chisq1 lt total(abs(y))/1e7/NFREE
;
;         Invert modified curvature matrix to find new parameters.

          repeat begin
             c = sqrt(alpha(diag) # alpha(diag))
             array = alpha/c
             array(diag) = array(diag)*(1.+flambda)              
             array = invert(array)
             b = a+ array/c # transpose(beta) ; New params
             call_procedure, function_name, x, b, yfit  ; Evaluate function
             chisqr = total(w*(y-yfit)^2)/nfree         ; New chisqr
	     if finite(chisqr) ne 1 then begin
		yfit=fltarr(n_elements(y))
		a=[0.,0.,1.,0.]
		goto,crash
	     endif
	     if all_done then goto, done
             flambda = flambda*10.                      ; Assume fit got worse
          endrep until chisqr le chisq1
;
          flambda = flambda/100.  ; Decrease flambda by factor of 10
          a=b                     ; Save new parameter estimate.
          if ((chisq1-chisqr)/chisq1) le tol then goto,done  ; Finished?

       endfor                        ;iteration loop
;
       message, 'Failed to converge', /INFORMATIONAL
       yfit=fltarr(n_elements(y))	
       a=[0.,0.,1.,0.]
;
done:  sigmaa = sqrt(array(diag)/alpha(diag)) ; Return sigma's
       chi2 = chisqr                          ; Return chi-squared
crash: return,yfit              ;return result
END