Numerical routines

program names, file names
variable names
prompts, commands, program code


Contents:

Overview

This page gives documentation for several numerical routines which I have found useful for a variety of applications. Most of them are only applicable to 1-D arrays of numbers, but the first (lims) is also useful for multi-dimensional input data.


lims

The lims function returns a 5-element vector giving certain statistical information about the input array (which can be a 1-D vector, a 2-D image, or even a higher-dimensional array):
IDLprompt> result_arr = lims( in_arr )
The fourth and fifth elements of this vector (   result_arr(3) and result_arr(4)   ) give the mean and standard deviation, respectively, of in_arr.

The original purpose of lims was to determine viewing parameters for images. As such, three optional parameters may be given to lims:
IDLprompt> result_arr = lims( in_arr , min , max , sig=3 )
In this example, min and max will, after runtime, contain the values which are three standard deviations below and above the mean of in_arr (the default value for sig is 5). These three numbers---min, max, and the value of sig used by lims) are, in order, the first three elements of result_arr.

lims may also be used on impacks, but to do so you MUST set either the im or the wt keyword, in the same context as used in filter:
IDLprompt> im_res_arr = lims( impack , /im ) IDLprompt> wt_res_arr = lims( impack , /wt )


fwhm

The fwhm function is used to compute the FWHM (full-width at half-maximum) of a 1-D array (it also computes the position of the centroid). Let's say you have an array y which contains the ordinate values of a singly-peaked function (we're assuming the abscissa values to be evenly spaced):

[Profile of a
singly-peaked function]
To find the FWHM (and the centriod), fwhm needs (1) to identify the peak of interest and (2) to determine the baseline of the function. In this example, (1) is easy: it's the only interesting peak in the function. However, if the function had more than one peak, such that the peak of interest were NOT the tallest feature, it would be necessary to manually identify the peak for fwhm. This is done by setting the pickpeak keyword. Otherwise, fwhm just uses the tallest peak in the input function (as for how to identify the peak, should that be necessary, set the help keyword and follow the directions).

Determining the baseline sometimes takes a bit of judgement. The baseline is an ordinate value which gives some indication of the zero-point of the function. Essentially, you need to pick a baseline value such that if you drew a horizontal line across the function at the height of the baseline value, that line would cross the peak on both sides near the base of the peak. Set the help keyword for directions on how to specify the baseline interactively. Or, if you have a good idea of the baseline value beforehand (a baseline of zero [0] would be appropriate for the above example), you can specify this with the zero keyword.

I recommend that you run fwhm with the help keyword set until you get used to how it works (when you do so, you'll likely recognize fwhm as being used by cuts). To find the FWHM of our function y, do one of the following:
IDLprompt> fwhm_val = fwhm( y , /help )
Or, to pick out the peak by hand (not necessary here, but illustrative),
IDLprompt> fwhm_val = fwhm( y , /pickpeak , /help )
Or, to specify a baseline value of zero,
IDLprompt> fwhm_val = fwhm( y , zero = 0 , /help )
In each case, the result fwhm_val will be a scalar giving the FWHM of the function. A plot will be generated, which will show the centroid position, the half-power value, and the positions of the half-power points:

[Plot generated by
fwhm.pro]

In some cases, you'll need to find the FWHM of a function whose abscissa values are not evenly spaced. If so, you'll need to provide fwhm with an array of abscissa values (let's call it x). You'll give it to fwhm as the second positional parameter, and then proceed as before:
IDLprompt> fwhm_val = fwhm( y , x , /help )

fwhm has several other features (reading input from a text file, or divering the output plot to a Postscript file), descriptions for which I refer you to the program's documentation (using dl).


fitgauss

The fitgauss procedure is used to fit Gaussian profiles to input functions (such as the profile discussed above, in the fwhm section). The following command is a typical call to fitgauss:
IDLprompt> fitgauss, ordinate, abscissa, ord_fit, fitparams
where fitgauss and ordinate are 1-D arrays giving the ordinate and abscissa values for which you want a Gaussian fit. After runtime, ord_fit will contain the ordinate values of the fit, and fitparams will be a vector containing the parameters of the fit (see the program documentation for details about fitparams).


findroot

findroot, as its name suggests, is used to find zero-crossings in a function. Its input is much like that of the fwhm function, in that the first positional parameter is the ordinate, and an optional second positional parameter is the abscissa (if no abscissa is given, unit spacing of the points is assumed). findroot returns the position of the first root found in the input function, this being a fractional array index or an interpolated abscissa value, depending on whether or not an abscissa was provided. The default is to begin looking for a root at the low-array-index end of the function and to proceed to higher array indices. However, you can specify a starting point with the keyword start (this is an array index, even if an abscissa is provided), and you can also move backwards through the input function (ie., to lower array indices) with the backwards keyword. The following command would find the first root of the function ordinate to the left of the abscissa value with array index 12:
IDLprompt> result = findroot( ordinate , abscissa , start=12 , /backwards )
The interpolated abscissa value of the root would be assigned to result. If no root is found, result will be given a value of -1.  .


polyfit

The polyfit procedure is used to fit polynomials to data points. polyfit uses the least-squares method of curve-fitting. Let's assume that we have two 1-D arrays x and y which hold the abscissa and ordinate values for a set of data points. To fit an Nth-order polynomial to this data, use the following command:
IDLprompt> polyfit, x, y, N, coefs=coefs, LS=LS
After runtime coefs will be an (N+1)-element vector containing the parameters of the fit (coefs(0) will be the constant term, coefs(1) will be the coefficient of the linear term, ..., and coefs(N) will be the coefficient of the Nth-order term), and LS (remember that IDL is case-insensitive) will have the ``least-squares'' value (this last quantity is computed just as a chi-squares value would be, but with unit uncertainty for each ordered pair).


Carl Welch