Changeset 8295

Timestamp:

Aug 11, 2006, 4:50:07 PM (20 years ago)

Author:

eugene

Message:

working on intermediate level between dec bands and gsc regions

Location:

trunk/Ohana/src/libdvo

Files:

• Makefile (modified) (1 diff)
• src/skyregion_gsc.c (modified) (5 diffs)
• src/tabletest.c (added)

trunk/Ohana/src/libdvo/Makefile

@@ -105 +105 @@
         cp $(INC)/$* $(DESTINC)/
 
+tabletest: install
+        gcc -L$(LLIB) -I$(LINC) -o tabletest $(SRC)/tabletest.c -ldvo -lFITS -lohana -lm
+
 clean:
         rm -f */*~

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

@@ -12 +12 @@
                            "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", "s5230", "s6000", "s6730", "s7500", "s8230", "none"};
 
+void SkyTableSort (SkyTable *table);
+
 SkyTable *SkyTableFromGSC (char *filename, int depth, int VERBOSE) {
 
-  int i, j, Nx, Ny, No, Nr, NR, Nlines;
+  int i, j, Nx, Ny, No, Nb, Nr, NR, Nlines;
   int nx, ny, Ne, Ns, Nbox;
   double RA0, RA1, DEC0, DEC1, dR, dD;
@@ -23 +25 @@
   SkyTable *skytable;
   SkyRegion *regions;
+  SkyTable bandregions;
   FILE *f;
   
@@ -95 +98 @@
   regions[No].childE = Nr;
 
-  /* level 2 : copy the data from the GSC Region files */
+  /* 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;
+
+    /* sort the band regions by Rmin */
+    SkyTableSort (&bandregions);
+    
+    { 
+      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++) {
@@ -152 +284 @@
 
   /* 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;
@@ -201 +334 @@
   return (skytable);
 }
+
+/* sort skytable by Rmin */
+void SkyTableSort (SkyTable *table) {
+
+  int i, j, k, l, ir, N;
+  SkyRegion *regions;
+  SkyRegion tempregion;
+  char **filename;
+  char *tempfile;
+  
+  N = table[0].Nregions;
+  regions = table[0].regions;
+  filename = table[0].filename;
+
+  if (N < 2) return;
+  l = N >> 1;
+  ir = N - 1;
+  for (;;) {
+    if (l > 0) {
+      l--;
+      tempregion = regions[l];
+      tempfile   = filename[l];
+    } else {
+      tempregion   = regions[ir];
+      tempfile     = filename[ir];
+      regions[ir]   = regions[0];
+      filename[ir] = filename[0];
+      if (--ir == 0) {
+        regions[0]   = tempregion;
+        filename[0] = tempfile;
+        return;
+      }
+    }
+    i = l;
+    j = (l << 1) + 1;
+    while (j <= ir) {
+      if (j < ir && regions[j].Rmin < regions[j+1].Rmin) ++j;
+      if (tempregion.Rmin < regions[j].Rmin) {
+        regions[i]   = regions[j];
+        filename[i] = filename[j];
+        j += (i=j) + 1;
+      } 
+      else j = ir + 1;
+    }
+    regions[i]   = tempregion;
+    filename[i] = tempfile;
+  }
+}