path: root/noao/onedspec/doc/scombine.hlp

.help scombine Sep97 noao.onedspec
.ih
NAME
scombine -- Combine spectra
.ih
USAGE
scombine input output
.ih
PARAMETERS
.ls input
List of input images containing spectra to be combined.  The spectra
in the images to be combined are selected with the \fIapertures\fR and
\fIgroup\fR parameters.  Only the primary spectrum is combined and
the associated band spectra are ignored.
.le
.ls output
List of output images to be created containing the combined spectra.
If the grouping option is "all"
or "apertures" then only one output image will be created.  In the
first case the image will contain only one spectrum and in the latter case
there will be a spectrum for each selected aperture.
If the grouping option is "images" then there will be one
output spectrum per input spectrum.
.le
.ls noutput = ""
List of output images to be created containing the number of spectra combined.
The number of images required is the same as the \fIoutput\fR list.
Any or all image names may be given as a null string, i.e. "", in which
case no output image is created.
.le
.ls logfile = "STDOUT"
File name for recording log information about the combining operation.
The file name "STDOUT" is used to write the information to the terminal.
If the null string is specified then no log information is printed or
recorded.
.le

.ls apertures = ""
List of apertures to be selected for combining.  If none is specified
then all apertures are selected.  The syntax is a blank or comma separated
list of aperture numbers or aperture ranges separated by a hyphen.
.le
.ls group = "apertures" (all|images|apertures)
Option for grouping input spectra for combining (after selection by aperture)
from one or more input images.  The options are:
.ls "all"
Combine all spectra from all images in the input list into a single output
spectrum.
.le
.ls "images"
Combine all spectra in each input image into a single spectrum in
separate output images.
.le
.ls "apertures"
Combine all spectra of the same aperture from all input images and put it
into a single output image with the other selected apertures.
.le
.le
.ls combine = "average" (average|median|sum)
Option for combining pixels at the same dispersion coordinate.  after any
rejection operation.  The options are to compute the  "average", "median",
or "sum" of the pixels.  The first two are applied after any pixel
rejection.  The sum option ignores the rejection and scaling parameters and
no rejection is performed.  In other words, the "sum" option is simply the
direct summation of the pixels.  The median uses the average of the two
central values when the number of pixels is even.
.le
.ls reject = "none" (none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip)
Type of rejection operation performed on the pixels which overlap at each
dispersion coordinate.  The algorithms are discussed in the
DESCRIPTION section.  The rejection choices are:

.nf
      none - No rejection
    minmax - Reject the nlow and nhigh pixels
   sigclip - Reject pixels using a sigma clipping algorithm
 avsigclip - Reject pixels using an averaged sigma clipping algorithm
   ccdclip - Reject pixels using CCD noise parameters
  crreject - Reject only positive pixels using CCD noise parameters
     pclip - Reject pixels using sigma based on percentiles
.fi

.le

.ls first = no
Use the first input spectrum of each set to be combined to define the
dispersion coordinates for combining and output?  If yes then all other
spectra to be combined will be interpolated to the dispersion of this
reference spectrum and that dispersion defines the dispersion of the
output spectrum.  If no, then all the spectra are interpolated to a linear
dispersion as determined by the following parameters.  The interpolation
type is set by the package parameter \fIinterp\fR.
.le
.ls w1 = INDEF, w2=INDEF, dw = INDEF, nw = INDEF, log = no
The output linear or log linear wavelength scale if the dispersion of the
first spectrum is not used.  INDEF values are filled in from the maximum
wavelength range and minimum dispersion of the spectra to be combined.  The
parameters are aways specified in linear wavelength even when the log
parameter is set to produce constant pixel increments in the log of the
wavelength.  The dispersion is interpreted in that case as the difference
in the log of the endpoints divided by the number of pixel increments.
.le

.ls scale = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>)
Multiplicative image scaling to be applied.  The choices are none,
multiply by the reciprocal of the mode , median, or mean of the specified
statistics section, scale by the exposure time in the image header, multiply
by the values in a specified file, or multiply by a specified image header
keyword.  When specified in a file the scales must be one per line in the
order of the input spectra.
.le
.ls zero = "none" (none|mode|median|mean|@<file>|!<keyword>)
Additive zero level image shifts to be applied.  The choices are none,
add the negative of the mode, median, or mean of the specified statistics
section, add the values given in a file, or add values given by an
image header keyword.  When specified in a file the zero values must be one
per line in the order of the input spectra. File or keyword zero offset
values do not allow a correction to the weights.
.le
.ls weight = "none" (none|mode|median|mean|exposure|@<file>|!<keyword>)
Weights to be applied during the final averaging.  The choices are none,
the mode, median, or mean of the specified statistics section, the exposure
time, values given in a file, or values given by an image header keyword.
When specified in a file the weights must be one per line in the order of
the input spectra.
.le
.ls sample = ""
Wavelength sample regions to use in computing spectrum statistics for
scaling and weighting.  If no sample regions are given then the entire
input spectrum is used.  The syntax is colon separated wavelengths
or a file containing colon separated wavelengths preceded by the
@ character; i.e. @<file>.
.le

.ce
Algorithm Parameters
.ls lthreshold = INDEF, hthreshold = INDEF
Low and high thresholds to be applied to the input pixels.  This is done
before any scaling, rejection, and combining.  If INDEF the thresholds
are not used.
.le
.ls nlow = 1,  nhigh = 1 (minmax)
The number of low and high pixels to be rejected by the "minmax" algorithm.
These numbers are converted to fractions of the total number of input spectra
so that if no rejections have taken place the specified number of pixels
are rejected while if pixels have been rejected by thresholding
or nonoverlap, then the fraction of the remaining pixels, truncated
to an integer, is used.
.le
.ls nkeep = 1
The minimum number of pixels to retain or the maximum number to reject
when using the clipping algorithms (ccdclip, crreject, sigclip,
avsigclip, or pclip).  When given as a positive value this is the minimum
number to keep.  When given as a negative value the absolute value is
the maximum number to reject.  This is actually converted to a number
to keep by adding it to the number of images.
.le
.ls mclip = yes (ccdclip, crreject, sigclip, avsigcliip)
Use the median as the estimate for the true intensity rather than the
average with high and low values excluded in the "ccdclip", "crreject",
"sigclip", and "avsigclip" algorithms?  The median is a better estimator
in the presence of data which one wants to reject than the average.
However, computing the median is slower than the average.
.le
.ls lsigma = 3., hsigma = 3. (ccdclip, crreject, sigclip, avsigclip, pclip)
Low and high sigma clipping factors for the "ccdclip", "crreject", "sigclip",
"avsigclip", and "pclip" algorithms.  They multiply a "sigma" factor
produced by the algorithm to select a point below and above the average or
median value for rejecting pixels.  The lower sigma is ignored for the
"crreject" algorithm.
.le
.ls rdnoise = "0.", gain = "1.", snoise = "0." (ccdclip, crreject)
Effective CCD readout noise in electrons, gain in electrons/DN, and
sensitivity noise as a fraction.  These parameters are used with the
"ccdclip" and "crreject" algorithms.  The values may be either numeric or
an image header keyword which contains the value.  Note that if the spectra
have been extracted from a 2D CCD image then the noise parameters must be
adjusted for background and the aperture summing.
.le
.ls sigscale = 0.1 (ccdclip, crreject, sigclip, avsigclip)
This parameter determines when poisson corrections are made to the
computation of a sigma for images with different scale factors.  If all
relative scales are within this value of unity and all relative zero level
offsets are within this fraction of the mean then no correction is made.
The idea is that if the images are all similarly though not identically
scaled, the extra computations involved in making poisson corrections for
variations in the sigmas can be skipped.  A value of zero will apply the
corrections except in the case of equal images and a large value can be
used if the sigmas of pixels in the images are independent of scale and
zero level.
.le
.ls pclip = -0.5 (pclip)
Percentile clipping algorithm parameter.  If greater than
one in absolute value then it specifies a number of pixels above or
below the median to use for computing the clipping sigma.  If less
than one in absolute value then it specifies the fraction of the pixels
above or below the median to use.  A positive value selects a point
above the median and a negative value selects a point below the median.
The default of -0.5 selects approximately the quartile point.
See the DESCRIPTION section for further details.
.le
.ls grow = 0
Number of pixels to either side of a rejected pixel
to also be rejected.  This applies only to pixels rejected by one of
the rejection algorithms and not the threshold rejected pixels.
.le
.ls blank = 0.
Value to use when there are no input pixels to combine for an output pixel.
.le
.ih
DESCRIPTION
\fBScombine\fR combines input spectra by interpolating them (if necessary)
to a common dispersion sampling, rejecting pixels exceeding specified low
and high thresholds, scaling them in various ways, applying a rejection
algorithm based on known or empirical noise statistics, and computing the
sum, weighted average, or median of the remaining pixels.  Note that
the "sum" option is the direct summation of the pixels and does not
perform any rejection or scaling of the data regardless of the parameter
settings.

The input spectra are specified using an image list in which each image
may contain multiple spectra.  The set of spectra may be restricted
by the \fIaperture\fR parameter to specific apertures.  The set of input
spectra may then be grouped using the \fIgroup\fR parameter and each
group combined separately into a final output spectrum.  The grouping
options are to select all the input spectra regardless of the input
image or aperture number, select all spectra of the same aperture,
or select all the spectra from the same input image.

The output consists of either a single image with one spectrum for each
combined group or, when grouping by image, an image with the single
combined spectra from each input image.  The output images and
combined spectra inherit the header parameters from the first spectrum
of the combined group.  In addition to the combined spectrum an associated
integer spectrum containing the number of pixels combined
and logfile listing the combined spectra, scaling, weights, etc, may
be produced.

The spectral combining is done using pixels at common dispersion
coordinates rather than physical or logical pixel coordinates.  If the
spectra to be combined do not have identical dispersion coordinates then
the spectra are interpolated to a common dispersion sampling before
combining.  The interpolation conserves pixel values rather pixel fluxes.
This means that flux calibrated data is treated correctly and that
spectra in counts are not corrected in the interpolation for changes
in pixel widths.  
The default interpolation function is a 5th order polynomial.  The
choice of interpolation type is made with the package parameter "interp".
It may be set to "nearest", "linear", "spline3", "poly5", or "sinc".
Remember that this applies to all tasks which might need to interpolate
spectra in the \fBonedspec\fR and associated packages.  For a discussion of
interpolation types see \fBonedspec\fR.

There are two choices for the common dispersion coordinate sampling. If the
\fIfirst\fR parameter is set then the dispersion sampling of the first
spectrum is used.  This dispersion system may be nonlinear.  If the
parameter is not set then the user specified linear or log linear
dispersion system is used.  Any combination of starting wavelength, ending
wavelength, wavelength per pixel, and number of output pixels may be
specified.  Unspecified values will default to reasonable values based on
the minimum or maximum wavelengths of all spectra, the minimum dispersion,
and the number of pixels needed to satisfy the other parameters.  If the
parameters overspecify the linear system then the ending wavelength is
adjusted based on the other parameters.  Note that for a log linear system
the wavelengths are still specified in nonlog units and the dispersion is
finally recalculated using the difference of the log wavelength endpoints
divided by the number pixel intervals (the number of pixels minus one).

There are several stages to combining a selected group of spectra.  The
first is interpolation to a common dispersion sampling as discussed
above.  The second stage is to eliminate any pixels outside the specified
thresholds.  Note that the thresholds apply to the interpolated
spectra.  Scaling and zero offset factors are computed and applied to the
spectra if desire.  The computation of these factors as well as weights is
discussed in the following section.  Next there is a choice of rejection
algorithms to identify and eliminate deviant pixels.  Some of these are
based on order statistics and some relative to the distance from an initial
median or average using a noise model cutoff.  A growing factor may be
applied to neighbors of rejected pixels to reject additional pixels.  The
various algorithms are described in detail in a following section.
Finally, the remaining pixels are combined by summing (which may not be
appropriate when pixels are rejected), computing a median, or computing a
weighted or unweighted average.  The combined spectrum is written to an
output image as well the number of pixels used in the final combining.

SCALES AND WEIGHTS

In order to combine spectra with rejection of pixels based on deviations
from some average or median they must be scaled to a common level.  There
are two types of scaling available, a multiplicative intensity scale and an
additive zero point shift.  The intensity scaling is defined by the
\fIscale\fR parameter and the zero point shift by the \fIzero\fR
parameter.  These parameters may take the values "none" for no scaling,
"mode", "median", or "mean" to scale by statistics of the spectrum pixels,
"exposure" (for intensity scaling only) to scale by the exposure time
keyword in the image header, any other image header keyword specified by
the keyword name prefixed by the character '!', and the name of a file
containing the scale factors for the input image prefixed by the
character '@'.

Examples of the possible parameter values are shown below where
"myval" is the name of an image header keyword and "scales.dat" is
a text file containing a list of scale factors.

.nf
	scale = none		No scaling
	zero = mean		Intensity offset by the mean
	scale = exposure	Scale by the exposure time
	zero = !myval		Intensity offset by an image keyword
	scale = @scales.dat	Scales specified in a file
.fi

The spectrum statistics factors are computed within specified sample
regions given as a series of colon separated wavelengths.  If no
regions are specified then all pixels are used.  If the
wavelength sample list is too long the regions can be defined in a file and
specified in the \fIsample\fR parameter using the syntax @<file> where file
is the filename.

The statistics are as indicated by their names.  In particular, the
mode is a true mode using a bin size which is a fraction of the
range of the pixels and is not based on a relationship between the
mode, median, and mean.  Also thresholded pixels are excluded from the
computations as well as during the rejection and combining operations.

The "exposure" option in the intensity scaling uses the value of the image
header keyword (EXPTIME, EXPOSURE, or ITIME).  Note that the exposure
keyword is also updated in the final image as the weighted average of the
input values.  If one wants to use a nonexposure time keyword and keep the
exposure time updating feature the image header keyword syntax is
available; i.e. !<keyword>.

Scaling values may be defined as a list of values in a text file.  The file
name is specified by the standard @file syntax.  The list consists of one
value per line.  The order of the list is assumed to be the same as the
order of the input spectra.  It is a fatal error if the list is incomplete
and a warning if the list appears longer than the number of input spectra.
Consideration of the grouping parameter must be included in
generating this list since spectra may come from different images,
some apertures may be missing, and, when there are multiple output spectra
or images, the same list will be repeatedly used.

If both an intensity scaling and zero point shift are selected the
multiplicative scaling is done first.  Use of both makes sense for images
if the intensity scaling is the exposure time to correct for
different exposure times and with the zero point shift allowing for
sky brightness changes.  This is less relevant for spectra but the option
is available.

The spectrum statistics and scale factors are recorded in the log file
unless they are all equal, which is equivalent to no scaling.  The
intensity scale factors are normalized to a unit mean and the zero
point shifts are adjusted to a zero mean.  When scal factors
or zero point shifts are specified by the user in an @file or by an
image header keyword, no normalization is done.

Scaling affects not only the mean values between spectra but also the
relative pixel uncertainties.  For example scaling an spectrum by a
factor of 0.5 will reduce the effective noise sigma of the spectrum
at each pixel by the square root of 0.5.  Changes in the zero
point also changes the noise sigma if the spectrum noise characteristics
are Poissonian.  In the various rejection algorithms based on
identifying a noise sigma and clipping large deviations relative to
the scaled median or mean, one may need to account for the scaling induced
changes in the spectrum noise characteristics.

In those algorithms it is possible to eliminate the "sigma correction"
while still using scaling.  The reasons this might be desirable are 1) if
the scalings are similar the corrections in computing the mean or median
are important but the sigma corrections may not be important and 2) the
spectrum statistics may not be Poissonian, either inherently or because the
spectra have been processed in some way that changes the statistics.  In the
first case because computing square roots and making corrections to every
pixel during the iterative rejection operation may be a significant
computational speed limit the parameter \fIsigscale\fR selects how
dissimilar the scalings must be to require the sigma corrections.  This
parameter is a fractional deviation which, since the scale factors are
normalized to unity, is the actual minimum deviation in the scale factors.
For the zero point shifts the shifts are normalized by the mean shift
before adjusting the shifts to a zero mean.  To always use sigma scaling
corrections the parameter is set to zero and to eliminate the correction in
all cases it is set to a very large number.

If the final combining operation is "average" then the spectra may be
weighted during the averaging.  The weights are specified in the same way
as the scale factors.  The weights, scaled to a unit sum, are printed in
the log output.

The weights are only used for the final weighted average and sigma image
output.  They are not used to form averages in the various rejection
algorithms.  For weights in the case of no scaling or only multiplicative
scaling the weights are used as given or determined so that images
with lower signal levels will have lower weights.  However, for
cases in which zero level scaling is used the weights are computed
from the initial weights (the exposure time, image statistics, or
input values) using the formula:

.nf
	weight_final = weight_initial / (scale * zero)
.fi

where the zero values are those before adjustment to zero mean over
all images.  The reasoning is that if the zero level is high the sky
brightness is high and so the S/N is lower and the weight should be lower.


THRESHOLD REJECTION

There is an initial threshold rejection step which may be applied.  The
thresholds are given by the parameters \fIlthreshold\fR and
\fIhthreshold\fR.  Values of INDEF mean that no threshold value is
applied.  Threshold rejection may be used to exclude very bad pixel values
or as a way of masking images.  The former case is useful to exclude very
bright cosmic rays.  Some of the rejection algorithms, such as "avsigclip",
can perform poorly if very strong cosmic rays are present.  For masking one
can use a task like \fBimedit\fR or \fBimreplace\fR to set parts of the
spectra to be excluded to some very low or high magic value.


REJECTION ALGORITHMS

The \fIreject\fR parameter selects a type of rejection operation to
be applied to pixels not thresholded.  If no rejection
operation is desired the value "none" is specified.  This task is
closely related to the image combining task \fBimcombine\fR and, in
particular, has the same rejection algorithms.
Some the algorithms are more appropriate to images but are available
in this task also for completeness.

MINMAX
.in 4
A specified fraction of the highest and lowest pixels are rejected.
The fraction is specified as the number of high and low pixels, the
\fInhigh\fR and \fInlow\fR parameters, when data from all the input spectra
are used.  If pixels are missing where there is no overlap or have been
rejected by thresholding then a matching fraction of the remaining pixels,
truncated to an integer, are used.  Thus,

.nf
	nl = n * nlow/nspectra + 0.001 
	nh = n * nhigh/nspectra + 0.001 
.fi

where n is the number of pixels to be combined, nspectra is the number
of input spectra, nlow and nhigh
are task parameters and nl and nh are the final number of low and
high pixels rejected by the algorithm.  The factor of 0.001 is to
adjust for rounding of the ratio.

As an example with 10 input spectra and specifying one low and two high
pixels to be rejected the fractions to be rejected are 0.1 and 0.2
and the number rejected as a function of n is:

.nf
	 n   0  1  2  3  4  5  6  7  8  9 10
	 nl  0  0  0  0  0  1  1  1  1  1  2
	 nh  0  0  0  0  0  0  0  0  0  0  1
.fi
.in -4
CCDCLIP
.in 4
If the noise characteristics of the spectra can be described by fixed
gaussian noise, a poissonian noise which scales with the square root of
the intensity, and a sensitivity noise which scales with the intensity,
the sigma in data values at a pixel with true value <I>,
as approximated by the median or average with the lowest and highest value
excluded, is given as:

.nf
	sigma = ((rn / g) ** 2 + <I> / g + (s * <I>) ** 2) ** 1/2
.fi

where rn is the read out noise in electrons, g is the gain in
electrons per data value, s is a sensitivity noise given as a fraction,
and ** is the exponentiation operator.  Often the sensitivity noise,
due to uncertainties in the pixel sensitivities (for example from the
flat field), is not known in which case a value of zero can be used.

This model is typically valid for CCD images.  During extraction of 
spectra from CCD images the noise parameters of the spectrum pixels
will be changed from those of the CCD pixels.  Currently it is up to
the user to determine the proper modifications of the CCD read noise
gain, and sensitivity noise.

The read out noise is specified by the \fIrdnoise\fR parameter.  The value
may be a numeric value to be applied to all the input spectra or an image
header keyword containing the value for spectra from each image.
Similarly, the parameter \fIgain\fR specifies the gain as either a value or
image header keyword and the parameter \fIsnoise\fR specifies the
sensitivity noise parameter as either a value or image header keyword.

The algorithm operates on each output pixel independently.  It starts by
taking the median or unweighted average (excluding the minimum and maximum)
of the unrejected pixels provided there are at least two input pixels.  The
expected sigma is computed from the CCD noise parameters and pixels more
that \fIlsigma\fR times this sigma below or \fIhsigma\fR times this sigma
above the median or average are rejected.  The process is then iterated
until no further pixels are rejected.  If the average is used as the
estimator of the true value then after the first round of rejections the
highest and lowest values are no longer excluded.  Note that it is possible
to reject all pixels if the average is used and is sufficiently skewed by
bad pixels such as cosmic rays.

If there are different CCD noise parameters for the input images
(as might occur using the image header keyword specification) then
the sigmas are computed for each pixel from each image using the
same estimated true value.

If the images are scaled and shifted and the \fIsigscale\fR threshold
is exceedd then a sigma is computed for each pixel based on the
spectrum scale parameters; i.e. the median or average is scaled to that of the
original image before computing the sigma and residuals.

After rejection the number of retained pixels is checked against the
\fInkeep\fR parameter.  If there are fewer pixels retained than specified
by this parameter the pixels with the smallest residuals in absolute
value are added back.  If there is more than one pixel with the same
absolute residual (for example the two pixels about an average
or median of two will have the same residuals) they are all added
back even if this means more than \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR parameter only applies to the pixels used
by the clipping rejection algorithm and does not apply to threshold
or bad pixel mask rejection.

This is the best clipping algorithm to use if the CCD noise parameters are
adequately known.  The parameters affecting this algorithm are \fIreject\fR
to select this algorithm, \fImclip\fR to select the median or average for
the center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, the CCD noise parameters \fIrdnoise, gain\fR and \fIsnoise\fR,
\fIlsigma\fR and \fIhsigma\fR to select the clipping thresholds,
and \fIsigscale\fR to set the threshold for making corrections to the sigma
calculation for different image scale factors.

.in -4
CRREJECT
.in 4
This algorithm is identical to "ccdclip" except that only pixels above
the average are rejected based on the \fIhsigma\fR parameter.  This
is appropriate for rejecting cosmic ray events and works even with
two spectra.

.in -4
SIGCLIP
.in 4
The sigma clipping algorithm computes at each output pixel the median or
average excluding the high and low values and the sigma about this
estimate.  There must be at least three input pixels, though for this method
to work well there should be at least 10 pixels.  Values deviating by more
than the specified sigma threshold factors are rejected.  These steps are
repeated, except that after the first time the average includes all values,
until no further pixels are rejected or there are fewer than three pixels.

After rejection the number of retained pixels is checked against the
\fInkeep\fR parameter.  If there are fewer pixels retained than specified
by this parameter the pixels with the smallest residuals in absolute
value are added back.  If there is more than one pixel with the same
absolute residual (for example the two pixels about an average
or median of two will have the same residuals) they are all added
back even if this means more than \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR parameter only applies to the pixels used
by the clipping rejection algorithm and does not apply to threshold
rejection.

The  parameters affecting this algorithm are \fIreject\fR to select
this algorithm, \fImclip\fR to select the median or average for the
center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, \fIlsigma\fR and \fIhsigma\fR to select the
clipping thresholds, and \fIsigscale\fR to set the threshold for
making corrections to the sigma calculation for different spectrum scale
factors.

.in -4
AVSIGCLIP
.in 4
The averaged sigma clipping algorithm assumes that the sigma about the
median or mean (average excluding the low and high values) is proportional
to the square root of the median or mean at each point.  This is
described by the equation:

.nf
	sigma(column,line) = sqrt (gain(line) * signal(column,line))
.fi

where the \fIestimated\fR signal is the mean or median (hopefully excluding
any bad pixels) and the gain is the \fIestimated\fR proportionality
constant having units of photons/data number.

This noise model is valid for spectra whose values are proportional to the
number of photons recorded.  In effect this algorithm estimates a
photon per data value gain for each spectrum.
The gain proportionality factor is computed
independently for each output spectrum by averaging the square of the residuals
(at points having three or more input values) scaled by the median or
mean.

Once the proportionality factor is determined, deviant pixels exceeding the
specified thresholds are rejected at each point by estimating the sigma
from the median or mean.  If any values are rejected the median or mean
(this time not excluding the extreme values) is recomputed and further
values rejected.  This is repeated until there are no further pixels
rejected or the number of remaining input values falls below three.  Note
that the proportionality factor is not recomputed after rejections.

If the spectra are scaled differently and the sigma scaling correction
threshold is exceedd then a correction is made in the sigma
calculations for these differences, again under the assumption that
the noise in an spectra scales as the square root of the mean intensity.

After rejection the number of retained pixels is checked against the
\fInkeep\fR parameter.  If there are fewer pixels retained than specified
by this parameter the pixels with the smallest residuals in absolute
value are added back.  If there is more than one pixel with the same
absolute residual (for example the two pixels about an average
or median of two will have the same residuals) they are all added
back even if this means more than \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR parameter only applies to the pixels used
by the clipping rejection algorithm and does not apply to threshold
rejection.

This algorithm works well for even a few input spectra.  It works better if
the median is used though this is slower than using the average.  Note that
if the spectra have a known read out noise and gain (the proportionality
factor above) then the "ccdclip" algorithm is superior.  However, currently
the CCD noise characteristics are not well propagated during extraction so
this empirical algorithm is the one most likely to be useful.  The two
algorithms are related in that the average sigma proportionality factor is
an estimate of the gain.

The  parameters affecting this algorithm are \fIreject\fR to select
this algorithm, \fImclip\fR to select the median or average for the
center of the clipping, \fInkeep\fR to limit the number of pixels
rejected, \fIlsigma\fR and \fIhsigma\fR to select the
clipping thresholds, and \fIsigscale\fR to set the threshold for
making corrections to the sigma calculation for different image scale
factors.

.in -4
PCLIP
.in 4
The percentile clipping algorithm is similar to sigma clipping using the
median as the center of the distribution except that, instead of computing
the sigma of the pixels from the CCD noise parameters or from the data
values, the width of the distribution is characterized by the difference
between the median value and a specified "percentile" pixel value.  This
width is then multipled by the scale factors \fIlsigma\fR and \fIhsigma\fR
to define the clipping thresholds above and below the median.  The clipping
is not iterated.

The pixel values at each output point are ordered in magnitude and the
median is determined.  In the case of an even number of pixels the average
of the two middle values is used as the median value and the lower or upper
of the two is the median pixel when counting from the median pixel to
selecting the percentile pixel.  The parameter \fIpclip\fR selects the
percentile pixel as the number (if the absolute value is greater
than unity) or fraction of the pixels from the median in the ordered set.
The direction of the percentile pixel from the median is set by the sign of
the \fIpclip\fR parameter with a negative value signifying pixels with
values less than the median.  Fractional values are internally converted to
the appropriate number of pixels for the number of input spectra.  A minimum
of one pixel and a maximum corresponding to the extreme pixels from the
median are enforced.  The value used is reported in the log output.  Note
that the same percentile pixel is used even if pixels have been rejected by
nonoverlap or thresholding; for example, if the 3nd pixel below
the median is specified then the 3rd pixel will be used whether there are
10 pixels or 5 pixels remaining after the preliminary steps.

After rejection the number of retained pixels is checked against the
\fInkeep\fR parameter.  If there are fewer pixels retained than specified
by this parameter the pixels with the smallest residuals in absolute
value are added back.  If there is more than one pixel with the same
absolute residual (for example the two pixels about an average
or median of two will have the same residuals) they are all added
back even if this means more than \fInkeep\fR pixels are retained.
Note that the \fInkeep\fR parameter only applies to the pixels used
by the clipping rejection algorithm and does not apply to threshold
or bad pixel mask rejection.

Some examples help clarify the definition of the percentile pixel.  In the
examples assume 10 pixels.  The median is then the average of the
5th and 6th pixels.  A \fIpclip\fR value of 2 selects the 2nd pixel
above the median (6th) pixel which is the 8th pixel.  A \fIpclip\fR
value of -0.5 selects the point halfway between the median and the
lowest pixel.  In this case there are 4 pixels below the median,
half of that is 2 pixels which makes the percentile pixel the 3rd pixel.

The percentile clipping algorithm is most useful for clipping small
excursions, such as the wings of bright lines when combining
disregistered observations, that are missed when using
the pixel values to compute a sigma.  It is not as powerful, however, as
using the CCD noise parameters (provided they are accurately known) to clip
about the median.  This algorithm is primarily used with direct images
but remains available for spectra.

The  parameters affecting this algorithm are \fIreject\fR to select this
algorithm, \fIpclip\fR to select the percentile pixel, \fInkeep\fR to limit
the number of pixels rejected, and \fIlsigma\fR and \fIhsigma\fR to select
the clipping thresholds.


.in -4
GROW REJECTION

Neighbors of pixels rejected by the rejection algorithms
may also be rejected.  The number of neighbors to be rejected on either
side is specified by the \fIgrow\fR parameter.

This rejection step is also checked against the \fInkeep\fR parameter
and only as many pixels as would not violate this parameter are
rejected.  Unlike it's application in the rejection algorithms at
this stage there is no checking on the magnitude of the residuals
and the pixels retained which would otherwise be rejected are randomly
selected.


COMBINING

After all the steps of offsetting the input images, masking pixels,
threshold rejection, scaling, and applying a rejection algorithms the
remaining pixels are combined and output.  The pixels may be combined
by computing the median or by computing a weighted average.
.ih
EXAMPLES
1.  Combine orders of echelle images.

.nf
	cl> scombine *.ec *%.ec%% group=images combine=sum
.fi

2.  Combine all spectra using range syntax and scale by the exposure times.

.nf
	cl> names irs 10-42 > irs.dat
	cl> scombine @irs.dat irscombine group=all scale=exptime
.fi

3.  Combine spectra by apertures using exposure time scaling and weighting.

.nf
	cl> scombine *.ms combine.ms nout=ncombine.ms \\
	>>> group=apertures scale=exptime weights=exptime
.fi
.ih
REVISIONS
.ls SCOMBINE V2.10.3
The weighting was changed from using the square root of the exposure time
or spectrum statistics to using the values directly.  This corresponds
to variance weighting.  Other options for specifying the scaling and
weighting factors were added; namely from a file or from a different
image header keyword.  The \fInkeep\fR parameter was added to allow
controlling the maximum number of pixels to be rejected by the clipping
algorithms.  The \fIsnoise\fR parameter was added to include a sensitivity
or scale noise component to the noise model.
.le
.ls SCOMBINE V2.10
This task is new.
.le
.ih
NOTES
The pixel uncertainties and CCD noise model are not well propagated.  In
particular it would be desirable to propagate the pixel uncertainties
and CCD noise parameters from the initial CCD images.
.ih
SEE ALSO
imcombine, odcombine, lscombine
.endhelp