Changeset 27840 for branches/simtest_nebulous_branches/Ohana/src/opihi/dvo/skyregion.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/opihi/dvo/skyregion.c (modified) (1 diff)

branches/simtest_nebulous_branches/Ohana/src/opihi/dvo/skyregion.c

@@ -44 +44 @@
   return TRUE;
 }
+
+/* find region which overlaps c at given depth (-1 : populated ) */
+int SkyRegionByPoint_r (SkyTable *table, SkyList *list, int depth, double ra, double dec) {
+  
+  int i, Ns, Ne, No;
+  SkyRegion *skyregion;
+
+  list[0].Nregions = 0;
+  list[0].ownElements = FALSE; // this list is only holding a view to the elements
+
+  skyregion = table[0].regions;
+
+  Ns = 0;
+  Ne = 1;
+
+  while (1) {
+    No = -1;
+    for (i = Ns; (i < Ne) && (i < table[0].Nregions); i++) {
+      if (ra  < skyregion[i].Rmin) continue;
+      if (ra  > skyregion[i].Rmax) continue;
+      if (dec < skyregion[i].Dmin) continue;
+      if (dec > skyregion[i].Dmax) continue;
+      No = i;
+      break;
+    }
+    if (No == -1) return (FALSE);
+    if ((depth == -1) && (skyregion[No].table)) break;
+    if (depth == skyregion[No].depth) break;
+    if ((depth > skyregion[No].depth) && !skyregion[No].child) return (FALSE);
+
+    /* need to check Ns, Ne, or guarantee valid range */
+    Ns = skyregion[No].childS;
+    Ne = skyregion[No].childE;
+  }
+
+  list[0].regions[0] = &skyregion[No];
+  list[0].filename[0] = table[0].filename[No];
+  list[0].Nregions = 1;
+  return (TRUE);
+}
+
+SkyList *SelectRegionsByCoordVectors (Vector *RA, Vector *DEC) {
+  
+  int i, j, Npts, Nout, NOUT, *found;
+  double ra, dec;
+  SkyList *new, *list;
+  SkyTable *sky;
+
+  sky = GetSkyTable ();
+
+  Npts = RA->Nelements;
+  ALLOCATE (found, int, Npts);
+  memset (found, 0, Npts*sizeof(int));
+
+  ALLOCATE (new, SkyList, 1);
+  ALLOCATE (new[0].regions,  SkyRegion *, 1);
+  ALLOCATE (new[0].filename,  char *, 1);
+
+  // output list
+  Nout = 0;
+  NOUT = 100;
+  ALLOCATE (list, SkyList, 1);
+  ALLOCATE (list[0].regions,  SkyRegion *, NOUT);
+  ALLOCATE (list[0].filename,  char *, NOUT);
+  list[0].Nregions = 0;
+  list[0].ownElements = FALSE; // this list is only holding a view to the elements
+
+  for (i = 0; i < Npts; i++) {
+    if (found[i]) continue;
+
+    SkyRegionByPoint_r (sky, new, -1, RA->elements.Flt[i], DEC->elements.Flt[i]);
+    if (new->Nregions == 0) continue;
+    assert (new->Nregions == 1);
+    // append new region to output list
+    list[0].regions[Nout] = new[0].regions[0];
+    list[0].filename[Nout] = new[0].filename[0];
+    Nout ++;
+    if (Nout >= NOUT) {
+        NOUT += 100;
+        REALLOCATE (list[0].regions, SkyRegion *, NOUT);
+        REALLOCATE (list[0].filename, char *, NOUT);
+    }
+    found[i] = TRUE;
+
+    // scan over the remaining points to find any that lie within this region
+    // if we sorted the ra,dec inputs we could break out of this more quickly
+    for (j = i + 1; j < Npts; j++) {
+      ra = RA->elements.Flt[j];
+      dec = DEC->elements.Flt[j];
+      if (ra  < new[0].regions[0][0].Rmin) continue;
+      if (ra  > new[0].regions[0][0].Rmax) continue;
+      if (dec < new[0].regions[0][0].Dmin) continue;
+      if (dec > new[0].regions[0][0].Dmax) continue;
+      found[j] = TRUE;
+    }
+  }
+  list[0].Nregions = Nout;
+
+  free (found);
+  SkyListFree (new);
+
+  return (list);
+}
+