Freudenreich AstroContrib Library

This page is a listing of the entire contents of this library for IDL. This listing is the long version. Viewing the much more compact listing may be handier.

[Go Back to Main IDL Libraries Search Page]


Last modified: Thu Dec 21 21:17:11 2000.

List of Routines


Routine Descriptions

AUTOHIST

[Next Routine] [List of Routines]
 NAME:
       AUTOHIST

 PURPOSE:
       Draws a histogram using automatic bin-sizing.
 
 EXPLANATION:
       AUTOHIST chooses a number of bins (initially, SQRT(2*N). If this leads 
       to a histogram in which > 1/5 of the central 50% of the bins are empty,
       it decreases the number of bins and tries again. The minimum # bins 
       is 5. The max=199.   Called by HISTOGAUSS and HALFAGAUSS.    Note that 
       the intrinsic HISTOGRAM function is *not* used.

 CALLING SEQUENCE:
       AUTOHIST, Sample, XLines, Ylines, XCenters, YCenters, [/NOPLOT]

 INPUT:
       Sample = the vector to be histogrammed

 OUTPUT:
       XLINES = the x coordinates of the points that trace the rectangular 
               histogram bins
       YLINES = the y coordinates. To draw the histogram plot YLINES vs XLINES.
       XCENTERS = the x values of the bin centers
       YCENTERS = the corresponding y values

 OPTIONAL INPUT KEYWORD:
       /NOPLOT  If set, nothing is drawn

 SUBROUTINE CALLS:
       MED, which calculates a median
 REVISION HISTORY:
       Written,   H. Freudenreich, STX, 1/91
       1998 March 17 - Changed shading of histogram.  RSH, RSTX

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/autohist.pro)


BBOOTSTRAP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BBOOTSTRAP

 PURPOSE:
	1-parameter bootstrap calculation of function FUNCT

CALLING SEQUENCE:
	biwt_means = bbootstrap( data, Num_Samp,[FUNCT='biweight_mean', /SLOW] )
		OR
	means = bbootstrap( data, Num_Samp )

 INPUT:
	DATA     = input vector 
	Num_Samp = the number of bootstrap samples to be taken. The more the 
		better.  Should be at least 30 if a standard deviation is to 
		be calculated.  At least 200 if 95% confidence limits are to 
		be calculated.
KEYWORD:
	FUNCT = the name of the function to be applied. If missing, AVERAGE is 
		assumed.
	SLOW    If present, the BALANCED bootstrap is performed. This is more
		accurate, especially for long-tailed distributions, but is also
		much slower and requires more memory.
 
 RETURNS:
         vector of NSAMP bootstrap answers. The mean, standard deviation
         and confidence intervals can be calculated from this

SUBROUTINES CALLED:
	'FUNCT'
	PERMUTE

NOTES:  
	This program randomly selects values to fill sets of the same size as Y,
	Then calculates the FUNCT of each. It does this NSAMP times. If the SLOW
	mode is requested, the balanced bootstrap method is used to obtain the 
	samples, hence the name BBOOTSTRAP. This method is preferable to that 
	used in BOOT_MEAN, but does require more virtual memory, and patience.

	The user should choose NSAMP large enough to get a good distribution.
	The sigma of that distribution is then the standard deviation of the 
	mean of ANSER.
	For example, if input X is normally distributed with a standard 
	deviation of 1.0, the standard deviation of the vector MEAN will be 
	~1.0/SQRT(N-1), where N is the number of values in X.

 WARNING:
	At least 5 points must be input. The more, the better.

 REVISION HISTORY:
	H.T. Freudenreich, HSTX, 2/95

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/bbootstrap.pro)


BIWEIGHT_MEAN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BIWEIGHT_MEAN 

 PURPOSE:
	Calculate the center and dispersion (like mean and sigma) of a 
	distribution using bisquare weighting.

 CALLING SEQUENCE:
	Mean = BIWEIGHT_MEAN( Vector, [ Sigma, Weights ] ) 

 INPUTS:
	Vector = Distribution in vector form

 OUTPUT:
	Mean - The location of the center.

 OPTIONAL OUTPUT ARGUMENTS:

	Sigma = An outlier-resistant measure of the dispersion about the 
	center, analogous to the standard deviation. The half-width of the 95%
	confidence interval = |STUDENT_T( .95, .7*(N-1) )*SIGMA/SQRT(N)|,
	where N = number of points.  

	Weights = The weights applied to the data in the last iteration.

SUBROUTINE CALLS:
	MED, which calculates a median

 REVISION HISTORY
	Written,  H. Freudenreich, STX, 12/89
	Modified 2/94, H.T.F.: use a biweighted standard deviation rather than
		median absolute deviation.
	Modified 2/94, H.T.F.: use the fractional change in SIGMA as the 
		convergence criterion rather than the change in center/SIGMA.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/biweight_mean.pro)


BOOT_BINDATA

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BOOT_BINDATA

 PURPOSE:
	Groups data into NBINS bins of equal size, and finds the robust mean of
	each bin. The X value per bin is not necessarily = bin-center. The boot-
	strap method is used. If a bin has fewer than 6 points, the biweighted
	mean is calculated; otherwise, a robust line is fitted to the data.

 CALLING SEQUENCE:
	BOOT_BINDATA, NBins, Xin, Yin, X1, X2, NSamp, X, Y, S, BinHits, [IType]

 INPUTS:
	NBins = number of bins
	Xin   = the input X values
	Yin   = the input Y values
	X1    = the low end of the range of X
	X2    = the high end of the range of X
	NSamp = the number of bootstrap samples
	IType (OPTIONAL) = use logarithmically equal bin-widths if present

 OUTPUTS:
	X       = the X values of the bins
	Y       = the y values (=-100000 if a bin is empty)
	S       = the robust standard deviations of the bins
	BinHits = the number of entries per bin.

 NOTES:

	If a bin has fewer than 6 entries, a biweighted mean of Y is
	calculated; X is weighted by the same weights, so the X value of the bin
	is not necessarily bin-center. If 6 points or more, a robust line is
	fitted to the data, and evaluated at bin-center. 

 REVISION HISTORY
	Written,  H.T. Freudenreich, Hughes STX, 12/3/92

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/boot_bindata.pro)


BOOT_MEAN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BOOT_MEAN

 PURPOSE:
	Calculate the Bootstrap Mean. (Also see BBOOTSTRAP.)

 CALLING SEQUENCE:
	BOOT_MEAN, y, N_Sample, Mean, [ /ROBUST ]
 INPUT:
	Y = vector to average
	N_SAMPLE = the number of bootstrap samples
 
 OUTPUT:
	MEAN = vector of N_SAMPLE means

 OPTIONAL INPUT KEYWORD:
	ROBUST - If this argument is present, an iterative biweighted mean is 
		returned

 SUBROUTINES CALLED:
	BIWEIGHT_MEAN, ROBUST_SIGMA, MED

NOTE:  
	This program randomly selects values to fill sets of the same size as Y,
	Then calculates the mean of each. It does this N_SAMPLE times.

	The user should choose N_SAMPLE large enough to get a good distribution.
	The sigma of that distribution is then the standard deviation of the 
	mean of MEAN.
	For example, if input Y is normally distributed with a standard 
	deviation of 1.0, the standard deviation of the vector MEAN will be 
	~1.0/SQRT(N-1), where N is the number of values in Y.

 WARNING:
	At least 5 points must be input. The more, the better.
 REVISION HISTORY:
	Written,   H.T. Freudenreich, HSTX, ?/92

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/boot_mean.pro)


BOOT_POLYFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BOOT_POLYFIT

 PURPOSE:
	Bootstrap polynomial fit to data. 

 CALLING SEQUENCE: 
	COEFF = BOOT_POLYFIT( X,Y, NDEG, N_SAMPLE, [ NUMOUT,  /ROBUST ] ) 

INPUT ARGUMENT:
	X = x vector
	Y = y vector
	NDEG = degree of polynomial fit
	N_SAMPLE = the number of bootstrap samples

 RETURNS:
	COEFF = array of coefficients. Dimensions=(NDEG+1,N_SAMPLE)
		First dimension in the order A0, A1, A2,...

 OUTPUT ARGUMENT:
	 NUMOUT = actual number of samples returned. If there is an error in
		the fit to any sample, NUMOUT is decreased by one.

 OPTIONAL INPUT KEYWORD:
	ROBUST  if present, an outlier-resistant fit is performed

 EXAMPLE:
	Takes 100 samples of a robust cubic fit to (X,Y)
	IDL> COEFF = BOOT_POLYFIT( X,Y, 3, 100, /ROBUST) 

 NOTE:  
	This program randomly selects (x,y) points and fits a curve to them. It
	does this N_SAMPLE times.

 WARNING:
	At least NDEG+1 points must be input. It is best to have at least twice
	that many for a robust fit.

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 3/16/94. Adapted from obsolete 
		BOOT_LINE.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/boot_polyfit.pro)


BOOT_XYFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	BOOT_XYFIT

 PURPOSE:
	Bootstrap linear fit to data in which there are errors in both Y and X,
	but the measurement errors are not available. Calculates the bisector 
	of the Y vs X and X vs Y regression.

 CALLING SEQUENCE:
	COEFF = BOOT_XYFIT( X,Y, NSAMPLE, YINT, SLOP, SYINT, SSLOP, NUMOUT, 
		[ /ROBUST ] ) 

 INPUT ARGUMENT:
	X = x vector
	Y = y vector
	N_SAMPLE = the number of bootstrap samples

 RETURNS:
	COEFF = array of coefficients. Dimensions=(NDEG+1,N_SAMPLE)
		First dimension in the order A0, A1

 OPTIONAL OUTPUT ARGUMENT:
	YINT = Y intercept (average of COEFF(0,*))
	SLOP = slope (average of (COEFF(1,*))
	SYINT = standard error of Y intercept (std. dev. of COEFF(0,*))
	SSLOP = standard error of slope (std. dev. of COEFF(1,*))
	NUMOUT = actual number of samples returned. If there is an error in
		the fit to any sample, NUMOUT is decreased by one.

 OPTIONAL INPUT KEYWORD:
	ROBUST  if present, an outlier-resistant fit is performed

 NOTES:  
	This program randomly selects (x,y) points and fits a line to them. The
	line is the bisector of the Y vs X and X vs Y regressions. It does this 
	N_SAMPLE times.

 WARNING:
	At least NDEG+1 points must be input. It is best to have at least twice
	that many for a robust fit.

 REVISION HISTORY:
	Written H.T. Freudenreich, HSTX, 3/16/94. Adapted from obsolete 
		BOOT_LINE.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/boot_xyfit.pro)


CHAINSAW

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	CHAINSAW

 PURPOSE:
	Remove isolated HIGH elements from a 2D array, replacing them with a 
	local surface fit. 

 CALLING SEQUENCE:
	NewMap = CHAINSAW( Map, Width,( FLOOR = ] )
 INPUT ARGUMENT:
	MAP = the array to de-point

INPUT ARGUMENT:
	WIDTH = the width of the square moving neighborhood centered on the 
		pixel in question. If there are too few non-empty pixels 
		within this neighborhood, nothing is done.

 OPTIONAL INPUT KEYWORS:
	FLOOR = the value representing an empty pixel. Default = 0.0

 OUTPUT:
	NEWMAP -  The input array, with point sources (or just exceptionally 
		bright pixels) removed.

 METHOD:
	Divides the neighborhood into 4 rectangular areas, finds the minimum 
	of each, and fits a plane to the 4 points. The central pixel is 
	replaced by the value of the plane at its location, unless it is 
	already fainter. The end result is much like that of EROSION+DILATION,
	but does not have mask-shaped artifacts. Edges of width WIDTH/2-1 are 
	unaffected.

 SUBROUTINES CALLED:
	Planefit, Med

 REVISIONS HISTORY:
	Written,  H.T. Freudenreich, HSTX, 10/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/chainsaw.pro)


CUBEROOT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	CUBEROOT
 PURPOSE:
	Real roots of a cubic equation. Complex roots set to -1.0e30.
	Called by HALFAGAUSS.
 CALLING SEQUENCE:
	roots = CUBEROOT(cc)
 INPUT:
	cc = 4 element vector, giving coefficients of cubic polynomial, 
		[c0,c1,c2,c3], where one seeks the roots of
		c0 + c1*x + c2*x^2 + c3*x^3
 OUTPUT:
	Function returns a vector of 3 roots. If only one root is real, then 
	that becomes the 1st element.
 EXAMPLE:
	Find the roots of the equation
		3.2 + 4.4*x -0.5*x^2 -x^3 = 0
	IDL> x = cuberoot([3.2,4.4,-0.5,-1])
	
	will return a 3 element vector with the real roots 
		-1.9228, 2.1846, -0.7618
 REVISION HISTORY:
	Henry Freudenreich, 1995

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/cuberoot.pro)


FITAGAUSS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	FITAGAUSS
 PURPOSE:
	A version of CURVEFIT that fits a Gaussian to the data, using the
	mean absolute deviation rather than the chi^2. DESIGNED ONLY FOR
	USE BY HISTOGAUSS AND HALFAGAUSS. EMPLOY ELSEWHERE AT USER'S RISK.
 CALLING SEQUNENCE: 
	value = FITAGAUSS( X, Y, A, SIGMAA )
 INPUT:
	X, Y = the points to fit
	A = the vector of parameters: [height,center,width,Y offset]. The
	offset is optional. If the A input has 3 element, no offset is returned.
       If A has four elements, the last parameter is a constant offset:
       y = G(x) + A(3)
 OUTPUT:
	Value -	The value of the Gaussian fit at the input points.
	SIGMAA = uncertainties of parameters
 REVISION HISTORY:
	Written,   ; H.T. Freudenreich, ?/?

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/fitagauss.pro)


FITXYERRS[1]

[Previous Routine] [Next Routine] [List of Routines]
 Part of FITXYERRS.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/fitxyerrs.pro)


FITXYERRS[2]

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	FITXYERRS

 PURPOSE:
	Linear fit to data in which there are errors in both Y and X, and 
	possibly outliers or points with true errors much larger than the 
	measurement errors.   It is outlier-resistant but not as resistant as 
	ROBUST_LINEFIT.   Uncertainties are calculated using the bootstrap 
	method.   

 CALLING SEQUENCE:
	Coeff = FITXYERRS( X, Y, Ex, Ey, N_sample, [ Yint, Slope, Sig_Yint, 
				Sig_Slope, /NONROBUST, Tolerance = ] )

 INPUT ARGUMENT:
	X = Input x vector
	Y = Input y vector
	Ex= vector of X measurement errors
	Ey= the vector of Y measurement errors
	N_SAMPLE = the number of bootstrap samples. 50 to 100 usually sufficient
			If N_SAMPLE = 1, standard errors are not returned.

 OPTIONAL INPUT KEYWORDS:
	NONROBUST -   If set, "robustness" weights are not calculated.
	TOLERANCE -   When the scatter about the fitted line, weighted by
		measurement errors, changes by less than this fraction,
		the fit is considered to have converged. Default=.0001

 OUTPUT:
	COEFF = array of coefficients. Dimensions=(NDEG+1,N_SAMPLE)
		First dimension in the order A0, A1

 OPTIONAL OUTPUT ARGUMENT:
	Yint = Y intercept (average of COEFF(0,*))
	Slop = slope (average of (COEFF(1,*))
	Sig_Yint = standard error of Y intercept (std. dev. of COEFF(0,*))
	Sig_Slop = standard error of slope (std. dev. of COEFF(1,*))

 METHOD:  
	Convert to standardized variables, calculate Y vs X and X vs Y fits,
	take a weighted mean of the two lines, the weights being a function of
	slope. The weights of the input points are functions of their 
	uncertainties and perpendicular distance to the fitted line. This is 
	an iterative procedure.   To the author's knowledge, neither this nor 
	a similar method has ever been, or may ever be, published.

 WARNING:
	At least 6 points must be input. It is better to have at least 10.
	It is best to have a large number.
	For large datasets, try it out with N_SAMPLE of 1 first to see how long 
	it takes.

 SUBROUTINES CALLED:
	XYFIT, MED (to calculate a median)

 REVSION HISTORY:
	Written,  H.T. Freudenreich, HSTX, 6/94. 
	Corrected serious bug in computation of the Y intercept 
		H.T. Freudenreich/Landsman HSTX, 11/95
	

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/fitxyerrs.pro)


GAUSSFUNC

[Previous Routine] [Next Routine] [List of Routines]
 Part of FITAGAUSS
 Evaluates a Gausian and its partial derivative for FITAGAUSS

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/gaussfunc.pro)


HALFAGAUSS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	HALFAGAUSS

 PURPOSE:
	Histogramming a distribution that can be characterized by a Gaussian
	with one heavy tail and returns the parameters of the Gaussian. Draws
	the histogram, the Gaussian, and the smoothed residuals. If a plausible
	fit to a Gaussian cannot be made (typically, if the distribution has a
	sharp edge or insufficient data near the mode) the parameters of the
	Gaussian are estimated. In this case an asterisk precedes the label of
	mean and width drawn on the plot.

 CALLLING SEQUENCE: 
	YFit = HALFAGAUSS( Data, Parm, X, Y, [ /NOPLOT, STITLE ] )

 INPUT:
	DATA = the data to histogram. Should be large. At least 500 points;
		10,000 is better.
 OUTPUT:
	YFit - HALFAGAUSS returns the vector of Gaussian fit to the histogram. 
		If the program fails, it will return a scalar, 0.0.

 OPTIONAL OUTPUT:
	PARM = the parameters of the Gaussian component of the distribution:
		[height,mode,sigma]
	X = the locations of bin-centers
	Y = the locations of the bin-means

 OPTIONAL INPUT KEYWORDS: 
	NOPLOT - if set, a histogram is not drawn
	STITLE - subtitle, if desired

SUBROUTINES NEEDED: 
	HISTOGAUSS,AUTOHIST,ROBUST_SIGMA,ROBUST_POLY_FIT,ROB_CHECKFIT,MED,
	LOWESS, POLY4PEAK,FITAGAUSS and GAUSSFUNC.

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 5/6/94
	2/10/95: Changed location of labels. HF

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/halfagauss.pro)


HISTOGAUSS

[Previous Routine] [Next Routine] [List of Routines]
NAME:
	HISTOGAUSS

 PURPOSE:
	Histograms data and overlays it with a Gaussian. Draws the mean, sigma,
	and number of points on the plot.

 CALLING SEQUENCE:
	HISTOGAUSS, Sample, A, XX, YY, GX, GY, [/NOPLOT, /NOFIT, CHARSIZE = ]

 INPUT:
	SAMPLE = Vector to be histogrammed

 OUTPUT ARGUMENTS:
	A = coefficients of the Gaussian fit: Height, mean, sigma
		A(0)= the height of the Gaussian
		A(1)= the mean
		A(2)= the standard deviation
		A(3)= the half-width of the 95% conf. interval 
		A(4)= 1/(N-1)*total( (y-mean)/sigma)^2 ) = a measure of 
			normality

	Below: superceded. The formula is not entirely reliable.
	A(4)= measure of the normality of the distribution. =1.0, perfectly
       normal. If no more than a few hundred points are input, there are
       formulae for the 90 and 95% confidence intervals of this quantity:
       M=ALOG10(N-1) ; N = number of points
       T90=ABS(.6376-1.1535*M+.1266*M^2)  ; = 90% confidence interval
       IF N LT 50 THEN T95=ABS(-1.9065-2.5465*M+.5652*M^2) $
                  ELSE T95=ABS( 0.7824-1.1021*M+.1021*M^2)   ;95% conf.
       (From Martinez, J. and Iglewicz, I., 1981, Biometrika, 68, 331-333.)

	XX = the x coordinates of the histogram bins (CENTER)
	YY = the y coordinates of the histogram bins

 OPTIONAL INPUT KEYWORDS:
	NOPLOT - If set, nothing is drawn
	FITIT   If set, a Gaussian is actually fitted to the distribution.
		By default, a Gaussian with the same mean and sigma is drawn; 
		the height is the only free parameter.
	CHARSIZE Size of the characters in the annotation. Default = 0.82.

 SUBROUTINE CALLS:
	BIWEIGHT_MEAN, which determines the mean and std. dev.
	AUTOHIST, which draws the histogram
	FITAGAUSS, which does just that
	MED, which calculates a median (called by AUTOHIST)

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 12/89
	Modified for IDL Version 2, 1/91, HF
	More quantities returned in A, 2/94, HF
	Added NOPLOT keyword and print if Gaussian, 3/94
	Stopped printing confidence limits on normality 3/31/94 HF
	Added CHARSIZE keyword, changed annotation format, 8/94 HF
	Simplified calculation of Gaussian height, 5/95 HF

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/histogauss.pro)


LOESS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	LOESS

 PURPOSE:
	Map-smoothing using the loess method (a local, weighted, polynomial fit
	in two dimensions). Allows fits only up to degree 2 in X and Y.

 CALLING SEUENCE:
 	smooth_map = LOESS( map, width, ndegree, [Noise, FLOOR = ,MINPTS =,
			/LOG, EDGE= ] )

 INPUT ARGUMENT:
	MAP = the array to smooth
	WIDTH = the (square) width of the local neighborhood over which 
		smoothing is performed.
	NDEGREE = the degree of the polynomial (in X and Y) used. 1 or 2 
		allowed.

 OUTPUT:
	SMOOTH_MAP - LOESS returns the smoothed map

 OPTIONAL OUTPUT:
	NOISE = the map of std. dev. (scaled from the median abs. deviation) of 
	the surface-fit residuals in the neighborhood centered on each pixel. 

 OPTIONAL INPUT KEYWORDS:
	FLOOR = the minimum value allowed. If a pixel < FLOOR it is considered
		 vacant.  The default FLOOR = 0.0
	MINPTS= the minimum # occupied pixels per neighborhood. 
		Default = total/2.  If fewer pixels are occupied, no fit is 
		made and the value of the central pixel is unchanged. MINPTS 
		must be > 3 for 1st degree fit and > 6 for 2nd degree fit.
	LOG  = if set, the LOG of the map is used in the local fits, then 
		exponentiated before returning (and before calculating noise)
	EDGE = if set, the edges of the map are not smoothed.

 SUBROUTINES CALLED:
	ROBUST_PLANEFIT, ROB_QUARTICFIT, ROB_MAPFIT, ROBUST_SIGMA, MED. 

 NOTE ON USAGE:
	If too many pixels within a a neighborhood are empty, the central pixel 
	retains its old value and noise, if desired, set to -1. for that pixel.

 METHOD:
	In each neighborhood: a polynomial (plane or quartic) is fitted robustly
	(resistant to outliers). Weights are taken as functions of the deviation
	from this fit (biweights) * a tri-cubic function of the distance of the
	neighborhood pixels from the center of the neighborhood. Then another
	surface is fitted using these weights. The edges are obtained from 
	robust surface fits to a neighborhood of width (WIDTH-2)^2; no 
	distance weighting is used on the edges.

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 5/27/93
	Keywords added by H.F., 3/94
	Added test on quadrant coverage, HF, 9/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/loess.pro)


LOWESS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	LOWESS

 PURPOSE:
	Robust smoothing of 1D data. A non-parametric way of drawing a smooth
	curve to represent data. (For 2D (map) data, use LOESS.)

 CALLING SEQUENCE:
	YSmooth = LOWESS( X, Y, Window, NDeg, Noise )

 INPUT ARGUMENTS:
	X = X values
	Y = Y values, to be smoothed
	WINDOW = width of smoothing window
	NDEG   = degree of polynomial to be used within the window (1 or 2 
		recommended)

 OUTPUT:
	Ysmooth - LOWESS returns the vector of smoothed Y values

 OPTIONAL OUTPUT ARGUMENT:
	Noise = the robust std. deviation w.r.t. the fit, at each point

 NOTE:
	This routine uses a least-squares fit within a moving window. The fit
	is weighted by statistical weights and weights that are a function of
	distance from the center of the window. This is a "local weighted
	polynomial regression smoother" (Cleveland 1979, Journal of the Amer.
	Statistical Association, 74, 829-836).
	A polynomial of degree NDEG+1 is fitted directly to the first and last 
	WINDOW/2 points. 

	This routine is fairly slow. 

 SUBROUTINES CALLED:
	ROBUST_LINEFIT
	ROBUST_POLY_FIT
	ROB_CHECKFIT
	ROBUST_SIGMA
	POLYFITW
	MED

 REVISION HISTORY:
	Written,  H.T. Freudenreich, HSTX, 1/8/93
	H.T. Freudenreich, 2/94 Return sigma rather than slope

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/lowess.pro)


MAPFITW

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	MAPFITW

 PURPOSE:
	Fit a 2D surface to a map through call to REGRESS. Each pixel in the map
	may be given a weight. (If pixels have equal weight, SURFACE_FIT 
	suffices.)

 CALLING SEQUENCE:
	surface = MAPFITW( map,ndeg,w )

 INPUT:
	map = array to fit, floating-point
 ndeg= degree in X and Y of polynomials to be fitted
       Maximum degree = 7 (35 terms)
 W   = array of same size and type as map containing the pixel weights.

RETURNS:
 The calculated fitted surface at each pixel: an array like the input array.

KEYWORDS:
 SINGLE  if set, single precision is used. This should be done only if there
         is not enough memory for the default double-precision calculations.

SUBROUTINES NEEDED:
 REGRESS (in USERLIB)

NOTE:
 This routine merely convert the pixel coordinates into a form REGRESS can 
 digest.

AUTHOR: H.T. Freudenreich, HSTX, 3/16/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/mapfitw.pro)


MED

[Previous Routine] [Next Routine] [List of Routines]
 NAME
	 MED
 PURPOSE
	Compute the median of an array, which may be of even length. 

 CALLING SEQUENCE:  
 	MID_VALUE = MED(A)

 OUTPUTS 
 	The median of array A
 REVISION HISTORY:
	H.T. Freudenreich, ?/89

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/med.pro)


MEDSMOOTH

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	MEDSMOOTH

 PURPOSE:
	Median smoothing of a vector, including pointe near its ends.

 CALLING SEQUENCE:
	SMOOTHED = MEDSMOOTH( VECTOR, WINDOW_WIDTH )

 INPUTS:
	VECTOR  = The vector to be smoothed
	WINDOW = The full width of the window over which the median is 
		determined for each point.

 OUTPUT:
	Function returns the smoothed vector

 SUBROUTINES CALLED: 
	MEDIAN, to find the median

 PROCEDURE:
	Each point is replaced by the median of the nearest WINDOW of points.
	The width of the window shrinks towards the ends of the vector, so that
	only the first and last points are not filtered. These points are 
	replaced by forecasting from smoothed interior points.

 REVISION HISTORY:
 	Written, H. Freudenreich, STX, 12/89
	H.Freudenreich, 8/90: took care of end-points by shrinking window.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/medsmooth.pro)


PERMUTE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	PERMUTE
 PURPOSE:
	Randomly permute the elements of a vector
 USAGE:
	NewIndex = PERMUTE( M, [ Seed] )
 INPUT:
	M = length of vector
 OPTIONAL
	SEED = random number seed to be passed to RANDOMU (optional)
 OUTPUT:
	PERMUTE returns M randomly shuffled integers between 0 and M-1
 EXAMPLE:
	To shuffle the elements of a vector V,

	V = V( PERMUTE(N_ELEMENTS(V))	

 WARNING: 
	This can be pretty slow. 

 REVISION HISTORY:
	Written,  H.T. Freudenreich, HSTX, 2/95

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/permute.pro)


PERMUTEST

[Previous Routine] [Next Routine] [List of Routines]
 NAME:  
	PERMUTEST
 PURPOSE: 
	Apply Fisher's Permutation Test for the equality of the means of two
	samples. This method is non-parametric and exact--and slow.

 CALLING SEQUENCE:
	asl = PERMUTEST( sample1, sample2, [ TD] )
 INPUT:	
	SAMPLE1, SAMPLE2 = vectors containing samples A and B

 OUTPUT:
	ASL - PERMUTEST returns the Achieved Significance Level, or the 
	probability the two distributions are the same. A fraction 0. to 1.0.

 OPTIONAL OUTPUT: 
	TD = the distribution of the difference of the means
 NOTE: 
	This is a SLOW routine. It may not be appropriate for large (~N > 1000) 
	samples, depending on the user's patience and the speed of his machine.
 REVISION HISTORY:
	H.T. Freudenreich, HSTX, 2/1/95

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/permutest.pro)


PLANEFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	PLANEFIT
 PURPOSE:
	Least-squares fit of a plane to a set of (X,Y,Z) points: 
		Z=R(0)+R(1)*X+R(2)*Y
 CALLING SEQUENCE:
	coef = PLANEFIT( x,y,z,0., yfit )  ; for equal weighting
	coef = PLANEFIT( x,y,z, weights_vector, yfit ) ;otherwise
 INPUT:
	XIN, YIN, ZIN = the points to be fit
	Weights_vector = their weights. (For equal weights, let W=0.)
 RETURNS:
	coef = the fit coefficients
 OPTIONAL OUTPUT:
	YFIT = the Z values at (X,Y), calculated from the fit

 NOTE:
	If a fit is not possible, R will be returned as a scalar, the average 
	of Z.
 REVISION HISTORY:
	Written by H.T. Freudenreich, HSTX, 5/18/93

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/planefit.pro)


POINT_REMOVER

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	POINT_REMOVER

 PURPOSE:
	Remove isolated HIGH SIGNAL/NOISE elements from a 2D array, replacing
	them with a local median. This routine iterates until convergence or
	a user-supplied maximum number. This is NOT a low-pass filter. Only
	points passing the S/N cut are affected. 

 CALLING SEQUENCE:
	NewMap = POINT_REMOVER( Map, Width, Sigma_Mult, ItMax,[ /BOTH, FLOOR= ])

 INPUT ARGUMENT:
	MAP = the array to de-point
	WIDTH = the width of the square moving neighborhood centered on the 
		pixel in question. If half the pixels within this neighborhood
		 are vacant, nothing is done.

 OPTIONAL INPUT ARGUMENT:
	SIGMA_MULT = the number of standard deviations by which a pixel must 
		exceed the local background for it to be replaced. DEFAULT = 2.0
	ITMAX = the maximum number of iterations. DEFAULT = 20

 OPTIONAL INPUT KEYWORDS:
	BOTH - if set, the absolute value of the difference between a pixel and 
		its neighborhood is used, so that low pixels may also be 
		replaced. Else, only + high points will be replaced.
	FLOOR = the value representing an empty pixel. Default = 0.0

 OUTPUT:
 	POINT_REMOVER returns the array, with point sources removed.

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, ?/90 or 91
	HF, 3/94, to include FLOOR and USE_ABS keywords

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/point_remover.pro)


POLY4PEAK

[Previous Routine] [Next Routine] [List of Routines]
 Finds the maximum of a 4th degree polynomial c0+c1x+c2x^2+c3x^3+c4x^4
 Called by HALFAGAUSS
 INPUT: COEFFICIENTS OF POLYNOMIAL. RETURNS: THE LOCATION OF THE MAX.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/poly4peak.pro)


QUARTICFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
      QUARTICFIT

 PURPOSE:
      Fit a general quadratic surface: Z=a+bX+cX^2+dXY+eY+fY^2

 CALLING SEQUENCE:
      coefficients = QUARTICFIT( X,Y,Z,W, YFIT )
 INPUT:
      X = independent variable vector, may be in double precision
      Y = independent variable vector, may be in double precision
      Z = dependent variable vector, may be in double precision
      W = weights (IF EQUAL, SET W=0.)

 RETURNS:
      The coefficients of the fit, in the order given above.

 OPTIONAL OUTPUT:
      YFIT = the calculated fit at points (X,Y)

 SIDE-EFFECTS:
      Sometimes nausea and dizziness, if taken in large doses.

 AUTHOR:
      H.T. Freudenreich, HSTX, 5/19/93

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/quarticfit.pro)


RESISTANT_MEAN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	RESISTANT_MEAN  

 PURPOSE:
	An outlier-resistant determination of the mean and its standard 
	deviation.  It trims away outliers using the median and the median 
	absolute deviation.

 CALLING SEQUENCE:
	RESISTANT_MEAN,VECTOR,SIGMA_CUT, MEAN,SIGMA,NUM_REJECTED

 INPUT ARGUMENT:
	VECTOR    = Vector to average
	SIGMA_CUT = Data more than this number of standard deviations from the
		median is ignored. Suggested values: 2.0 and up.

 OUTPUT ARGUMENT:
	MEAN  = the mean
	SIGMA = the standard deviation of the mean
	NUM_REJECTED = the number of points trimmed

 SUBROUTINE CALLS:
	MED, which calculates a median

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 1989; Second iteration added 5/91.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/resistant_mean.pro)


ROBUST_BIN2D

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_BIN2D

 PURPOSE:
	Bins data in a 2-dimensional array suitable for surface or contour 
	plots.  Cells with more than 2 entries are robustly averaged to remove
	outliers.

 CALLING SEQUENCE:
	ROBUST_BIN2D,NX,NY,X0,X1,Y0,Y1,X,Y,DATA,Z,NUM,SIGGMA

 INPUTS:
	NX = Number of bins in the X direction
	NY = Number of bins in the Y direction
	X0 = Lower bound on X
	X1 = Higher bound on X
	Y0 = Lower bound on Y
	Y1 = Higher bound on Y
	X = Independent variable vector, floating-point
	Y = Dependent variable vector
	DATA = Vector of quantity to be binned. X, Y, DATA must have same 
		length.

 OUTPUTS:
	Z   = 2-dimensional array containing the average Z per bin
	NUM = 2-dimensional array containing the number of entries per bin.
	SIGGMA = (OPTIONAL) standard deviation per pixel (not OF THE MEAN).

 SUBROUTINES CALLED:
	RESISTANT_MEAN, MED

 REVISION HISTORY:
	Written by H.T. Freudenreich, ?/1990
	Modified: H.F., 3/94 to return SIGGMA

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_bin2d.pro)


ROBUST_BINDATA

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_BINDATA

 PURPOSE: 
	Groups data into NBINS bins of equal size, and finds the robust mean of
	each bin. If there are enough pts in a bin, a line is fitted to the 
	data.   The bootstrap method is used to determine the uncertainty of 
	the bin averages.

 CALLING SEQUENCE:
	ROBUST_BINDATA, NBINS, XIN, YIN, X1, X2, X, Y, S, BINHITS,itype
 INPUT:
	NBINS = number of bins
	XIN   = the input X values
	YIN   = the input Y values
	X1    = the low end of the range of X
	X2    = the high end of the range of X
	ITYPE (OPTIONAL) = use logarithmically equal bin-widths if present

 OUTPUTS:
	X       = the X values of the bins
	Y       = the y values (=-100000 if a bin is empty)
	S       = the robust standard deviations of the bins
	BINHITS = the number of entries per bin.

 NOTES:   

	If a bin has fewer than 5 entries, a biweighted mean of Y is
	calculated; X is weighted by the same weights, so the X value of the bin
	is not necessarily bin-center. If 5 points or more, a robust line is
	fitted to the data, and evaluated at bin-center.

 REVISION HISTORY:
	Written, H.T. Freudenreich, Hughes STX, 12/3/92

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_bindata.pro)


ROBUST_BOXCAR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_BOXCAR 

 PURPOSE:
	Calculate robust boxcar averages of fixed number of points. One average
	is calculated per NPER points. (Remaining points go into the last 
	average.)   Biweights are used. The same weights are applied to the X 
	variable in calculating its average, so that the (X,Y) correspondence 
	is not lost.

 CALLING SEQUENCE:
	ROBUST_BOXCAR, X, Y, NPER,  U, V, SIG, [ CONF_INT = ]

 INPUT ARGUMENTS:
	X = X coordinates vector
	Y = Distribution in vector form
	NPER = # points per average

 OUTPUT ARGUMENT:
	U = vector of average X
	V = vector of average Y

 OPTIONAL OUTPUT ARGUMENT:
	SIG = standard error of V (a robust std. dev. of the mean) [OPTIONAL]

 OPTINAL INPUT KEYWORD:
	CONF_INT = confidence interval in percent. If this keyword is present
		then output SIG = the confidence interval

 SUBROUTINE CALLS:
	BIWEIGHT_MEAN, which calculates the averages

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 1/92
	Modifications: H.F. 3/94 -- added confidence intervals option

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_boxcar.pro)


ROBUST_CORR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_CORR

 PURPOSE:
	Derive an outlier-resistant measure of the correlation coefficient of
	variables X and Y.

 CALLING SEQUENCE:
	Correl_coeff = ROBUST_CORR( X, Y )

 INPUT ARGUMENTS:
	X = Vector of quantity X
	Y = Vector of quantity Y

 RETURNS:
	Estimate of correlation coefficient. In the absence of outliers this
	equals the true correlation coefficient.

 CALLS:
	Function ROBUST_LINEFIT to perfom an outlier-resistant fit. 
	Function ROBUST_SIGMA to calculate a resistant analog to the
	standard deviation. 
	Also: ROB_CHECKFIT, MED

 REVISION HISTORY:
	Written,   H. Freudenreich, STX, 8/90

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_corr.pro)


ROBUST_LINEFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_LINEFIT

 PURPOSE:
	An outlier-resistant two-variable linear regression. Either Y on X
	or, for the case in which there is no true independent variable, the
	bisecting line of Y vs X and X vs Y is calculated. No knowledge of
	the errors of the input points is assumed.

 CALLING SEQUENCE:
	COEFF = ROBUST_LINEFIT( X, Y, YFIT, SIG, COEF_SIG, [ BISECT = ,
			BiSquare_Limit = , Close_factor = , NumIT = ] )

 INPUTS:
	X = Independent variable vector, floating-point or double-precision
	Y = Dependent variable vector

 OUTPUTS:
	Function result = coefficient vector. 
	If = 0.0 (scalar), no fit was possible.
	If vector has more than 2 elements (the last=0) then the fit is dubious.

 OPTIONAL OUTPUT PARAMETERS:
	YFIT = Vector of calculated y's
	SIG  = The "standard deviation" of the fit's residuals. If BISECTOR 
		is set, this will be smaller by ~ sqrt(2).
	COEF_SIG  = The estimated standard deviations of the coefficients. If 
		BISECTOR is set, however, this becomes the vector of fit 
		residuals measured orthogonal to the line.

 OPTIONAL INPUT KEYWORDS:
	NUMIT = the number of iterations allowed. Default = 25
	BISECT  if set, the bisector of the "Y vs X" and "X vs Y" fits is 
		determined.  The distance PERPENDICULAR to this line is used 
		in calculating weights. This is better when the uncertainties 
		in X and Y are comparable, so there is no true independent 
		variable.  Bisquare_Limit  Limit used for calculation of 
		bisquare weights. In units of outlier-resistant standard 
		deviations. Default: 6.
		Smaller limit ==>more resistant, less efficient
 Close_Factor  - Factor used to determine when the calculation has converged.
		Convergence if the computed standard deviation changes by less 
		than Close_Factor * ( uncertainty of the std dev of a normal
		distribution ). Default: 0.03.
 SUBROUTINE CALLS:
	MED, to calculate the median
	ROB_CHECKFIT
	ROBUST_SIGMA, to calculate a robust analog to the std. deviation

 PROCEDURE:
	For the initial estimate, the data is sorted by X and broken into 2
	groups. A line is fitted to the x and y medians of each group.
	Bisquare ("Tukey's Biweight") weights are then calculated, using the 
	a limit of 6 outlier-resistant standard deviations.
	This is done iteratively until the standard deviation changes by less 
	than CLOSE_ENOUGH = CLOSE_FACTOR * {uncertainty of the standard 
		deviation of a normal distribution}

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 4/91.
	4/13/93 to return more realistic SS's HF
	2/94 --more error-checking, changed convergence criterion HF
	5/94 --added BISECT option. HF.
       8/94 --added Close_Factor and Bisquare_Limit options  Jack Saba.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_linefit.pro)


ROBUST_PLANEFIT

[Previous Routine] [Next Routine] [List of Routines]
NAME:
	ROBUST_PLANEFIT

PURPOSE:
	An outlier-resistant fit to Z = a + bX + cY

 CALLING SEQUENCE:
	COEFF = ROBUST_PLANEFIT(X,Y,Z, ZFIT,SIG)

 INPUTS:
 X = Independent variable vector. DOUBLE-PRECISION RECOMMENDED
 Y = Second independent variable vector
 Z = Dependent variable vector

OUTPUTS:
 Function result = coefficient vector [a,b,c,d,e]
 Either floating point or double precision.

OPTIONAL OUTPUT PARAMETERS:
 ZFIT = Vector of calculated y's
 SIG  = the std. dev. of the residuals

KEYWORDS:
 NUMIT = the number of iterations allowed. Default = 20.
 FLOOR = the minimum value accepted. Z < FLOOR is ignored. Default=-1.0e20

RESTRICTIONS:
 This routine works best when the number of points >> 6. Minimum=6

PROCEDURE:
 For the initial estimate, a regular least-squares fit is performed.
 Bisquare ("Tukey's Biweight") weights are then calculated and the
 fit it iterated until convergence (see ROB_CHECKFIT).

AUTHOR: H. Freudenreich, HSTX, 5/19/93. 
 Minor modifications: H.F., 3/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_planefit.pro)


ROBUST_POLY_FIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_POLY_FIT

 PURPOSE:
	An outlier-resistant polynomial fit.

 CALLING SEQUENCE:
	COEFF = ROBUST_POLY_FIT(X,Y,NDEGREE  ,[ YFIT,SIG, NUMIT =] )

 INPUTS:
	X = Independent variable vector, floating-point or double-precision
	Y = Dependent variable vector

 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

 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)).

 REVISION HISTORY
	Written, H. Freudenreich, STX, 8/90. Revised 4/91.
	2/94 -- changed convergence criterion

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_poly_fit.pro)


ROBUST_REGRESS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_REGRESS

 PURPOSE:
	An outlier-resistant linear regression. Calls REGRESS.

 CALLING SEQUENCE:
	COEFF = ROBUST_REGRESS(X,Y, YFIT,SIG, [ NUMIT = , FLOOR = ] )

 INPUTS:
	X = Matrix of independent variable vectors, dimensioned 
		(NUM_VAR,NUM_PTS), as in REGRESS
	Y = Dependent variable vector. All Y> NVAR+1 points for the fit to be truly outlier-resistant.

 PROCEDURE:
	This iteratively: calls REGRESS, checks the dispersion of the
	residuals, and computes weights (biweights). 
	This is done until the standard deviation, also calculated using 
	biweights, begins to grow or changes by less than CLOSE_ENOUGH, now set
	to .03 X [uncertainty of the standard deviation of a normal 
	distribution]

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 5/19/93
	H.F. 3/94 --add FLOOR keyword and related modifications.

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_regress.pro)


ROBUST_SIGMA

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROBUST_SIGMA  

 PURPOSE:
	Calculate a resistant estimate of the dispersion of a distribution.
	For an uncontaminated distribution, this is identical to the standard
	deviation.

 CALLING SEQUENCE:
	result = ROBUST_SIGMA( Y, [ /ZERO ] )

 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.

 OUTPUT:
	ROBUST_SIGMA returns the dispersion. In case of failure, returns 
	value of -1.0

 SUBROUTINE CALLS:
	MED, which calculates the median

 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.

 REVSION HISTORY: 
	H. Freudenreich, STX, 8/90

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/robust_sigma.pro)


ROB_CHECKFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROB_CHECKFIT
 PURPOSE:
	Used by ROBUST_... routines to determine the quality of a fit and to
	return biweights.
 CALLING SEQUENCE:
	status = ROB_CHECKFIT( Y, YFIT, EPS, DEL, SIG, FRACDEV, NGOOD, W, B
				BISQUARE_LIMIT = )
 INPUT:
	Y     = the data
	YFIT  = the fit to the data
	EPS   = the "too small" limit
	DEL   = the "close enough" for the fractional median abs. deviations
 RETURNS:
	Integer status. if =1, the fit is considered to have converged

 OUTPUTS:
	SIG   = robust standard deviation analog
	FRACDEV = the fractional median absolute deviation of the residuals
	NGOOD   = the number of input point given non-zero weight in the 
		calculation
	W     = the bisquare weights of Y
	B     = residuals scaled by sigma

 OPTIONAL INPUT KEYWORD:
	BISQUARE_LIMIT = allows changing the bisquare weight limit from 
			default 6.0
 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 1/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/rob_checkfit.pro)


ROB_MAPFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	 ROB_MAPFIT 

 PURPOSE:
	Robustly fit a polynomial surface to a 2D floating-point array

 CALLING SEQUENCE:
	surface = ROB_MAPFIT( map,ndeg, coefficients, Sigma, [/SINGLE, FLOOR= ])

 INPUTS:
	map = array to fit, floating-point
	ndeg= degree of polynomials to be fitted. 1: Z=a+bX+cY
						2: Z=a+bX+cY+eX^2+eY^2+fXY
				3: Z=a+bX+cY+eX^2+eY^2+fX^3+gY^3+hX^2Y+iY^2X
						and so on...
		Maximum degree = 7 (35 terms)

 OUTPUT:
	ROB_MAPFIT returns the calculated fitted surface at each pixel: 
	an array like the input array.

 OPTIONAL INPUT KEYWORDS:
	FLOOR = the minimum value allowed. Pixels <= this are considered vacant.
		The default = -1.0e20
	SINGLE  if set, single precision is used. This should be done only if 
		there is not enough memory for the default double-precision 
		calculations.

 OPTIONAL OUTPUT:
	coefficients = the coefficients of the fit, in the order [a,b,c,...] 
		as above.
	sigma = a robust analog of the std. deviation of the residuals

 SUBROUTINES NEEDED:
	Planefit, Robust_planefit, Quarticfit, Rob_Quarticfit, Rob_Checkfit,
	Robust_Sigma, Med, REGRESS (in USERLIB)

 NOTES:
	This routine will usually give excellent results as long as no more than
	25% of the pixels are contaminated. When more than 25% are affected, it
	becomes harder to obtain the underlying shape. However, it will still do
	a pretty good job even if the data is a mess. But try to keep the degree
	of the fit down. Best if #pts >> degree
 
 REVISION HISTORY: 
	Written,  H.T. Freudenreich, HSTX, 5/20/93

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/rob_mapfit.pro)


ROB_QUARTICFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	ROB_QUARTICFIT

 PURPOSE:
	An outlier-resistant fit to Z = a + bX + cX^2 + dY +eY^2 + fX*Y.

 CALLING SEQUENCE:
	COEFF = ROB_QUARTICFIT(X, Y, Z, ZFIT, SIG, [NUMIT = , FLOOR = ] )

INPUTS:
	X = Independent variable vector. DOUBLE-PRECISION RECOMMENDED
	Y = Second independent variable vector
	Z = Dependent variable vector

 RETURNS:
	Function result = coefficient vector [a,b,c,d,e]
	Either floating point or double precision.

 OPTIONAL INPUT KEYWORDS:
	NUMIT = the number of iterations permitted. Default = 20
	FLOOR = the minimum Z value allowed. Z> 6. Minimum=6

 PROCEDURE:
	For the initial estimate, a regular least-squares fit is performed.
	Bisquare ("Tukey's Biweight") weights are then calculated, using 
	a limit of 6 outlier-resistant standard deviations.
	This is done iteratively until the standard deviation, also calculated
	using biweights, changes by less than CLOSE_ENOUGH, now set to
	.03 X [uncertainty of the standard deviation of a normal distribution]

 REVISION HISTORY:
	Written, H. Freudenreich, STX, 5/19/93. 
	Minor modifications by H.F. 3/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/rob_quarticfit.pro)


RPLOT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	RPLOT
 PURPOSE:  
	Autoscaled plot of glitchy data. Short horizontal bars on the left of 
	the plot denote the allowed Y range.

 CALLING SEQUENCE:
	RPLOT, X, Y, [ CUT= , _EXTRA= ] 
			or
	RPLOT, Y, [ CUT = , _EXTRA = ]

 INPUT ARGUMENTS:
	X   = the independent variable [OPTIONAL]
	Y   = the dependent variable

 OPTIONAL INPUT KEYWORDS:
	CUT = the number of sigmas from the mean at which the data are to be 
	clipped.  (Actually uses median and interquartile range.) Default = 3.0

	RPLOT will accept any keyword accepted by the PLOT command.
 NOTES:  
	The interquartile range is used to estimate the dispersion of the data
       Points out of range are given the value of the ceiling or floor. If
       there is a trend in the data that is large compared to the noise, the
       ends may be cut off.

 EXAMPLES:
	RPLOT,X,Y,CUT=4.0, PSYM=1,SYMSIZ=.5
		or
	RPLOT,Y,CUT=2.5, PSYM=1,SYMSIZ=.5

 REVISION HISTORY:
	H.T. Freudenreich, HSTX, 3/91. _EXTRA added 2/94

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/rplot.pro)


R_DILATE

[Previous Routine] [Next Routine] [List of Routines]
NAME:
	R_DILATE

 PURPOSE:
	Replace pixel values with the neighborhood maximum. Like DILATE, but
	1. takes floating-point maps
	2. ignores pixels <= a minimum value/

 CALLING SEQUENCE:
	dilatedmap = R_DILATE( map, [ kernel, FLOOR = ])

 INPUT ARGUMENT:
	MAP = the array to dilate
 OPTIONAL INPUT:
	KERNEL = the byte-sized kernel that describes the shape of the 
		neighborhood.    For example, for a 3x3 neighborhood, 
		KERNEL=[[1,1,1],[1,1,1],[1,1,1]].    Or you could define a 
		cross-shaped neighborhood by KERNEL = [[0,1,0],[1,1,1],[0,1,0]]
		It should be square and of odd dimension.
 
          The default kernel = 7x7 with the 24 corner pixels zeroed.

 OPTIONAL INPUT KEYWORD:
	 FLOOR = the minimum value of an occupied pixel. Default=0.

 OUTPUT:
	DILATEDMAP - The floating-point dilated map.

NOTE: 
	The edges of the map are not filtered.
	If too few pixels are occupied within a neighborhood, no filtering is 
	done.   The companion routine is R_ERODE

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 5/19/93

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/r_dilate.pro)


R_ERODE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	R_ERODE

 PURPOSE:
	Replace pixel values with the neighborhood minimum. Like ERODE, but
	1. takes floating-point maps
	2. ignores pixels <= MINVAL

 CALLING SEQUENCE:
	erodedmap = r_erode(map,kernel)

 INPUT ARGUMENT:
	MAP = the array to erode
 OPTIONAL INPUT:
	KERNEL = the byte-sized kernel that describes the shape of the 
		neighborhood
          For example, for a 3x3 neighborhood, KERNEL=[[1,1,1],[1,1,1],[1,1,1]]
          Or you could define a cross-shaped neighborhood by
          KERNEL = [[0,1,0],[1,1,1],[0,1,0]]. It should be square and of odd
          dimension.
 
          The default kernel = 7x7 with the 24 corner pixels zeroed.
 OPTIONAL INPUT KEYWORD:
	FLOOR = the minimum valid value. Default = 0.

 OUTPUT: 
	R_ERODE returns the floating-point eroded map.

 NOTES: 
	The edges of the map are not filtered.  If too few pixels are occupied
	 within a neighborhood, no filtering is done.   The companion routine 
	is R_DILATE

 REVISION HISTORY:
	Written, H.T. Freudenreich, HSTX, 5/19/93

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/r_erode.pro)


XYFIT

[Previous Routine] [List of Routines]
 NAME:
	XYFIT
 PURPOSE:
	Returns the bisecting line of the regressions: Y on X and X on Y.
	OR the mean of the two lines, with FX and FY the weights.
	CALLED BY FITXYERRS.
 CALLING SEQUENCE:
	YFIT = XYFIT( X, Y, W, CCX, CCY, DIST, FX, FY )
 INPUTS:
	X,Y,W = points and weights
	FX, FY = weight given to X|Y and Y|X regressions in their average
 OUTPUT:
	CCX = coefficients of Y on X regression
	CCY = same for X on Y
	DIST = perpendicular distance of points to line
 AUTHOR:
	H.T. Freudenreich

(See /host/bluemoon/usr2/idllib/astron/contrib/freudenreich/xyfit.pro)