Introduction to Spectroscopy

by Andrew Walker
modified by S. Yang


The data reduction consists of several steps:

Pre-Processing

As usual we will be using IRAF to work with the data frames, which should be converted from FITS (a standard format for output from detectors) to IRAF format, using the dataio.rfits task.

First, if you have more than one calibration frame of any given type ( i.e. arc, dark, or flat-field) taken consecutively you should average or median them together (using images.imsum or images.imcombine), to increase the signal-to-noise ratio.

The data is pre-processed by first removing the signature of the detector. A dark frame is subtracted from all of the arc, lamp, and stellar spectra, using images.imarith, and the star-dark spectrum is then divided by the lamp-dark spectrum. For the arc spectrum it is sufficient to simply subtract the dark. You should make sure that the pixtype and calctype are both set to real when using imarith, as integer results will not be useful to you.

Converting to 1-D

If you are using a 2-dimensional detector you will now want to convert your spectra to just 1-dimension. This can be done in the following rather primitive manner (or use the apextract):

You should do this conversion for each of the pre-processed stellar spectra and arcs in the same way.

Continuum fitting

Now that we have a pre-processed 1-D spectrum we can fit the continuum. The task noao.onedspec.continuum has many parameters but nearly all of them can be left unchanged. Other than the input and output images name you will probably only want to alter function and order, which determine the function used to fit the continuum. The best function and order to use will depend on the shape of the continuum, and are best determined by trial and error.
input = input_name
output = output_name
interac = yes
function = spline3
order = 1
While inside the continuum task there are many options available to you to obtain a well fitted continuum. Some of the more useful are listed below:
? - gives help
f - fits the continuum with the currently set function type and order
s - sets sample range
t - resets sample range to full spectrum
z - deletes sample region nearest cursor
q - leaves the continuum task and creates the output file
:function
 - sets the function type where function is one of
{legendre,spline3,spline1,chebyshev}
:order - sets the fitting function order where order is a numeral

Wavelength Calibration

Before your stellar spectrum is really useful you need to calibrate the pixel position as a function of wavelength. This is done using the arc spectra and the noao.onedspec.identify task. Again, this is a complex task with many hidden parameters.

The parameter values given below are only used as examples. You should read through the help page on this task and set appropriate values.
images = image_name
fwidth = 4.0
cradius = 5.0
threshold = 10.0
function = spline3
order = 1
As with continuum there are numerous options inside the identify task, the possibly more useful being:
? - gives help
f - fits the dispersion relation
m - enter wavelength for feature nearest cursor
d - delete feature nearest cursor
z - zoom feature near cursor
r - redraws graph
q - leaves the identify task and saves the fitted dispersion relation
Now that you have determined the dispersion relation for the arc spectrum you have to apply it to the stellar spectrum. First you need to run noao.onedspec.refspec to inform IRAF that the stellar spectrum is to have the same dispersion relation as the arc. You should change the parameters with epar in the following way:
input = input_name
reference = arc_name
ignoreaps = yes
sort =
group =
verbose = yes
You can now tell IRAF to convert your spectrum as a function of pixel number to a spectrum versus wavelength. This is done using noao.onedspec.dispcor with the following parameters set:
input = input_name
output = output_name
database = database
table =
flux = no
ignoreaps = yes
The result should be a wavelength calibrated spectrum of a star. You will need to follow this procedure for each of the stellar spectra that you have taken. There are several ways within IRAF to perform the same task on several images at once, in order to achieve this more efficiently.

Velocity Calculation

For this lab we will calculate only relative radial velocities (i.e. we will determine only the change in R.V. relative to a reference spectrum and not the absolute R.V.). The first thing to do is to select a reference spectrum, which should have a good signal-to-noise (S/N) ratio.

It is then simply a case of taking a cross-correlation between this reference stellar spectrum and every other stellar spectrum that you have taken. This can be done with the noao.rv.fxcor task. This task has numerous parameters, and only the changes from the default values are listed below:
objects = spectrum_name
output = reference_spectrum_name
continuum = none
background = INDEF
There are also an abundance of interactive options available with this task, the more critical of which are listed below:
? - gives help
b - fixes the background level for the Gaussian fit
g - marks correlation peak lag limits and fits
q - quits task
w - writes current correlation results to the log file
y - marks correlation peak lower limit and fits
:function - sets the function type where function is one of
{gaussian,lorentzian,center1d,parabola}
You should read through the help page for the fxcor task carefully and make sure that you understand what it is doing, and how to achieve the most accurate results possible.

Heliocentric Correction

When making spectrographic observations with the intention of calculating radial velocities you have to consider that you are sitting on a moving platform, namely the Earth. The motion of the Earth with respect to the fixed barycentre (or the heliocentre which is essentially equivalent) must be calculated and accounted for. This correction must account for the rotation and orbital motion of the Earth, as well as numerous second order effects for a precise correction.

This can be done using the astutil.rvcorrect task. First you have to let IRAF know which observatory the observations were taken from (as this will effect the heliocentric velocity correction), which is done by typing in the parameters to a file, which could look like this: 

#
# observatory parameters
#
observatory  = "dao"
        name = "Dominion Astrophysical Observatory"
        longitude = 123:25.0
        latitude  =  48:31.2
        altitude  = 238.0
        timezone  = 8
Assuming you called the file obsdb.dat you would then type set obsdb = "obsdb.dat" in IRAF. To use the rvcorrect one of the simpler methods is to generate a file containing the year, month, day, and u.t. of each observation, together with the right ascension and declination of the star with one observation per line. A typical file might look like: 
      1987 10 21 11:00:24  3:36:15   0:22:04
      1987 10 21 11:08:00  8:19:35  -0:51:35
      1987 10 21 11:15:47  8:35:12   6:40:29
      1987 10 21 12:12:10  9:13:20  61:28:49
      1987 10 21 12:16:03  9:27:48   9:07:08
After entering the observatory identification (which is ubc in the example above) and assuming your file to be called rv.obs you would then type rvcorrect f=rv.obs > rv.dat, which would generate a file called rv.dat with all the heliocentric velocity corrections in.

Making the R.V. Curve

All that is now left to do is to apply the heliocentric corrections and plot your data. You should plot the radial velocity versus time, with whatever constitutes your favourite plotting package.