FUNCTION ROBUST_POLY_FIT,X,Y,NDEG,YFIT,SIG, NUMIT=THIS_MANY, DOUBLE=DOUBLE
;+
; NAME:
;	ROBUST_POLY_FIT
;
; PURPOSE:
;	An outlier-resistant polynomial fit.
;
; CALLING SEQUENCE:
;	COEFF = ROBUST_POLY_FIT(X,Y,NDEGREE, [ YFIT,SIG, /DOUBLE, NUMIT=] )
;
; INPUTS:
;	X = Independent variable vector, floating-point or double-precision
;	Y = Dependent variable vector
;       NDEGREE - integer giving degree of polynomial to fit, maximum = 6
; OUTPUTS:
;	Function result = coefficient vector, length NDEGREE+1. 
;	IF COEFF=0.0, NO FIT! If N_ELEMENTS(COEFF) > degree+1, the fit is poor
;	(in this case the last element of COEFF=0.)
;	Either floating point or double precision.
;
; OPTIONAL OUTPUT PARAMETERS:
;	YFIT = Vector of calculated y's
;	SIG  = the "standard deviation" of the residuals
;
; OPTIONAL INPUT KEYWORD:
;       /DOUBLE - If set, then force all computations to double precision.
;       NUMIT - Maximum number of iterations to perform, default = 25
; RESTRICTIONS:
;	Large values of NDEGREE should be avoided. This routine works best
;	when the number of points >> NDEGREE.
;
; PROCEDURE:
;	For the initial estimate, the data is sorted by X and broken into 
;	NDEGREE+2 sets. The X,Y medians of each set are fitted to a polynomial
;	 via POLY_FIT.   Bisquare ("Tukey's Biweight") weights are then 
;	calculated, using a limit  of 6 outlier-resistant standard deviations.
;	The fit is repeated iteratively until the robust standard deviation of 
;	the residuals changes by less than .03xSQRT(.5/(N-1)).
;
; PROCEDURES CALLED:
;        POLY(), POLY_FIT()
;       ROB_CHECKFIT()
; REVISION HISTORY
;	Written, H. Freudenreich, STX, 8/90. Revised 4/91.
;	2/94 -- changed convergence criterion
;        Added /DOUBLE keyword, remove POLYFITW call  W. Landsman  Jan 2009
;-

ON_ERROR,2
COMPILE_OPT IDL2

EPS   = 1.0E-20
DEL   = 5.0E-07
DEGMAX= 6

IF N_ELEMENTS(THIS_MANY) GT 0 THEN ITMAX=THIS_MANY ELSE ITMAX=25

BADFIT=0

NPTS = N_ELEMENTS(X)
MINPTS=NDEG+1
IF (NPTS/4*4) EQ NPTS THEN NEED2 = 1 ELSE NEED2 = 0
N3 = 3*NPTS/4  &  N1 = NPTS/4

; If convenient, move X and Y to their centers of gravity:
IF NDEG LT DEGMAX THEN BEGIN
   X0=TOTAL(X)/NPTS  &  Y0=TOTAL(Y)/NPTS
   U=X-X0            &  V=Y-Y0
ENDIF ELSE BEGIN
   U=X               &  V=Y
ENDELSE

; The initial estimate.

; Choose an odd number of segments:
NUM_SEG = NDEG+2
IF (NUM_SEG/2*2) EQ NUM_SEG THEN NUM_SEG =NUM_SEG+1
MIN_PTS = NUM_SEG*3
IF NPTS LT 10000 THEN BEGIN ;MIN_PTS THEN BEGIN
;  Settle for least-squares:
   LSQFIT = 1
   CC = POLY_FIT( U, V, NDEG, YFIT , DOUBLE=DOUBLE)
ENDIF ELSE BEGIN
;  Break up the data into segments:
   LSQFIT = 0
   Q = SORT(U)
   U = U[Q]  &  V = V[Q]
   N_PER_SEG = REPLICATE( NPTS/NUM_SEG, NUM_SEG)

;  Put the leftover points in the middle segment:
   N_LEFT = NPTS - N_PER_SEG[0]*NUM_SEG
   N_PER_SEG[NUM_SEG/2] = N_PER_SEG[NUM_SEG/2] + N_LEFT
   R = DBLARR(NUM_SEG)  &  S = DBLARR(NUM_SEG)
   R[0]=MEDIAN( U[0:N_PER_SEG[0]-1],/EVEN ) 
   S[0]=MEDIAN( V[0:N_PER_SEG[0]-1],/EVEN )
   I2 = N_PER_SEG[0]-1
   FOR I=1,NUM_SEG-1 DO BEGIN
     I1 = I2 + 1
     I2 = I1 + N_PER_SEG[I] - 1
     R[I] = MEDIAN( U[I1:I2], /EVEN)      &  S[I] = MEDIAN( V[I1:I2],/EVEN )
   ENDFOR
;  Now fit:
   CC = POLY_FIT( R,S, NDEG, DOUBLE=DOUBLE )
   YFIT = POLY(U,CC)  
ENDELSE

ISTAT = ROB_CHECKFIT(V,YFIT,EPS,DEL,  SIG,FRACDEV,NGOOD,W,S)

IF ISTAT EQ 0 THEN GOTO,AFTERFIT

IF NGOOD LT MINPTS THEN BEGIN
   IF LSQFIT EQ 0 THEN BEGIN
      ;  Try a least-squares:
      CC = POLY_FIT( U, V, NDEG, YFIT, DOUBLE=DOUBLE )
      ISTAT = ROB_CHECKFIT(V,YFIT,EPS,DEL,  SIG,FRACDEV,NGOOD,W,S)
      IF ISTAT EQ 0 THEN GOTO,AFTERFIT
      NGOOD = NPTS-COUNT
   ENDIF
   IF NGOOD LT MINPTS THEN BEGIN
      PRINT,'ROBUST_POLY_FIT: No Fit Possible!'
      RETURN,0.
   ENDIF
ENDIF

; Now iterate until the solution converges:
CLOSE_ENOUGH = .03*SQRT(.5/(NPTS-1)) > DEL 
DIFF= 1.0E10
SIG_1= (100.*SIG) < 1.0E20
NIT = 0
WHILE( (DIFF GT CLOSE_ENOUGH) AND (NIT LT ITMAX) ) DO BEGIN
  NIT=NIT+1
  SIG_2=SIG_1
  SIG_1=SIG
; We use the "obsolete" POLYFITW routine because it allows us to input weights
; rather than measure errors
  g = where(W gt 0, Ng)
  if Ng LT N_elements(w) then begin    ;Throw out points with zero weight
    u = u[g]
    v = v[g]
    w = w[g]
  endif  
  CC = POLY_FIT( U, V, NDEG, YFIT, MEASURE_ERRORS = 1/W^2, DOUBLE=DOUBLE )
  ISTAT = ROB_CHECKFIT(V,YFIT,EPS,DEL,  SIG,FRACDEV,NGOOD,W,S)
  IF ISTAT EQ 0 THEN GOTO,AFTERFIT
  IF NGOOD LT MINPTS THEN BEGIN
     PRINT,'ROBUST_POLY_FIT: Questionable Fit!'
     BADFIT=1
     GOTO,AFTERFIT
  ENDIF
  DIFF = (ABS(SIG_1-SIG)/SIG) < (ABS(SIG_2-SIG)/SIG)
ENDWHILE

;IF NIT GE ITMAX THEN PRINT,'ROBUST_POLY_FIT: Did not converge in',ITMAX,$
;' iterations!'

AFTERFIT:
CC=REFORM(CC)

IF NDEG LT DEGMAX THEN BEGIN
CASE NDEG OF
 1: CC[0] = CC[0]-CC[1]*X0 + Y0
 2: BEGIN   
   CC[0] = CC[0]-CC[1]*X0+CC[2]*X0^2 + Y0
   CC[1] = CC[1]-2.*CC[2]*X0
    END
 3: BEGIN
   CC[0] = CC[0]-CC[1]*X0+CC[2]*X0^2-CC[3]*X0^3 + Y0
   CC[1] = CC[1]-2.*CC[2]*X0+3.*CC[3]*X0^2
   CC[2] = CC[2]-3.*CC[3]*X0
    END
 4: BEGIN
   CC[0] = CC[0]-   CC[1]*X0+CC[2]*X0^2-CC[3]*X0^3+CC[4]*X0^4+ Y0
   CC[1] = CC[1]-2.*CC[2]*X0+3.*CC[3]*X0^2-4.*CC[4]*X0^3
   CC[2] = CC[2]-3.*CC[3]*X0+6.*CC[4]*X0^2
   CC[3] = CC[3]-4.*CC[4]*X0
    END
 5: BEGIN
   CC[0] = CC[0]-  CC[1]*X0+CC[2]*X0^2-CC[3]*X0^3+CC[4]*X0^4-CC[5]*X0^5+ Y0
   CC[1] = CC[1]-2.*CC[2]*X0+ 3.*CC[3]*X0^2- 4.*CC[4]*X0^3+5.*CC[5]*X0^4
   CC[2] = CC[2]-3.*CC[3]*X0+ 6.*CC[4]*X0^2-10.*CC[5]*X0^3
   CC[3] = CC[3]-4.*CC[4]*X0+10.*CC[5]*X0^2
   CC[4] = CC[4]-5.*CC[5]*X0
    END
 ENDCASE
ENDIF

; Calculate the fit at points X:
IF( N_PARAMS(0) GT 3 )THEN YFIT=POLY(X,CC)

IF BADFIT EQ 1 THEN CC=[CC,0.]

RETURN,CC
END