Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 32695)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 33652)
@@ -38,4 +38,6 @@
 }  
 
+// This function operates on both Measure and MeasureTiny.  In the big stages, this should
+// be called with just MeasureTiny set and Measure == NULL
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
@@ -72,11 +74,14 @@
   T2000 = ohana_date_to_sec ("2000/01/01");
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+
   Nsecfilt = GetPhotcodeNsecfilt ();
-  assert (catalog[0].Nsecfilt == Nsecfilt);
+  if (Ncatalog) {
+    assert (catalog[0].Nsecfilt == Nsecfilt);
+  }
 
   NaveSum = NparSum = NpmSum = NskipSum = 0;
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
+    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
 
     Nave = Npar = Npm = Nskip = 0;
@@ -86,13 +91,4 @@
       XVERB = FALSE;
 
-      // 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
-
       if (catalog[i].average[j].Nmeasure == 0) {
 	  continue;
@@ -101,50 +97,56 @@
       N = 0;
       m = catalog[i].average[j].measureOffset;
-      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
+      MeasureTiny *measure = &catalog[i].measureT[m];
+      Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL;
+      // when we update the output measure values, we need to do it here
+
       mode = FIT_MODE;
 
       // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
-      Tmean = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-
-	//does the measurement pass the supplied filtering constraints?
-	if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) {
-	  catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+
+	// does the measurement pass the supplied filtering constraints?
+	if (!MeasFilterTestTiny(&measure[k], FALSE)) {
+	  measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+	  if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
 	  continue;
 	}
 
 	//outlier rejection
-	if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
-	  catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+	if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
+	  measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+	  if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
 	  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;
+	measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+	if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
+	R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+
+	// XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects
+	// to get Time in years relative to UNIX Tzero (1970/01/01)
+	T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
+
+	// dX, dY : error in arcsec -- 
+	dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
+	dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
+
+	// allow a given photcode or measurement to be
+	// ignored if the error is NAN (for photcode, set astromErrSys to NaN)
+	if (isnan(dX[N])) {
+	  measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
 	  continue;
 	}
-# endif
-
-	catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
-
-	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 - T2000) / (86400*365.25) ; // time relative to J2000 in years
-
-	Tmin = MIN(Tmin, T[N]);
-	Tmax = MAX(Tmax, T[N]);
-	Tmean += T[N];
-
-	// dX, dY : error in arcsec -- 
-	dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
-	dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+	if (isnan(dY[N])) {
+	  measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+	  continue;
+	}
 
 	// add systematic error in quadrature, if desired
 	// only do this after the fit has converged (or you will never improve the poor images)
 	// if (INCLUDE_SYS_ERR) {
-	// float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
+	// float dRsys = FromShortPixels(measure[k].dRsys);
 	// dX[N] = hypot(dX[N], dRsys);
 	// dY[N] = hypot(dY[N], dRsys);
@@ -154,5 +156,5 @@
 	// dY[N] = 0.1;
 
-	dT[N] = catalog[i].measure[m].dt;
+	dT[N] = measure[k].dt;
 
 	// XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
@@ -172,8 +174,22 @@
       catalog[i].average[j].flags &= ~ID_STAR_FEW;
 
+      // find Tmin & Tmax from the list of accepted measurements
+      Tmean = 0;
+      Tmin = Tmax = T[0];
+      for (k = 0; k < N; k++) {
+	Tmin = MIN(Tmin, T[k]);
+	Tmax = MAX(Tmax, T[k]);
+	Tmean += T[k];
+      }
       // XXX add the parallax factor range as a criterion as well
       Trange = Tmax - Tmin;
       if (Trange < PM_DT_MIN) mode = FIT_AVERAGE;
       if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
+
+      if (FIT_TARGET == TARGET_HIGH_SPEED) {
+	  Tmean = 0.5*(Tmax - Tmin);
+      } else {
+	  Tmean /= (float) N;
+      }
 
       // too few measurements for average position (require 2 values)
@@ -192,14 +208,7 @@
       coords.crval2 = D[0];
 
-      if (FIT_TARGET == TARGET_HIGH_SPEED) {
-	  Tmean = 0.5*(Tmax - Tmin);
-      } else {
-	  Tmean /= (float) N;
-      }
-      
       // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
-      // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
 
       // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -232,11 +241,25 @@
 	// fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__);
 
+	float pXmin = +2.0;
+	float pXmax = -2.0;
+	float pYmin = +2.0;
+	float pYmax = -2.0;
 	for (k = 0; k < N; k++) {
 	  ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
-	}
-	FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
-	XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
-	catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
-	Npar ++;
+	  pXmin = MIN (pXmin, pX[k]);
+	  pXmax = MAX (pXmax, pX[k]);
+	  pYmin = MIN (pYmin, pY[k]);
+	  pYmax = MAX (pYmax, pY[k]);
+	}
+	float dXRange = pXmax - pXmin;
+	float dYRange = pYmax - pYmin;
+	float parRange = hypot (dXRange, dYRange);
+	
+	if (parRange >= PAR_FACTOR_MIN) {
+	  FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
+	  XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+	  catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
+	  Npar ++;
+	}
       }	  
 
@@ -307,8 +330,11 @@
 
       // the measure fields must be updated before the average fields
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	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]);
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+	setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	if (measureBig) {
+	  setMeanR_Big (fit.Ro, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	  setMeanD_Big (fit.Do, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	}
       }