FUNCTION HSI_PMTRAS_QUADFIT,  x, y, error_est, yfit, okj
;
; Returns the parameters of a double-precision quadratic fit of y vs x, with outlier rejection.
; 	Also sets 	yfit 		= fitted values of y,
;				est_error 	= (rms difference between yfit and y) / SQRT(nok-2)
;				okj			= 'nok' indices of acceptable data points, referenced to the input data array
; 	Requires at least 3 values of x.
;	If only   3 values are available, then error-est is set to 0.
;	If only   3 values are available, then no outlier rejection is performed.  (N_ELEMENTS(okj) will be >=2.)
;
; 10-Sep-03 gh	Initial version adapted from hsi_pmtras_linfit.
; 13-May-05 gh	Accomodate name change of hsi_outlier routine.
;
nok 		= N_ELEMENTS(x)							; number of ok input data points (>2)
IF nok LE 2 THEN STOP, 'HSI_PMTRAS_QUADFIT input vectors should have length >=3.'
ok			= INTARR(nok) + 1						; initialized to 1
okj			= WHERE (ok EQ 1)						; array of indices to input x,y arrays
;
; Begin iterative loop, rejecting one outlier at a time, if necessary.
;
done		= 0										; completion flag, set true after last iteration
REPEAT BEGIN
; In special case of just 3 ok values left, just determine the fitparm parameters directly and quit.
	IF nok EQ 3 THEN BEGIN
		fitparm 	= DBLARR(3)
		fitparm[2]	= ((y[okj[2]]-y[okj[1]]) / (x[okj[2]]-x[okj[1]])  - 	$
					   (y[okj[1]]-y[okj[0]]) / (x[okj[1]]-x[okj[0]])) / 	$
					   (x[okj[2]]-x[okj[0]])
		yp			= y - fitparm[2]*x^2
		fitparm[1] 	= (yp[okj[1]]-yp[okj[0]]) / (x[okj[1]]-x[okj[0]])
		fitparm[0] 	= yp[okj[0]] - fitparm[1]*x[okj[0]]
		error_est	= 0.
		yfit		= x^2 * fitparm[2] + x * fitparm[1] + fitparm[0]	; returned yfit array includes rejected points
		RETURN, fitparm
	ENDIF
	fitparm = POLY_FIT(x[okj], y[okj], 2, /DOUBLE, YFIT=yfit)	; quadratic fit using ok values only.
	diff	= y[okj] - yfit
	diff2	= diff^2
	outs	= hsi_outlier(diff)							; indices of outliers, or -1 if no outlier
	IF outs[0] EQ -1 THEN done = 1					; done when no outliers
	IF done EQ 0 THEN BEGIN
		dummy 			= MAX(diff2, jmax)			; jmax is referenced to truncated yfit array
		nok 			= nok -1
		ok[okj[jmax]]	= 0
		okj				= WHERE (ok EQ 1)
	ENDIF
ENDREP UNTIL (done EQ 1)
yfit 		= x^2 * fitparm[2] + x * fitparm[1] + fitparm[0]		; returned yfit array includes rejected points
error_est 	= SQRT(TOTAL(diff2)/(nok*(nok-1)))
;PRINT, MEAN(x), fitparm, error_est, N_ELEMENTS(x), nok , format='(4f10.5, 2I5)
RETURN, fitparm
END