Changeset 23138

Timestamp:

Mar 3, 2009, 3:52:18 PM (17 years ago)

Author:

eugene

Message:

adding the rings mode (not quite done)

Location:

trunk/Ohana/src/addstar/src

Files:

• args_skycells.c (modified) (2 diffs)
• sky_tessalation.c (modified) (3 diffs)

trunk/Ohana/src/addstar/src/args_skycells.c

@@ -35 +35 @@
       MODE = LOCAL;
     }
+    if (!strcasecmp (argv[N], "rings")) {
+      MODE = RINGS;
+    }
     remove_argument (N, &argc, argv);
   }
@@ -140 +143 @@
 
   LEVEL = 8;
-  if ((MODE != LOCAL) && (N = get_argument (argc, argv, "-level"))) {
-    remove_argument (N, &argc, argv);
-    LEVEL = strtol (argv[N], &ptr, 10);
-    remove_argument (N, &argc, argv);
-    if ((*ptr != 0) || (LEVEL < 0)) {
-      fprintf (stderr, "-level requires an integer (>= 0) argument\n");
+  if ((MODE == SQUARES) || (MODE == TRIANGLES)) {
+    if ((N = get_argument (argc, argv, "-level"))) {
+      remove_argument (N, &argc, argv);
+      LEVEL = strtol (argv[N], &ptr, 10);
+      remove_argument (N, &argc, argv);
+      if ((*ptr != 0) || (LEVEL < 0)) {
+        fprintf (stderr, "-level requires an integer (>= 0) argument\n");
+        help ();
+      }  
+    }
+  }
+
+  CELLSIZE = 4.0;
+  if ((MODE == RINGS) && (N = get_argument (argc, argv, "-cellsize"))) {
+    remove_argument (N, &argc, argv);
+    CELLSIZE = strtod (argv[N], &ptr);
+    if ((*ptr != 0) || (CELLSIZE < 0.0)) {
+      fprintf (stderr, "-level requires a floating-point argument\n");
       help ();
     }  
+    remove_argument (N, &argc, argv);
   }
 

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

@@ -11 +11 @@
   sky_tessellation_init (scale);
 
-  if (mode == SQUARES) {
-    sky_tessellation_squares (db, level, Nmax);
-    return TRUE;
-  }
-
-  if (mode == TRIANGLES) {
-    sky_tessellation_triangles (db, level, Nmax);
-    return TRUE;
-  }
-
-  if (mode == LOCAL) {
-    sky_tessellation_local (db, level, Nmax);
-    return TRUE;
+  switch (mode) {
+    case SQUARES:
+      sky_tessellation_squares (db, level, Nmax);
+      return TRUE;
+    case TRIANGLES:
+      sky_tessellation_triangles (db, level, Nmax);
+      return TRUE;
+    case LOCAL:
+      sky_tessellation_local (db, level, Nmax);
+      return TRUE;
+    case RINGS:
+      sky_tessellation_rings (db, level, Nmax);
+      return TRUE;
+    default:
+      break;
   }
 
@@ -222 +224 @@
 
   free (image);
+  return (TRUE);
+}
+
+// the RINGS tessellation uses the declination zones proposed by Tamas Budavari
+// we generate projects on uniform rings of constant dec height
+int sky_tessellation_rings (FITS_DB *db, int level, int Nmax) {
+
+  int j, nDEC, Nimage, Nring;
+  float dec, dDEC;
+  SkyRectangle *ring;
+  Image *image;
+
+  // The tessellation has one input parameter: the approximate cell size.  Starting with
+  // the cell size, determine the optimal projection cell height (dDEC) that results in an
+  // integer number of dec zones between -90 and +90
+
+  nDEC = 180.0 / CELLSIZE;
+  dDEC = 180.0 / nDEC;
+
+  // generate the a collection of rectangles for each ring
+  for (dec = -90.0 + 0.5*dDEC; dec < +90.0; dec += dDEC) {
+
+    ring = sky_rectangle_ring (dec, dDEC, &Nring);
+
+    // subdivide each image (Nx x Ny subcells)
+    Nimage = NX_SUB*NY_SUB*Nring;
+    ALLOCATE (image, Image, Nimage);
+    for (j = 0; j < Nring; j++) {
+      // convert the SkyRectangles to Images for output
+      sky_subdivide_image (&image[j*NX_SUB*NY_SUB], &ring[j], NX_SUB, NY_SUB);
+    }
+
+    /* add the new images and save */
+    dvo_image_addrows (db, image, Nimage);
+    SetProtect (TRUE);
+    dvo_image_update (db, VERBOSE);
+    SetProtect (FALSE);
+    dvo_image_clear_vtable (db);
+    
+    free (ring);
+    free (image);
+  }    
   return (TRUE);
 }
@@ -480 +524 @@
 
   return (TRUE);
+}
+
+// define the parameters of a single sky projection center
+SkyRectangle *sky_rectangle_ring (float dec, float dDEC, int *nring) {
+
+  int i, nRA, NX, NY;
+  float dRA, decLower;
+  SkyRectangle *ring;
+
+  // 'dec' is a guess at the center of the cell; in fact, we need to choose decLower and
+  // decUpper to ensure complete overlap of the cells
+
+  // we can determine the 'lower' bound (bound closest to the equator):
+  decLower = (dec > 0.0) ? dec - 0.5*dDEC : dec + 0.5*dDEC;
+
+  // Subdivide the 'lower' bound into an integer number of segments:
+  nRA = cos(dec*RAD_DEG) * 360.0 / CELLSIZE; // CELLSIZE is a projection size
+  dRA = 360.0 / nRA;                         // dRA is a size in RA degrees
+
+  // I think we need to return the value of dec for the next ring, but I am not sure...
+
+  ALLOCATE (ring, SkyRectangle, nRA);
+
+  for (i = 0; i < nRA; i++) {
+    memset (&ring[i], 0, sizeof(SkyRectangle));
+    memset (&ring[i].coords, 0, sizeof(Coords));
+    ring[i].coords.crval1 = i*dRA;
+    ring[i].coords.crval2 = dec;
+
+    ring[i].coords.pc1_1 = +1.0;
+    ring[i].coords.pc1_2 = +0.0;
+    ring[i].coords.pc2_1 = -0.0;
+    ring[i].coords.pc2_2 = +1.0;
+  
+    // range values are in projected degrees
+    NX = cos(dec*RAD_DEG) * dRA  * 3600.0 / SCALE;
+    NY =                    dDEC * 3600.0 / SCALE;
+
+    // crpix1,crpix2 is the projection center
+    ring[i].coords.crpix1 = 0.5*NX;
+    ring[i].coords.crpix2 = 0.5*NY;
+
+    ring[i].coords.cdelt1 = SCALE / 3600.0;
+    ring[i].coords.cdelt2 = SCALE / 3600.0;
+
+    strcpy (ring[i].coords.ctype, "DEC--TAN");
+
+    ring[i].NX = NX;
+    ring[i].NY = NY;
+    ring[i].photcode = 1; // this needs to be set more sensibly
+
+
+    fprintf (stderr, "%f %f  : %f %f\n", 
+             ring[i].coords.crval1, ring[i].coords.crval2, 
+             ring[i].coords.crpix1, ring[i].coords.crpix2);
+  }
+
+  *nring = nRA;
+  return (ring);
 }