************************ * Reducing Gemini Data * ************************ 16 Jan 2003 I haven't done a full, careful reduction like this one in a while. I think the terse treatment in the subsequent link is just as good and like 1.5 million times faster.. Maybe you should start there, poncho. 05.09.01 In particular, I'm going to describe the data reduction that went into my project to get IR imaging of faint radio sources. For the record, these were all steep spectrum sources in the NDWFS. We imaged each target field in J and K' (simultaneously) for 2 x 18 minutes, where each 18 minute exposure was a 2 minute, 9 point dither. None of the targets had optical counterparts in the DSS. Since the whole idea is to see if we can find counterparts in J and K', to these radio sources, I intend to cover quick astrometry and photometry at the end of these notes. INTRODUCTORY TIPS: ----------------- a) The Gemini naming convention is: 08maa036.fts - frame 36 from channel a, taken 08may2001 08mab036.fts - frame 36 from channel b, taken 08may2001 For us, channel a is J; channel b is K'. Note that you will have to change the .fts to .fits in order to use IRAF. b) To look at a directory made on shane from orthrus, you need to precede the shane path with "/net/shane". That is: /data/08ma on shane is /net/shane/data/08ma on orthrus. c) The claim in the manual is that Gemini reaches J = 19.0 K' = 17.2 in a one minute integration. A SAMPLE REDUCTION ------------------ 1. Change .fts extension to .fits extenstion. There are definitely better ways to do this (e.g. with awk & sed). My way is inelegant, but it works just fine. I assume your directory contains nothing but the night's .fts files. a) Make list of file names: > ls > filenames b) Make first-pass .csh script "runme" > awk '{ printf("mv %s %s.fits\n",$1,$1)}' filenames > runme c) "Runme" won't look quite right. Open with VI and type: :.,$s/\.fts\.fits/\.fits/g Also, remove the line "mv filenames filenames.fits" from the end of the file. d) Now execture .csh script "runme" > source runme 2. Make dark frames. You'll need a dark frame for each combination of exposure times and co-adds for each filter. We had three combinations for each of J and K' (for objects, twiflats, and standards). For the record, the combinations were: J K' ------ ------ 10 x 1 5 x 2 (for twiflats) -> dark_J_10x1 / dark_K_5x2 2 x 5 2 x 5 (for standards) -> dark_J_2x5 / dark_K_2x5 40 x 3 6 x 20 (for objects) -> dark_J_40x3 / dark_K_6x20 Make the dark frames with imcombine set to "average" with min/max rejection and no scaling. The naming convention I inherited from Mark is shown above. To speed things up, I made lists like the following: :::::::::::::: :::::::::::::: list.dark_J_40x3 list.dark_K_2x5 :::::::::::::: :::::::::::::: 08maa013.fits 08mab197.fits 08maa014.fits 08mab198.fits 08maa015.fits 08mab199.fits 08maa016.fits 08mab200.fits 08maa017.fits 08mab201.fits 08maa018.fits 08mab202.fits 08maa019.fits 08mab203.fits 08maa020.fits 08mab204.fits 08maa021.fits 08mab205.fits 08maa022.fits 08mab206.fits 08maa023.fits and then made all the darks with this IRAF script: imcombine @list.dark_J_10x1 dark_J_10x1 combine=average reject=minmax scale=none imcombine @list.dark_J_2x5 dark_J_2x5 combine=average reject=minmax scale=none imcombine @list.dark_J_40x3 dark_J_40x3 combine=average reject=minmax scale=none imcombine @list.dark_K_2x5 dark_K_2x5 combine=average reject=minmax scale=none imcombine @list.dark_K_5x2 dark_K_5x2 combine=average reject=minmax scale=none imcombine @list.dark_K_6x20 dark_K_6x20 combine=average reject=minmax scale=none When you're done, look at the darks and check the counts. The object darks should have about 10 counts (J) and 10^2 counts (K'). Note that frames with longer exposure times actually have fewer counts in the darks. IR detectors work by having the charge in the capacitor decrease as photons come along and knock them out. 3. Make twilight flats. This is a bit complicated owing to the nature of IR observing. You can treat the J twiflats essentially like you would in the optical -- subtract off the appropriate dark frame and then imcombine. The K' side must be treated more carefully, since it's sensitive to both the sky and to the thermal signature of the dome... I treat each separately. Incidentally, we're employing the Bunker method, which is to take twiflats of fixed exposure time knowing that your first flats will be saturated. This gives you a big range of count levels, which will be useful when we get to the K' band twiflats, below. J-band twilight flats --------------------- a) Weed out saturated frames. A complication (which we actually don't encounter here, because we had only 1 coadd) is that the output frames have counts which are equal the number of counts in the individual integration times the number of coadds. On the other hand, the saturation values are quoted in counts per coadd: channel A: 0 < DN < 5000 (counts / coadd) channel B: 0 < DN < 4000 (counts / coadd) Hence, if you had two coadds, your channel A saturation would be 10,000 counts. To actually do this, I made a list of all of the twiflats and ran imstat on each. Be careful! When the data numbers get high, they wrap! Don't be fooled by the fact that your first, definitely saturated frames may have low numbers -- watch them unwrap, perhaps several times. Since our integration was 10x1, I kept only frames with means < 5000. Check it: imstat 08maa034.fits 08maa034.fits 65536 246.5 350.6 imstat 08maa035.fits --> 08maa035.fits 65536 4022. 4187. imstat 08maa036.fits 08maa036.fits 65536 4645. 4781. These are the first frames -- don't be fooled! They are saturated! b) Subtract the appropriate dark frame from the raw (not saturated) twiflat frame: imarith 08maa053.fits - dark_J_10x1 08maa053_ds imarith 08maa054.fits - dark_J_10x1 08maa054_ds imarith 08maa055.fits - dark_J_10x1 08maa055_ds imarith 08maa056.fits - dark_J_10x1 08maa056_ds imarith 08maa057.fits - dark_J_10x1 08maa057_ds imarith 08maa058.fits - dark_J_10x1 08maa058_ds c) Now combine dark-subtracted twiflats with imcombine; "scale" and "weight" by "mode". (This is to compensate for the fact that the twiflats have enormously different numbers of counts, given the fixed exposure time throughout.) Accordingly, the result should have on the order of 10^3 counts. PACKAGE = immatch TASK = imcombine input = @list.imcombineJ List of images to combine output = twi_J_10x1 List of output images (rejmask= ) List of rejection masks (optional) (plfile = ) List of pixel list files (optional) (sigma = ) List of sigma images (optional) (logfile= STDOUT) Log file (combine= average) Type of combine operation (reject = minmax) Type of rejection (project= no) Project highest dimension of input images? (outtype= real) Output image pixel datatype (offsets= none) Input image offsets (masktyp= none) Mask type (maskval= 0.) Mask value (blank = 0.) Value if there are no pixels (scale = mode) Image scaling (zero = none) Image zero point offset (weight = mode) Image weights (statsec= ) Image section for computing statistics (expname= ) Image header exposure time keyword (lthresh= INDEF) 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= 6) ccdclip: CCD readout noise (electrons) (gain = 3.8) 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) d) Finally, normalize the twiflat by dividing by the median (which you determine with imstat). Check it: "imstat twi_J_10x1" returns, # IMAGE NPIX MIDPT MEAN STDDEV MIN MAX MIDPT twi_J_10x1 65536 3404. 3286. 585.6 -115.1 5804. 3404. So to normalize: imarith twi_J_10x1 / 3404. twi_J_10x1_norm K'-band twilight flats ---------------------- The idea here is that the K' band twilight flats have a contribution from the sky and from the thermal properties of the dome. The trick to handling this is to take flats starting just after the point of saturation and then keep going until the counts level off. Those final frames should have just the thermal component, so you subtract them from the initial frames (which have both). Note that this means it is not necessary to subtract a dark frame, as the dark current should be present in both! a) Again, weed out saturated twiflats. Keep in mind that there were 2 coadds here, so the effective saturation value is 8000. Here's the result of running imstat on the lot of 'em: # IMAGE NPIX MIDPT MEAN STDDEV MIN MAX MIDPT 08mab034.fits 65536 2982. 2993. 601.1 -843. 6744. 2982. 08mab035.fits 65536 12074. 11945. 1007. -1394. 14967. 12074. 08mab036.fits 65536 10904. 10820. 905.6 -1122. 13811. 10904. 08mab037.fits 65536 9941. 9840. 816. -928. 12619. 9941. 08mab038.fits 65536 9063. 8978. 738.2 -786. 11957. 9063. 08mab039.fits 65536 8237. 8174. 667.5 -778. 11028. 8237. 08mab040.fits 65536 7531. 7472. 606.9 -781. 10101. 7531. 08mab041.fits 65536 6926. 6867. 555.2 -673. 9330. 6926. 08mab042.fits 65536 6421. 6363. 512.7 -415. 8634. 6421. 08mab043.fits 65536 5980. 5937. 477.9 -467. 8029. 5980. 08mab044.fits 65536 5644. 5604. 450.2 -379. 7781. 5644. 08mab045.fits 65536 5348. 5314. 426.9 -366. 7630. 5348. 08mab046.fits 65536 5144. 5105. 410. -276. 7496. 5144. 08mab047.fits 65536 4953. 4917. 394.6 -232. 7386. 4953. 08mab048.fits 65536 4779. 4750. 381.8 -216. 7353. 4779. 08mab049.fits 65536 4645. 4613. 370.7 -203. 7323. 4645. 08mab050.fits 65536 4550. 4523. 364.1 -261. 7175. 4550. 08mab051.fits 65536 4480. 4451. 358.2 -181. 7146. 4480. 08mab052.fits 65536 4438. 4410. 355.1 -180. 7173. 4438. 08mab053.fits 65536 4400. 4374. 352.6 -178. 7162. 4400. 08mab054.fits 65536 4384. 4356. 351. -277. 7165. 4384. 08mab055.fits 65536 4369. 4347. 350.1 -205. 7157. 4369. 08mab056.fits 65536 4377. 4353. 350.7 -173. 7219. 4377. 08mab057.fits 65536 4381. 4354. 350.6 -240. 7352. 4381. 08mab058.fits 65536 4375. 4350. 350.5 -179. 7205. 4375. So these guys get unsaturated starting at frame 08mab040. I actually plotted the mean counts per coadd in IDL and saw that leveling-off trend seems to start at 08mab049. So, we'll ditch 08mab034 - 08mab039; we'll consider 08mab040 - 08mab048 our true twiflats, and we'll use 08mab049 - 08mab058.fits to establish our thermal component. b) Make this "thermal twi" now. Don't dark subract it! I ran with (note that we're scaling and weighting with "mode" again): PACKAGE = immatch TASK = imcombine input = @list.therm List of images to combine output = twi_K_5x2therm List of output images (rejmask= ) List of rejection masks (optional) (plfile = ) List of pixel list files (optional) (sigma = ) List of sigma images (optional) (logfile= STDOUT) Log file (combine= average) Type of combine operation (reject = minmax) Type of rejection (project= no) Project highest dimension of input images? (outtype= real) Output image pixel datatype (offsets= none) Input image offsets (masktyp= none) Mask type (maskval= 0.) Mask value (blank = 0.) Value if there are no pixels (scale = mode) Image scaling (zero = none) Image zero point offset (weight = mode) Image weights (statsec= ) Image section for computing statistics (expname= ) Image header exposure time keyword (lthresh= INDEF) 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= 6) ccdclip: CCD readout noise (electrons) (gain = 3.8) 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) :::::::::::::: list.therm :::::::::::::: 08mab049.fits 08mab050.fits 08mab051.fits 08mab052.fits 08mab053.fits 08mab054.fits 08mab055.fits 08mab056.fits 08mab057.fits 08mab058.fits c) Subtract the "thermal twi" from each "true twi". Here's the script: imarith 08mab040 - twi_K_5x2therm 08mab040_ts imarith 08mab041 - twi_K_5x2therm 08mab041_ts imarith 08mab042 - twi_K_5x2therm 08mab042_ts imarith 08mab043 - twi_K_5x2therm 08mab043_ts imarith 08mab044 - twi_K_5x2therm 08mab044_ts imarith 08mab045 - twi_K_5x2therm 08mab045_ts imarith 08mab046 - twi_K_5x2therm 08mab046_ts imarith 08mab047 - twi_K_5x2therm 08mab047_ts imarith 08mab048 - twi_K_5x2therm 08mab048_ts d) Now combine your thermal-component-subtracted twilight frames. Use the same parameters for the imcombine on the J twiflat (i.e. weight/scale by mode). imcombine @list.thermsubK -> twi_K_5x2 e) Normalize. Imstat gives: # IMAGE NPIX MIDPT MEAN STDDEV MIN MAX MIDPT twi_K_5x2 65536 1434. 1416. 123. -246.7 1962. 1434. Hence: imarith twi_K_5x2 / 1434. twi_K_5x2_norm 4. Dark subtract all your object and standard frames. Be sure and use the appropriate integration-coadd combination for each, and separate according to band. imarith @list.obj_J - dark_J_40x3 @list.obj_J//_ds imarith @list.obj_K - dark_K_6x20 @list.obj_K//_ds imarith @list.std_J - dark_J_2x5 @list.std_J//_ds imarith @list.std_K - dark_K_2x5 @list.std_K//_ds 5. Divide all object and standard frames by appropriate twiflats. imarith @list.obj_J//_ds / twi_J_10x1_norm @list.obj_J//_dsff imarith @list.obj_K//_ds / twi_K_5x2_norm @list.obj_K//_dsff imarith @list.std_J//_ds / twi_J_10x1_norm @list.std_J//_dsff imarith @list.std_K//_ds / twi_K_5x2_norm @list.std_K//_dsff 6. Construct bad pixel masks. You can make bad pixel masks out of your flatfields with imreplace. Since they've been normalized to unity, you can figure than any pixel in the flat with values < 0.5 or > 1.5 are bad, and you can flag 'em with imreplace. Here's the mapping you'll use along with sample imreplace commands. Note that you must copy your normalized flatfield to a temporary file before acting on it with imreplace! Also, when you're done, change the .fits bad pixel masks to .pl. The mapping: pix < 0.5 -> 0 e.g. imreplace file 0 lower=INDEF upper=0.5 0.5 < pix < 1.5 -> 1 e.g. imreplace file 1 lower=0.5 upper=1.5 pix > 1.5 -> 0 e.g. imreplace file 0 lower=1.5 upper=INDEF Here's the script I used to achieve all this: imcopy twi_J_10x1_norm bpm_J imcopy twi_K_5x2_norm bpm_K imreplace bpm_J 0 lower=INDEF upper=0.5 imreplace bpm_J 1 lower=0.5 upper=1.5 imreplace bpm_J 0 lower=1.5 upper=INDEF imreplace bpm_K 0 lower=INDEF upper=0.5 imreplace bpm_K 1 lower=0.5 upper=1.5 imreplace bpm_K 0 lower=1.5 upper=INDEF imcopy bpm_J.fits bpm_J.pl imcopy bpm_K.fits bpm_K.pl 7. Start DIMSUM. a) Mark recommends making a separate directory for object. Put in that directory the flatfielded object frames and the bad pixel masks. DIMSUM's smart enough to handle the two separate 9-point integrations, so throw it all in there! Keep the bands separate, though. b) Now use the task "reduce". Some of these parameters are a bit obscure. Where confused, accept the defaults. Here are the important ones: img_lst = duh out_img = duh ref_img = use the first image in your mosaic nmed = takes 1/2 this many frames on either side of the frame being flattened bpimage = make sure to use .pl file mfac = use 2 to allow for 1/2 pixel shifts PACKAGE = dimsum TASK = reduce img_lst = @list.obj List of images out_img = J142543+33_J Base name for output firstpass and maskpass images (ref_img= 08maa064_dsff) Name of a reference image (WITHOUT prefix) in image list (prefix = ss_) Prefix for sky-subtracted images (chk = yes) Do you want chances to check shifts and masks? kpchking= no Keep checking mask? (fp_xslm= yes) Do firstpass xslm? (fp_fixp= yes) Do firstpass fixpix? (fp_xzap= yes) Do firstpass xzap? (fp_badp= yes) Do firstpass bad pixel file update? (mk_shif= yes) Determine shifts? (fp_regi= yes) Do firstpass registar? (masking= yes) Make masks? (maskder= yes) Deregister masks? (mp_xslm= yes) Do maskpass xslm? (mp_fixp= yes) Do maskpass fixpix? (mp_xzap= yes) Do maskpass xzap? (mp_badp= yes) Do maskpass bad pixel file update? (mp_regi= yes) Do maskpass registar? (nmed = 6) Number of images to use in sky frame (bpimage= bpm_J.pl) Bad pixel mask image (mfac = 2) Mag factor for maskpass image (secpexp= 1.) Seconds per unit exposure time (ncoavg = 1) Number of internal coaverages per frame (rjct = 1) Number of pixels for xslm minmax reject (nskymin= 3) Minimum number of frames to use for sky (forcesc= yes) Force recalculation of image medians in sky scaling? (premask= yes) Mask working sky images in maskpass? (pixin = yes) Are input coords measured pixels of ref object? (scale_x= 1.) x pixels per ab coord unit (scale_y= 1.) y pixels per ab coord unit (ang_a2x= 0.) Angle in degrees from a CCW to x (senseab= yes) Is a through b counterclockwise? (ang_y2N= 0.) Angle in degrees from y to N dir N through E (rotatio= yes) Is N through E counterclockwise? (nsig_ob= 1.1) factor x suggested threshold for object masking (nsig_cr= 1.5) factor x suggested threshold for cr masking (ilist = ) (mode = ql) ------------------------------------------------------ The Interactive Bits ------------------------------------------------------ (i) Give it coarse shift: DIMSUM wants you to click on the the same bright star in each image. So: a) Center the cursor on some star that you'll recognize from image to image, b) hit "a" on it, c) hit "n" to advance to next image. Watch the file names so that you only scroll through your set of images once. When you're done scrolling through, hit "q" to quit. This will bring up a VI window. Check that each image is represented. (ii) Give it a bunch of stars to look for: DIMSUM wants you to click on a bunch of stars in just one image. It'll use these stars to make a nice map between images. So: a) Center the cursor on a bunch of bright stars, b) hit "a" on them, c) hit "q" to quit. Aim for > 5 objects. After quitting, you'll get a VI window. Delete an imexam lines that didn't converge. (iii) DIMSUM now brings up another VI window. This has all the image frames with the calculated shifts. Make sure that every frame is represented. (iv) Now DIMSUM goes through two rounds of mask-making. In each, it's going to suggest a threshold value. Make sure that the sky is below the threshold and the objects are above it. (DIMSUM generally chooses well.) After you're satisfied from imexam'ing around, hit "q", and then enter the suggested value at the prompt. Subsequently, DIMSUM will ask if you want to keep checking the masks. Say no. ------------------------------------------------------ Trouble Shooting ------------------------------------------------------ (i) Crashed on first-pass crzap ("can't divide by zero"): Epar reduce and turn off first-pass crzap. Start again. (This doesn't seem to affect me when running at Berkeley -- only when on the mountain.) (ii) Crashed on maskpass sxlm ("won't overwrite ss_*"): Those frames were from first-pass sky-subtraction and should have been deleted. Delete them; epar DIMSUM and switch off steps it has already done; start again. You may also have to delete some _crap. (iii) See a trailed image? Or an image you want to exclude? Give it a zero weight when DIMSUM opens the list of shifts. (Last column of shift list.) (iv) When I tried this at home, dimsum complained about not understanding the parameter "rej" in imcombine. Mark said this had something to do with setting the version parameter in dimsum to the current version of IRAF (V2.11.3). That didn't work. Nor did anything. In the end, I just got all pissed and copied down the dimsum files from Mark's observer directory on orthrus. That got me as far as the next bug... Stupid. (v) After making the firstpass registar (_fp), I got this weird "value expression" prompt which was clearly looking at the current value of "value" in hedit. I set the hidden parameters to: (add = no) add rather than edit fields (delete = no) delete rather than edit fields (verify = yes) verify each edit operation (show = yes) print record of each edit operation (update = yes) enable updating of the image header (mode = al) and the problem went away. Stupid. (vi) Now it's crashing at the "Make masks?" stage with: ERROR on line 34: Pixel storage file is truncated (J142907+33_K_fp) sigmanorm (image=J142907+33_K_fp, prefix=norm) mkmask (image=J142907+33_K_fp, nthresh=1.5, nsmooth=3, subsample=2, filtsize=15, ngrow=0, interact=yes) reduce (mode=h) This was really weird. I won't go into the details, but the short story is I fucked things up when dark-subtracting and twi-flattening. Dimsum didn't crash on data from other nights; it didn't crash on the raw data from this night. Weird. (vii) Sometimes you get caught in this infinite loop at the mask checking stage. For some reason, I think this happes only if DIMSUM crashed at some prior stage, and you restarted the reduction. Whatever. The solution is: a) CONTROL-C your ass out of the loop, b) \rm _* in the working directory, c) epar reduce, d) and start back up with everything turned off PRIOR to "Make Masks". ------------------------------------------------------ c) So you made it through! The final image is the sum of the input object frames, which have been mapped up to arrays of 2x their original size. (This means that you have correlated noise of correlation length 4 pix in your final image. Hence, its best to get your RMS noise on the from the "_fp" image, which is unexpanded.) Final image: ROOT_mp = "mask pass" 8. Do astrometry. We want to see if we can see radio sources with accurate FIRST positions in these IR images. So, we need to put some astrometry in the header of the IR images. Then we can tool around with imexam and see if there's something at the radio position. STRATEGY: We're most interested in the K band, so you should do your astrometry there first. Having done this a bit, it looks like you need to do J and K separately... So, do the Ks and look for detections. You might be able to spot it in J without doing actual astrometry. At this point, it would be good to look at a FIRST image of the thing. a) Get a DSS finding chart for your field. A good size is 4 minutes square. We're going to use the astrometry in the header of this thing. (Use the "dss" link off my homepage, if you like.) b) Now determine an astrometric solution in your IR image. The first thing to do is get RA and Dec for objects in reference image. This is two steps: (i) Get X and Y for objects in your reference image with imexam. Pipe it to a file to make life easier. Call if ref.xy. (ii) Now covert those X's and Y's to RAs and Decs with xyeq in stsdas.analysis.gasp. Do it at the command line with: xyeq yes dss_plate_finder "" ref.xy iraf 1 2 0 > ref.radec c) Now get X and Y for those same objects in your IR image. It'll help to mark the objects you used in the reference image with: tvmark 2 ref.xy mark=circle radii=5 color=206 label=no number=yes Now just blink 'em, and imexam the same objects. Pipe it to image.xy d) Now make a file with the X and Y from image and the RA and Dec from the reference. If you've used "ref.xy", "ref.radec", and "image.xy" for your file names, just run "matchpos.pro": pro matchpos,xyfile=xyfile,posfile=posfile,outfile=outfile if NOT keyword_set(xyfile) then xyfile = 'image.xy' if NOT keyword_set(posfile) then posfile = 'ref.radec' if NOT keyword_set(outfile) then outfile = 'image.xyradec' fmt = 'A,A,A,A' readcol,posfile,refx,refy,refra,refdec,f=fmt fmt = 'A,A' readcol,xyfile,imagex,imagey,f=fmt openw,unit,outfile,/get_lun for i = 0,n_elements(imagex)-1 do begin printf,unit,imagex(i),imagey(i),refra(i+2),refdec(i+2), $ format='(A10,A10,A20,A20)' endfor free_lun,unit end e) Finally, run ccmap in the imcoords package. Recall that this bring up an interactive window: g - show you the bird's eye view of the stars being fitted f - make a new fit x - plots some residuals r - plots some residuals s - plots some residuals y - plots some residuals q - quit, and get prompted for update Basically, I cycled through x, r, s, and y and deleted outliers until my residuals sank below some fun level. The parameters I used for ccmap were: PACKAGE = imcoords TASK = ccmap input = image.xyradec The input coordinate files database= image.db The output database file (solutio= ) The database plate solution names (images = J142543+33_J_mp.fits) The input images (results= image.summary) The optional results summary files (xcolumn= 1) Column containing the x coordinate (ycolumn= 2) Column containing the y coordinate (lngcolu= 3) Column containing the ra / longitude (latcolu= 4) Column containing the dec / latitude (xmin = INDEF) Minimum logical x pixel value (xmax = INDEF) Maximum logical x pixel value (ymin = INDEF) Minimum logical y pixel value (ymax = INDEF) Maximum logical y pixel value (lngunit= hours) Input ra / longitude units (latunit= degrees) Input dec / latitude units (insyste= j2000) Input celestial coordinate system (refpoin= coords) Source of the reference point definition (lngref = INDEF) Reference point ra / longitude telescope coordinate (latref = INDEF) Reference point dec / latitude telescope coordinate (refsyst= INDEF) Reference point telescope coordinate system (lngrefu= ) Reference point ra / longitude units (latrefu= ) Reference point dec / latitude units (project= tan) Sky projection geometry (fitgeom= general) Fitting geometry (functio= polynomial) Surface type (xxorder= 2) Order of xi fit in x (xyorder= 2) Order of xi fit in y (xxterms= half) Xi fit cross terms type (yxorder= 2) Order of eta fit in x (yyorder= 2) Order of eta fit in y (yxterms= half) Eta fit cross terms type (reject = INDEF) Rejection limit in sigma units (update = yes) Update the image world coordinate system ? (pixsyst= logical) Input pixel coordinate system (verbose= yes) Print messages about progress of task ? (interac= yes) Fit the transformation interactively ? (graphic= stdgraph) Default graphics device (cursor = ) Graphics cursor (mode = ql) f) To find your radio source, fire up imexam and start hittin' "x" until you arrive at the radio source position. See anything? -------------------------------------------------------------------------------------------------------------- Sample Run (J142748+3, K band) -------------------------------------------------------------------------------------------------------------- display J142748+34_K_mp 1 # load image display dss_plate_finder 2 # load reference imexam > ref.xy # get stars X,Y off of reference vi ref.xy # edit ref.xy (want x and y column only) xyeq yes dss_plate_finder "" ref.xy iraf 1 2 0 > ref.radec # convert ref. pix to celestial tvmark 2 ref.xy mark=circle radii=5 color=206 label=no number=yes # indicate ref. stars imexam > image.xy vi image.xy # edit image.xy (want x and y column only) ## IN IDL ## matchpos.pro # make imagex imagey refra refdec ccmap image.xyradec image.db images=J142748+34_K_mp results=image.summary # do astrom.; update header -------------------------------------------------------------------------------------------------------------- 9. Do photometry. I lied -- no photometry here. See the later (terse) iteration of these notes.