Extended IDL Help

This page was created by the IDL library routine mk_html_help. For more information on this routine, refer to the IDL Online Help Navigator or type:

     ? mk_html_help

at the IDL command line prompt.

Last modified: Tue Aug 18 15:38:05 2015.


List of Routines


Routine Descriptions

ACOSH

[Next Routine] [List of Routines]
 NAME:
     ACOSH
 PURPOSE:
     Return the inverse hyperbolic sine of the argument
 EXPLANATION:
     The inverse hyperbolic sine is used for the calculation of acosh 
     magnitudes, see Lupton et al. (1999, AJ, 118, 1406)

 CALLING SEQUENCE
     result = acosh( x) 
 INPUTS:
     X - hyperbolic sine, numeric scalar or vector or multidimensional array 
        (not complex) 

 OUTPUT:
     result - inverse hyperbolic sine, same number of elements as X
              double precision if X is double, otherwise floating pt.

 METHOD:
     Expression given in  Numerical Recipes, Press et al. (1992), eq. 5.6.8 

 REVISION HISTORY:
   MODIFICATION HISTORU FOR ASINH
     Written W. Landsman                 February, 2001
     Work for multi-dimensional arrays  W. Landsman    August 2002
     Simplify coding, and work for scalars again  W. Landsman October 2003
   MODIFICATION HISTORU FOR ACOSH
     Modified from asinh, Leslie Young, Aug 2006

(See ../math/acosh.pro)


AMOEBA

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

 PURPOSE:
	Multidimensional minimization of a function FUNC(X), where
	X is an N-dimensional vector, using the downhill simplex
	method of Nelder and Mead, 1965, Computer Journal, Vol 7, pp 308-313.

    This routine is modified from IDL's amoeba routine, which
	in turn is based on the AMOEBA routine, Numerical
	Recipes in C: The Art of Scientific Computing (Second Edition), Page
	411.

 CATEGORY:
	Function minimization/maximization. Simplex method.

 CALLING SEQUENCE:
	Result = AMOEBAFIT(Ftol, data, ....)
 INPUTS:
    FTOL:  the fractional tolerance to be achieved in the function
	value.  e.g. the fractional decrease in the function value in the
	terminating step.  This should never be less than the
	machine's single or double precision.
    DATA: The data values
 KEYWORD PARAMETERS:
    FUNCTION_NAME: a string containing the name of the function to
	be minimized.  If omitted, the function FUNC is minimized.
	This function must accept an Ndim vector as its only parameter and
	return a scalar single or double precision floating point value as its
	result. 
    *************************************************************
    This function is now the same form as for curvefit, eg
    function_name(coords, param, values, deriv)
    *************************************************************

    COORD: The dependent variable, [0,1,2,3...] if not passed
    WEIGHT: The weights for the data, [1,1,1,1,] if not passed

    FUNCTION_VALUE: (output) on exit, an Ndim+1 element vector
	containing the function values at the simplex points.  The first
	element contains the function minimum. 
    
    NCALLS: (output) the of times the function was evaluated. 
    NMAX: the maximum number of function evaluations allowed
	before terminating.  Default = 5000.
    P0: Initial starting point, an Ndim element vector.  The starting
	point must be specified using either the keyword SIMPLEX, or P0 and
	SCALE.  P0 may be either single or double precision floating.
	For example, in a 3-dimensional problem, if the initial guess
	is the point [0,0,0], and it is known that the function's
	minimum value occurs in the interval: -10 <
	X(0) < 10, -100 < X(1) < 100, -200 < X(2) < 200, specify: P0=[0,0,0],
	SCALE=[10, 100, 200]. 
    SCALE: a scalar or Ndim element vector contaiing the problem's
	characteristic length scale for each dimension.
	SCALE is used with P0 to form an initial (Ndim+1) point simplex.
	If all dimensions have the same	scale, specify a scalar.
    SIMPLEX: (output and/or optional input) On input, if P0 and SCALE
	are not set, SIMPLEX contains the Ndim+1 vertices, each of
	Ndim elements, of starting simplex, in either single or double
	precision floating point, in an (Ndim, Ndim+1) array. On output,
	SIMPLEX contains the simplex, of dimensions (Ndim, Ndim+1), enclosing
	the function minimum.  The first point, Simplex(*,0), corresponds to
	the function's minimum.
	METHOD: What to minimize.  Choices are 
       'sumsqres' to minimize the sum of (data-model)^2, 
       'sumabres' to minimize the sum of |data-model|, 
       'maxabdev' to minimize the maximum of |data-model|

 OUTPUTS:
   Result: If the minimum is found, an Ndim vector, corresponding to
	the Function's minimum value is returned.  If a function minimum
	within the given tolerance, is NOT found in the given number of
	evaluations, a scalar value of -1 is returned.

 COMMON BLOCKS:
	None.

 SIDE EFFECTS:
	None.

 PROCEDURE:
	This procedure implements the Simplex method, described in
	Numerical Recipes, Section 10.4.  See also the POWELL procedure.

	Advantages:  requires only function evaluations, not
	derivatives, may be more reliable than the POWELL method.
	Disadvantages: not as efficient as Powell's method, and usually
	requires more function evaluations.

	Results are performed in the mode (single or double precision)
	returned by the user-supplied function.  The mode of the inputs P0,
	SCALE, or SIMPLEX, should match that returned by the function. The
	mode of the input vector supplied to the user-written function, is
	determined by P0, SCALE, or SIMPLEX.

 EXAMPLE:
	Use Amoeba to find the slope and intercept of a straight line fitting
	a given set of points minimizing the maximum error:

	The function to be minimized returns the maximum error,
	given p(0) = intercept, and p(1) = slope:
 FUNCTION FUNC, x, p, y
 RETURN, (p[0] + p[1] * x)
 END

	Define the data points:
   x = findgen(17)*5
   y = [ 12.0,  24.3,  39.6,  51.0,  66.5,  78.4,  92.7, 107.8, 120.0, $
        135.5, 147.5, 161.0, 175.4, 187.4, 202.5, 215.4, 229.9]

	Call the function.  Fractional tolerance = 1 part in 10^5, 
	Initial guess = [0,0], and the minimum should be found within
	a distance of 100 of that point: 
   r = AMOEBAFIT(1.0e-5, y, COORD=x, SCALE=1.0e2, P0 = [0, 0], FUNCTION_VALUE=fval, METHOD='maxabres')

	Check for convergence:
   if n_elements(r) eq 1 then message,'AMOEBA failed to converge'
	Print results.
   print, 'Intercept, Slope:', r, 'Function value (max error): ', fval(0)
Intercept, Slope:      11.4100      2.72800
Function value:       1.33000



 MODIFICATION HISTORY:
	LAY, Oct, 2002.	Modified Amoeba.

(See ../math/amoebafit.pro)


COMPLEXPHASE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  complexphase
 PURPOSE: (one line)
  Return the phase of a complex number
 DESCRIPTION:
  Return the phase of a complex number
 CATEGORY:
  Math
 CALLING SEQUENCE:
  p = complexphase(c)
 INPUTS:
  c - scalar or array; int, real, double or complex
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  atan(im/real)
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000, by Leslie Young, SwRI

(See ../math/complexphase.pro)


COR2COV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  cor2cov
 PURPOSE: (one line)
  Covariance matrix to matric of correlation coefficients
 CATEGORY:
  Math
 CALLING SEQUENCE:
  cor2cov, cor, sigma, cov
 INPUTS:
  cov = covariance matrix
 OUTPUTS:
  sigma = arror array
  cor = correlation coefficients
 RESTRICTIONS:
  cor must be square, symmetric, with positive diagonals
 MODIFICATION HISTORY:
  Oct 2006

(See ../math/cor2cov.pro)


COV2COR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  cov2cor
 PURPOSE: (one line)
  Covariance matrix to matric of correlation coefficients
 CATEGORY:
  Math
 CALLING SEQUENCE:
  cov2cor, cov, sigma, cor
 INPUTS:
  cov = covariance matrix
 OUTPUTS:
  sigma = arror array
  cor = correlation coefficients
 RESTRICTIONS:
  cor must be square, symmetric, with positive diagonals
 MODIFICATION HISTORY:
  Oct 2006

(See ../math/cov2cor.pro)


CUBIC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  cubic
 PURPOSE: (one line)
  evaluate a cubic given values and derivatives at two points
 DESCRIPTION:
  evaluate a cubic given values and derivatives at two points
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = cubic(x0, x1, y0, y1, dy0, dy1, x, dy, dx0, dx1, a)
 INPUTS:
  x0, x1 - the fixed points
  y0, y1 - values at the fixed points
  dy0, dy1 - derivatives at the fixed points, same
             units as (y1-y0)/(x1-x0)
  x - new point at which to evaluate the cubic
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  returns y = value of the cubic at x
 OPTIONAL OUTPUT PARAMETERS:
  dy = derivative of the cubic at x
  dx0 = (x-x0)/(x1-x0) and dx1 = (x-x1)/(x1-x0) - for debug only
  a = coefficients for y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
     (for debug only)
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
  Given a cubic function with the following endpoints 
  y(x0) = y0
  y(x1) = y1
  dy(x0)/dx = dy0
  dy(x1)/dx = dy0
  calculate the coefficients so that
  y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
  dx0 = (x-x0)/(x1-x0)
  dx1 = (x-x1)/(x1-x0)
 MODIFICATION HISTORY:
  Written 2000 October, by Leslie Young, SwRI
  Nov 2000.  Allow non-scalar argument, LAY
  11 Mar 2006.  LAY.  Moved to $idl/layoung/math, check nparams
    before calculating dy
     

(See ../math/cubic.pro)


CURVEFIT

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

 PURPOSE:
       Non-linear least squares fit to a function of an arbitrary
       number of parameters.  The function may be any non-linear
       function.  If available, partial derivatives can be calculated by
       the user function, else this routine will estimate partial derivatives
       with a forward difference approximation.

 CATEGORY:
       E2 - Curve and Surface Fitting.

 CALLING SEQUENCE:
       Result = CURVEFIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
                         ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)

 INPUTS:
       X:  A row vector of independent variables.  This routine does
           not manipulate or use values in X, it simply passes X
           to the user-written function.

       Y:  A row vector containing the dependent variable.

  Weights:  A row vector of weights, the same length as Y.
            For no weighting,
                 Weights(i) = 1.0.
            For instrumental (Gaussian) weighting,
                 Weights(i)=1.0/sigma(i)^2
            For statistical (Poisson)  weighting,
                 Weights(i) = 1.0/y(i), etc.

            For no weighting, set Weights to an undefined variable.

       A:  A vector, with as many elements as the number of terms, that
           contains the initial estimate for each parameter.  IF A is double-
           precision, calculations are performed in double precision,
           otherwise they are performed in single precision. Fitted parameters
           are returned in A.

 KEYWORDS:
	FITA:   A vector, with as many elements as A, which contains a zero for
   		each fixed parameter, and a non-zero value for elements of A to
   		fit. If not supplied, all parameters are taken to be non-fixed.

       FUNCTION_NAME:  The name of the function (actually, a procedure) to
       fit.  IF omitted, "FUNCT" is used. The procedure must be written as
       described under RESTRICTIONS, below.

       ITMAX:  Maximum number of iterations. Default = 20.
       ITER:   The actual number of iterations which were performed
       TOL:    The convergence tolerance. The routine returns when the
               relative decrease in chi-squared is less than TOL in an
               interation. Default = 1.e-3.
       CHI2:   The value of chi-squared on exit (obselete)

       CHISQ:   The value of reduced chi-squared on exit
       NODERIVATIVE:   IF this keyword is set THEN the user procedure will not
               be requested to provide partial derivatives. The partial
               derivatives will be estimated in CURVEFIT using forward
               differences. IF analytical derivatives are available they
               should always be used.

       DOUBLE = Set this keyword to force the calculation to be done in
                double-precision arithmetic.

   STATUS: Set this keyword to a named variable in which to return
           the status of the computation. Possible values are:
           STATUS = 0: The computation was successful.
           STATUS = 1: The computation failed. Chi-square was
                       increasing without bounds.
           STATUS = 2: The computation failed to converge in ITMAX
                       iterations.

   YERROR: The standard error between YFIT and Y.

 OUTPUTS:
       Returns a vector of calculated values.
       A:  A vector of parameters containing fit.

 OPTIONAL OUTPUT PARAMETERS:
       Sigma:  A vector of standard deviations for the parameters in A.

     Note: if Weights is undefined, then you are assuming that
           your model is correct. In this case, SIGMA is multiplied
           by SQRT(CHISQ/(N-M)), where N is the number of points
           in X and M is the number of terms in the fitting function.
           See section 15.2 of Numerical Recipes in C (2nd ed) for details.

 COMMON BLOCKS:
       NONE.

 SIDE EFFECTS:
       None.

 RESTRICTIONS:
       The function to be fit must be defined and called FUNCT,
       unless the FUNCTION_NAME keyword is supplied.  This function,
       (actually written as a procedure) must accept values of
       X (the independent variable), and A (the fitted function's
       parameter values), and return F (the function's value at
       X), and PDER (a 2D array of partial derivatives).
       For an example, see FUNCT in the IDL User's Libaray.
       A call to FUNCT is entered as:
       FUNCT, X, A, F, PDER
 where:
       X = Variable passed into CURVEFIT.  It is the job of the user-written
           function to interpret this variable.
       A = Vector of NTERMS function parameters, input.
       F = Vector of NPOINT values of function, y(i) = funct(x), output.
       PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
               PDER(I,J) = DErivative of function at ith point with
               respect to jth parameter.  Optional output parameter.
               PDER should not be calculated IF the parameter is not
               supplied in call. IF the /NODERIVATIVE keyword is set in the
               call to CURVEFIT THEN the user routine will never need to
               calculate PDER.

 PROCEDURE:
       Copied from "CURFIT", least squares fit to a non-linear
       function, pages 237-239, Bevington, Data Reduction and Error
       Analysis for the Physical Sciences.  This is adapted from:
       Marquardt, "An Algorithm for Least-Squares Estimation of Nonlinear
       Parameters", J. Soc. Ind. Appl. Math., Vol 11, no. 2, pp. 431-441,
       June, 1963.

       "This method is the Gradient-expansion algorithm which
       combines the best features of the gradient search with
       the method of linearizing the fitting function."

       Iterations are performed until the chi square changes by
       only TOL or until ITMAX iterations have been performed.

       The initial guess of the parameter values should be
       as close to the actual values as possible or the solution
       may not converge.

 EXAMPLE:  Fit a function of the form f(x) = a * exp(b*x) + c to
           sample pairs contained in x and y.
           In this example, a=a(0), b=a(1) and c=a(2).
           The partials are easily computed symbolicaly:
           df/da = exp(b*x), df/db = a * x * exp(b*x), and df/dc = 1.0

           Here is the user-written procedure to return F(x) and
           the partials, given x:

       pro gfunct, x, a, f, pder      ; Function + partials
         bx = exp(a(1) * x)
         f= a(0) * bx + a(2)         ;Evaluate the function
         IF N_PARAMS() ge 4 THEN $   ;Return partials?
         pder= [[bx], [a(0) * x * bx], [replicate(1.0, N_ELEMENTS(f))]]
       end

         x=findgen(10)                  ;Define indep & dep variables.
         y=[12.0, 11.0,10.2,9.4,8.7,8.1,7.5,6.9,6.5,6.1]
         Weights=1.0/y            ;Weights
         a=[10.0,-0.1,2.0]        ;Initial guess
         yfit=curvefit(x,y,Weights,a,sigma,function_name='gfunct')
         print, 'Function parameters: ', a
         print, yfit
       end

 MODIFICATION HISTORY:
       Written, DMS, RSI, September, 1982.
       Does not iterate IF the first guess is good.  DMS, Oct, 1990.
       Added CALL_PROCEDURE to make the function's name a parameter.
              (Nov 1990)
       12/14/92 - modified to reflect the changes in the 1991
            edition of Bevington (eq. II-27) (jiy-suggested by CreaSo)
       Mark Rivers, U of Chicago, Feb. 12, 1995
           - Added following keywords: ITMAX, ITER, TOL, CHI2, NODERIVATIVE
             These make the routine much more generally useful.
           - Removed Oct. 1990 modification so the routine does one iteration
             even IF first guess is good. Required to get meaningful output
             for errors.
           - Added forward difference derivative calculations required for
             NODERIVATIVE keyword.
           - Fixed a bug: PDER was passed to user's procedure on first call,
             but was not defined. Thus, user's procedure might not calculate
             it, but the result was THEN used.

      Steve Penton, RSI, June 1996.
            - Changed SIGMAA to SIGMA to be consistant with other fitting
              routines.
            - Changed CHI2 to CHISQ to be consistant with other fitting
              routines.
            - Changed W to Weights to be consistant with other fitting
              routines.
            _ Updated docs regarding weighing.

      Chris Torrence, RSI, Jan,June 2000.
         - Fixed bug: if A only had 1 term, it was passed to user procedure
           as an array. Now ensure it is a scalar.
         - Added more info to error messages.
         - Added /DOUBLE keyword.
      CT, RSI, Nov 2001: If Weights is undefined, then assume no weighting,
           and boost the Sigma error estimates according to NR Sec 15.2
           Added YERROR keyword.
      CT, RSI, May 2003: Added STATUS keyword.
	        Added FITA keyword (courtesy B. LaBonte)
      CT, RSI, August 2004: Added ON_ERROR, 2

      Sep 15 2006 Leslie Young. added keywords alpha, beta, sigmaalt, derivdelta
      Jan 26 2009 Leslie Young. increase parameter increment for
                                numerical derivatives

(See ../math/curvefitlay.pro)


CURVEFITLAYTEST

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

 PURPOSE:
       tester/driver for curvefitlay

 CATEGORY:
       E2 - Curve and Surface Fitting.

 CALLING SEQUENCE:
       CURVEFITLAYTEST

 INPUTS:
 KEYWORDS:
 OUTPUTS:

 OPTIONAL OUTPUT PARAMETERS:

 COMMON BLOCKS:

 SIDE EFFECTS:

 RESTRICTIONS:
 PROCEDURE:


 MODIFICATION HISTORY:
      Leslie Young 2004 Apr 30.
           

(See ../math/curvefitlaytest.pro)


DBRENT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	dbrent
 PURPOSE:
	minimum finding with derivatives using Brent's method
 DESCRIPTION:
  Given a function f and its derivative function df, and given a
  bracketing triplet of abscissas ax, bx, cx [such that bx is between ax
  and cx, and f(bx) is less than both f(ax) and f(cx)], this routine
  isolates the minimum to a fractional precision of about tol using a modi
  cation of Brent's method that uses derivatives. The abscissa of the
  minimum is returned as xmin, and the minimum function value is returned
  as dbrent, the returned function value.
 CALLING SEQUENCE:
	min = dbrent(func, ax,bx,cx, xmin, tol=tol, extras=extras)

 INPUT PARAMETERS
   func: a string containing the name of the function to be minimized. 
     func, x, f, df,   OR  func,x,f,df, p
   ax, bx, cx: a bracketing triplet of abscissas ax, bx, cx 
    [such that bx is between ax and cx, and f(bx) is less than both f(ax) and f(cx)]

 OPTIONAL INPUT PARAMETERS
   tol - tolerance
   funcp - optional other parameters to be passed to f
     

 OUTPUT PARAMETERS
	xmin abscissa of the minimum
   dbrent - minimum function value, returned

 PROCEDURE:
   Based on dbrent in Numerical Recipes

 REVISION HISTORY

(See ../math/dbrent.pro)


EXPINTEI

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  expintei
 PURPOSE: (one line)
  Return the exponential integral, Ei
 DESCRIPTION:
  Return the exponential integral, Ei
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = expintei(x)
 INPUTS:
  x - argument, real number > 0
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  Ei(x)
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  x is a scalar
 PROCEDURE:
  Based on ei.c, Numerical Recipes 2nd Ed, pp 225.
 MODIFICATION HISTORY:
  Written 2006 Sept, by Leslie Young, SwRI

(See ../math/expintei.pro)


EXPM1

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  expm1
 PURPOSE: (one line)
  Return exp(x)-1, even for small x.
 DESCRIPTION:
  Return exp(x)-1, even for small x.
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = expm1(x)
 INPUTS:
  x - exponent
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  exp(x)-1
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000 October, by Leslie Young, SwRI
  Nov 2000.  Allow non-scalar argument, LAY

(See ../math/expm1.pro)


FMAXLOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	fMaxLoc
 PURPOSE: (one line)
	Returns the location of the maximum in a 1-D array
 DESCRIPTION:
	This function returns finds the maximum point in a 1-D array, then fits
	a parabola to the maximum and its left and right neighbors. The center of
	the parabola is assumed to be the fractional location of the maximum. 
	NOTICE that this function is useless if the FWHM is small compared to a pixel.
 CATEGORY:
	Spectral extraction
 CALLING SEQUENCE:
	fMax = fMaxLoc(v)
 INPUTS:
	myVec	- Input 1-D array.
 OPTIONAL INPUT PARAMETERS:
 KEYWORD PARAMETERS:
 OUTPUTS:
	fMax	- the location of the maximum
 COMMON BLOCKS:
 SIDE EFFECTS:
 RESTRICTIONS:
 PROCEDURE:
 MODIFICATION HISTORY:
	Written June 6, 1996, Eliot Young, NASA Ames Research Center
	Modified June 22, 1999, Jason Cook & Leslie Young, BU.
		Found y1=0, y3=0, and arrWidth bugs.
	Modified Dec 21, 1999, Leslie Young, SwRI.
		Found bracketing bug

(See ../math/fmaxloc.pro)


GCF

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  gcf
 PURPOSE: (one line)
  Returns the incomplete gamma function Q(a; x)
 DESCRIPTION:
 Returns the incomplete gamma function Q(a; x) evaluated by its continued fraction representation
 as gammcf. Also returns ln Gamma
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = gcf(a,x, gln=gln)
 INPUTS:
  a - first parameter for Q
  x - second parameter for Q
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  gln = ln(Gamma(a))
 OUTPUTS:
  Q(a; x)
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2003 November, by Leslie Young, SwRI
  based on Numerical Recipes gcf, Section 6.2

(See ../math/gcf.pro)


IM

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
     im
 PURPOSE:
     Return the imaginary part of a number
 EXPLANATION:
 CALLING SEQUENCE
     r = im(x)
 INPUTS:
     X - number or array of type complex, dcomplex
         (other types returned as 0)
 OUTPUT:
     result - imaginary part of x
        float for x is complex
        double for x is dcomplex
        original type for all others
 METHOD:
     IDL routines float and dfloat;
 REVISION HISTORY:
     2012-01-08 Leslie Young, SwRI

(See ../math/im.pro)


LEGPOLY

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

 PURPOSE:
	Evaluate a sum of legendre polynomial functions of a variable.

 CATEGORY:
	C1 - Operations on polynomials.

 CALLING SEQUENCE:
	Result = LEGPOLY(X,C)

 INPUTS:
	X:	The variable.  This value can be a scalar, vector or array.

	C:	The vector of legendre polynomial coefficients.  The degree of
		of the polynomial is N_ELEMENTS(C) - 1.

 OUTPUTS:
	LEGPOLY returns a result equal to:
		 C[0] + c[1] * LEGENDRE(X,1) + c[2]*LEGENDRE(x,2) + ...

 COMMON BLOCKS:
	None.

 SIDE EFFECTS:
	None.

 RESTRICTIONS:
	Designed for X between -1 and 1.

 PROCEDURE:
	Straightforward.

 MODIFICATION HISTORY:
       LAY, Written, Dec 2010

	BASED ON RSI FUNCTION POLY
 $Id: //depot/idl/IDL_70/idldir/lib/poly.pro#2 $

 Copyright (c) 1983-2008, ITT Visual Information Solutions. All
       rights reserved. Unauthorized reproduction is prohibited.
	DMS, Written, January, 1983.
   CT, RSI, Nov 2004: Special case for zero-order polynomial.
       Be sure to still return an array.

(See ../math/legpoly.pro)


LEGPOLY_FIT

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

 PURPOSE:
   Perform a least-square fit with optional error estimates to
   legendre coefficients

   This routine uses matrix inversion.  A newer version of this routine,
   SVDFIT, uses Singular Value Decomposition.  The SVD technique is more
   flexible, but slower.

 CATEGORY:
   Curve fitting.

 CALLING SEQUENCE:
   Result = LEGPOLY_FIT(X, Y, Degree)

 INPUTS:
   X:  The independent variable vector.

   Y:  The dependent variable vector, should be same length as x.

   Degree: The degree of the polynomial to fit.

 OUTPUTS:
   LEGPOLY_FIT returns a vector of coefficients with a length of NDegree+1.

 KEYWORDS:
   CHISQ:   Sum of squared errors divided by MEASURE_ERRORS if specified.

   COVAR:   Covariance matrix of the coefficients.

   DOUBLE:  if set, force computations to be in double precision.

   MEASURE_ERRORS: Set this keyword to a vector containing standard
       measurement errors for each point Y[i].  This vector must be the same
       length as X and Y.

     Note - For Gaussian errors (e.g. instrumental uncertainties),
        MEASURE_ERRORS should be set to the standard
        deviations of each point in Y. For Poisson or statistical weighting
        MEASURE_ERRORS should be set to sqrt(Y).

   SIGMA:   The 1-sigma error estimates of the returned parameters.

     Note: if MEASURE_ERRORS is omitted, then you are assuming that
           your model is correct. In this case, SIGMA is multiplied
           by SQRT(CHISQ/(N-M)), where N is the number of points
           in X and M is the number of terms in the fitting function.
           See section 15.2 of Numerical Recipes in C (2nd ed) for details.

   STATUS = Set this keyword to a named variable to receive the status
          of the operation. Possible status values are:
          0 for successful completion, 1 for a singular array (which
          indicates that the inversion is invalid), and 2 which is a
          warning that a small pivot element was used and that significant
          accuracy was probably lost.

    Note: if STATUS is not specified then any error messages will be output
          to the screen.

   YBAND:  1 standard deviation error estimate for each point.

   YERROR: The standard error between YFIT and Y.

   YFIT:   Vector of calculated Y's. These values have an error
           of + or - YBAND.

 COMMON BLOCKS:
   None.

 SIDE EFFECTS:
   None.

 MODIFICATION HISTORY:
   2010 Dec Leslie Young.  Modify from polyfit.
   
 $Id: //depot/idl/IDL_70/idldir/lib/poly_fit.pro#1 $

 Distributed by ITT Visual Information Solutions.

   Written by: George Lawrence, LASP, University of Colorado,
       December, 1981.

   Adapted to VAX IDL by: David Stern, Jan, 1982.
       Modified:    GGS, RSI, March 1996
                    Corrected a condition which explicitly converted all
                    internal variables to single-precision float.
                    Added support for double-precision inputs.
                    Added a check for singular array inversion.
            SVP, RSI, June 1996
                     Changed A to Corrm to match IDL5.0 docs.
                    S. Lett, RSI, December 1997
                     Changed inversion status check to check only for
                     numerically singular matrix.
                    S. Lett, RSI, March 1998
                     Initialize local copy of the independent variable
                     to be of type DOUBLE when working in double precision.
       CT, RSI, March 2000: Changed to call POLYFITW.
       CT, RSI, July-Aug 2000: Removed call to POLYFITW,
                   added MEASURE_ERRORS keyword,
                   added all other keywords (except DOUBLE),
                   made output arguments obsolete.
       CT, RSI, Jan 2003: Combine some vector math expressions,
                   general code cleanup. About 50% faster.

(See ../math/legpoly_fit.pro)


MADLINE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  madline
 PURPOSE: (one line)
  Fits y = a + bx by the criterion of minimum absolute deviations (MAD).
 DESCRIPTION:
   Fits y = a + bx by the criterion of least absolute deviations. The arrays x[0..ndata-1] and
   y[0..ndata-1] are the input experimental points. The fitted parameters a and b are output,
   along with abdev, which is the mean absolute deviation (in y) of the experimental points from
   the fitted line. 
 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  [a,b] = madline(x, y, yfit, abdev)
 INPUTS:
  x  - An n-element vector of independent variables.
  y  - A vector of dependent variables, the same length as X.
 OPTIONAL INPUT PARAMETERS:
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
 OUTPUTS:
  returns [a,b], the fitted parameters for y = a + b * x
  Yfit - A named variable that will contain the vector of calculated Y values. 
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:


 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
     based on Numerical Recipes medline.

(See ../math/madline.pro)


MADLINEW

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  madlinew
 PURPOSE: (one line)
  Fits y = a + bx by the criterion of minimum absolute deviations (MAD).
 DESCRIPTION:
   Fits y = a + bx by the criterion of least absolute deviations. The arrays x[0..ndata-1] and
   y[0..ndata-1] are the input experimental points. The fitted parameters a and b are output,
   which minimize the sum of (w * |y - a - bx| ). 
 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  [a,b] = madlinew(x, y, w, yfit)
 INPUTS:
  x  - An n-element vector of independent variables.
  y  - A vector of dependent variables, the same length as X.
  w  - weights
 OPTIONAL INPUT PARAMETERS:
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
 OUTPUTS:
  returns [a,b], the fitted parameters for y = a + b * x
  Yfit - A named variable that will contain the vector of calculated Y values. 
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:


 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
     based on Numerical Recipes medfit.

(See ../math/madlinew.pro)


MADSCALE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  madscale
 PURPOSE: (one line)
  minimize the average deviation (MAD) for y = a * x.
 DESCRIPTION:
  This routine estimates the scaling factor, a, such that
  the average deviation, total(abs(y-a*x) ), is minimized.

  The statistics are robust in that the result is insensitive
  to outliers.  

  Note that for y[i] = a * x[i], x can be thought of as
  either a coordinate, or a template (e.g., a line profile).

 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  a = madscale(x, y)
 INPUTS:
  x   - dependent variable.
  y   - data, proportional to x.
 OPTIONAL INPUT PARAMETERS:
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
 OUTPUTS:
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:

 The goal is to minimize f = total(abs(y-a*x)) = total( abs(x) * abs(y/x - a) )
 This is equivalent to finding the weighted median of y/x, with
 abs(x) as the weights.

 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.

(See ../math/madscale.pro)


MADSCALEW

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  madscalew
 PURPOSE: (one line)
  minimize the average deviation (MAD) for y = a * x.
 DESCRIPTION:
  This routine estimates the scaling factor, a, such that
  the average deviation, total(w * abs(y-a*x) ), is minimized.

  The statistics are robust in that the result is insensitive
  to outliers.  

  Note that for y[i] = a * x[i], x can be thought of as
  either a coordinate, or a template (e.g., a line profile).

 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  a = madscale(x, y)
 INPUTS:
  x   - dependent variable.
  y   - data, proportional to x.
 OPTIONAL INPUT PARAMETERS:
  w   - weights for y
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
 OUTPUTS:
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:

 The goal is to minimize f = total(w * abs(y-a*x)) = total( w * abs(x) * abs(y/x - a) )
 This is equivalent to finding the weighted median of y/x, with
 abs(x) as the weights.

 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Institute, 2002 Oct 14.

(See ../math/madscalew.pro)


MEDIANW

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  medianw
 PURPOSE: (one line)
  minimize the weighted average deviation.
 DESCRIPTION:
  This routine estimates the average data value, xmed, such that
  the average deviation, total(weight * abs(x-xmed) ), is minimized.
  For equally-weighted points, this is the median.

  The statistics are robust in that the result is insensitive
  to outliers.  
 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  mean = medianw(x, w)
 INPUTS:
  x   - Input data to be analyzed.
 OPTIONAL INPUT PARAMETERS:
  w - weights.  Assumed equally weighted if not passed.
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
 OUTPUTS:
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:


 Minimizing f(xmed) = total(weight * abs(x-xmed)) is equivalent to finding
 df(xmed)/d xmed = 0 = -total(weight * sgn(x-xmed)).  
 If x is sorted, and xmed=x[imed], this can be written as
 0 = total(weight[0:imed-1] * sgn(x[0:imed-1]-xmed)) 
   + total(weight[imed+1:n-1] * sgn(x[imed+1:n-1]-xmed))
 or
 0 = - total(weight[0:imed-1]) + total(weight[imed+1:n-1])

 or
 total(weight[0:imed-1]) = total(weight[imed+1:n-1])


 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Institute, 2002 Jul 23.
  Modified LAY Oct 2002.  Better medianwTEST; 
            Avoid eq for comparing floats in test for one median value
            (equivalent to an odd number of equally weighted values)
            Average two neighbors if xmed is not one of the listed x's
            (equivalent to an even number of equally weighted values)
 Modified LAY Oct 2002.  Changed name to medianw; 

(See ../math/medianw.pro)


NORMV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  normv
 PURPOSE: (one line)
  Return a normalized vector
 DESCRIPTION:
  Return a normalized vector
 CATEGORY:
  Math
 CALLING SEQUENCE:
  n = normv(v)
 INPUTS:
  v - a vector of arbitrary length
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  v/vabs(v)
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000, by Leslie Young, SwRI

(See ../math/normv.pro)


PIECECUB_INTERP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:   
       PIECECUB_INTERP  
 PURPOSE: 
       interpolate 1-d data and derivatives with piece-wise cubic functions
 EXPLANATION:
       Interpolation into a function with continuous values and
       derivatives, with tabulated values.  Between grid points, the
       function is defined by a cubic polynomial.  Unlike cubic
       splines, the second derivative is not necessarily continuous.

       Calculations are done in double precision.

 CALLING SEQUENCE:
       piececub_interp, Xtab, Ytab, dYtab, Xint, Yint, dYint

 INPUT PARAMETERS: 
       Xtab -  Vector containing the current independent variable grid.
               Must be monotonic increasing or decreasing
       Ytab -  Vector containing the current dependent variable values at 
               the XTAB grid points.
       dYtab - Vector containing the current derivative of the dependent 
               variable values at the XTAB grid points.
       Xint -  Scalar or vector containing the new independent variable grid 
               points for which interpolated value(s) of the dependent 
               variable are sought.    Note that -- due to a limitation of the
               intrinsic INTERPOLATE() function -- Xint is always converted to
               floating point internally.

 OUTPUT PARAMETERS:
       Yint  - Scalar or vector with the interpolated value(s) of the 
               dependent variable at the XINT grid points.
               YINT is double precision if Ytab is, float otherwise.
       dYint - Scalar or vector with the interpolated value(s) of the 
               derivative at the XINT grid points
               DYINT is double precision if DYtab is, float otherwise.


 OPTIONAL INPUT KEYWORD:

 EXAMPLE:

 PROCEDURES CALLED:
       TABINV, ZPARCHECK
 PROCEEDURE


 ________ DEFINITIONS ___________________
 Start with some definitions:
 
 x   the independent variable
 x0  x at node 0 (or node n)
 x1  x at node 1 (or node n+1)
 y   the dependent variable
 y0  y at node 0
 dy0 dy/dx at node 0
 y1  y at node 0
 dy1 dy/dx at node 0
 
 Let's define some things to make our lives easier:
 d0  = x-x0 
 d1  = x-x1 
 del = x1-x0
 so d1 and d0 are functions of x.
 
 d0(x0) = 0      d0(x1) = del
 d1(x0) = -del   d1(x1) = 0
 
 Using d0, d1, and del, define four functions:
 
 P0(x) =  (3*d0-d1) * d1^2 / del^3
 P1(x) = -(3*d1-d0) * d0^2 / del^3
 Q0(x) = d0 * d1^2 / del^2
 Q1(x) = d1 * d0^2 / del^2
 
 with the following derivatives 
 
 dP0(x)/dx =  6*d0*d1/del^3
 dP1(x)/dx = -6*d0*d1/del^3
 dQ0(x)/dx =  d1*(d1+2*d0)/del^2
 dQ1(x)/dx =  d0*(d0+2*d1)/del^2
 
 These functions have very tidy values and derivatives at the nodes:
 
 P0(x0) = 1   P0(x1) = 0   dP0(x0)/dx = 0   P0(x1)/dx = 0
 P1(x0) = 0   P1(x1) = 1   dP1(x0)/dx = 0   P1(x1)/dx = 0
 Q0(x0) = 0   Q0(x1) = 0   dQ0(x0)/dx = 1   Q0(x1)/dx = 0
 Q1(x0) = 0   Q1(x1) = 0   dQ1(x0)/dx = 0   Q1(x1)/dx = 1
 
 ________ CALCULATING THE CUBIC FUNCTION ___________________
 This makes it trivial to define the cubic function that
 matches the nodes:
 
 y(x) = y0 * P0(x) + y1 * P1(x) + dy0 * Q0(x) + dy1 * Q1(x)  [CUBIC]
 

 MODIFICATION HISTORY:

(See ../math/piececub_interp.pro)


QINTERP[1]

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  qinterp
 PURPOSE: (one line)
 given a quadratic through xx, yy, and a pt (or array) x, give val & 2 derive
 DESCRIPTION:
 given a quadratic through xx, yy, and a pt (or array) x, give val & 2 derive
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = qinterp(xx, yy, x, dydx, d2ydx2
 INPUTS:
  xx - absissa (3 elements long)
  yy - value (3 elements long)
  x - absissa for evaluation
 OUTPUTS:
  returns y(x)
  dydx - first derivative
  d2ydx2 - second derivative
 RESTRICTIONS:
  None
 PROCEDURE:
  Lagrange formula (e.g., Numerical Recipies, 2nd ed, 3.1.1)
 MODIFICATION HISTORY:
  Written 2009 Feb 15, by Leslie Young, SwRI

(See ../math/qinterp.pro)


QINTERP[2]

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  qinterp
 PURPOSE: (one line)
 given sorted arrays x and y, find 2 derivatives at x
 DESCRIPTION:
 CATEGORY:
  Math
 CALLING SEQUENCE:
  qinterp_arr,x, y, dydx, d2ydx2
 INPUTS:
  x - absissa (sorted)
  y - value (3 elements long)
 OUTPUTS:
  dydx - first derivative at x
  d2ydx2 - second derivative at x
 RESTRICTIONS:
  None
 PROCEDURE:
  Lagrange formula (e.g., Numerical Recipies, 2nd ed, 3.1.1)
 MODIFICATION HISTORY:
  Written 2009 Feb 15, by Leslie Young, SwRI

(See ../math/qinterp_arr.pro)


RE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
     re
 PURPOSE:
     Return the real part of a number
 EXPLANATION:
 CALLING SEQUENCE
     r = re(x)
 INPUTS:
     X - number or array of type complex, dcomplex
         (other types returned unchanged)
 OUTPUT:
     result - real part of x
        float for x is complex
        double for x is dcomplex
        original type for all others
 METHOD:
     IDL routines float and dfloat;
 REVISION HISTORY:
     2012-01-08 Leslie Young, SwRI

(See ../math/re.pro)


ROBOCUBE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	robocube
 PURPOSE: (one line)
	Combine arrays with a median average.
 DESCRIPTION:
       This will combine a series of arrays into a single array by filling each
       pixel in the output array with the robust mean of the corresponding pixels
       in the input arrays. SAME AS BUIE'S AVCLIP, but probably slower.
       It calls
	robomean,data,thresh,eps,avg,avgdev,stddev,var,skew,kurt,nfinal
 CATEGORY:
	CCD data processing
 CALLING SEQUENCE:
       robocube, inarr, outarr
 INPUTS:
       inarr  -- A three dimensional array containing the input arrays to
                 combine together.  Each of the input arrays must be two
                 dimensional and must have the same dimensions.  These arrays
                 should then be stacked together into a single 3-D array,
                 creating INARR.
 OPTIONAL INPUT PARAMETERS:
	None.
 KEYWORD PARAMETERS:
	None.
 OUTPUTS:
       outarr -- The output array.  It will have dimensions equal to the
                 first two dimensions of the input array.
 COMMON BLOCKS:
	None.
 SIDE EFFECTS:
	None.
 RESTRICTIONS:
	This will run VERY slow if inarr and outarr won't fit into real memory
	on your computer.  Don't try this using virtual memory.
 PROCEDURE:
       The output array is created and then each pixel is extracted from the
	cube.  Once extracted, the pixel stack is sorted and the middle value
	is put into the output array.
 MODIFICATION HISTORY:
       Written by Leslie Young, Boston University, 28 June 1996
       based on medarr, which was
	Written by Marc W. Buie, Lowell Observatory, 17 January 1992

(See ../math/robocube.pro)


ROBOFIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  robofit
 PURPOSE: (one line)
  Least-squares fitting with the identification of outliers.

 DESCRIPTION:
       Robust non-linear least squares fit to a function of an arbitrary 
       number of parameters.  The function may be any non-linear 
       function.  If available, partial derivatives can be calculated by 
       the user function, else this routine will estimate partial derivatives
       with a forward difference approximation.

       This routine differs from curvefit in that it assumes that
       the distribution may be non-normal -- having an "outerlier"
       tail.  The function is first fit to the data minimizing
       the sum of absolute deviations.  From this solution,
       outliers are identified, and expunged from the final
       least-squares fit.

       To maintain compatibility with CURVEFIT, the calling
       sequence is the same (with the addition of THRESH).
 CATEGORY:
       E2 - Curve and Surface Fitting.
 CALLING SEQUENCE:
       Result = ROBOFIT(X, Y, Weights, A, SIGMA, FUNCTION_NAME = name, $
                         ITMAX=ITMAX, ITER=ITER, TOL=TOL, /NODERIVATIVE)
 INPUTS:
 INPUTS:
       X:  A row vector of independent variables.  This routine does
           not manipulate or use values in X, it simply passes X
           to the user-written function.

       Y:  A row vector containing the dependent variable.

  Weights:  A row vector of weights, the same length as Y.
            For no weighting,
                 Weights(i) = 1.0.
            For instrumental (Gaussian) weighting,
                 Weights(i)=1.0/sigma(i)^2
            For statistical (Poisson)  weighting,
                 Weights(i) = 1.0/y(i), etc.

       A:  A vector, with as many elements as the number of terms, that 
           contains the initial estimate for each parameter.  IF A is double-
           precision, calculations are performed in double precision, 
           otherwise they are performed in single precision. Fitted parameters
           are returned in A.
 KEYWORDS:
       FUNCTION_NAME:  The name of the function (actually, a procedure) to 
       fit.  IF omitted, "FUNCT" is used. The procedure must be written as
       described under RESTRICTIONS, below.

       ITMAX:  Maximum number of iterations. Default = 20.
       ITER:   The actual number of iterations which were performed
       TOL:    The convergence tolerance. The routine returns when the
               relative decrease in chi-squared is less than TOL in an 
               interation. Default = 1.e-3.
       CHI2:   The value of chi-squared on exit (obselete)
     
       CHISQ:   The value of reduced chi-squared on exit
       NODERIVATIVE:   IF this keyword is set THEN the user procedure will not
               be requested to provide partial derivatives. The partial
               derivatives will be estimated in CURVEFIT using forward
               differences. IF analytical derivatives are available they
               should always be used.
 OUTPUTS:
       Returns a vector of calculated values.
       A:  A vector of parameters containing fit.
 OPTIONAL OUTPUT PARAMETERS:
       Sigma:  A vector of standard deviations for the parameters in A.
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
       The function to be fit is defined identically to that used by CURVEFIT.

       The function to be fit must be defined and called FUNCT,
       unless the FUNCTION_NAME keyword is supplied.  This function,
       (actually written as a procedure) must accept values of
       X (the independent variable), and A (the fitted function's
       parameter values), and return F (the function's value at
       X), and PDER (a 2D array of partial derivatives).
       For an example, see FUNCT in the IDL User's Libaray.
       A call to FUNCT is entered as:
       FUNCT, X, A, F, PDER
 where:
       X = Variable passed into CURVEFIT.  It is the job of the user-written
           function to interpret this variable.
       A = Vector of NTERMS function parameters, input.
       F = Vector of NPOINT values of function, y(i) = funct(x), output.
       PDER = Array, (NPOINT, NTERMS), of partial derivatives of funct.
               PDER(I,J) = DErivative of function at ith point with
               respect to jth parameter.  Optional output parameter.
               PDER should not be calculated IF the parameter is not
               supplied in call. IF the /NODERIVATIVE keyword is set in the
               call to CURVEFIT THEN the user routine will never need to
               calculate PDER.
 PROCEDURE:
  Call AMOEBA

 MODIFICATION HISTORY:
  Written by Leslie A. Young, Soutwest Research Instituter, 2002 Oct 29.
     based on Numerical Recipes medfit.

(See ../math/robofit.pro)


ROBOLINE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  roboline
 PURPOSE: (one line)
  Robust fitting of a line
 DESCRIPTION:
   Robust fitting of y = a = b * x
 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  [a,b] = roboline(x, y, sigma, thresh, bplist, GDINDEX = gdindex, coeferr = coefErr)
 INPUTS:
  x - the indices
  y - the data array, same length as x
  sigma - estimated noise for each pixel in y
  thresh - threshold for outliers 
 OPTIONAL INPUT PARAMETERS:
  bplist - optional bad pixel list (-1 for all OK)
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
  gdindex - named variable that will contain an array of indices used
 OUTPUTS:
  returns scaleFac such that y ~= a + b * x
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:
   First iteration finds minimum absolute deviation.
   Outliers are identified after the first iteration.
   The second iteration minimizes least squares 

 MODIFICATION HISTORY:
  Written July 22, 2002 Eliot Young, SwRI
  Rewritten by Leslie A. Young, SwRI, 2002 Oct 29.
  - first iteration minimizes absolute deviation
  - exchanged order of x and y in argument list
  - only two iterations
  - optionally returns the error and the list of good pixel indices

(See ../math/roboline.pro)


ROBOMARG

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	robomarg
 PURPOSE: (one line)
	Return the robust averages of each row or each column.
 DESCRIPTION:
	This routine calls robomean on a row-by-row or col-by-col basis
	to get the robust mean of each row or each column. The /GetRowMeans
	keyword returns a column containing the robust average of all of
	the rows. The /GetColMeans keyword returns a row with the robust
	averages of all the columns. /GetRowMeans is the default.
 CATEGORY:
	spectral extraction
 CALLING SEQUENCE:
	means = robomarg(array)
 INPUTS:
	array - Input array to be scanned.
 OPTIONAL INPUT PARAMETERS:
 KEYWORD INPUT PARAMETERS:
	GetRowMeans	- FLAG, return a column which contains the 
					robust means of each row.
	GetColMeans	- FLAG, return a row which contains the 
					robust means of each column.
     bad            - bad pixel mask (1 = bad, 0 = good)
 OUTPUTS:
	means	- The 1-D column or row with averages
 KEYWORD OUTPUT PARAMETERS:
 COMMON BLOCKS:
 SIDE EFFECTS:
 RESTRICTIONS:
 PROCEDURE:
 MODIFICATION HISTORY:
	96/06/06,	EFY, NASA Ames Research Center
	96/10/17,	EFY, changed to a function that returns means
	01/11/02,	LAY, add bad pixel mask option

(See ../math/robomarg.pro)


ROBOSCALE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  roboscale
 PURPOSE: (one line)
  Robust scaling of a template to match an array.
 DESCRIPTION:
   Robust fitting of y = scaleFac * x
 CATEGORY:
  Statistics
 CALLING SEQUENCE:
  scaleFac = roboscale(x, y, sigma, thresh, bplist, )
 INPUTS:
  x - the template
  y - the data array, same length as x
  sigma - estimated noise for each pixel in y
  thesh - threshold for outliers 
  bplist - optional bad pixel list (-1 for all OK)
 OPTIONAL INPUT PARAMETERS:
 INPUT KEYWORD PARAMETERS:
 OUTPUT KEYWORD PARAMETERS:
  gdindex - named variable that will contain an array of indices used
 OUTPUTS:
  returns scaleFac such that y ~= scaleFac * x
 COMMON BLOCKS:
  None.
 SIDE EFFECTS:
  None.
 RESTRICTIONS:
  None.
 PROCEDURE:
   First iteration finds minimum absolute deviation.
   Outliers are identified after the first iteration.
   The second iteration minimizes least squares 

 MODIFICATION HISTORY:
  Written July 22, 2002 Eliot Young, SwRI
  Rewritten by Leslie A. Young, SwRI, 2002 Oct 29.
  - first iteration minimizes absolute deviation
  - exchanged order of x and y in argument list
  - only two iterations
  - optionally returns the error and the list of good pixel indices

(See ../math/roboscale.pro)


SGN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  sgn
 PURPOSE: (one line)
  return the sign of a number or list
 DESCRIPTION:
 CATEGORY:
  Math
 CALLING SEQUENCE:
  y = sgn(x)
 INPUTS:
  x - argument
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  exp(x)-1
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2002 August, by Leslie Young, SwRI

(See ../math/sgn.pro)


VABS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vabs
 PURPOSE: (one line)
  Return the length of a vector
 DESCRIPTION:
  Return the length of a vector
 CATEGORY:
  Math
 CALLING SEQUENCE:
  len = abs(v)
 INPUTS:
  v - a vector of arbitrary length
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  sqrt(total(v^2))
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000, by Leslie Young, SwRI
  2006 Jan 12 LAY change to double

(See ../math/vabs.pro)


VANG

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vang
 PURPOSE: (one line)
  Return the angle between two vectors
 DESCRIPTION:
  Return the angle (in radians) between two vectors
 CATEGORY:
  Math
 CALLING SEQUENCE:
  d = vang(v1,v2)
 INPUTS:
  v1, v2 - two vector of arbitrary length
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  angle between v1 and v2
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000, by Leslie Young, SwRI

(See ../math/vang.pro)


VDOT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vdot
 PURPOSE: (one line)
  Return the dot product of two vectors
 DESCRIPTION:
  Return the dot product of two vectors
 CATEGORY:
  Math
 CALLING SEQUENCE:
  d = vdot(v1,v2)
 INPUTS:
  v1, v2 - two vector of arbitrary length
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  v1 . v2
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2000, by Leslie Young, SwRI

(See ../math/vdot.pro)


WMEAN

[Previous Routine] [List of Routines]
 NAME:
  wmean
 PURPOSE: (one line)
  weighted mean.
 DESCRIPTION:
  Return weighted mean and the error of the mean.
 CATEGORY:
  Math
 CALLING SEQUENCE:
  meanx = wmean(x, xerr, meanerr)
 INPUTS:
  x - array of values
  xerr - array of errors for x
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  exterr - calculate the "external error" using scatter
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  return weighted mean.
  manxerr - error in the mean
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Added to layoung/math 2004 Feb, by Leslie Young, SwRI
  fixed typo 2005 Aug 4 LAY
  2009 Oct 20 LAY - mult. residual by w (weight), to wtot

(See ../math/wmean.pro)