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

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
```

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

[Previous Routine] [Next Routine] [List of Routines]
``` NAME:
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 )

while (conv_factor GT 0) do begin
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).
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:
PROCEDURE:
Algorithm adapted from Numerical Recipes, sec.10.6 (p.305).
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.
```

## 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,
(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 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
```

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

[Previous Routine] [Next Routine] [List of Routines]
``` NAME:
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
```

## 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 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,
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)