**
Make a log which lists, for each exposure: the filter, exposure time, and seeing. Also note when the telescope was dithered.
**
eliminate images with awful seeing.
**
copy over a valid OPTIC wcs with wcscopy. Laura has one, if you don't. This is important! If you try to do the next step without doing this, very strange and inexplicable things will happen.
PACKAGE = immatch
TASK = wcscopy
images = *.fits[4] List of input images
refimage= msc066[4] List of reference images
(verbose= yes) Print messages about actions taken ?
(mode = ql)
(must do this four times, once for each extension)
**
edit login.cl file to include this line:
set disable_wcs_maps=""
**
Run 'mscsplit' on all images. This will split your images into five single-extension fits files. The first (_0) is just a master header. Mscsplit is smart enough to add the extension numbers to the outputted files.
ms> lpar mscsplit
input = "msc3086" List of input MEF files
(output = "ext3086") List of output root names
(mefext = ".fits") MEF filename extension
(delete = no) Delete MEF file after splitting?
(verbose = yes) Verbose?
(fd1 = "")
(fd2 = "")
(mode = "q")
**
Run KOORDS from the karma package on one image for each dither position. By one image, I actually mean you will have to run koords four times for each dither, because you'll have to run it once for each extension. I use the DSS image of my galaxy for the 'reference image' in koords. Running koords on a single extension can be a bit tricky-- sometimes it's hard to even find three reasonably bright stars! It helps if you display the entire mosaic in some ds9 window to keep you oriented.
**
Copy over the koords wcs to the rest of the images taken at that dither position using 'wcscopy'. Again, you'll have to do this once per extension.
**
Re-join the single extensions back into an MEF file, using mscjoin.
ms> lpar mscjoin
input = "@extlist" List of input root names
(output = "@msclist") List of output MEF names
(delete = no) Delete input images after joining?
(verbose = yes) Verbose?
(fd1 = "")
(fd2 = "")
(mode = "q")
**
Get a list of RAs/Decs of stars in your region of the sky using 'mscgetcat'
ms> lpar mscgetcat
input = "ext.086_1" List of Mosaic files
output = "usno2.cat" Output file of sources
(magmin = 0.) Minimum magnitude
(magmax = 26.) Maximum magnitude
(catalog = "CADC:USNO-A2") Catalog
(rmin = 16.) Minimum radius (arcmin)
(mode = "q")
Sometimes the server goes down, so if you are getting a blank file, try both the CADC and NOAO server.
This gives me a list of all USNO stars within a 16 arcminute radius of the input image.
**
Now, improve the wcs using msccmatch. I like to do this once for the entire MEF frame, to check that everything is all right, and then once for each extension independently.
When you run msccmatch, it should center your field at x=0, y=0. If it's centering it at some crazy big numbers, there's probably something wrong.
ms> lpar msccmatch
input = "@fitslist" List of input mosaic exposures
coords = "usno2.cat" Coordinate file (ra/dec)
accept = yes Accept solution?
(outcoords = "") List of updated coordinate files
(usebpm = no) Use bad pixel masks?
(verbose = yes) Verbose?\n\n# Coarse Search
(nsearch = 50) Maximum number of positions to use in search
(search = 0.) Translation search radius (arcsec)
(rsearch = 0.) Rotation search radius (deg)\n\n# Fine Centroid
(cbox = 11) Centering box (pixels)
(maxshift = 5.) Maximum centering shift to accept (arcsec)
(csig = 0.1) Maximum centering uncertainty to accept (arcsec
(cfrac = 0.5) Minimum fraction of accepted centers
(listcoords = yes) List centered coordinates in verbose mode?\n\n#
(nfit = 4) Min for fit (>0) or max not found (<=0)
(rms = 2.) Maximum fit RMS to accept (arcsec)
(fitgeometry = "general") Fitting geometry
(reject = 3.) Fitting rejection limit (sigma)
(update = yes) Update coordinate systems?
(interactive = yes) Interactive?
(fit = yes) Interactive fitting?
(graphics = "stdgraph") Graphics device
(cursor = "") Graphics cursor\n
(mode = "ql")
This will cycle through all images in the list. I like to run this program interactively-- type 'x' to see x residuals, 'y' to see y residuals, 'r' and 's' to see cross-term residuals. 'd' to delete a point, 'f' to refit, 'u' to undelete a point. 'g' will take you back to the original screen that shows the distribution of stars.
Then I rerun it four more times, once for each extension. Each time I give it a list of images, with the extension number in square brackets. This seems to help a bit with some of the crazy higher order wcs stuff, and leads to better stacks.
**
Now, you're going to have to fix the cosmic rays before you turn your MEF file into a single image. Here's what I do:
--load package crutil
**
Crmedian seems to work the best for identifying the cosmic rays. The below parameters have so far worked well for my OPTIC images, and don't seem to think that stars are cosmic rays, or anything nasty like that. But definitely make sure the below parameters are reasonable for your data (you can overlay a cr mask using 'display')
PACKAGE = crutil
TASK = crmedian
input = msc063[1] Input image
output = Output image
(crmask = ) Output cosmic ray mask
(median = ) Output median image
(sigma = ) Output sigma image
(residua= ) Output residual image
(var0 = 0.) Variance coefficient for DN^0 term
(var1 = 0.) Variance coefficient for DN^1 term
(var2 = 0.) Variance coefficient for DN^2 term
(lsigma = 10.) Low clipping sigma factor
(hsigma = 4.5) High clipping sigma factor
(ncmed = 5) Column box size for median level
calculation
(nlmed = 5) Line box size for median level calculation
(ncsig = 25) Column box size for sigma calculation
(nlsig = 25) Line box size for sigma calculation
(mode = ql)
you have to run crmed on each extension, so make a file that looks like this:
crmedian msc064[1] output="" crmask="crmsk064_1"
crmedian msc064[2] output="" crmask="crmsk064_2"
crmedian msc064[3] output="" crmask="crmsk064_3"
crmedian msc064[4] output="" crmask="crmsk064_4"
#
crmedian msc065[1] output="" crmask="crmsk065_1"
crmedian msc065[2] output="" crmask="crmsk065_2"
crmedian msc065[3] output="" crmask="crmsk065_3"
crmedian msc065[4] output="" crmask="crmsk065_4"
#
crmedian msc066[1] output="" crmask="crmsk066_1"
crmedian msc066[2] output="" crmask="crmsk066_2"
crmedian msc066[3] output="" crmask="crmsk066_3"
crmedian msc066[4] output="" crmask="crmsk066_4"
(if there are no empty lines at bottom iraf might crash)
save this to something.cl
run this script by typing
cl < something.cl
at iraf prompt
this will be sorta slow.
**
--NOW-- run crgrow to 'cover' more of the cosmic rays.
run on each mask (give it a list of input mask files. make output the same list)
PACKAGE = crutil
TASK = crgrow
input = testy[1] Input cosmic ray masks
output = testy[1] Output cosmic ray masks
(radius = 1.) Grow radius
(inval = INDEF) Input mask value to grow (INDEF for any)
(outval = INDEF) Output grown mask value (INDEF for any)
(mode = ql)
**
apply the masks with fixpix
PACKAGE = proto
TASK = fixpix
images = List of images to be fixed
masks = List of bad pixel masks
(linterp= INDEF) Mask values for line interpolation
(cinterp= INDEF) Mask values for column interpolation
(verbose= no) Verbose output?
(pixels = no) List pixels?
(mode = ql)
make a goat.cl file that looks like this:
fixpix msc063.fits[1] masks="crmsk063_1.fits[1]"
fixpix msc063.fits[2] masks="crmsk063_2.fits[1]"
fixpix msc063.fits[3] masks="crmsk063_3.fits[1]"
fixpix msc063.fits[4] masks="crmsk063_4.fits[1]"
#
fixpix msc064.fits[1] masks="crmsk064_1.fits[1]"
fixpix msc064.fits[2] masks="crmsk064_2.fits[1]"
fixpix msc064.fits[3] masks="crmsk064_3.fits[1]"
fixpix msc064.fits[4] masks="crmsk064_4.fits[1]"
#
fixpix msc065.fits[1] masks="crmsk065_1.fits[1]"
fixpix msc065.fits[2] masks="crmsk065_2.fits[1]"
fixpix msc065.fits[3] masks="crmsk065_3.fits[1]"
fixpix msc065.fits[4] masks="crmsk065_4.fits[1]"
#
fixpix msc066.fits[1] masks="crmsk066_1.fits[1]"
fixpix msc066.fits[2] masks="crmsk066_2.fits[1]"
fixpix msc066.fits[3] masks="crmsk066_3.fits[1]"
fixpix msc066.fits[4] masks="crmsk066_4.fits[1]"
#
fixpix msc067.fits[1] masks="crmsk067_1.fits[1]"
fixpix msc067.fits[2] masks="crmsk067_2.fits[1]"
fixpix msc067.fits[3] masks="crmsk067_3.fits[1]"
fixpix msc067.fits[4] masks="crmsk067_4.fits[1]"
run this:
cl < goat.cl
**
--Now, turn MEF images into a single-extension FITS images using 'mscimage'.
PACKAGE = mscred
TASK = mscimage
input = @mimage_in List of input mosaic exposures
output = @mimage_out List of output images
(format = image) Output format (image|mef)
(pixmask= yes) Create pixel mask?
(verbose= )_.verbose) Verbose output?
# Output WCS parameters
(wcssour= image) Output WCS source (image|parameters)
(referen= ) Reference image
(ra = INDEF) RA of tangent point (hours)
(dec = INDEF) DEC of tangent point (degrees)
(scale = INDEF) Scale (arcsec/pixel)
(rotatio= INDEF) Rotation of DEC from N to E (degrees)
# Resampling parmeters
(blank = 0.) Blank value
(interpo= sinc) Interpolant for data
(minterp= linear) Interpolant for mask
(boundar= reflect) Boundary extension
(constan= 0.) Constant boundary extension value
(fluxcon= no) Preserve flux per unit area?
(ntrim = 8) Edge trim in each extension
(nxblock= 2048) X dimension of working block size in pixels
(nyblock= 2048) Y dimension of working block size in pixels
# Geometric mapping parameters
(interac= no) Fit mapping interactively?
(nx = 10) Number of x grid points
(ny = 20) Number of y grid points
(fitgeom= general) Fitting geometry
(xxorder= 4) Order of x fit in x
(xyorder= 4) Order of x fit in y
(xxterms= half) X fit cross terms type
(yxorder= 4) Order of y fit in x
(yyorder= 4) Order of y fit in y
(yxterms= half) Y fit cross terms type
(fd_in = )
(fd_ext = )
(fd_coor= )
(mode = ql)
First i ran this with just msc066, to make sing066. This will be the reference image for all future runs of mscimage! I'm not sure how important it is to have the same reference image for all images, but it seems like it has potential to be significant.
This is very slow! The slowest step of all.
********************************************************************************************
From now on, split the images up by filter. Do the following steps for each filter separately!
***********************************************************************************************
--Subtract the sky. If your galaxy is not too big, the mode of the image should work just fine.
--Use mscimatch to write scaling factors into header
ms> lpar mscimatch
input = "ssingj066,ssingj3117" List of images
coords = "usno2.cat" File of coordinates
accept = yes Accept scaling and update images?
(bpm = "") List of bad pixel masks
(measured = "") Measurment file
(scale = yes) Determine scale?
(zero = no) Determine zero offset?
(box1 = 21) Inner box size for statistics
(box2 = 51) Outer box size for statistics
(lower = -100.) Lower limit for good data
(upper = INDEF) Upper limit for good data
(niterate = 3) Number of sigma clipping iterations
(sigma = 11.) Sigma clipping factor
(interactive = yes) Interactive?
(verbose = yes) Verbose?
(mode = "ql")
Curiously, mscimatch gives much better results if you only run it on one image at a time with a given reference image. If you give it a list of images, it does not match up every frame to a reference-- instead it progresses through the list, and matches the current image to the previous image, thus propagating errors. I really reccomend that you only give mscimatch two files at a time.
**
Combine the images using mscstack.
You can apply the bad pixel masks which came from mscimage using mscstack. In your image headers, there should be a 'BPM' entry, and it should say the bad pixel mask which corresponds to your image. If it doesn't, use hedit to correct it. Mscstack is smart enough to use the bpm that is listed in your header, as long as masktype does NOT equal 'none'. In the below parameters mscstack masks everything that has a value of 1 in your bad pixel masks.
PACKAGE = mscred
TASK = mscstack
input = @w16_ssing.list List of images to combine
output = stack_w16 Output image
(headers= ) List of header files (optional)
(bpmasks= ) List of bad pixel masks (optional)
(rejmask= ) List of rejection masks (optional)
(nrejmas= ) List of number rejected masks (optional)
(expmask= ) List of exposure masks (optional)
(sigmas = ) List of sigma images (optional)
(combine= average) Type of combine operation (median|average)
(reject = avsigclip) Type of rejection
(masktyp= badvalue) Mask type
(maskval= 1.) Mask value
(blank = 0.) Value if there are no pixels
(scale = !mscscale) Image scaling
(zero = !msczero) Image zero point offset
(weight = exposure) Image weights
(statsec= ) Image section for computing statistics
(lthresh= 1.) Lower threshold
(hthresh= 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)
(sigscal= 0.1) Tolerance for sigma clipping scaling corrections
(pclip = -0.5) pclip: Percentile clipping parameter
(grow = 0.) Radius (pixels) for neighbor rejection
(mode = ql)
Hopefully this will produce a nice deep image!
