FUNCTION CRVEFIT,X,Y,W,A,SIGMAA,TOLER,MINLIMIT,MAXLIMIT
;+
; NAME:
;	CURVEFIT
; PURPOSE:
;	Non-linear least squares fit to a function of an
;	arbitrary number of parameters.
;	Function may be any non-linear function where
;	the partial derivatives are known or can be approximated.
; CATEGORY:
;	E2 - Curve and Surface Fitting
; CALLING SEQUENCE:
;	YFIT = CRVEFIT(X,Y,W,A,SIGMAA)
; INPUTS:
;	X = Row vector of independent variables.
;	Y = Row vector of dependent variable, same length as x.
;	W = Row vector of weights, same length as x and y.
;		For no weighting
;		w(i) = 1., instrumental weighting w(i) =
;		1./y(i), etc.
;	A = Vector of nterms length containing the initial estimate
;		for each parameter.  If A is double precision, calculations
;		are performed in double precision, otherwise in single prec.
;	Sigmaa = sigmaa(i) ne 0 --> fit parameter i
;	         sigmaa(i) eq 0 --> use supplied parameter i, and fix it
;       Toler  = Vector of n terms containing the tolerance for each
;		 parameter.  If a parameter gets out of the tolerance
; 		 range, CRVEFIT will return -1 as its result.  If tolerance
;		 for a parameter is negative, CRVEFIT will exit, returning
;		 -1 if the corresponding parameter is negative or outside
;		 the tolerance range.
;	      ** The tolerance parameter is optional.
;	Minlimit, maxlimit = two vectors of n terms containing absolute 
;		  minimum and maximum values that the respective parameters
;		  can have.  If the parameter is out of this range, CRVEFIT
;		  will exit, returning -1.
;	      ** The minlimit and maxlimit parameters are optional, but both
;		  must be specified if they are used.
; 
; OUTPUTS:
;	A = Vector of parameters containing fit.
;	Function result = YFIT = Vector of calculated
;		values.  -1 if parameter outside of tolerance.
;	Sigmaa = Vector of standard deviations for parameters
;		A.
;	
; COMMON BLOCKS:
;	NONE.
; SIDE EFFECTS:
;	The function to be fit must be defined and called FUNCT.
;	For an example see FUNCT in the IDL User's Library.
;	Call to FUNCT is:
;	FUNCT,X,A,F,PDER
; where:
;	X = Vector of NPOINT independent variables, input.
;	A = Vector of NTERMS function parameters, input.
;	F = Vector of NPOINT values of function, y(i) = funct(x(i)), output.
;	PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
;		PDER(I,J) = Derivative of function at ith point with
;		respect to jth parameter.  Optional output parameter.
;		PDER should not be calculated if parameter is not
;		supplied in call (Unless you want to waste some time).
; RESTRICTIONS:
;	NONE.
; PROCEDURE:
;	Copied from "CURFIT", least squares fit to a non-linear
;	function, pages 237-239, Bevington, Data Reduction and Error
;	Analysis for the Physical Sciences.
;
;	"This method is the Gradient-expansion algorithm which
;	compines the best features of the gradient search with
;	the method of linearizing the fitting function."
;
;	Iterations are perform until the chi square changes by
;	only 0.1% or until 40 iterations have been performed.
;
;	The initial guess of the parameter values should be
;	as close to the actual values as possible or the solution
;	may not converge.
;
; MODIFICATION HISTORY:
;	Written, DMS, RSI, September, 1982.
;                CJB       July       1989 (input switches in sigmaa)
;-
	ON_ERROR,2		;RETURN TO CALLER IF ERROR
;        print,'one'
        if(total(abs(sigmaa)) eq 0.) then begin
           funct,x,a,yfit
           return,yfit
        endif
;        print,'two'

        mask = sigmaa
	AA = 1. * A	&    bb = aa	;MAKE PARAMS FLOATING
	initial = a
;        print,'three'

	NTERMS = N_ELEMENTS(where(mask ne 0.))	;# OF PARAMS.
	NFREE = (N_ELEMENTS(Y)<N_ELEMENTS(X))-NTERMS ;Degs of freedom
	IF NFREE LE 0 THEN STOP,'Crvefit - not enough data points.'
	FLAMBDA = 0.001		;Initial lambda
	DIAG = INDGEN(NTERMS)*(NTERMS+1) ;SUBSCRIPTS OF DIAGONAL ELEMENTS
;        print,'four'

;  check for optional tolerance parameter

	if ( n_params(0) eq 6 ) then begin
             checktol = 1 
             inittol = toler
        endif else begin
             checktol = 0
        endelse
;        print,'five'

;  check for optional minlimit and maxlimit parameters

	if ( n_params(0) eq 7 ) then begin
	  print,'ERROR : minlimit and maxlimit must BOTH be specified.'
	  return, -1
	endif
;        print,'six'

	if ( n_params(0) eq 8 ) then begin
             checklimits = 1 
             inittol = toler
             initmin = minlimit
             initmax = maxlimit
        endif else begin
             checklimits = 0
        endelse

        nparam=n_elements(a)
        a=fltarr(nterms)
        pder=fltarr(n_elements(x),nterms)
	parnum = intarr( nparam )
;        print,'seven'

	FOR ITER = 1,40 DO BEGIN	;Iteration loop

;		EVALUATE ALPHA AND BETA MATRICIES.
	FUNCT,X,AA,YFIT,PDERA	;COMPUTE FUNCTION AT A.
;        print,'total((y-yfit)^2)',total((y-yfit)^2)
        index=0
        for i=0,nparam-1 do begin
           if(mask(i) ne 0.) then begin
             pder(0,index)=pdera(*,i)
             a(index)=aa(i)
	     if ( n_params(0) ge 6 ) then begin	     
                toler(index) = inittol(i)
                minlimit(index) = initmin(i)
                maxlimit(index) = initmax(i)
             endif
	     parnum(index) = i
             index=index+1
           endif
        endfor
;        print,'eight'

	BETA = (Y-YFIT)*W # PDER
	ALPHA = TRANSPOSE(PDER) # (W # (FLTARR(NTERMS)+1)*PDER)

	CHISQ1 = TOTAL(W*(Y-YFIT)^2)/NFREE ;PRESENT CHI SQUARED.
;        print,' TOTAL(W*(Y-YFIT)^2)/NFREE,TOTAL(W),TOTAL((Y-YFIT)^2),nfree',$
;                TOTAL(W*(Y-YFIT)^2)/NFREE,TOTAL(W),TOTAL((Y-YFIT)^2),nfree


;	INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS.

	REPEAT BEGIN
;        print,'nine'

		C = SQRT(ALPHA(DIAG) # ALPHA(DIAG))
		ARRAY = ALPHA / C
		ARRAY(DIAG) = 1.+FLAMBDA
		ARRAY = INVERT(ARRAY)
		B = A + ARRAY / C # TRANSPOSE(BETA) ;NEW PARAMS

		if ( checktol ) then begin
		  for i = 0, nterms-1 do begin
		    if ( ( toler(i) lt 0. ) and ( b(i) lt 0. ) ) then begin
			print,'** Parameter ',parnum(i),' is negative'
			a = initial
		        toler = inittol
			minlimit = initmin
			maxlimit = initmax
			return, -1
		    endif
		    if ( abs( b(i) - initial(i) ) gt toler(i) ) then begin
			print,'** Tolerance exceeded for parameter ',parnum(i),$
			 ' Value = ',b(i)
			a = initial
			toler = inittol
			minlimit = initmin
			maxlimit = initmax
			return, -1
		    endif
		  endfor
		endif
;        print,'ten'

		if ( checklimits eq 1 ) then begin
		 for i = 0, nterms-1 do begin
		  if ((b(i) lt minlimit(i)) or (b(i) gt maxlimit(i))) then begin
			print,'** Limits exceeded for parameter ',parnum(i),$
			  ' Value = ',b(i)
			a = initial
			toler = inittol
			minlimit = initmin
			maxlimit = initmax
			return, -1
		   endif
		 endfor
		endif
;        print,'eleven'

                index=0
                for i=0,nparam-1 do begin
                  if ( mask(i) ne 0. ) then begin
                    bb(i) = b(index)
                    index = index + 1
                  endif
                endfor
		FUNCT,X,BB,YFIT		;EVALUATE FUNCTION
;        print,'twelve'
		CHISQR = TOTAL(W*(Y-YFIT)^2)/NFREE ;NEW CHISQR
		FLAMBDA = FLAMBDA*10.	;ASSUME FIT GOT WORSE
	ENDREP UNTIL CHISQR LE CHISQ1

	FLAMBDA = FLAMBDA/100.	;DECREASE FLAMBDA BY FACTOR OF 10
	Aa=Bb & a=b    ;SAVE NEW PARAMETER ESTIMATE.
	PRINT,'ITERATION =',ITER,' ,CHISQR =',CHISQR,'DIFF=',$
           (chisq1-chisqr)/chisq1
;	PRINT,A
	IF ((CHISQ1-CHISQR)/CHISQ1) LE .001 THEN GOTO,DONE ;Finished?
	ENDFOR			;ITERATION LOOP


	PRINT,'CRVEFIT - Failed to converge'

DONE:	
;        print,'thirteen'

	a=aa
        SIGMAAa = SQRT(ARRAY(DIAG)/ALPHA(DIAG)) ;RETURN SIGMA'S
                index=0
                for i=0,nparam-1 do begin
                  if(mask(i) ne 0.) then begin
                    sigmaa(i)=sigmaaa(index)
                    index=index+1
                  endif
                endfor
;                 print,'fourteen'

	RETURN,YFIT		;RETURN RESULT
END