FUNCTION HSI_PMTRAS_MULTIFIT,  x, y, error_est, yfit, okj
;
; Returns the parameters of a double-precision linear or 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-1)
;				okj			= 'nok' indices of acceptable data points, referenced to the input data array
; Tries both a linear and quadratic fit and uses whichever worked best.
;
; CALLS: 	hsi_pmtras_linfit
;			hsi_pmtras_quadfit
;
; 12-Sep-03 gh	Initial working version
;
nok 		= N_ELEMENTS(x)							; number of ok input data points (>2)
IF nok LE 1 THEN STOP, 'HSI_PMTRAS_MULTIFIT input vectors should have length >=2.'
linflag = 1											; preset flag to indicate that linear fit was best
linfit 	= hsi_pmtras_linfit(x,y, error_est, yfit, okj)
linfit	= REFORM(linfit, 2, /OVERWRITE)
linfit 	= [linfit, 0.D]								; Set quadratic parameter to zero for linear fit.
IF nok LE 4 THEN RETURN, linfit						; Require some redundancy to attempt quadratic fit.
quadfit	= hsi_pmtras_quadfit(x,y, quad_error_est, quad_yfit, quad_okj)
;
; Use linear fit if there is no valid quadratic error estimate or if it is larger than the linear value.
IF (quad_error_est EQ 0. OR quad_error_est GE error_est) THEN RETURN, linfit
;
; It seems that quadratic worked best - substitute quadratic versions of parameters and return
error_est 	= quad_error_est
yfit		= quad_yfit
okj			= quad_okj
RETURN, quadfit
END