FUNCTION HSI_PMTRAS_LINFIT,  x, y, error_est, yfit, okj
;
; Returns the parameters of a double-precision linear 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
; 	Requires at least 2 values of x.
;	If only   2 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.)
;
; 12-Nov-02 gh	Initial working version
; 15-Jul-03 gh	Add okj to list of returned parameters
;				Add provision for special case of 2-element input arrays.
; 22-Jul-03 gh	Minor correction to calculation of error_est
; 10-Sep-03 gh	Accomodate name change of hsi_outlier function.
;
nok 		= N_ELEMENTS(x)							; number of ok input data points (>2)
IF nok LE 1 THEN STOP, 'HSI_PMTRAS_LINFIT input vectors should have length >=2.'
ok			= INTARR(nok) + 1						; initialized to 1
okj			= WHERE (ok EQ 1)						; array of indices to input x,y arrays
;
; In special case of 2 available values, just determine the fitparm parameters directly.
;
IF nok EQ 2 THEN BEGIN
	fitparm 	= DINDGEN(2)
	fitparm[1] 	= (y[1]-y[0]) / (x[1]-x[0])
	fitparm[0] 	= y[0] - fitparm[1]*x[0]
	error_est	= 0.
	yfit		= y
	RETURN, fitparm
ENDIF
;
; Normal case of 3 or more available values...
;
done		= 0										; completion flag, set true after last iteration
REPEAT BEGIN
	fitparm = POLY_FIT(x[okj], y[okj], 1, /DOUBLE, YFIT=yfit)	; linear 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 nok EQ 3 OR outs[0] EQ -1 THEN done = 1		; done when only 3 ok points left or 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 * 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