Calibration using IRAF

/!\ To all the students currently working on data reduction: a (new) major revision was done on 29-07-15, 00:00. Explanations are a bit more clear (especially for those which never worked with IRAF or command lines), some tricks are detailed, and the workflow has been modified (bias subtraction has been added). If you are experiencing trouble with your reduced data (results of insufficient quality, etc), it may be worth to have a look here again.

The text in verbatim are IRAF commands used, or their output. Lines starting with "#" are comments. Be sure to insert the correct filenames, those of your files! (Make use of the linux wildcards * and ?) The workflow follows the analysis sequence described in the Quick Analysis Guide. The official documentation for CCD reduction within IRAF is e.g. the guide by Philip Massey.

1. Preparation

In the example below, we assume that all our data is in a directory called "Data", and the calibration files (bias, dark, and flat frames) are in a subdirectory called "Calib", where we want to prepare the master calibration frames.

If you are working on a virtual machine as those provided by the supervisors of the lab course, then it is better to start iraf within the "iraf_work" directory, where a "login.cl" file is located. (As described in the following example). Otherwise, you can go directly to your data directory. In that case, skip the first line of the example.

First, we go to the iraf directory and we start iraf by typing: 

cd iraf_work
cl

We are now within the IRAF environment. Usual shell commands like "cd", "ls" are available to browse files and folders. (You can leave the IRAF environment by typing "logout".) We can now start the reduction of the data. At first, we need to set some parameters and file lists for our convenience. 

   1 # Set the correct image size (4096x4096 pixels):
   2 reset stdimage=imt4096
   3 
   4 # To load the necessary packages type in one after another:
   5 noao
   6 imred
   7 ccdred
   8 # (NOTE: You must type these three commands every time you restart IRAF if you want to do data reduction. IRAF packages are organized in a tree structure, and browsing through packages is like browsing through directories, except you come back to the beginning at each restart.)
   9 
  10 # Go to your calibration folder (replace <calib directory path> with the path of the calibration directory, in this case, it is "Data/Calib")
  11 cd <calib directory path>
  12 
  13 # Create a list of the calibration files (you can display the files directly on the screen by typing "files *fit")
  14 files *fit > calibrationframes.list
  15 
  16 # Check the list
  17 cat calibrationframes.list
  18 
  19 # GENERAL REMARK: the IRAF software does not like names starting with a "-". It is a special character used for options in command lines. If you have image files named "-001G60.fit" or something, rename these files sothat their names start with everything but "-". It will avoid you lots of errors.
  20 
  21 # Create separate lists for flats, bias and darks (assuming the prefixes "DARK", "FLAT" and "BIAS")
  22 cat calibrationframes.list | grep "DARK" > darkframes.list
  23 cat calibrationframes.list | grep "FLAT" > flatframes.list
  24 cat calibrationframes.list | grep "BIAS" > biasframes.list
  25 
  26 # (If necessary, you can edit the lists using any external text editor outside IRAF. Useful if your bias frames are named like "DARK_000s"...)
  27 
  28 # Set image type in IRAF-compliant name (zero, dark, flat, object...) for every calibration frame. These are usually already set, but better checking if they are correct. Otherwise, they will not be compiled by IRAF.
  29 # You can read-out the header using: imheader <name> l+
  30 
  31 # GENERAL REMARK: at any step where filenames are required, you can type them separately (separated with commas, no spaces) or you can provide a list you created using the character "@". In this example, we prefer the latter for convenience.
  32 
  33 
  34 hedit @flatframes.list imagetyp "flat" update+ verify-
  35 hedit @biasframes.list imagetyp "zero" update+ verify-
  36 hedit @darkframes.list imagetyp "dark" update+ verify-
  37 
  38 # Check image types
  39 ccdlist @calibrationframes.list ccdtype=""

We here have a bifurcation in the workflow according to what kind of calibration frames you have. If you have dark frames for all of your observations (i.e., calibration frames with exactly the same exposure time of every frame), then you can apply the direct master dark subtraction: bias and dark current will be accurately subtracted, leaving no flaw. If you don't (for instance, if you did not take any calibration for your flat fields), then you have to apply separate bias and dark current subtractions in order to achieve the same quality of your final results.

If you can offer to perform direct dark subtraction, please refer to Section 2. If you need to perform separate bias and dark subtractions, please refer to Sections 2' and 3. (We describe both procedures with the most general case as an example, i.e., a set of scaleable dark frames. you can easily adapt the reduction sequence for every exposure time you have.)

Afterwards, go back to the standard data reduction sequence, from Section 4.

2. Direct master dark subtraction: prepare the master darks

First, create several dark frame lists according to exposure time and camera temperature (after having checked what you had). Name those lists with self-explanatory names. Then, combine the dark frames with the task "darkcombine". 

   1 # Create the dark lists. Here is a first manual way to do it.
   2 
   3 files <insert filenames of darks with 1s> > dark_1s_-25deg.list
   4 files <insert filenames of darks with 2s> > dark_2s_-25deg.list
   5 files <insert filenames of darks with 3s> > dark_3s_-25deg.list
   6 files <insert filenames of darks with 5s> > dark_5s_-25deg.list
   7 files <insert filenames of darks with 10s> > dark_10s_-25deg.list
   8 files <insert filenames of darks with 20s> > dark_20s_-25deg.list
   9 files <insert filenames of darks with 60s> > dark_60s_-25deg.list
  10 files <insert filenames of darks with 120s> > dark_120s_-25deg.list
  11 files <insert filenames of darks with 180s> > dark_180s_-25deg.list
  12 files <insert filenames of darks with 240s> > dark_240s_-25deg.list
  13 
  14 
  15 # Another, more automatic way (assuming the filenames contain information about exposure time):
  16 
  17 cat darkframes.list | grep "001s" > dark_1s_-25deg.list
  18 cat darkframes.list | grep "002s" > dark_2s_-25deg.list
  19 cat darkframes.list | grep "003s" > dark_3s_-25deg.list
  20 cat darkframes.list | grep "005s" > dark_5s_-25deg.list
  21 cat darkframes.list | grep "010s" > dark_10s_-25deg.list
  22 cat darkframes.list | grep "020s" > dark_20s_-25deg.list
  23 cat darkframes.list | grep "060s" > dark_60s_-25deg.list
  24 cat darkframes.list | grep "120s" > dark_120s_-25deg.list
  25 cat darkframes.list | grep "180s" > dark_180s_-25deg.list
  26 cat darkframes.list | grep "240s" > dark_240s_-25deg.list
  27 
  28 
  29 # Set the parameters of darkcombine:
  30 epar darkcombine
  31 
  32 # You are now in an interactive interface showing a list of all the parameters for the task "darkcombine". You can navigate using the arrow keys, type and press enter to edit a parameter value (error messages are shown at the bottom). Set the parameter list to the following:
  33 
  34 PACKAGE = ccdred
  35    TASK = darkcombine
  36 
  37 input   =                     ) List of dark images to combine
  38 (output =                     ) Output dark image root name
  39 (combine=               median) Type of combine operation
  40 (reject =              sigclip) Type of rejection
  41 (ccdtype=                 dark) CCD image type to combine
  42 (process=                   no) Process images before combining?
  43 (delete =                   no) Delete input images after combining?
  44 (clobber=                   no) Clobber existing output image?
  45 (scale  =             exposure) Image scaling
  46 (statsec=                     ) Image section for computing statistics
  47 (nlow   =                    0) minmax: Number of low pixels to reject
  48 (nhigh  =                    1) minmax: Number of high pixels to reject
  49 (nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
  50 (mclip  =                  yes) Use median in sigma clipping algorithms?
  51 (lsigma =                   3.) Lower sigma clipping factor
  52 (hsigma =                   3.) Upper sigma clipping factor
  53 (rdnoise=                   0.) ccdclip: CCD readout noise (electrons)
  54 (gain   =                   1.) ccdclip: CCD gain (electrons/DN)
  55 (snoise =                   0.) ccdclip: Sensitivity noise (fraction)
  56 (pclip  =                 -0.5) pclip: Percentile clipping parameter
  57 (blank  =                   0.) Value if there are no pixels
  58 (mode   =                   ql)
  59 
  60 # Exit the interface by typing ":q".
  61 
  62 # Run darkcombine to combine dark frames of equal exposure times into master darks. Choose self-explanatory names for the resulting, combined frames.
  63 # NOTE: If some lists have only a few dark frames (3 or less), use "average" instead of "median" for the combine parameter. (In this example, the longest exposure times have the least dark frames.)
  64 
  65 darkcombine @dark_1s_-25deg.list output="master_dark_1s_-25deg" combine="median" &
  66 darkcombine @dark_2s_-25deg.list output="master_dark_2s_-25deg" combine="median" &
  67 darkcombine @dark_3s_-25deg.list output="master_dark_3s_-25deg" combine="median" &
  68 darkcombine @dark_5s_-25deg.list output="master_dark_5s_-25deg" combine="median" &
  69 darkcombine @dark_10s_-25deg.list output="master_dark_10s_-25deg" combine="median" &
  70 darkcombine @dark_20s_-25deg.list output="master_dark_20s_-25deg" combine="median" &
  71 darkcombine @dark_60s_-25deg.list output="master_dark_60s_-25deg" combine="median" &
  72 darkcombine @dark_120s_-25deg.list output="master_dark_120s_-25deg" combine="median" &
  73 darkcombine @dark_180s_-25deg.list output="master_dark_180s_-25deg" combine="average" &
  74 darkcombine @dark_240s_-25deg.list output="master_dark_240s_-25deg" combine="average" &
  75 
  76 # If some individual frames are giving you trouble (error messages), you can exclude them from the lists.
  77 
  78 #NOTE: the "&" parameter is only meant to run the task in the background. It is not mandatory.
  79

Now jump on Section 4 to continue with the standard reduction sequence.

2'. Separate subtraction: prepare the master bias frame

If you didn't, create a list containing the names of all the bias frames you have (as explained in Section 1). Then, combine the bias frames with the task "zerocombine". 

   1 # Set the parameters of zerocombine:
   2 epar zerocombine
   3 
   4 # You are now in an interactive interface showing a list of all the parameters for the task "zerocombine". You can navigate using the arrow keys, type and press enter to edit a parameter value (error messages are shown at the bottom). Set the parameter list to the following:
   5 
   6 PACKAGE = ccdred
   7    TASK = zerocombine
   8 
   9 input   =                       List of zero level images to combine
  10 (output =                     ) Output zero level name
  11 (combine=               median) Type of combine operation
  12 (reject =              sigclip) Type of rejection
  13 (ccdtype=                 zero) CCD image type to combine
  14 (process=                   no) Process images before combining?
  15 (delete =                   no) Delete input images after combining?
  16 (clobber=                   no) Clobber existing output image?
  17 (scale  =                 none) Image scaling
  18 (statsec=                     ) Image section for computing statistics
  19 (nlow   =                    0) minmax: Number of low pixels to reject
  20 (nhigh  =                    1) minmax: Number of high pixels to reject
  21 (nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
  22 (mclip  =                  yes) Use median in sigma clipping algorithms?
  23 (lsigma =                   3.) Lower sigma clipping factor
  24 (hsigma =                   3.) Upper sigma clipping factor
  25 (rdnoise=                   0.) ccdclip: CCD readout noise (electrons)
  26 (gain   =                   1.) ccdclip: CCD gain (electrons/DN)
  27 (snoise =                   0.) ccdclip: Sensitivity noise (fraction)
  28 (pclip  =                 -0.5) pclip: Percentile clipping parameter
  29 (blank  =                   0.) Value if there are no pixels
  30 (mode   =                   ql)
  31 
  32 # Exit the interface by typing ":q".
  33 
  34 # Run zerocombine to combine your bias frames into a master bias. Give it a self-explanatory name as output.
  35 # NOTE: If you have 3 bias frames or less, use "average" instead of "median" for the combine parameter. You can set it directly in the zerocombine parameters or specify it in the command line for one run (as examplified below).
  36 
  37 zerocombine @biasframes.list output="master_bias" combine="average" &
  38 
  39 # If some individual frames are giving you trouble (error messages), you can exclude them from the list.
  40 
  41 #NOTE: the "&" parameter is only meant to run the task in the background. It is not mandatory.
  42

Now you have your master bias, picturing the bias signal caused by the voltage of the CCD detector during the readout. Because it affects every image taken with the camera, we need to subtract it from absolutely every frame: darks, flats and scientific frames. In this example, we perform this task separately in different Sections for clarity.

The subtraction is performed by the task imarith (it does not need any parameter settings, you can use it as is). 

   1 # Do subtraction. The syntax for imarith is: imarith [frame_1] [operation] [frame_2] [frame_result]. You can provide lists instead of individual frames.
   2 # In this example, we subtract a whole list by a single frame, and we save the results in separate files by giving them a suffix. (This is done by providing a filelist plus a suffix as an output. You can also provide a new list you made containing custom filenames.) 
   3 
   4 imarith @darkframes.list - master_bias.fits  bias_sub_@darkframes.list &

3. Separate subtraction: prepare the master dark frames

First, create several dark frame lists (with bias-subtracted darks only) according to exposure time and camera temperature (after having checked what you had). Name those lists with self-explanatory names. Then, combine the dark frames with the task "darkcombine". 

   1 # Create the dark lists. Here is a first manual way to do it.
   2 
   3 files <insert filenames of bias-subtracted darks with 1s> > dark_1s_-25deg.list
   4 files <insert filenames of bias-subtracted darks with 2s> > dark_2s_-25deg.list
   5 files <insert filenames of bias-subtracted darks with 3s> > dark_3s_-25deg.list
   6 files <insert filenames of bias-subtracted darks with 5s> > dark_5s_-25deg.list
   7 files <insert filenames of bias-subtracted darks with 10s> > dark_10s_-25deg.list
   8 files <insert filenames of bias-subtracted darks with 20s> > dark_20s_-25deg.list
   9 files <insert filenames of bias-subtracted darks with 60s> > dark_60s_-25deg.list
  10 files <insert filenames of bias-subtracted darks with 120s> > dark_120s_-25deg.list
  11 files <insert filenames of bias-subtracted darks with 180s> > dark_180s_-25deg.list
  12 files <insert filenames of bias-subtracted darks with 240s> > dark_240s_-25deg.list
  13 
  14 
  15 # Another, more automatic way (assuming the filenames contain information about exposure time):
  16 
  17 cat darkframes.list | grep "bias_sub_" | grep "001s" > dark_1s_-25deg.list
  18 cat darkframes.list | grep "bias_sub_" | grep "002s" > dark_2s_-25deg.list
  19 cat darkframes.list | grep "bias_sub_" | grep "003s" > dark_3s_-25deg.list
  20 cat darkframes.list | grep "bias_sub_" | grep "005s" > dark_5s_-25deg.list
  21 cat darkframes.list | grep "bias_sub_" | grep "020s" > dark_20s_-25deg.list
  22 cat darkframes.list | grep "bias_sub_" | grep "060s" > dark_60s_-25deg.list
  23 cat darkframes.list | grep "bias_sub_" | grep "120s" > dark_120s_-25deg.list
  24 cat darkframes.list | grep "bias_sub_" | grep "180s" > dark_180s_-25deg.list
  25 cat darkframes.list | grep "bias_sub_" | grep "240s" > dark_240s_-25deg.list
  26 
  27 
  28 # Set the parameters of darkcombine to the following:
  29 epar darkcombine
  30 
  31 PACKAGE = ccdred
  32    TASK = darkcombine
  33 
  34 input   =                     ) List of dark images to combine
  35 (output =                     ) Output dark image root name
  36 (combine=               median) Type of combine operation
  37 (reject =              sigclip) Type of rejection
  38 (ccdtype=                 dark) CCD image type to combine
  39 (process=                   no) Process images before combining?
  40 (delete =                   no) Delete input images after combining?
  41 (clobber=                   no) Clobber existing output image?
  42 (scale  =             exposure) Image scaling
  43 (statsec=                     ) Image section for computing statistics
  44 (nlow   =                    0) minmax: Number of low pixels to reject
  45 (nhigh  =                    1) minmax: Number of high pixels to reject
  46 (nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
  47 (mclip  =                  yes) Use median in sigma clipping algorithms?
  48 (lsigma =                   3.) Lower sigma clipping factor
  49 (hsigma =                   3.) Upper sigma clipping factor
  50 (rdnoise=                   0.) ccdclip: CCD readout noise (electrons)
  51 (gain   =                   1.) ccdclip: CCD gain (electrons/DN)
  52 (snoise =                   0.) ccdclip: Sensitivity noise (fraction)
  53 (pclip  =                 -0.5) pclip: Percentile clipping parameter
  54 (blank  =                   0.) Value if there are no pixels
  55 (mode   =                   ql)
  56 
  57 
  58 # Run darkcombine to combine dark frames of equal exposure times into master darks. Choose self-explanatory names for the resulting, combined frames. Again, use 'average' instead of 'median' if some lists have 3 exposures or less (which is generally the case for longer exposure times). In this example, we don't have enough frames for 3min and 4min exposure time:
  59 
  60 darkcombine @dark_1s_-25deg.list output="master_dark_1s_-25deg" combine="median" &
  61 darkcombine @dark_2s_-25deg.list output="master_dark_2s_-25deg" combine="median" &
  62 darkcombine @dark_3s_-25deg.list output="master_dark_3s_-25deg" combine="median" &
  63 darkcombine @dark_5s_-25deg.list output="master_dark_5s_-25deg" combine="median" &
  64 darkcombine @dark_10s_-25deg.list output="master_dark_10s_-25deg" combine="median" &
  65 darkcombine @dark_20s_-25deg.list output="master_dark_20s_-25deg" combine="median" &
  66 darkcombine @dark_60s_-25deg.list output="master_dark_60s_-25deg" combine="median" &
  67 darkcombine @dark_120s_-25deg.list output="master_dark_120s_-25deg" combine="median" &
  68 darkcombine @dark_180s_-25deg.list output="master_dark_180s_-25deg" combine="average" &
  69 darkcombine @dark_240s_-25deg.list output="master_dark_240s_-25deg" combine="average" &

You now have a master bias frame (containing only the bias noise) and several master dark frames (containing only dark current). You may now combine those frames to create new ones with custom exposure times using imarith. We explain here how to proceed. (Let's suppose we want to make a master dark frame with an exposure time of 1.4 seconds: 

   1 # First, locate the closest master dark (according to exposure time) among your master dark set. It also needs to have a LONGER exposure time. In this case, our closest (longer) dark frame is 2s. We then use the imarith command to scale the dark current down to 1.4s:
   2 
   3 imarith master_dark_2s_-25deg.fits * 0.7 master_dark_1s4_-25deg.fits &
   4 
   5 # (to know the converson factor to apply, use your mind or your calculator.)
   6

Proceed the same way to produce master dark frames with any exposure time you need.

The reason why you must always use a master dark frame with a longer exposure time to create a new, shorter one is the risk of noise amplification: when you increase the dark current by multiplying its intensity a constant value, the standard deviation is increased the same, and you artifically produce a larger noise than it could be in reality. This is also true for the other way around: reducing the dark current intensity decreases its scatter and weakens the noise, which is less damaging since a lower noise will be hidden behind other, larger kinds of noise anyway. So NEVER mutiply the dark current by a >1 constant value!

If you have no dark frame with a longer exposure time that the one you need, then it is preferable to sum up different frames of shorter exposure times with imarith. In this example, let's suppose we need an exposure time of 6 minutes (300s); the best would be to add the 240s frame to the 60s frame: 

   1 imarith master_dark_240s_-25deg.fits + master_dark_60s_-25deg.fits  master_dark_300s_-25deg.fits &

4. Prepare the master flat-fields

First, subtract the master dark frame having an exposure time similar to that of each flat-field frame (the exposure time is listed in the header, check it with imheader as before). Subtract the master bias frame too if you have one. The task used is imarith. 

   1 # If you have separate masterdarks and masterbias (in case you didn't perform direct dark subtraction), subtract the masterbias from all your flats first. (Skip this step otherwise.)
   2 imarith @flatframes.list - master_bias.fits bias_sub_@flatframes.list &
   3 
   4 # Now prepare list of files as explained for the masterdarks. (If you have performed the masterbias subtraction at the previous step, these lists must contain subtracted flats only.)
   5 files <insert filenames of flats in red filter> > flatred.list
   6 files <insert filenames of flats in green filter> > flatgreen.list
   7 files <insert filenames of flats in blue filter> > flatblue.list
   8 files <insert filenames of flats in Halpha filter> > flathalpha.list
   9 files <insert filenames of flats in luminance filter> > flatluminance.list
  10 
  11 # Do the masterdark subtraction separately for each filter and each exposure time.
  12 imarith @flatred.list - master_dark_1s_-25deg.fits  dark_sub_@flatred.list &
  13 imarith @flatblue.list - master_dark_2s5_-25deg.fits  dark_sub_@flatblue.list &
  14 imarith @flatluminance.list - master_dark_0s2_-25deg.fits  dark_sub_@flatluminance.list &
  15 imarith @flatgreen.list - master_dark_1s4_-25deg.fits  dark_sub_@flatgreen.list &

Then, combine the dark-subtracted flat-fields separately for each filter. We do not normalise the flat field by the mean value now (the flats are scaled by their mean value when used by ccdproc later). 

   1 # Edit parameters of flatcombine as the following:
   2 epar flatcombine
   3 
   4 input   =                       List of flat field images to combine
   5 (output =                     ) Output flat field root name
   6 (combine=              average) Type of combine operation
   7 (reject =              sigclip) Type of rejection
   8 (ccdtype=                 flat) CCD image type to combine
   9 (process=                   no) Process images before combining?
  10 (subsets=                   no) Combine images by subset parameter?
  11 (delete =                   no) Delete input images after combining?
  12 (clobber=                   no) Clobber existing output image?
  13 (scale  =                 mode) Image scaling
  14 (statsec=                     ) Image section for computing statistics
  15 (nlow   =                    1) minmax: Number of low pixels to reject
  16 (nhigh  =                    1) minmax: Number of high pixels to reject
  17 (nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
  18 (mclip  =                  yes) Use median in sigma clipping algorithms?
  19 (lsigma =                   3.) Lower sigma clipping factor
  20 (hsigma =                   3.) Upper sigma clipping factor
  21 (rdnoise=                   0.) ccdclip: CCD readout noise (electrons)
  22 (gain   =                   1.) ccdclip: CCD gain (electrons/DN)
  23 (snoise =                   0.) ccdclip: Sensitivity noise (fraction)
  24 (pclip  =                 -0.5) pclip: Percentile clipping parameter
  25 (blank  =                   1.) Value if there are no pixels
  26 (mode   =                   ql)
  27 
  28 # Run flatcombine to produce the master flat:
  29 flatcombine dark_sub_@flatred.list output=masterflat_red.fits &
  30 
  31 # Here is an example of the output of flatcombine (in the "logfile" file created in the current directory):
  32 May  8 16:57: IMCOMBINE
  33   combine = average, scale = mode, zero = none, weight = none
  34   reject = sigclip, mclip = yes, nkeep = 1
  35   lsigma = 3., hsigma = 3.
  36   blank = 1.
  37                 Images     Exp    Mode  Scale
  38   darksub_AutoFlat-PANoRot-Red-Bin1-001-NoRot.fits    3.7  16360.  1.035
  39   darksub_AutoFlat-PANoRot-Red-Bin1-002-NoRot.fits    4.1  17095.  0.990
  40   darksub_AutoFlat-PANoRot-Red-Bin1-003-NoRot.fits    4.4  17113.  0.989
  41   darksub_AutoFlat-PANoRot-Red-Bin1-004-NoRot.fits    4.7  17089.  0.991
  42   darksub_AutoFlat-PANoRot-Red-Bin1-005-NoRot.fits    5.0  17001.  0.996
  43 
  44   Output image = combineflat_red.fits, ncombine = 5
  45 
  46 # Proceed similarly for flats in other filters.
  47 
  48 flatcombine dark_sub_@flatgreen.list output=masterflat_green.fits &
  49 flatcombine dark_sub_@flatblue.list output=masterflat_blue.fits &
  50 flatcombine dark_sub_@flatluminance.list output=masterflat_luminance.fits &

5. Process the raw data

We now go to a directory with the raw data of an object. In that example we have 10*180 second exposures of M53, in the R, G, and B filters. 

   1 # Go to the object directory
   2 cd ../M53
   3 
   4 # Display the files on the screen
   5 ls
   6 
   7 #Example of output:
   8 M53_-003_G.fit         M53_-006_B.fit  M53_-008_R.fit
   9 M53_-001_B.fit         M53_-003_R.fit  M53_-006_G.fit  M53_-009_B.fit
  10 M53_-001_G.fit         M53_-004_B.fit  M53_-006_R.fit  M53_-009_G.fit
  11 M53_-001_R.fit         M53_-004_G.fit  M53_-007_B.fit  M53_-009_R.fit
  12 M53_-002_B.fit         M53_-004_R.fit  M53_-007_G.fit  M53_-010_B.fit
  13 M53_-002_G.fit         M53_-005_B.fit  M53_-007_R.fit  M53_-010_G.fit
  14 M53_-002_R.fit         M53_-005_G.fit  M53_-008_B.fit  M53_-010_R.fit
  15 M53_-003_B.fit         M53_-005_R.fit  M53_-008_G.fit
  16 
  17 # make list of files (all of them, and filter-separated):
  18 files M53_-*.fit > M53-all.list
  19 files M53_*R.fit > M53-R.list
  20 files M53_*G.fit > M53-G.list
  21 files M53_*B.fit > M53-B.list
  22 
  23 # Set image type in IRAF-compliant name (object) for each frame.
  24 hedit @M53-all.list imagetyp "object" update+ verify-
  25 
  26 # Make a copy of your raw data, because some task (like ccdproc) will overwrite the input files. (Here, we give a "raw" prefix to all of them.)
  27 imcopy @M53-all.list raw//@M53-all.list &
  28 
  29 # Copy the necessary calibration files (master dark frame and master flats) into the current directory.
  30 imcopy ../Calibration/master*.fits . &
  31 
  32 # It is always a good idea to have a look at your files. You can open some of them to check them by eye. To do this, open another terminal and open DS9 by typing "ds9". Back to the IRAF terminal, type for instance:
  33 
  34 display M53_-007_R.fit 1
  35 display M53_-007_G.fit 2
  36 display M53_-007_B.fit 3


   1 # First, subtract the masterbias frame from all your scientific images. We named the output files as the input files, adding a "biassub_" prefix:
   2 imarith @M53-all.list - master_bias.fits  biassub_@M53-all.list &
   3 # Now, subtract the masterdark frames from all your scientific images (in this example, all of them have an exposure time of 180s, so we can do it in a single command line). We named the output files as the input files, adding a "darksub_" prefix:
   4 imarith biassub_@M53-all.list - master_dark_180s_-25deg.fits  darksub_@M53-all.list &
   5 
   6 
   7 # divide by flat field with ccdproc
   8 # prepare the task; leave the 'flat' parameter free
   9 epar ccdproc
  10 
  11 PACKAGE = ccdred
  12    TASK = ccdproc
  13     
  14 images  =                       List of CCD images to correct
  15 (output =                     ) List of output CCD images
  16 (ccdtype=               object) CCD image type to correct
  17 (max_cac=                    0) Maximum image caching memory (in Mbytes)
  18 (noproc =                   no) List processing steps only?
  19 
  20 (fixpix =                   no) Fix bad CCD lines and columns?
  21 (oversca=                   no) Apply overscan strip correction?
  22 (trim   =                   no) Trim the image?
  23 (zerocor=                   no) Apply zero level correction?
  24 (darkcor=                   no) Apply dark count correction?
  25 (flatcor=                  yes) Apply flat field correction?
  26 (illumco=                   no) Apply illumination correction?
  27 (fringec=                   no) Apply fringe correction?
  28 (readcor=                   no) Convert zero level image to readout correction?
  29 (scancor=                   no) Convert flat field image to scan correction?
  30 
  31 (readaxi=                 line) Read out axis (column|line)
  32 (fixfile=                     ) File describing the bad lines and columns
  33 (biassec=                     ) Overscan strip image section
  34 (trimsec=                     ) Trim data section
  35 (zero   =                     ) Zero level calibration image
  36 (dark   =                     ) Dark count calibration image
  37 (flat   =                     ) Flat field images
  38 (illum  =                     ) Illumination correction images
  39 (fringe =                     ) Fringe correction images
  40 (minrepl=                   1.) Minimum flat field value
  41 (scantyp=            shortscan) Scan type (shortscan|longscan)
  42 (nscan  =                    1) Number of short scan lines
  43 
  44 (interac=                   no) Fit overscan interactively?
  45 (functio=             legendre) Fitting function
  46 (order  =                    1) Number of polynomial terms or spline pieces
  47 (sample =                    *) Sample points to fit
  48 (naverag=                    1) Number of sample points to combine
  49 (niterat=                    1) Number of rejection iterations
  50 (low_rej=                   3.) Low sigma rejection factor
  51 (high_re=                   3.) High sigma rejection factor
  52 (grow   =                   0.) Rejection growing radius
  53 (mode   =                   ql)
  54 
  55 
  56 # Flat-field division with ccdproc and correct parameters, for each filter :
  57 ccdproc darksub_@M53-R.list flat="masterflat_red.fits" &
  58 ccdproc darksub_@M53-G.list flat="masterflat_green.fits" &
  59 ccdproc darksub_@M53-B.list flat="masterflat_blue.fits" &
  60 
  61 
  62 # Check the header of one file:
  63 imheader darksub_M53_-001_R.fit l+
  64 
  65 # At the end of the header, you will notice that some keywords have been added by ccdproc (if you don't know anymore wether you reduced a frame or not, you can use them for checking):
  66 FLATCOR = 'May  8 17:08 Flat field image is masterflat_red.fits with scale=1750'
  67 CCDSEC  = '[1:4096,1:4096]'
  68 CCDPROC = 'May  8 17:09 CCD processing done'

6. Align and stack the calibrated images 

   1 #            SECTION     NPIX     MEAN   MEDIAN   STDDEV      MIN      MAX
   2 [2073:2077,2922:2926]       25    1535.    1540.    24.39    1479.    1572.
   3 #            SECTION     NPIX     MEAN   MEDIAN   STDDEV      MIN      MAX
   4 [3193:3197,1827:1831]       25    1537.    1546.    24.71    1493.    1580.
   5 #            SECTION     NPIX     MEAN   MEDIAN   STDDEV      MIN      MAX
   6 [1493:1497,1415:1419]       25    1548.    1544.    23.35    1509.    1610.

So let us use a value of 24 for this example. Leave the imexamine mode by moving the cursor to the window and stroking 'q'. The IRAF window should be back to normal command line mode.

We now set necessary parameters for the star finding algorithm (you can find a thorough description of DAOPHOT and its main parameters in the Data Analysis section): 

   1 # to load the necessary packages type in one after another:
   2 digiphot
   3 daophot
   4 
   5 # Edit the main parameters of daofind (defined in the datapars and findpars parameter sets):
   6 datapars
   7 
   8 PACKAGE = daophot
   9    TASK = datapars
  10 
  11 (scale  =                   1.) Image scale in units per pixel
  12 (fwhmpsf=                   6.) FWHM of the PSF in scale units
  13 (emissio=                  yes) Features are positive ?
  14 (sigma  =                  24.) Standard deviation of background in counts
  15 (datamin=                INDEF) Minimum good data value
  16 (datamax=                INDEF) Maximum good data value
  17 (noise  =              poisson) Noise model
  18 (ccdread=              rdnoise) CCD readout noise image header keyword
  19 (gain   =                egain) CCD gain image header keyword
  20 (readnoi=                   0.) CCD readout noise in electrons
  21 (epadu  =                   1.) Gain in electrons per count
  22 (exposur=              exptime) Exposure time image header keyword
  23 (airmass=              airmass) Airmass image header keyword
  24 (filter =               filter) Filter image header keyword
  25 (obstime=             date-obs) Time of observation image header keyword
  26 (itime  =                   1.) Exposure time
  27 (xairmas=                INDEF) Airmass
  28 (ifilter=                INDEF) Filter
  29 (otime  =                INDEF) Time of observation
  30 (mode   =                   ql)
  31 ($nargs =                    0)
  32 
  33 findpars
  34 
  35 PACKAGE = daophot
  36    TASK = findpars
  37 
  38 (thresho=                  100.) Threshold in sigma for feature detection
  39 (nsigma =                  1.5) Width of convolution kernel in sigma
  40 (ratio  =                   1.) Ratio of minor to major axis of Gaussian kernel
  41 (theta  =                   0.) Position angle of major axis of Gaussian kernel
  42 (sharplo=                  0.2) Lower bound on sharpness for feature detection 
  43 (sharphi=                   1.) Upper bound on sharpness for  feature detection
  44 (roundlo=                  -1.) Lower bound on roundness for feature detection
  45 (roundhi=                   1.) Upper bound on roundness for feature detection
  46 (mkdetec=                   no) Mark detections on the image display ?
  47 (mode   =                   ql)
  48 ($nargs =                    0)
  49 
  50 # Run daofind on your image (in this example, darksub_M53_-005_R.fit). Some user input is prompted, but all critical parameters have been set already. You can just press enter:
  51 daofind darksub_M53_-005_R.fit default
  52 FWHM of features in scale units (6.) (CR or value): 
  53         New FWHM of features: 6. scale units  6. pixels
  54 Standard deviation of background in counts (24.) (CR or value): 
  55         New standard deviation of background: 24. counts
  56 Detection threshold in sigma (100.) (CR or value): 
  57         New detection threshold: 100. sigma 2400. counts
  58 Minimum good data value (INDEF) (CR or value): 
  59         New minimum good data value: INDEF counts
  60 Maximum good data value (INDEF) (CR or value): 
  61         New maximum good data value: INDEF counts
  62 
  63 # A starlist with all the star detections has been created and named <original filename>.coo.1. Sort the star detections by brightness with the command:
  64 psort darksub_M53_-005_R.fit.coo.1 MAG
  65 
  66 # Now mark the stars on the image in ds9 by running the command:
  67 tvmark 1 darksub_M53_-005_R.fit.coo.1 mark=circle radii=20 color=204
  68 
  69 # GENERAL REMARK: The naming convention to avoid overwriting is that the .coo starlists are named with a number as a suffix. This number denotes the file version. If you create a new starlist later on, it will be named "darksub_M53_-005_R.fit.coo.2", and so on.
  70

We now need to select well identifiable stars for alignment. The task pstselect will choose suitable stars automatically, e.g. the brightest neighbour in any region with a size set by psfrad: 

   1 # Set the daopars parameters:
   2 daopars
   3 
   4 PACKAGE = daophot
   5    TASK = daopars
   6 
   7 (functio=                gauss) Form of analytic component of psf model
   8 (varorde=                    2) Order of empirical component of psf model
   9 (nclean =                    5) Number of cleaning iterations for computing psf model
  10 (saturat=                   no) Use wings of saturated stars in psf model computation
  11 (matchra=                   3.) Object matching radius in scale units
  12 (psfrad =                  25.) Radius of psf model in scale units
  13 (fitrad =                  10) Fitting radius in scale units
  14 (recente=                  yes) Recenter stars during fit ?
  15 (fitsky =                   no) Recompute group sky value during fit ?
  16 (groupsk=                  yes) Use group rather than individual sky values ?
  17 (sannulu=                   0.) Inner radius of sky fitting annulus in scale units
  18 (wsannul=                  11.) Width of sky fitting annulus in scale units
  19 (flaterr=                 0.75) Flat field error in percent
  20 (proferr=                   5.) Profile error in percent
  21 (maxiter=                  100) Maximum number of fitting iterations
  22 (clipexp=                    6) Bad data clipping exponent
  23 (clipran=                  2.5) Bad data clipping range in sigma
  24 (mergera=                INDEF) Critical object merging radius in scale units
  25 (critsnr=                   1.) Critical S/N ratio for group membership
  26 (maxnsta=                50000) Maximum number of stars to fit
  27 (maxgrou=                  100) Maximum number of stars to fit per group
  28 (mode   =                   ql)
  29 ($nargs =                    0)
  30 
  31 # Select alignment stars. They will be dropped in a new starlist (the 3rd argument). The syntax is
  32 # pstselect [frame name] [input starlist] [output starlist] [number of objects in output list]
  33 pstselect darksub_M53_-005_R.fit darksub_M53_-005_R.fit.coo.1 darksub_M53_-005_R_alignstar.coo 50
  34 
  35 # Check the selected stars on the display:
  36 display darksub_M53_-005_R.fit 1
  37 tvmark 1 darksub_M53_-005_R.fit.coo.1 mark=circle radii=30 color=205  
  38 pdump darksub_M53_-005_R_alignstar.coo xcenter,ycenter,id yes | tvmark 1 STDIN mark=circle radii="50" col=204 label+

If the red-circle stars shown on the image are not acceptable (too faint, too few, not stars, spurious, too close from other stars) repeat the star detection with other parameters.

NOTE: you may select stars for alignment in a manual way. Have a look at the red-circled stars on the image and take note of the id number of those you prefer. Then in a new terminal, create a new "xxx_alignstar.coo" file and manually copy-paste them from the old starlist to the new starlist. (10 stars are sufficient to obtain a good alignment.) Here is how to proceed in a convenient way: 

   1 # Outside the ds9 environment:
   2 
   3 # Make a copy the original starlist without the stars:
   4 cat darksub_M53_-005_R_alignstar.coo | grep '#' > darksub_M53_-005_R_alignstar2.coo
   5 
   6 # Copy-paste the stars you picked up one by one into the new list using their index. In this example, indexes are 336,8,74. (Put five spaces after the index to be sure to select indexes and not other numbers.)
   7 
   8 cat darksub_M53_-005_R_alignstar.coo | grep '336     ' >> sub_M13-010R60_alignstar2.coo
   9 cat darksub_M53_-005_R_alignstar.coo | grep '8     ' >> sub_M13-010R60_alignstar2.coo
  10 cat darksub_M53_-005_R_alignstar.coo | grep '74     ' >> sub_M13-010R60_alignstar2.coo
  11 
  12 # Then check if everything went fine in the new starlist.
  13

Once you are sure the selected stars ar OK, use their coordinates to align all other frames in the same filter to the reference image (here, we chose darksub_M53_-005_R.fit): 

   1 # Select x and y colums of your favorite starlist (here, darksub_M53_-005_R_alignstar.coo) and create a new, reduced one:
   2 pdump darksub_M53_-005_R_alignstar.coo xcenter,ycenter,id yes > darksub_M53_-005_R_alignstar.coo.1
   3 # The resulting file will contain 3 columns: position (x,y) and index number of every star.
   4 
   5 # Align all red frames to the reference image (darksub_M53_-005_R.fit) and name the resulting files with a "align" prefix:
   6 imalign darksub_@M53-R.list darksub_M53_-005_R.fit darksub_M53_-005_R_alignstar.coo.1 align//@M53-R.list boxsize=15 bigbox=51 trimimages=no
   7 
   8 
   9 # Exemple of output: Check the X- and Y-shift. They are small but would produce artefacts if ignored. N is the number of stars used to compute the shift. Almost all of them were used, meaning that the alignment stars were correctly selected. A very small number of used stars (e.g., 2 or 3) would have led to huge positional errors.
  10 
  11 #Shifts        Image    X-shift   Err      Y-shift   Err      N      Internal
  12 darksub_M53_-001_R.fit     -5.369 (0.003)      9.666 (0.004)   48   (0.009,0.013)
  13 darksub_M53_-002_R.fit     -4.648 (0.003)      6.147 (0.003)   48   (0.006,0.009)
  14 darksub_M53_-003_R.fit     -1.384 (0.003)      5.500 (0.004)   48   (0.005,0.009)
  15 darksub_M53_-004_R.fit     -0.120 (0.003)      2.605 (0.004)   48   (0.006,0.014)
  16 darksub_M53_-005_R.fit      0.000 (0.003)      0.000 (0.004)   48   (0.000,0.000)
  17 darksub_M53_-006_R.fit      2.095 (0.003)     -1.438 (0.004)   48   (0.008,0.016)
  18 darksub_M53_-007_R.fit      1.396 (0.003)     -3.046 (0.003)   48   (0.008,0.019)
  19 darksub_M53_-008_R.fit      3.854 (0.003)     -3.082 (0.004)   48   (0.011,0.024)
  20 darksub_M53_-009_R.fit      3.915 (0.003)     -1.297 (0.004)   48   (0.013,0.026)
  21 darksub_M53_-010_R.fit      7.158 (0.003)     -2.472 (0.004)   48   (0.009,0.018)
  22 
  23 # If you notice some frames are incorrectly aligned (they are visibly off compared to the others), you can retry the alignment for these frames with another starlist. If the image is glitched and definitely cannot be aligned, you can exclude it.
  24

The aligned, calibrated images taken with the red filters are finally combined using the task imcombine: 

   1 # first create a list with the aligned frames:
   2 files alignM53_-*R.fit > alignM53-R.list
   3 
   4 # Set parameters of imcombine as the following:
   5 epar imcombine
   6 
   7 PACKAGE = immatch
   8    TASK = imcombine
   9     
  10 input   =                       List of images to combine
  11 output  =                       List of output images
  12 (headers=                     ) List of header files (optional)
  13 (bpmasks=                     ) List of bad pixel masks (optional)
  14 (rejmask=                     ) List of rejection masks (optional)
  15 (nrejmas=                     ) List of number rejected masks (optional)
  16 (expmask=                     ) List of exposure masks (optional)
  17 (sigmas =                     ) List of sigma images (optional)
  18 (logfile=              logfile) Log file
  19 
  20 (combine=               median) Type of combine operation
  21 (reject =              sigclip) Type of rejection
  22 (project=                   no) Project highest dimension of input images?
  23 (outtype=                 real) Output image pixel datatype
  24 (outlimi=                     ) Output limits (x1 x2 y1 y2 ...)
  25 (offsets=                 none) Input image offsets
  26 (masktyp=                 none) Mask type
  27 (maskval=                    0) Mask value
  28 (blank  =                   0.) Value if there are no pixels
  29 
  30 (scale  =             exposure) Image scaling
  31 (zero   =                 none) Image zero point offset
  32 (weight =                 none) Image weights
  33 (statsec=                     ) Image section for computing statistics
  34 (expname=             exposure) Image header exposure time keyword
  35 
  36 (lthresh=                INDEF) Lower threshold
  37 (hthresh=                INDEF) Upper threshold
  38 (nlow   =                    1) minmax: Number of low pixels to reject
  39 (nhigh  =                    1) minmax: Number of high pixels to reject
  40 (nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
  41 (mclip  =                  yes) Use median in sigma clipping algorithms?
  42 (lsigma =                   3.) Lower sigma clipping factor
  43 (hsigma =                   3.) Upper sigma clipping factor
  44 (rdnoise=                   0.) ccdclip: CCD readout noise (electrons)
  45 (gain   =                   1.) ccdclip: CCD gain (electrons/DN)
  46 (snoise =                   0.) ccdclip: Sensitivity noise (fraction)
  47 (sigscal=                  0.1) Tolerance for sigma clipping scaling corrections
  48 (pclip  =                 -0.5) pclip: Percentile clipping parameter
  49 (grow   =                   0.) Radius (pixels) for neighbor rejection
  50 (mode   =                   ql)
  51 
  52 # Combine your scientific, aligned images:
  53 imcombine @alignM53-R.list output=align_M53red &
  54 
  55 # Make a backup copy of all your final images:
  56 imcopy align_M53*fits %align_%final_%M53*fits &

Finally, repeat the last step for green and blue images.

NOTES:

So you can simply do 

   1 imalign darksub_@M53-all.list darksub_M53_-005_R.fit darksub_M53_-005_R_alignstar.coo.1 align//@M53-G.list boxsize=15 bigbox=51 trimimages=no
   2

And then combine the frames per filter to produce your final reduced frames.

You can now critically examine the results. Repeat any necessary step. Compare the well calibrated, stacked images, to the raw data, and be proud!

Here are the files used in that example:

/!\ T.B.D.: Include intermediate images

[ATTACH]

[ATTACH]

[ATTACH]

Nice image of M53 with R, V, and B filter!