Changeset 8328

Timestamp:

Aug 14, 2006, 1:14:22 PM (20 years ago)

Author:

eugene

Message:

modifying dvo_catalog APIs

Location:

trunk/Ohana/src

Files:

• addstar/src/addstar.c (modified) (2 diffs)
• addstar/src/find_matches.c (modified) (5 diffs)
• addstar/src/find_matches_closest.c (modified) (8 diffs)
• addstar/src/gcatalog.c (modified) (1 diff)
• addstar/src/load_pt_catalog.c (modified) (2 diffs)
• libdvo/doc/dvo-catalogs.txt (added)
• libdvo/include/dvo.h (modified) (4 diffs)
• libdvo/src/dvo_catalog.c (modified) (9 diffs)
• libdvo/src/dvo_catalog_create.c (added)
• libdvo/src/dvo_catalog_mef.c (modified) (2 diffs)
• libdvo/src/dvo_catalog_raw.c (modified) (2 diffs)
• libdvo/src/dvo_catalog_split.c (modified) (6 diffs)
• libdvo/src/skyregion_gsc.c (modified) (7 diffs)
• libdvo/src/skyregion_ops.c (modified) (3 diffs)
• opihi/dvo/skycat.c (modified) (2 diffs)
• relphot/src/load_catalogs.c (modified) (1 diff)

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

@@ -95 +95 @@
 
     catalog.filename = skylist[0].filename[i];
-    load_pt_catalog (&catalog, skylist[0].regions[i]);
+    catalog.catformat = dvo_catalog_format (CATFORMAT);  // set the default catformat from config data
+    catalog.catmode   = dvo_catalog_mode (CATMODE);      // set the default catmode from config data
+    catalog.lockmode  = LCK_XCLD;
+    dvo_catalog_open (&catalog, skylist[0].regions[i], Nsecfilt, "w");
 
     /* for only_match, skip empty catalogs XXX EAM : this leaves behind empty files */
     if ((catalog.Nave_disk == 0) && options.only_match) {
-      unlock_catalog (&catalog);
-      free (catalog.filename);
-      free_catalog (&catalog);
+      dvo_catalog_unlock (&catalog);
+      free (catalog.filename);  // can this be done within the dvo_catalog_free function?? 
+      dvo_catalog_free (&catalog);
       continue;
     }
@@ -136 +139 @@
 
     if (Nsubset == 0) {
+      // XXX remove empty catalogs
       unlock_catalog (&catalog);
     } else {
       SetProtect (TRUE);
-      if (!options.only_images) wcatalog (&catalog);
+      if (!options.only_images) dvo_catalog_save (&catalog);
       SetProtect (FALSE);
-      unlock_catalog (&catalog);
+      dvo_catalog_unlock (&catalog);
     }
-    free_catalog (&catalog);
-    free (catalog.filename);
+    dvo_catalog_free (&catalog);
+    free (catalog.filename); // XXX ????
 
     if (options.mode == M_REFCAT) free (stars);

trunk/Ohana/src/addstar/src/find_matches.c

@@ -1 +1 @@
 # include "addstar.h"
 
-void find_matches (SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options) {
-
-  int i, j, n, N, J;
+void find_matches (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options) {
+
+  int i, j, n, N, J, status, Nstars;
   double X, Y, RADIUS, RADIUS2, secz;
   float *X1, *Y1, *X2, *Y2;
@@ -23 +23 @@
 
   /** allocate local arrays (stars) **/
-  ALLOCATE (X1, float, Nstars);
-  ALLOCATE (Y1, float, Nstars);
-  ALLOCATE (N1, int,   Nstars);
+  ALLOCATE (X1, float, NstarsIn);
+  ALLOCATE (Y1, float, NstarsIn);
+  ALLOCATE (N1, int,   NstarsIn);
 
   /** allocate local arrays (catalog) **/
@@ -55 +55 @@
     tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
     tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    tcoords.Npolyterms = 1;
     strcpy (tcoords.ctype, "RA---TAN");
   }
 
   /* build spatial index (RA sort) */
-  for (i = 0; i < Nstars; i++) {
-    fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords);
-    N1[i] = i;
+  Nstars = 0;
+  for (i = 0; i < NstarsIn; i++) {
+    status = fRD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].R, stars[i].D, &tcoords);
+    if (!status) continue;
+    N1[Nstars] = i;
+    Nstars ++;
+  }
+  if (Nstars < 1) {
+    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
+    free (catalog[0].found);
+    free (X1);
+    free (Y1);
+    free (N1);
+    free (X2);
+    free (Y2);
+    free (N2);
+    return;
   }
   if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars);
@@ -333 +348 @@
 
   /* note stars which have been found in this catalog */
-  for (i = 0; i < Nstars; i++) {
+  for (i = 0; i < NstarsIn; i++) {
     if (stars[i].found > -1) {
       stars[i].found = -2;
@@ -346 +361 @@
   catalog[0].Nmissing = Nmiss;
   if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
+
+  free (catalog[0].found);
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (X2);
+  free (Y2);
+  free (N2);
   return;
 }

trunk/Ohana/src/addstar/src/find_matches_closest.c

@@ -1 +1 @@
 # include "addstar.h"
-# define DEBUG 1
-
-void find_matches_closest (SkyRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options) {
-
-  int i, j, n, N, J, Jmin;
+
+void find_matches_closest (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, Image *image, Image *overlap, int Noverlap, Coords *mosaic, AddstarClientOptions options) {
+
+  int i, j, n, N, J, Jmin, status, Nstars;
   double X, Y, RADIUS, RADIUS2, Rmin, secz;
   float *X1, *Y1, *X2, *Y2;
@@ -24 +23 @@
 
   /** allocate local arrays (stars) **/
-  ALLOCATE (X1, float, Nstars);
-  ALLOCATE (Y1, float, Nstars);
-  ALLOCATE (N1, int,   Nstars);
+  ALLOCATE (X1, float, NstarsIn);
+  ALLOCATE (Y1, float, NstarsIn);
+  ALLOCATE (N1, int,   NstarsIn);
 
   /** allocate local arrays (catalog) **/
@@ -40 +39 @@
   NMEAS = Nmeas = catalog[0].Nmeasure;
   NMISS = Nmiss = catalog[0].Nmissing;
-  
+
   /* project onto rectilinear grid with 1 arcsec pixels */
   /* we keep the original crpix1,2 and crref1,2 */
@@ -56 +55 @@
     tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
     tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
+    tcoords.Npolyterms = 1;
     strcpy (tcoords.ctype, "RA---TAN");
   }
 
-  /* build spatial index (RA sort) */
-  for (i = 0; i < Nstars; i++) {
-    fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords);
-    N1[i] = i;
+  /* build spatial index (RA sort) referencing input array sequence */
+  Nstars = 0;
+  for (i = 0; i < NstarsIn; i++) {
+    status = fRD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].R, stars[i].D, &tcoords);
+    if (!status) continue;
+    N1[Nstars] = i;
+    Nstars ++;
+  }
+  if (Nstars < 1) {
+    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
+    free (catalog[0].found);
+    free (X1);
+    free (Y1);
+    free (N1);
+    free (X2);
+    free (Y2);
+    free (N2);
+    return;
   }
   if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars);
-  
+
   /* build spatial index (RA sort) */
   for (i = 0; i < Nave; i++) {
@@ -149 +163 @@
     n = N2[Jmin];
     N = N1[i];
-
-    if (DEBUG) fprintf (stderr, "matched %f,%f and %f,%f\n", 
-                        catalog[0].average[n].R, catalog[0].average[n].D,
-                        stars[N].R, stars[N].D);
 
     /* add to end of measurement list */
@@ -329 +339 @@
   }
       
-  free (catalog[0].found);
   REALLOCATE (catalog[0].average, Average, Nave);
   REALLOCATE (catalog[0].measure, Measure, Nmeas);
@@ -344 +353 @@
 
   /* note stars which have been found in this catalog */
-  for (i = 0; i < Nstars; i++) {
+  for (i = 0; i < NstarsIn; i++) {
     if (stars[i].found > -1) {
       stars[i].found = -2;
@@ -356 +365 @@
   catalog[0].Nmeasure = Nmeas;
   catalog[0].Nmissing = Nmiss;
-  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", NstarsIn, Nave, Nmeas, Nmiss, Nmatch);
+
+  free (catalog[0].found);
+  free (X1);
+  free (Y1);
+  free (N1);
+  free (X2);
+  free (Y2);
+  free (N2);
   return;
 }

trunk/Ohana/src/addstar/src/gcatalog.c

@@ -8 +8 @@
 
   /* read catalog header */
-  if (!load_catalog (catalog, VERBOSE)) {
+  if (!dvo_catalog_load (catalog, VERBOSE)) {
     fprintf (stderr, "ERROR: failure loading catalog\n");
     exit (1);
   }
 
-  /* should this be moved into save_catalog?? */
-  status = gfits_scan (&catalog[0].header, "SORTED", "%t", 1, &catalog[0].sorted);
-  if (!status) catalog[0].sorted = TRUE;
-  /* XXX EAM - is this a good choice?  should the default be 'FALSE'? */
-
   /* check Nsecfile value, update if needed */
   Nsecfilt = GetPhotcodeNsecfilt ();
-  if (catalog[0].Nsecfilt < Nsecfilt) {
-
-    int i, j, Nextra, in, out;
-    SecFilt *insec, *outsec;
-
-    Nextra = Nsecfilt - catalog[0].Nsecfilt;
-    insec = catalog[0].secfilt;
-    ALLOCATE (outsec, SecFilt, catalog[0].Naverage * Nsecfilt);
-    for (in = out = i = 0; i < catalog[0].Naverage; i++) {
-      for (j = 0; j < catalog[0].Nsecfilt; j++, in++, out++) {
-        outsec[out].M_PS  = insec[in].M_PS;
-        outsec[out].dM_PS = insec[in].dM_PS;
-        outsec[out].Xm = insec[in].Xm;
-      }
-      for (j = 0; j < Nextra; j++, out++) {
-        outsec[out].M_PS  = NO_MAG;
-        outsec[out].dM_PS = NO_MAG;
-        outsec[out].Xm    = NO_MAG;
-      }
-    }
-    free (catalog[0].secfilt);
-    catalog[0].secfilt = outsec;
-    catalog[0].Nsecfilt = Nsecfilt;
-  }
-
-  if (catalog[0].Nsecfilt > Nsecfilt) {
+  if (!dvo_catalog_check (catalog, Nsecfilt, TRUE)) {
     fprintf (stderr, "ERROR: can't reduce number of secondary filters\n");
-    exit (1);
+    exit (2);
   }
   return (TRUE);

trunk/Ohana/src/addstar/src/load_pt_catalog.c

@@ -1 +1 @@
 # include "addstar.h"
 
+// XXX can I replace this with dvo_catalog_open ()
+// XXX set the catalog properties before the function call (not as args)
 int load_pt_catalog (Catalog *catalog, SkyRegion *region) {
   
-  if (!check_file_access (catalog[0].filename, TRUE, TRUE)) {
-    exit (1);
-  }
+  if (!check_file_access (catalog[0].filename, TRUE, TRUE)) exit (1);
 
   if (VERBOSE) fprintf (stderr, "adding to %s\n", catalog[0].filename);
     
+  /* check Nsecfile value, update if needed */
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
   switch (lock_catalog (catalog, LCK_XCLD)) {
   case 0:
@@ -14 +17 @@
     exit (1);
   case 1:
-    gcatalog (catalog); /* load from disk */
-    if (VERBOSE) fprintf (stderr, "loading existing file %s\n", catalog[0].filename);
+    if (!dvo_catalog_load (catalog, VERBOSE)) {
+      fprintf (stderr, "ERROR: failure loading catalog\n");
+      exit (1);
+    }
+    if (!dvo_catalog_check (catalog, Nsecfilt, TRUE)) {
+      fprintf (stderr, "ERROR: can't reduce number of secondary filters\n");
+      exit (2);
+    }
+    if (VERBOSE) fprintf (stderr, "loaded existing file %s\n", catalog[0].filename);
     break;
   case 2:
-    mkcatalog (region, catalog); /* fills in new header info */
+    dvo_catalog_create (region, catalog, Nsecfilt, CATFORMAT, CATMODE); /* fills in new header info */
     if (VERBOSE) fprintf (stderr, "creating new file %s\n", catalog[0].filename);
     break;

trunk/Ohana/src/libdvo/include/dvo.h

@@ -197 +197 @@
 void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch);
 
-int lock_catalog (Catalog *catalog, int lockmode);
-int unlock_catalog (Catalog *catalog);
-int load_catalog (Catalog *catalog, int VERBOSE);
-int save_catalog (Catalog *catalog, char VERBOSE);
-int update_catalog (Catalog *catalog, char VERBOSE);
-
 int FindMosaicForImage (Image *images, int Nimages, int entry);
 int FindMosaicForImage_TableSearch (Image *images, int Nimages, int entry);
@@ -221 +215 @@
 PhotCode *GetPhotcodebyNsec (int Nsec);
 PhotCode *GetPhotcodeEquivbyCode (int code);
-
 char     *GetPhotcodeNamebyCode (int code);
 
@@ -274 +267 @@
 int   ConvertStruct (char *buffer, int size, int Nbytes, char *type);
 
+/** dvo_catalog APIs */
+void dvo_catalog_init (Catalog *catalog);
+void dvo_catalog_create (SkyRegion *region, Catalog *catalog, int Nsecfilt, char *catformat, char *catmode);
+int dvo_catalog_lock (Catalog *catalog, int lockmode);
+int dvo_catalog_unlock (Catalog *catalog);
+int dvo_catalog_load (Catalog *catalog, int VERBOSE);
+int dvo_catalog_save (Catalog *catalog, char VERBOSE);
+int dvo_catalog_update (Catalog *catalog, char VERBOSE);
+
+/* catmode-specific APIs */
+int dvo_catalog_load_raw (Catalog *catalog, int VERBOSE);
+int dvo_catalog_save_raw (Catalog *catalog, char VERBOSE);
+int dvo_catalog_load_mef (Catalog *catalog, int VERBOSE);
+int dvo_catalog_save_mef (Catalog *catalog, char VERBOSE);
+int dvo_catalog_load_split (Catalog *catalog, int VERBOSE);
+int dvo_catalog_save_split (Catalog *catalog, char VERBOSE);
+int dvo_catalog_update_split (Catalog *catalog, char VERBOSE);
+
 /*** conversion functions / I/O conversions ***/
-
 Average *ReadRawAverage (FILE *f, int Naverage, int format);
 Measure *ReadRawMeasure (FILE *f, int Nmeasure, int format);
@@ -312 +322 @@
 int MeasureToFtable (FTable *ftable, Measure *measure, int Nmeasure, int format);
 int SecFiltToFtable (FTable *ftable, SecFilt *secfilt, int Nsecfilt, int format);
-
-int load_catalog_raw (Catalog *catalog, int VERBOSE);
-int save_catalog_raw (Catalog *catalog, char VERBOSE);
-int load_catalog_mef (Catalog *catalog, int VERBOSE);
-int save_catalog_mef (Catalog *catalog, char VERBOSE);
-int load_catalog_split (Catalog *catalog, int VERBOSE);
-int save_catalog_split (Catalog *catalog, char VERBOSE);
-int update_catalog_split (Catalog *catalog, char VERBOSE);
 
 /*** DVO image db I/O Functions ***/

trunk/Ohana/src/libdvo/src/dvo_catalog.c

@@ -1 +1 @@
 # include <ohana.h>
 # include <dvo.h>
+# define DEBUG 1
+
+void dvo_catalog_init (Catalog *catalog) {
+
+  catalog[0].f = NULL;
+  catalog[0].filename = NULL;
+
+  gfits_init_header (&catalog[0].header);
+
+  catalog[0].average = NULL;
+  catalog[0].measure = NULL; 
+  catalog[0].missing = NULL; 
+  catalog[0].secfilt = NULL;
+  
+  catalog[0].Naverage = 0;
+  catalog[0].Nmeasure = 0;
+  catalog[0].Nmissing = 0;
+  catalog[0].Nsecfilt = 0;
+
+  catalog[0].Nave_disk  = 0;
+  catalog[0].Nmeas_disk = 0;
+  catalog[0].Nmiss_disk = 0;
+  catalog[0].Nmeas_off  = 0;
+
+  /* pointers to SPLIT data files */
+  catalog[0].measure_catalog = NULL;
+  catalog[0].missing_catalog = NULL;
+  catalog[0].secfilt_catalog = NULL;
+
+  /* extra catalog information */
+  catalog[0].lockmode = 0;
+  catalog[0].catmode  = 0;
+  catalog[0].catformat = 0;
+  catalog[0].catflags = 0;
+  catalog[0].sorted = 0;
+  
+  /* pointers for data manipulation */
+  catalog[0].found = NULL;
+  catalog[0].image = NULL;
+  catalog[0].mosaic = NULL;
+  catalog[0].X = NULL;
+  catalog[0].Y = NULL;
+}
 
 /* possible exit status for lock_catalog: 
@@ -7 +50 @@
    2 - empty file (file may be open or closed!) 
 */
-
-int lock_catalog (Catalog *catalog, int lockmode) {
+int dvo_catalog_lock (Catalog *catalog, int lockmode) {
 
   int dbstate;
@@ -16 +58 @@
 
   /* set lock on database, create stream f */
+  // fprintf (stderr, "locking: %s\n", catalog[0].filename);
   catalog[0].f = fsetlockfile (catalog[0].filename, 3600.0, catalog[0].lockmode, &dbstate);
 
@@ -26 +69 @@
 }
 
-int unlock_catalog (Catalog *catalog) {
+int dvo_catalog_unlock (Catalog *catalog) {
 
   int dbstate;
 
   if (catalog[0].f == (FILE *) NULL) return (2);
+  // fprintf (stderr, "unlocking: %s\n", catalog[0].filename);
   fclearlockfile (catalog[0].filename, catalog[0].f, catalog[0].lockmode, &dbstate);
+
+  if (catalog[0].catmode == DVO_MODE_SPLIT) {
+    if (catalog[0].measure_catalog != NULL) dvo_catalog_unlock (catalog[0].measure_catalog);
+    if (catalog[0].missing_catalog != NULL) dvo_catalog_unlock (catalog[0].missing_catalog);
+    if (catalog[0].secfilt_catalog != NULL) dvo_catalog_unlock (catalog[0].secfilt_catalog);
+  }
   return (1);
 }
 
-int load_catalog (Catalog *catalog, int VERBOSE) {
+int dvo_catalog_open (Catalog *catalog, int NSECFILT, int VERBOSE, int MODE) {
+
+  if (!check_file_access (catalog[0].filename, TRUE, VERBOSE)) {
+    if (VERBOSE) fprintf (stderr, "no permission to access %f\n", catalog[0].filename);
+    return (FALSE);
+  }
+
+  if (VERBOSE) fprintf (stderr, "adding to %s\n", catalog[0].filename);
+    
+  switch (dvo_catalog_lock (catalog, catalog[0].lockmode)) {
+  case 0:
+    if (VERBOSE) fprintf (stderr, "can't lock file %s\n", catalog[0].filename);
+    return (FALSE);
+  case 1:
+    if (!dvo_catalog_load (catalog, VERBOSE)) {
+      if (VERBOSE) fprintf (stderr, "failure loading catalog\n");
+      return (FALSE);
+    }
+    if (!dvo_catalog_check (catalog, Nsecfilt, TRUE)) {
+      if (VERBOSE) fprintf (stderr, "can't reduce number of secondary filters\n");
+      return (FALSE);
+    }
+    if (VERBOSE) fprintf (stderr, "loaded existing file %s\n", catalog[0].filename);
+    break;
+  case 2:
+    if (MODE == DVO_READ_ONLY) return (FALSE);
+    dvo_catalog_create (region, catalog, NSECFILT); /* fills in new header info */
+    if (VERBOSE) fprintf (stderr, "creating new file %s\n", catalog[0].filename);
+    break;
+  }
+  return (TRUE);
+}
+
+int dvo_catalog_load (Catalog *catalog, int VERBOSE) {
   
   int Naxis, split;
@@ -53 +136 @@
      an old table will keep its mode
   */
+
+  /* check if the table has been sorted or not */
+  status = gfits_scan (&catalog[0].header, "SORTED", "%t", 1, &catalog[0].sorted);
+  if (!status) catalog[0].sorted = TRUE;
 
   catalog[0].catmode = DVO_MODE_MEF;
@@ -68 +155 @@
     case DVO_MODE_RAW:
       if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_RAW)\n");
-      load_catalog_raw (catalog, VERBOSE);
+      dvo_catalog_load_raw (catalog, VERBOSE);
       break;
     case DVO_MODE_MEF:
       if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_MEF)\n");
-      load_catalog_mef (catalog, VERBOSE);
+      dvo_catalog_load_mef (catalog, VERBOSE);
       break;
     case DVO_MODE_SPLIT:
       if (VERBOSE) fprintf (stderr, "reading catalog (mode DVO_MODE_SPLIT)\n");
-      load_catalog_split (catalog, VERBOSE);
+      dvo_catalog_load_split (catalog, VERBOSE);
       break;
     default:
@@ -85 +172 @@
 }
 
-int save_catalog (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_save (Catalog *catalog, char VERBOSE) {
 
   switch (catalog[0].catmode) {
     case DVO_MODE_RAW:
-      save_catalog_raw (catalog, VERBOSE);
+      dvo_catalog_save_raw (catalog, VERBOSE);
       break;
     case DVO_MODE_MEF:
-      save_catalog_mef (catalog, VERBOSE);
+      dvo_catalog_save_mef (catalog, VERBOSE);
       break;
     case DVO_MODE_SPLIT:
-      save_catalog_split (catalog, VERBOSE);
+      dvo_catalog_save_split (catalog, VERBOSE);
       break;
     default:
@@ -104 +191 @@
 }
 
-int update_catalog (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_update (Catalog *catalog, char VERBOSE) {
 
   /* update is only valid for catmode SPLIT */
   switch (catalog[0].catmode) {
     case DVO_MODE_RAW:
-      save_catalog_raw (catalog, VERBOSE);
+      dvo_catalog_save_raw (catalog, VERBOSE);
       break;
     case DVO_MODE_MEF:
-      save_catalog_mef (catalog, VERBOSE);
+      dvo_catalog_save_mef (catalog, VERBOSE);
       break;
     case DVO_MODE_SPLIT:
       /* new file needs to use save_catalog_split */
       if (catalog[0].Nave_disk == 0) {
-        save_catalog_split (catalog, VERBOSE);
+        dvo_catalog_save_split (catalog, VERBOSE);
       } else {
-        update_catalog_split (catalog, VERBOSE);
+        dvo_catalog_update_split (catalog, VERBOSE);
       }
       break;
@@ -125 +212 @@
       fprintf (stderr, "invalid catalog mode\n");
       exit (2);
+  }
+  return (TRUE);
+}
+
+int dvo_catalog_check (Catalog *catalog, int Nsecfilt, int extend) {
+
+  int i, j, Nextra, in, out;
+  SecFilt *insec, *outsec;
+
+  if (catalog[0].Nsecfilt == Nsecfilt) return (TRUE);
+
+  if (catalog[0].Nsecfilt > Nsecfilt) return (FALSE);
+
+  if ((catalog[0].Nsecfilt < Nsecfilt) && !extend) return (FALSE);
+
+  if ((catalog[0].Nsecfilt < Nsecfilt) && extend) {
+
+    Nextra = Nsecfilt - catalog[0].Nsecfilt;
+    insec = catalog[0].secfilt;
+    ALLOCATE (outsec, SecFilt, catalog[0].Naverage * Nsecfilt);
+    for (in = out = i = 0; i < catalog[0].Naverage; i++) {
+      for (j = 0; j < catalog[0].Nsecfilt; j++, in++, out++) {
+        outsec[out].M_PS  = insec[in].M_PS;
+        outsec[out].dM_PS = insec[in].dM_PS;
+        outsec[out].Xm = insec[in].Xm;
+      }
+      for (j = 0; j < Nextra; j++, out++) {
+        outsec[out].M_PS  = NO_MAG;
+        outsec[out].dM_PS = NO_MAG;
+        outsec[out].Xm    = NO_MAG;
+      }
+    }
+    free (catalog[0].secfilt);
+    catalog[0].secfilt = outsec;
+    catalog[0].Nsecfilt = Nsecfilt;
   }
   return (TRUE);

trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c

@@ -1 +1 @@
 # include <dvo.h>
 
-int load_catalog_mef (Catalog *catalog, int VERBOSE) {
+int dvo_catalog_load_mef (Catalog *catalog, int VERBOSE) {
   
   int Nitems, Nbytes, Nexpect, Naverage, Nmeasure, Nmissing, Nsecfilt;
@@ -131 +131 @@
 /* save_catalog_mef writes a complete new file from scratch */
 
-int save_catalog_mef (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_save_mef (Catalog *catalog, char VERBOSE) {
 
   int Nitems;

trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c

@@ -3 +3 @@
 /* read data from raw-style catalog file; set data format values based on file */
 
-int load_catalog_raw (Catalog *catalog, int VERBOSE) {
+int dvo_catalog_load_raw (Catalog *catalog, int VERBOSE) {
   
   int Nitems, nitems;
@@ -189 +189 @@
 }
 
-int save_catalog_raw (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_save_raw (Catalog *catalog, char VERBOSE) {
 
   int Nitems, nitems;

trunk/Ohana/src/libdvo/src/dvo_catalog_split.c

@@ -1 +1 @@
 # include <dvo.h>
 
-int load_catalog_split (Catalog *catalog, int VERBOSE) {
+int dvo_catalog_load_split (Catalog *catalog, int VERBOSE) {
 
   int Nitems, Nexpect, Naverage, Nmeasure, Nmissing, Nsecfilt;
@@ -68 +68 @@
   if (catalog[0].catflags & LOAD_MEAS) {
     ALLOCATE (measure, Catalog, 1);
+    measure[0].measure_catalog = NULL;
+    measure[0].missing_catalog = NULL;
+    measure[0].secfilt_catalog = NULL;
 
     /* get split filename from main header (paths relative to cpt file) */
@@ -145 +148 @@
   if (catalog[0].catflags & LOAD_MISS) {
     ALLOCATE (missing, Catalog, 1);
+    missing[0].measure_catalog = NULL;
+    missing[0].missing_catalog = NULL;
+    missing[0].secfilt_catalog = NULL;
 
     /* get split filename from main header (paths relative to cpt file) */
@@ -192 +198 @@
   if (catalog[0].catflags & LOAD_SECF) {
     ALLOCATE (secfilt, Catalog, 1);
-  
+    secfilt[0].measure_catalog = NULL;
+    secfilt[0].missing_catalog = NULL;
+    secfilt[0].secfilt_catalog = NULL;
+
     /* get split filename from main header (paths relative to cpt file) */
     if (!gfits_scan (&catalog[0].header, "SECFILT",  "%s", 1, string)) return (FALSE);
@@ -244 +253 @@
 /* save_catalog_split writes complete new files from scratch */
 
-int save_catalog_split (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_save_split (Catalog *catalog, char VERBOSE) {
 
   int Nitems;
@@ -434 +443 @@
  */
 
-int update_catalog_split (Catalog *catalog, char VERBOSE) {
+int dvo_catalog_update_split (Catalog *catalog, char VERBOSE) {
 
   int i, Nx, Ny, Nlines;

trunk/Ohana/src/libdvo/src/skyregion_gsc.c

@@ -1 +1 @@
 # include "dvo.h"
-# define NBANDS 24
+# define NDECBANDS 24
 # define NDIV 4
 
-static int DecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 25, 0,
-                         597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48, 0};
-
-static 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};
-
-static 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"};
+static int DecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 25, 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48};
+static double DecMin[] = {0.0, +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, -52.5, -60.0, -67.5, -75.0, -82.5, -90.0};
+static double DecMax[] = {+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};
+static char *DecNames[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", "n5230", "n6000", "n6730", "n7500", "n8230", "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", "s5230", "s6000", "s6730", "s7500", "s8230"};
+
+// L0 : full sky
+// L1 : Dec bands
+// L2 : RA segments
+// L3 : GSC regions
+// L4 : GSC subdivisions
+
+// a zone contains a set of L3 regions of the same size
+typedef struct {
+  int L1band;   // parent band
+  int L3start;  // start of zone in L1 band
+  int L3end;    // end of zone in L1 band (last + 1)
+  float dDec;   // height of L3 region in zone
+  float Rmin;
+  float Rmax;
+  float Dmin;
+  float Dmax;
+  int Nset;
+} SkyRegionZone;
+
+SkyTable *SkyRegionForDecBand (char *buffer, int Nregions, char *DecName, float Dmin, float Dmax);
+SkyRegionZone *SkyRegionFindZones (SkyTable *band, int *Nzones, int parent);
+void SkyTableL2fromZone (SkyTable *L2, SkyTable *L3, SkyTable *L4, SkyTable *band, SkyRegionZone *zone, int parent);
+void SkyTableL3fromL2 (SkyRegion *L2, SkyTable *L3, SkyTable *L4, SkyTable *band, int Ns, int Ne);
+void SkyTableL4fromL3 (SkyRegion *L3, SkyTable *L4);
 
 void SkyTableSort (SkyTable *table);
+void SkyTableAppend (SkyTable *old, SkyTable *new, int Nprev);
+int NsetForDecRange (float dDec);
+void SkyRegionPrint (SkyRegion *region);
 
 SkyTable *SkyTableFromGSC (char *filename, int depth, int VERBOSE) {
 
-  int i, j, Nx, Ny, No, Nb, Nr, NR, Nlines;
-  int nx, ny, Ne, Ns, Nbox;
-  double RA0, RA1, DEC0, DEC1, dR, dD;
+  int i, j, skipLines, Nzones;
   char temp[80], name[80];
 
@@ -25 +47 @@
   SkyTable *skytable;
   SkyRegion *regions;
-  SkyTable bandregions;
+  SkyTable *band;
+  SkyTable L0, L1, L2, L3, L4;
+  SkyRegionZone *zones;
+
   FILE *f;
   
@@ -34 +59 @@
   }
 
-  /* load in table data */
+/* load in table data */
   ftable.header = &theader;
   if (!gfits_fread_ftable (f, &ftable, "REGIONS")) {
@@ -41 +66 @@
     return (NULL);
   }
-
-  gfits_scan (ftable.header, "NAXIS1", "%d", 1, &Nx);
-  gfits_scan (ftable.header, "NAXIS2", "%d", 1, &Ny);
-  
   gfits_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 : identify the blocks in the DEC bands */
-  No = 1;
-  for (i = 0; i < 25; i++) {
-    if (i == 12) continue;
-
-    Nlines = 0;
-    for (j = 0; j < i; j++) Nlines += DecLines[j];
-
-    /* find and save all GSC Regions in this band */
-    Nb = 0;
-    ALLOCATE (bandregions.regions, SkyRegion, DecLines[i]);
-    ALLOCATE (bandregions.filename, char *, DecLines[i]);
-    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;
-
-      bandregions.regions[Nb].Dmin = DEC0;
-      bandregions.regions[Nb].Dmax = DEC1;
-      bandregions.regions[Nb].Rmin = RA0;
-      bandregions.regions[Nb].Rmax = RA1;
-      bandregions.filename[Nb] = NULL;
-      Nb ++;
-    }
-    bandregions.Nregions = Nb;
+/* generate the regions for each level independently. indices (parent,childE,childS) 
+   refer to the value within the independent group.  at the end, we combine and adjust
+   the indices */
+
+/* L0 : full sky */
+  L0.Nregions = 1;
+  ALLOCATE (L0.regions, SkyRegion, L0.Nregions);
+  L0.regions[0].Rmin    =   0;
+  L0.regions[0].Rmax    = 360;
+  L0.regions[0].Dmin    = -90;
+  L0.regions[0].Dmax    = +90;
+  L0.regions[0].index   =  0;
+  L0.regions[0].depth   =  0;
+  L0.regions[0].parent  = -1;
+  L0.regions[0].child   =  TRUE;
+  L0.regions[0].table   = -1;
+  L0.regions[0].childS  =  0;
+  L0.regions[0].childE  =  NDECBANDS;
+  strcpy (L0.regions[0].name, "fullsky");
+  // SkyRegionPrint (&L0.regions[0]);
+
+  /* allocate space for all levels */
+  L1.Nregions = L2.Nregions = L3.Nregions = L4.Nregions = 0;
+  ALLOCATE (L1.regions, SkyRegion, 1);
+  ALLOCATE (L2.regions, SkyRegion, 1);
+  ALLOCATE (L3.regions, SkyRegion, 1);
+  ALLOCATE (L4.regions, SkyRegion, 1);
+
+  /* L1 : dec bands */
+  skipLines = 0;
+  for (i = 0; i < NDECBANDS; i++) {
+    L1.regions[i].Rmin     = 0;
+    L1.regions[i].Rmax     = 360;
+    L1.regions[i].Dmin     = DecMin[i];
+    L1.regions[i].Dmax     = DecMax[i];
+    L1.regions[i].index    = i;
+    L1.regions[i].depth    = 1;
+    L1.regions[i].parent   = 0;
+    L1.regions[i].child    = TRUE;
+    L1.regions[i].table    = -1;
+    strcpy (L1.regions[i].name, DecNames[i]);
+    // SkyRegionPrint (&L1.regions[i]);
+
+    /* build the L2 regions for this L1 region (based on zones) */
+    L1.regions[i].childS   = L2.Nregions;
+
+    /* load all GSC Regions in this band */
+    band = SkyRegionForDecBand (&ftable.buffer[skipLines*48], DecLines[i], DecNames[i], L1.regions[i].Dmin, L1.regions[i].Dmax);
+    skipLines += DecLines[i];
 
     /* sort the band regions by Rmin */
-    SkyTableSort (&bandregions);
+    SkyTableSort (band);
     
-    { 
-      int k, Nt, Nrange, Ndiv, Nset;
-      int *transitions;
-      float dDec, dDcur;
-
-      /* search for transitions in the region dDec */
-      ALLOCATE (transitions, int, 50);
-      Nt = 0;
-      dDec = bandregions.regions[0].Dmax - bandregions.regions[0].Dmin;
-      dDcur = dDec;
-      transitions[0] = 0;
-      Nt++;
-      DEC0 = bandregions.regions[0].Dmin;
-      DEC1 = bandregions.regions[0].Dmax;
-      for (j = 0; j < Nb; j++) {
-        DEC0 = MIN (bandregions.regions[j].Dmin, DEC0);
-        DEC1 = MAX (bandregions.regions[j].Dmax, DEC1);
-        dDec = bandregions.regions[j].Dmax - bandregions.regions[j].Dmin;
-        if (dDec != dDcur) {
-          dDcur = dDec;
-          transitions[Nt] = j;
-          Nt++;
-        }
-      }
-      transitions[Nt] = Nb;
-      Nt++;
-
-      /* subdivide each transition range */
-      for (j = 1; j < Nt; j++) {
-        Ne = transitions[j];
-        Ns = transitions[j-1];
-        dDec = bandregions.regions[Ns].Dmax - bandregions.regions[Ns].Dmin;
-        fprintf (stderr, "trans: %d - %d : ", Ns, Ne);
-        fprintf (stderr, "%f,%f - %f,%f\n", 
-                 bandregions.regions[Ns].Rmin, bandregions.regions[Ne-1].Rmax,
-                 bandregions.regions[Ns].Dmin, bandregions.regions[Ns].Dmax);
-
-        Nrange = transitions[j] - transitions[j-1];
-        Ndiv = (int) (0.5 + 7.5 / dDec);
-        switch (Ndiv) {
-          case 2:
-            Nset = 4*2;
-            break;
-          case 3:
-            Nset = 6*3;
-            break;
-          case 4:
-            Nset = 8*4;
-            break;
-          default:
-            fprintf (stderr, "programming error\n");
-            exit (2);
-        }
-        for (k = Ns; k < Ne; k += Nset) {
-          regions[Nr].Rmin = bandregions.regions[k].Rmin;
-          if (k + Nset < Ne) {
-            regions[Nr].Rmax = bandregions.regions[k+Nset-1].Rmax;
-          } else {
-            regions[Nr].Rmax = bandregions.regions[Ne-1].Rmax;
-          }     
-          regions[Nr].Dmin = DEC0;
-          regions[Nr].Dmax = DEC1;
-          fprintf (stderr, "new: %f - %f :%f - %f\n", 
-                   regions[Nr].Rmin, regions[Nr].Rmax, 
-                   regions[Nr].Dmin, regions[Nr].Dmax);
-        }
-
-        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;
-      }
-      regions[No].childS = Nr;
-      regions[No].childE = Nr;
-      No ++;
-    }
-  }
-
-  /* level 3 : copy the data from the GSC Region files */
-  /* XXX fix this: level here is now 3 */
-  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 */
-  /* XXX level here is now 4 */
-  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;
-  }
-
+    zones = SkyRegionFindZones (band, &Nzones, i);
+
+    /* subdivide each zone */
+    for (j = 0; j < Nzones; j++) {
+      SkyTableL2fromZone (&L2, &L3, &L4, band, &zones[j], i);
+    }
+    free (zones);
+    SkyTableFree (band);
+
+    /* at end of loop, L2.Nregions has been updated to end of L2.regions for L1.region */
+    L1.regions[i].childE   = L2.Nregions;
+
+    L1.Nregions ++;
+    REALLOCATE (L1.regions, SkyRegion, L1.Nregions + 1);
+  }
+
+  /* merge L1 into L0 */
   ALLOCATE (skytable, SkyTable, 1);
-  skytable[0].regions = regions;
-  skytable[0].Nregions = Nr;
-
+  ALLOCATE (skytable[0].regions, SkyRegion, 1);
+  skytable[0].Nregions = 0;
+
+  SkyTableAppend (skytable, &L0, 0);
+  SkyTableAppend (skytable, &L1, skytable[0].Nregions - L0.Nregions);
+  SkyTableAppend (skytable, &L2, skytable[0].Nregions - L1.Nregions);
+  SkyTableAppend (skytable, &L3, skytable[0].Nregions - L2.Nregions);
+  SkyTableAppend (skytable, &L4, skytable[0].Nregions - L3.Nregions);
+
+  /* fix the depth elements and create the filename array */
   ALLOCATE (skytable[0].filename, char *, skytable[0].Nregions);
   for (i = 0; i < skytable[0].Nregions; i++) {
     skytable[0].filename[i] = NULL;
-  }
-
+    skytable[0].regions[i].table = (skytable[0].regions[i].depth == depth);
+  }
   return (skytable);
+}
+
+SkyTable *SkyRegionForDecBand (char *buffer, int Nregions, char *DecName, float DminBand, float DmaxBand) {
+
+  int i;
+  double Rmin, Rmax, Dmin, Dmax;
+  char temp[80], name[80];
+  SkyTable *band;
+  SkyRegion *regions;
+
+  ALLOCATE (band, SkyTable, 1);
+  ALLOCATE (regions, SkyRegion, Nregions);
+  for (i = 0; i < Nregions; i++) {
+    strncpy (temp, &buffer[i*48], 48);
+    temp[49] = 0;
+    hstgsc_hms_to_deg (&Rmin, &Rmax, &Dmin, &Dmax, &temp[7]);
+    if (Dmax < Dmin) SWAP (Dmin, Dmax);
+
+    /* check that we are in correct Dec band */
+    if ((Dmax < DminBand) || (Dmin > DmaxBand)) {
+      fprintf (stderr, "error with raw table: table mis-match\n");
+      fprintf (stderr, "line: %s\n", temp);
+      fprintf (stderr, "region %d: %f - %f vs %f - %f\n", i, Dmin, Dmax, DminBand, DmaxBand);
+      exit (2);
+    }
+
+    /* regions near 0,360 boundary have Rmax == 0.0 */
+    if (Rmax < Rmin) Rmax += 360.0;
+
+    regions[i].Dmin = Dmin;
+    regions[i].Dmax = Dmax;
+    regions[i].Rmin = Rmin;
+    regions[i].Rmax = Rmax;
+
+    temp[5] = 0;
+    sprintf (name, "%s/%s", DecName, &temp[1]);
+    strcpy (regions[i].name, name);
+  }
+  band[0].filename = NULL;
+  band[0].regions = regions;
+  band[0].Nregions = Nregions;
+  return (band);
 }
 
@@ -355 +219 @@
       l--;
       tempregion = regions[l];
-      tempfile   = filename[l];
+      if (filename) tempfile   = filename[l];
     } else {
       tempregion   = regions[ir];
-      tempfile     = filename[ir];
+      if (filename) tempfile     = filename[ir];
       regions[ir]   = regions[0];
-      filename[ir] = filename[0];
+      if (filename) filename[ir] = filename[0];
       if (--ir == 0) {
         regions[0]   = tempregion;
-        filename[0] = tempfile;
+        if (filename) filename[0] = tempfile;
         return;
       }
@@ -373 +237 @@
       if (tempregion.Rmin < regions[j].Rmin) {
         regions[i]   = regions[j];
-        filename[i] = filename[j];
+        if (filename) filename[i] = filename[j];
         j += (i=j) + 1;
       } 
@@ -379 +243 @@
     }
     regions[i]   = tempregion;
-    filename[i] = tempfile;
-  }
-}
+    if (filename) filename[i] = tempfile;
+  }
+}
+
+/* give a complete set of regions in a Dec band, subdivide into zones of equal-sized regions */
+SkyRegionZone *SkyRegionFindZones (SkyTable *band, int *Nzones, int parent) {
+  
+  float Dmin, Dmax, dDec;
+  int i, Nz, NZ, Ndiv, Nregions;
+  SkyRegion *regions;
+  SkyRegionZone *zones;
+  
+  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);
+
+  /* search for transitions in the region dDec, find the max Dec range */
+  Dmin = regions[0].Dmin;
+  Dmax = regions[0].Dmax;
+  for (i = 0; i < Nregions; i++) {
+    Dmin = MIN (regions[i].Dmin, Dmin);
+    Dmax = MAX (regions[i].Dmax, Dmax);
+    dDec = regions[i].Dmax - regions[i].Dmin;
+    if (dDec != zones[Nz].dDec) {
+      /* we've found the end of the current zone */
+      zones[Nz].L3end = i;
+      zones[Nz].Rmax = regions[i-1].Rmax;
+      Nz++;
+      CHECK_REALLOCATE (zones, SkyRegionZone, NZ, Nz, 10);
+
+      /* start info for the new zone */
+      zones[Nz].Rmin = regions[i].Rmin;
+      zones[Nz].dDec = dDec;
+      zones[Nz].L3start = i;
+      zones[Nz].L1band = parent;
+      zones[Nz].Nset = NsetForDecRange (zones[Nz].dDec);
+    }
+  }
+  zones[Nz].L3end = i;
+  zones[Nz].Rmax = regions[i-1].Rmax;
+  Nz ++;
+
+  for (i = 0; i < Nz; i++) {
+    zones[i].Dmin = Dmin;
+    zones[i].Dmax = Dmax;
+  }
+
+  *Nzones = Nz;
+  return (zones);
+}
+
+int NsetForDecRange (float dDec) {
+
+  int Ndiv, Nset;
+  
+  /* max segment size is ~7.5 x 15 degrees */
+  Ndiv = (int) (0.5 + 7.5 / dDec);
+  Nset = 2*Ndiv*Ndiv;
+  return (Nset);
+}
+  
+void SkyTableL2fromZone (SkyTable *L2, SkyTable *L3, SkyTable *L4, SkyTable *band, SkyRegionZone *zone, int parent) {
+
+  int i, Nr, Ns, Ne;
+  char *p, name[80];
+
+  Nr = L2[0].Nregions;
+  REALLOCATE (L2[0].regions, SkyRegion, Nr + 1);
+  
+  /* divide this zone into L2 regions with Nset L3 regions each (fewer on ends) */
+  for (i = zone[0].L3start; i < zone[0].L3end; i += zone[0].Nset) {
+    Ns = i; 
+    Ne = MIN (i + zone[0].Nset, zone[0].L3end);
+
+    L2[0].regions[Nr].Rmin = band[0].regions[Ns].Rmin;
+    L2[0].regions[Nr].Rmax = band[0].regions[Ne - 1].Rmax;
+    L2[0].regions[Nr].Dmin = zone[0].Dmin;
+    L2[0].regions[Nr].Dmax = zone[0].Dmax;
+    
+    L2[0].regions[Nr].index    =  Nr;
+    L2[0].regions[Nr].depth    =  2;
+    L2[0].regions[Nr].table    =  -1;
+    L2[0].regions[Nr].parent   =  parent;
+    L2[0].regions[Nr].child    =  FALSE;
+    L2[0].regions[Nr].childS   =  0;
+    L2[0].regions[Nr].childE   =  0;
+
+    /* define names for L2 regions */
+    p = strchr (band[0].regions[Ns].name, '/');
+    if (p == NULL) {
+      fprintf (stderr, "programming error in SkyTableL2fromZone\n");
+      exit (2);
+    }
+    *p = 0;
+    sprintf (name, "%s/z%03d", band[0].regions[Ns].name, Nr);
+    *p = '/';
+    strcpy (L2[0].regions[Nr].name, name);
+    // SkyRegionPrint (&L2[0].regions[Nr]);
+
+    /* childS and childE are set in SkyTableL3fromL2 */
+    SkyTableL3fromL2 (&L2[0].regions[Nr], L3, L4, band, Ns, Ne);
+
+    Nr++;
+    REALLOCATE (L2[0].regions, SkyRegion, Nr + 1);
+  }
+  L2[0].Nregions = Nr;
+}
+
+void SkyTableL3fromL2 (SkyRegion *L2, SkyTable *L3, SkyTable *L4, SkyTable *band, int Ns, int Ne) {
+
+  int i, Nr;
+
+  Nr = L3[0].Nregions;
+  L3[0].Nregions += Ne - Ns;
+  REALLOCATE (L3[0].regions, SkyRegion, L3[0].Nregions);
+
+  L2[0].child  = TRUE;
+  L2[0].childS = Nr;
+  L2[0].childE = L3[0].Nregions;
+
+  /* copy the band entries corresponding to this region in the L3 table */
+  for (i = Ns; i < Ne; i++) {
+
+    L3[0].regions[Nr] = band[0].regions[i];
+    
+    L3[0].regions[Nr].index    =  Nr;
+    L3[0].regions[Nr].depth    =  3;
+    L3[0].regions[Nr].table    =  -1;
+    L3[0].regions[Nr].parent   =  L2[0].index;
+    L3[0].regions[Nr].child    =  FALSE;
+    L3[0].regions[Nr].childS   =  0;
+    L3[0].regions[Nr].childE   =  0;
+
+    // SkyRegionPrint (&L3[0].regions[Nr]);
+    /* name is set for the band in SkyRegionForDecBand */
+
+    /* childS and childE are set in SkyTableL3fromL3 */
+    SkyTableL4fromL3 (&L3[0].regions[Nr], L4);
+    Nr++;
+  }
+
+  return;
+}
+
+void SkyTableL4fromL3 (SkyRegion *L3, SkyTable *L4) {
+
+  int nx, ny, Nr, Nbox;
+  double Rmin, Dmin, dR, dD;
+  char name[80];
+
+  Nr = L4[0].Nregions;
+  L4[0].Nregions += NDIV*NDIV;
+  REALLOCATE (L4[0].regions, SkyRegion, L4[0].Nregions);
+
+  L3[0].child  = TRUE;
+  L3[0].childS = Nr;
+  L3[0].childE = L4[0].Nregions;
+
+  /* subdivide L3 into NDIV boxes */
+  Rmin = L3[0].Rmin;
+  Dmin = L3[0].Dmin;
+  dR = (L3[0].Rmax - L3[0].Rmin) / NDIV;
+  dD = (L3[0].Dmax - L3[0].Dmin) / NDIV;
+
+  Nbox = 0;
+  for (ny = 0; ny < NDIV; ny ++) {
+    for (nx = 0; nx < NDIV; nx ++) {
+      L4[0].regions[Nr].Rmin     = Rmin  + (nx + 0)*dR;
+      L4[0].regions[Nr].Rmax     = Rmin  + (nx + 1)*dR;
+      L4[0].regions[Nr].Dmin     = Dmin + (ny + 0)*dD;
+      L4[0].regions[Nr].Dmax     = Dmin + (ny + 1)*dD;
+
+      L4[0].regions[Nr].index    =  Nr;
+      L4[0].regions[Nr].depth    =  4;
+      L4[0].regions[Nr].table    =  -1;
+      L4[0].regions[Nr].parent   =  L3[0].index;
+      L4[0].regions[Nr].child    =  FALSE;
+      L4[0].regions[Nr].childS   =  0;
+      L4[0].regions[Nr].childE   =  0;
+
+      sprintf (name, "%s.%02d", L3[0].name, Nbox);
+      strcpy (L4[0].regions[Nr].name, name);
+      // SkyRegionPrint (&L4[0].regions[Nr]);
+
+      Nr ++;
+      Nbox ++;
+    }
+  }
+  return;
+}
+
+void SkyTableAppend (SkyTable *old, SkyTable *new, int Nprev) {
+
+  int i, Nold;
+
+  Nold = old[0].Nregions;
+
+  old[0].Nregions += new[0].Nregions;
+  REALLOCATE (old[0].regions, SkyRegion, old[0].Nregions);
+
+  for (i = Nprev; i < Nold; i++) {
+    old[0].regions[i].childS += Nold;
+    old[0].regions[i].childE += Nold;
+  }
+
+  for (i = 0; i < new[0].Nregions; i++) {
+    old[0].regions[i + Nold] = new[0].regions[i];
+    old[0].regions[i + Nold].parent += Nprev;
+  }
+  return;
+}
+
+void SkyRegionPrint (SkyRegion *region) {
+
+  int i;
+
+  fprintf (stderr, "L%d:", region[0].depth);
+  for (i = 0; i < region[0].depth; i++) {
+    fprintf (stderr, " ");
+  } 
+  fprintf (stderr, "%s : %f - %f : %f - %f\n", region[0].name, region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax);
+  return;
+}
+
+/***
+    notes and questions:
+    1) is the regions.index value used?
+
+
+***/

trunk/Ohana/src/libdvo/src/skyregion_ops.c

@@ -327 +327 @@
   int i;
 
+  /* XXX do we need to free the filename array as well? */
   if (list == NULL) return (TRUE);
   if (list[0].regions != NULL) {
@@ -349 +350 @@
 
   if (table == NULL) return (TRUE);
-  if (table[0].regions != NULL) {
+  if (table[0].filename != NULL) {
     for (i = 0; i < table[0].Nregions; i++) {
       if (table[0].filename[i] != NULL) {
@@ -355 +356 @@
       }
     }
+    free (table[0].filename);
+  }
+  if (table[0].regions != NULL) {
     free (table[0].regions);
   }

trunk/Ohana/src/opihi/dvo/skycat.c

@@ -6 +6 @@
   
   double Radius;
-  int i, j, N, Nregions, ShowAll, NPTS, Npts, leftside, Depth, TableDepth, VERBOSE;
+  int i, j, N, Nregions, ShowAll, NPTS, Npts, leftside, Depth, VERBOSE;
   struct stat filestat;
   Vector Xvec, Yvec;
@@ -39 +39 @@
   }
   SetGraph (graphmode);
-
-  TableDepth = (Depth == 3) ? 3 : 2;
 
   Radius = MAX (fabs(graphmode.xmax), fabs(graphmode.ymax));

trunk/Ohana/src/relphot/src/load_catalogs.c

@@ -15 +15 @@
   for (i = 0; i < skylist[0].Nregions; i++) {
     tcatalog.filename = skylist[0].filename[i];
+
+    tcatalog.catformat = dvo_catalog_format (CATFORMAT);  // set the default catformat from config data
+    tcatalog.catmode   = dvo_catalog_mode (CATMODE);      // set the default catmode from config data
+    tcatalog.lockmode  = LCK_SOFT;
+    dvo_catalog_open (&tcatalog, skylist[0].regions[i], Nsecfilt, "r");
+
     switch (lock_catalog (&tcatalog, LCK_SOFT)) {
     case 0: