Astrometry

   Comments in italics are from Lindsey Davis, from the NOAO staff.

Draft Astrometry Cookbook

1. You need the J2000.0 (important! NOT B1950!) coordinates of the approximate center of your field and the width of the field. You should also have a rough idea of the magnitudes of the faintest well-detected stars.
2. Go to the US Naval Observatory PMM website, http://ftp.nofs.navy.mil/projects/pmm/ . From there, go to the "Finding Chart Program". Use the Basic or Detailed form as you see fit. (Basic is good enough unless you are determined to get the orientation of the map exactly right on the page.) Enter the coordinates of your field and the size. Set the faint magnitude limit to something corresponding to your actual data. You don't need the pixels. Other parameters are according to taste.
3. When you submit the request you get a sky plot of your field. Get the postscript file and print it. Then get the USNO A2.0 Star List file. Save it as plain text.
4. Make a copy of the star list and edit it so that it will be readable by IRAF. Delete the text at the top, leaving only the data. Delete everything to the right of the declination in each line. Be sure the equatorial coordinates are delimited by colons, i.e. hh:mm:ss.ss dd:mm:ss.ss . Don't round off anything, you need ALL the decmial points. Let's say this file is called fred.usno .

Steps 1-4 look fine. My only additional comment is that it may not be necessary to delete everything to the right of the declination coordinate as ccfind will append the x / y columns to the end of each input coordinate line and ccmap can be programmed to use any ra / dec / x / y columns. It is necessary to delete this data if any of the columns trailing the declination column may be blank as this will confuse the simple text file decoding routines.

5. Load the images package, then the imcoords package.
6. Check the header of your image and make sure that the keyword EPOCH is set to the date of the exposure (expressed in years, to the nearest hundredth; e.g. June 29, 2000 is 2000.49.) If it isn't, do (for example)

           hedit fred epoch 2000.49
   

   If the EPOCH keyword isn't in the header at all, you have to add it. If this is the case, append "add+" to the end of the above command.

If I understand your cookbook correctly I think step 6 should only be neccessary if:

1. you need to define the epoch (actually the equinox) of prexisting CRVAL[1/2] keywords which are not in the J2000 system and which are used by the ccfind task to transform the catalog coordinate system to the image coordinate system in order to locate objects if usewcs=yes.
2. you set the ccmap parameter refpoint to "user" and set refsystem, lngref, latref to user values. You are using refpoint=coords which is usually ok so this may not matter but see note under step 11.

It is probably better to use equinox (new FITS convention) rather than epoch (old FITS convention above) in step 6 since IRAF will use equinox in preference to epoch if both are present.

7. Check the header to see if the wcs coordinate system parameters are present. These are things like CTYPE1, CTYPE2, CRVAL1, CRVAL2, CRPIX1, CRPIX2, CD1_1, CD1_2, CD2_1, and CD2_2. If most of them are there then press ahead. If none of them are there, do

           wcsreset fred physical
           wcsreset fred world
   

Resetting the physical coordinate system is a good idea if you don't want to retain information about the coordinate system of any parent image. People dealing with mosaic data might not want to do with this however. Resetting the world coordinate system should not necessary as this should get overwritten by the ccmap or ccsetwcs tasks ...

8. If you don't already know this, use imexamine or some other task to measure the FWHM of the PSF (measured in pixels) from your stellar images. You don't need an exact number; 10% accuracy should be good enough unless you are concerned about discriminating between faint stars and faint galaxies. The 'a' command in imexamine will also give you total counts in your brighter stars; take note of this too.
9. Use the starfind task to find pixel coordinates of some of the brighter stars in the field. (More than 3 is essential, more than 6 is overkill at this point.) Give the task the HWHM (half of the FWHM) of your PSF. Set the detection threshold to somewhere around the total counts you measured above. (You can adjust this to get the number of stars you want.) Call the output file something appropriate like fred.stars .
10. Using the sky map and the star list from USNO, the stars file you just created, and your image, match up your few bright stars with their true coordinates. Use your favorite text editor to write a file that contains one line for each of these stars, and 4 columns containing RA, Dec, Pixel X, and Pixel Y coordinates. Be sure that the equatorial coordinates are delimited by colons, as described in step 4, and that nothing is rounded off. Let's imagine this file is called fred.coords1 .
11. Use the ccmap task to do a first approximation to the plate solution. In this example, the input coordinate file would be fred.coords1 and the input image would be fred.imh. Call the output database something like fred.db1. Plate solution names and summary file names can be left blank. Set xcolumn, ycolumn, lngcolumn, and latcolumn to 3, 4, 1, and 2. Make sure update is set to "yes". Other parameters should be OK at their default values. (Check that insystem = "j2000", refoint = "coords", projection = "tan", fitgeometry = "general", and function = "polynomial", the orders are all 2, and the cross terms are "half".) When you run the task you get put in an interactive fitting program, which works like other IRAF fitting programs. (Mechanics of fitting are not described here. Type ? for help.) All you need now is for the RMS residuals to be smaller than your HWHM.

A quick note about refpoint. If you use refpoint=coords then the tangent point is set to the mean of the input ra and dec columns. If these quantities are scattered more or less uniformly around the true tangent point (usually the center of the image) then everything is probably ok. This is not ok if the tangent point is actually off the image as it may be for mosaic data, or all the coordinates are in one part of a large image. Then you should set refpoint=user and set lngref, latref and maybe refsystem to known quantities.

12. Assuming that you got a satisfactory fit with all the orders = 2, then the current plate solution is recorded in the image header, and you can use the ccfind task to find more stars in the image and match them against the USNO star list. This is a little tricky. The list input celestial coordinate file is your edited usno list (fred.usno). The input image is of course the image (fred) and you can call the output matched list something like fred.coords2 . Make sure usewcs is yes. Most of the other defaults should work; check that insystem and refsystem are both j2000. The parameters to play with are sbox, cbox, datamin, and background. Make cbox about the diameter of the full image (as seen on the display) of one of the stars you used above, and sbox about twice that. Check the typical pixel values in the sky and set background to be in that ballpark. Finally, datamin should be around the peak counts in the faintest star that you want to try to match up. Your output should be a 4-column list in the same format as your fred.coords1 file.
13. You can see what stars you have identified by doing the following. (There SHOULD be a better way to do this but I haven't found it.) Make a copy of the output file from step 12 (call it something like fred.mark). Remove everything except the list of X and Y pixel coordinates. Display your image in, for instance, frame 1. Then run the tvmark task with frame = 1, coords = "fred.mark", mark = "plus", color = 0, and pointsize = 20. You'll get TEENY black plus signs in the centers of all the stars that you've matched. (Well, sort of. tvmark puts its marks 1 pixel to the left of where they ought to go.) If you don't like what you see, edit your fred.coords2 file or try fiddling with the parameters in ccfind. (E.g., if you've set datamin too low, ccfind will "find" stars that don't exist.) Hopefully you'll be able to get a dozen stars or more matched up.

Alternate method to steps 9-13 ? I have a couple of suggestions here which may make it unnecessary to run ccfind at all and requires only one run of ccmap. In step 9 set the starfind detection threshold so as to detect most of the USNO stars. In step 10 determine the ra / dec / x / y of 3 well separated stars, create a tie point file which looks like

ra1 dec1 ra2 dec2 ra3 dec3
x1 y1 x2 y2 x3 y3

and match up the celestial and pixel coordinates using the ccxymatch and the method illustrated in example 4. The reference coordinate file in this case would be fred.usno and the input coordinate file the output of starfind. Then go directly to step 14 and run ccmap on the output of ccxymatch. The starfind coordinate should be accurate enough to give reasonable astrometry ...

14. When you're happy with your "coords2" file, go back and run ccmap again, using this file as the new input, and changing the output database name, say, to fred.db2 . Use the interactive fitting routine to increase the order of the fit (remember there are FOUR order parameters!) and/or delete stars until the fit is good. (One ought to be able to get an RMS substantially smaller than 1/10 of the HWHM if all goes well.) The output will be your final plate solution.

My experience with the USNO catalog is that the rms is usually around 0.20-0.30" and that the main source of scatter is in the catalog rather than the x / y centers assuming reasonable sampling.

15. Now the whole point: to measure coordinates of other objects in the image, use the cctran task. The input coordinate file is a list of X and Y pixel coordinates; other parameters should be fairly obvious. Make sure geometry = "geometric" and forward = "yes". As a check, run this on your output file from starfind; the results ought to do a pretty decent job of recovering the coordinates from the USNO star list. (Most stars' positions ought to be recovered to within hundredths of an arcsec. But don't expect every star to be perfect; some stars have large proper motions and are not in the same places they were when the plates the USNO catalog was built from were taken.) You may want to use starfind to get the pixel coordinates of your new objects. The elongation limits may be useful in picking out the objects (e.g. galaxies, asteroid trails) you want.

If you want the higher order terms of any fit recorded in the image header by ccmap or ccsetwcs then set projection=tnx in ccmap and this will be done. Any IRAF task will understand tnx although non-IRAF software will not. This avoids the necessity of keeping the database file around. The tasks wcsctran / skyctran can transfrom from pixel to celesitial coordinates and vice versa using the image header ...