Changeset 3389 for trunk/Ohana/src/addstar

Timestamp:

Mar 7, 2005, 6:15:46 PM (22 years ago)

Author:

eugene

Message:

adding mosaic astrometry concepts

Location:

trunk/Ohana/src/addstar

Files:

• Makefile (modified) (2 diffs)
• include/addstar.h (modified) (1 diff)
• src/args.c (modified) (5 diffs)
• src/find_matches.c (modified) (5 diffs)
• src/gimages.c (modified) (6 diffs)
• src/gregion_image.c (modified) (1 diff)
• src/gstars.c (modified) (1 diff)

trunk/Ohana/src/addstar/Makefile

@@ +1 @@
+#install:
+#       @echo "addstar under re-construction"
 default: addstar
 help:
@@ -110 +112 @@
 # utilities #################################################
 
-install:
+# install:
         for i in $(INSTALL); do make $$i.install; done
 

trunk/Ohana/src/addstar/include/addstar.h

@@ -62 +62 @@
 PhotCode *thiscode;
 GSCRegion UserPatch;
+Coords *MOSAIC;
 
 char DateKeyword[64];

trunk/Ohana/src/addstar/src/args.c

@@ -14 +14 @@
 int args (int argc, char **argv) {
   
-  int i, N;
+  int N;
 
   /* check for help request */
@@ -61 +61 @@
   /* provide a time for dataset */
   TIMEREF = 0; 
-  if ((i = get_argument (argc, argv, "-time"))) {
-    remove_argument (i, &argc, argv);
-    if (!str_to_time (argv[i], &TIMEREF)) { 
+  if ((N = get_argument (argc, argv, "-time"))) {
+    remove_argument (N, &argc, argv);
+    if (!str_to_time (argv[N], &TIMEREF)) { 
       fprintf (stderr, "syntax error in time\n");
       exit (1);
     }
-    remove_argument (i, &argc, argv);
+    remove_argument (N, &argc, argv);
+  }
+  /* provide a mosaic for distortion */
+  MOSAIC = NULL;
+  if ((N = get_argument (argc, argv, "-mosaic"))) {
+    Header header;
+    ALLOCATE (MOSAIC, Coords, 1);
+
+    remove_argument (N, &argc, argv);
+    if (!fits_read_header (argv[N], &header)) {
+      fprintf (stderr, "ERROR: can't read header for mosaic %s\n", argv[N]);
+      exit (1);
+    }
+    if (!GetCoords (MOSAIC, &header)) {
+      fprintf (stderr, "ERROR: no astrometric solution in header\n");
+      exit (1);
+    }
+    if (strcmp(&MOSAIC[0].ctype[4], "-DIS")) {
+      fprintf (stderr, "ERROR: not a mosaic distortion header\n");
+      exit (1);
+    }
+    remove_argument (N, &argc, argv);
+    fits_free_header (&header);
   }
   
@@ -73 +95 @@
   /* only add to existing objects */
   ONLY_MATCH = FALSE;
-  if ((i = get_argument (argc, argv, "-only-match"))) {
+  if ((N = get_argument (argc, argv, "-only-match"))) {
     ONLY_MATCH = TRUE;
-    remove_argument (i, &argc, argv);
+    remove_argument (N, &argc, argv);
   }
   /* don't add missed pts to Missed table (image only) */
@@ -85 +107 @@
   /* replace measurement, don't duplicate (ref/cat only) */
   REPLACE = FALSE;
-  if ((i = get_argument (argc, argv, "-replace"))) {
+  if ((N = get_argument (argc, argv, "-replace"))) {
     REPLACE = TRUE;
-    remove_argument (i, &argc, argv);
+    remove_argument (N, &argc, argv);
   }
   /* only add image potion to image table */
@@ -133 +155 @@
     fprintf (stderr, "USAGE: addstar -ref (filename)\n");
     fprintf (stderr, "  cur: %s ", argv[0]);
-    for (i = 0; i < argc; i++) { fprintf (stderr, "%s ", argv[i]); }
+    for (N = 0; N < argc; N++) { fprintf (stderr, "%s ", argv[N]); }
     fprintf (stderr, "\n");  
     exit (2);

trunk/Ohana/src/addstar/src/find_matches.c

@@ -49 +49 @@
   REALLOCATE (catalog[0].missing, Missing, NMISS);
   
-  /* project onto rectilinear grid with 1 arcsec pixels, sort by X */
-  /* reference for coords is this image */
-  tcoords = image[0].coords;
-  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
-  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
-  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
-  strcpy (tcoords.ctype, "RA---TAN");
-  
+  /* project onto rectilinear grid with 1 arcsec pixels */
+  /* we keep the original crpix1,2 and crref1,2 */
+  /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */
+  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
+    tcoords = MOSAIC[0];
+    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    tcoords.Npolyterms = 1;
+    strcpy (tcoords.ctype, "RA---TAN");
+  } else {
+    tcoords = image[0].coords;
+    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    strcpy (tcoords.ctype, "RA---TAN");
+  }
+
   for (i = 0; i < Nstars; i++) {
     fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords);
@@ -134 +144 @@
       /* calculate accurate per-star airmass */
       secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude);
-
       
       /** add measurements for this star **/
@@ -162 +171 @@
       if (*Mval == NO_MAG) *Mval = Mrel;
 
+      /** this is a little tricky: the new measures are not in the measure sequence
+          this call can run into the next star if we have had two measures found at this location. **/
       if (CALIBRATE) {
-        /** this is a little tricky: the new measures are not in the measure sequence
-            this call can run into the next star if we have had two measures found at this location. **/
         AddToCalibration (&catalog[0].average[n], &catalog[0].measure[m], &catalog[0].measure[Nmeas], N);
       }
@@ -202 +211 @@
 
   /* add reference for undetected catalog stars */
+  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) RegisterMosaic (MOSAIC);
   for (j = 0; (j < Nave) && !SKIP_MISSED; j++) {
     n = N2[j];
@@ -287 +297 @@
     /** now add references from all previous non-detection observations of this spot on the sky */
     for (j = 0; (j < Noverlap) && !SKIP_MISSED; j++) {
+      if (!FindMosaicForImage (overlap, Noverlap, j)) continue;
       if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue;
       if (catalog[0].average[Nave].Nn < 1) {

trunk/Ohana/src/addstar/src/gimages.c

@@ -38 +38 @@
   /* project onto rectilinear grid with 1 arcsec pixels */
   /* we keep the original crpix1,2 and crref1,2 */
-  tcoords = image[0].coords;
-  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
-  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
-  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
-  strcpy (tcoords.ctype, "RA---TAN");
+  /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */
+  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
+    tcoords = MOSAIC[0];
+    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    tcoords.Npolyterms = 1;
+    strcpy (tcoords.ctype, "RA---TAN");
+    RegisterMosaic (MOSAIC);
+    /* register so image->sky conversions below have correct mosaic */
+  } else {
+    tcoords = image[0].coords;
+    tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+    tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
+    tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    strcpy (tcoords.ctype, "RA---TAN");
+  }
 
   /* define original corners */
@@ -51 +63 @@
   Xo[4] = 0;           Yo[4] = 0;  /* so we can make a loop easily */
   for (j = 0; j < 5; j++) {
+    /* XY-to_RD is two-level if ctype == WRP */
     XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords);
     RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords);
@@ -77 +90 @@
   } 
 
-  NTIMAGE = 100;
+  /* kind of dangerous, see note below */
+  NTIMAGE = MAX (MIN (Nimage, 1000), 1);
   ALLOCATE (timage, Image, NTIMAGE);
 
@@ -87 +101 @@
   for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) {
     ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image");
+    if (ntimage == 0) Shutdown ("expected %d images, only found %d", Nimage, Ntimage);
+    BuildChipMatch (timage, ntimage);
+
     for (i = 0; i < ntimage; i++) {
+
+      /* if any of these images are WRP images, need to find matching DIS */
+      if (!FindMosaicForImage (timage, ntimage, i)) {
+        fprintf (stderr, "!");
+        continue;
+      }
+
       addtolist = FALSE;
       /* define image corners */
@@ -116 +140 @@
       }
     }
+  }
 
-    if (ntimage == 0) {
-      Shutdown ("expected %d images, only found %d", Nimage, Ntimage);
-    }
-  }
+  BuildChipMatch (pimage, npimage);
 
   if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage);
@@ -129 +151 @@
 }
   
+/* this function tries to find overlapping images with the current image.
+   it loads the image table in blocks (up to 1000 at a time) and looks
+   for overlaps within the set.  the search is not very robust. in 
+   addition, it needs to look for mosaic DIS images matching any WRP images.
+   if the WRP image is loaded near the beginning of a block, then the corresponding DIS 
+   may be in the previous block, and thus missed 
+*/

trunk/Ohana/src/addstar/src/gregion_image.c

@@ -12 +12 @@
   impatch.DEC[0] = +90;
   impatch.DEC[1] = -90;
+
+  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
+    if (MOSAIC == NULL) Shutdown ("no mosaic for WRP image (use -mosaic)");
+    RegisterMosaic (MOSAIC);
+  }    
 
   for (x = 0; x <= 1; x++) {

trunk/Ohana/src/addstar/src/gstars.c

@@ -52 +52 @@
     exit (1);
   }
+  if (!strcmp (&image[0].coords.ctype[4], "-WRP")) {
+    if (MOSAIC == NULL) {
+      fprintf (stderr, "ERROR: no mosaic for WRP image (use -mosaic)\n");
+      exit (1);
+    }
+    RegisterMosaic (MOSAIC);
+  }    
     
   /* get ST (used for airmass calculation) */