Index: /trunk/Ohana/src/libdvo/src/skyregion_gsc.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/skyregion_gsc.c	(revision 12391)
+++ /trunk/Ohana/src/libdvo/src/skyregion_gsc.c	(revision 12392)
@@ -128,5 +128,5 @@
     /* subdivide each zone */
     for (j = 0; j < Nzones; j++) {
-      if (DEBUG) fprintf (stderr, "zone: %d : %f - %f : %f - %f\n", j, zones[j].Rmin, zones[j].Rmax, zones[j].Dmin, zones[j].Dmax);
+      if (DEBUG >= 1) fprintf (stderr, "zone: %d : %f - %f : %f - %f\n", j, zones[j].Rmin, zones[j].Rmax, zones[j].Dmin, zones[j].Dmax);
       SkyTableL2fromZone (&L2, &L3, &L4, band, &zones[j], i);
     }
@@ -251,42 +251,66 @@
 }
 
-/* give a complete set of regions in a Dec band, subdivide into zones of equal-sized regions */
+/* given a complete set of regions in a Dec band, subdivide into zones of equal-sized regions.
+ * this function assumes the DEC band is sliced first by lines of constant RA, then subdivided
+ * with lines of constant DEC.  Groups ('zones') of these smaller boxes have the same size in
+ * delta-Dec.  We are trying to fine the boundaries of these zones.  The initial assumption
+ * breaks down for the pole region (86.25 < DEC < 90), where the entire range is divided into a
+ * single circular patch.  we need to treat it differently.
+ */
+
 SkyRegionZone *SkyRegionFindZones (SkyTable *band, int *Nzones, int parent) {
   
-  float Dmin, Dmax, dDec;
-  int i, Nz, NZ, Nregions;
+  float Dmin, Dmax, dDec, dDecZone;
+  int i, j, Nz, NZ, Nregions, poleRegion;
   SkyRegion *regions;
+  SkyRegion tempregion;
   SkyRegionZone *zones;
+  char basename[64];
   
   regions = band[0].regions;
   Nregions = band[0].Nregions;
 
-  Nz = 0;
   NZ = 10;
   ALLOCATE (zones, SkyRegionZone, NZ);
  
-  zones[0].dDec = regions[0].Dmax - regions[0].Dmin;
-  zones[0].Rmin = regions[0].Rmin;
-  zones[0].L3start = 0;
-  zones[0].L1band = parent;
-  zones[0].Nset = NsetForDecRange (zones[0].dDec);
+  /* the pole regions need to be separately handled.  skip in analysis below */
+  poleRegion = 0;
 
   /* search for transitions in the region dDec, find the max Dec range */
-  Dmin = regions[0].Dmin;
-  Dmax = regions[0].Dmax;
+  Nz = -1;
+  dDecZone = -1;
+  Dmin = +90;
+  Dmax = -90;
   for (i = 0; i < Nregions; i++) {
+    if ((regions[i].Dmin >= +86.00) || (regions[i].Dmax <= -86.00)) {
+      /* drop the pole region (re-generated below) */
+      poleRegion = SIGN(regions[i].Dmin);
+      Nregions --;
+      tempregion = regions[i];
+      for (j = i; j < Nregions; j++) {
+	regions[j] = regions[j+1];
+      }
+      i--;
+      continue;
+    }
+    /* is this region in the current zone? */
     Dmin = MIN (regions[i].Dmin, Dmin);
     Dmax = MAX (regions[i].Dmax, Dmax);
     dDec = regions[i].Dmax - regions[i].Dmin;
     if (fabs(dDec - zones[Nz].dDec) > 0.001) {
-      /* we've found the end of the current zone */
-      zones[Nz].L3end = i;
-      zones[Nz].Rmax = regions[i-1].Rmax;
+      /* we've found the end of the current zone; set the ending info */
+      if (Nz >= 0) {
+	zones[Nz].L3end = i;
+	zones[Nz].Rmax = regions[i-1].Rmax;
+      }
+
+      /* go to the next zone */
       Nz++;
       CHECK_REALLOCATE (zones, SkyRegionZone, NZ, Nz, 10);
 
       /* start info for the new zone */
+      dDecZone = dDec;
       zones[Nz].Rmin = regions[i].Rmin;
-      zones[Nz].dDec = dDec;
+      zones[Nz].dDec = dDecZone;
       zones[Nz].L3start = i;
       zones[Nz].L1band = parent;
@@ -301,4 +325,56 @@
     zones[i].Dmin = Dmin;
     zones[i].Dmax = Dmax;
+  }
+  
+  /* the pole region is mis-allocated in the gsc table to L3.  elevate it 
+     to L2 (zone) and subdivide */
+  if (poleRegion) {
+    CHECK_REALLOCATE (zones, SkyRegionZone, NZ, Nz, 10);
+    /* pole region @ L2 */
+    zones[Nz].Rmin = 0.0;
+    zones[Nz].Rmax = 360.0;
+    if (poleRegion > 0) {
+      zones[Nz].Dmin = +86.25;
+      zones[Nz].Dmax = +90.00;
+      strcpy (basename, "n8230/pole");
+    } else {
+      zones[Nz].Dmin = -90.00;
+      zones[Nz].Dmax = -86.25;
+      strcpy (basename, "s8230/pole");
+    }
+    zones[Nz].dDec = 90.00 - 86.25;
+    zones[Nz].L3start = Nregions;
+    zones[Nz].L3end = Nregions + 5;
+    zones[Nz].L1band = parent;
+
+    band[0].Nregions += 4;
+    REALLOCATE (band[0].regions, SkyRegion, band[0].Nregions);
+    for (i = 0; i < 4; i++) {
+      band[0].regions[Nregions + i].Rmin = i * 90.0;
+      band[0].regions[Nregions + i].Rmax = i * 90.0 + 90.0;
+      if (poleRegion > 0) {
+	band[0].regions[Nregions + i].Dmin = +86.250;
+	band[0].regions[Nregions + i].Dmax = +88.125;
+      } else {
+	band[0].regions[Nregions + i].Dmin = -88.125;
+	band[0].regions[Nregions + i].Dmax = -86.250;
+      }
+      sprintf (band[0].regions[Nregions + i].name, "%s.%d", basename, i);
+    }
+    band[0].regions[Nregions + i].Rmin = 0.0;
+    band[0].regions[Nregions + i].Rmax = 360.0;
+    if (poleRegion > 0) {
+      band[0].regions[Nregions + i].Dmin = +88.125;
+      band[0].regions[Nregions + i].Dmax = +90.000;
+    } else {
+      band[0].regions[Nregions + i].Dmin = -90.000;
+      band[0].regions[Nregions + i].Dmax = -88.125;
+    }
+    sprintf (band[0].regions[Nregions + i].name, "%s.%d", basename, i);
+
+    zones[0].Nset = 6;
+    zones[1].Nset = 5;
+
+    Nz ++;
   }
 
@@ -353,5 +429,5 @@
     *p = '/';
     strcpy (L2[0].regions[Nr].name, name);
-    if (DEBUG) SkyRegionPrint (&L2[0].regions[Nr]);
+    if (DEBUG >= 2) SkyRegionPrint (&L2[0].regions[Nr]);
 
     /* childS and childE are set in SkyTableL3fromL2 */
@@ -389,5 +465,5 @@
     L3[0].regions[Nr].childE   =  0;
 
-    if (DEBUG) SkyRegionPrint (&L3[0].regions[Nr]);
+    if (DEBUG >= 3) SkyRegionPrint (&L3[0].regions[Nr]);
     /* name is set for the band in SkyRegionForDecBand */
 
@@ -413,4 +489,6 @@
   L3[0].childS = Nr;
   L3[0].childE = L4[0].Nregions;
+
+  // XXX handle the pole regions just a bit differently...
 
   /* subdivide L3 into NDIV boxes */
@@ -438,5 +516,5 @@
       sprintf (name, "%s.%02d", L3[0].name, Nbox);
       strcpy (L4[0].regions[Nr].name, name);
-      if (DEBUG) SkyRegionPrint (&L4[0].regions[Nr]);
+      if (DEBUG >= 4) SkyRegionPrint (&L4[0].regions[Nr]);
 
       Nr ++;