Changeset 34844 for trunk/Ohana/src/addstar

Timestamp:

Dec 18, 2012, 10:04:35 AM (14 years ago)

Author:

eugene

Message:

merging changes from eam_branches/ipp-20121130

Location:

trunk/Ohana

Files:

• . (modified) (1 prop)
• src/addstar/src/findskycell.c (modified) (10 diffs)
• src/addstar/src/sky_tessalation.c (modified) (3 diffs)
• src/addstar/test/boundarytree.sh (copied) (copied from branches/eam_branches/ipp-20121130/Ohana/src/addstar/test/boundarytree.sh )

trunk/Ohana/src/addstar/src/findskycell.c

@@ -15 +15 @@
 
 void usage (void) {
-  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir)\n");
+  fprintf (stderr, "USAGE: findcell -mktree (tree) (catdir) [-nx Nx] [-ny Ny]\n");
   fprintf (stderr, "USAGE: findcell -tree (tree) (datafile)\n");
   fprintf (stderr, "   (datafile) should contain a list of RA,DEC pairs\n");
@@ -24 +24 @@
 int apply_tree (char *treefile, char *datafile);
 
+float SCALE = 1.0;
+int NX_SUB = 1;
+int NY_SUB = 1;
+
 int main (int argc, char **argv) {
 
@@ -37 +41 @@
   if (get_argument (argc, argv, "-help")) usage ();
   if (get_argument (argc, argv, "-h")) usage ();
+
+  NX_SUB = NY_SUB = 1;
+  if ((N = get_argument (argc, argv, "-nx"))) {
+    remove_argument (N, &argc, argv);
+    NX_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-ny"))) {
+    remove_argument (N, &argc, argv);
+    NY_SUB = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* pixel scale (arcsec/pixel) */
+  SCALE = 1.0;
+  if ((N = get_argument (argc, argv, "-scale"))) {
+    remove_argument (N, &argc, argv);
+    SCALE = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   /* extra error messages */
@@ -109 +133 @@
 
   tree.Nzone = 46;
+  tree.NX_SUB = NX_SUB;
+  tree.NY_SUB = NY_SUB;
+  tree.dPix = SCALE/3600.0;
 
   ALLOCATE (tree.Nband, int, tree.Nzone);
@@ -122 +149 @@
   ALLOCATE (tree.ra,   double *, tree.Nzone);
   ALLOCATE (tree.dec,  double *, tree.Nzone);
+  ALLOCATE (tree.Xo,   double *, tree.Nzone);
+  ALLOCATE (tree.Yo,   double *, tree.Nzone);
+  ALLOCATE (tree.dX,      int *, tree.Nzone);
+  ALLOCATE (tree.dY,      int *, tree.Nzone);
   ALLOCATE (tree.cell,    int *, tree.Nzone);
   ALLOCATE (tree.name,  char **, tree.Nzone);
@@ -152 +183 @@
     ALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
     ALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
+    ALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+    ALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
+    ALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
     ALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
     ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
@@ -157 +192 @@
       tree.ra[zone][band] = NAN;
       tree.dec[zone][band] = NAN;
+      tree.Xo[zone][band] = NAN;
+      tree.Yo[zone][band] = NAN;
+      tree.dX[zone][band] = -1;
+      tree.dY[zone][band] = -1;
       tree.cell[zone][band] = -1;
       ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH);
@@ -179 +218 @@
       REALLOCATE (tree.ra[zone],   double, tree.NBAND[zone]);
       REALLOCATE (tree.dec[zone],  double, tree.NBAND[zone]);
-      REALLOCATE (tree.cell[zone], int,    tree.NBAND[zone]);
+      REALLOCATE (tree.Xo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.Yo[zone],   double, tree.NBAND[zone]);
+      REALLOCATE (tree.dX[zone],      int, tree.NBAND[zone]);
+      REALLOCATE (tree.dY[zone],      int, tree.NBAND[zone]);
+      REALLOCATE (tree.cell[zone],    int, tree.NBAND[zone]);
       REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]);
       for (j = start; j < tree.NBAND[zone]; j++) {
         tree.ra[zone][j] = NAN;
         tree.dec[zone][j] = NAN;
+        tree.Xo[zone][band] = NAN;
+        tree.Yo[zone][band] = NAN;
+        tree.dX[zone][band] = -1;
+        tree.dY[zone][band] = -1;
         tree.cell[zone][j] = -1;
         ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH);
@@ -190 +237 @@
     tree.ra[zone][band] = ra;
     tree.dec[zone][band] = dec;
+    tree.Xo[zone][band] = x;
+    tree.Yo[zone][band] = y;
+    tree.dX[zone][band] = image[i].NX / NX_SUB;
+    tree.dY[zone][band] = image[i].NY / NY_SUB;
+    
     tree.cell[zone][band] = i;
 
@@ -325 +377 @@
     }
 
-    fprintf (stdout, "%10.6f %10.6f  %3d %3d  %s\n", ra, dec, zone, band, tree->name[zone][band]);
+    // I know the projection cell (band,zone), but I need to find the skycell within that projection cell.
+    // the proj cell is divided into Nx, Ny bits.
+    // (ra,dec) for (Ro,Do) -> (x,y).  given (Xo,Yo),(dX,dY) I can find ix,iy
+    // I currently track Ro,Do (tree->ra[zone][band], tree->dec[zone][band])
+    // I need to have the skycell center in pixels (Xo,Yo) and the scale dX,dY
+    
+    // convert R,D to X,Y with hard-wired projection and scale, orientation?
+    
+    double x = 0.0;
+    double y = 0.0;
+    BoundaryTreeProjection (&x, &y, ra, dec, tree, zone, band);
+
+    int xi = x / tree->dX[zone][band];
+    int yi = y / tree->dY[zone][band];
+    
+    char format[24], skycellname[128];
+    int Ndigit = (int)(log10(tree->NX_SUB*tree->NY_SUB)) + 1 ;
+    snprintf (format, 24, "%s.%%0%dd", tree->name[zone][band], Ndigit);
+
+    int N = xi + tree->NX_SUB * yi;
+    snprintf (skycellname, 128, format, N);
+
+    fprintf (stdout, "%10.6f %10.6f  %8.3f %8.3f  %3d %3d  %s\n", ra, dec, x, y, zone, band, skycellname);
   }
 
   exit (0);
 }
+

trunk/Ohana/src/addstar/src/sky_tessalation.c

@@ -275 +275 @@
       // convert the SkyRectangles to Images for output
       sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB);
-      printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2);
+      // printf("%s %8.2f %8.2f\n", ring[j].name, ring[j].coords.crval1, ring[j].coords.crval2);
     }
 
@@ -691 +691 @@
     // fprintf (stdout, "%f %f  %f  %f  %f %f  %f %f  %f %f %f\n", dec, decUpper, dRA, arg, root1, root2, theta1, theta2, r1, r2, k*r2);
   }
-  fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
+  // fprintf (stdout, "%f %f  %f x %f (%d)\n", dec, decUpper, dRA, theta, nRA);
 
   // I think we need to return the value of dec for the next ring, but I am not sure...
@@ -776 +776 @@
       ring[N].coords.crval2 = (j == 0) ? dec / d2r : -dec / d2r;
 
-      printf(" \t %d   %25.20f   %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1);
+      // printf(" \t %d   %25.20f   %25.20f\n", i, ring[N].coords.crval2, ring[N].coords.crval1);
 
       ring[N].coords.pc1_1 = +1.0 * X_PARITY;