Changeset 3389

Timestamp:

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

Author:

eugene

Message:

adding mosaic astrometry concepts

Location:

trunk/Ohana/src

Files:

• addstar/Makefile (modified) (2 diffs)
• addstar/include/addstar.h (modified) (1 diff)
• addstar/src/args.c (modified) (5 diffs)
• addstar/src/find_matches.c (modified) (5 diffs)
• addstar/src/gimages.c (modified) (6 diffs)
• addstar/src/gregion_image.c (modified) (1 diff)
• addstar/src/gstars.c (modified) (1 diff)
• libohana/Makefile (modified) (1 diff)
• libohana/doc/ChangeLog.txt (modified) (1 diff)
• libohana/doc/coordtrans.txt (modified) (1 diff)
• libohana/include/loneos.h (modified) (2 diffs)
• libohana/src/coordops.c (modified) (8 diffs)
• libohana/src/mosaic_astrom.c (added)

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) */

trunk/Ohana/src/libohana/Makefile

@@ -28 +28 @@
 $(SRC)/Fread.$(ARCH).o          $(SRC)/glockfile.$(ARCH).o       \
 $(SRC)/LoadPhotcodes.$(ARCH).o  $(SRC)/photfits.$(ARCH).o        \
-$(SRC)/phot_catalog.$(ARCH).o   $(SRC)/imreg_datatypes.$(ARCH).o
+$(SRC)/phot_catalog.$(ARCH).o   $(SRC)/imreg_datatypes.$(ARCH).o \
+$(SRC)/mosaic_astrom.$(ARCH).o
 
 $(DESTLIB)/libohana.a: $(LIB)/libohana.$(ARCH).a

trunk/Ohana/src/libohana/doc/ChangeLog.txt

@@ +1 @@
+
+  changed meaning of PLY / WRP / DIS (see coordtrans.txt)
+  added check on ALLOCATE/REALLOCATE to avoid failures on Nvalues < 1
+
+libohana-1-2
+  added 'STG' to allowed coordinate transforms (inaccurate & stop-gap)
+  set default Npolyterms to 1
+  initialized filemode in glockfile.c
 
 libohana-1-1

trunk/Ohana/src/libohana/doc/coordtrans.txt

@@ +1 @@
+
+I have changed the meaning of coord type PLY from the historical
+meaning.  The only previous use that is in a DVO-style database is in
+the loneos data (all cfht data needed to go to linear terms).  these
+can be easily converted by changing from PLY to DIS.
+
+There are now three types of higher-order transformation which may be
+used depending on the circumstance.  The new naming scheme is:
+
+PLY : a cartesian transformation (ie, no projection implied)
+DIS : equivalent to TAN projection with up-to-3rd order polynomial
+WRP : a second-level cartesian transformation, implying the presence
+      of a DIS coord frame (must be registered with libohana function)
 
 { 

trunk/Ohana/src/libohana/include/loneos.h

@@ -410 +410 @@
 int  fRD_to_XY (float *x,  float *y,   double ra, double dec, Coords *coords);
 int  GetCoords (Coords *coords, Header *header);
+int  PutCoords (Coords *coords, Header *header);
+void RegisterMosaic (Coords *coords);
+void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch);
+
 short int putMi (double value);
 double getMi (short int value);
@@ -458 +462 @@
 void ParseColorTerms (char *terms, float *X, int *N);
 
-void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch);
-
 int get_image_type (char *name);
 char *get_type_name (int type);

trunk/Ohana/src/libohana/src/coordops.c

@@ -2 +2 @@
 # include <loneos.h>
 
+static Coords *mosaic = NULL;
+
+void RegisterMosaic (Coords *coords) {
+  mosaic = coords;
+}
+
 int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) {
 
-  int Polynomial, Zenithal, Cartesian, PseudoCyl;
+  int Zenith1, Zenith2, Zenithal, Polynomial, Cartesian, PseudoCyl;
   char *type;
   double L, M, X, Y, T, Z;
@@ -15 +21 @@
   type = &coords[0].ctype[4];
   
-  /* DIS and PLY are identical & apply projection; WRP is cartesian */
+  /* DIS is equiv to TAN with higher order terms
+     PLY is cartesian with higher order terms
+     WRP is equiv to PLY, with implied mosaic */
+
   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
-  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
-  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+  Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
+  Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
+  Zenithal   = !strcmp(type, "-ZEA") || Zenith1 || Zenith2;
   if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
 
@@ -45 +56 @@
     *ra  = L + coords[0].crval1;
     *dec = M + coords[0].crval2;
+
+    /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */
+    if (!strcmp(type, "-WRP")) {
+      if (mosaic == NULL) return (FALSE);
+      XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic);
+    }
     return (TRUE);
   }
@@ -60 +77 @@
 
     /* this is wrong : STG is not TAN - need to put in correct relationships.  but is a close approx */
-    if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) {
+    if (Zenith1) {
       if (R == 0) {
         stht = 1.0;
@@ -70 +87 @@
       }
     }
-    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
+    if (Zenith2) {
       ctht = RAD_DEG * R;
       stht = sqrt (1 - ctht*ctht);
@@ -124 +141 @@
 
   char *type;
-  int i, status, Polynomial, Zenithal, Cartesian, PseudoCyl;
+  int i, status, Polynomial, Zenith1, Zenith2, Zenithal, Cartesian, PseudoCyl;
   double phi, theta;
   double determ;
@@ -139 +156 @@
 
   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
-  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-STG") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
-  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
+  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
+  Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
+  Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
+  Zenithal   = !strcmp(type, "-ZEA") || Zenith1 || Zenith2;
   if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
 
   /**** Locally Cartesian Projections ****/
   if (Cartesian) {
-    L = (ra  - coords[0].crval1);
-    M = (dec - coords[0].crval2);
+    if (!strcmp(type, "-WRP")) {
+      if (mosaic == NULL) return (FALSE);
+      RD_to_XY (&Lo, &Mo, ra, dec, mosaic);
+      L = (Lo - coords[0].crval1);
+      M = (Mo - coords[0].crval2);
+    } else {
+      L = (ra  - coords[0].crval1);
+      M = (dec - coords[0].crval2);
+    }
   }
   
@@ -164 +190 @@
     if (stht < 0) status = FALSE;
 
-    if (!strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG")) {
+    if (Zenith1) {
       L =  DEG_RAD * sphi / stht;
       M = -DEG_RAD * cphi / stht;
     }
-    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
+    if (Zenith2) {
       L =  DEG_RAD * sphi;
       M = -DEG_RAD * cphi;