First decide on the y-axis values over which there is a useful amount of spectrum. If vignetting causes the signal to trail away gradually at one or both ends, then don't try to keep too much of the low-signal bit: this will usually be more difficult to continuum-subtract (or flux-calibrate) and will have a lower signal to noise ratio anyway. On the other hand, it may be useful and can always be ignored later, so you may not want throw away all of the vignetted parts of the spectrum.
The first thing to do is to track the position of the spectrum through the image. This position is needed to define the regions from which the spectrum and the sky are extracted. The basic concept is to calculate precisely where the spectrum lies on the CCD and adjust the regions from which the spectrum and the sky are extracted based on where the spectrum is at each row of the image. For this we must fit a polynomial to the x-axis position of the spectrum on each row of the image and store the result in a file. We use the pamela command track:
track [image (without .sdf)]
This will prompt you for many parameters, because it is a complicated algorithm which can been fooled in many ways by different datasets. I normally call the output file zzz_trpoly (which will result in a file called zzz_trpoly.sdf). The parameter NPOLY (order of the polynomial to fit to the spectrum position) is important: start with 3 but increase this value if the fit is poor. Also set PICK to true unless the algorithm is having difficulty locating the spectrum. The parameter NOBJ exists to select different spectra if there are more than one on each image.
track will plot up the image with the calculated spectrum track on it, and the residuals of the fit. After checking these, look to see that the chi-squared and rms deviation values are reasonable too. Record the x-axis pixel value of the midpoint of the region, which is outputted to screen by track, as this will be needed later.
We need to define the columns from which the spectra and the sky brightnesses should be extracted. For this we use the pamela command regpic:
regpic [image (without .sdf)]
At the prompts, give the balance image (the inverse of the flat field), don't bother with a dark image, give the spectrum distortion trace file (I use zzz_trpoly), and name the sky region output file zzz_sregion (which will create sregion.sdf). Ask for a greyscale plot, give an x-axis range which concentrates on the spectrum to be extracted, and maximise the y-axis range.
Now, with cursor or keyboard (depending on which you chose), define two regions over which the sky should be extracted: one on each side of the spectrum. Both regions should be the same size and reasonably large (so photon noise isn't a big problem) but not too large (because then an overall okay fit might happen to be bad in the small region around the spectrum). At blue optical wavelengths you should try to get quite a large region (60 to 80 pixels each side) to minimise the Poisson noise introduced by the sky fit. At redder wavelengths sky emission lines become problematic and it is best to choose a smaller region (40 to 50 pixels each side). Then define the region over which the spectrum is to be extracted. You should choose the smallest region which covers all the spectrum adequately. Note that the distortions of the spectrum have been taken into account when the regions have been plotted, because we have given the track polynomial file.
Now we must calculate the sky level in the region from which we will extract the spectrum, by fitting polynomials to the counts in the adjacent sky regions using the pamela command skyfit:
skyfit [image] [balance image] mastersky
Give zzz_trpoly as a map of the spectrum distortions. The parameters XSTART and XEND should contain the plot limits you entered previously for regfit. For YSTART and YEND put the limits over which to extract the spectrum (see above). For NPOLY you should give the order of the simplest polynomial which gives a good fit to the sky. Normally 3 will do, but 5 or 7 may be needed sometimes, particularly at redder wavelengths where there are strong night-sky emission lines to fit. Then enter 3 for the sky fit threshold and enter the readout noise you found previously using picstat (see here). PHOTON is the gain of the CCD system and can usually be found in the documentation or the website for the system (or determined from your observations if you took an appropriate set of calibration images).
skyfit now outputs some diagnostics for each hundredth row. Check that the number of rejects is low (0.5% is good, and less than 1% is the maximum, depending on the characteristics of your data).
Now we can actually extract a spectrum from an image. There are three different pamela commands to do this: extnor, extopt and optext, all of which require several parameters (most mentioned previously) to work.
extnor does a simple and robust normal extraction of the spectrum (analagous to aperture photometry). It then reports a noise increase, which is a measure of how much noise the sky fit is adding to the extracted spectrum. A good value for this is close to 1, and you should worry if it is similar to or larger than 2. Changes in the sky regions or polynomial order can decrease or increase this number. The output .sdf file can be plotted using pplot to see the result.
Using Keith Horne's method
The extopt command does an optimal extraction of a spectrum according to the procedures laid down by Keith Horne (1986, PASP, 98, 609). It has now been superseded by Tom Marsh's method for the extraction of highly distorted spectra (see below) and the command is retained in case it may be useful.
This does an optimal extraction, which maximises the signal to noise ratio of the resulting spectrum. NPOLY here asks for the coefficients of y-direction polynomials which provide some compensation for tilts in the spectra. It will need quite a high order (roughly 8) to be useful.
This method was introduced by Tom Marsh (1989, PASP, 101, 1032) to properly extract grating spectra which are significantly tilted on the CCD image. Once track and regpic and skyfit have been run, you will need to use the commands profit and optext to complete the extraction.
The profit command performs a fit to the profile of the spectrum. It runs a series of parallel polynomials along the spectrum (not along columns of pixels) and fits them to the spectrum profile. Somewhat confusingly, the resulting plots are given along the CCD columns, not along the polynomials. Type:
and follow the prompts. I normally set the output file to zzz_fract. A third order polynomial (i.e. set NPOLY to 3) should be enough to properly fit the fractions. Start off with a spacing of 0.5 pixels between polynomials, but it is good to make this figure smaller if you have undersampled data. Make sure you put PLOT to true, so profit will plot up the fit for each column. Monitor these fit plots to ensure that the fits are good. If not, then you may have to increase NPOLY or the pixel polynomial separation.
Now we have everything we need to optimally extract our curved spectrum using the modified algorithm by Tom Marsh, optext:
Enter parameter values at the prompts like you did for extopt. The final result is your (painfully) extracted spectrum (on a pixel scale). Next we can script this to make it all automatic (see next page).