path: root/noao/twodspec/multispec/doc/MSalgo.ms

.de FX
.nf
.ps -2
.ss 25
.cs R 25
..
.de EX
.ps +2
.ss
.cs R
.fi
..
.EQ
delim $$
.EN
.RP
.TL
Algorithms for the Multi-Spectra Extraction Package
.AU
Francisco Valdes
.K2
.TU
.AB
The algorithms for the Multi-Spectra Extraction Package (\fBmultispec\fR)
in the Image Reduction and Analysis Facility (\fBIRAF\fR) is described.
The basic aspects of the general two dimensional aperture spectra extraction
problem are first discussed.
The specific algorithms for extraction of multi-aperture plate and
Echelle digital data are presented.  Results of the authors experiments
with this type of data are included.
The detailed specification of the package is given in a second document,
\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB.
.AE
.NH
Introduction
.PP
There are an increasing number of astronomical instruments which produce
multiple spectra on a two dimensional detector.
The basic concept is to use one dimension for wavelength,
the dispersion dimension, and the other, the cross dimension, for
packing additional information during a single exposure.
For example, the cross dimension can be different objects or
different spectral orders.  The classic multi-spectra instrument is
the Echelle spectrograph.  New instruments are the aperture plate and
Medusa spectrographs.
.PP
There is an additional aspect of the multi-spectra format; namely,
the individual spectra can contain spatial data.  An example of
this would be multiple slit spectra in which each slit spectra contains
sky signal and object signal.  In the following
discussion we limit the spectra to be simple aperture spectra in
which we only desire to sum the intensities at each wavelength to form
a one dimensional spectrum.
.PP
The analysis of multi-spectra aperture data consists of two steps; the
separation and extraction into individual aperture spectra
and the calibration and measurement of the spectra.  These steps can
either be incorporated into one analysis package or two separate
packages.  There are advantages to the first approach since some
aspects of the individual spectra are directly related by the physical
geometry of the multi-spectra format.  However, because packages for
the analysis of individual spectra exist we begin by dividing the
reduction into separate extraction and analysis tasks.  It is
important to realize, however, that the existing analysis tools are not well
suited to reducing the larger number of spectra and treating sets of
spectra together.
.PP
The latter part of this paper describes the algorithms for the
extraction of two types of data; the multi-aperture plate (MAP)
and Echelle used with digital detectors.  However,
it is important to keep the more general problem in mind
and the remainder of this introduction considers the different
conceptual aspects of the multi-spectra extraction task.
Table 1 lists many of the general properties of multi-spectra aperture data.
The other two columns give possible alternatives that each property may take.

.TS
center box;
c s s
c s s
c c s
= = =
c || c | c.
Table 1:  Aspects of Multi-Spectral Data

Property	Alternatives
detector	digital	photographic
alignment	aligned	skewed
blending	blended	unblended
aperture	holes	slits
spectral resolution	low	high
.TE

.PP
The detector determines what kind of calibration procedures are
required to produce intensity from the measurements.
A digital detector requires sensitivity calibrations on all scales.
This is the "flat field" problem.  There are also corrections for
bias and dark current.  Photographic detectors require
intensity calibration.  Data which are not aligned with the natural
dimensions of the digital image require extra procedures.  Two types
of non-alignment are a rotation of the dispersion dimension relative
to the pixel dimension and a "wiggling" or "snaking" of the dispersion
dimension.  Blending refers to the degree of contamination along the
cross dimension.  Blended data requires extra effort to correct for
the overlap between different spectra and to determine the background.
The aperture defines the extent of the spectra in the cross dimension.
The two most relevant choices are holes and slits.  In some
instruments, like the Echelle, the size of the aperture can be varied
at the time of the observations.  Finally, the spectral resolution is
important in conjunction with digital detectors.  If the resolution is
high then quartz flat field calibrations are relatively easy because
the spectral
signature need not be considered.  Otherwise, the flat field problem
is more difficult because the gain variations of the detector
must be separated from the natural spectral intensity variation of the
quartz.
.PP
There is always some confusion of terms when talking about multi-spectra
data; in particular, the terms x, y, line, and band.
The image pixel dimensions are refered to as x and y.  We assume
for the moment that the alignment of the multi-spectra format is such
that x corresponds to the cross dimension and y to the dispersion
dimension.  If this is not the case a rotation or interpolation
program can be used to approximately orient the data in this way.
A line is the set of intensity values as a function of x at constant y.
In other words, a line is a cut across the dispersion dimension.
A band is the average of more than one line.
The image residing on disk will generally be organized
such that x varies more rapidly and a line of the image is easily
obtained.  In this form a display of the image will have the spectra
running vertically.  In the Cyber extraction package the data is
organized with x corresponding to the dispersion dimension.
.NH
Multi-Spectra Image Formats
.PP
The remainder of this paper will refer to two specfic and very
different multi-spectra formats; the Kitt Peak Multi-Aperture Plate
System and the Kitt Peak Echelle Spectrograph.
.NH 2
Kitt Peak Multi-Aperture Plate System
.PP
The reduction of data from multi-aperture plate observations is the 
driving force for the development of a multi-spectra extraction
package.  This application turns out to have most of the worse aspects
of the properties listed in Table 1.  The multi-aperture plate spectrograph uses
digital dectectors with low resolution, the spectra are blended and
change alignment along the pixel dimension.  Furthermore, the camera
has a variable point-spread function and focus,
suffers from flexture problems, has a different illumination for
the quartz than object exposures, and unexplained background level
variations (in the CRYOCAM).  There are two detectors which have been
used with the multi-aperture plate system, the Cryogenic Camera
(CRYOCAM) and the High Gain Video Spectrometer (HGVS).
.NH 2
Echelle
.PP
As some of the algorithms were developed the Echelle data was brought
to my attention.  It is considerably simpler than the MAP data because
it is unblended and of high spectral resolution.
Furthermore, the quartz exposure
can be made wider than the object exposures and flexture is not a
problem.  The principle problem in this data was the
prevalence of cosmic rays.  It pointed to the need to maintain generality
in dealing with both the MAP data and other types of
multi-spectra data which have different profiles, may or may not be
merged, and may or may not have different widths in quartz and object.
Dealing with the cosmic ray problem lead to a very effective solution
usable in both the Echelle and multi-aperture plate data.
.NH
User Level Extraction Logic
.PP
The user should generally only be concerned with the logical steps of
extracting the individual spectra from the multi-spectra image.  This
means that apart from specifying the detector system and the format
he should not deal with details of the detector and the format.
In the paper, 
\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB,
the \fBIRAF\fR extraction package design and program specifications
are described.
.NH
Flat Fields
.PP
There are two types of flat field situations depending on the spectral
resolution.  When the resolution is high then the spectral signature of
the continum calibration source, a quartz exposure, will be unimportant
and variations in the signal will be due to detector sensitivity variations.
In this case the quartz frame, or average of several frames, is the flat
field and division of the object frames by the quartz frame is all that
is required.  However, a special
image division program is desirable to handle the region of low or absent
signal between the the spectra.  This is described in section 4.2.
.PP
In the alternate case of lower resolution the quartz spectral signature is
larger than the detector response variations.  A flat
field in which the intrinsic quartz spectrum is removed is produced by
assuming that the true value of a pixel is given by the smoothed average
of the pixels near that point in position and wavelength and taking
the ratio of the data value to the smoothed value.
This requires a special smoothing program described in section 4.1.
After the flat field is generated then the same image division
program used for the Echelle data is applied.
The image division and smoothing programs are general image operators and
not specific to the Multi-Spectra Extraction Package.
.NH 2
MULTISPEC_FLAT
.PP
The multi-aperture plate data varies in both dimensions.  Thus, any averaging
to smooth the image must take this variation into account.  In the Cyber
a flat field for the multi-aperture plate data smooths across the dispersion
by modeling the spectra.  This is a difficult task to do accurately because
the true shape of the spectra is not known and the counts vary greatly
and rapidly in this dimension.  This approach has the further difficulty
that it is not possible to average several quartz exposures directly.
.PP
The alternate approach to modeling is statistical averaging.
Averaging across the dispersion requires very high order polynomials
because of the rapid variations;
the spectra are typically spaced about 8 pixels apart and there are on
the order of 50 spectra.  On the other hand, variations along the dispersion
are much slower even if the spectra are slightly skewed; a bad case would
have two peaks in 800 pixels along the y dimension.  This kind
of variation is tractable with relatively simple averaging polynomials
and is the one adopted for the multi-aperture plate data.
.PP
The flat fields are produced by a quadratic moving average along the
y direction.  This means that the region centered at a given pixel
is fitted by a least-squares quadratic polynomial and the value of the
polynomial at that point is the appropriate statistical average.
The width of the moving average is an adjustable parameter.
At the edges of the frame where it is not possible to center a region of
the specified width about the desired pixel the polynomial fit is used to
extrapolate the average value to the edge.
.PP
Because the quadratic fit will
be influenced by bad pixels an attempt is made to detect and smooth over
the bad pixels.  This is accomplished by comparing the smoothed values to
the observed values and ignoring pixels with a value of

.EQ (1)
	chi = | observed - smoothed | / sqrt smoothed
.EN

greater than a specified value.  Then the smoothing is recalculated and tested
again for bad pixels.  This iteration continues until no further pixels
are rejected.
.PP
Following the smoothing the flat field is produced by ratioing the raw quartz
to the smoothed quartz.  Pixels of low signal (specified by the
parameter \fIconstant\fR )
are treated by the equation

.EQ
	r = (data + (constant - smoothed) ) / constant  .
.EN

The resultant flat field image is then divided into the object frames in
the manner described in the next section.
.PP
Experience with data from the Cryogenic Camera has proved very good.
The flat field which is produced can be examined on a display.  It
shows fringing at red wavelengths and is not too strongly affected
by bad pixels.  Some further effort, however, could go into smoothing
over the bad pixels.
.PP
The smoothing operation on data from the Cryogenic Camera actually
consists of four steps.  The quartz exposures are first averaged.
The average quartz is rotated so that the dispersion
direction is the most rapidly varying or x dimension.  Then the
smoothing is performed along x followed by another rotation to return
the flat field image to its original orientation.  The reason for the
rotations is that they can be done quickly and efficiently whereas
smoothing along the y dimension is very slow and coding an efficient
version is much more complicated than doing a single line at a time.
.NH 2
FLAT_DIVIDE
.PP
The Echelle data has quartz frames which can be used directly as flat fields.
One just has to divide the object frames by the quartz or average of several
quartz.  However, in image division consideration has to be given the
problem of division by zero or very small numbers.  In direct imaging this
may not be much of a problem but in multi-spectra data the region between
the spectra and near the edges of the spectra will have very low counts.
Another aspect of image division for making flat field corrections is the
scaling of the result.  The flat field integer image data must be large
to give accurate relative response values.  However, one wants to divide
an object frame by values near unity.
This section describes a special image division operator allowing the user
to specify how to handle these cases.
.PP
The parameters are a \fIdivision threshold\fR
(default of zero) and a \fIthreshold violation value\fR.  Values of the
denominator above the \fIthreshold\fR are treated separatedly from those
below the \fIthreshold\fR.  The denominator image is scaled to have an
average of one for pixels above the \fIthreshold\fR.  The pixel by pixel
division is then performed for those points for which the denominator
is above the \fIthreshold\fR.  Pixels for which the denominator is below the
\fIthreshold\fR are set to the \fIthreshold violation value\fR in the resultant
image if the \fIviolation value\fR is specified.  If the value is not
specified then the numerator value is taken as the resultant value.
The divisions can be done in place or the result put into a new image file.
.PP
For the multi-spectra situation where the object spectra have a
smaller width than the quartz, as in the Echelle, one can either
set the \fIthreshold
violation value\fR to zero or not set it at all resulting in either
exactly zero or background values between the spectra while still flattening
the spectra.  This allows looking at the flattened spectra without the
annoying "grass" between the spectra caused by dividing by small
values.
.NH
Extraction
.NH 2
MULTIAP_EXTRACT
.PP
The extraction of spectra from multi-aperture plate images consists of
a series of steps.  The steps are executed from a script.
The command

.FX
ms> multiap_extract "ap165.*", "", 165, 50
.EX

will take the flattened images, ap165.*, from aperture plate 165 with 50
spectra and automatically locate the spectra, model the profiles, and
extract the one dimensional spectra.  The script consists of
a number of steps as described below.
.PP
\fBFind_spectra\fR (section 6) initializes the \fBmultispec\fR data file
and does a peak search to determine the initial positions of the
spectra.
\fBFind_bckgrnd\fR fits a polynomial of order 1 (or more) for the pixels which
are not near the spectra as defined by \fBfind_spectra\fR.
.PP
The spectra are then modeled in bands of 32 lines by the model profiles
described in section 8.1.  The first \fBmodel_fit\fR uses three Gaussian
parameters for
each spectra measuring the peak intensity, peak position, and width.
The second \fBmodel_fit\fR adds a fourth parameter to modify the wings of the
profile.
.PP
The \fBmodel_extract\fR program extracts the spectra line by line and also
detects and removes cosmic ray events which do not fit the model
profiles (see section 9).
In outline, the extraction of blended spectral data uses the
model profiles to determine the fraction of light
from each of the neighboring spectra at the pixel in question.  The
appropriate fraction of the
.ul
observed
pixel intensity (minus the background) is
assigned to the luminosities of the spectra.  There are two versions
of the \fBmodel_extract\fR extraction.  The first simultaneously fits the
peak intensity of all the spectra and the second uses the
data value at the peak of each spectra to normalize the model.  The
first method is slow and accurate and the second is fast and approximate.  
Because the models are used in extraction only to define the relative
contributions of neighboring spectra to the total observed pixel luminosity
the speed of the approximate method far outweighs the need for
accuracy.  However, cleaning the spectra of cosmic rays is a different
matter and is discussed further in section 12.
.PP
After the extraction the spectra are correlated with the aperture plate
description using \fBap_plate\fR (see Section 10) to determine the
relative wavelength offsets and assign identification information to
the spectra.
.PP
For successive frames it is not necessary to resort to the initial
steps of finding the spectra and fitting from scratch.  The \fBcopy_params\fR
routine makes a new copy of the fitting database.  Small shifts in positions
of the spectra and the peak intensities are determined by doing a two
parameter fit for the peak and position using the previously determined
shape parameters.
Changes in the shape of the spectra are then determined by the three
and four parmater fits.  Because the solution is likely to be close to
the previously determined shape the transfering of one solution from a
previously solved image is faster than starting from scratch.
Note that the shapes as well as the positions and peak intensities
of all frames including the object exposures are allowed to change.
.PP
The spectra are then extracted from the image by \fBmodel_extract\fR and the
process repeats for the succeeding images.
.PP
One useful feature is the ability to specify the bands or lines to be
modeled or extracted.
This feature is useful for diagnosising the programs quickly.
The default is all bands or lines.
.NH 2
ECHELLE_EXTRACT
.PP
The extraction of the unblended Echelle spectra is performed
begins in a similar way with \fBfind_spectra\fR and \fBfind_bckgrnd\fR.
The extraction and cleaning, however, uses \fBstrip_extract\fR which
adds up the instrumental counts for each unblended spectra at each
wavelength to get the total luminosity.
.NH
FIND_SPECTRA -- Finding the Spectra
.PP
The first step in the extraction and processing of multi-spectra data is
to locate the spectra.  This can be done interactively by
the user but it is far preferable to automate the process;
particularly since multi-spectra data can have a large number of
spectra and frames.  The approach is to find the peaks in a line, or
average of lines, sort the peaks found in some manner, such as by
strength, and select the expected number of peaks from the top of the
list.
Beyond this simple outline are several algorithmic details such as how
to define and locate valid peaks and how to sort the list of peaks.
Peak finding is a general problem and a subroutine for peak finding is
described below.  The \fBfind_spectra\fR program provides an
interface between the \fBmultispec\fR data file and the
general peak finding algorithm.
.PP
The \fBpeaks\fR function takes arrays of x (position) and y (value) points
and the number of
points in the arrays and returns the number of peaks found.  It also
returns the estimated positions of the peaks in the x array and the
extimated peak values in the y array in order of peak likelihood.
There is one user parameter, the smoothing \fIwidth\fR.
The choice of the \fIwidth\fR parameter is dicatated by how closely and how
wide the peaks are to be expected.
The algorithm takes a region of \fIwidth\fR points
centered on each x point and fits a quadratic;

.EQ
y sub fit = a + b x + c x sup 2~.
.EN

A peak is defined
when the slopes, $b sub 1$ and $b sub 2$, of two neighboring points
$x sub 1$ and $x sub 2$ change
sign from positive to negative and the curvatures, $c sub 1$ and $c
sub 2$, are less than -0.001 for both points.
The quadratic polynomials define two estimated peak positions

.EQ
x sub 1 sub peak = x sub 1 - b sub 1 / (2 * c sub 1 ),~~
x sub 2 sub peak = x sub 2 - b sub 2 / (2 * c sub 2 )~.
.EN

The offsets are then normalized to give a linear interpolation
fraction
$f = ( x sub 1 sub peak - x sub 1 ) / ( x sub 2 sub peak - x sub 1 sub
peak )$ in the interval between the two points.
The estimated position of the peak is then

.EQ
x sub peak = f * ( x sub 1 - x sub 2 )
.EN

and the estimated peak value is the average value of the two quadratic
polynomials at $x sub peak$.  The curvature at the peak is
estimated by $c sub peak = c sub 1 + f * (c sub 1 - c sub 2 )$.
Finally, the peaks are sorted by the magnitude of the peak curvature.
.PP
This peak finding algorithm works quite well.  I have also used it to
automatically locate peaks in the extracted one dimensional spectra
and then do peak correlations between spectra to find a relative
wavelength solution.  Some such use of this program may be implemented
in either future additions to the Multi-Spectra Extraction Package or
the Spectral Reduction Package.
.PP
In \fBfind_spectra\fR the number of spectra to be found is specified by
the user.  The user should have previously looked at an image
on a display or done a profile plot across the
dispersion to count the observed spectra.
Additional parameters specify the columns in which the spectra
are to be found and the minimum separation and width of the spectra.
The column specification allows the elimination of problems with defective
areas of the detector such as the LED in the Cryogenic Camera.  The minimum
width and separation provide algorithmic constraints to the spectra finding
procedure.
.PP
The peaks are found at two or more points in the
multi-spectra image for a band of 32 lines using a
\fBpeaks\fR \fIwidth\fR parameter of 5.  After the peaks are found
at a number of bands in the image a linear fit is made to determine any small
slope of the spectra relative to the columns.
The reason for specifying only a few bands is that the process of
finding the peaks is moderately slow and only two bands are needed for
determining the initial position angle of the spectra to the y
dimension.  Furthermore, some bands do not give a satisfactory result
because of extraneous data such as the LED in the CRYOCAM or bad
focus.  Another possibility is that a spectrum may go off the edge
and in order to "find" the spectrum for partial extraction bands that
include the on image part of the spectrum must be specified.
.NH
FIND_BCKGRND -- Background
.PP
The background on a multi-spectra image is the result of very broad
scattering as opposed to the narrower scattering which produces
distinguishable wings on individual spectra.
Modeling of the background in a Cryogenic Camera multi-aperture plate
image shows that the background is well explained by a broad
scattering function.
It is not reasonable, however, to model the scattering to this detail
in actual extractions.
Instead a smooth polynomial is fitted to the pixels not covered by
spectra.  The order of the polynomial is a specified parameter.
For the CRYOCAM MAP data a quadratic is appropriate.
.PP
The algorithm is the same for all multi-spectra data except for the
choice of parameters.  First, the location of the spectra must be
determined.  This is done by the \fBfind_spectra\fR program.  There
are two principle parameters, a buffer distance and the order of the
fitting polynomial.  Each line, or average of several lines, is fitted
by least-squares for the points lying farther than the buffer
distance from any spectra.  If there are no points which completely
stradle the spectra, i.e. located on each side of the spectra, then
the order of the fitting polynomial is ignored and a constant, or
first order polynomial, is determined.
A hidden parameter specifying the columns allowed for searching for
background points is available so that bad parts of the image can be
ignored.
.PP
A difference in philosophy with the Cyber programs is that the
determined background is not subtracted from the image data.  It is
instead kept in the database for the image.  Generally, it is better to
modify the basic image data as little as possible.  This is the approach
used in the Multi-Spectra Extraction Package.
.NH
Spectra Profiles
.NH 2
MODEL_FIT -- Models for Multi-Spectra Images
.PP
The object of modeling is to separate blended spectra for extraction
and to remove artifacts, such as cosmic rays, which do not fit
the model.  The models should have the minimum number of parameters
which give residuals approaching the detector statistics, they
should incorporate constraints  based on the physics of the
detector/camera system,  and the models must be ammenable to a
statistical fitting algorithm which is stable.
There are a large number of possibilities.
.PP
An important point to bear in mind during the following discussion is
the necessary accuracy of the model fitting.  In the design proposed
here the model fitting is not used for determining the smooth quartz.
Use of a model for making a flat field would require a very accurate
model and using an average quartz is not possible.  However, for
extraction the model is used only to indicate the
relative fraction of light for each spectrum when the spectra are
blended.  The cleaning application is more critical but not nearly so
much as in the flat field modeling.  Thus, though we do a good job of
model fitting (better the the Cyber modeling) some additional features
such as smoothing along the spectra are not included.
Also, though some improvement can be gained by the additional shape
parameters in the fit, they are not necessary for the required purpose
and can be left out leading to a faster extraction.
.PP
During the course of my investigation I tried more than one hundred
models and combinations of constraints.  Some general results of this
study follow.
The model which I find gives the best results has six parameters not
counting the background.  The model is defined by the following
equations where x is the cross dimension.

.EQ (1)
I = I sub 0 exp (-s * ( DELTA x sup 2 ))
.EN
.EQ
DELTA x = (x - x sub 0 )
.EN
.EQ
s = s sub 0 + s sub 1 y sup 2 + s sub 2 y sup 3
.EN
.EQ
y = DELTA x / sqrt { DELTA x sup 2 + x sub 1 sup 2 }
.EN

The model consists of a intensity scale parameter, $I sub 0$,
and a profile which is
written in a Gaussian form.  The center of the profile is given by
the parameter $x sub 0$.  The profile is not exactly Gaussian because the
scale, $s$, is not a constant but depends on $DELTA x$.  The scale
function has three terms; a constant term, $s sub 0$, which is the scale
near the center of the profile, and even and odd terms, $s sub 1$
and $s sub 2$,
which change the scale in the wings of the profile.
.PP
The characteristic of the profile which must be satisfied is that at
large distances from the profile center the scale is positive.  This
requirement means that the profile will be monotonically decreasing at
large distances and will have a finite luminosity.  This point was
crucial in determining the form of the scale function.  A straight
power series in $DELTA x$ does not work because power series diverge.
Instead, the scale function is defined in terms of a separation
variable $y$ which is bounded by -1 and 1 at infinite separation and is
zero at zero separation.  The parameter $x sub 1$ defines a characteristic
distance where the character of $y$ changes from going as $DELTA x$ to
asymptotic to 1.  The parameters are, thus, $I sub 0$, $x sub 0$, $s sub 0$,
$s sub 1$, $s sub 2$, $x sub 1$.
.PP
An important property of this model is that the terms have a physical
interpretation.  The profile scale and profile center are obvious and
any model must include them.  It is the remaining terms, $s sub 0$,
$s sub 1$, $s sub 2$,
and $x sub 1$, which are called the shape parameters, which are interesting.
In an ideal aperture plate system the shape of a profile would be
given by the projection of the circular aperture into the cross dimension:

.EQ
P( DELTA x ) = sqrt {1 - a DELTA x sup 2}
.EN

where the constant a is related to the size of the hole by 

.EQ
a = 1 / r sup 2
.EN

For small $DELTA x$ the profile can be expressed in the Gaussian form with
a scale

.EQ
s = a( 1/2  + a DELTA x sup 2 + ...)
.EN

Thus, even in a perfect aperture plate system a Gaussian form shows the
scale increasing from a central value determined by the size of the hole.
In other words, the profile decreases more sharply than a Gaussian.
.PP
However, no aperture plate system is ideal because the thickness of
the aperture plate is finite and there is scattering and changes in
the focus of the system.  One must
convolve the profile above with a scattering/focus function.  One can show
that for reasonable functions, exponential and Gaussian,
with some scale b the final profile is a function of the ratio b/a.
If the ratio is less than 1 then the profile will be more like that of
the hole and the profile will be sharper than a Gaussian in the wings.
If the ratio is much greater than 1 then the profile will become the
scattering profile at large separations.  Simulations using Gaussian
and exponential scattering profiles show behaviors very much like the
profile (1) with $s sub 1$ greater than zero when b/a < 1
meaning the profile becomes sharper (than a Gaussian) in the wings
and $s sub 1$ < 0 when b/a > 1.
Thus, $s sub 1$ defines the scale of the scattering profile relative
to the hole size.
The size of the hole is incorporated into the parameter $x sub 1$.
The parameter $s sub 2$ allows an asymmetry in the profile.
.PP
An interesting property of the scale function is that it is all right
for it to be negative at small distances from the profile center.  This
occurs when $s sub 0$ is negative.  The effect of this, provided $s$
becomes positive at large distances, is to give a two horned profile.
This, in fact, is observed when the focus of the system becomes very
poor.
.PP
The best fits (least chi-square or rms residual) are
obtained when each spectrum at each wavelength has independent
parameters.  However, this sometimes gives rise to unphysical results.
If left entirely unconstrained the parameter fitting algorithm can
make one line broad and dominant and a neighboring line weak and
sharp.
This is not, of course, a property of the camera or detector.
Thus, constraints based on the physics of the
camera/detector are used.  This means that the shape
parameters $s sub 0$, $s sub 1$, $s sub 2$, and $x sub 1$
are coupled locally by making them vary as a polynomial of position
across the dispersion.  One might also
constrain the variation of the shape along the spectrum as is done in
the Cyber.  This is not needed because there are no drastic differences
between the fitted parameters at neighboring points along the spectra.
.PP
My experience with the Cyrogenic Camera system has shown the
following.  The focus ripples twice across the CCD with the
propagation angle being approximately 30 degrees from the long dimension.
The change in focus is partly just a scale change.  This is seen in
the correlation of $s sub 0$ with the image scale found by \fBap_plate\fR.
The shape parameter $s sub 1$ changes sign from positive to
negative indicating that when the focus is good the profile
decreases faster than a Gaussian and when the focus is bad it decreases
slower.  Occassionally the focus is very bad and $s sub
0$ is negative and $s sub 1$ is small and positive causing a broad two
horned profile.  The
assymetry parameter, $s sub 2$, is useful only when the signal is strong near
the peak of a quartz exposure.  It is not really necessary to include
it in the model fits.  The assymetry parameter was dominant, however,
in some Echelle data which were clearly asymmetric.  The value of
$x sub 1$ is
not highly sensitive and can be fixed for a given hole size.  Large
changes in the hole size would require resetting $x sub 1$.
The use of the four parameters, $I sub 0$, $x sub 0$, $s sub 0$,
and $s sub 1$, allow good fits
to all the data I've examined including those in which the peak/valley
intensity ratio across the spectra is about 1.1.  It is the importance
of the parameter $s sub 1$ which improves the fitting dramatically over the
Cyber three parameter fitting (in addition to a different fitting
algorithm).
.PP
The heart of profile fitting is the solution of the multi-parameter
least-squares problem.  In a blended multi-spectra image the profile
parameters of one spectra are affected by its neighbors which are,
in turn, affected by their neighbors and so forth.  The key to this
type of problem is to realize that only nearest neighbors affect the
profile parameters and this leads to a "banded" least-squares matrix.
A banded matrix is one in which cross terms far from the diagonal are
zero.  Solution of banded matrices is much more efficient than solving
the entire matrix.  This allows solution for more than 100 parameters
simultaneously in a short time.
.PP
Use of the banded multi-parameter solution has the restriction, however,
that there can be no parameters in the model which are not local to
the profiles.  This affects the way
global constraints are applied to the parameters.  In particular,
the way the shape parameters are constrained to vary smoothly across the
detector.
The shape parameters are first found as independent parameters by the
banded matrix solution and then smoothed by a polynomial in x.
.PP
An area which was extensively investigated was the appropriate
weighting to use for the model fitting.  The most likely choices are
weighting by $1 / sqrt data$ and unit weight corresponding to
$chi sup 2$
and least squares fitting.  It was found that the two methods
agreed fairly closely but that the least squares fitting was more
appropriate because the blending correction depends largely on the
value of the peak intensity and less on the exact fit of the wings.
With $chi sup 2$ the peak is fit with less accuracy in order to improve
the fit in the wings of the profile.  In some cases this gave clear
errors in estimating the peak intensity and, hence, the proper contributions
between the blended spectra were not made.
.PP
Now follows the details of the fitting algorithm.
The algorithm is a series of script steps in \fBmultiap_extract\fR
which call the model fitting program \fBmodel_fit\fR with different
parameters.  In the script all bands are fit, $x sub 1$ is fixed,
and the asymmetry shape parameter $s sub 2$ is ignored.
The four parameter fit is applied to bands of 32 lines.  The band
solutions are linearly interpolated to the full image and then only
the intensity scale parameter is calculated for each line during the
extraction of the spectra with \fBmodel_extract\fR.
.PP
The general fitting scheme proceeds as follows:
.LP
1.  Fit the three parameters $I sub 0$, $x sub 0$, $s sub 0$ with
$x sub 1$ fixed and $s sub 1$ and $s sub 2$
zero.  This is precisely a Gaussian fit.  The three parameters are
determined simultaneously for all the lines at once using the banded
matrix method.  Thus for 50 lines the solution has 150 variables.
After each fit the scale
parameter $s sub 0$ is smoothed by a polynomial in x.  The polynomial is
taken with seven terms.
.LP
2.  Once the improvement in each iteration becomes less than a
specified amount (2% in rms residual) the next parameter $s sub 1$ is added.
The solution has two steps: fit for $s sub 0$ and $s sub 1$ with $I sub 0$
and $x sub 0$ fixed and
then fit $I sub 0$ and $x sub 0$ with $s sub 0$ and $s sub 1$ fixed.  As before the scale terms
are smoothed by a seventh order polynomial.  Attempts to solve for all
four parameters a once gave unstable results for reasons I don't
understand.
.LP
3.  If desired, the last shape parameter $s sub 2$ can be added by solving
for $s sub 0$, $s sub 1$, and $s sub 2$ while holding $I sub 0$ and
$x sub 0$ fixed and then solving for
$I sub 0$ and $x sub 0$.  This provides some improvement when the signal is very
strong but is generally not needed in the multi-aperture plate data.
It can be an important term in the Echelle data.
.LP
4.  It is possible to then adjust $x sub 1$ followed by steps 2 or 3.
However, this gives very little improvement and $x sub 1$ should be fixed for
each type of data.
.LP
5.  During the final extraction when individual lines are evaluated a one
parameter fit is used to find $I sub 0$ for each spectra.  This is
rather slow, however, on the order of 3 hours per frame.  By using
the pixel value near $x sub 0$ as the value for $I sub 0$ the extraction is reduced
to 13 minutes per frame (see section 12).
.PP
In addition to the preceeding steps the fitting algorithm applies some
heuristic constraints.  These constraints limit how far the peak positions can
shift in each iteration, require the peak intensity to remain positive, and
limit the scale function to be positive at large values of y.
.NH 2
STRIP_EXTRACT -- Unblended Profiles
.PP
For unblended multi-spectra data the profiles can be anything.  The profiles
are obtained by averaging a number of lines (say 32) and normalizing
at some point like the peak value.  These profiles are then used for
detecting bad pixels, such as cosmic rays, and correcting for them as
discussed in the section on cleaning.  Modeling using the \fBmodel_fit\fR
program is only used on Echelle data to find peak positions
accurately in order to follow any curvature of the spectra.
.NH
Extraction and Cleaning
.PP
The extraction of spectra are done separately from the modeling.  It is
possible to extract spectra without any modeling at all using
\fBstrip_extract\fR.  The extraction step also allows the user to specify
if cleaning of the spectra for cosmic rays is desired.  Also modifying
the image is an option.
.NH 2
MODEL_EXTRACT
.PP
Extraction and cleaning using a model fit is described here.
First the $I sub 0$ values for the model profiles are determined for
all the spectra in a line either by multi-parameter fitting or by
taking the peak value.  The pixel values are then compared to the
model in a chi-squared way:

.EQ
r = (data - model) / sqrt model
.EN

If the value of r is larger than a set amount, say 5, then the pixel
value is set to that of the model.  Since the "bad" pixel may affect
the intensity scale $I sub 0$ the cleaning is iterated until no further
pixels are changed.
.PP
The fitting of the data from an individual line of data to the model profiles
is the key element in this scheme.  The best method is to use all the
data in a multi-parameter least square fit.  This minimizes the effect
of bad pixels on the estimated profile which is the essence of this
cleaning method.  However, while the time required to do this for one
line is not great, it adds up to nearly three hours for the 800 lines
in a CRYOCAM frame.  A quick alternative is to scale the model profile
by the data value at the peak position.  This is
quite fast.  However, if the peak has a cosmic ray event or is
otherwise bad then the estimated profile will not correspond closely
to the data profile and the cleaning procedure will make gross errors.
The limited experience I've had with the Echelle and MAP data
has worked well with using the peak estimate.  However, the possible
problems make me nervous and some compromise based on using more than
the peak to estimate the intensity scale of the profile to the data
needs to be found.  This is important because much of the feedback on
the \fBmultispec\fR package from Paul Hintzen and Caty Pilachowski
have dealt with
the particular usefulness of a good cosmic ray cleaning algorithm in
extracting multi-spectra data.
.NH 2
STRIP_EXTRACT
.PP
Removing cosmic rays is the major part of Echelle extraction.
Because these are unblended spectra of arbitrary shape a strip
extraction is all that is needed.
The cleaning is done by the same algorithm used for the multi-aperture
plates except that the profiles are found, as described earlier, by
averaging a number of lines.
The intensity scaling is determined from either a least-square fit
or the peak value.
The same question about the appropriate way to
determine the fit of the profiles to the data discussed previously
applies here except since the spectra are not blended the spectra
can be treated separately in any least square fitting.
.NH
AP_PLATE -- Aperture Plate Correlation
.PP
The final step in the extraction of a multi-aperture plate image is to
correlate the spectra with the on-line database description of the 
drilled hole positions.  This allows for estimates of relative wavelength
offsets and the identification of the spectra with the ID, RA, and DEC
parameters.
.PP
The spectra are fitted to the known aperture plate drilled positions, given in
millimeters, to find an \fIangle\fR for the aperture plate relative to the
detector x dimension and the image \fIscale\fR in pixels / millimeter,

.EQ
x sub fit = a + scale (x sub drill cos (angle) + y sub drill sin (angle))~.
.EN

If the number of spectra is less than that given by the aperture plate drilled
positions then a correlation is done leaving out sequences of
consecutive holes until the fit residual is minimized.  If the number of
spectra is greater than that supposedly drilled then sequences of
consecutive peaks are left out of the fit to minimize the residual.
The missing holes or extra peaks are printed out and, if allowed, the aperture
plate description file is modified, otherwise the program terminates.
In all cases if the final fit residual is greater than 1
pixel the program will terminate.
The program prints out the \fIangle\fR of the aperture plate and the \fIscale\fR
which is also stored in the database.
.PP
An indication that a large part of the focus variation is purely a
scale change is that the derived image \fIscale\fR correlates very well with
the width of the spectra as derived from the profile fitting.  I
estimate that at least 50% of the focus variation is a scale
variation.  This is good news in the sense that a scale variation will
be taken out in the dispersion solution and lines in different parts
of the detector will become more similiar after the solution.
However, the scale variation does not account for all the profile
shape changes and there is indeed a change in the point spread function
across the detector.
.NH
Problems
.PP
There a few problems which I have not been able to resolve or have not
had the time to consider.  The problems which are largely intractable
without a great deal of effort are the unexplained background
variations and deconvolving the spectra for the variation in the
point-spread-function.  The background variations are abrupt increases
in the background in parts of the CRYOCAM detector.  The step edge sometimes
occurs under the spectra and so any smooth polynomial fit to the
background will not be very good.  The modeling of the multi-aperture
plate profiles provides information about the point-spread function
but a deconvolution of the variation in the PSF is a difficult problem
and not warrented at this time.
.PP
I had expected that the large scale response of the CRYOCAM could be
corrected by determining an overall average quartz spectrum from all the
extracted quartz spectra and then dividing the object spectra in each
hole by the ratio of the average quartz spectra from that hole to the
overall average quartz spectrum.  This was attempted and it was found
to work only partially.  Specifically, while there might be a 20%
difference between a spectrum on the edge and one near the center of
the detector the quartz correction left a 10% difference in the object
spectra.  This is apparently due to a poor illumination by the quartz
light source which does not correspond to the telescope illumination.
This quartz correction technique may be made available to users if
desired.
.NH
Comparison with the Cyber Extraction Package
.PP
The discussion of this section must be qualified by the fact that I
have not used the Cyber Extraction Package and I base my understanding on the
algorithms from the Multi-Aperture Plate Data Reduction Manual and
conversations with knowledgable people.  There are many differences
both major and minor and this section only seeks to mention the
some of the important differences.  In the Cyber package:

The background is subtracted from the images as a preliminary process.

The background is either constant or linear across the spectra.

The flat fields are produced by modeling the quartz and data from
several quartz exposures cannot be easily combined.

The initial peak finding and aperture plate correlation algorithm is less
automated in determining missing or additional holes.

The model fitting uses only a three parameter Gaussian model
and the algorithms do not yield results when the focus becomes poor.

The fitting algorithm is neighbor subtraction rather than full
simultaneous solution for all the profiles.

The model fitting is applied only to a quartz and the model is transfered to
object exposures.  This does not allow the shape of the profiles to
change with time as the telescope moves.

The modelling does not couple solutions for neighboring spectra
across the dispersion as is suggested in this design and it does smooth
along the spectra which is not done in this proposal.

The extraction is only to some specified sigma in the model profile and
there is no attempt to correct for blending.

There is no cleaning of the spectra.
.NH
Discussion
.PP
The only data which has passed beyond the extraction phase using the
algorithms described here was that of Paul Hintzen.
Comparison of data reduced by the TV package for
spectra extracted by both the Cyber package and the techniques of the
suggested \fBmultispec\fR package were quite comparable.  To the level he
examined
the spectra there was no clear increase in accuracy though the \fBmultispec\fR
extractions generally had higher counts due to the full extraction of
the blended spectra.  The big advantages found were
the ability to extract all the data even when the focus
became very poor and the success of the cosmic ray cleaning
algorithm.  Thus, Hintzen feels that the need for speed in the extraction
(primarily dependent on the cleaning algorithm)
is modified significantly by the difficulty of dealing with cosmic
rays in the TV spectral analysis programs.  More exploration
of techniques for determining the profile intensity scale from the
model without the full multi-parameter solution is warrented for this
reason.
.PP
I have extracted some Echelle data including field flattening.  The
data had a considerable number of cosmic rays which were removed
quite well.  The extracted spectra were put into a CAMERA format
for further analysis.
.PP
The programs were recently applied to a long slit analysis problem
being studied by Vesa Junkkarinen.  The image was already flat fielded.
The data had two closely spaced and very faint diffuse objects and scattered
light from a nearby QSO.
The three spectra were so weak and closely spaced
that the automatic finding was not used.  However, the rest of the modeling
and extraction were applied directly.
The \fBfind_bckgrnd\fR program, whose original purpose was to correct for
scattered light, worked well to extrapolate the sky across the
image.  The model fitting accurately followed
the peaks of the spectra but the shape fitting was only moderately accurate
since the model shape parameters are not suited to modeling galaxies.
It successfully extracted spectra with a minimum of effort on my part.
Analysis of the extracted spectra and comparison with other techniques
must still be done.  The conclusions to be drawn from this experiment are
that with sufficiently general multi-spectra tools multiple objects in
long slit spectroscopy can be handled.
.PP
One area in which I do not have practical experience is
the extraction of HGVS data.  I believe
the proposed design will work on this type of data.
.PP
A point which needs to be considered in the final design are the
formats of the data files.  The currently used one dimensional spectra
formats are an IIDS format and a CAMERA image format.
The formating of data files for the current spectral analysis packages by
\fBto_iids\fR starts from the \fBmultispec\fR database and throws away a lot 
of information about the spectra.  
Some refinement of this database should focus on the format
to be used by a new \fBIRAF\fR spectral analysis package.
.PP
It should be pointed out that many of the operations can have
alternate algorithms substituted.  In particular, the smoothing
algorithm for the multi-aperture plate flat fields can be replaced by
some other scheme.  The links between the multi-parameter fitting
program and the model have been made very general for investigating
a broad range of models.  Thus, it is also possible to substitute
additional model profiles with relative ease.
.PP
Estimates of excution time are taken from the experimental C programs
implementing the algorithms of this design and they are only
approximate estimates.  The steps corresponding
to \fBdebias\fR, \fBmultispec_flat\fR, and \fBflat_divide\fR for
the multi-aperture data from the CRYOCAM take
about 1 hour for a typical set of frames, say 5 to 15.  This includes
debiasing, triming, computing a flat field from several quartz frames
and dividing the quartz into the object frames.
.PP
The CRYOCAM \fBmultiap_extract\fR phase takes about 40 minutes for the modeling of a frame using 32 lines per band and either 3 hours for an extraction
using the profile fitting
method or 14 minutes for extraction using the peak profile scaling
method.
.PP
Finally, the \fBto_iids\fR takes about 3 minutes per frame.  It takes
this long because it has to convert the \fBmultispec\fR database organized across
the dispersion into formats in which the data is stored as consecutive
spectra; i.e. a type of rotation operation.