ISOPHOT Interactive Analysis (PIA): Astrolib

MATH


List of Routines


Routine Descriptions

FITEXY

[Next Routine] [List of Routines]
 NAME:
	fitexy
 PURPOSE:
	Linear Least-squares approximation in one-dimension (y = a + b*x),
		when both x and y data have errors (e.g. Gaussian noise).

 CALLING EXAMPLE:
	fitexy, x, y, A, B, X_SIG=sigx, Y_SIG=sigy, [sigma_A_B, chi_sq, q ]

 INPUTS:
	x = array of values for independent variable.
	y = array of data values assumed to be linearly dependent on x.

 REQUIRED INPUT KEYWORDS:
	X_SIGMA = scalar or array specifying the standard deviation of x data.
	Y_SIGMA = scalar or array specifying the standard deviation of y data.

 OPTIONAL INPUT KEYWORD:
	TOLERANCE = desired accuracy of minimum & zero location, default=1.e-3.

 OUTPUTS:
	A_intercept = constant parameter result of linear fit,
	B_slope = slope parameter, so that:
			( A_intercept + B_slope * x ) approximates the y data.
 OPTIONAL OUTPUT:
	sigma_A_B = two element array giving standard deviation of 
		 A_intercept and B_slope parameters, respectively.
                The standard deviations are not meaningful if (i) the
                fit is poor (see parameter q), or (ii) b is so large that
                the data are consistent with a vertical (infinite b) line.
                If the data are consistent with *all* values of b, then
                sigma_A_B = [1e33,e33]  
	chi_sq = resulting minimum Chi-Square of Linear fit, scalar
       q - chi-sq probability, scalar (0-1) giving the probability that
              a correct model would give a value equal or larger than the
              observed chi squared.   A small value of q indicates a poor
              fit, perhaps because the errors are underestimated.

 COMMON:
	common fitexy, communicates the data for computation of chi-square.

 CALLS:
	function chisq_fitexy
	pro minf_bracket
	pro minf_parabolic
	function zbrent
       function chi_sqr1   (from Statistics library)
 PROCEDURE:
	From "Numerical Recipes" column by Press and Teukolsky: 
       in "Computer in Physics",  May, 1992 Vol.6 No.3
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC  September 1992.
       Now returns q rather than 1-q   W. Landsman  December 1992

(See /usr/local/idl/lib/zastron/math/fitexy.pro)


GAMMLN

[Previous Routine] [Next Routine] [List of Routines]
 NAME
	GAMMLN
 PURPOSE:
	Return the natural log of the gamma function.   Adapted from the
	algorithm GAMMLN in Section 6.1 in "Numerical Recipes " (Second
	Edition) by Press  et al. 1992.    This function became obsolete 
	in IDL 2.4.0 when the equivalent LNGAMMA intrinsic function was 
	introduced.

 CALLING SEQUENCE:
	result = gammln ( xx )

 INPUTS:
	xx - numeric scalar or vector for which the log of the gamma function
		will be evaluated.    Must be > 0

 OUTPUT:
	result = alog ( gamma(xx) ).   The result will double precision if xx
		is double, otherwise the result is floating point.
 NOTES:
	IDL has an intrinsic gamma function, GAMMA, but overflow occurs in
	GAMMA for X > 34.5.    By computing the log of the gamma function, one
	can deal with much larger input values.   GAMMLN also allows double 
	precision computation, not available with GAMMA.

 EXAMPLE:
	Compare the output of GAMMA with GAMMLN

	IDL> x = findgen(15)+0.5
	IDL> print, alog(gamma(x))
	IDL> print,gammln(x)

 METHOD:
	Uses the expansion of Lanczos as described in Press et al. (1986)
 REVISION HISTORY:
	Written,   W. Landsman            July, 1992
	Double Precision update           October, 1992

(See /usr/local/idl/lib/zastron/math/gammln.pro)


GAUSSIAN

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

 PURPOSE:
	Compute the 1-D Gaussian function and optionally the derivative 
	at an array of points.

 CALLING:
	y = gaussian( xi, parms,[ pderiv ])

 INPUTS:
	xi = array, independent variable of Gaussian function.

	parms = parameters of Gaussian, 2 or 3 element array:
		parms(0) = maximum value (factor) of Gaussian,
		parms(1) = mean value (center) of Gaussian,
		parms(2) = standard deviation (sigma) of Gaussian.
		(if parms has only 2 elements then sigma taken from common).

 OPTIONAL OUTPUT:
	pderiv = optional output of partial derivatives,
		computed only if parameter is present in call.

		pderiv(*,i) = partial derivative at all xi absisca values
		with respect to parms(i), i=0,1,2.

	Function returns array of Gaussian evaluated at xi.

 EXAMPLE:
	Evaulate a Gaussian centered at x=0, with sigma=1, and a peak value
	of 10 at the points 0.5 and 1.5.   Also compute the derivative

	IDL> f = gaussian( [0.5,1.5], [10,0,1], DERIV )
	==> f= [8.825,3.25].   DERIV will be a 2 x 3 array containing the
	numerical derivative at the two point with respect to the 3 parameters.
 
 COMMON BLOCKS:
	common gaussian, sigma
 HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.

(See /usr/local/idl/lib/zastron/math/gaussian.pro)


KSONE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	KSONE
 PURPOSE:
	Return the Kolmogorov-Smirnov statistic and associated probability for 
	for an array of data values and a user-supplied cumulative distribution
	function (CDF) of a single variable.   Algorithm from the procedure of
	the same name in "Numerical Recipes" by Press et al. 2nd edition (1992)

 CALLING SEQUENCE:
	ksone, data, func_name, D, prob, [ /PLOT ]

 INPUT PARAMATERS:
	data -  vector of data values, must contain at least 4 elements for the
		K-S statistic to be meaningful 
	func_name - scalar string giving the name of the cumulative distribution
		function

 OUTPUT PARAMETERS:
	D - floating scalar giving the Kolmogorov-Smirnov statistic.   It 
		specified the maximum deviation between the cumulative 
		distribution of the data and the supplied function 
	prob - floating scalar between 0 and 1 giving the significance level of
		the K-S statistic.   Small values of PROB show that the 
       	cumulative distribution function of DATA is significantly 
		different from FUNC_NAME.

 OPTIONAL INPUT KEYWORD:
	PLOT - If this keyword is set and non-zero, then KSONE will display a
		plot of the CDF of the data with the supplied function 
		superposed.   The data value where the K-S statistic is 
		computed (i.e. at the maximum difference between the data CDF 
		and the function) is indicated by a vertical line.

 EXAMPLE:
	Determine if a vector created by the RANDOMN function is really 
	consistent with a gaussian distribution.
	The CDF of a gaussian is the error function except that a factor
	of 2 is included in the error function.   So we must create a special
	function:

	function gauss_cdf
	return, errorf( x/sqrt(2) )
	end

	IDL> data = randomn(seed, 50)          ;create data array to be tested
	IDL> ksone, abs(data), 'gauss_cdf', D, prob, /PLOT     ;Use K-S test
      
	PROB gives the probability that DATA came from a gaussian distribution

 NOTES:
	Note that the 2nd (1992) edition of Numerical Recipes includes
	a more accurate computation of the K-S significance for small 
	values of N.

 PROCEDURE CALLS
	procedure prob_ks - computes significance of K-S distribution

 REVISION HISTORY:
	Written     W. Landsman                August, 1992

(See /usr/local/idl/lib/zastron/math/ksone.pro)


KSTWO

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	KSTWO
 PURPOSE:
	Return the Kolmogorov-Smirnov statistic and associated probability 
	that two arrays of data values are drawn from the same distribution
	Algorithm taken from procedure of the same name in "Numerical
	Recipes" by Press et al., 2nd edition (1992)

 CALLING SEQUENCE:
	kstwo, data1, data2, D, prob

 INPUT PARAMATERS:
	data1 -  vector of data values, at least 4 data values must be included
		for the K-S statistic to be meaningful
	data2 -  second set of data values, does not need have the same number
		of elements as data1

 OUTPUT PARAMETERS:
	D - floating scalar giving the Kolmogorov-Smirnov statistic.   It
		specifies the maximum deviation between the cumulative 
		distribution of the data and the supplied function 
	prob - floating scalar between 0 and 1 giving the significance level of
		the K-S statistic.   Small values of PROB show that the 
		cumulative distribution function of DATA1 is significantly 
		different from DATA2

 EXAMPLE:
	Test whether two vectors created by the RANDOMN function likely came
	from the same distribution

	IDL> data1 = randomn(seed,40)        ;Create data vectors to be 
	IDL> data2 = randomn(seed,70)        ;compared
	IDL> kstwo, data1, data2, D, prob   & print,D,prob

 PROCEDURE CALLS
	procedure prob_ks - computes significance of K-S distribution

 REVISION HISTORY:
	Written     W. Landsman                August, 1992

(See /usr/local/idl/lib/zastron/math/kstwo.pro)


LINTERP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:   
	LINTERP  
 PURPOSE: 
	Linearly interpolate tabulated data from one data grid to another.
	Use QUADTERP for quadratic interpolation.

 CALLING SEQUENCE:
	LINTERP, Xtab, Ytab, Xint, Yint, [MISSING = ]   

 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.
	Xint -  Scalar or vector containing the new independent variable grid 
		points for which interpolated value(s) of the dependent 
		variable are sought.

 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 XTAB or YTAB are double,
		otherwise YINT is REAL*4

 OPTIONAL INPUT KEYWORD:
	MISSING - Scalar specifying YINT value(s) to be assigned, when Xint
		value(s) are outside of the range of Xtab.     Default is to
		truncate the out of range YINT value(s) to the nearest value 
		of YTAB.   See the help for the INTERPOLATE function.

 EXAMPLE:
	To linearly interpolate from an IUE spectrum wavelength grid
	WAVE, FLUX to another grid defined as:
	WGRID = [1540., 1541., 1542., 1543., 1544, 1545.]
   
	IDL>  LINTERP, WAVE, FLUX, WGRID, FGRID  

	FGRID will be a 6 element vector containing the values of FLUX 
	linearly interpolated onto the WGRID wavelength scale

 PROCEDURE: 
	Uses TABINV to calculate the effective index of the values
	in Xint in the table Xtab.  The resulting index is used
	with the intrinsic INTERPOLATE function to find the corresponding 
	Yint value in Ytab.  Unless the MISSING keyword is supplied, out
	of range Yint values are truncated to the nearest value of Ytab.

 NOTES:
	Users with IDL versions before V2.2.2 need to replace the call to
	INTERPOLATE with the commented lines at the end of the procedure

 PROCEDURES CALLED:
	TABINV, ZPARCHECK
 MODIFICATION HISTORY:
	Adapted from the IUE RDAF,  W. Landsman      October, 1988
	Modified to use the new INTERPOLATE function        June, 1992
	Modified to always return REAL*4             October, 1992
	Added MISSING keyword                        August, 1993

(See /usr/local/idl/lib/zastron/math/linterp.pro)


MINF_BRACKET

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	minF_bracket
 PURPOSE:
	Bracket a local minimum of a 1-D function with 3 points,
	thus ensuring that a minimum exists somewhere in the interval.
	This routine assumes that the function has a minimum somewhere....
	Routine can also be applied to a scalar function of many variables,
	for such case the local minimum in a specified direction is bracketed,
	This routine is called by minF_conj_grad, to bracket minimum in the 
	direction of the conjugate gradient of function of many variables
 CALLING EXAMPLE:
	xa=0  & xb=1					
	minF_bracket, xa,xb,xc, fa,fb,fc, FUNC_NAME="name"	;for 1-D func.
  or:
	minF_bracket, xa,xb,xc, fa,fb,fc, FUNC="name",     $
					  POINT=[0,1,1],   $
					  DIRECTION=[2,1,1]	;for 3-D func.
 INPUTS:
	xa = scalar, guess for point bracketing location of minimum.
	xb = scalar, second guess for point bracketing location of minimum.
 KEYWORDS:
	FUNC_NAME = function name (string)
		Calling mechanism should be:  F = func_name( px )
		where:
			px = scalar or vector of independent variables, input.
			F = scalar value of function at px.
	POINT_NDIM = when working with function of N variables,
		use this keyword to specify the starting point in N-dim space.
		Default = 0, which assumes function is 1-D.
	DIRECTION = when working with function of N variables,
		use this keyword to specify the direction in N-dim space
		along which to bracket the local minimum, (default=1 for 1-D).
		(xa,xb,xc) are then relative distances from POINT_NDIM.
 OUTPUTS:
	xa,xb,xc = scalars, 3 points which bracket location of minimum,
		that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
		When working with function of N variables
		(xa,xb,xc) are then relative distances from POINT_NDIM,
		in the direction specified by keyword DIRECTION,
		with scale factor given by magnitude of DIRECTION.
 OPTIONAL OUTPUT:
	fa,fb,fc = value of function at 3 points which bracket the minimum,
			again note that fb < fa and fb < fc if minimum exists.
 PROCEDURE:
	algorithm from Numerical Recipes (by Press, et al.), sec.10.1 (p.281).
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.

(See /usr/local/idl/lib/zastron/math/minf_bracket.pro)


MINF_CONJ_GRAD

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	minF_conj_grad
 PURPOSE:
	Find the local minimum of a scalar function of several variables using 
	the Conjugate Gradient method (Fletcher-Reeves-Polak-Ribiere algorithm).
	Function may be anything with computable partial derivatives.
	Each call to minF_conj_grad performs one iteration of algorithm,
	and returns an N-dim point closer to the local minimum of function.
 CALLING EXAMPLE:
	p_min = replicate( 1, N_dim )
	minF_conj_grad, p_min, f_min, conv_factor, FUNC_NAME="name",/INITIALIZE

	while (conv_factor GT 0) do begin
		minF_conj_grad, p_min, f_min, conv_factor, FUNC_NAME="name"
	endwhile
 INPUTS:
	p_min = vector of independent variables, location of minimum point
		obtained from previous call to minF_conj_grad, (or first guess).
 KEYWORDS:
	FUNC_NAME = function name (string)
		Calling mechanism should be:  F = func_name( px, gradient )
	  where:
		F = scalar value of function at px.
		px = vector of independent variables, input.
		gradient = vector of partial derivatives of the function
			with respect to independent variables, evaluated at px.
			This is an optional output parameter:
			gradient should not be calculated if parameter is not
			supplied in call (Unless you want to waste some time).
      /INIT must be specified on first call (whenever p_min is a guess),
			to initialize the iteration scheme of algorithm.
      /USE_DERIV causes the directional derivative of function to be used
			in the 1-D minimization part of algorithm
			(default is not to use directional derivative).
	TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).
      /QUADRATIC runs simpler version which works only for quadratic function.
 OUTPUTS:
	p_min = vector giving improved solution for location of minimum point.
	f_min = value of function at p_min.
	conv_factor = gives the current rate of convergence (change in value),
			iteration should be stopped when rate gets near zero.
 EXTERNAL CALLS:
	pro minF_bracket,  to find 3 points which bracket the minimum in 1-D.
	pro minF_parabolic,  to find minimum point in 1-D.
	pro minF_parabol_D,  to find minimum point in 1-D, using derivatives.
 COMMON BLOCKS:
	common minf_conj_grad, grad_conj, grad_save, gs_norm
	(to keep conjugate gradient, gradient and norm from previous iteration)
 PROCEDURE:
	Algorithm adapted from Numerical Recipes, sec.10.6 (p.305).
	Conjugate gradient is computed from gradient, which then gives
	the best direction (in N-dim space) in which to proceed to find
	the minimum point. The function is then minimized along
	this direction of conjugate gradient (a 1-D minimization).
	The algorithm is repeated starting at the new point by calling again.
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.

(See /usr/local/idl/lib/zastron/math/minf_conj_grad.pro)


MINF_PARABOLIC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	minF_parabolic
 PURPOSE:
	Find a local minimum of a 1-D function up to specified tolerance.
	This routine assumes that the function has a minimum nearby.
	(recommend first calling minF_bracket, xa,xb,xc, to bracket minimum).
	Routine can also be applied to a scalar function of many variables,
	for such case the local minimum in a specified direction is found,
	This routine is called by minF_conj_grad, to locate minimum in the 
	direction of the conjugate gradient of function of many variables.

 CALLING EXAMPLES:
	minF_parabolic, xa,xb,xc, xmin, fmin, FUNC_NAME="name"	;for 1-D func.
  or:
	minF_parabolic, xa,xb,xc, xmin, fmin, FUNC="name", $
					  POINT=[0,1,1],   $
					  DIRECTION=[2,1,1]	;for 3-D func.
 INPUTS:
	xa,xb,xc = scalars, 3 points which bracket location of minimum,
		that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
		When working with function of N variables
		(xa,xb,xc) are then relative distances from POINT_NDIM,
		in the direction specified by keyword DIRECTION,
		with scale factor given by magnitude of DIRECTION.
 KEYWORDS:
	FUNC_NAME = function name (string)
		Calling mechanism should be:  F = func_name( px )
		where:
			px = scalar or vector of independent variables, input.
			F = scalar value of function at px.

	POINT_NDIM = when working with function of N variables,
		use this keyword to specify the starting point in N-dim space.
		Default = 0, which assumes function is 1-D.
	DIRECTION = when working with function of N variables,
		use this keyword to specify the direction in N-dim space
		along which to bracket the local minimum, (default=1 for 1-D).
		(xa, xb, xc, x_min are then relative distances from POINT_NDIM)
	MAX_ITER = maximum allowed number iterations, default=100.
	TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).
 OUTPUTS:
	xmin = estimated location of minimum.
		When working with function of N variables,
		xmin is the relative distance from POINT_NDIM,
		in the direction specified by keyword DIRECTION,
		with scale factor given by magnitude of DIRECTION,
		so that min. Loc. Pmin = Point_Ndim + xmin * Direction.
	fmin = value of function at xmin (or Pmin).
 PROCEDURE:
	Brent's method to minimize a function by using parabolic interpolation,
	from Numerical Recipes (by Press, et al.), sec.10.2 (p.285).
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.

(See /usr/local/idl/lib/zastron/math/minf_parabolic.pro)


MINF_PARABOL_D

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	minF_parabol_D
 PURPOSE:
	Find a local minimum of a 1-D function up to specified tolerance,
	using the first derivative of function in the algorithm.
	This routine assumes that the function has a minimum nearby.
	(recommend first calling minF_bracket, xa,xb,xc, to bracket minimum).
	Routine can also be applied to a scalar function of many variables,
	for such case the local minimum in a specified direction is found,
	This routine is called by minF_conj_grad, to locate minimum in the 
	direction of the conjugate gradient of function of many variables.

 CALLING EXAMPLES:
	minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC_NAME="name"	;for 1-D func.
  or:
	minF_parabol_D, xa,xb,xc, xmin, fmin, FUNC="name", $
					  POINT=[0,1,1],   $
					  DIRECTION=[2,1,1]	;for 3-D func.
 INPUTS:
	xa,xb,xc = scalars, 3 points which bracket location of minimum,
		that is, f(xb) < f(xa) and f(xb) < f(xc), so minimum exists.
		When working with function of N variables
		(xa,xb,xc) are then relative distances from POINT_NDIM,
		in the direction specified by keyword DIRECTION,
		with scale factor given by magnitude of DIRECTION.
 KEYWORDS:
	FUNC_NAME = function name (string)
		Calling mechanism should be:  F = func_name( px, gradient )
		where:
			px = scalar or vector of independent variables, input.
			F = scalar value of function at px.
			gradient = derivative of function, a scalar if 1-D,
				a gradient vector if N-D,
				(should only be computed if arg. is present).

	POINT_NDIM = when working with function of N variables,
		use this keyword to specify the starting point in N-dim space.
		Default = 0, which assumes function is 1-D.
	DIRECTION = when working with function of N variables,
		use this keyword to specify the direction in N-dim space
		along which to bracket the local minimum, (default=1 for 1-D).
		(xa, xb, xc, x_min are then relative distances from POINT_NDIM)
	MAX_ITER = maximum allowed number iterations, default=100.
	TOLERANCE = desired accuracy of minimum location, default=sqrt(1.e-7).

 OUTPUTS:
	xmin = estimated location of minimum.
		When working with function of N variables,
		xmin is the relative distance from POINT_NDIM,
		in the direction specified by keyword DIRECTION,
		with scale factor given by magnitude of DIRECTION,
		so that min. Loc. Pmin = Point_Ndim + xmin * Direction.
	fmin = value of function at xmin (or Pmin).
 PROCEDURE:
	Brent's method to minimize a function by using parabolic interpolation
	and using first derivative of function,
	from Numerical Recipes (by Press, et al.), sec.10.3 (p.287),
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.

(See /usr/local/idl/lib/zastron/math/minf_parabol_d.pro)


PCA

[Previous Routine] [Next Routine] [List of Routines]
	PCA

 PURPOSE:
	Carry out a Principal Components Analysis (Karhunen-Loeve Transform)
	Results can be directed to the screen, a file, or output variables

 CALLING SEQUENCE:
	PCA, data, eigenval, eigenvect, percentages, proj_obj, proj_atr, 
		[MATRIX =, TEXTOUT = ,/COVARIANCE, /SSQ, /SILENT ]

 INPUT PARAMETERS:
	data -  2-d data matrix, data(i,j) contains the jth attribute value
		for the ith object in the sample.    If N_OBJ is the total
		number of objects (rows) in the sample, and N_ATTRIB is the 
		total number of attributes (columns) then data should be
		dimensioned N_OBJ x N_ATTRIB.

 OPTIONAL INPUT KEYWORD PARAMETERS:
	/COVARIANCE - if this keyword set, then the PCA will be carried out
		on the covariance matrix (rare), the default is to use the
		correlation matrix
	/SSQ - if this keyword is set, then the PCA will be carried out on
		on the sums-of-squares & cross-products matrix (rare)
	TEXTOUT - Controls print output device, defaults to !TEXTOUT

		textout=1	TERMINAL using /more option
		textout=2	TERMINAL without /more option
		textout=3	.prt
		textout=4	laser.tmp
		textout=5      user must open file
		textout = filename (default extension of .prt)

 OPTIONAL OUTPUT PARAMETERS:
	eigenval -  N_ATTRIB element vector containing the sorted eigenvalues
	eigenvect - N_ATRRIB x N_ATTRIB matrix containing the corresponding 
		eigenvectors
	percentages - N_ATTRIB element containing the cumulative percentage 
		variances associated with the principal components
	proj_obj - N_OBJ by N_ATTRIB matrix containing the projections of the 
		objects on the principal components
	proj_atr - N_ATTRIB by N_ATTRIB	matrix containing the projections of 
		the attributes on the principal components

 OPTIONAL OUTPUT PARAMETER
 	MATRIX   = analysed matrix, either the covariance matrix if /COVARIANCE
		is set, the "sum of squares and cross-products" matrix is
		/SSQ is set, or the (by default) correlation matrix.    Matrix
		will have dimensions N_ATTRIB x N_ATTRIB

 NOTES:
	This procedure performs Principal Components Analysis (Karhunen-Loeve
	Transform) according to the method described in "Multivariate Data 
	Analysis" by Murtagh & Heck [Reidel : Dordrecht 1987], pp. 33-48. 

	Keywords /COVARIANCE and /SSQ are mutually exclusive.

	The printout contains only (at most) the first seven principle 
	eigenvectors.    However, the output variables EIGENVECT contain 
	all the eigenvectors
	
       Different authors scale the covariance matrix in different ways.
	The eigenvalues output by PCA may have to be scaled by 1/N_OBJ or
	1/(N_OBJ-1) to agree with other calculations when /COVAR is set.

	PCA uses the non-standard system variables !TEXTOUT and !TEXTUNIT.
	These can be added to one's session using the procedure ASTROLIB.

	Users of V3.5.0 or later could change the calls to TQLI and TRED2
	to NR_TQLI and NR_TRED2

 EXAMPLE:
	Perform a PCA analysis on the covariance matrix of a data matrix, DATA,
	and write the results to a file

	IDL> PCA, data, /COVAR, t = 'pca.dat'

 	Perform a PCA analysis on the correlation matrix.   Suppress all 
	printing, and save the eigenvectors and eigenvalues in output varaibles

	IDL> PCA, data, eigenval, eigenvect, /SILENT

 PROCEDURES CALLED:
	Procedures TEXTOPEN, TEXTCLOSE

 COPYRIGHT:
	Copyright 1993, Hughes STX Corporation, Lanham MD 20706.

 REVISION HISTORY:
	Immanuel Freedman (after Murtagh F. and Heck A.).     December 1993
	Wayne Landsman, modified I/O              December 1993

(See /usr/local/idl/lib/zastron/math/pca.pro)


POIDEV

[Previous Routine] [Next Routine] [List of Routines]
 NAME
	POIDEV
 PURPOSE:
	Return an integer random deviate drawn from a Poisson distribution with
	a specified mean.    Adapted from procedure of the same name in 
	"Numerical Recipes" by Press et al. (1986), Section 7.3

 CALLING SEQUENCE:
	result = POIDEV( xm, [ SEED = ] )

 INPUTS:
	xm - integer scalar or vector, specifying mean of the Poisson 
		distribution

 OUTPUT:
	result - integer scalar or vector, same size as xm

 OPTIONAL KEYWORD INPUT-OUTPUT:
	SEED - floating point scalar, used as the seed for the random 
		distribution.     This keyword can be used to have POIDEV
		give identical results on consecutive runs.     

 EXAMPLE:
	(1) Add Poisson noise to an integral image array, im
		IDL> imnoise = POIDEV( im)

	(2) Verify the expected mean  and sigma for an input value of 81
		IDL> p = POIDEV( intarr(10000) + 81)   ;Test for 10,000 points
		IDL> print,avg(p),sigma(p)
	Average and sigma of the 10000 points should be close to 81 and 9

 METHOD: 
	For small values (< 20) independent exponential deviates are generated 
	until their sum exceeds the specfied mean, the number of events 
	required is returned as the Poisson deviate.   For large (> 20) values,
	uniform random variates are compared with a Lorentzian distribution 
	function.

 NOTES:
	Negative values in the input array will be returned as zeros.  

 PROCEDURES CALLED:
	GAMMLN - returns log of the Gamma function.    Users with IDL V3.0.0
		or later may wish to replace this call with the intrinisc 
		function LNGAMMA which was introduced in IDL 2.4.0.
             
 REVISION HISTORY:
	Version 1               Wayne Landsman        July  1992
	Added SEED keyword                            September 1992

(See /usr/local/idl/lib/zastron/math/poidev.pro)


POLINT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	POLINT
 PURPOSE:
	Interpolate a set of N points by fitting a polynomial of degree N-1
	Adapted from algorithm in Numerical Recipes, Press et al. (1986) 
	Section 3.1.

 CALLING SEQUENCE
	POLINT, xa, ya, x, y, [ dy ]
 INPUTS:
	XA - X Numeric vector, all values must be distinct.   The number
		of values in XA should rarely exceed 10 (i.e. a 9th order
		polynomial)
	YA - Y Numeric vector, same number of elements
	X - Numeric scalar specifying value to be interpolated 

 OUTPUT:
	Y - Scalar, interpolated value in (XA,YA) corresponding to X

 OPTIONAL OUTPUT
	DY - Error estimate on Y, scalar

 EXAMPLE:
	Find sin(2.5) by polynomial interpolation on sin(indgen(10))

		IDL> xa = indgen(10)
		IDL> ya = sin( xa )
		IDL> polint, xa, ya, 2.5, y ,dy
             
	The above method gives  y = .5988 & dy = 3.1e-4  a close
	approximation to the actual sin(2.5) = .5985

 METHOD:
	Uses Neville's algorithm to iteratively build up the correct
	polynomial, with each iteration containing one higher order.

 REVISION HISTORY:
	Written W. Landsman                 January, 1992

(See /usr/local/idl/lib/zastron/math/polint.pro)


POLY_SMOOTH

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

 PURPOSE:
	Reduce noise in 1-D data (e.g. time-series, spectrum) 
	but retain dynamic range of variations in the data by
	applying a least squares smoothing polynomial filter,

	also called the Savitzky-Golay smoothing filter, cf. Numerical
	Recipes (Press et al. 1992, Sec.14.8)

	The low-pass filter coefficients are computed by effectively
	least-squares fitting a polynomial in moving window,
	centered on each data point, so the new value will be the
	zero-th coefficient of the polynomial. Approximate first derivates
	of the data can be computed by using first degree coefficient of
	each polynomial, and so on. The filter coefficients for a specified
	polynomial degree and window width are computed independent of any
	data, and stored in a common block. The filter is then convolved
	with the data array to result in smoothed data with reduced noise,
	but retaining higher order variations (better than SMOOTH).

 CALLING SEQUENCE:

	spectrum = poly_smooth( data, [ width, DEGREE = , NLEFT = , NRIGHT = 
					DERIV_ORDER = ,COEFF = ]

 INPUTS:
	data = 1-D array, such as a spectrum or time-series.

	width = total number of data points to use in filter convolution,
		(default = 5, using 2 past and 2 future data points),
		must be larger than DEGREE of polynomials, and a guideline is 
		make WIDTH between 1 and 2 times the FWHM of desired features.

 KEYWORDS:

	DEGREE = degree of polynomials to use in designing the filter
		via least squares fits, (default DEGREE = 2), and
		the higher degrees will preserve sharper features.

	NLEFT = # of past data points to use in filter convolution,
		excluding current point, overrides width parameter,
		so that width = NLEFT + NRIGHT + 1.  (default = NRIGHT)

	NRIGHT = # of future data points to use (default = NLEFT).

	DERIV_ORDER = order of derivative desired (default = 0, no derivative).

	COEFFICIENTS = optional output of the filter coefficients applied,
		but they are all stored in common block for reuse, anyway.
 RESULTS:
	Function returns the data convolved with polynomial filter coefs.

 EXAMPLE:

	Given a wavelength - flux spectrum (w,f), apply a 31 point quadratic
	smoothing filter and plot

	IDL> plot, w, poly_smooth(f,31)	
 COMMON BLOCKS:
	common poly_smooth, degc, nlc, nrc, coefs, ordermax

 PROCEDURE:
	As described in Numerical Recipies, 2nd edition sec.14.8, 
	Savitsky-Golay filter.
	Matrix of normal eqs. is formed by starting with small terms
	and then adding progressively larger terms (powers).
	The filter coefficients of up to derivative ordermax are stored
	in common, until the specifications change, then recompute coefficients.
	Coefficients are stored in convolution order, zero lag in the middle.
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1993.

(See /usr/local/idl/lib/zastron/math/poly_smooth.pro)


PROB_KS

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	PROB_KS
 PURPOSE:
	Returns the significance level of an observed value of the 
	Kolmogorov-Smirnov statistic D for an effective number of data points
	N_eff.   Called by KSONE and KSTWO

 CALLING SEQUENCE:
	prob_ks, D, N_eff, probks

 INPUT PARAMATERS:
	D -  Kolmogorov statistic, floating scalar, always non-negative
	N_eff - Effective number of data points, scalar.   For a 2 sided test 
		this is given by (N1*N2)/(N1+N2) where N1 and N2 are the number 
		of points in each data set.

 OUTPUT PARAMETERS:
	probks - floating scalar between 0 and 1 giving the significance level of
		the K-S statistic.   Small values of PROB suggest that the 
		distribution being tested are not the same

 REVISION HISTORY:
	Written     W. Landsman                August, 1992

(See /usr/local/idl/lib/zastron/math/prob_ks.pro)


QSIMP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	QSIMP
 PURPOSE:
	Integrate a function to specified accuracy using the extended 
	trapezoidal rule.   Adapted from algorithm in Numerical Recipes, 
	by Press et al. (1986) p. 111.     This procedure became mostly
	obsolete in IDL V3.5 with the introduction of the intrinsic 
	procedure NR_QSIMP.

 CALLING SEQUENCE:
	QSIMP, func, A, B, S, [ EPS = , MAX_ITER = ]

 INPUTS:
	func - scalar string giving name of function of one variable to 
		be integrated
	A,B  - numeric scalars giving the lower and upper bound of the 
		integration

 OUTPUTS:
	S - Scalar giving the approximation to the integral of the specified
		function between A and B.

 OPTIONAL KEYWORD PARAMETERS:
	EPS - scalar specify the fractional accuracy before ending the 
		iteration.  Default = 1E-6
	MAX_ITER - Integer specifying the total number iterations at which 
		QSIMP will terminate even if the specified accuracy has not yet
		been met.   The maximum number of function evaluations will be
		2^(MAX_ITER).    Default value is MAX_ITER = 20

 NOTES:
	The function QTRAP is robust way of doing integrals that are not very 
	smooth.  However, if the function has a continuous 3rd derivative then 
	QSIMP will likely be more efficient at performing the integral.

 EXAMPLE:
	Compute the integral of sin(x) from 0 to !PI/3.
    
	IDL> QSIMP, 'sin', 0, !PI/3, S   & print, S
   
	The value obtained should be cos(!PI/3) = 0.5

 PROCEDURES CALLED:
	TRAPZD, ZPARCHECK

 REVISION HISTORY:
	W. Landsman         ST Systems Co.         August, 1991

(See /usr/local/idl/lib/zastron/math/qsimp.pro)


QTRAP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	QTRAP
 PURPOSE:
	Integrate a function to specified accuracy using the extended trapezoidal
	rule.   Adapted from Numerical Recipes, p. 111. 

 CALLING SEQUENCE:
	QTRAP, func, A, B, S, [EPS = , MAX_ITER = ]

 INPUTS:
	func - scalar string giving name of function of one variable to 
		be integrated
	A,B  - numeric scalars giving the lower and upper bound of the 
		integration

 OUTPUTS:
	S - Scalar giving the approximation to the integral of the specified
		function between A and B.

 OPTIONAL KEYWORD PARAMETERS:
	EPS - scalar specify the fractional accuracy before ending the 
		iteration.    Default = 1E-6
	MAX_ITER - Integer specifying the total number iterations at which 
		QTRAP will terminate even if the specified accuracy has not yet
		 been met.    The maximum number of function evaluations will 
		be 2^(MAX_ITER).   Default value is MAX_ITER = 20

 NOTES:
	QTRAP is robust way of doing integrals that are not very smooth.  If the
	function has a continuous 3rd derivative then the function QSIMP will 
	   likely be more efficient at performing the integral.
 EXAMPLE:
	Compute the integral of sin(x) from 0 to !PI/3.
    
	IDL> QTRAP, 'sin', 0, !PI/3, S   & print,S
   
	The value obtained should be cos(!PI/3) = 0.5

 PROCEDURES CALLED:
	TRAPZD, ZPARCHECK
 REVISION HISTORY:
	W. Landsman         ST Systems Co.         August, 1991

(See /usr/local/idl/lib/zastron/math/qtrap.pro)


QUADTERP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	QUADTERP     
 PURPOSE:
	Quadratically interpolate (3 point Lagrangian) a function Y = f(X)
	at specified grid points.  Use LINTERP for linear interpolation

 CALLING SEQUENCE:
	QUADTERP, Xtab, Ytab, Xint, Yint, [ MISSING = ]

 INPUT: 
	Xtab - Vector (X TABle) containing the current independent variable 
		Must be either monotonic increasing or decreasing
	Ytab - Vector (Y TABle) containing the dependent variable defined
		at each of the points of XTAB.
	Xint - Scalar or vector giving the values of X for which interpolated 
		Y values are sought

 OUTPUT: 
	Yint - Interpolated value(s) of Y, same number of points as Xint

 OPTIONAL INPUT KEYWORD:
	MISSING - Scalar specifying Yint value(s) to be assigned, when Xint
		value(s) are outside of the range of Xtab.     Default is to
		truncate the out of range Yint value(s) to the nearest value 
		of Ytab.   See the help for the INTERPOLATE function.
 METHOD:
	3-point Lagrangian interpolation.  The average of the two quadratics 
	derived from the four nearest  points is returned in YTAB
	TABINV is used to locate center point of the interpolation.

 RESTRICTIONS:
	Unless MISSING keyword is set, points outside the range of Xtab in 
	which two valid quadratics can be computed are returned at the value 
	of the next to last point of Ytab (i.e. Ytab(1) and Ytab(NPTS-2) ).

 EXAMPLE:
	A spectrum has been defined using a wavelength vector WAVE and a
	flux vector FLUX.  Interpolate onto a new wavelength grid, e.g. 

	IDL> wgrid = [1540.,1541.,1542.,1543.,1544.,1545.]
	IDL> quadterp, wave, flux, wgrid, fgrid 
     
	FGRID will be a 5 element vector containing the quadratically
	interpolated values of FLUX at the wavelengths given in WGRID.

  EXTERNAL ROUTINES:
	TABINV, ZPARCHECK, DATATYPE
  REVISION HISTORY:
	31 October 1986 by B. Boothman, adapted from the IUE RDAF
	12 December 1988 J. Murthy, corrected error in Xint
	September 1992, W. Landsman, fixed problem with double precision
       August 1993, W. Landsman, added MISSING keyword

(See /usr/local/idl/lib/zastron/math/quadterp.pro)


SIGMA

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	SIGMA
 PURPOSE:
	Calculate the standard deviation value of an array, or calculate the 
	standard deviation over one dimension of an array as a function of all
	the other dimensions.
 CALLING SEQUENCE:
	RESULT = SIGMA(ARRAY)
	RESULT = SIGMA(ARRAY,N_PAR)
	RESULT = SIGMA(ARRAY,N_PAR,DIM)
 INPUTS:
	ARRAY = Input array.  May be any type except string.
 OPTIONAL INPUT PARAMETERS:
	N_PAR = Number of parameters.  Default value is zero.  The number of
		degrees of freedom is N_FREE = N_POINTS - N_PAR, where N_POINTS
		is either N_ELEMENTS(ARRAY) or the size of dimension DIM.  The
		value of SIGMA varies as 1 / SQRT(N_FREE).  If N_PAR is
		negative, then the absolute value of N_PAR is used, and an
		additional factor of 1 / SQRT(N_POINTS) is included, yielding
		the reduced sigma.
	DIM   = Optional dimension to do standard deviation over.
 OUTPUTS:
	The standard deviation value of the array when called with one
	parameter.

	If DIM is passed, then the result is an array with all the dimensions
	of the input array except for the dimension specified, each element of
	which is the standard deviation of the corresponding vector in the
	input array.

	For example, if A is an array with dimensions of (3,4,5), then the
	command B = SIGMA(A,N,1) is equivalent to

			B = FLTARR(3,5)
			FOR J = 0,4 DO BEGIN
			    FOR I = 0,2 DO BEGIN
			        B(I,J) = SIGMA( A(I,*,J) , N )
			    ENDFOR
			ENDFOR

 RESTRICTIONS:
	Dimension specified must be valid for the array passed; otherwise the
	input array is returned as the output array.
 PROCEDURE:
	Uses the function AVG.
	When DIM is passed, then the function SUM is used.
 MODIFICATION HISTORY:
	William Thompson	Applied Research Corporation
	July, 1986		8201 Corporate Drive
				Landover, MD  20785

(See /usr/local/idl/lib/zastron/math/sigma.pro)


SIXLIN

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	SIXLIN
 PURPOSE:
	Compute linear regression coefficients by six different methods.
	Adapted from the FORTRAN program (Rev. 1.1)  supplied by Isobe, 
	Feigelson, Akritas, and Babu Ap. J. Vol. 364, p. 104 (1990).   
	Suggested when there is no understanding about the nature of the 
	scatter about a linear relation, and NOT when the errors in the 
	variable are calculable.

 CALLING SEQUENCE:
		SIXLIN, xx, yy, a, siga, b, sigb   

 INPUTS:
	XX - vector of X values
	YY - vector of Y values, same number of elements as XX

 OUTPUTS:
	A - Vector of 6 Y intercept coefficients
	SIGA - Vector of standard deviations of 6 Y intercepts
	B - Vector of 6 slope coefficients
	SIGB - Vector of standard deviations of slope coefficients

	The output variables are computed using linear regression for each of 
	the following 6 cases:
		(0) Ordinary Least Squares (OLS) Y vs. X
		(1) Ordinary Least Squares  X vs. Y
		(2) Ordinary Least Squares Bisector
		(3) Orthogonal Reduced Major Axis
		(4) Reduced Major-Axis 
		(5) Mean ordinary Least Squares

 NOTES:
	Isobe et al. make the following recommendations

	(1) If the different linear regression methods yield similar results
 		then quoting OLS(Y|X) is probably the most familiar.

	(2) If the linear relation is to be used to predict Y vs. X then
		OLS(Y|X) should be used.   

	(3) If the goal is to determine the functional relationship between
		X and Y then the OLS bisector is recommended.

 REVISION HISTORY:
	Written   Wayne Landsman          February, 1991         
	Corrected sigma calculations      February, 1992

(See /usr/local/idl/lib/zastron/math/sixlin.pro)


SPLINE_SMOOTH

[Previous Routine] [Next Routine] [List of Routines]
	SPLINE_SMOOTH

 PURPOSE:

	Construct cubic smoothing spline (or give regression solution) to given
	data with minimum "roughness" (measured by the energy in the second 
	derivatives) whilst restricting the weighted mean square distance
	of the approximation from the data.  The results may be written to
	the screen or a file or both and are optionally returned in the 
	parameters.  The results may be optionally displayed graphically.

 CATEGORY:
	Math and Statistics [PRO.MATH].

 CALLING SEQUENCE:
	SPLINE_SMOOTH,X,Y,Yerr,distance,[coefficients,smoothness,xplot,yplot XTITLE=xtitle,YTITLE=ytitle, INTERP=interp, TEXTOUT=,/SILENT,/PLOT,/ERRBAR]

 INPUT PARAMETERS:
	X - N_POINT element vector containing the data abcissae
	Y - N_POINT element vector containing the data ordinates
	Yerr -     estimated uncertainty in ordinates ( positive)
	distance - upper bound on the weighted mean square distance
		   of the approximation from the data


 OPTIONAL INPUT PARAMETERS

      xplot -    vector of spline evaluation abcissae

 OPTIONAL INPUT KEYWORD PARAMETERS:
	TEXTOUT - Controls print output device, defaults to !TEXTOUT
	
		textout=1	TERMINAL using /more option
		textout=2	TERMINAL without /more option
		textout=3	.prt
		textout=4	laser.tmp
		textout=5	user must open file
		textout = filename (default extension of .prt)

 OPTIONAL OUTPUT PARAMETERS:

	coefficients - N_POINT x 4 element array containing the sequence of
                      spline coefficients including the smoothed ordinates.

	smoothness  - N_POINT element vector containing the energy in second 
		      derivatives of approximated function.
       yplot       - vector of evaluated spline ordinates.

 OPTIONAL OUTPUT KEYWORD PARAMETERS
	/SILENT	    - suppress all printing.
	/PLOT	    - display smooth curve, data ordinates and error bars
       /ERRBAR     - display error bars
	XTITLE      - optional title for X-axis
	YTITLE      - optional title for Y-axis
       INTERP      - optionally returned interpolated smooth spline
 NOTES:
	This procedure constructs a smoothing spline according to the method
	described in "Fundamentals of Image Processing" by A. Jain  [Prentice-
	Hall : New Jersey 1989].
	If the distance parameter is sufficiently large a linear regression
	is performed, otherwise a cubic smoothing spline is constructed.

	This procedure assumes regular sampling and independent identically 
	distributed normal errors without missing data.  The data are sorted.
	
	SPLINE_SMOOTH uses the non-standard system variables !TEXTOUT and
       !TEXTUNIT.
	These can be added to one's session using the procedure ASTROLIB.

 COMMON BLOCKS:
	None.
 EXAMPLE:
	Obtain coefficients of a univariate smoothing spline fitted to data
       X,Y assuming normally distributed errors Yerr and write the results to
       a file.

	IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, smoothness,
	     t='spline.dat'

	Fit a smoothing spline to observational data.  Suppress all printing 
	and save the smoothed ordinates in output variables. Display results.

	IDL> SPLINE_SMOOTH, X, Y, Yerr, distance, coefficients, /SILENT, /PLOT
	
 PROCEDURES CALLED:
	Procedures TEXTOPEN, TEXTCLOSE, PLOT, PLOTERR

 COPYRIGHT:
	Copyright 1993, Hughes STX Corporation, Lanham MD 20706.
 AUTHOR:
	Immanuel Freedman (after A. Jain).	December, 1993
 REVISIONS
       January 12, 1994    I. Freedman (HSTX)  Adjusted formats
       March   14, 1994    I. Freedman (HSTX)  Improved convergence
       March   15, 1994    I. Freedman (HSTX)  User-specified interpolates

(See /usr/local/idl/lib/zastron/math/spline_smooth.pro)


TABINV

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	TABINV     
 PURPOSE:  
	To find the effective index of a function value in
	an ordered vector.

 CALLING SEQUENCE:
	TABINV, XARR, X, IEFF
 INPUTS:
	XARR - the vector array to be searched, must be monotonic
		increasing or decreasing
	X    - the function value(s) whose effective
		index is sought (scalar or vector)

 OUTPUT:
	IEFF - the effective index or indices of X in XARR
		real or double precision, same # of elements as X

 RESTRICTIONS:
	TABINV will abort if XARR is not monotonic.  (Equality of 
	neighboring values in XARR is allowed but results may not be
	unique.)  This requirement may mean that input vectors with padded
	zeroes could cause routine to abort.

 PROCEDURE:
	A binary search is used to find the values XARR(I)
	and XARR(I+1) where XARR(I) < X < XARR(I+1).
	IEFF is then computed using linear interpolation 
	between I and I+1.
		IEFF = I + (X-XARR(I)) / (XARR(I+1)-XARR(I))
	Let N = number of elements in XARR
		if x < XARR(0) then IEFF is set to 0
		if x > XARR(N-1) then IEFF is set to N-1

 EXAMPLE:
	Set all flux values of a spectrum (WAVE vs FLUX) to zero
	for wavelengths less than 1150 Angstroms.
         
	IDL> tabinv, wave, 1150.0, I
	IDL> flux( 0:fix(I) ) = 0.                         

 FUNCTIONS CALLED:
	ISARRAY
 REVISION HISTORY:
	Adapted from the IUE RDAF                     January, 1988         
	More elegant code  W. Landsman                August, 1989
	Mod to work on 2 element decreasing vector    August, 1992

(See /usr/local/idl/lib/zastron/math/tabinv.pro)


TRAPZD

[Previous Routine] [Next Routine] [List of Routines]
 NAME
	TRAPZD
 PURPOSE:
	Compute the nth stage of refinement of an extended trapezoidal rule.
	This procedure is called by QSIMP and QTRAP.   Algorithm from Numerical
	Recipes, Section 4.2.   TRAPZD is meant to be called iteratively from
	a higher level procedure.

 CALLING SEQUENCE:
	TRAPZD, func, A, B, S, step

 INPUTS:
	func - scalar string giving name of function to be integrated.   This
		must be a function of one variable.
	A,B -  scalars giving the limits of the integration

 INPUT-OUTPUT:
	S -    scalar giving the total sum from the previous interations on 
		input and the refined sum after the current iteration on output.

	step - LONG scalar giving the number of points at which to compute the
		function for the current iteration.   If step is not defined on
		input, then S is intialized using the average of the endpoints
		of limits of integration.

 NOTES:
	TRAPZD will check for math errors when computing the function at the
	endpoints, but not on subsequent iterations.

 REVISION HISTORY:
	Written         W. Landsman                 August, 1991

(See /usr/local/idl/lib/zastron/math/trapzd.pro)


TSUM

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
	TSUM
 PURPOSE:
	Trapezoidal summation of the area under a curve.   This procedure
	was formerly know as INTEG

 CALLING SEQUENCE:
       Result = TSUM(y)
              or
       Result = TSUM( x, y, [ imin, imax ] )  
 INPUTS:
	x = array containing independent variable.  If omitted, then
		x is assumed to contain the index of the y variable.
		x = indgen( N_elements(y) ).
	y = array containing dependent variable y = f(x)

 OPTIONAL INPUTS:
	imin = index of x array at which to begin the integration, integer
		scalar.  If omitted, then summation starts at x(0).
	imax = index of x value at which to end the integration, integer 
		scalar.  If omitted then the integration ends at x(npts).

 OUTPUTS:
	result = area under the curve y=f(x) between x(imin) and x(imax).

 PROCEDURE:
	The area is determined of indivdual trapezoids defined by x(i),
	x(i+1), y(i) and y(i+1).

 MODIFICATION HISTORY:
	Written, W.B. Landsman, STI Corp. May 1986
	Modified so X is not altered in a one parameter call Jan 1990

(See /usr/local/idl/lib/zastron/math/tsum.pro)


ZBRENT

[Previous Routine] [List of Routines]
 NAME:
	ZBRENT
 PURPOSE:
	Find the zero of a 1-D function up to specified tolerance.
	This routine assumes that the function is known to have a zero.

 CALLING:
	x_zero = ZBRENT( x1, x2, FUNC_NAME="name" )

 INPUTS:
	x1, x2 = scalars, 2 points which bracket location of function zero,
						that is, F(x1) < 0 < F(x2).
	Note: computations are performed with
	same precision (single/double) as the inputs and user supplied function.

 REQUIRED INPUT KEYWORD:
	FUNC_NAME = function name (string)
		Calling mechanism should be:  F = func_name( px )
		where:	px = scalar independent variable, input.
			F = scalar value of function at px,
			    should be same precision (single/double) as input.

 OPTIONAL INPUT KEYWORD:
	MAX_ITER = maximum allowed number iterations, default=100.
	TOLERANCE = desired accuracy of minimum location, default = 1.e-3.

 OUTPUTS:
	Returns the location of zero, with accuracy of specified tolerance.

 PROCEDURE:
	Brent's method to find zero of a function by using bracketing,
	bisection, and inverse quadratic interpolation,
	from Numerical Recipes (by Press, et al.), sec.9.3 (p.251).

 EXAMPLE:
       Find the root of the COSINE function between 1. and 2.  radians

        IDL> print, zbrent( 1, 2, FUNC = 'COS')

       and the result will be !PI/2 within the specified tolerance
 MODIFICATION HISTORY:
	Written, Frank Varosi NASA/GSFC 1992.
	FV.1994, mod to check for single/double prec. and set zeps accordingly.

(See /usr/local/idl/lib/zastron/math/zbrent.pro)