Changeset 17211

Timestamp:

Mar 28, 2008, 3:41:56 PM (18 years ago)

Author:

eugene

Message:

init uR, uD, PA; clip by photcode, photflags

File:

• trunk/Ohana/src/relastro/src/UpdateObjects.c (modified) (10 diffs)

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -40 +40 @@
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
-  int i, j, k, m, N;
+  int i, j, k, m, N, Nsecfilt, found, kp;
   StatType statsR, statsD;
   Coords coords;
@@ -49 +49 @@
   float errorScale;
   float mag;
+  int mask;
+  PhotCode *code;
 
   initObjectData (catalog, Ncatalog);
@@ -67 +69 @@
   Nave = Npar = Npm = 0;
 
+  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+  Nsecfilt = GetPhotcodeNsecfilt ();
+  assert (catalog[0].Nsecfilt == Nsecfilt);
+
   for (i = 0; i < Ncatalog; i++) {
 
@@ -73 +79 @@
     Nskip = 0;
     for (j = 0; j < catalog[i].Naverage; j++) {
-
       /* calculate the average value of R,D for a single star */
+
+      // skip objects which are known to be problematic
+      // XXX include this code or not?
+      # if (0)
       if (catalog[i].average[j].code & STAR_BAD) {
         Nskip ++;
         continue;  
       }
+      # endif
 
       N = 0;
@@ -88 +98 @@
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 
-        /* exclude measurements by previous outlier detection */
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+        // exclude measurements by previous outlier detection
+        // XXX include this code or not?
+        # if (0)
+        if (catalog[i].measure[m].dbFlags & MEAS_BAD) { 
+          catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+          continue;
+        }
+        # endif
 
         /* exclude measurements by quality */
-        if (PhotFlagSelect && (catalog[i].measure[m].photFlags & PhotFlagValue)) continue;
+        if (PhotFlagSelect) {
+          if (PhotFlagBad) {
+            mask = PhotFlagBad;
+          } else {
+            code = GetPhotcodebyCode (catalog[i].measure[m].photcode);
+            mask = code[0].astromBadMask;
+          }
+          if (mask & catalog[i].measure[m].photFlags) {
+            catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+            continue;
+          }
+        }
         
         /* exclude measurements by mag limit */
         if (ImagSelect) {
           mag = PhotInst (&catalog[i].measure[m]);
-          if (mag < ImagMin) continue;
-          if (mag > ImagMax) continue;
-        }
-
-        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+          if (mag < ImagMin) {
+            catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+            continue;
+          }
+          if (mag > ImagMax) {
+            catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+            continue;
+          }
+        }
+
+        /* select or exclude measurements by photcode, or equiv photcode, if specified */
+        if (NphotcodesKeep > 0) {
+          found = FALSE;
+          for (kp = 0; (kp < NphotcodesKeep) && !found; kp++) {
+            if (photcodesKeep[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
+            if (photcodesKeep[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
+          }
+          if (!found) continue;
+        }
+        if (NphotcodesSkip > 0) {
+          found = FALSE;
+          for (kp = 0; (kp < NphotcodesSkip) && !found; kp++) {
+            if (photcodesSkip[kp][0].code == catalog[i].measure[m].photcode) found = TRUE;
+            if (photcodesSkip[kp][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode)) found = TRUE;
+          }
+          if (found) continue;
+        }
+
+        R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
         T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years
 
@@ -119 +170 @@
       // position, consider including the lower-quality detections
 
+      catalog[i].average[j].code &= ~ID_STAR_FEW;
+
       // XXX add the parallax factor range as a criterion as well
       if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE;
       if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
 
-      // XXX This criterion needs to be better considered: adjust to match Ndof
-      if (N < POS_TOOFEW) { /* too few measurements */
-        catalog[i].average[j].code |= ID_STAR_FEW;
+      // too few measurements for average position (require 2 values)
+      if (N < 2) {
+        // XXX need to define PHOTOM and ASTROM object flags
+        // catalog[i].average[j].code |= ID_STAR_FEW;
         continue;
-      } else {
-        catalog[i].average[j].code &= ~ID_STAR_FEW;
-      } 
+      }
 
       /* we need to do the fit in a locally linear space; choose a ref coordinate */
@@ -157 +209 @@
           fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
           fit.Nfit = N;
+
+          fit.uR = fit.duR = 0.0;
+          fit.uD = fit.duD = 0.0;
+          fit.p  = fit.dp  = 0.0;
+
           Nave ++;
           break;
@@ -167 +224 @@
           // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
           // continue;
+
+          fit.p  = fit.dp  = 0.0;
+
           Npm ++;
           break;
 
         case FIT_PAR_ONLY:
+          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+          exit (2);
+
           for (k = 0; k < N; k++) {
             ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
@@ -178 +241 @@
           // project Ro, Do back to RA,DEC
           XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
+
+          fit.uR = fit.duR = 0.0;
+          fit.uD = fit.duD = 0.0;
+
           Npar ++;
           break;
 
         case FIT_PM_AND_PAR:
+          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+          exit (2);
+
           for (k = 0; k < N; k++) {
             ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
@@ -207 +277 @@
       m = catalog[i].average[j].measureOffset;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+        // XXX why was this here?? if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
       }