FUNCTION  ROBUST_SIGMA,Y, ZERO=REF, GOODVEC = Q
;
;+
; NAME:
;	ROBUST_SIGMA  
;
; PURPOSE:
;	Calculate a resistant estimate of the dispersion of a distribution.
; EXPLANATION:
;	For an uncontaminated distribution, this is identical to the standard
;	deviation.
;
; CALLING SEQUENCE:
;	result = ROBUST_SIGMA( Y, [ /ZERO, GOODVEC = ] )
;
; INPUT: 
;	Y = Vector of quantity for which the dispersion is to be calculated
;
; OPTIONAL INPUT KEYWORD:
;	/ZERO - if set, the dispersion is calculated w.r.t. 0.0 rather than the
;		central value of the vector. If Y is a vector of residuals, this
;		should be set.
;
; OPTIONAL OUPTUT KEYWORD:
;       GOODVEC = Vector of non-trimmed indices of the input vector
; OUTPUT:
;	ROBUST_SIGMA returns the dispersion. In case of failure, returns 
;	value of -1.0
;
; PROCEDURE:
;	Use the median absolute deviation as the initial estimate, then weight 
;	points using Tukey's Biweight. See, for example, "Understanding Robust
;	and Exploratory Data Analysis," by Hoaglin, Mosteller and Tukey, John
;	Wiley & Sons, 1983, or equation 9 in Beers et al. (1990, AJ, 100, 32)
;
; REVSION HISTORY: 
;	H. Freudenreich, STX, 8/90
;       Replace MED() call with MEDIAN(/EVEN)  W. Landsman   December 2001
;       Don't count NaN values  W.Landsman  June 2010
;
;-
 On_error,2
 compile_opt idl2
 
 EPS = 1.0E-20
 IF KEYWORD_SET(REF) THEN Y0=0. ELSE Y0  = MEDIAN(Y,/EVEN)

; First, the median absolute deviation MAD about the median:

 MAD = MEDIAN( ABS(Y-Y0), /EVEN )/0.6745

; If the MAD=0, try the MEAN absolute deviation:
 IF MAD LT EPS THEN MAD = MEAN( ABS(Y-Y0) )/.80
 IF MAD LT EPS THEN RETURN, 0.0

; Now the biweighted value:
 U   = (Y-Y0)/(6.*MAD)
 UU  = U*U
 Q   = WHERE(UU LE 1.0, COUNT)
 IF COUNT LT 3 THEN BEGIN
   PRINT,'ROBUST_SIGMA: This distribution is TOO WEIRD! Returning -1'
   SIGGMA = -1.
   RETURN,SIGGMA
 ENDIF

 N = TOTAL(FINITE(Y),/INT)      ;In case Y has NaN values          ;
 NUMERATOR = TOTAL( (Y[Q]-Y0)^2 * (1-UU[Q])^4 )
 DEN1  = TOTAL( (1.-UU[Q])*(1.-5.*UU[Q]) )
 SIGGMA = N*NUMERATOR/(DEN1*(DEN1-1.))
 
 IF SIGGMA GT 0. THEN RETURN, SQRT(SIGGMA) ELSE RETURN, 0.

 END