Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 17210)
+++ /trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 17211)
@@ -40,5 +40,5 @@
 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,4 +49,6 @@
   float errorScale;
   float mag;
+  int mask;
+  PhotCode *code;
 
   initObjectData (catalog, Ncatalog);
@@ -67,4 +69,8 @@
   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,10 +79,14 @@
     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,19 +98,60 @@
       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,15 +170,16 @@
       // 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,4 +209,9 @@
 	  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,8 +224,14 @@
 	  // 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,8 +241,15 @@
 	  // 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,7 +277,7 @@
       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]);
       }