Changeset 37261

Timestamp:

Aug 14, 2014, 2:25:44 PM (12 years ago)

Author:

eugene

Message:

fix stack R,D in relastro

Location:

trunk/Ohana/src/relastro

Files:

• Makefile (modified) (2 diffs)
• include/relastro.h (modified) (2 diffs)
• src/ParFactor.c (modified) (1 diff)
• src/UpdateObjects.c (modified) (4 diffs)
• src/extra.c (added)

trunk/Ohana/src/relastro/Makefile

@@ -59 +59 @@
 $(SRC)/save_catalogs.$(ARCH).o       \
 $(SRC)/CoordOps.$(ARCH).o            \
+$(SRC)/extra.$(ARCH).o            \
 $(SRC)/FixProblemImages.$(ARCH).o    \
 $(SRC)/StarMaps.$(ARCH).o    \
@@ -117 +118 @@
 $(SRC)/save_catalogs.$(ARCH).o       \
 $(SRC)/CoordOps.$(ARCH).o            \
+$(SRC)/extra.$(ARCH).o            \
 $(SRC)/high_speed_catalogs.$(ARCH).o  \
 $(SRC)/high_speed_objects.$(ARCH).o  \

trunk/Ohana/src/relastro/include/relastro.h

@@ -133 +133 @@
   int    Nmeas;
 } StatType;
+
+typedef struct {
+  off_t Nave;
+  off_t Npm;
+  off_t Npar;
+  off_t Nskip;
+  off_t Noffset;
+} FitStats;
 
 /* global variables set in parameter file */
@@ -535 +543 @@
 int areImagesLoaded ();
 int areImagesMatched ();
+
+int isGPC1chip (int photcode);
+int isGPC1stack (int photcode);
+int isGPC1warp (int photcode);

trunk/Ohana/src/relastro/src/ParFactor.c

@@ -63 +63 @@
   double jd, lambda, beta, epsilon, Radius;
 
-  /* given a time T in UNIX seconds, determine the solar longitude S */
+  /* given a Time relative to Tmean, Tmean in years since J2000, determine the solar
+    longitude S */
 
   // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));

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

@@ -1 +1 @@
 # include "relastro.h"
+# define PAR_TOOFEW 5
+
+int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m);
+int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
 
 static off_t Nmax;
@@ -13 +17 @@
 static double *C_red;
 
+static Coords coords;
+
+static time_t T2000;
+
+void initFitStats (FitStats *fitStats) {
+  fitStats->Nave = 0;  
+  fitStats->Npm = 0;   
+  fitStats->Npar = 0;  
+  fitStats->Nskip = 0; 
+  fitStats->Noffset = 0;
+  return;
+}
+
+void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) {
+  tgtFitStats->Nave    += srcFitStats->Nave    ;  
+  tgtFitStats->Npm     += srcFitStats->Npm     ;   
+  tgtFitStats->Npar    += srcFitStats->Npar    ;  
+  tgtFitStats->Nskip   += srcFitStats->Nskip   ; 
+  tgtFitStats->Noffset += srcFitStats->Noffset ;
+  return;
+}
+
 void initObjectData (Catalog *catalog, int Ncatalog) {
 
@@ -41 +67 @@
   ALLOCATE (C_blue, double, MAX (1, Nmax));
   ALLOCATE (C_red,  double, MAX (1, Nmax));
-}  
-
-void freeObjectData () {
-
-  free (R);
-  free (D);
-  free (T);
-  free (X);
-  free (Y);
-
-  free (dR);
-  free (dD);
-  free (dT);
-  free (dX);
-  free (dY);
-
-  free (pX);
-  free (pY);
-
-  free (C_blue);
-  free (C_red);
-}  
-
-// 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) {
-
-  off_t j, k, m;
-  int i, N, NcBlue, NcRed, Nsecfilt, mode, result, status, XVERB;
-  StatType statsR, statsD;
-  Coords coords;
-  PMFit fitAve, fitPM, fitPAR, fit;
-  time_t T2000;
-  off_t Nave, Npm, Npar, Nskip, Noffset;
-  off_t NaveSum, NpmSum, NparSum, NskipSum, NoffSum;
-  double Tmin, Tmax, Tmean, Trange;
-
-  memset (&fit,    0, sizeof(fit));
-  memset (&fitAve, 0, sizeof(fitAve));
-  memset (&fitPM,  0, sizeof(fitPM));
-  memset (&fitPAR, 0, sizeof(fitPAR));
-  initObjectData (catalog, Ncatalog);
-
-  int setRefColor = areImagesMatched();
 
   /* project coordinates to a plane centered on the object with units of arcsec */
@@ -97 +79 @@
   strcpy (coords.ctype, "DEC--SIN");
 
-  XVERB = FALSE;
-
   // use J2000 as a reference time
-  T2000 = ohana_date_to_sec ("2000/01/01");
+  T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
+}  
+
+void freeObjectData () {
+
+  free (R);
+  free (D);
+  free (T);
+  free (X);
+  free (Y);
+
+  free (dR);
+  free (dD);
+  free (dT);
+  free (dX);
+  free (dY);
+
+  free (pX);
+  free (pY);
+
+  free (C_blue);
+  free (C_red);
+}  
+
+// 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) {
+
+  initObjectData (catalog, Ncatalog);
+
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
-
-  Nsecfilt = GetPhotcodeNsecfilt ();
+  int Nsecfilt = GetPhotcodeNsecfilt ();
   if (Ncatalog) {
     assert (catalog[0].Nsecfilt == Nsecfilt);
   }
 
-  NaveSum = NparSum = NpmSum = NoffSum = NskipSum = 0;
+  FitStats sumStatsChips; initFitStats (&sumStatsChips);
+  FitStats sumStatsStack; initFitStats (&sumStatsStack);
+
+  int i;
   for (i = 0; i < Ncatalog; i++) {
 
     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 = Noffset = 0;
+    FitStats fitStatsChips; initFitStats (&fitStatsChips);
+    FitStats fitStatsStack; initFitStats (&fitStatsStack);
+
+    off_t j;
     for (j = 0; j < catalog[i].Naverage; j++) {
       /* calculate the average value of R,D for a single star */
-
-      XVERB = FALSE;
-      fitAve.chisq = NAN;
-      fitPM.chisq = NAN;
-      fitPAR.chisq = NAN;
-
-      // if we fail to fit the astrometry for some reason, we need to set/reset these
-      catalog[i].average[j].flags |= ID_STAR_NO_ASTROM;
-      catalog[i].average[j].ChiSqAve  = NAN;
-      catalog[i].average[j].ChiSqPM   = NAN;
-      catalog[i].average[j].ChiSqPar  = NAN;
-      catalog[i].average[j].Npos = 0;
-
-      if (catalog[i].average[j].Nmeasure == 0) {
+      off_t m = catalog[i].average[j].measureOffset;
+      MeasureTiny *measure = &catalog[i].measureT[m];
+      Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL;
+      Average *average = &catalog[i].average[j];
+      SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt];
+
+      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack);
+      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m);
+    }
+    if (VERBOSE) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
+    if (VERBOSE) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
+    sumFitStats (&fitStatsChips, &sumStatsChips);
+    sumFitStats (&fitStatsStack, &sumStatsStack);
+  }
+  freeObjectData ();
+
+  if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
+  if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
+  return (TRUE);
+}
+
+// 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_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m) {
+
+  int setRefColor = areImagesMatched();
+
+  /* calculate the average value of R,D for a single star */
+
+  PMFit fit;    memset (&fit,    0, sizeof(fit));
+  PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN;
+  PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq = NAN;
+  PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;
+
+  // if we fail to fit the astrometry for some reason, we need to set/reset these
+  average[0].flags |= ID_STAR_NO_ASTROM;
+  average[0].ChiSqAve  = NAN;
+  average[0].ChiSqPM   = NAN;
+  average[0].ChiSqPar  = NAN;
+  average[0].Npos = 0;
+
+  // an object with no measurements is externally supplied
+  if (average[0].Nmeasure == 0) return TRUE;
+
+  int NcBlue = 0;
+  int NcRed = 0;
+  int N = 0;
+
+  int mode = FIT_MODE; // start with the globally-defined fit mode
+
+  int XVERB = FALSE;
+  XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC);
+  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
+
+  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
+  off_t k;
+  for (k = 0; k < average[0].Nmeasure; k++) {
+
+    if (XVERB) {
+      char *date = ohana_sec_to_date (measure[k].t);
+      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
+      fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
+      free (date);
+    }
+
+    // SKIP gpc1 stack data
+    if (isGPC1stack(measure[k].photcode)) continue;
+
+    // SKIP gpc1 forced-warp data
+    if (isGPC1warp(measure[k].photcode)) continue;
+
+    // reset the bit to note that a detection was used (or not)
+    measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+    if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
+
+    // does the measurement pass the supplied filtering constraints?
+    // MeasFilterTestTiny does not test psfQF
+    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
+    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
+    if (!keepMeasure) {
+      continue;
+    }
+
+    double Ri = getMeanR (&measure[k], average, secfilt);
+    double Di = getMeanD (&measure[k], average, secfilt);
+
+    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
+    if (ExcludeBogus) {
+      coords.crval1 = average[0].R;
+      coords.crval2 = average[0].D;
+      double Xi, Yi;
+      RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
+      double radius = hypot(Xi, Yi);
+      if (radius > ExcludeBogusRadius) {
+        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
+        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
         continue;
       }
-
-      NcBlue = 0;
-      NcRed = 0;
-      N = 0;
-      m = catalog[i].average[j].measureOffset;
-      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;
-
-      // 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);
-
-      // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
-
-        if (XVERB) {
-          char *date = ohana_sec_to_date (measure[k].t);
-          int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
-          fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
-          free (date);
-        }
-
-        // reset the bit to note that a detection was used (or not)
-        measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
-        if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
-
-        // does the measurement pass the supplied filtering constraints?
-        // MeasFilterTestTiny does not test psfQF
-        // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
-        int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
-        if (!keepMeasure) {
-          continue;
-        }
-
-        // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
-        if (ExcludeBogus) {
-          double Ri = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-          double Di = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-          coords.crval1 = catalog[i].average[j].R;
-          coords.crval2 = catalog[i].average[j].D;
-          double Xi, Yi;
-          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
-          double radius = hypot(Xi, Yi);
-          if (radius > ExcludeBogusRadius) {
-              measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
-              if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
-              continue;
-          }
-          measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
-          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
-        }
-
-        // outlier rejection
-        if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
-          continue;
-        }
-
-        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])) continue;
-        if (isnan(dY[N])) 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(measure[k].dRsys);
-        // dX[N] = hypot(dX[N], dRsys);
-        // dY[N] = hypot(dY[N], dRsys);
-        // }
-
-        // dX[N] = 0.1;
-        // dY[N] = 0.1;
-
-        dT[N] = measure[k].dt;
-
-        // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
-        // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
-        // non-circular errors (different values for X and Y), then dR and dD will be
-        // incorrect: they would need to be rotated to take out the position angle
-        dR[N] = dX[N] / 3600.0;
-        dD[N] = dY[N] / 3600.0;
-
-        if (setRefColor) {
-          float colorBlue = getColorBlue (m+k, i);
-          if (!isnan(colorBlue)) {
-            C_blue[NcBlue] = colorBlue;
-            NcBlue++;
-          }
-          float colorRed = getColorRed (m+k, i);
-          if (!isnan(colorRed)) {
-            C_red[NcRed] = colorRed;
-            NcRed++;
-          }
-        }
-
-        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
-        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
-
-        N++;
-      } // loop over measurements : catalog[i].average[j].Nmeasure 
-
-      // if we have too few good detections for the desired fit, or too limited a
-      // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
-      // for even the average position, consider including the lower-quality detections?
-
-      // 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];
+      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+      if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
+    }
+
+    // outlier rejection
+    if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
+      continue;
+    }
+
+    R[N] = Ri;
+    D[N] = Di;
+
+    // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
+    // T[] is time in years since J2000 (jd = 2451545)
+    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])) continue;
+    if (isnan(dY[N])) 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(measure[k].dRsys);
+    // dX[N] = hypot(dX[N], dRsys);
+    // dY[N] = hypot(dY[N], dRsys);
+    // }
+
+    // dX[N] = 0.1;
+    // dY[N] = 0.1;
+
+    dT[N] = measure[k].dt;
+
+    // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
+    // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
+    // non-circular errors (different values for X and Y), then dR and dD will be
+    // incorrect: they would need to be rotated to take out the position angle
+    dR[N] = dX[N] / 3600.0;
+    dD[N] = dY[N] / 3600.0;
+
+    if (setRefColor) {
+      float colorBlue = getColorBlue (m+k, i);
+      if (!isnan(colorBlue)) {
+        C_blue[NcBlue] = colorBlue;
+        NcBlue++;
       }
-      // 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) || (mode == FIT_PM_AND_PAR)) && (N <= PM_TOOFEW)) mode = FIT_AVERAGE;
-
-      if (RELASTRO_OP == OP_HIGH_SPEED) {
-          Tmean = 0.5*(Tmax - Tmin);
-      } else {
-          Tmean /= (float) N;
+      float colorRed = getColorRed (m+k, i);
+      if (!isnan(colorRed)) {
+        C_red[NcRed] = colorRed;
+        NcRed++;
       }
-
-      // too few measurements for average position (require 2 values)
-      if (N < SRC_MEAS_TOOFEW) {
-        if (N < 2) continue;
+    }
+
+    measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+    if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
+    N++;
+  } // loop over measurements : average[0].Nmeasure 
+
+  if (N < 1) { 
+    if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
+      average[0].R  = average[0].Rstk;
+      average[0].D  = average[0].Dstk;
+      average[0].dR = average[0].dRstk;
+      average[0].dD = average[0].dDstk;
+    }
+    return FALSE;
+  }
+
+  // if we have too few good detections for the desired fit, or too limited a
+  // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
+  // for even the average position, consider including the lower-quality detections?
+
+  // find Tmin & Tmax from the list of accepted measurements
+  double Tmean = 0.0;
+  double Tmin = T[0];
+  double Tmax = T[0];
+  for (k = 0; k < N; k++) {
+    Tmin = MIN(Tmin, T[k]);
+    Tmax = MAX(Tmax, T[k]);
+    Tmean += T[k];
+  }
+  double Trange = Tmax - Tmin;
+
+  if (RELASTRO_OP == OP_HIGH_SPEED) {
+    Tmean = 0.5*(Tmax - Tmin);
+  } else {
+    Tmean /= (float) N;
+  }
+
+  /* we need to do the fit in a locally linear space; choose a ref coordinate */
+  coords.crval1 = R[0];
+  coords.crval2 = D[0];
+
+  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
+
+  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
+  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
+    if (Trange < PM_DT_MIN) {
+      mode = FIT_AVERAGE;
+      goto skipPM;
+    }
+    if (N <= PM_TOOFEW) {
+      mode = FIT_AVERAGE;
+      goto skipPM;
+    }
+
+    // project all of the R,D coordinates to a plane centered on this coordinate. set
+    // the times to be relative to Tmean (this is required for parallax as well)
+    for (k = 0; k < N; k++) {
+      RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
+      T[k] -= Tmean;
+      if (XVERB) {
+        fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
       }
-
-      /* we need to do the fit in a locally linear space; choose a ref coordinate */
-      coords.crval1 = R[0];
-      coords.crval2 = D[0];
-
-      // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
-
-      if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
-        // project all of the R,D coordinates to a plane centered on this coordinate. set
-        // the times to be relative to Tmean (this is required for parallax as well)
-        for (k = 0; k < N; k++) {
-          RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
-          T[k] -= Tmean;
-          if (XVERB) {
-            fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
-          }
-        }         
-
-        FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
-
-        if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
-
-        // project Ro, Do back to RA,DEC
-        XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
-        if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
-
-        fitPM.p  = fitPM.dp  = 0.0;
-        catalog[i].average[j].flags |= ID_STAR_FIT_PM;
-        Npm ++;
-
-        // XXX a hard-wired hack...
-        if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
-          mode = FIT_AVERAGE;
-          catalog[i].average[j].flags |= ID_STAR_BAD_PM;
-        }
-      }
-
-      if (mode == FIT_PM_AND_PAR) {
-        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);
-          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);
+    }     
+
+    FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
+
+    if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
+
+    // project Ro, Do back to RA,DEC
+    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
+    if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
+
+    fitPM.p  = fitPM.dp  = 0.0;
+    average[0].flags |= ID_STAR_FIT_PM;
+    fitStats->Npm ++;
+
+    // XXX a hard-wired hack...
+    if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
+      mode = FIT_AVERAGE;
+      average[0].flags |= ID_STAR_BAD_PM;
+    }
+  }
+  
+skipPM:
+  // fit the parallax + proper-motion model
+  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
+  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
+  if (mode == FIT_PM_AND_PAR) {
+    if (Trange < PM_DT_MIN) {
+      mode = FIT_PM_ONLY;
+      goto skipPAR;
+    }
+    if (N <= PAR_TOOFEW) {
+      mode = FIT_PM_ONLY;
+      goto skipPAR;
+    }
+    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);
+      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);
         
-# define PAR_TOOFEW 5
-        if ((parRange >= PAR_FACTOR_MIN) && (N > PAR_TOOFEW)) {
-          FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
-          if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);
-
-          XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
-          catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
-          Npar ++;
-
-          // XXX a hard-wired hack...
-          if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
-            mode = FIT_AVERAGE;
-            catalog[i].average[j].flags |= ID_STAR_BAD_PM;
-          }
-        } else {
-          // need to set mode = FIT_PM_ONLY if we do not fit for parallax
-          mode = FIT_PM_ONLY;
-        }
-      }   
-
-      // fit the average model
-      if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
-        liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
-        liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
-
-        fitAve.Ro = statsR.mean;
-        fitAve.dRo = 3600.0*statsR.sigma;
-
-        fitAve.Do = statsD.mean;
-        fitAve.dDo = 3600.0*statsD.sigma;
-
-        fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
-        fitAve.Nfit = N;
-
-        fitAve.uR = fitAve.duR = 0.0;
-        fitAve.uD = fitAve.duD = 0.0;
-        fitAve.p  = fitAve.dp  = 0.0;
-        catalog[i].average[j].flags |= ID_STAR_FIT_AVE;
-        Nave ++;
-      }
-
-      if (setRefColor) {
-        float colorMedian;
-        dsort (C_blue, NcBlue);
-        colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
-        catalog[i].average[j].refColorBlue = colorMedian;
-        dsort (C_red, NcRed);
-        colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
-        catalog[i].average[j].refColorRed = colorMedian;
-      }
-
-      /* choose the result based on the chisq values */
-      // XXXX for now, just use the mode as the result:
-      result = mode;
-
-      switch (result) {
-      case FIT_AVERAGE:
-        catalog[i].average[j].flags |= ID_STAR_USE_AVE;
-        fit = fitAve;
-        break;
-      case FIT_PM_ONLY:
-        catalog[i].average[j].flags |= ID_STAR_USE_PM;
-        fit = fitPM;
-        break;
-      case FIT_PM_AND_PAR:
-        catalog[i].average[j].flags |= ID_STAR_USE_PAR;
-        fit = fitPAR;
-        break;
-      }
-      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
-                          catalog[i].average[j].R, 
-                          catalog[i].average[j].D, 
-                          fit.Ro, fit.Do, 
-                          3600*(catalog[i].average[j].R - fit.Ro), 
-                          3600*(catalog[i].average[j].D - fit.Do),
-                          fit.uR, fit.uD, fit.p, fit.dp);
-
-      // make sure that the fit succeeded
-      status = TRUE;
-      status &= finite(fit.Ro);
-      status &= finite(fit.Do);
-      status &= finite(fit.dRo);
-      status &= finite(fit.dDo);
-      status &= finite(fit.uR);
-      status &= finite(fit.uD);
-      status &= finite(fit.duR);
-      status &= finite(fit.duD);
-      status &= finite(fit.p);
-      status &= finite(fit.dp);
-      if (!status) {
-        Nskip ++;
-        continue;
-      }
-
-      // what is the offset relative to the mean fit position?
-      coords.crval1 = catalog[i].average[j].R;
-      coords.crval2 = catalog[i].average[j].D;
-
-      double dXoff, dYoff;
-      RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
-      float dPos = hypot (dXoff, dYoff);
-      if (dPos > MaxMeanOffset) {
-        if (Noffset < 100) {
-          fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
-        }
-        Noffset ++;
-        continue;
-      }
-
-
-      // the measure fields must be updated before the average fields
-      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]);
-        }
-      }      
-
-      catalog[i].average[j].R   = fit.Ro; // RA in degrees
-      catalog[i].average[j].D   = fit.Do; // DEC in degrees
-      catalog[i].average[j].dR  = fit.dRo; // RA scatter in arcsec
-      catalog[i].average[j].dD  = fit.dDo; // DEC scatter in arcsec
-
-      catalog[i].average[j].uR  = fit.uR; // RA proper motion in arcsec/year
-      catalog[i].average[j].uD  = fit.uD; // DEC proper motion in arcsec/year
-      catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year
-      catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year
-
-      catalog[i].average[j].P   = fit.p; // parallax in arcsec
-      catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec
-
-      catalog[i].average[j].ChiSqAve  = fitAve.chisq;
-      catalog[i].average[j].ChiSqPM   = fitPM.chisq;
-      catalog[i].average[j].ChiSqPar  = fitPAR.chisq;
-      catalog[i].average[j].Tmean = (Tmean * 86400 * 365.25) + T2000;
-      catalog[i].average[j].Trange = (Trange * 86400 * 365.25);
-      catalog[i].average[j].Npos = fit.Nfit;
-
-      // XXX EAM 20140812: for a test, set average.Rstk,Dstk, etc to match R,D
-      catalog[i].average[j].Rstk  = fit.Ro; // RA in degrees
-      catalog[i].average[j].Dstk  = fit.Do; // DEC in degrees
-      catalog[i].average[j].dRstk = fit.dRo; // RA scatter in arcsec
-      catalog[i].average[j].dDstk = fit.dDo; // DEC scatter in arcsec
-
-      // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
-      catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM;
-
-      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
-                          catalog[i].average[j].R,
-                          catalog[i].average[j].D,
-                          fit.Ro, fit.Do,
-                          3600*(catalog[i].average[j].R - fit.Ro),
-                          3600*(catalog[i].average[j].D - fit.Do),
-                          catalog[i].average[j].uR,
-                          catalog[i].average[j].uD,
-                          fitAve.chisq, fitPM.chisq, fitPAR.chisq);
-    }
-
-    NaveSum += Nave;
-    NpmSum += Npm;
-    NparSum += Npar;
-    NskipSum += Nskip;
-    NoffSum += Noffset;
-    if (VERBOSE) fprintf (stderr, "catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  Nave,  Npm,  Npar,  Nskip, Noffset);
-  }
-
-  freeObjectData ();
-
-  if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (NaveSum + NpmSum + NparSum),  NaveSum,  NpmSum,  NparSum,  NskipSum, NoffSum);
+    if (parRange < PAR_FACTOR_MIN) {
+      mode = FIT_PM_ONLY;
+      goto skipPAR;
+    }
+
+    FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
+    if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);
+
+    XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+    average[0].flags |= ID_STAR_FIT_PAR;
+    fitStats->Npar ++;
+
+    // XXX a hard-wired hack...
+    if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
+      mode = FIT_PM_ONLY;
+    }
+  }       
+
+skipPAR:
+  {
+    // ALWAYS fit the average model
+    StatType statsR, statsD;
+    liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
+    liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
+
+    fitAve.Ro = statsR.mean;
+    fitAve.dRo = 3600.0*statsR.sigma;
+
+    fitAve.Do = statsD.mean;
+    fitAve.dDo = 3600.0*statsD.sigma;
+
+    fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
+    fitAve.Nfit = N;
+
+    fitAve.uR = fitAve.duR = 0.0;
+    fitAve.uD = fitAve.duD = 0.0;
+    fitAve.p  = fitAve.dp  = 0.0;
+    average[0].flags |= ID_STAR_FIT_AVE;
+    fitStats->Nave ++;
+  }
+
+  if (setRefColor) {
+    float colorMedian;
+    dsort (C_blue, NcBlue);
+    colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
+    average[0].refColorBlue = colorMedian;
+    dsort (C_red, NcRed);
+    colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
+    average[0].refColorRed = colorMedian;
+  }
+
+  /* choose the result based on the chisq values */
+  // XXXX for now, just use the mode as the result:
+  int result = mode;
+
+  switch (result) {
+    case FIT_AVERAGE:
+      average[0].flags |= ID_STAR_USE_AVE;
+      fit = fitAve;
+      break;
+    case FIT_PM_ONLY:
+      average[0].flags |= ID_STAR_USE_PM;
+      fit = fitPM;
+      break;
+    case FIT_PM_AND_PAR:
+      average[0].flags |= ID_STAR_USE_PAR;
+      fit = fitPAR;
+      break;
+  }
+  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
+                      average[0].R, 
+                      average[0].D, 
+                      fit.Ro, fit.Do, 
+                      3600*(average[0].R - fit.Ro), 
+                      3600*(average[0].D - fit.Do),
+                      fit.uR, fit.uD, fit.p, fit.dp);
+
+  // make sure that the fit succeeded
+  int status = TRUE;
+  status &= finite(fit.Ro);
+  status &= finite(fit.Do);
+  status &= finite(fit.dRo);
+  status &= finite(fit.dDo);
+  status &= finite(fit.uR);
+  status &= finite(fit.uD);
+  status &= finite(fit.duR);
+  status &= finite(fit.duD);
+  status &= finite(fit.p);
+  status &= finite(fit.dp);
+  if (!status) {
+    fitStats->Nskip ++;
+    return FALSE; // XXX ??
+  }
+
+  // what is the offset relative to the mean fit position?
+  coords.crval1 = average[0].R;
+  coords.crval2 = average[0].D;
+
+  double dXoff, dYoff;
+  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
+  float dPos = hypot (dXoff, dYoff);
+  if (dPos > MaxMeanOffset) {
+    if (fitStats->Noffset < 100) {
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
+    }
+    fitStats->Noffset ++;
+    return FALSE; // XXX ??
+  }
+
+  average[0].R          = fit.Ro; // RA in degrees
+  average[0].D          = fit.Do; // DEC in degrees
+  average[0].dR         = fit.dRo; // RA scatter in arcsec
+  average[0].dD         = fit.dDo; // DEC scatter in arcsec
+
+  average[0].uR         = fit.uR; // RA proper motion in arcsec/year
+  average[0].uD         = fit.uD; // DEC proper motion in arcsec/year
+  average[0].duR        = fit.duR; // RA proper motion error in arcsec/year
+  average[0].duD        = fit.duD; // DEC proper motion error in arcsec/year
+
+  average[0].P          = fit.p; // parallax in arcsec
+  average[0].dP         = fit.dp; // parallax error in arcsec
+
+  average[0].ChiSqAve   = fitAve.chisq;
+  average[0].ChiSqPM    = fitPM.chisq;
+  average[0].ChiSqPar   = fitPAR.chisq;
+
+  average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
+  average[0].Trange     = (Trange * 86400 * 365.25);
+  average[0].Npos       = fit.Nfit;
+
+  // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
+  average[0].flags &= ~ID_STAR_NO_ASTROM;
+
+  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
+                      average[0].R,
+                      average[0].D,
+                      fit.Ro, fit.Do,
+                      3600*(average[0].R - fit.Ro),
+                      3600*(average[0].D - fit.Do),
+                      average[0].uR,
+                      average[0].uD,
+                      fitAve.chisq, fitPM.chisq, fitPAR.chisq);
+
   return (TRUE);
 }
+
+// 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_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {
+
+  off_t k;
+
+  // set the default values
+  average[0].Rstk  = NAN; // RA in degrees
+  average[0].Dstk  = NAN; // DEC in degrees
+  average[0].dRstk = NAN; // RA scatter in arcsec
+  average[0].dDstk = NAN; // DEC scatter in arcsec
+
+  /* calculate the average value of R,D for a single star */
+  PMFit fitAve;
+  memset (&fitAve, 0, sizeof(fitAve));
+  fitAve.chisq = NAN;
+
+  if (average[0].Nmeasure == 0) return TRUE;
+
+  int N = 0;
+
+  int XVERB = FALSE;
+  XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC);
+  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
+
+  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
+  for (k = 0; k < average[0].Nmeasure; k++) {
+
+    if (XVERB) {
+      char *date = ohana_sec_to_date (measure[k].t);
+      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
+      fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
+      free (date);
+    }
+
+    // SKIP everything except gpc1 stack data
+    if (!isGPC1stack(measure[k].photcode)) continue;
+
+    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
+    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
+    if (!keepMeasure) {
+      continue;
+    }
+
+    R[N] = getMeanR (&measure[k], average, secfilt);
+    D[N] = getMeanD (&measure[k], average, secfilt);
+
+    // 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])) continue;
+    if (isnan(dY[N])) continue;
+
+    // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
+    // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
+    // non-circular errors (different values for X and Y), then dR and dD will be
+    // incorrect: they would need to be rotated to take out the position angle
+    dR[N] = dX[N] / 3600.0;
+    dD[N] = dY[N] / 3600.0;
+
+    // XXX use a different flag for stack measurements?
+    // measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+    // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
+    N++;
+  } // loop over measurements : average[0].Nmeasure 
+
+  // if we have too few good detections for the desired fit, or too limited a
+  // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
+  // for even the average position, consider including the lower-quality detections?
+
+  // too few measurements for average position (require 2 values)
+  if (N < 1) return FALSE; // XXX ?? 
+
+  // find the mean position
+  StatType statsR, statsD;
+  liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
+  liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
+
+  fitAve.Ro = statsR.mean;
+  fitAve.dRo = 3600.0*statsR.sigma;
+
+  fitAve.Do = statsD.mean;
+  fitAve.dDo = 3600.0*statsD.sigma;
+
+  fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
+  fitAve.Nfit = N;
+
+  // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
+  fitStats->Nave ++;
+
+  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
+                      average[0].R, 
+                      average[0].D, 
+                      fitAve.Ro, fitAve.Do, 
+                      3600*(average[0].R - fitAve.Ro), 
+                      3600*(average[0].D - fitAve.Do));
+
+  // make sure that the fit succeeded
+  int status = TRUE;
+  status &= finite(fitAve.Ro);
+  status &= finite(fitAve.Do);
+  status &= finite(fitAve.dRo);
+  status &= finite(fitAve.dDo);
+  if (!status) {
+    fitStats->Nskip ++;
+    return FALSE;
+  }
+
+  // what is the offset relative to the mean fit position?
+  coords.crval1 = average[0].R;
+  coords.crval2 = average[0].D;
+
+  double dXoff, dYoff;
+  RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords);
+  float dPos = hypot (dXoff, dYoff);
+  if (dPos > MaxMeanOffset) {
+    if (fitStats->Noffset < 100) {
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
+    }
+    fitStats->Noffset ++;
+    return FALSE;
+  }
+
+  // set the stack position values
+  average[0].Rstk  = fitAve.Ro; // RA in degrees
+  average[0].Dstk  = fitAve.Do; // DEC in degrees
+  average[0].dRstk = fitAve.dRo; // RA scatter in arcsec
+  average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec
+
+  return (TRUE);
+}
+
+
 
 /* fitting proper-motion and parallax: