BENCH-MOUNTED ECHELLE REDUCTION MANUAL
nicholas b suntzeff 27 May 1998
The following is a data reduction session for a single night's worth of data taken with the BME in May 1998. The data are of reasonably bright A stars, where the S/N is expected to be about 100. For your data reduction, there may be different specific tasks you will include, but the following reduction is a pretty standard set. This reduction will take your data from raw CCD format into flattened, wavelength calibrated echelle data. Once you are comfortable with the reduction tasks, you will find it takes a full day to reduce one night's worth of data. These reduction notes assume you are well acquainted with IRAF. If you only know IDL, my heart goes out to you. OBSERVATION DETAILS: 1.5m telescope CTIO Bench Mounted Echelle Observer: Carol Grady (NASA Goddard) 17-8 - 22-3 May 1998 31.6 g/mm echelle, tilt 23.5 Gr400 cross disperser, tilt 191.6, centered at 5800A 750mm folded schmidt camera no filters CCD read out in full format (2122,2048) through a single amplifier. Site2K_6 CCD The following chip parameters were measured with "nproto.findgain" on 30 data sections. The CCD info at the telescope gave similar values.nsuntzeff@noao.edusigma biweight N gain 0.94 0.05 0.94 30 (e-/ADU) ron 4.42 0.37 4.37 30 (e-) To calculate this, I used reduced data and ran the following. It is not necessary to reduce the data. Just make sure that the data sections are away from the overscan region. After it runs, edit the "gain" file for wild points, and take an average. cl < fgain > gain fgain: findgain pflat060 pflat061 zero001 zero002 sec=[0100:0125,0100:0125] findgain pflat060 pflat061 zero001 zero002 sec=[0126:0150,0126:0150] findgain pflat060 pflat061 zero001 zero002 sec=[0150:0175,0150:0175] findgain pflat060 pflat061 zero001 zero002 sec=[0176:0199,0176:0199] findgain pflat060 pflat061 zero001 zero002 sec=[0200:0225,0200:0225] findgain pflat060 pflat061 zero001 zero002 sec=[0226:0250,0226:0250] findgain pflat060 pflat061 zero001 zero002 sec=[0250:0275,0250:0275] findgain pflat060 pflat061 zero001 zero002 sec=[0276:0299,0276:0299] findgain pflat060 pflat061 zero001 zero002 sec=[0300:0325,0300:0325] findgain pflat060 pflat061 zero001 zero002 sec=[0326:0350,0326:0350] findgain pflat060 pflat061 zero001 zero002 sec=[0350:0375,0350:0375] findgain pflat060 pflat061 zero001 zero002 sec=[0476:0499,0476:0499] findgain pflat060 pflat061 zero001 zero002 sec=[0400:0425,0400:0425] findgain pflat060 pflat061 zero001 zero002 sec=[0426:0450,0426:0450] findgain pflat060 pflat061 zero001 zero002 sec=[0450:0475,0450:0475] findgain pflat060 pflat061 zero001 zero002 sec=[0476:0499,0476:0499] findgain pflat060 pflat061 zero001 zero002 sec=[0500:0525,0500:0525] findgain pflat060 pflat061 zero001 zero002 sec=[0526:0550,0526:0550] findgain pflat060 pflat061 zero001 zero002 sec=[0550:0575,0550:0575] findgain pflat060 pflat061 zero001 zero002 sec=[0576:0599,0576:0599] findgain pflat060 pflat061 zero001 zero002 sec=[0600:0625,0600:0625] findgain pflat060 pflat061 zero001 zero002 sec=[0626:0650,0626:0650] findgain pflat060 pflat061 zero001 zero002 sec=[0650:0675,0650:0675] findgain pflat060 pflat061 zero001 zero002 sec=[0676:0699,0676:0699] findgain pflat060 pflat061 zero001 zero002 sec=[0700:0725,0700:0725] findgain pflat060 pflat061 zero001 zero002 sec=[0726:0750,0726:0750] findgain pflat060 pflat061 zero001 zero002 sec=[0750:0775,0750:0775] findgain pflat060 pflat061 zero001 zero002 sec=[0776:0799,0776:0799] findgain pflat060 pflat061 zero001 zero002 sec=[0800:0825,0800:0825] findgain pflat060 pflat061 zero001 zero002 sec=[0826:0850,0826:0850] findgain pflat060 pflat061 zero001 zero002 sec=[0850:0875,0850:0875] findgain pflat060 pflat061 zero001 zero002 sec=[0876:0899,0876:0899] findgain pflat060 pflat061 zero001 zero002 sec=[0900:0925,0900:0925] findgain pflat060 pflat061 zero001 zero002 sec=[0926:0950,0926:0950] findgain pflat060 pflat061 zero001 zero002 sec=[0950:0975,0950:0975] findgain pflat060 pflat061 zero001 zero002 sec=[0976:0999,0976:0999] CCD REDUCTIONS Introduction: We pulled off the second night to reduce as a data reduction tutorial. We set up a script file in the local directory which has the basic packages and some definitions. This we called "setup" and ran it everytime we logged in. You could also put these commands in your "loginuser.cl" if you wish. cl < setup setup: noao imred echelle ccdr type /usr/iraf/dev/imtoolrc set stdimage = imt512 set imdir = HDR$pixv2/ set node = ctios3! set raw = /uw53/v2/raw/ keep Some comments about the "setup" file. The word "keep" is necessary otherwise the definitions are only local to the script and will be forgotten on termination of the script. The "noao", "imred", "echelle", and "ccdr" loads these packages. We set the imtool to "imt512" since we were using the 512x512 ximtool format. The command "type /usr/iraf/dev/imtoolrc" types out the full file for your information in case you want to chose another image size. We set "set node = ctios3!" since we were sitting at the "ctios3" computer but using the cpu of ctiow5. If you didn't do this, if you tried to display an image, it would appear on the ctiow5 window. Finally, we set "set raw = /uw53/v2/raw" as a quick way to jump around. We could now move to the directory by merely saying "cd raw$". The basic tasks are in the package "imred" in "noao". We will use the "ccdred" package to reduce the CCD data and the "echelle" package to do the spectral extractions. Note that "ccdred" only works if the CCD was read using one amplifier. If you are using 2 amps, you need to use the "ared.quad" package. If you need help on reducing 2 or 4 amp data, let me know. You may want to run "setinstrument" to set up the basic parameters for the "ccdred" and "ccdproc" tasks. The "subsets" file is not important here. This file is used to keep track of the equivalence of names between filter names in the image header, and the names you really want to call the filters. For the echelle, it can be ignored. ccdred: (pixeltype = "real real") Output and calculation pixel datatypes (verbose = yes) Print log information to the standard output? (logfile = "logfile") Text log file (plotfile = "") Log metacode plot file (backup = "") Backup directory or prefix (instrument = "echccd.dat") CCD instrument file (ssfile = "subsets") Subset translation file (graphics = "stdgraph") Interactive graphics output device (cursor = "") Graphics cursor input (version = "2: October 1987") (mode = "ql") ccdproc: images = "comp*,dark*,obj*,pflat*" List of CCD images to correct (ccdtype = "") CCD image type to correct (max_cache = 0) Maximum image caching memory (in Mbytes) (noproc = no) List processing steps only?\n (fixpix = no) Fix bad CCD lines and columns? (overscan = yes) Apply overscan strip correction? (trim = yes) Trim the image? (zerocor = yes) Apply zero level correction? (darkcor = no) Apply dark count correction? (flatcor = yes) Apply flat field correction? (illumcor = no) Apply illumination correction? (fringecor = no) Apply fringe correction? (readcor = no) Convert zero level image to readout correction? (scancor = no) Convert flat field image to scan correction?\n (readaxis = "line") Read out axis (column|line) (fixfile = "") File describing the bad lines and columns (biassec = "[1:62,390:1800]") Overscan strip image section (trimsec = "[160:1880,390:1800]") Trim data section (zero = "biasred") Zero level calibration image (dark = "") Dark count calibration image (flat = "flat") Flat field images (illum = "") Illumination correction images (fringe = "") Fringe correction images (minreplace = 1.) Minimum flat field value (scantype = "shortscan") Scan type (shortscan|longscan) (nscan = 1) Number of short scan lines\n (interactive = no) Fit overscan interactively? (function = "spline3") Fitting function (order = 5) Number of polynomial terms or spline pieces (sample = "*") Sample points to fit (naverage = 1) Number of sample points to combine (niterate = 4) Number of rejection iterations (low_reject = 2.5) Low sigma rejection factor (high_reject = 2.5) High sigma rejection factor (grow = 0.) Rejection growing radius (mode = "ql") The instrument file is needed to allow the programs to recognize key header information. It was copied from ccddb$ctio/echccd.dat. If you want, you can just point "ccdred" to that file as: ccdred: (instrument = "ccddb$ctio/echccd.dat") CCD instrument file echccd.dat: # ECHCCD.DAT -- Instrument file to be used with ccdred when reducing echelle # spectroscopic data obtained with ArCon. subset none exptime exptime darktime darktime imagetyp imagetyp biassec biassec datasec datasec trimsec trimsec fixfile fixfile OBJECT object DARK dark "PROJECTOR FLAT" flat "SKY FLAT" flat COMPARISON comp ZERO zero # New software BIAS zero # Old software "DOME FLAT" flat MASK other FOCUS object If everything is set up okay, you will get all the useful information using ccdlist: ccdlist *.imh temp3.imh[1721,46][real][object][][OTZF]:HD 163296 thar.ec.imh[1721,46][real][comp][][OTZF]:wavecal Th Ar source zero001.imh[2122,2048][ushort][zero][]:bias frames 5/18 evening zero002.imh[2122,2048][ushort][zero][]:bias frames 5/18 evening etc. If it is not set up correctly, ccdlist will give lots of "unknown"'s. The zero data are raw images written in the fav Arcon integer type called "ushort" or unsigned integer, which can range from 0 to 2^16 (~65000). The "echccd.dat" file tells the "ccdred" program that the image type is located in the FITS image keyword "IMAGETYP" For instance, the image "thar.ec" has the FITS header line: IMAGETYP= 'COMPARISON' / Type of picture (object, dark, etc.) According to the "echccd.dat" table, this converts to "comp" which is a parameter that "ccdred" understands. I often change the "other" to "object" in the instrument file, to allow me to reduce the focus frames like any other frame. Reductions: Do some simple bookkeeping stuff first. Put the JD information in the header (needed for the dispersion correction) and correct the airmass to mid-exposure. Run: observatory (set to ctio) setjd *.imh setairmass *.imh setjd: images = Images (observatory = )_.observatory) Observatory of observation (date = "date-obs") Date of observation keyword (time = "ut") Time of observation keyword (exposure = "exptime") Exposure time keyword (ra = "ra") Right ascension (hours) keyword (dec = "dec") Declination (degrees) keyword (epoch = "epoch") Epoch (years) keyword\n (jd = "jd") Output Julian date keyword (hjd = "hjd") Output Helocentric Julian date keyword (ljd = "ljd") Output local Julian date keyword\n (utdate = yes) Is observation date UT? (uttime = yes) Is observation time UT? (listonly = no) List only without modifying images? (mode = "ql") setairmass: images = "*.imh" Input images (observatory = )_.observatory) Observatory for images (intype = "beginning") Input keyword time stamp (outtype = "effective") Output airmass time stamp\n (date = "date-obs") Observation date keyword (exposure = "exptime") Exposure time keyword (seconds) (airmass = "airmass") Airmass keyword (output) (utmiddle = "utmiddle") Mid-observation UT keyword (output)\n (show = yes) Print the airmasses and mid-UT? (update = yes) Update the image header? (override = yes) Override previous assignments? (mode = "ql") The zeros and the milky flats were combined with "zerocomb" and "flatcomb". It is important at this point to chose carefully the trim section and overscan region using "implot" and also deciding the wavelength coverage. For the x direction (along the echelle dispersion), chose an x trim that brings you down to rougly the same level in the echelle blaze - at about the 30-40% level is fine. Make sure, though, that you are not cutting out important spectral information in the red where the free spectral range (FSR) gets to be the size of the chip. On the other hand, including too large an x range makes later data reduction more difficult for the blue orders where there is little light either in the spectra or the milky flats. You will beging to lose traces and get divide by zero problems if you chose a huge x range. If this happens, the data reduction becomes much less automatic and much more labor intensive. The y-trim is chosen based on how much spectral coverage you want. For our data, we chose the y minimum in the blue where we just began to have signal (about 4500A) and stopped it where there was noticable blue 2nd order leaking into the red first order(about 7500A). This is easy to see as faint "interorder" orders appearing in the red and rapidly getting stronger. Zeros are combined with the task "zerocomb". We combined the raw images which should be good enough rather than reducing the zeros to [OT] and then combining. You probably would not want to combine the raw images if there were large DC variations from image to image, but this not the case with this chip. For the zeros, we used "minmax" rejection in the task "zerocomb" and just threw out the single high pixel to get rid of cr's. We then reduced this combined zero image with "ccdproc" to [OT] and called this "biasred." zerocomb: input = "zero*" List of zero level images to combine (output = "biasred") Output zero level name (combine = "median") Type of combine operation (reject = "minmax") Type of rejection (ccdtype = "zero") CCD image type to combine (process = no) Process images before combining? (delete = no) Delete input images after combining? (clobber = no) Clobber existing output image? (scale = "none") Image scaling (statsec = "") Image section for computing statistics (nlow = 0) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 3.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise = "0.") ccdclip: CCD readout noise (electrons) (gain = "1.") ccdclip: CCD gain (electrons/DN) (snoise = "0.") ccdclip: Sensitivity noise (fraction) (pclip = -0.5) pclip: Percentile clipping parameter (blank = 0.) Value if there are no pixels (mode = "ql") Next we reduced all the milky flat images to [OTZ] with "ccdpr" and combined them using "flatcombine" to an image called "mflat." The flats were scaled using a small section for the scale. We chose the section in a region of the milky flat where the intensity was roughly constant. We did a simple minmax reject throwing out the one high pixel. This is okay since the pflat exposures were short. If they were long (30min) or so, the average may still contain cr's - I would use a median in that case. flatcomb: input = "pflat*" List of flat field images to combine (output = "mflat") Output flat field root name (combine = "average") Type of combine operation (reject = "minmax") Type of rejection (ccdtype = "flat") CCD image type to combine (process = no) Process images before combining? (subsets = no) Combine images by subset parameter? (delete = no) Delete input images after combining? (clobber = no) Clobber existing output image? (scale = "mode") Image scaling (statsec = "[800:900,900:1000]") Image section for computing statistics (nlow = 0) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 3.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise = "0.") ccdclip: CCD readout noise (electrons) (gain = "1.") ccdclip: CCD gain (electrons/DN) (snoise = "0.") ccdclip: Sensitivity noise (fraction) (pclip = -0.5) pclip: Percentile clipping parameter (blank = 1.) Value if there are no pixels (mode = "ql") The image mflat has light everywhere, but the intesity varies tremendously across the image. To make the mflat into a good 2d flat field takes some art. There is no since best way of doing this. The idea is to bring the flat to an average of 1.0 everywhere while leaving the pix-to-pix variations in. We used a simple script as follows: in2: imdel temp* #mflat[1721,1411][real]: milky flats 1e/ADU blkavg mflat temp1 1 1411 blkrep temp1 temp2 1 1411 imarith mflat / temp2 temp3 blkavg temp3 temp4 1721 1 blkrep temp4 temp5 1721 1 imar temp3 / temp5 temp6 fmedian temp6 temp7 xw=51 yw=51 imar temp6 / temp7 flat imdel temp* The script does the following: The image size is [1721,1411]. First we average all the rows together ("blkavg"), create a new image by replicating the averaged row, and then dividing the original by the averaged row. In other words, you are normalizing the data by dividing the image by the data marginalized in the y-direction. You then do the same in the other direction. The image "temp6" is prettly close to flat, with about a 20% variation in the corners. The only thing to watch out for here is that we are dividing the data by single rows and columns that were created by averaging the whole data along one dimension. If there happened to be a bad column or row, this would produce a "bad" pixel in this 1d image. To avoid this, you can fit the 1d image with a curve and use the curve in "blkrep". Use the task "continuum" to do this. To do the final flattening, we create a medianed image with a 51x51 box, and divide the data by that. This final flat is not perfect - you will see a bit of ringing around large "holes" in the CCD chip - but it is pretty good. Note that any residual wiggle that was put into the flat by this technique, will be divided out later if you divide the data by the extracted quartz spectrum. You are now ready to flatten all your data. If you took quartz spectral data during the night, note that it will not be flattened by the 2d image if it its "imagetyp" is "flat." In our case, we had to change the "imagetyp" header from "PROJECTOR FLAT" to "OBJECT" in order to process the quartz images (the ones taken during the night). We used the command "hselect" to do this as: in3: hedit pflat011.imh imagetyp "OBJECT" up+ ver- hedit pflat012.imh imagetyp "OBJECT" up+ ver- hedit pflat013.imh imagetyp "OBJECT" up+ ver- hedit pflat014.imh imagetyp "OBJECT" up+ ver- hedit pflat015.imh imagetyp "OBJECT" up+ ver- hedit pflat023.imh imagetyp "OBJECT" up+ ver- hedit pflat030.imh imagetyp "OBJECT" up+ ver- hedit pflat037.imh imagetyp "OBJECT" up+ ver- hedit pflat044.imh imagetyp "OBJECT" up+ ver- hedit pflat051.imh imagetyp "OBJECT" up+ ver- We reduced all the comp*,obj*,pflat* ( the quartz data) to [OTZF]. DARKS It is good to check the dark count level. The dark signal for these chips is due to scattered light and not a very small structures on the CCD emitting light. Because the dark signal is a very smooth low-frequency signal across the chip, the scattered light correction which we do next will also remove the dark. To reduce the darks, we combined the 3 darks using "combine." combine: input = "dark*" List of images to combine output = "darkn2" List of output images (plfile = "") List of output pixel list files (optional) (sigma = "") List of sigma images (optional)\n (ccdtype = "") CCD image type to combine (optional) (subsets = no) Combine images by subset parameter? (delete = no) Delete input images after combining? (clobber = no) Clobber existing output image?\n (combine = "median") Type of combine operation (reject = "none") Type of rejection (project = no) Project highest dimension of input images? (outtype = "real") Output image pixel datatype (offsets = "none") Input image offsets (masktype = "none") Mask type (maskvalue = 0.) Mask value (blank = 0.) Value if there are no pixels\n (scale = "none") Image scaling (zero = "none") Image zero point offset (weight = "none") Image weights (statsec = "") Image section for computing statistics\n (lthreshold = INDEF) Lower threshold (hthreshold = INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 3.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise = "0.") ccdclip: CCD readout noise (electrons) (gain = "1.") ccdclip: CCD gain (electrons/DN) (snoise = "0.") ccdclip: Sensitivity noise (fraction) (sigscale = 0.1) Tolerance for sigma clipping scaling correction (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0) Radius (pixels) for 1D neighbor rejection (mode = "ql") There was noticeable dark counts at the level of about 1000 ADU in 9000 sec or about 6.6 e-/min. This is rather high, and looks like a light leak in the spectrograph. Should not effect the data too much, but this is too high. 1e-/min is the typical best value we should achieve and 2e-/min is acceptable. SCATTERED LIGHT There is not much scattered light in the spectrograph, but it is easy and rather simple to remove it, so why not? To do this you must first trace the orders to define the region "in between" the orders. Load the echelle package. echelle: (extinction = "onedstds$ctioextinct.dat") Extinction file (caldir = "onedstds$spec16cal/") Standard star calibration directory (observatory = "observatory") Observatory of data (interp = "poly5") Interpolation type (dispaxis = 1) Image axis for 2D/3D images (nsum = "1") Number of lines/columns/bands to sum for 2D/3D (database = "database") Database (verbose = yes) Verbose output? (logfile = "logfile") Text log file (plotfile = "") Plot file\n (records = "") Record number extensions (version = "ECHELLE V3: July 1991") (mode = "ql") We then traced all the object spectra with apall. To run apall, you need to decide a few things as follows. 1. The number of orders. You can quickly figure this out by letting the program find the orders on the quartz spectra. Run "apall" with "nfind=100" or so and find the number of orders. Playing with the "thresh" command makes the aperture search much easier. Set the final number of apertures in the "nfind" parameter. 2. The line width. Find the FWHM of a typical line and use at least 2x this for the "width" parameter. The "width" parameter is really the FW0M. 3. Chose a typical "lower" and "upper" aperture window. You should probably use a fixed lower and upper window size, although this is not required. Run "apall" once on the quartz spectrum and allow the program to adjust the size. The size is adjusted to the "ylevel" value, which is 0.1 or 10% here. Look at all the apertures and take a reasonable average of them. I ususally take the lower and upper values to be roughly the largest of all the apertures at the 10% level. Put these values into the "lower" and "upper" parameters and set "resize=no." 4. Set "nsum" to some negative value like -20 or so. This will take the 20 central columns and median them together (that's what the negative value means). The median will remove the cosmic rays. 5. If you have lots of weak spectra, you may want to set "npeaks=0.5" to take only the strongest 50% of the lines for the shift calculation. 6. Take some sort of reasonable sum for the trace. If you have 2000 pix along the dispersion, taking "t_nsum=20" will mean you are fitting a curve with 100 points. 40 would also be fine here. Generally "t_nsum" = "t_step". Set "t_funct"="spline3" and "t_order"= 1-4 or so. It is important to set "t_niter" = 1 or greater. Don't leave it with the default value of 0. 7. For the traces here, keep clean and variance weighting turned off. Run through a session of "apall" to make sure you are comfortable with the parameters. We did not write out the extracted images at this point. We ran apall as: apall: input = "obj027" List of input images nfind = 46 Number of apertures to be found automatically (output = "") List of output spectra (apertures = "") Apertures (format = "echelle") Extracted spectra format (references = "") List of aperture reference images (profiles = "") List of aperture profile images\n (interactive = yes) Run task interactively? (find = yes) Find apertures? (recenter = yes) Recenter apertures? (resize = no) Resize apertures? (edit = yes) Edit apertures? (trace = yes) Trace apertures? (fittrace = yes) Fit the traced points interactively? (extract = no) Extract spectra? (extras = no) Extract sky, sigma, etc.? (review = no) Review extractions?\n (line = INDEF) Dispersion line (nsum = -20) Number of dispersion lines to sum or median DEFAULT APERTURE PARAMETERS (lower = -5.7) Lower aperture limit relative to center (upper = 4.8) Upper aperture limit relative to center (apidtable = "") Aperture ID table (optional)\n\n# DEFAULT BACKGROUND PARAMETERS (b_function = "chebyshev") Background function (b_order = 1) Background function order (b_sample = "-20:-10,10:20") Background sample regions (b_naverage = -3) Background average or median (b_niterate = 0) Background rejection iterations (b_low_reject = 3.) Background lower rejection sigma (b_high_rejec = 3.) Background upper rejection sigma (b_grow = 0.) Background rejection growing radius\n\n# APERTURE CENTERING PARAMETERS (width = 12.) Profile centering width (radius = 10.) Profile centering radius (threshold = 0.) Detection threshold for profile centering\n\n# AUTOMATIC FINDING AND ORDERING PARAMETERS nfind = 46 Number of apertures to be found automatically (minsep = 10.) Minimum separation between spectra (maxsep = 1000.) Maximum separation between spectra (order = "increasing") Order of apertures\n\n# RECENTERING PARAMETERS (aprecenter = "") Apertures for recentering calculation (npeaks = 30.) Select brightest peaks (shift = yes) Use average shift instead of recentering?\n\n# RESIZING PARAMETERS (llimit = INDEF) Lower aperture limit relative to center (ulimit = INDEF) Upper aperture limit relative to center (ylevel = 0.1) Fraction of peak or intensity for automatic width (peak = yes) Is ylevel a fraction of the peak? (bkg = yes) Subtract background in automatic width? (r_grow = 0.) Grow limits by this factor (avglimits = no) Average limits over all apertures?\n\n# TRACING PARAMETERS (t_nsum = 40) Number of dispersion lines to sum (t_step = 40) Tracing step (t_nlost = 5) Number of consecutive times profile is lost before quit (t_function = "spline3") Trace fitting function (t_order = 3) Trace fitting function order (t_sample = "*") Trace sample regions (t_naverage = 1) Trace average or median (t_niterate = 4) Trace rejection iterations (t_low_reject = 2.5) Trace lower rejection sigma (t_high_rejec = 2.5) Trace upper rejection sigma (t_grow = 0.) Trace rejection growing radius\n\n# EXTRACTION PARAMETERS (background = "none") Background to subtract (skybox = 1) Box car smoothing length for sky (weights = "none") Extraction weights (none|variance) (pfit = "fit1d") Profile fitting type (fit1d|fit2d) (clean = no) Detect and replace bad pixels? (saturation = INDEF) Saturation level (readnoise = "12") Read out noise sigma (photons) (gain = "0.94") Photon gain (photons/data number) (lsigma = 4.) Lower rejection threshold (usigma = 4.) Upper rejection threshold (nsubaps = 1) Number of subapertures per aperture (mode = "ql") We fully traced one bright star (obj038) and called this our reference image. For the rest of the tracing, all the tracing was done relative to this reference image. You can run the tracing for the rest of the images in background if you want, but beware that the traces may get lost in the blue orders for weak S/N spectra. You are advised to do all the traces by hand, at least through the blue orders to make sure the traces are okay. Use the "d" key to delete bad points in the trace, and "a" to force a point where you think the curve should go - if the trace curve gets lost. We looked at all the spectra during "apall" in the lower orders to verify the tracing went well. "Apall" was used as: apall obj039 reference=obj038 We then ran "apscatter" on a test image to get the parameters right. The main parameters to worry about are "buffer" and the fitting parameters for the curves. We found that "buf=1" still left adequate interorder pixels to gauge the scattered light. If you go much further to the blue, this may not work. apscatter: input = "obj038" List of input images to subtract scattered ligh output = "s038" List of output corrected images (apertures = "") Apertures (scatter = "") List of scattered light images (optional) (references = "") List of aperture reference images\n (interactive = yes) Run task interactively? (find = no) Find apertures? (recenter = no) Recenter apertures? (resize = no) Resize apertures? (edit = no) Edit apertures? (trace = no) Trace apertures? (fittrace = no) Fit the traced points interactively? (subtract = yes) Subtract scattered light? (smooth = yes) Smooth scattered light along the dispersion? (fitscatter = yes) Fit scattered light interactively? (fitsmooth = yes) Smooth the scattered light interactively?\n (line = INDEF) Dispersion line (nsum = 10) Number of dispersion lines to sum or median (buffer = 1.) Buffer distance from apertures (apscat1 = "") Fitting parameters across the dispersion (apscat2 = "") Fitting parameters along the dispersion (mode = "ql") Note that apscat1 and apscat2 are also parameter files. You can edit them by the usual way: epar apscat1 etc, or you can also access them from within "epar apscatter" by typing ":e" (which means edit in vi) on the line for the parameter. We then ran the apscatter in script mode: in4: apscatter obj017.imh s017 interactive- apscatter obj018.imh s018 interactive- apscatter obj019.imh s019 interactive- apscatter obj020.imh s020 interactive- apscatter obj021.imh s021 interactive- apscatter obj024.imh s024 interactive- apscatter obj025.imh s025 interactive- apscatter obj026.imh s026 interactive- apscatter obj027.imh s027 interactive- apscatter obj028.imh s028 interactive- apscatter obj031.imh s031 interactive- apscatter obj032.imh s032 interactive- apscatter obj033.imh s033 interactive- apscatter obj034.imh s034 interactive- apscatter obj035.imh s035 interactive- #apscatter obj038.imh s038 interactive- apscatter obj039.imh s039 interactive- apscatter obj040.imh s040 interactive- apscatter obj041.imh s041 interactive- apscatter obj042.imh s042 interactive- apscatter obj045.imh s045 interactive- apscatter obj046.imh s046 interactive- apscatter obj047.imh s047 interactive- apscatter obj048.imh s048 interactive- apscatter obj049.imh s049 interactive- apscatter obj050.imh s050 interactive- EXTRACTIONS These were easy because we have already done the traces. We ran these as background. Note that we did cleaning on the quartz and program objects. With the 2048x2048 format, the cr cleaning works exceedingly well. in5: #quartz flats apall pflat012 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat013 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat014 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat015 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat023 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat030 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat037 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat044 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- apall pflat051 refer=obj038 intera- find- recent+ resiz- edit- trac+ clea+ wei- #th-ar spectra apall comp016 refer=pflat015 intera- find- recent- resize- edit- trac- apall comp022 refer=pflat023 intera- find- recent- resize- edit- trac- apall comp029 refer=pflat030 intera- find- recent- resize- edit- trac- apall comp036 refer=pflat037 intera- find- recent- resize- edit- trac- apall comp043 refer=pflat044 intera- find- recent- resize- edit- trac- apall comp052 refer=pflat051 intera- find- recent- resize- edit- trac- #program spectra apall s017 refer=obj017 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s018 refer=obj018 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s019 refer=obj019 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s020 refer=obj020 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s021 refer=obj021 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s024 refer=obj024 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s025 refer=obj025 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s026 refer=obj026 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s027 refer=obj027 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s028 refer=obj028 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s031 refer=obj031 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s032 refer=obj032 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s033 refer=obj033 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s034 refer=obj034 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s035 refer=obj035 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s038 refer=obj038 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s039 refer=obj039 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s040 refer=obj040 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s041 refer=obj041 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s042 refer=obj042 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s045 refer=obj045 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s046 refer=obj046 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s047 refer=obj047 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s048 refer=obj048 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s049 refer=obj049 intera- find- recen- resiz- edit- trac- clea+ wei=var apall s050 refer=obj050 intera- find- recen- resiz- edit- trac- clea+ wei=var apall: input = "obj027" List of input images nfind = 46 Number of apertures to be found automatically (output = "") List of output spectra (apertures = "") Apertures (format = "echelle") Extracted spectra format (references = "") List of aperture reference images (profiles = "") List of aperture profile images\n (interactive = yes) Run task interactively? (find = yes) Find apertures? (recenter = yes) Recenter apertures? (resize = no) Resize apertures? (edit = yes) Edit apertures? (trace = yes) Trace apertures? (fittrace = yes) Fit the traced points interactively? (extract = yes) Extract spectra? (extras = yes) Extract sky, sigma, etc.? (review = no) Review extractions?\n (line = INDEF) Dispersion line (nsum = -10) Number of dispersion lines to sum or median\n\n# DEFAUL (lower = -5.7) Lower aperture limit relative to center (upper = 4.8) Upper aperture limit relative to center (apidtable = "") Aperture ID table (optional)\n\n# DEFAULT BACKGROUND PA (b_function = "chebyshev") Background function (b_order = 1) Background function order (b_sample = "-20:-10,10:20") Background sample regions (b_naverage = -3) Background average or median (b_niterate = 0) Background rejection iterations (b_low_reject = 3.) Background lower rejection sigma (b_high_rejec = 3.) Background upper rejection sigma (b_grow = 0.) Background rejection growing radius\n\n# APERTURE CENTE (width = 12.) Profile centering width (radius = 10.) Profile centering radius (threshold = 0.) Detection threshold for profile centering\n\n# AUTOMATI (minsep = 10.) Minimum separation between spectra (maxsep = 1000.) Maximum separation between spectra (order = "increasing") Order of apertures\n\n# RECENTERING PARAMETERS\n (aprecenter = "") Apertures for recentering calculation (npeaks = 30.) Select brightest peaks (shift = yes) Use average shift instead of recentering?\n\n# RESIZING (llimit = INDEF) Lower aperture limit relative to center (ulimit = INDEF) Upper aperture limit relative to center (ylevel = 0.1) Fraction of peak or intensity for automatic width (peak = yes) Is ylevel a fraction of the peak? (bkg = yes) Subtract background in automatic width? (r_grow = 0.) Grow limits by this factor (avglimits = no) Average limits over all apertures?\n\n# TRACING PARAMET (t_nsum = 20) Number of dispersion lines to sum (t_step = 20) Tracing step (t_nlost = 3) Number of consecutive times profile is lost before quit (t_function = "spline3") Trace fitting function (t_order = 2) Trace fitting function order (t_sample = "*") Trace sample regions (t_naverage = 1) Trace average or median (t_niterate = 4) Trace rejection iterations (t_low_reject = 2.5) Trace lower rejection sigma (t_high_rejec = 2.5) Trace upper rejection sigma (t_grow = 0.) Trace rejection growing radius\n\n# EXTRACTION PARAMETE (background = "none") Background to subtract (skybox = 1) Box car smoothing length for sky (weights = "none") Extraction weights (none|variance) (pfit = "fit1d") Profile fitting type (fit1d|fit2d) (clean = no) Detect and replace bad pixels? (saturation = INDEF) Saturation level (readnoise = "12") Read out noise sigma (photons) (gain = "0.94") Photon gain (photons/data number) (lsigma = 4.) Lower rejection threshold (usigma = 4.) Upper rejection threshold (nsubaps = 1) Number of subapertures per aperture (mode = "ql") Note that I have increased the read noise to 12. This is because we had an excessive amount of dark (6e-/min/pix). Since your typical exposure was 20 min, this means the effective read noise in the scattered light subtracted image was about 120e-. There is a bug in version 2.11EXPORT for variance extracted data using the echelle package. The WCS headers are not correctly written. The extracted data are okay - you just can't get to band=2 (the aperture extracted spectra) or band=3 (the variance image). If you run the following on all extracted data extracted via variance extraction (in our case, all images with [1721,46,3] ), it will fix the problem. in8; hedit @in9 wcsdim 3 ver- up+ show+ hedit @in9 cd3_3 1. add+ ver- show+ hedit @in9 ltm3_3 1. add+ ver- show+ in9; d012 d013 d014 s050.ec etc. ... Notice that we did the extractions for the object spectra using the "clean" algorithm. To check how well this did, use "splot" to look at a given order. Overplot the band=2 data (you can quickly do this by using the "o" and "%" keys). You will see how many cr's are removed. You may want to check the data for the shifts that were found. Do: grep shift logfile The shifts during the night should be much less than 1 pix. You may see a discontinuity in the shifts at the time of the dewar fill. WAVELENGTH CALIBRATIONS We checked the extracted th-ar spectra for shifts. We chose a subsample of the th-ar spectra that were shifted by less than 1pix. We then combined the comparison spectra into one image. In our case, the first comp image was shifted by about 1/2 FWHM and was not included in the combined image. combine @stuff thar.ec combime: input = "@stuff" List of images to combine output = "thar.ec" List of output images (plfile = "") List of output pixel list files (optional) (sigma = "") List of sigma images (optional)\n (ccdtype = "") CCD image type to combine (optional) (subsets = no) Combine images by subset parameter? (delete = no) Delete input images after combining? (clobber = no) Clobber existing output image?\n (combine = "median") Type of combine operation (reject = "none") Type of rejection (project = no) Project highest dimension of input images? (outtype = "real") Output image pixel datatype (offsets = "none") Input image offsets (masktype = "none") Mask type (maskvalue = 0.) Mask value (blank = 0.) Value if there are no pixels\n (scale = "none") Image scaling (zero = "none") Image zero point offset (weight = "none") Image weights (statsec = "") Image section for computing statistics\n (lthreshold = INDEF) Lower threshold (hthreshold = INDEF) Upper threshold (nlow = 0) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg) (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 3.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise = "0.") ccdclip: CCD readout noise (electrons) (gain = "1.") ccdclip: CCD gain (electrons/DN) (snoise = "0.") ccdclip: Sensitivity noise (fraction) (sigscale = 0.1) Tolerance for sigma clipping scaling corrections (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0) Radius (pixels) for 1D neighbor rejection (mode = "ql") The wavelength solution is done by ecid: images = "comp052.ec" Images containing features to be identified (database = "database") Database in which to record feature data (coordlist = "linelists$thar.dat") User coordinate list (units = "") Coordinate units (match = 0.6) Coordinate list matching limit in user units (maxfeatures = 1000) Maximum number of features for automatic identification (zwidth = 10.) Zoom graph width in user units (ftype = "emission") Feature type (fwidth = 5.) Feature width in pixels (cradius = 5.) Centering radius in pixels (threshold = 10.) Feature threshold for centering (minsep = 2.) Minimum pixel separation (function = "chebyshev") Coordinate function (xorder = 3) Order of coordinate function along dispersion (yorder = 2) Order of coordinate function across dispersion (niterate = 1) Rejection iterations (lowreject = 3.) Lower rejection sigma (highreject = 3.) Upper rejection sigma (autowrite = no) Automatically write to database? (graphics = "stdgraph") Graphics output device (cursor = "") Graphics cursor input (mode = "ql") ecreid: images = "@junk" Spectra to be reidentified reference = "thar.ec" Reference spectrum (shift = 0.) Shift to add to reference features (cradius = 5.) Centering radius (threshold = 10.) Feature threshold for centering (refit = yes) Refit coordinate function? (database = "database") Database (logfiles = "STDOUT,logfile") List of log files (mode = "ql") The wavelength solution only needs to be done for one full echelle spectrum. There are a few tricks to note to get a quick convergence. To start the wavelength solution ,first id a bunch of lines in one order across the dispersion - say 15 lines. Then go to the next order, and id just afew. Set "xor=3" and "yor=2". Do a solution with the "f" key on these two orders. Go to the next order and id a few lines. Do another solution. Now, with the solution from 3 sequential orders, the program should have a good enough solution to id the lines in nearby orders correctly. Now jump 3 orders and check a few lines. Make sure you include lines at the edges of the dispersion also. After jumping 10 orders or so, start increasing the yorder higher - to "yor=3 or 4". For these spectra, we found that "yor=5" was needed for all 46 orders. Also "xor=4 to 6" was the best. It is also important to constantly write out the intermediate solutions using ":w" so you don't lose your work if the program crashes. The program can crash if you have id'ed a line which is not in the table and the program puts in INDEF for the wavelength. This usually does not crash the program but it does sometimes. Beware! To do the final wavelength id's, type "l" for the automatic line id's. With a good prelimiary solution, it will find all the lines up to "maxfeatures." Do the fit, edit out some of the wild points, and adjust the "xor" and "yor" for the final fit. Other things to note: 1. You want to have about 20 lines per order or so. If you have 50 orders this means that you will have 1000 lines in your solution. Make sure that the parameter "maxfeatures=1000" or so. 2. Make sure the "fwidth" parameter is 2x the FWHM. 3. Make sure "niterate" is 1 or greater and not the default 0. 4. After you do the final solution, quickly jump through the orders with the "j" and "k" keys. Did the program id the brighest lines? If it is finding only faint lines and missing the bright ones, you have a lousy solution and must start over. 5. In the fitting section of the program, you can remove the wild points with the "d" key. I would recommend removing them. You can then adjust the "xorder" and "yorder" parameters to figure out what is the minimum order to get the rms down to a mininum. 6. You may want to note down, for your records, the order number for aperture 1 (the blue aperture ususally). In our case it was 123. The echelle equation can be recovered using this value where lambda_c = m/order where m=order(aper 1) * lambda_c (aper 1) and lambda_c is the wavelength in the center of the order. In our case, it is: lambda_c = 563990/order (A) 7. Notice that I have set the "match" parameter to 0.6A. This allows me to enter the line to the nearest integer. If you want to use this solution for another night, in another directory, you can read it in. From within ecid do: :database ../n1/database (you are now taking to the database from the previous night) :read thar.ec :database ../n2/database (go back to the database for your night) a s The "a" says apply the following command to all lines. The "s" key says shift the solution line centers (that were read in the database) to the actual line centers. Once the final solution is done, run the "ecreid" task. This task does a good job at taking the first solution and using this to find the wavelength solutions for other spectra. You now have a solution for the wavelengths on your next night. Typical shifts; ECREIDENTIFY: NOAO/IRAF V2.11EXPORT v2@ctiow5 Wed 11:05:15 27-May-98 Reference image = thar.ec, Refit = yes Image Found Fit Pix Shift User Shift Z Shift RMS comp016.ec 981/981 979/981 -0.952 -4.09 -7.3E-6 0.00261 comp022.ec 980/981 957/980 -0.521 -2.24 -4.0E-6 0.00202 comp029.ec 981/981 954/981 -0.228 -0.981 -1.7E-6 0.00198 comp036.ec 981/981 975/981 -6.54E-4 -0.00283 -5.0E-9 0.0025 comp043.ec 981/981 965/981 0.545 2.34 4.15E-6 0.00218 comp052.ec 981/981 954/981 0.876 3.76 6.67E-6 0.00199 During the night we had shifts of up to 1 pixel, which is large than I normally see. DISPERSION CORRECTIONS First we assign the wavelengths solutions using "refspec". Here we use interpolation of the wavelength solutions using the JD parameter. refspec: input = "s*.ec.imh,pflat*.ec.imh" List of input spectra answer = Accept assignment? (references = "*.imh") List of reference spectra (apertures = "") Input aperture selection list (refaps = "") Reference aperture selection list (ignoreaps = yes) Ignore input and reference apertures? (select = "interp") Selection method for reference spectra (sort = "jd") Sort key (group = "ljd") Group key (time = no) Is sort key a time? (timewrap = 17.) Time wrap point for time sorting (override = no) Override previous assignments? (confirm = no) Confirm reference spectrum assignments? (assign = yes) Assign the reference spectra to the input spectrum? (logfiles = "STDOUT,logfile") List of logfiles (verbose = yes) Verbose log output? (mode = "ql") We used a reference table called "table". Normally you can enter all the comparison data as the second entry on a line, like: pflat012.ec comp016.ec,comp022.ec,comp029.ec,comp036.ec,comp043.ec,comp052.ec pflat013.ec etc. This will default to interpolating the wavelength solutions for a particular program object using the two comparison frames that bracket in time the exposure in question. You may not want to do this though. Generally the dewar is filled during the night, in our case around frame 23. Since this will cause a small jump in the data, it is best to divide the data into two separate interpolation groups: before and after the fill. Thus our table becomes: table: pflat012.ec comp016.ec,comp022.ec pflat013.ec pflat014.ec pflat015.ec s017.ec s018.ec s019.ec s020.ec s021.ec pflat023.ec s024.ec comp029.ec,comp036.ec,comp043.ec,comp052.ec s025.ec s026.ec s027.ec s028.ec pflat030.ec s031.ec s032.ec s033.ec s034.ec s035.ec pflat037.ec s038.ec s039.ec s040.ec s041.ec s042.ec pflat044.ec s045.ec s046.ec s047.ec s048.ec s049.ec s050.ec pflat051.ec Note that if the second item in the table is blank, it defaults to the last second entry. The output from "refspec" looks like the following. It is very important to check this output quickly for obvious assignment errors. Note that if it did not have two comparison spectra that bracketed the exposure in time, it takes the closest one in time. [s017.ec] refspec1='comp016.ec 0.903' [s017.ec] refspec2='comp022.ec 0.097' [s018.ec] refspec1='comp016.ec 0.799' [s018.ec] refspec2='comp022.ec 0.201' [s019.ec] refspec1='comp016.ec 0.599' [s019.ec] refspec2='comp022.ec 0.401' [s020.ec] refspec1='comp016.ec 0.387' [s020.ec] refspec2='comp022.ec 0.613' [s021.ec] refspec1='comp016.ec 0.158' [s021.ec] refspec2='comp022.ec 0.842' [s024.ec] refspec1='comp029.ec' [s025.ec] refspec1='comp029.ec' [s026.ec] refspec1='comp029.ec' [s027.ec] refspec1='comp029.ec' [s028.ec] refspec1='comp029.ec' [pflat044.ec] refspec2='comp052.ec 0.034' [pflat051.ec] refspec1='comp043.ec 0.034' [pflat051.ec] refspec2='comp052.ec 0.966' etc. We ran "dispcor" using "global=yes" defaults. This will put all the data for a given order on the same beginning wavelength and dispersion. The downside of this is that echelle calculates the beginning wavelength in the global solution as the MININUM beginning wavelength for all the data acorss the night. This means that for the most of the data, the first pixel will be 0. Similarly it calculates the MAXIMUM ending wavelength. These zero pixels are a pain if you try to combine spectra across orders or divide the data by the extracted quartz pflats. If you want to avoid this, you have to calculate the MAXIMUM beginning wavelength and the minimum ending wavelength by running "dispcor" with "list+" and "global-". You then will have to specify the wavelength table using these values. dispcor: input = "@in6" List of input spectra output = "@in7" List of output spectra (linearize = yes) Linearize (interpolate) spectra? (database = "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? (samedisp = no) Same dispersion in all apertures? (global = yes) Apply global defaults? (ignoreaps = no) Ignore apertures? (confirm = no) Confirm dispersion coordinates? (listonly = no) List the dispersion coordinates only? (verbose = yes) Print linear dispersion assignments? (logfile = "") Log file (mode = "ql") FINAL FLATTENING This is done to your taste. Do you want to keep the echelle blaze in, which tells you exactly the number of counts (may be needed for EQW chi-sq calculations) or do you want to flatten the spectrum and remove any last vestige of the flat-field wiggles? If you want to merely remove the echelle blaze and retain the rough number of counts, do the following. Here "d051" is an extracted quartz pflat and "d050" is a program object. cont: input = "d051" Input images output = "temp1" Output images ask = "NO" (lines = "*") Image lines to be fit (bands = "1") Image bands to be fit (type = "ratio") Type of output (replace = no) Replace rejected points by fit? (wavescale = yes) Scale the X axis with wavelength? (logscale = no) Take the log (base 10) of both axes? (override = no) Override previously fit lines? (listonly = no) List fit but don't modify any images? (logfiles = "logfile") List of log files (interactive = yes) Set fitting parameters interactively? (sample = "*") Sample points to use in fit (naverage = 1) Number of points in sample averaging (function = "legendre") Fitting function (order = 1) Order of fitting function (low_reject = 4.) Low rejection in sigma of fit (high_reject = 4.) High rejection in sigma of fit (niterate = 1) Number of rejection iterations (grow = 1.) Rejection growing radius in pixels (markrej = yes) Mark rejected points? (graphics = "stdgraph") Graphics output device (cursor = "") Graphics cursor input (mode = "ql") and divide the data as: imar d050 / temp1 c050. c050 is your flattened data. Here you have fit a constant to each order and divided the order by that constant. By dividing your object data by this "normalized" quartz pflat, the number of counts will be approximately conserved and the echelle blaze will be rmoved to first order. Confused as to why you are flat-fielding your data twice? The first flat-fielding (the 2d one in "ccdproc") takes out all the pixel-to-pixel variations in all your data, be it object, comparison, or quartz pflats. It does not perfectly take out lower order wiggles, such as dust "donuts" on the CCD or the echelle blaze. The division of your extracted program data by the extracted quartz data will do this. This "double" flat-fielding is especially robust to spectral shifts. If the spectra shift during the night, this technique will still take out the pixel-to-pixel variations to high order. You can prove this to yourself. Use "splot" and the "m" key to measure the S/N in the extracted quartz pflat. You should find that the S/N is given by sqrt(Ne-) up to S/N of at least 300. On the other hand, you may want to directly continuum flat your data right now without the quartz. If your data is of an absorption line spectrum, do the following: contin; input = "d050" Input images output = "c050"" Output images ask = "NO" (lines = "*") Image lines to be fit (bands = "1") Image bands to be fit (type = "ratio") Type of output (replace = no) Replace rejected points by fit? (wavescale = yes) Scale the X axis with wavelength? (logscale = no) Take the log (base 10) of both axes? (override = no) Override previously fit lines? (listonly = no) List fit but don't modify any images? (logfiles = "logfile") List of log files (interactive = yes) Set fitting parameters interactively? (sample = "*") Sample points to use in fit (naverage = 1) Number of points in sample averaging (function = "spline3") Fitting function (order = 3) Order of fitting function (low_reject = 1.5) Low rejection in sigma of fit (high_reject = 3.) High rejection in sigma of fit (niterate = 3) Number of rejection iterations (grow = 1.) Rejection growing radius in pixels (markrej = yes) Mark rejected points? (graphics = "stdgraph") Graphics output device (cursor = "") Graphics cursor input (mode = "ql") This will exclude the absorption lines and to a robust fit to the continuum. OTHER RANDOM STUFF OF PROBABLY NO INTEREST If you want to get rid of the 0 pixels, try using the task "fixpix" to do so. Make a badpix file with the bad columns marked. badpix: 1 3 1 46 1718 1721 1 46 fixpix: images = "temp3" List of images to be fixed masks = "badpix" List of bad pixel masks (linterp = "INDEF") Mask values for line interpolation (cinterp = "INDEF") Mask values for column interpolation (verbose = yes) Verbose output? (pixels = no) List pixels? (mode = "ql") You can combine a continuum flattened image into a hellishly long single spectrum using the "scombine" command. This is especially useful for emission line objects. Make sure though that you have fixed the 0 pixel problem first. scomb: input = "temp3" List of input spectra output = "temp" List of output spectra (noutput = "") List of output number combined spectra (logfile = "STDOUT") Log file\n (apertures = "") Apertures to combine (group = "all") Grouping option (combine = "average") Type of combine operation (reject = "none") Type of rejection\n (first = no) Use first spectrum for dispersion? (w1 = INDEF) Starting wavelength of output spectra (w2 = INDEF) Ending wavelength of output spectra (dw = INDEF) Wavelength increment of output spectra (nw = INDEF) Length of output spectra (log = no) Logarithmic increments?\n (scale = "none") Image scaling (zero = "none") Image zero point offset (weight = "none") Image weights (sample = "") Wavelength sample regions for statistics\n (lthreshold = INDEF) Lower threshold (hthreshold = INDEF) Upper threshold (nlow = 1) minmax: Number of low pixels to reject (nhigh = 1) minmax: Number of high pixels to reject (nkeep = 1) Minimum to keep (pos) or maximum to reject (neg (mclip = yes) Use median in sigma clipping algorithms? (lsigma = 3.) Lower sigma clipping factor (hsigma = 3.) Upper sigma clipping factor (rdnoise = "0.") ccdclip: CCD readout noise (electrons) (gain = "1.") ccdclip: CCD gain (electrons/DN) (snoise = "0.") ccdclip: Sensitivity noise (fraction) (sigscale = 0.1) Tolerance for sigma clipping scaling correction (pclip = -0.5) pclip: Percentile clipping parameter (grow = 0) Radius (pixels) for 1D neighbor rejection (blank = 0.) Value if there are no pixels (mode = "ql") You can combine a range of apertures as: scomb temp3 temp aper="35-40" For those who are really brave, you can try a flux calibration. Steve Heathcote (sheathcote@noao.edu) has flux tables for the echelle spectrophotometric standards binned every 5A, which are quite useful for this. You can use the "splot" command to look at the regions of spectral overlap. Run "splot" and plot a given order. Then use the "o" and "(" or ")" keys to plot the next higher or lower order.