Index: /trunk/Ohana/src/libdvo/Makefile
===================================================================
--- /trunk/Ohana/src/libdvo/Makefile	(revision 5440)
+++ /trunk/Ohana/src/libdvo/Makefile	(revision 5441)
@@ -40,5 +40,8 @@
 $(SRC)/dvo_convert_elixir.$(ARCH).o    \
 $(SRC)/dvo_convert_loneos.$(ARCH).o    \
-$(SRC)/dvo_convert_panstarrs.$(ARCH).o
+$(SRC)/dvo_convert_panstarrs.$(ARCH).o \
+$(SRC)/skyregion_io.$(ARCH).o    \
+$(SRC)/skyregion_gsc.$(ARCH).o    \
+$(SRC)/skyregion_ops.$(ARCH).o
 
 AOBJ = \
Index: /trunk/Ohana/src/libdvo/include/dvo.h
===================================================================
--- /trunk/Ohana/src/libdvo/include/dvo.h	(revision 5440)
+++ /trunk/Ohana/src/libdvo/include/dvo.h	(revision 5441)
@@ -106,4 +106,15 @@
 } FITS_DB;
 
+/* SkyRegion : better implementation than GSCRegion */
+typedef struct {
+  int Nregions;
+  SkyRegion *regions;
+} SkyTable;
+
+typedef struct {
+  int Nregions;
+  SkyRegion **regions;
+} SkyList;
+
 # if (0)
 /* structure for data on a catalog region */
@@ -192,4 +203,5 @@
 int FindMosaicForImage_MatchSearch (Image *images, int Nimages, int entry);
 int BuildChipMatch (Image *images, int Nimages);
+void SetImageCorners (double *X, double *Y, Image *image);
 
 short int putMi (double value);
@@ -322,3 +334,17 @@
 ImagePanstarrs *ImageInternalToPanstarrs (Image *in, int Nvalues);
 
+/* skyregion APIs */
+int SkyTableSave (SkyTable *table, char *filename);
+SkyTable *SkyTableLoad (char *filename, int VERBOSE);
+SkyTable *SkyTableFromGSC (char *filename, int depth, int VERBOSE);
+SkyRegion *SkyRegionByPoint (SkyTable *table, int depth, double ra, double dec);
+SkyList *SkyListByPoint (SkyTable *table, double ra, double dec);
+SkyList *SkyListByRadius (SkyTable *table, int depth, double RA, double DEC, double radius);
+SkyList *SkyListByPatch (SkyTable *table, int depth, SkyRegion *patch);
+SkyList *SkyListByImage (SkyTable *table, int depth, Image *image, Coords *mosaic);
+SkyList *SkyListByBounds (SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax);
+SkyList *SkyListChildrenByBounds (SkyTable *table, int No, int depth, double Rmin, double Rmax, double Dmin, double Dmax);
+int SkyListFree (SkyList *list);
+int SkyTableFree (SkyTable *table);
+
 # endif
Index: /trunk/Ohana/src/libdvo/src/skyregion_gsc.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/skyregion_gsc.c	(revision 5441)
+++ /trunk/Ohana/src/libdvo/src/skyregion_gsc.c	(revision 5441)
@@ -0,0 +1,197 @@
+# include "dvo.h"
+# define NBANDS 24
+# define NDIV 4
+
+int DecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, 0,
+                  597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48, 0};
+
+double DecBands[] = {0.0, +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, +90.0,
+                     0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, -52.5, -60.0, -67.5, -75.0, -82.5, -90.0};
+
+char *DecNames[] = {
+  "n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", "n5230", "n6000", "n6730", "n7500", "n8230", "none",
+  "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", "s5230", "s6000", "s6730", "s7500", "s8230", "none",
+};
+
+SkyTable *SkyTableFromGSC (char *filename, int depth, int VERBOSE) {
+
+  int i, j, Nx, Ny, No, Nr, NR, Nlines;
+  int nx, ny, Ne, Ns, Nbox;
+  double RA0, RA1, DEC0, DEC1, dR, dD;
+  char temp[80], name[80];
+
+  Header theader;
+  FTable ftable;
+  SkyTable *skytable;
+  SkyRegion *regions;
+  FILE *f;
+  
+  f = fopen (filename, "r");
+  if (f == NULL) {
+    if (VERBOSE) fprintf (stderr, "can't find GSC Region file %s\n", filename);
+    return (NULL);
+  }
+
+  /* load in table data */
+  ftable.header = &theader;
+  if (!fits_fread_ftable (f, &ftable, "REGIONS")) {
+    if (VERBOSE) fprintf (stderr, "can't read GSC Region table\n");
+    fclose (f);
+    return (NULL);
+  }
+
+  fits_scan (ftable.header, "NAXIS1", "%d", 1, &Nx);
+  fits_scan (ftable.header, "NAXIS2", "%d", 1, &Ny);
+  
+  fits_free_header (&theader);
+
+  /* build supporting level 0 and 1 regions */
+  Nr = 0;
+  NR = 100;
+  ALLOCATE (regions, SkyRegion, 100);
+  
+  /* level 0 : full sky */
+  regions[Nr].Rmin 	=   0;
+  regions[Nr].Rmax 	= 360;
+  regions[Nr].Dmin 	= -90;
+  regions[Nr].Dmax 	= +90;
+  regions[Nr].index  =  0;
+  regions[Nr].depth  =  0;
+  regions[Nr].parent = -1;
+  regions[Nr].child  =  TRUE;
+  regions[Nr].table  =  (depth == 0);
+  strcpy (regions[Nr].name, "fullsky");
+  
+  No = Nr;
+  Nr ++;
+
+  /* level 1 : add the dec bands */
+  regions[No].childS = Nr;
+  /* first north */
+  for (i = 0; i < 12; i++, Nr++) {
+    regions[Nr].Rmin   =   0;
+    regions[Nr].Rmax   = 360;
+    regions[Nr].Dmin   = DecBands[i];
+    regions[Nr].Dmax   = DecBands[i+1];
+    regions[Nr].index  =  i+1;
+    regions[Nr].depth  =  1;
+    regions[Nr].parent =  0;
+    regions[Nr].child  =  TRUE;
+    regions[Nr].table  =  (depth == 1);
+    strcpy (regions[Nr].name, DecNames[i]);
+  }
+  /* now south */
+  for (i = 0; i < 12; i++, Nr++) {
+    regions[Nr].Rmin   =   0;
+    regions[Nr].Rmax   = 360;
+    regions[Nr].Dmin   = DecBands[i+14];
+    regions[Nr].Dmax   = DecBands[i+13];
+    regions[Nr].index  =  i+1;
+    regions[Nr].depth  =  1;
+    regions[Nr].parent =  0;
+    regions[Nr].child  =  TRUE;
+    regions[Nr].table  =  (depth == 1);
+    strcpy (regions[Nr].name, DecNames[i+13]);
+  }
+  regions[No].childE = Nr;
+
+  /* level 2 : copy the data from the GSC Region files */
+  No = 1;
+  for (i = 0; i < 25; i++) {
+    if (i == 12) continue;
+
+    Nlines = 0;
+    for (j = 0; j < i; j++) Nlines += DecLines[j];
+
+    /* find all GSC Regions in this band */
+    regions[No].childS = Nr;
+    for (j = 0; j < DecLines[i]; j++) {
+      strncpy (temp, &ftable.buffer[(j + Nlines)*48], 48);
+      temp[49] = 0;
+      hstgsc_hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]);
+      if (DEC1 < DEC0) SWAP (DEC0, DEC1);
+
+      /* check that we are in correct Dec band */
+      if (DEC1 < regions[No].Dmin) {
+	fprintf (stderr, "error: table mis-match (1)\n");
+	fprintf (stderr, "line: %s\n", temp);
+	fprintf (stderr, "region %d: %f - %f\n", No, regions[No].Dmin, regions[No].Dmax);
+	return (NULL);
+      }
+      if (DEC0 > regions[No].Dmax) {
+	fprintf (stderr, "error: table mis-match (2)\n");
+	fprintf (stderr, "line: %s\n", temp);
+	fprintf (stderr, "region %d: %f - %f\n", No, regions[No].Dmin, regions[No].Dmax);
+	return (NULL);
+      }
+
+      /* regions with RA0 < 360 have RA1 == 0.0 */
+      if (RA1 < RA0) RA1 += 360.0;
+
+      regions[Nr].Dmin = DEC0;
+      regions[Nr].Dmax = DEC1;
+      regions[Nr].Rmin = RA0;
+      regions[Nr].Rmax = RA1;
+
+      regions[Nr].index  =  Nr;
+      regions[Nr].depth  =  2;
+      regions[Nr].parent =  No;
+      regions[Nr].child  =  TRUE;
+      regions[Nr].table  =  (depth == 2);
+      regions[Nr].childS =  0;
+      regions[Nr].childE =  0;
+
+      temp[5] = 0;
+      sprintf (name, "%s/%s", DecNames[i], &temp[1]);
+      strcpy (regions[Nr].name, name);
+      Nr ++;
+      CHECK_REALLOCATE (regions, SkyRegion, NR, Nr, 100);
+    }
+    regions[No].childE = Nr;
+    No ++;
+  }
+
+  /* subdivide level 2 into NDIV boxes */
+  Ns = regions[regions[0].childS].childS;
+  Ne = regions[regions[0].childE - 1].childE;
+  for (i = Ns; i < Ne; i++) {
+
+    Nbox = 0;
+    RA0 = regions[i].Rmin;
+    DEC0 = regions[i].Dmin;
+    dR = (regions[i].Rmax - regions[i].Rmin) / NDIV;
+    dD = (regions[i].Dmax - regions[i].Dmin) / NDIV;
+
+    regions[i].childS = Nr;
+    for (ny = 0; ny < NDIV; ny ++) {
+      for (nx = 0; nx < NDIV; nx ++) {
+	regions[Nr].Rmin = RA0  + (nx + 0)*dR;
+	regions[Nr].Rmax = RA0  + (nx + 1)*dR;
+	regions[Nr].Dmin = DEC0 + (ny + 0)*dD;
+	regions[Nr].Dmax = DEC0 + (ny + 1)*dD;
+
+	regions[Nr].index  =  Nr;
+	regions[Nr].depth  =  3;
+	regions[Nr].parent =  i;
+	regions[Nr].child  =  FALSE;
+	regions[Nr].table  =  (depth == 3);
+	regions[Nr].childS =  0;
+	regions[Nr].childE =  0;
+
+	temp[5] = 0;
+	sprintf (name, "%s.%02d", regions[i].name, Nbox);
+	strcpy (regions[Nr].name, name);
+
+	Nr ++;
+	Nbox ++;
+	CHECK_REALLOCATE (regions, SkyRegion, NR, Nr, 100);
+      }
+    }
+    regions[i].childE = Nr;
+  }
+
+  ALLOCATE (skytable, SkyTable, 1);
+  skytable[0].regions = regions;
+  skytable[0].Nregions = Nr;
+  return (skytable);
+}
Index: /trunk/Ohana/src/libdvo/src/skyregion_io.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/skyregion_io.c	(revision 5441)
+++ /trunk/Ohana/src/libdvo/src/skyregion_io.c	(revision 5441)
@@ -0,0 +1,54 @@
+# include "dvo.h"
+
+SkyTable *SkyTableLoad (char *filename, int VERBOSE) {
+
+  Header header;
+  Matrix matrix;
+  Header theader;
+  FTable ftable;
+  SkyTable *skytable;
+  FILE *f;
+  
+  f = fopen (filename, "r");
+  if (f == NULL) {
+    if (VERBOSE) fprintf (stderr, "can't find Sky Region file %s\n", filename);
+    return (NULL);
+  }
+
+  /* load in table data */
+  ftable.header = &theader;
+  if (!fits_fread_header (f, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read Sky Region header\n");
+    fclose (f);
+    return (NULL);
+  }
+  if (!fits_fread_matrix (f, &matrix, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read Sky Region matrix\n");
+    fits_free_header (&header);
+    fclose (f);
+    return (NULL);
+  }
+  if (!fits_fread_ftable (f, &ftable, "SKY_REGIONS")) {
+    if (VERBOSE) fprintf (stderr, "can't read Sky Region table\n");
+    fits_free_header (&header);
+    fits_free_matrix (&matrix);
+    fclose (f);
+    return (NULL);
+  }
+
+  ALLOCATE (skytable, SkyTable, 1);
+  skytable[0].regions = fits_table_get_SkyRegion (&ftable, &skytable[0].Nregions, NULL);
+  
+  fits_free_header (&header);
+  fits_free_matrix (&matrix);
+  fits_free_header (&theader);
+
+  return (skytable);
+}
+
+int SkyTableSave (SkyTable *table, char *filename) {
+
+  fprintf (stderr, "not yet defined\n");
+
+}
+
Index: /trunk/Ohana/src/libdvo/src/skyregion_ops.c
===================================================================
--- /trunk/Ohana/src/libdvo/src/skyregion_ops.c	(revision 5441)
+++ /trunk/Ohana/src/libdvo/src/skyregion_ops.c	(revision 5441)
@@ -0,0 +1,394 @@
+# include "dvo.h"
+
+/* each region is bounded on the sky by lines of constant RA & DEC.
+   the region represents a specified depth, and it may or may not have
+   children.  The start and end of the children in the array are childS and childE 
+   If the region at the given depth is populated with an object table, then 'object' is TRUE.
+   If the region at the given depth is populated with an image table, then 'image' is TRUE.
+   I have defined no accelerators other than the table hierarchy */
+
+/* find region which overlaps c at given depth (-1 : populated ) */
+SkyRegion *SkyRegionByPoint (SkyTable *table, int depth, double ra, double dec) {
+  
+  int i, Ns, Ne, No;
+  SkyRegion *region;
+
+  region = table[0].regions;
+
+  Ns = 0;
+  Ne = 1;
+
+  while (1) {
+    No = -1;
+    for (i = Ns; (i < Ne) && (i < table[0].Nregions); i++) {
+      if (ra  < region[i].Rmin) continue;
+      if (ra  > region[i].Rmax) continue;
+      if (dec < region[i].Dmin) continue;
+      if (dec > region[i].Dmax) continue;
+      No = i;
+      break;
+    }
+    if (No == -1) return (NULL);
+    if ((depth == -1) && (region[No].table)) return (&region[No]);
+    if (depth == region[No].depth) return (&region[No]);
+    if ((depth > region[No].depth) && !region[No].child) return (NULL);
+
+    /* need to check Ns, Ne, or guarantee valid range */
+    Ns = region[No].childS;
+    Ne = region[No].childE;
+  }
+}
+
+/* find regions at all levels which overlap c */
+SkyList *SkyListByPoint (SkyTable *table, double ra, double dec) {
+
+  int i, Ns, Ne, No, N, NREGIONS;
+  SkyList *list;
+  SkyRegion *region;
+  
+  N = 0;
+  NREGIONS = 10;
+  ALLOCATE (list, SkyList, 1);
+  ALLOCATE (list[0].regions,  SkyRegion *, NREGIONS);
+  list[0].Nregions = N;
+
+  region = table[0].regions;
+
+  Ns = 0;
+  Ne = 1;
+
+  while (1) {
+    No = -1;
+    for (i = Ns; (i < Ne) && (i < table[0].Nregions); i++) {
+      if (ra  < region[i].Rmin) continue;
+      if (ra  > region[i].Rmax) continue;
+      if (dec < region[i].Dmin) continue;
+      if (dec > region[i].Dmax) continue;
+      No = i;
+      break;
+    }
+    if (No == -1) return (list);
+
+    list[0].regions[N] = &region[No];
+    N++;
+    if (N == NREGIONS) {
+      NREGIONS += 10;
+      REALLOCATE (list[0].regions, SkyRegion *, NREGIONS);
+    }      
+    list[0].Nregions = N;
+
+    /* need to check Ns, Ne, or guarantee valid range */
+    if (region[No].child) {
+      Ns = region[No].childS;
+      Ne = region[No].childE;
+    } else {
+      return (list);
+    }
+  }
+}
+
+SkyList *SkyListByRadius (SkyTable *table, int depth, double RA, double DEC, double radius) {
+
+  double rad;
+  double Rmin, Rmax, Dmin, Dmax;
+  SkyList *list;
+
+  Dmin = DEC - radius;
+  Dmax = DEC + radius;
+
+  if ((Dmin <= -89) || (Dmax >= 89)) {
+    Rmin = 0;
+    Rmax = 360;
+  } else {
+    rad = MAX (radius / (cos(Dmin*RAD_DEG)), radius / (cos(Dmax*RAD_DEG)));
+    Rmin = RA - rad;
+    Rmax = RA + rad;
+  }
+
+  list = SkyListByBounds (table, depth, Rmin, Rmax, Dmin, Dmax);
+  return (list);
+}
+
+SkyList *SkyListByPatch (SkyTable *table, int depth, SkyRegion *patch) {
+
+  SkyList *list;
+
+  list = SkyListByBounds (table, depth, patch[0].Rmin, patch[0].Rmax, patch[0].Dmin, patch[0].Dmax);
+  return (list);
+}
+
+/* user must be careful about mosaic registration */
+SkyList *SkyListByImage (SkyTable *table, int depth, Image *image, Coords *mosaic) {
+
+  int i, j;
+  SkyList *list;
+  double r, d, X[4], Y[4];
+  double Rmin, Rmax, Dmin, Dmax;
+  
+  RegisterMosaic (mosaic);
+  SetImageCorners (X, Y, image);
+
+  Rmin = 360.0;
+  Rmax =   0.0;
+  Dmin = +90.0;
+  Dmax = -90.0;
+
+  /* does this work at 0,360 boundary? XY_to_RD must return 
+     coord offsets relative to CRVAL1,2 (ie, NOT renormalize) */
+  for (j = 0; j < 4; j++) {
+    /* XY_to_RD is two-level if ctype == WRP */
+    XY_to_RD (&r, &d, X[j], Y[j], &image[0].coords);
+    Rmin = MIN (Rmin, r);
+    Rmax = MAX (Rmax, r);
+    Dmin = MIN (Dmin, d);
+    Dmax = MAX (Dmax, d);
+  }
+
+  list = SkyListByBounds (table, depth, Rmin, Rmax, Dmin, Dmax);
+  return (list);
+}
+
+SkyList *SkyListByBounds (SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax) {
+
+  int i, Ns;
+  SkyList *list, *extra;
+
+  /* rationalize Rmin, Rmax */
+  while (Rmin <   0.0) Rmin += 360;
+  while (Rmin > 360.0) Rmin -= 360;
+  while (Rmax <   0.0) Rmax += 360;
+  while (Rmax > 360.0) Rmax -= 360;
+
+  /* handle 0,360 boundary requests */
+  /* this is probably wrong: I will get duplicates for all Dec bands... */
+  if (Rmin > Rmax) {
+    list = SkyListChildrenByBounds (table, -1, depth, 0.0, Rmax, Dmin, Dmax);
+
+    extra = SkyListChildrenByBounds (table, -1, depth, Rmin, 360.0, Dmin, Dmax);
+    REALLOCATE (list[0].regions, SkyRegion *, list[0].Nregions + extra[0].Nregions);
+    Ns = list[0].Nregions;
+    for (i = 0; i < extra[0].Nregions; i++) {
+      list[0].regions[i + Ns] = extra[0].regions[i];
+    }
+    list[0].Nregions += extra[0].Nregions;
+    SkyListFree (extra);
+  } else {
+    list = SkyListChildrenByBounds (table, -1, depth, Rmin, Rmax, Dmin, Dmax);
+  }
+
+  return (list);
+}
+
+SkyList *SkyListChildrenByBounds (SkyTable *table, int No, int depth, double Rmin, double Rmax, double Dmin, double Dmax) {
+
+  int i, j, Ns, Ne, Nnew, NNEW;
+  int append;
+  SkyList *children;
+  SkyList *list;
+  SkyRegion *region;
+
+  Nnew = 0;
+  NNEW = 50;
+  ALLOCATE (list, SkyList, 1);
+  ALLOCATE (list[0].regions, SkyRegion *, NNEW);
+
+  region = table[0].regions;
+
+  if (No == -1) {
+    Ns = 0;
+    Ne = 1;
+  } else {
+    Ns = region[No].childS;
+    Ne = region[No].childE;
+  }
+
+  for (i = Ns; (i < Ne) && (i < table[0].Nregions); i++) {
+    if (Rmax < region[i].Rmin) continue;
+    if (Rmin > region[i].Rmax) continue;
+    if (Dmax < region[i].Dmin) continue;
+    if (Dmin > region[i].Dmax) continue;
+
+    if ((depth > region[i].depth) && (region[i].child == FALSE)) continue;
+
+    append = FALSE;
+    append |= (depth > region[i].depth) && (region[i].child == TRUE);
+    append |= (depth ==             -1) && (region[i].table == FALSE);
+
+    if (append) {
+      /* append children to new */
+      children = SkyListChildrenByBounds (table, i, depth, Rmin, Rmax, Dmin, Dmax);
+      if (Nnew + children[0].Nregions >= NNEW) {
+	NNEW = Nnew + children[0].Nregions + 50;
+	REALLOCATE (list[0].regions, SkyRegion *, NNEW);
+      }
+      for (j = 0; j < children[0].Nregions; j++) {
+	list[0].regions[Nnew + j] = children[0].regions[j];
+      }
+      Nnew += children[0].Nregions;
+      SkyListFree (children);
+    } else {
+      list[0].regions[Nnew] = &region[i];
+      Nnew ++;
+      CHECK_REALLOCATE (list[0].regions, SkyRegion *, NNEW, Nnew, 50);
+    }
+  }
+
+  list[0].Nregions = Nnew;
+  return (list);
+}
+
+/* XXX EAM : this should go elsewhere (libdvo?) */
+void SetImageCorners (double *X, double *Y, Image *image) {
+
+  if (!strcmp(&image[0].coords.ctype[4], "-DIS")) {
+    X[0] = -0.5*image[0].NX; Y[0] = -0.5*image[0].NY;
+    X[1] = +0.5*image[0].NX; Y[1] = -0.5*image[0].NY;
+    X[2] = +0.5*image[0].NX; Y[2] = +0.5*image[0].NY;
+    X[3] = -0.5*image[0].NX; Y[3] = +0.5*image[0].NY;
+  } else {
+    X[0] = 0;           Y[0] = 0;
+    X[1] = image[0].NX; Y[1] = 0;
+    X[2] = image[0].NX; Y[2] = image[0].NY;
+    X[3] = 0;           Y[3] = image[0].NY;
+  }
+}
+
+int SkyListFree (SkyList *list) {
+
+  if (list == NULL) return (TRUE);
+  if (list[0].regions != NULL) {
+    free (list[0].regions);
+  }
+  free (list);
+  list = NULL;
+  return (TRUE);
+}
+
+int SkyTableFree (SkyTable *table) {
+
+  if (table == NULL) return (TRUE);
+  if (table[0].regions != NULL) {
+    free (table[0].regions);
+  }
+  free (table);
+  table = NULL;
+  return (TRUE);
+}
+
+# if (0)
+
+/* find regions contained within rectangular region  c1 - c2 */
+SkyRegion **SkyFindArea (SkyRegion *db, SkyCoord c1, SkyCoord c2, int *nlist) {
+
+  int Nlist;
+  SkyRegion *list, *new, *sub;
+
+  while (c1.r > 360.0) c1.r -= 360.0;
+  while (c1.r <   0.0) c1.r += 360.0;
+  while (c2.r > 360.0) c2.r -= 360.0;
+  while (c2.r <   0.0) c2.r += 360.0;
+
+  /* check on c1.r > c2.r : cross boundary */
+  if (c1.r > c2.r) {
+    c0 = c2; c0.r = 360.0;
+    list = SkyFindArea (db, c1, c0, &Nlist);
+    c0 = c1; c0.r = 0.0;
+    new  = SkyFindArea (db, c0, c2, &Nnew);
+    REALLOCATE (list, SkyRegion *, MAX (1, Nlist + Nnew));
+    memcpy (&list[Nlist], new, Nnew*sizeof(SkyRegion *));
+    free (new);
+    Nlist += Nnew;
+    *nlist = Nlist;
+    return (list);
+  }    
+
+  Ns = 0;
+  Ne = 1;
+  Nlist = 1;
+  ALLOCATE (list, SkyRegion *, Nlist);
+  list[0] = &region[0];
+  getchild = region[0].child;
+
+  while (getchild) {
+
+    getchild = FALSE;
+    Nnew = 0;
+    NNEW = 100;
+    ALLOCATE (new, SkyRegion *, NNEW);
+
+    for (i = Ns; i < Ne; i++) {
+      children = SkyFindChildren(db, list[i], c1, c2, &Nchildren);
+      if (children == NULL) continue;
+      if (Nnew + Nchildren >= NNEW) {
+	NNEW += 100;
+	REALLOCATE (new, SkyRegion *, NNEW);
+      }
+      for (i = 0; i < Nchildren; i++) {
+	getchild |= children[i][0].child;
+	new[Nnew] = children[i];
+	Nnew++;
+      }
+      free (children);
+    }
+
+    REALLOCATE (list, SkyRegion *, Nlist + Nnew);
+    memcpy (&list[Nlist], new, Nnew*sizeof(SkyRegion *));
+    free (new);
+    Nlist += Nnew;
+  }
+  return (list);
+  *nlist = Nlist;
+}
+
+/* find all children of the given sky region */
+SkyRegion **SkyFindChildren (SkyRegion *db, SkyRegion *parent, SkyCoord c1, SkyCoord c2, int *Nchildren) {
+
+  int i, Ns, Ne, Nregion, NREGIONS;
+  SkyRegion **region;
+  
+  *Nchildren = 0;
+  if (!parent[0].child) return (NULL);
+
+  Ns = parent[0].childS;
+  Ne = parent[0].childE;
+
+  Nregion = 0;
+  NREGIONS = 100;
+  ALLOCATE (region, SkyRegion *, NREGIONS);
+
+  for (i = Ns; i < Ne; i++) {
+    if (region[i].Rmax < c1.r) continue;
+    if (region[i].Rmin > c2.r) continue;
+    if (region[i].Dmax < c1.d) continue;
+    if (region[i].Dmin > c2.d) continue;
+    region[Nregion] = &region[i];
+    Nregion++;
+    if (Nregion == NREGIONS) {
+      NREGIONS += 100;
+      REALLOCATE (region, SkyRegion *, NREGIONS);
+    }      
+  }
+  *Nchildren = Nregion; 
+  return (region);
+}
+
+/* region is pointer to entry in db */
+SkyRegion *SkyExtend (SkyRegion *db, SkyRegion *region) {
+
+  fprintf (stderr, "not implemented\n");
+
+}
+
+/* region is pointer to entry in db */
+SkyRegion *SkyContract (SkyRegion *db, SkyRegion *region) {
+
+  fprintf (stderr, "not implemented\n");
+
+}
+
+SkyRegion *SkyFindGCircle (SkyRegion *db, SkyCoord c1, SkyCoord c2) {
+
+  fprintf (stderr, "not implemented\n");
+
+}
+# endif