******************************** * Make an Astrometric Solution * ******************************** 05.16.00 ----------------------- The Verbose Version ----------------------- 1. Get a Digital Sky Survey image of the same field. e.g. hdf_xdss.fits, 10' x 10' field centered on hdf from http://cadcwww.dao.nrc.ca/cadcbin/getdss/$s 2. Find stars common to DSS field and image of interest. You'll use daofind to find the stars (noao.digiphot.daophot), but you must set the threshold # of sigma, and feed it an actual value for sigma from iterstat. First do the DSS field, then do the same on the image of interest. a) Set daofind sigma threshold. In noao.digiphot.daophot, epar daofind, then scroll down to the parameter findpar. Type ":e", this will bring up a sublist of parameters. Change thresho appropriately, then ":q" to quit and save out of the sublist, and ":q" to quit and save out of the daofind parameter list. e.g. thresho = 4 for DSS thresho = 20 for hawaii I-band b) Run iterstat to get sigma. e.g. iterstat hdf_xdss c) Run daofind (having set the threshold) with the sigma returned by iterstat. e.g. daofind hdf_xdss datapars.sigma=iterstat.sigma You will be prompted for the following: Output coordinate file(s): [take default] FWHM of features in scale units: [run imexam on some stars] Standard deviation of background in counts: [could find, I suppose] Detection threshold in sigma: [set here => skip step (a)] Minimum good data value: [leave blank, or use zero] Maximum good data value: [where CCD goes non-linear; not sure for plates The output will be the inputfile.coo.#, where # is the number of times you've run daofind. d) To show which objects it found, use tvmark. e.g. tvmark 1 hdf_xdss.coo.1 3. Now we match the objects found in the image of interest to the reference objects object found in the DSS image. First, you need to make a file listing objects you know to be the same between the images. The format is: x1 y1 x2 y2 x3 y3 (reference image) x1 y1 x2 y2 x3 y3 (input image) Call this dsshaw_match.ref, and set the following at the command line: input = hawaii.coo.4 The input lists referenc= hdf_xdss.coo.3 The reference lists output = dsshaw_match.dat The output matched coordinate lists toleranc= 10. The matching tolerance in pixels (refpoin= dsshaw_match.ref) Optional list of reference points e.g. xyxy hawaii.coo.1 hdf_xdss.coo.1 dsshaw_match.dat 10 refpo=dsshaw_match.ref 4. Now we have our paired objects, and we want to generate a solution between them. Feed the DSS star xy values to xyeq to get their RA/decs: e.g. xyeq yes hdf_xdss "" dsshaw_match.dat iraf 1 2 0 > radec.dat then match these reference RA/Decs to the input XYs (with match.f, for instance); put results in "xyradec.dat". This file looks like: 1794.955 104.281 12:36:40.09 62:08:40.31 1207.862 204.746 12:36:55.98 62:08:58.00 2873.201 151.322 12:36:11.07 62:08:51.45 2102.547 240.967 12:36:31.55 62:09:08.62 5. From the file containg the input xy and the reference RA/dec, generate a plate solution for the image of interest with images.imcoords.ccmap. PACKAGE = imcoords TASK = ccmap input = xyradec.dat The input coordinate files database= dsshaw.db The output database file (solutio= ) The database plate solution names (images = hawaii) The input images (results= dsshaw.res) 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 = 2.) 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) You could probably be smarter in the interactive mode than I have been, but this seems to do an okay job without any finessing. This outputs a database file which contains the fit parameters and a residual file which contains the original points plus their transformed brethren along with the residuals. Since I have "update" set to "yes", this automatically updates the image header so that imexam reads out RA/dec, and saotng displays it automatically when you scroll. Sweet. --------------------- The Terse Example --------------------- % Display DSS, find stars, and mark % sigma threshold = 4 display hdf_xdss.fits iterstat hdf_xdss daofind hdf_xdss datapars.sigma=iterstat.sigma tvmark 1 hdf_xdss.coo.1 % Display hawaii, find stars, and mark % sigma threshold = 25 display hawaii iterstat hawaii daofind hawaii datapars.sigma=iterstat.sigma tvmark 1 hawaii.coo.1 % Match DSS (reference) to Hawaii (input) % Make file of reference points: % 31.13 495.32 159.80 270.39 462.75 523.65 % 68.05 2572.48 713.87 1350.54 2382.37 2650.20 xyxy hawaii.coo.1 hdf_xdss.coo.1 dsshaw_match.dat 10 refpo=dsshaw_match.ref % Get equatorial coordinates for matched DSS objects... xyeq yes hdf_xdss "" dsshaw_match.dat iraf 1 2 0 > radec.dat % ...then match those equatorial coordinates to the Hawaii pixel % positions (use match.f to make the file xyradec.dat), and generate % the plate solution. You'll need to epar images.imcoords.ccmap; see above. ----------- match.f ----------- program match implicit none c To be used when generating astrometic solutions for an c input image based on a reference image. Program reads c the output of xyeq, which contains the reference image c xy coordinates matched to equatorial coordinates, and c matches those equatorial coordinates with the INPUT image c xy coordinates so that a plate solution for the input c image can be generated. integer p, i, j, k parameter (p=72) real x(p), y(p), junk(p) integer ra1(p), ra2(p), dec1(p), dec2(p) real ra3(p), dec3(p) i=1 j=1 k=1 c open(unit=10, file='temp1.dat', status='old') c open(unit=11, file='temp2.dat', status='old') open(unit=10, file='dsshaw_match.dat',status='old') open(unit=11, file='radec.dat', status='old') open(unit=12, file='xyradec.dat', status='unknown') c Want the input image xy from match file. Want the reference c image RA/Dec from xyeq output. read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) read(10,*) do 2000 i=1,p read(10,*,end=1000) junk(i), junk(i), x(i), y(i), junk(i), junk(i) write(*,*) junk(i), junk(i), x(i), y(i), junk(i), junk(i) 1000 continue 2000 continue read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) read(11,*) do 4000 j=1,p read(11,2500,end=3000) junk(j), junk(j), ra1(j), ra2(j), ra3(j), & dec1(j), dec2(j), dec3(j) write(*,*) ra1(j), ra2(j), ra3(j), dec1(j), dec2(j), dec3(j) 2500 format(3x, f5.2, 4x, f5.2, 2x, I2, x, I2, x, f5.2, 2x, I2, x, I2, x, f5.2) 3000 continue 4000 continue c Write data to file in form: c input X input Y ref RA ref Dec do 5000 k=1,p write(12,4500) x(k), y(k), ra1(k), ra2(k), ra3(k), dec1(k), dec2(k), dec3(k) 4500 format(2x, f8.3, 2x, f8.3, 2x, I2,':',I2,':',f5.2, 2x, I2,':',I2,':',f5.2) 5000 continue close(unit=10) close(unit=11) close(unit=12) end