FUNCTION ROBUST_POLY_FIT,X,Y,NDEG,YFIT,MEANABSDEV
;+
; Name:
; ROBUST_POLY_FIT
; Purpose:
; An outlier-resistant polynomial fit.
;
; Calling Sequence:
; COEFF = ROBUST_POLY_FIT(X,Y,NDEGREE  ,YFIT,MEANABSDEV)
; Inputs:
; X = Independent variable vector, floating-point or double-precision
; Y = Dependent variable vector
;
; Outputs:
; Function result = coefficient vector, length NDEGREE+1.
; Either floating point or double precision.
;
; Optional output parameters:
; YFIT = Vector of calculated y's
; MEANABSDEV = the mean absolute deviation of the y's from the fit
; curve
;
; Restrictions:
; Large values of NDEGREE should be avoided. This routine works best
; when the number of points >> NDEGREE.
;
; Procedure:
; IDL routine POLY_FIT is called for the initial estimate. Points with
; the largests mean absolute deviation are rejected, and the fit repeated.
; This goes on until the mean absolute fractional deviation stops decreasing,
; or only NDEGREE+1 points are left. Usually only a few iterations are 
; necessary.
;
;-

  ON_ERROR,2

  SIGMA_CUT = 3.0

  U = X
  V = Y
  NQ=N_ELEMENTS(U)
; The initial estimate:
  CC= POLY_FIT( U, V, NDEG, YFIT )
; Calculate the mean abs. deviation:
  DEV   = ABS(V-YFIT) 
  SIGMA = AVG(DEV) 
  LASTSIGMA=SIGMA
  LAST_NQ = NQ

; Loop until we converge, or run out of points:
  NUM_IT = 0
  WHILE( NQ GT NDEG ) DO BEGIN
     NUM_IT = NUM_IT + 1
;     PRINT,NQ,SIGMA
;    Reject points more than SIGMA_CUT mean absolute deviations:
     DEVLIM = SIGMA*SIGMA_CUT
     Q = WHERE(DEV LT DEVLIM)
     NQ=N_ELEMENTS(Q)
     IF( NQ EQ LAST_NQ )THEN Q=WHERE(DEV LT MAX(DEV))
     NQ=N_ELEMENTS(Q)
     IF( NQ GT NDEG )THEN BEGIN
      ; Re-fit:
        U=U(Q) & V=V(Q)
        CC= POLY_FIT( U, V, NDEG )
      ; If the mean abs. deviation didn't decrease, we're done:
        YFIT  = POLY(U,CC)
        DEV   = ABS(V-YFIT)
        SIGMA = AVG(DEV)
        IF( LASTSIGMA GE SIGMA )THEN GOTO,ALLDONE  ; EXIT!
        LASTSIGMA = SIGMA
     ENDIF ELSE BEGIN
;       PRINT,' ROBUST_POLY_FIT: Fit did not converge in',NUM_IT,' iterations.'
     ENDELSE
;    PRINT,NQ,SIGMA
  ENDWHILE
  ALLDONE:
  IF( N_PARAMS(0) GT 3 )THEN YFIT = POLY(X,CC)
  IF( N_PARAMS(0) GT 4 )THEN MEANABSDEV = AVG(ABS(Y-YFIT))
  RETURN,CC
END