reducing apo dis spectra with iraf

General Tips

The canonical references for reducing long-slit spectra in IRAF are below:

Some general information:

Use the following to get a quick look at the exposure times in an image header:
```
IRAF> imhead [image.fits] long+ | match EXPTIME
```
IRAF multispec (*.ms.fits and resulting dispersion/calibrated/normalized FITS) file image band contents:
- the optimally extracted, sky-subtracted spectrum
- same, but without optimal extraction
- the sky background spectrum
- an estimate of the standard deviation of the first spectrum
If optimal extraction was not used, then the multispec file contains:
- the extracted, sky-subtracted spectrum
- the sky background spectrum

Reducing APO Spectra in IRAF

CCD processing steps

Put all the r and b images in separate directories. This keeps things simple.

Set IRAF to spectrum mode:

IRAF> noao
IRAF> imred
IRAF> ccdred
IRAF> setinstrument specphot

Edit the parameters for ccdproc. Set the following:

biassec = image (to pick up the bias section from the FITS header)
trimsec = (for 2x2 binning, set this to [1:1024,103:481] for R channel,
                                        [27:1023,41:394] for B channel;
           basically get rid of everything not illuminated by flat lamp)
trim = yes
oversca = yes
zerocor = yes

Combine the bias files into a superbias:

IRAF> epar zerocombine
IRAF> zerocombine

At this point, we can also choose to interpolate over bad pixels if they are obviously visible in the flat images. Make a file with the coordinates of the bad pixels to interpolate over:
```
[bad x1] [bad x2] [bad y1] [bad y2]
```

Combine the flats into a superflat:

IRAF> epar flatcombine
IRAF> flatcombine

In ccdproc, set fixpix = yes, and fixfile = [name of bad pixel coordinate file]. Edit the parameters for ccdproc and turn on flat field correction and (if desired) bad pixel interpolation. The task ccdlist is especially helpful for keeping track of what's been done for each frame.

Load the packages twodspec and longslit. Then fit the dispersion function to the superflat to make a normalized superflat:

IRAF> epar response

--> calibrat=[name of combined flat image]
--> normaliz=[name of combined flat image]
--> response=[OUTPUT name of reponse corrected normalized combined flat]

IRAF> response

Use ccdproc to process all files with the generated superbias and superflat, and the trimsec settings. Set zero = [name of superbias] and flat = [name of normalized superflat].
Use imcombine in the sum mode to combine target frames if desirable, and in the median mode to combine comparison arc and standard frames as necessary. An odd number of frames appear to give the best rejection of cosmic ray hits. Use either the crreject or avsigclip methods for rejecting outliers when combining images.

Spectrum extraction

Load twodspec and then apextract to get to the apall routine.

Edit the parameters for apall:

 b_naver = [some negative number to get the median value]
 b_sampl = [the offsets from the center of the aperture to use when fitting background]
 backgro = fit [to subtract the background, set to none otherwise]
 weights = none [for nonoptimal extraction, variance for optimal extraction]
 saturat = INDEF [for nonoptimal, set to saturated DN value for optimal extraction]
 readnoi = rdnoise (tells IRAF to get this from image header keyword RDNOISE)
 gain = gain (tells IRAF to get this from image header keyword GAIN)

Run apall to extract spectra:

IRAF> apall [input fits] output=[output fits]

First extract the aperture by using the interactive fit routines. l sets the lower limit of the aperture, u sets the upper limit of the aperture at the respective cursor locations. b goes to the background fit mode, where t clears the current fit, and using s key in conjunction with the cursor marks the lower/upper limit of background sampling on either side of the aperture. Get the new fit by using the f key. If the background you want isn't visible, use the command: sample -AA:-BB,CC:DD where AA, BB, CC, DD are offsets from the center of the aperture. You can also use w and then a to set the view window to the maximum size.
Then fit the trace to the aperture. :order [x] will set the function order to [x], :function [legendre/spline3] will set the function appropriately. The f key runs the fit and displays it. Hit q to go to the next step.
Review the extracted spectrum.

Wavelength calibration

For each object exposure, prepare a comparison exposure by running apall with identical parameters. We do not subtract the background, and forgo recentering and fitting the trace of the aperture since it must be the same in spatial extent and shape as the target science frame spectrum extraction aperture.
```
IRAF> apall [comparison frame] ref=[object frame] recen- trace- back- intera-
```
Now we need to inspect the comparison arc frame and pick out the wavelengths. DIS has several wavelength plots available at the DIS website. See the plots in the Wavelength Calibration section. It's generally a good idea to download these plots and display them alongside the IRAF windows when carrying out wavelength calibrations. Load the kpnoslit or the twodspec -> longslist packages and then run the identify task:
```
IRAF> epar identify
IRAF> identify
```
For several arcs of the same type, use the reidentify task.
Using the plot above, mark four or five lines using m, and type in their associated wavelengths. Hit f to see the fit, using :funct cheb and :order 2 to fit to a straight line. If the residuals are low, then hit q to return to the wavelength plot, and hit l to automatically identify all the other lines. Try the fit again: f and j to see residuals, then see if a spline fits the data better: :funct spline3 and :order 3. Delete outliers by mousing over them and hitting d. Get the residuals to below 0.1 A if possible (hit f to fit again). Once satisfied, hit q to return to the plot. Hit enter to accept the wavelengths and write them to the database.

Assign this comparison arc to be this object frame's reference spectrum by editing the header for the object frame:

IRAF> hedit [object frame]
fields to be edited (REFSPEC1): [type in REFSPEC1]
value expression (Ne-comp): [type in the comparison arc filename]
add [object frame filename],REFSPEC1 = [comparison arc filename]
update [object frame filename] ? (yes): [hit enter]
[object frame filename] updated

Or more directly:

hedit [object-frame.ms.fits] REFSPEC1 "[arc-comp-frame.ms.fits]" add+

Run the dispersion correction task. This will finally add wavelengths to the x-axis of the spectrum plots:
```
 IRAF> dispcor [name of object frame] output=[name of output object spectrum file]
```

Looking at the spectrum using `splot`

Use splot to look at the final spectrum:

IRAF> splot [name of output object spectrum file]

The following splot commands may be useful:

a = autoexpand between cursor positions
d = deblend several lines by fitting
c = clear and redraw plot
e = equivalent width, integrated flux, and line center between cursor positions
k = profile fit to one line between cursor positions
m = mean, RMS, and SNR in selected region
o = next image will be overplotted
g = get the next image to plot (use in combination with o to overplot two spectra)
s = smooth using boxcar
t = fit continuum (can subtract/normalize, etc)
w = use a custom view window (hit e to mark one corner, then e again to mark the other corner of the window)
( = look at previous spectrum in image
) = look at next spectrum in image
[space] = coordinates and flux at position

When using deblending, mark the continuum first by hitting d on both sides of the lines in question. Then mark the centers of the lines and the line fit type by hitting g (Gaussian), l (Lorentzian), or v (Voigt) at the centers of each of the lines in question. After marking the lines type q to go on to the next step of selecting what to fit. The first question is whether to use the marked positions with f for fixed, fit a single position (the relative separations will be preserved) with s, or fit all the centers independently with a. Use a. The last question is whether to use the continuum entered or fit it simultaneously with the profiles. Hit y to fit the continuum. The deblended fits will then be shown. Use + or - to scroll through the different lines and see their centers, flux values, equivalent widths, and FWHMs.

Flux calibration

First run the setairmass task to correct the observation dates to middle exposure. Edit the parameters and correct the FITS header keywords for APO:
```
equinox = equinox
st = lst
ut = ut1 <-- that's a 'one', not an 'ell'
```
Get the APO extinction curve data file and save it somewhere local (as a .dat text file). Point the kpnoslit/longslit package to it by editing its parameters:
```
extinct = /path/to/apo/extinction/file.dat
```

Run the standard task. Give it the FITS filename of the dispersion corrected spectrophotometric standard you observed. It must be one of the stars known to the twodspec -> longslit or kpnoslit packages. If you can't find your observed standard in the list returned by standard, try loading the ctioslit package and using the standard stars found by the standard task in that package. Execute

ls
onedstds$ctionewcal

in IRAF to see these standard stars more suited for the southern hemisphere (thus ctioslit). The list below is the standard stars from the twodspec -> longslit package.

Standard stars in onedstds$spec50cal/  (3200A - 8100 A)

bd284211   eg247      feige34    hd192281   pg0205134  pg0934554  wolf1346
cygob2no9  eg42       feige66    hd217086   pg0216032  pg0939262
eg139      eg71       feige67    hilt600    pg0310149  pg1121145
eg158      eg81       g191b2b    hz14       pg0823546  pg1545035
eg20       feige110   gd140      hz44       pg0846249  pg1708602

Standard stars in onedstds$spec50cal/  (3200A - 10200A)

bd284211   eg247      feige34    g191b2b    hz44
eg139      eg71       feige66    gd140      pg0823546
eg158      feige110   feige67    hilt600    wolf1346

Fit the sensitivity function by using the sensfunc task. Use the interactive fitting keywords and commands outlined earlier to get the RMS as low as possible. For slit sizes that are well-matched to the seeing, the fit should be very good already.
Now use the calibrate task to apply the flux calibration and extinction corrections to each dispersion corrected spectrum.

Continuum subtraction (optional)

To subtract the continuum from the spectra and normalize them, use the continuum function:

continuum [input spectrum] output=[output spectrum]

To continuum-subtract the sky background as well, use bands=* in the parameters for continuum.