Changeset 27838 for branches/tap_branches/Ohana/src/addstar/src/find_matches.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/addstar/src/find_matches.c (modified) (14 diffs)

branches/tap_branches/Ohana/src/addstar/src/find_matches.c

@@ -1 +1 @@
 # include "addstar.h"
 
-int find_matches (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, AddstarClientOptions options) {
-
-  int i, j, n, N, J, status, Nstars;
+int find_matches (SkyRegion *region, Stars *stars, unsigned int NstarsIn, Catalog *catalog, AddstarClientOptions options) {
+
+  off_t i, j, n, N, J, status, Nstars;
   double RADIUS, RADIUS2;
   double *X1, *Y1, *X2, *Y2;
   double dX, dY, dR;
-  int *N1, *N2,  *next_meas;
-  int Nave, NAVE, Nmeas, NMEAS, Nmatch;
+  off_t *N1, *N2, *next_meas;
+  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
   int Nsecfilt, Nsec;
   unsigned int objID, catID;
   Coords tcoords;
+
+  if (NSTAR_GROUP <= 0) {
+      fprintf (stderr, "ERROR: NSTAR_GROUP NOT SET!\n");
+      exit (1);
+  }
 
   /* photcode data - must by of type DEP; options.photcode is equiv photcode for all input
@@ -24 +29 @@
   ALLOCATE (X1, double, NstarsIn);
   ALLOCATE (Y1, double, NstarsIn);
-  ALLOCATE (N1, int,   NstarsIn);
+  ALLOCATE (N1, off_t,  NstarsIn);
 
   /** allocate local arrays (catalog) **/
@@ -30 +35 @@
   ALLOCATE (X2, double, NAVE);
   ALLOCATE (Y2, double, NAVE);
-  ALLOCATE (N2, int,   NAVE);
-  ALLOCATE (catalog[0].found, int, NAVE);
+  ALLOCATE (N2, off_t,  NAVE);
+  ALLOCATE (catalog[0].found, off_t, NAVE);
   /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
 
@@ -103 +108 @@
 
   /** find matched stars **/
+  // XXX could use NSTAR_GROUP to do this match more quicky if NSTAR_GROUP > 1
   for (i = j = 0; (i < Nstars) && (j < Nave); ) {
     if (!finite(X1[i]) || !finite(Y1[i])) { 
@@ -115 +121 @@
     /* negative dX: j is too large */
     dX = X1[i] - X2[j];
-    if (dX <= -2*RADIUS) {
+    if (dX <= -1.02*RADIUS) {
       i++;
       continue;
     }
     /* positive dX, i is too large */
-    if (dX >= 2*RADIUS) {
+    if (dX >= 1.02*RADIUS) {
       j++;
       continue;
@@ -126 +132 @@
 
     /* within match range; look for matches */
-    for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) {
+    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
       dX = X1[i] - X2[J];
       dY = Y1[i] - Y2[J];
@@ -132 +138 @@
       if (dR > RADIUS2) continue;
 
+      /*** a match is found, add to average, measure ***/
+      Nmatch ++;
+      n = N2[J];
+      N = N1[i];
+
       /* make sure there is space for next entry */
       if (Nmeas >= NMEAS) {
         NMEAS = Nmeas + 1000;
-        REALLOCATE (next_meas, int, NMEAS);
+        REALLOCATE (next_meas, off_t, NMEAS);
+
         REALLOCATE (catalog[0].measure, Measure, NMEAS);
       }
-
-      Nmatch ++;
-      n = N2[J];
-      N = N1[i];
 
       /* add to end of measurement list */
@@ -226 +234 @@
   /* incorporate unmatched image stars, if this star is in field of this catalog */
   /* these new entries are all written out in UPDATE mode */ 
-  for (i = 0; (i < Nstars) && !options.only_match; i++) {
+  for (i = 0; (i < Nstars) && !options.only_match; i += NSTAR_GROUP) {
     /* make sure there is space for next entry */
-    if (Nmeas >= NMEAS) {
+    if (Nmeas >= NMEAS - NSTAR_GROUP) {
       NMEAS = Nmeas + 1000;
-      REALLOCATE (next_meas, int, NMEAS);
+      REALLOCATE (next_meas, off_t, NMEAS);
       REALLOCATE (catalog[0].measure, Measure, NMEAS);
     }
@@ -248 +256 @@
     catalog[0].average[Nave].dD            = 0;
 
-    catalog[0].average[Nave].Nmeasure      = 1;
+    catalog[0].average[Nave].Nmeasure      = NSTAR_GROUP;
     catalog[0].average[Nave].Nmissing      = 0;
     catalog[0].average[Nave].Nextend       = 0;
@@ -264 +272 @@
 
     catalog[0].average[Nave].Xp            = 0;
-    catalog[0].average[Nave].ChiSq         = 0.0;
+    catalog[0].average[Nave].ChiSqAve      = 0.0;
+    catalog[0].average[Nave].ChiSqPM       = 0.0;
+    catalog[0].average[Nave].ChiSqPar      = 0.0;
+    catalog[0].average[Nave].Tmean         = 0;
+    catalog[0].average[Nave].Trange        = 0;
     catalog[0].average[Nave].Npos          = 0;
 
@@ -270 +282 @@
     catalog[0].average[Nave].catID         = catID;
     catalog[0].average[Nave].flags         = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
 
     objID ++;
@@ -283 +301 @@
     }
 
-    // supply the measurments from this detection
-    catalog[0].measure[Nmeas]           = stars[N].measure;
-
-    // the following measure elements cannot be set until here:
-    catalog[0].measure[Nmeas].dR        = 0.0;
-    catalog[0].measure[Nmeas].dD        = 0.0;
-    catalog[0].measure[Nmeas].dbFlags   = 0;
-    catalog[0].measure[Nmeas].averef    = Nave; // XXX EAM : must be absolute Nave if partial read
-    catalog[0].measure[Nmeas].objID     = catalog[0].average[Nave].objID;
-    catalog[0].measure[Nmeas].catID     = catalog[0].catID;
-
-    /* set the average magnitude if not already set and the photcode.equiv is not 0 */
-    /* in UPDATE mode, this value is not saved; use relphot to recalculate */
-    if (Nsec > -1) { 
-        catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
-    }
-
-    /* next[Nmeas] should always be -1 in this context (it is always the only
-       measurement for the star) */
-    stars[N].found = Nmeas;
-    next_meas[Nmeas] = -1;
-    Nmeas ++;
+    for (j = 0; j < NSTAR_GROUP; j++) {
+      // supply the measurments from this detection
+      catalog[0].measure[Nmeas]           = stars[N].measure;
+      N = N1[i + j];
+
+      // the following measure elements cannot be set until here:
+      catalog[0].measure[Nmeas].dR              = 0.0;
+      catalog[0].measure[Nmeas].dD              = 0.0;
+      catalog[0].measure[Nmeas].dbFlags         = 0;
+      catalog[0].measure[Nmeas].averef          = Nave; // XXX EAM : must be absolute Nave if partial read
+      catalog[0].measure[Nmeas].objID     = catalog[0].average[Nave].objID;
+      catalog[0].measure[Nmeas].catID     = catalog[0].catID;
+
+      /* set the average magnitude if not already set and the photcode.equiv is not 0 */
+      /* in UPDATE mode, this value is not saved; use relphot to recalculate */
+      if (Nsec > -1) { 
+          catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
+      }
+
+      /* next[Nmeas] should always be -1 in this context (it is always the only
+         measurement for the star) */
+      stars[N].found = Nmeas;
+      next_meas[Nmeas] = -1;  // initial value here update below
+      Nmeas ++;
+    }
+    for (j = 0; j < NSTAR_GROUP - 1; j++) {
+      next_meas[Nmeas - NSTAR_GROUP + j] = Nmeas - NSTAR_GROUP + j + 1;
+    }
     Nave ++;
   }
@@ -332 +356 @@
   catalog[0].Nmeasure = Nmeas;
   catalog[0].Nsecf_mem = Nave*Nsecfilt;
-  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmatch);
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %lld %lld %lld, (%lld matches)\n", (long long) Nstars, (long long) Nave, (long long) Nmeas, (long long) Nmatch);
 
   free (catalog[0].found);
@@ -338 +362 @@
   free (Y1);
   free (N1);
+  free (N2);
   free (X2);
   free (Y2);
-  free (N2);
+  free (next_meas);
+
   return (Nmatch);
 }