Index: /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/average.d	(revision 38599)
@@ -69,5 +69,9 @@
 FIELD catID,          CAT_ID,      unsigned int,    unique ID for table in which object was first realized
 FIELD extID,          EXT_ID,      uint64_t,        external ID for object (eg PSPS objID)
-FIELD extIDgc,        EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
+
+# replace extIDgc (unused) with uRgal, uDgal:
+# FIELD extIDgc,      EXT_ID_GC,   uint64_t,        external ID for object in galactic coords
+FIELD uRgal,          U_RA_GAL,    float,           modeled proper motion based on galactic motion
+FIELD uDgal,          U_DEC_GAL,   float,           modeled proper motion based on galactic motion
 
 # this structure should only be used for internal representations
Index: /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/measure.d
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/measure.d	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/libautocode/def/measure.d	(revision 38599)
@@ -42,4 +42,7 @@
 FIELD XoffCAM,        X_OFF_CAM,     float,          X offset from correction,     pixels
 FIELD YoffCAM,        Y_OFF_CAM,     float,          Y offset from correction,     pixels
+
+# XXX I can deprecate these as I am going to apply the correct uR,uD offset
+# XXX not sure how to use this yet...
 FIELD RoffGAL,        R_OFF_GAL,     float,          RA offset from correction,    arcsec
 FIELD DoffGAL,        D_OFF_GAL,     float,          DEC offset from correction,   arcsec
Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/include/dvo.h	(revision 38599)
@@ -171,7 +171,9 @@
   ID_PROPER            = 0x00000010, // star with large proper motion
   ID_TRANSIENT         = 0x00000020, // identified as a non-periodic (stationary) transient
-  ID_VARIABLE          = 0x00000040, // identified as a period variable
+  ID_VARIABLE          = 0x00000040, // identified as a periodic variable
   ID_ASTEROID          = 0x00000080, // identified with a known solar-system object (asteroid or other)
-  // bits 0x00000100 - 0x00008000 are currently unused
+  ID_STACK_ASTROM      = 0x00000100, // stack position used for astrometry
+  
+  // bits 0x00000200 - 0x00008000 are currently unused
   ID_STAR_FIT_AVE      = 0x00010000, // average position fitted
   ID_STAR_FIT_PM       = 0x00020000, // proper motion fitted
Index: /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_catalog_create.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_catalog_create.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/libdvo/src/dvo_catalog_create.c	(revision 38599)
@@ -33,5 +33,5 @@
     catalog[0].lensobj_catalog 	= dvo_catalog_create_subcat (catalog, "cpy", "LENSOBJ");
     catalog[0].starpar_catalog 	= dvo_catalog_create_subcat (catalog, "cpz", "STARPAR");
-    catalog[0].galphot_catalog = dvo_catalog_create_subcat (catalog, "cpq", "GALPHOT");
+    catalog[0].galphot_catalog  = dvo_catalog_create_subcat (catalog, "cpq", "GALPHOT");
 
     // lock the additional split files
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h	(revision 38599)
@@ -115,13 +115,54 @@
   double Ro, dRo;
   double Do, dDo;
-
   double uR, duR;
   double uD, duD;
-
-  double p, dp;
-
+  double  p, dp;
+
+  int getChisq;
   double chisq;
   int Nfit;
-} PMFit;
+} FitAstromResult;
+
+typedef struct {
+  double **A;
+  double **B;
+  int Nterms;
+} FitAstromData;
+
+typedef struct {
+  double X, dX;
+  double Y, dY;
+  double R, dR;
+  double D, dD;
+  double T, dT;
+  double pX;
+  double pY;
+  double C_blue;
+  double C_red;
+} FitAstromPoints;
+
+typedef struct {
+  off_t Nave;
+  off_t Npm;
+  off_t Npar;
+  off_t Nskip;
+  off_t Noffset;
+
+  FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats
+  int Nfit;
+  int NfitAlloc;
+
+  FitAstromPoints *points;
+  FitAstromPoints *sample;
+  int Npoints;
+  int NpointsAlloc;
+
+  FitAstromData *fitPos;
+  FitAstromData *fitPM;
+  FitAstromData *fitPar;
+
+  Coords coords;
+  time_t T2000;
+} FitStats;
 
 typedef struct {
@@ -148,12 +189,4 @@
   int    Nmeas;
 } StatType;
-
-typedef struct {
-  off_t Nave;
-  off_t Npm;
-  off_t Npar;
-  off_t Nskip;
-  off_t Noffset;
-} FitStats;
 
 typedef struct {
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c	(revision 38599)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c	(revision 38599)
@@ -0,0 +1,147 @@
+# include "relastro.h"
+
+// I am modifying FitPM with an eye to (a) threaded operations and (b) bootstrap resampling tests.
+
+FitStats *FitStatsInit (int Nmax, int Nboot) {
+
+  // counters to record successful fits or failures
+  fitStats->Nave = 0;  
+  fitStats->Npm = 0;   
+  fitStats->Npar = 0;  
+  fitStats->Nskip = 0; 
+  fitStats->Noffset = 0;
+
+  // container to hold the fit results (Nboot > 1 for bootstrap resampling)
+  fitStats->fit = NULL;
+  fitStats->Nfit = 0;
+  fitStats->NfitAlloc = Nboot;
+  if (Nboot > 0) {
+    ALLOCATE (fitStats->fit, FitAstromResult, Nboot);
+  }
+
+  // containers to hold the measurements for a given star
+  fitStats->points = NULL;
+  fitStats->sample = NULL;
+  fitStats->Npoints = 0;
+  fitStats->NpointsAlloc = Nmax;
+  if (Nmax > 0) {
+    ALLOCATE (fitStats->points, FitAstromPoints, Nmax);
+    ALLOCATE (fitStats->sample, FitAstromPoints, Nmax);
+  }
+
+  // pre-allocated fit matrices for the 3 fit options
+  fitStats->fitPos = NULL;
+  fitStats->fitPM  = NULL;
+  fitStats->fitPar = NULL;
+
+  if (Nmax > 0) {
+    fitStats->fitPos = FitAstromDataInit (2);
+    fitStats->fitPM  = FitAstromDataInit (4);
+    fitStats->fitPar = FitAstromDataInit (5);
+  }
+
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  InitCoords (&fitStats->coords, "DEC--SIN");
+  fitStats->coords.cdelt1 = fitStats->coords.cdelt2 = 1.0 / 3600.0;
+
+  // use J2000 as a reference time
+  fitStats->T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
+  return;
+}
+
+void FitStatsReset (FitStats *tgt) {
+  tgt->Nave    = 0;  
+  tgt->Npm     = 0;   
+  tgt->Npar    = 0;  
+  tgt->Nskip   = 0; 
+  tgt->Noffset = 0;
+  return;
+}
+
+void FitStatsSum (FitStats *src, FitStats *tgt) {
+  tgt->Nave    += src->Nave    ;  
+  tgt->Npm     += src->Npm     ;   
+  tgt->Npar    += src->Npar    ;  
+  tgt->Nskip   += src->Nskip   ; 
+  tgt->Noffset += src->Noffset ;
+  return;
+}
+
+void FitStatsFree (FitStats *fitStats) {
+  if (!fitStats) return;
+
+  FREE (fitStats->fit);
+  FREE (fitStats->points);
+  FREE (fitStats->sample);
+
+  FitAstromDataFree (fitStats->fitPos);
+  FitAstromDataFree (fitStats->fitPM);
+  FitAstromDataFree (fitStats->fitPar);
+
+  free (fitStats);
+}
+
+FitAstromData *FitAstromDataInit (int Nterms) {
+
+  FitAstromData *fit = NULL;
+  ALLOCATE (fit, FitAstromData, 1);
+
+  /* do I need to do this as 2 2x2 matrix equations? */
+  fit->A = array_init (Nterms, Nterms);
+  fit->B = array_init (Nterms, 1);
+  fit->Nterms = Nterms;
+
+  return fit;
+}
+
+void FitAstromDataFree (FitAstromData *fit) {
+
+  if (!fit) return;
+
+  array_free (fit->A, fit->Nterms);
+  array_free (fit->B, fit->Nterms);
+  free (fit);
+  return;
+}
+
+}  
+
+void FitAstromObjectInit (FitAstromObject *object) {
+  object->X      = 0.0;
+  object->Y      = 0.0;
+  object->R      = 0.0;
+  object->D      = 0.0;
+  object->T      = 0.0;
+  object->dX     = 0.0;
+  object->dY     = 0.0;
+  object->dR     = 0.0;
+  object->dD     = 0.0;
+  object->dT     = 0.0;
+  object->pX     = 0.0;
+  object->pY     = 0.0;
+  object->C_blue = 0.0;
+  object->C_red  = 0.0;
+  return;
+}
+
+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);
+}  
+
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c	(revision 38599)
@@ -1,88 +1,81 @@
 # include "relastro.h"
 
-/* do we want an init function which does the alloc and a clear function to free? */
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB) {
+int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
 
   int i;
 
-  double **A, **B;
   double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
-  double chisq, Xf, Yf;
 
-  /* do I need to do this as 2 2x2 matrix equations? */
-  A = array_init (4, 4);
-  B = array_init (4, 1);
+  myAssert (data->Nterms == 4, "invalid fit arrays");
 
   Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
-  for (i = 0; i < Npts; i++) {
+
+  for (i = 0; i < Npoints; i++) {
     /* handle case where dX or dY = 0.0 */
-    wx = 1.0 / SQ(dX[i]);
-    wy = 1.0 / SQ(dY[i]);
+    wx = 1.0 / SQ(points[i].dX);
+    wy = 1.0 / SQ(points[i].dY);
 
     Wx += wx;
     Wy += wy;
 
-    Tx += T[i]*wx;
-    Ty += T[i]*wy;
+    Tx += points[i].T*wx;
+    Ty += points[i].T*wy;
     
-    Tx2 += SQ(T[i])*wx;
-    Ty2 += SQ(T[i])*wy;
+    Tx2 += SQ(points[i].T)*wx;
+    Ty2 += SQ(points[i].T)*wy;
     
-    Xs += X[i]*wx;
-    Ys += Y[i]*wy;
+    Xs += points[i].X*wx;
+    Ys += points[i].Y*wy;
 
-    XT += X[i]*T[i]*wx;
-    YT += Y[i]*T[i]*wy;
+    XT += points[i].X*points[i].T*wx;
+    YT += points[i].Y*points[i].T*wy;
   }
 
-  A[0][0] = Wx;
-  A[0][1] = Tx;
+  data->A[0][0] = Wx;
+  data->A[0][1] = Tx;
 
-  A[1][0] = Tx;
-  A[1][1] = Tx2;
+  data->A[1][0] = Tx;
+  data->A[1][1] = Tx2;
 
-  A[2][2] = Wy;
-  A[2][3] = Ty;
+  data->A[2][2] = Wy;
+  data->A[2][3] = Ty;
 
-  A[3][2] = Ty;
-  A[3][3] = Ty2;
+  data->A[3][2] = Ty;
+  data->A[3][3] = Ty2;
 
-  B[0][0] = Xs;
-  B[1][0] = XT;
-  B[2][0] = Ys;
-  B[3][0] = YT;
+  data->B[0][0] = Xs;
+  data->B[1][0] = XT;
+  data->B[2][0] = Ys;
+  data->B[3][0] = YT;
 
-  dgaussjordan (A, B, 4, 1);
+  dgaussjordan (data->A, data->B, 4, 1);
 
-  fit[0].Ro = B[0][0];
-  fit[0].uR = B[1][0];
-  fit[0].Do = B[2][0];
-  fit[0].uD = B[3][0];
-  fit[0].p  = 0.0;
+  fit->Ro = data->B[0][0];
+  fit->uR = data->B[1][0];
+  fit->Do = data->B[2][0];
+  fit->uD = data->B[3][0];
+  fit->p  = 0.0;
   
-  fit[0].dRo = sqrt(A[0][0]);
-  fit[0].duR = sqrt(A[1][1]);
-  fit[0].dDo = sqrt(A[2][2]);
-  fit[0].duD = sqrt(A[3][3]);
-  fit[0].dp  = 0.0;
+  fit->dRo = sqrt(data->A[0][0]);
+  fit->duR = sqrt(data->A[1][1]);
+  fit->dDo = sqrt(data->A[2][2]);
+  fit->duD = sqrt(data->A[3][3]);
+  fit->dp  = 0.0;
   
-  array_free (A, 4);
-  array_free (B, 4);
+  fit->Nfit = Npoints;
 
   // add up the chi square for the fit
-  chisq = 0.0;
-  for (i = 0; i < Npts; i++) {
-    Xf = fit[0].Ro + fit[0].uR*T[i];
-    Yf = fit[0].Do + fit[0].uD*T[i];
-    chisq += SQ(X[i] - Xf) / SQ(dX[i]);
-    chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
-    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+  if (fit->getChisq) {
+    double chisq = 0.0;
+    for (i = 0; i < Npoints; i++) {
+      double Xf = fit->Ro + fit->uR*points[i].T;
+      double Yf = fit->Do + fit->uD*points[i].T;
+      chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
+      chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
+    }
+    // the reduced chisq is divided by (Ndof = 2*Npts - 4)
+    fit->chisq = chisq / (2.0*Npts - 4.0);
   }
-  fit[0].Nfit = Npts;
 
-  // the reduced chisq is divided by (Ndof = 2*Npts - 4)
-  fit[0].chisq = chisq / (2.0*Npts - 4.0);
   return (TRUE);
 }
-
-// XXX this function should (optionally?) iterate and clip outlier detections
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c	(revision 38599)
@@ -1,113 +1,106 @@
 # include "relastro.h"
 
-/* do we want an init function which does the alloc and a clear function to free? */
-int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB) {
+int FitPMandPAR (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
 
   int i;
 
-  double **A, **B;
   double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
   double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
-  double chisq, Xf, Yf;
 
-  A = array_init (5, 5);
-  B = array_init (5, 1);
+  myAssert (data->Nterms == 5, "invalid fit arrays");
 
   PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
   Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
-  for (i = 0; i < Npts; i++) {
+
+  for (i = 0; i < Npoints; i++) {
     /* handle case where dX or dY = 0.0 */
-    wx = 1.0 / SQ(dX[i]);
-    wy = 1.0 / SQ(dY[i]);
+    wx = 1.0 / SQ(points[i].dX);
+    wy = 1.0 / SQ(points[i].dY);
 
     Wx += wx;
     Wy += wy;
 
-    Tx += T[i]*wx;
-    Ty += T[i]*wy;
+    Tx += points[i].T*wx;
+    Ty += points[i].T*wy;
     
-    Tx2 += SQ(T[i])*wx;
-    Ty2 += SQ(T[i])*wy;
+    Tx2 += SQ(points[i].T)*wx;
+    Ty2 += SQ(points[i].T)*wy;
     
-    PR += pR[i]*wx;
-    PD += pD[i]*wy;
+    PR += points[i].pR*wx;
+    PD += points[i].pD*wy;
     
-    PRT += pR[i]*T[i]*wx;
-    PDT += pD[i]*T[i]*wy;
+    PRT += points[i].pR*points[i].T*wx;
+    PDT += points[i].pD*points[i].T*wy;
     
-    PRX += pR[i]*X[i]*wx;
-    PDY += pD[i]*Y[i]*wy;
+    PRX += points[i].pR*points[i].X*wx;
+    PDY += points[i].pD*points[i].Y*wy;
     
-    PR2 += SQ(pR[i])*wx;
-    PD2 += SQ(pD[i])*wy;
+    PR2 += SQ(points[i].pR)*wx;
+    PD2 += SQ(points[i].pD)*wy;
 
-    Xs += X[i]*wx;
-    Ys += Y[i]*wy;
+    Xs += points[i].X*wx;
+    Ys += points[i].Y*wy;
 
-    XT += X[i]*T[i]*wx;
-    YT += Y[i]*T[i]*wy;
+    XT += points[i].X*points[i].T*wx;
+    YT += points[i].Y*points[i].T*wy;
   }
 
-  A[0][0] = Wx;
-  A[0][1] = Tx;
-  A[0][4] = PR;
+  data->A[0][0] = Wx;
+  data->A[0][1] = Tx;
+  data->A[0][4] = PR;
 
-  A[1][0] = Tx;
-  A[1][1] = Tx2;
-  A[1][4] = PRT;
+  data->A[1][0] = Tx;
+  data->A[1][1] = Tx2;
+  data->A[1][4] = PRT;
 
-  A[2][2] = Wy;
-  A[2][3] = Ty;
-  A[2][4] = PD;
+  data->A[2][2] = Wy;
+  data->A[2][3] = Ty;
+  data->A[2][4] = PD;
 
-  A[3][2] = Ty;
-  A[3][3] = Ty2;
-  A[3][4] = PDT;
+  data->A[3][2] = Ty;
+  data->A[3][3] = Ty2;
+  data->A[3][4] = PDT;
 
-  A[4][0] = PR;
-  A[4][1] = PRT;
-  A[4][2] = PD;
-  A[4][3] = PDT;
-  A[4][4] = PR2 + PD2;
+  data->A[4][0] = PR;
+  data->A[4][1] = PRT;
+  data->A[4][2] = PD;
+  data->A[4][3] = PDT;
+  data->A[4][4] = PR2 + PD2;
 
-  B[0][0] = Xs;
-  B[1][0] = XT;
-  B[2][0] = Ys;
-  B[3][0] = YT;
-  B[4][0] = PRX + PDY;
+  data->B[0][0] = Xs;
+  data->B[1][0] = XT;
+  data->B[2][0] = Ys;
+  data->B[3][0] = YT;
+  data->B[4][0] = PRX + PDY;
 
-  dgaussjordan (A, B, 5, 1);
+  dgaussjordan (data->A, data->B, 5, 1);
 
-  fit[0].Ro = B[0][0];
-  fit[0].uR = B[1][0];
-  fit[0].Do = B[2][0];
-  fit[0].uD = B[3][0];
-  fit[0].p  = B[4][0];
+  fit->Ro = data->B[0][0];
+  fit->uR = data->B[1][0];
+  fit->Do = data->B[2][0];
+  fit->uD = data->B[3][0];
+  fit->p  = data->B[4][0];
   
-  fit[0].dRo = sqrt(A[0][0]);
-  fit[0].duR = sqrt(A[1][1]);
-  fit[0].dDo = sqrt(A[2][2]);
-  fit[0].duD = sqrt(A[3][3]);
-  fit[0].dp  = sqrt(A[4][4]);
+  fit->dRo = sqrt(data->A[0][0]);
+  fit->duR = sqrt(data->A[1][1]);
+  fit->dDo = sqrt(data->A[2][2]);
+  fit->duD = sqrt(data->A[3][3]);
+  fit->dp  = sqrt(data->A[4][4]);
   
-  array_free (A, 5);
-  array_free (B, 5);
-
-  /* get the chisq from the matrix values */
+  fit->Nfit = Npoints;
 
   // add up the chi square for the fit
-  chisq = 0.0;
-  for (i = 0; i < Npts; i++) {
-    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
-    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
-    chisq += SQ(X[i] - Xf) / SQ(dX[i]);
-    chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
-    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+  if (fit->getChisq) {
+    double chisq = 0.0;
+    for (i = 0; i < Npoints; i++) {
+      double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR;
+      double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD;
+      chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
+      chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
+    }
+    // the reduced chisq is divided by (Ndof = 2*Npts - 5)
+    fit->chisq = chisq / (2.0*Npts - 5.0);
+  }
 
-  }
-  fit[0].Nfit = Npts;
-
-  // the reduced chisq is divided by (Ndof = 2*Npts - 5)
-  fit[0].chisq = chisq / (2.0*Npts - 5.0);
   return (TRUE);
 }
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c	(revision 38599)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c	(revision 38599)
@@ -0,0 +1,53 @@
+# include "relastro.h"
+
+int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
+
+  int i;
+
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+
+  myAssert (data->Nterms == 2, "invalid fit arrays");
+
+  Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
+
+  for (i = 0; i < Npoints; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(points[i].dX);
+    wy = 1.0 / SQ(points[i].dY);
+
+    Wx += wx;
+    Wy += wy;
+
+    Tx += points[i].T*wx;
+    Ty += points[i].T*wy;
+    
+    Xs += points[i].X*wx;
+    Ys += points[i].Y*wy;
+  }
+
+  fit->Ro = (Xs - fit->uR*Tx) / Wx;
+  fit->Do = (Ys - fit->uD*Ty) / Wy;
+
+  fit->p  = 0.0;
+  
+  // should be ~ 1.0 / Wx
+  fit->dRo = sqrt(1.0 / Wx);
+  fit->dDo = sqrt(1.0 / Wy);
+  
+  fit->Nfit = Npoints;
+
+  // add up the chi square for the fit
+  if (fit->getChisq) {
+    double chisq = 0.0;
+    for (i = 0; i < Npoints; i++) {
+      double Xf = fit->Ro + fit->uR*points[i].T;
+      double Yf = fit->Do + fit->uD*points[i].T;
+      chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
+      chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
+    }
+    // the reduced chisq is divided by (Ndof = 2*Npts - 4)
+    fit->chisq = chisq / (2.0*Npts - 4.0);
+  }
+
+  return (TRUE);
+}
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/GetAstromError.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/GetAstromError.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/GetAstromError.c	(revision 38599)
@@ -1,4 +1,8 @@
 # include "relastro.h"
 # define WEIGHTED_ERRORS 1
+
+// XXX hard-wire the trends identified by CZW
+static float BrightMo[] = {-15.6, -16.8, -17.0, -16.7, -16.0};
+static float BrightMs[] = {1.3, 1.3, 1.3, 1.8, 2.0}; 
 
 float GetAstromErrorTiny (MeasureTiny *measure, int mode) {
@@ -40,4 +44,15 @@
   dM    = measure[0].dM;
   dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
+
+  // for GPC1 data, we have a bright end model:
+  if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) {
+    int Np = ((int) (measure[0].photcode / 100)) % 100;
+    myAssert (Np >= 0, "oops");
+    myAssert (Np <= 4, "oops");
+
+    float Minst = measure[0].M - measure[0].dt - 25.0;
+    float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np])));
+    dPtotal = hypot(dPtotal, dPbright);
+  }
 
   dPtotal = MAX (dPtotal, MIN_ERROR);
@@ -84,4 +99,15 @@
   dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
 
+  // for GPC1 data, we have a bright end model:
+  if ((measure[0].photcode > 10000) && (measure[0].photcode < 10480)) {
+    int Np = ((int) (measure[0].photcode / 100)) % 100;
+    myAssert (Np >= 0, "oops");
+    myAssert (Np <= 4, "oops");
+
+    float Minst = measure[0].M - measure[0].dt - 25.0;
+    float dPbright = 0.335 / (1.0 + exp(BrightMs[Np]*(Minst - BrightMo[Np])));
+    dPtotal = hypot(dPtotal, dPbright);
+  }
+
   dPtotal = MAX (dPtotal, MIN_ERROR);
   return (dPtotal);
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ImageOps.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ImageOps.c	(revision 38599)
@@ -748,7 +748,14 @@
     n = measure[0].averef;
 
+    // apply proper-motion from average position to measure epoch:
+    float dTime = (measure[k].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years
+
     /* apply the current image transformation or use the current value of R+dR, D+dD? */
     ref[i].R = catalog[c].average[n].R;
     ref[i].D = catalog[c].average[n].D;
+    
+    // XXX do this in a better way?
+    ref[i].R += dTime * catalog[c].average[n].uR / 3600.0 / cos(ref[i].D*RAD_DEG);
+    ref[i].D += dTime * catalog[c].average[n].uD / 3600.0;
 
     // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ParFactor.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ParFactor.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/ParFactor.c	(revision 38599)
@@ -2,49 +2,9 @@
 # define J2000 2451545.       /* Julian date at standard epoch */
 
-# if (0)
-/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
-/* ra and dec are returned as decimal hours and decimal degrees. */
-void lpsun (double jd, double *ra, double *dec) {
+/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
+// jdoff is days since J2000
+int sun_ecliptic (double jdoff, double *lambda, double *beta, double *epsilon, double *Radius) {
 
-  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
-
-  n = jd - J2000;
-  L = 280.460 + 0.9856474 * n;
-  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
-  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
-  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
-
-  // this is the conversion from ecliptic to celestial coords
-  x = cos(lambda);
-  y = cos(epsilon)*sin(lambda);
-  z = sin(epsilon)*sin(lambda);
-
-  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
-  *dec = (asin(z))*DEG_IN_RADIAN;
-}
-# endif
-
-# if (0)
-/* code borrowed from Skycalc : fix this stuff XXX */
-/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
-int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
-
-  double n, L, g;
-
-
-  n = jd - J2000;
-  L = 280.460 + 0.9856474 * n;
-  g = (357.528 + 0.9856003 * n)*RAD_DEG;
-  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
-  *beta = 0.0;					  // approx latitude
-  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
-  return TRUE;
-}
-# endif
-
-/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
-int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) {
-
-  double n = jd - J2000;	      // day number
+  double n = jdoff;	      // day number
   double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
   double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly
@@ -58,17 +18,14 @@
 
 /* given RA, DEC, Time, calculate the parallax factor */
-// Time is relative to Tmean, Tmean is years relative to J2000
-int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) {
+// Time is years since J2000
+int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) {
 
-  double jd, lambda, beta, epsilon, Radius;
+  double lambda, beta, epsilon, Radius;
 
-  /* given a Time relative to Tmean, Tmean in years since J2000, determine the solar
-    longitude S */
+  /* given a Time in years since J2000, determine the solar longitude S */
 
-  // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
-  jd = 365.25*(Time + Tmean) + J2000;
-  // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd);
+  double jdoff = 365.25*Time;
 
-  sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius);
+  sun_ecliptic (jdoff, &lambda, &beta, &epsilon, &Radius);
 
   double lambda_rad = lambda*RAD_DEG;
@@ -99,2 +56,42 @@
   return TRUE;
 }
+
+# if (0)
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+/* ra and dec are returned as decimal hours and decimal degrees. */
+void lpsun (double jd, double *ra, double *dec) {
+
+  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
+  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
+  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
+
+  // this is the conversion from ecliptic to celestial coords
+  x = cos(lambda);
+  y = cos(epsilon)*sin(lambda);
+  z = sin(epsilon)*sin(lambda);
+
+  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
+  *dec = (asin(z))*DEG_IN_RADIAN;
+}
+
+/* code borrowed from Skycalc : fix this stuff XXX */
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
+
+  double n, L, g;
+
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)*RAD_DEG;
+  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
+  *beta = 0.0;					  // approx latitude
+  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
+  return TRUE;
+}
+# endif
+
Index: /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c	(revision 38598)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c	(revision 38599)
@@ -1,105 +1,12 @@
 # include "relastro.h"
 # define PAR_TOOFEW 5
+# define NBOOT 100
 
 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
 int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
-
-static off_t Nmax;
-static double *X, *dX;
-static double *Y, *dY;
-static double *R, *dR;
-static double *D, *dD;
-static double *pX;
-static double *pY;
-static double *T;
-static double *dT;
-static double *C_blue;
-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) {
-
-  off_t i, j;
-  
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
-    }
-  }
-
-  ALLOCATE (R, double, MAX (1, Nmax));
-  ALLOCATE (D, double, MAX (1, Nmax));
-  ALLOCATE (T, double, MAX (1, Nmax));
-  ALLOCATE (X, double, MAX (1, Nmax));
-  ALLOCATE (Y, double, MAX (1, Nmax));
-
-  ALLOCATE (dR, double, MAX (1, Nmax));
-  ALLOCATE (dD, double, MAX (1, Nmax));
-  ALLOCATE (dT, double, MAX (1, Nmax));
-  ALLOCATE (dX, double, MAX (1, Nmax));
-  ALLOCATE (dY, double, MAX (1, Nmax));
-
-  ALLOCATE (pX, double, MAX (1, Nmax));
-  ALLOCATE (pY, double, MAX (1, Nmax));
-
-  ALLOCATE (C_blue, double, MAX (1, Nmax));
-  ALLOCATE (C_red,  double, MAX (1, Nmax));
-
-  /* project coordinates to a plane centered on the object with units of arcsec */
-  InitCoords (&coords, "DEC--SIN");
-  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
-
-  // use J2000 as a reference time
-  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, int Nloop) {
-
-  initObjectData (catalog, Ncatalog);
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -109,6 +16,12 @@
   }
 
-  FitStats sumStatsChips; initFitStats (&sumStatsChips);
-  FitStats sumStatsStack; initFitStats (&sumStatsStack);
+  int NmeasureMax = CatalogMaxNmeasure (catalog, Ncatalog);
+
+  // allocate summary stats with Nmax = 0, Nboot = 0
+  FitStats *sumStatsChips = FitStatsInit (0, 0);
+  FitStats *sumStatsStack = FitStatsInit (0, 0);
+
+  FitStats *fitStatsChips = FitStatsInit (NmeasureMax, NBOOT);
+  FitStats *fitStatsStack = FitStatsInit (NmeasureMax, NBOOT);
 
   int i;
@@ -117,6 +30,6 @@
     if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
 
-    FitStats fitStatsChips; initFitStats (&fitStatsChips);
-    FitStats fitStatsStack; initFitStats (&fitStatsStack);
+    FitStatsReset (fitStatsChips);
+    FitStatsReset (fitStatsStack);
 
     off_t j;
@@ -129,16 +42,22 @@
       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, Nloop);
+      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack);
+      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m, Nloop);
     }
     if (VERBOSE2) 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 (VERBOSE2) 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 ();
+    FitStatsSum (fitStatsChips, sumStatsChips);
+    FitStatsSum (fitStatsStack, sumStatsStack);
+  }
+
+  FitStatsFree (fitStatsChips);
+  FitStatsFree (fitStatsStack);
 
   if (VERBOSE && (Ncatalog > 1)) 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 && (Ncatalog > 1)) 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);
+
+  FitStatsFree (sumStatsChips);
+  FitStatsFree (sumStatsStack);
+
   return (TRUE);
 }
@@ -146,5 +65,5 @@
 // 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 applyGalaxyOffset) {
+int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff) {
 
   int setRefColor = areImagesMatched();
@@ -154,5 +73,5 @@
   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 fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq  = NAN;
   PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;
 
@@ -169,5 +88,4 @@
   int NcBlue = 0;
   int NcRed = 0;
-  int N = 0;
 
   int mode = FIT_MODE; // start with the globally-defined fit mode
@@ -177,123 +95,9 @@
   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);
-
-    // if we are correcting for the Galaxy Motion Model, only should apply it here 
-    // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and
-    // (b) if we think the measure R,D is already at the image epoch position
-    if (USE_GALAXY_MODEL && applyGalaxyOffset) {
-      Ri -= measure[k].RoffGAL / 3600.0;
-      Di -= measure[k].DoffGAL / 3600.0;
-    }
-
-    // XXX add in dR,dD GAL here
-
-    // 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;
-      }
-      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++;
-      }
-      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 : average[0].Nmeasure 
-
-  if (N < 1) { 
+  // select the measurements to be used in this analysis
+  FitAstromObject *points = UpdateObjects_SelectMeasures (average, secfilt, measure, measureBig, cat, measOff, &Npoints);
+
+  // if there are no exposure detections, use the stack position
+  if (Npoints < 1) { 
     if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
       average[0].R  = average[0].Rstk;
@@ -301,44 +105,16 @@
       average[0].dR = average[0].dRstk;
       average[0].dD = average[0].dDstk;
+      average[0].flags |= ID_STACK_ASTROM;
     }
     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];
+  double Tmean, Trange, parRange;
+  UpdateObjects_Project (points, Npoints, &Tmean, &Trange, &parRange);
 
   // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
 
-  // 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]);
-    }
-  }	  
+  // if we have too few good detections for the desired fit, or too limited a baseline,
+  // use a fit with fewer parameters.
 
   // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
@@ -346,12 +122,18 @@
     if (Trange < PM_DT_MIN) {
       mode = FIT_AVERAGE;
-      goto skipPM;
-    }
-    if (N <= PM_TOOFEW) {
+      goto justPosition;
+    }
+    if (Npoints <= PM_TOOFEW) {
       mode = FIT_AVERAGE;
-      goto skipPM;
-    }
-
-    FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
+      goto justPosition;
+    }
+
+    fitStats->Nfit = 0;
+    for (k = 0; k < fitStats->NfitAlloc; k++) {
+      BootstrapResample (sample, points, Npoints);
+      if (!FitPM (&fitStats->fit[k], fitdata, sample, Npoints)) continue;
+      fitStats->Nfit ++;
+    }
+    // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
 
     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);
@@ -374,38 +156,25 @@
   }
   
-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);
-	
+      goto justPosition;
+    }
     if (parRange < PAR_FACTOR_MIN) {
       mode = FIT_PM_ONLY;
-      goto skipPAR;
-    }
-
-    FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
+      goto justPosition;
+    }
+
+    fitStats->Nfit = 0;
+    for (k = 0; k < fitStats->NfitAlloc; k++) {
+      BootstrapResample (sample, points, Npoints);
+      FitPMandPar (&fitStats->fit[k], fitsdata, sample, Npoints);
+      fitStats->Nfit ++;
+    }
+    // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
+
     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);
 
@@ -414,5 +183,5 @@
     fitStats->Npar ++;
 
-    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
+    if (fabs(fitPAR.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
 
     // XXX a hard-wired hack...
@@ -422,29 +191,43 @@
   }	  
 
-skipPAR:
+justPosition:
   {
-    // ALWAYS fit the average model
-    StatType statsR, statsD;
-    liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
-    liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
+    // use bootstrap resampling to check the error distribution
+    // if we only have one point, this is silly...
+    
+    fitStats->Nfit = 0;
+    for (k = 0; k < fitStats->NfitAlloc; k++) {
+      BootstrapResample (sample, points, Npoints);
+      FitPosPMfixed (&fitStats->fit[k], fitsdata, sample, Npoints);
+      fitStats->Nfit ++;
+    }
+
+    if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
 
     // project Ro, Do back to RA,DEC
-    XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords);
-    if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
-
-    fitAve.dRo = statsR.sigma;
-    fitAve.dDo = 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;
+    XY_to_RD (&fitAve.Ro, &fitAve.Do, fitAve.Ro, fitAve.Do, &coords);
     average[0].flags |= ID_STAR_FIT_AVE;
     fitStats->Nave ++;
+
+    // XXX fitAve.dRo = statsR.sigma;
+    // XXX fitAve.dDo = statsD.sigma;
+    // XXX 
+    // XXX fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
+    // XXX fitAve.Nfit = N;
+    // XXX 
+    // XXX fitAve.uR = fitAve.duR = 0.0;
+    // XXX fitAve.uD = fitAve.duD = 0.0;
+    // XXX fitAve.p  = fitAve.dp  = 0.0;
+  }
+
+  // XXX update the bit flags of which points were used
+  for (k = 0; k < Npoints; k++) {
+    int Nm = points[k].measure;
+    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
+    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
   }
 
   if (setRefColor) {
+    // need to reassign here if isfinite()
     float colorMedian;
     dsort (C_blue, NcBlue);
@@ -696,5 +479,182 @@
 }
 
-
+int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) {
+  int i;
+
+  // I need to draw Npoints random entries from 'points' with replacement:
+  for (i = 0; i < Npoints; i++) {
+    int N = Npoints * drand48();
+    sample[i] = points[N];
+  }
+  return TRUE;
+}
+
+FitAstromObject *UpdateObjects_SelectMeasures (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int *npoints) {
+
+  // XXX move this above, but just for re-working now:
+  int Npoints = 0;
+  int NPOINTS = average->Nmeasure;
+  FitAstromObject *points = NULL;
+  ALLOCATE (points, FitAstromObject, NPOINTS);
+
+  // 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 = measure[k].R;
+    double Di = measure[k].D;
+
+    // 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;
+      }
+      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;
+    }
+
+    FitAstromObjectInit (&points[Npoints]);
+
+    points[Npoints].R = Ri;
+    points[Npoints].D = Di;
+
+    // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
+    // T[] is time in years since J2000 (jd = 2451545)
+    points[Npoints].T = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
+
+    // add measured systematic error in quadrature?  only do this after the fit has
+    // converged (or you will never improve the poor images)
+
+    // 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
+
+    // dX, dY : error in arcsec:
+    points[Npoints].dX = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
+    points[Npoints].dY = 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(points[Npoints].dX)) continue;
+    if (isnan(points[Npoints].dY)) continue;
+
+    points[Npoints].dT = measure[k].dt;
+
+    if (setRefColor) {
+      points[Npoints].C_blue = getColorBlue (measOff+k, cat);
+      points[Npoints].C_red  = getColorRed (measOff+k, cat);
+    }
+
+    points[Npoints].measure = k;
+    Npoints++;
+  } // loop over measurements : average[0].Nmeasure 
+
+  *npoints = Npoints;
+  return points;
+}
+
+int FitAstromObject_Project (FitAstromObject *points, int Npoints, double *Tmean, double *Trange, double *parRange) {
+
+  int k;
+
+  // find Tmin & Tmax from the list of accepted measurements
+  double Tmin  = points[0].T;
+  double Tmax  = points[0].T;
+  double pXmin = +2.0;
+  double pXmax = -2.0;
+  double pYmin = +2.0;
+  double pYmax = -2.0;
+
+  *Tmean = 0.0;
+  for (k = 0; k < Npoints; k++) {
+    Tmin = MIN(Tmin, points[k].T);
+    Tmax = MAX(Tmax, points[k].T);
+    Tmean += points[k].T;
+
+    // at this point, T is in years since J2000
+    ParFactor (&points[k].pX, &points[k].pY, points[k].R, points[k].D, points[k].T);
+    pXmin = MIN (pXmin, points[k].pX);
+    pXmax = MAX (pXmax, points[k].pX);
+    pYmin = MIN (pYmin, points[k].pY);
+    pYmax = MAX (pYmax, points[k].pY);
+  }
+  *Trange = Tmax - Tmin;
+
+  // mean epoch
+  *Tmean /= (float) Npoints;
+
+  // for HIGH_SPEED, just use the center of the range
+  if (RELASTRO_OP == OP_HIGH_SPEED) {
+    *Tmean = 0.5*(Tmax - Tmin);
+  }
+
+  double dXRange = pXmax - pXmin;
+  double dYRange = pYmax - pYmin;
+  *parRange = hypot (dXRange, dYRange);
+
+  /* we need to do the fit in a locally linear space; choose a ref coordinate */
+  coords.crval1 = points[0].R;
+  coords.crval2 = points[0].D;
+
+  // project all of the R,D coordinates to a plane centered on this coordinate. set
+  // the times to be relative to Tmean
+  for (k = 0; k < Npoints; k++) {
+    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &coords);
+    points[k].T -= Tmean;
+  }	  
+  return TRUE;
+}
+
+int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
+
+  int i, j;
+
+  Nmax = 0;
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
+    }
+  }
+
+  return Nmax;
+}
 
 /* fitting proper-motion and parallax: