ASTR 4880 Computer Lab Session 8 Version 2

The goals of this session are

  1. Adding header keywords to our stellar spectrum
  2. Calculating the airmass for the stellar observation
  3. Calculating the heliocentric corrections for the stellar observation
  4. Constructing a wavelength calibration for the stellar spectrum

    1. Identifying the lines in the Th-Ar lamp spectrum
    2. Applying the wavelength calibration to the stellar spectrum

Log in to your astro1 account, change to your iraf directory, open an xgterm window, and type cl to start IRAF. Change to the directory where your data from 20091113 are stored. We'll continue working with the same data.

Editing the header of the stellar spectrum

Proceed in the same way as you did for the lamp spectra, but use different keyword and command files. In addition to the data you have already copied, the directory /mnt/vol01/ndm/CCD3/20091113/ now contains files named:

Copy them to your current directory and type page cmds.asth

Unlike the lamp spectra, this stellar one will be archived, and for this reason the full name of the original raw data file is to be included: 200911130017.FIT . Keywords for the star's right ascension and declination and their epoch are provided for. The sidereal time of the observation will be calculated by means of asthedit's built-in function, mst, and the effective airmass of the observation is calculated by means of another built-in function. The effective airmass is a weighted mean over the time of observation; it will differ from, for example, the time at mid-exposure if the exposure is long and the airmass is large, because the airmass is a nonlinear function function of time.

The new keywords file is similar to the one you made for the lamp, but with only one line and with more keywords on that line. Notice that, because the UT of the observation is earlier than 0:00, the observation date is given as 12 November. Because the new keywords file contains additional keywords, we will have to edit the colname parameter for asthedit accordingly. Load the astutil package and type epar asthedit .

PACKAGE = astutil
TASK = asthedit

images = 200911130017.fits Images to be operated upon
commands= cmds.asth File of commands
(table = keywds.tbl) File of values
(colname= fname ut exptime datatype dateobs object ra dec epoch observer) Column names in table
(prompt = asthedit> ) Prompt for STDIN commands
(update = no) Update image header?
(verbose= yes) Verbose output?
(oldstyl= no) Use old style format?
(mode = ql)

Inspect the output for reasonableness. When you are satisfied, type asthedit update+ to enable updating of the image header.

Calculating the heliocentric corrections

Although we often speak of just "the heliocentric correction," in fact there are two corrections, which are calculated by the task rvcorrect. Both carry out a transformation from the moving refence frame of the Earth to the quasi-stationary frame at the center of mass of the solar system. These calculations require the time and date, the latitude and longitude of the observatory, and the right ascension and declination of the star at a specified epoch. For some types of research, these calculations must be highly precise. Procedure
  1. Store the observatory location information. This has to be done only once. This information can be found, for many if not all observatories, in Section J of the Astronomical Almanac. Other methods of storing the data exist, but this is the simplest.
  2. epar observatory

    PACKAGE = noao
    TASK = observatory

    command = set Command (set|list|images)
    obsid = ritter Observatory to set, list, or image default
    images = List of images
    (verbose= yes) Verbose output?

    (observa= obspars) Observatory identification
    (name = Ritter Observatory) Observatory name
    (longitu= 83.61) Observatory longitude (degrees)
    (latitud= 41.65) Observatory latitude (degrees)
    (altitud= 200.) Observatory altitude (meters)
    (timezon= 5.) Observatory time zone
    override= obspars Observatory identification
    (mode = ql)

    In any task requiring observatory information, the parameter value obspars directs IRAF to obtain the parameters from the observatory task.
  3. epar rvcorrect

    PACKAGE = astutil
    TASK = rvcorrect

    (files = ) List of files containing observation data
    (images = 200911130017.fits) List of images containing observation data
    (header = yes) Print header?
    (input = yes) Print input data?
    (imupdat= no) Update image header with corrections?

    (epoch = INDEF) Epoch of observation coordinates (years)
    (observa= obspars) Observatory
    (vsun = 20.) Solar velocity (km/s)
    (ra_vsun= 18.) Right ascension of solar velocity (hours)
    (dec_vsu= 30.) Declination of solar velocity (degrees)
    (epoch_v= 1900.) Epoch of solar coordinates (years)

    (year = ) Year of observation
    (month = ) Month of observation (1-12)
    (day = ) Day of observation
    (ut = ) UT of observation (hours)
    (ra = ) Right ascension of observation (hours)
    (dec = ) Declination of observation (degrees)
    (vobs = 0.) Observed radial velocity
    (hjd = ) Helocentric Julian Day (output)
    (vhelio = ) Helocentric radial velocity (km/s) (output)
    (vlsr = ) Local standard or rest radial velocity (km/s) (output)
    (mode = ql)

    This task permits data entry through a text file (first parameter), through an image header (second parameter, the method demonstrated here), through the third set of task parameters, or by keyboard input (if files = stdin). The solar motion is used to calculate the object's velocity relative to the local standard of rest, if a velocity relative to the Earth is given. In this class, we will not use this feature.

    The image name in the second line should be the most recently processed version of the star image.

    The verbose output specified in the given parameters will result in all the input information being printed. The output will include the value of the VHELIO and HJD keywords. With imupdat = yes these keywords will be added to the image header. When comfortable with the results, type: rvcorrect imupdat+

Extract the Th-Ar spectra

In the night of data we are processing, the stellar observation, number 17, is flanked by Th-Ar spectra, image numbers 16 and 18. In a previous session, you should already have subtracted the median bias from all three images. Now, extract the Th-Ar images with apsum as you did the stellar image, with the flat as a reference and no interactive options set.

PACKAGE = echelle
TASK = apsum

input = 16b.fits List of input images
(output = 16bx.fits) List of output spectra
(apertur= ) Apertures
(format = echelle) Extracted spectra format
(referen= flatredb) List of aperture reference images
(profile= ) List of aperture profile images

(interac= no) Run task interactively?
(find = no) Find apertures?
(recente= no) Recenter apertures?
(resize = no) Resize apertures?
(edit = no) Edit apertures?
(trace = no) Trace apertures?
(fittrac= no) Fit the traced points interactively?
(extract= yes) Extract apertures?
(extras = no) Extract sky, sigma, etc.?
(review = yes) Review extractions?

(line = INDEF) Dispersion line
(nsum = 10) Number of dispersion lines to sum or median

(backgro= fit) Background to subtract (none|average|fit)
(weights= variance) Extraction weights (none|variance)
(pfit = fit1d) Profile fitting type (fit1d|fit2d)
(clean = yes) Detect and replace bad pixels?
(skybox = 1) Box car smoothing length for sky
(saturat= INDEF) Saturation level
(readnoi= 8) Read out noise sigma (photons)
(gain = 6) Photon gain (photons/data number)
(lsigma = 4.) Lower rejection threshold
(usigma = 4.) Upper rejection threshold
(nsubaps= 1) Number of subapertures per aperture
(mode = ql)

Your extracted spectrum should include 4 summed apertures. Printed copies of similar spectra with selected lines identified will be provided, or a larger spectral atlas of the Th-Ar lamp can be found here (PDF, 605 KB). In that atlas, our apertures correspond to apertures 1 - 4. Use the given line identifications as a guide in what follows.


Wavelength calibration (1) Establishing a dispersion solution with the Th-Ar spectrum

This session has been revised from this point onward by NDM.


In the onedspec package: epar identify

PACKAGE = onedspec
TASK = identify

images = 16bx.fits Images containing features to be identified
(section= first line) Section to apply to two dimensional images
(databas= database) Database in which to record feature data
(coordli= linelists$thar.dat) User coordinate list
(units = ) Coordinate units
(nsum = 1) Number of lines/columns/bands to sum in 2D image
(match = 1.) Coordinate list matching limit
(maxfeat= 0) Maximum number of features for automatic identif
(zwidth = 100.) Zoom graph width in user units
(ftype = emission) Feature type
(fwidth = 10.) Feature width in pixels
(cradius= 30.) Centering radius in pixels
(thresho= 100.) Feature threshold for centering
(minsep = 2.) Minimum pixel separation
(functio= chebyshev) Coordinate function
(order = 4) Order of coordinate function
(sample = *) Coordinate sample regions
(niterat= 0) Rejection iterations
(low_rej= 3.) Lower rejection sigma
(high_re= 3.) Upper rejection sigma
(grow = 0.) Rejection growing radius
(autowri= no) Automatically write to database
(graphic= stdgraph) Graphics output device
(cursor = ) Graphics cursor input
crval = Approximate coordinate (at reference pixel)
cdelt = Approximate dispersion
(aidpars= ) Automatic identification algorithm parameters
(mode = ql)

This task uses the interactive curve fitting package to fit a polynomial - here, a cubic - to the feature wavelength as a function of feature position in pixels, in each aperture separately. In using this method, we are sacrificing the information contained in the wavelength relationship among the orders, which is used with ecidentify (session 8.1). However, the last stage of this process does not work when only 4 échelle orders have been calibrated with ecidentify.

Choose an identified line, zoom in on it with the w menu if desired, and type m to mark the feature. A prompt will appear in the yellow bar at the bottom of the graphics window, at which you should type the wavelength given in the printout, to the nearest Å. Do the same for 3 more features. Type f, and check that there are no obvious problems with the fit. With 4 features and 4 polynomical coefficients to determine, the fit should be essentially perfect. Then type q to quit the curve fitting module and move on to the next aperture with k (j takes you to the previous aperture) and repeat.

When you have satisfactory polynomial fits for all the apertures, return to aperture 1 and mark the remaining features that are identified in the printed chart. If your initial fit was correct, the task should provide you with the correct wavelength for each feature, and all you have to do is hit the Enter key to accept it. When you have marked them all, type f to perform a new fit.

Your final rms for each fit should be in the neighborhood of 0.01 Å or less in each aperture. When satisfied, type q and agree to write the results to the database. Your database directory will now contain a text file whose name starts with 'id' and that contains a feature list and the result of the polynomial fitting.

The easiest way to identify the features in the second comparison spectrum is to use the task reidentify, as follows.

PACKAGE = onedspec
TASK = reidentify

referenc= 16bx Reference image
images = 18bx.hhh Images to be reidentified
(interac= yes) Interactive fitting?
(section= first line) Section to apply to two dimensional images
(newaps = no) Reidentify apertures in images not in reference?
(overrid= no) Override previous solutions?
(refit = yes) Refit coordinate function?

(trace = no) Trace reference image?
(step = 10) Step in lines/columns/bands for tracing an image
(nsum = 10) Number of lines/columns/bands to sum
(shift = 0.) Shift to add to reference features (INDEF to search)
(search = 0.) Search radius
(nlost = 0) Maximum number of features which may be lost

(cradius= 30.) Centering radius
(thresho= 0.) Feature threshold for centering
(addfeat= no) Add features from a line list?
(coordli= linelists$thar.dat) User coordinate list
(match = 1.) Coordinate list matching limit
(maxfeat= 0) Maximum number of features for automatic identification
(minsep = 2.) Minimum pixel separation

(databas= database) Database
(logfile= logfile) List of log files
(plotfil= ) Plot file for residuals
(verbose= no) Verbose output?
(graphic= stdgraph) Graphics output device
(cursor = ) Graphics cursor input

answer = yes Fit dispersion function interactively?
crval = Approximate coordinate (at reference pixel)
cdelt = Approximate dispersion
(aidpars= ) Automatic identification algorithm parameters
(mode = ql)

The task will open the first aperture with identified lines marked. Examine each feature to check the identification, and then type f to check the quality of the fit. There could be discrepant points caused by cosmic rays or other problems; check and delete them with d if necessary. When satisfied with the fit in the first aperture, type q and the task will query you (in the text window) whether to move on to the next aperture. (The j and k keys are not available in reidentify.) After the last q, the results will be automatically writeen to the database.

Wavelength calibration (2) Dispersion correction of the stellar spectrum

Unload the onedspec package by typing bye and return to the echelle package, reloading it if necessary.

First, we assign the Th-Ar spectrum as the wavelength reference spectrum. For this purpose, we add a keyword to the image header of the stellar spectrum: refspec1 and give it as a value the name of the comparison lamp spectrum. Use your current name for the stellar spectrum.

PACKAGE = imutil
TASK = hedit

images = 17bxf.fits images to be edited
fields = refspec1 fields to be edited
value = 16bx.hhh value expression
(add = yes) add rather than edit fields
(delete = no) delete rather than edit fields
(verify = yes) verify each edit operation
(show = yes) print record of each edit operation
(update = yes) enable updating of the image header
(mode = ql)

It is possible to use as references two Th-Ar spectra, one taken before and one taken after the stellar spectrum. In that case, the second spectrum will be labeled with the keyword refspec2 in the same way as above. The dispersion correction task will linearly interpolate the dispersion solution between the two comparison spectra. In this session, for simplicity we will use only one comparison spectrum.

Now to apply the dispersion correction:

PACKAGE = echelle
TASK = dispcor

input = 17bxf.fits List of input spectra
output = 17bxfl.fits List of output spectra
(lineari= yes) Linearize (interpolate) spectra?
(databas= database) Dispersion solution database
(table = ) Wavelength table for apertures
(w1 = INDEF) Starting wavelength
(w2 = INDEF) Ending wavelength
(dw = INDEF) Wavelength interval per pixel
(nw = INDEF) Number of output pixels
(log = no) Logarithmic wavelength scale?
(flux = yes) Conserve flux?
(blank = 0.) Output value of points not in input
(samedis= no) Same dispersion in all apertures?
(global = no) Apply global defaults?
(ignorea= no) Ignore apertures?
(confirm= no) Confirm dispersion coordinates?
(listonl= no) List the dispersion coordinates only?
(verbose= yes) Print linear dispersion assignments?
(logfile= ) Log file
(mode = ql)

Plot your wavelength calibrated spectrum and confirm that your calibration is reasonable. Here are the wavelengths of some selected lines in the spectrum of Deneb (alpha Cyg):

6562.817 H alpha
6347.091 Si II
6371.359 Si II

If time permits, we will demonstrate Doppler correction with dopcor. First, with scopy (a task designed to copy wavelength-calibrated spectra), make a fresh copy of your spectrum with a 'd' appended to the file name:

scopy 17tcbx.hhh 17tcbxd.hhh

Then, open the parameter editor for the Doppler correction task. Because of a bug, it does not work for making a new output file, but only for applying the correction in place to a spectrum.

PACKAGE = onedspec
TASK = dopcor

input = 17tcbxd.hhh List of input spectra
output = List of output spectra
redshift= -vhelio Redshift or velocity (Km/s)
(isveloc= yes) Is the redshift parameter a velocity?
(add = no) Add to previous dispersion correction?
(dispers= yes) Apply dispersion correction?
(flux = no) Apply flux correction?
(factor = 3.) Flux correction factor (power of 1+z)
(apertur= ) List of apertures to correct
(verbose= no) Print corrections performed?
(mode = ql)

To confirm that a Doppler shift has been applied, plot aperture 3 of the original spectrum, zoom in on the strong absorption lines, and overplot the Doppler corrected spectrum.

If time permits, continuum normalization will be demonstrated. Otherwise, measure radial velocities and equivalent widths with splot as we have done previously.

Return to ASTR 4880 home page

Latest update (by NDM): 12/9/09