Changeset 38986 for trunk/Ohana/src/relastro

trunk/Ohana

Property svn:mergeinfo set to
/branches/eam_branches/ipp-20150625/Ohana merged eligible

trunk/Ohana/src/relastro/Makefile

-              r38062
+              r38986
 relastro_client: $(BIN)/relastro_client.$(ARCH)
+install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client $(DESTBIN)/testparallax
+install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client
+# $(DESTBIN)/testparallax
 RELASTRO = \
 …
 $(SRC)/FitChip.$(ARCH).o             \
 $(SRC)/FitMosaic.$(ARCH).o           \
+$(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/FitAstromOps.$(ARCH).o           \
 $(SRC)/FitPM.$(ARCH).o               \
-$(SRC)/FitPar.$(ARCH).o              \
 $(SRC)/FitPMandPar.$(ARCH).o         \
 $(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/FitPosPMfixed.$(ARCH).o       \
 $(SRC)/ImageOps.$(ARCH).o            \
 $(SRC)/MosaicOps.$(ARCH).o           \
 …
 RELASTRO_CLIENT = \
 $(SRC)/ConfigInit.$(ARCH).o          \
+$(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/FitAstromOps.$(ARCH).o           \
 $(SRC)/FitPM.$(ARCH).o               \
-$(SRC)/FitPar.$(ARCH).o              \
 $(SRC)/FitPMandPar.$(ARCH).o         \
 $(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/FitPosPMfixed.$(ARCH).o       \
 $(SRC)/ImageOps.$(ARCH).o            \
 $(SRC)/MosaicOps.$(ARCH).o           \

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

-              r38441
+              r38986
 typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
+typedef enum {
+  FIT_RESULT_RA,
+  FIT_RESULT_DEC,
+  FIT_RESULT_uR,
+  FIT_RESULT_uD,
+  FIT_RESULT_PLX,
+} FitAstromResultMode;
 typedef enum {
 …
   double Ro, dRo;
   double Do, dDo;
   double uR, duR;
   double uD, duD;
+  double p, dp;
+  double  p, dp;
   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 pR;
+  double pD;
+  double C_blue;
+  double C_red;
+  int measure;
+} FitAstromPoint;
+typedef struct {
+  off_t Nave;
+  off_t Npm;
+  off_t Npar;
+  off_t Nskip;
+  off_t Noffset;
+  double *values;
+  FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats
+  int Nfit;
+  int NfitAlloc;
+  FitAstromPoint *points;
+  FitAstromPoint *sample;
+  int Npoints;
+  int NpointsAlloc;
+  FitAstromData *fitdataPos;
+  FitAstromData *fitdataPM;
+  FitAstromData *fitdataPar;
+  Coords coords;
+  time_t T2000;
+} FitStats;
 typedef struct {
 …
 typedef struct {
-  off_t Nave;
-  off_t Npm;
-  off_t Npar;
-  off_t Nskip;
-  off_t Noffset;
-} FitStats;
-typedef struct {
   double scale;
   double **Roff;
 …
   int Nicrfobj;
 } ICRFobj;
+# define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM
 /* global variables set in parameter file */
 …
 int    PARALLEL_OUTPUT;
+int    VERBOSE_IMAGE;
 int    VERBOSE;
 int    VERBOSE2;
 …
 int    USE_ICRF_SHFIT;
 int    USE_ICRF_POLE;
+int    USE_ALL_IMAGES;
 int    RESET;
 …
 int    CHIPORDER;
 int    CHIPMAP;
+int    N_BOOTSTRAP_SAMPLES;
 int MaxDensityUse;
 …
 int           liststats_pos       PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB));
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile));
 int           load_images         PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap));
+int           load_images         PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap, int UseAllImages));
 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap));
 …
 int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop);
 int UpdateSimple (Catalog *catalog, int Ncatalog);
 int UpdateChips (Catalog *catalog, int Ncatalog);
+int UpdateChips (Catalog *catalog, int Ncatalog, int Nloop);
 int UpdateMosaic (Catalog *catalog, int Ncatalog);
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
 …
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius);
 int ParFactor (double *pR, double *pD, double RA, double Dec, double Time, double Tmean);
 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB);
 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB);
+int ParFactor (double *pR, double *pD, double RA, double Dec, double Time);
+int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
+int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
 Mosaic *getMosaicForImage (off_t N);
 …
 int client_logger_init (char *dirname);
 int client_logger_message (char *format,...);
+int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode);
+double VectorFractionInterpolate (double *values, float fraction, int Npts);
+int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode);
+int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints);
+int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
+int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);
+int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack);
+int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
+int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff);
+void FitAstromResultInit (FitAstromResult *fit);
+void FitAstromPointInit (FitAstromPoint *object);
+void FitAstromDataFree (FitAstromData *fit);
+FitAstromData *FitAstromDataInit (int Nterms);
+void FitStatsFree (FitStats *fitStats);
+void FitStatsSum (FitStats *src, FitStats *tgt);
+void FitStatsReset (FitStats *tgt);
+FitStats *FitStatsInit (int Nmax, int Nboot);
+int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average);
+void AstromErrorSetLoop (int Nloop, int isImageMode);

trunk/Ohana/src/relastro/src/FitChip.c

-              r37807
+              r38986
+    }
+    if (VERBOSE2) fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
+    if (VERBOSE2) {
+      fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
+    }
     // when fitting the map, first fit a linear model (below? change Npolyterms to -1)
 …
     for (i = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
       fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
+      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, ref[i].dPos);
+    }
 …
   int i, N;
   float *x, *y, *dX, *dY;
+  float *x, *y, *dX, *dY, *dP;
   ALLOCATE (x,  float, Npts);
   ALLOCATE (y,  float, Npts);
   ALLOCATE (dX, float, Npts);
   ALLOCATE (dY, float, Npts);
+  ALLOCATE (dP, float, Npts);
   N = 0;
 …
     dX[N] = ref[i].X - raw[i].X;
     dY[N] = ref[i].Y - raw[i].Y;
+    dP[N] = ref[i].dPos;
     N++;
+  }
 …
   // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY)
   AstromOffsetMapFit (map, x, y, dX, N, TRUE);
   AstromOffsetMapFit (map, x, y, dY, N, FALSE);
+  AstromOffsetMapFit (map, x, y, dX, dP, N, TRUE);
+  AstromOffsetMapFit (map, x, y, dY, dP, N, FALSE);
   AstromOffsetMapRepair (map, TRUE);
 …
   free (dX);
   free (dY);
+  free (dP);
   return TRUE;

trunk/Ohana/src/relastro/src/FitPM.c

-              r32695
+              r38986
 # 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) {
+// initial values of *fit are ignored
+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);
+  }
-  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

trunk/Ohana/src/relastro/src/FitPMandPar.c

-              r32695
+              r38986
 # 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) {
+// initial values of *fit are ignored
+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);
+  fit->Nfit = Npoints;
-  /* get the chisq from the matrix values */
-  // 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);
+  }
-  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);
+}

trunk/Ohana/src/relastro/src/FrameCorrection.c

-              r38062
+              r38986
+    }
     AstromOffsetMapFit (map, Xfit, Yfit, dXfit, Nkeep, TRUE);
     AstromOffsetMapFit (map, Xfit, Yfit, dYfit, Nkeep, FALSE);
+    AstromOffsetMapFit (map, Xfit, Yfit, dXfit, NULL, Nkeep, TRUE);
+    AstromOffsetMapFit (map, Xfit, Yfit, dYfit, NULL, Nkeep, FALSE);
     for (i = 0; i < Npts; i++) {

trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c

-              r38553
+              r38986
   float *buffer = (float *)matrix->buffer;
   float **value = (raDirection) ? map->dXv : map->dYv;
+  float *value = (raDirection) ? map->dXv : map->dYv;
   int ix, iy;
   for (ix = 0; ix < map->Nx; ix++) {
     for (iy = 0; iy < map->Ny; iy++) {
       buffer[ix + map->Nx*iy] = value[ix][iy];
+      buffer[ix + map->Nx*iy] = value[ix + map->Nx*iy];
+    }
+  }
 …
   float *buffer = (float *) matrix->buffer;
   float **value = (raDirection) ? map->dXv : map->dYv;
+  float *value = (raDirection) ? map->dXv : map->dYv;
   int ix, iy;
   for (ix = 0; ix < map->Nx; ix++) {
     for (iy = 0; iy < map->Ny; iy++) {
       value[ix][iy] = buffer[ix + map->Nx*iy];
+      value[ix + map->Nx*iy] = buffer[ix + map->Nx*iy];
+    }
+  }

trunk/Ohana/src/relastro/src/GetAstromError.c

-              r36833
+              r38986
 # 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};
+static int Nloop = -1;
+static int isImage = FALSE;
+// Nloop is used to modify the per detection errors
+void AstromErrorSetLoop (int N, int isImageMode) {
+  Nloop = N;
+  isImage = isImageMode;
+}
 float GetAstromErrorTiny (MeasureTiny *measure, int mode) {
 …
   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);
+  // early on, we want 2MASS and Tycho to have a very high weight.  This will force images
+  // to match the 2MASS / Tycho / ICRS reference frame.  As Nloop gets higher, the weight
+  // needs to drop to allow the ps1 measurements to drive the solution
+  int is2MASS  = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2011) && (measure[0].photcode <= 2013);
+  int isTycho  = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2020) && (measure[0].photcode <= 2021);
+  int has2MASS = USE_GALAXY_MODEL &&  isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS);
+  // serious hack: if the object has 2MASS, we set this internal bit and adjust the
+  // weight to ensure the image is tied down to the 2mass frame
+  if (has2MASS) {
+    switch (Nloop) {
+      case 0:
+        dPtotal = dPtotal / 1000.0;
+        break;
+      case 1:
+        dPtotal = dPtotal / 300.0;
+        break;
+      case 2:
+        dPtotal = dPtotal / 100.0;
+        break;
+      case 3:
+        dPtotal = dPtotal / 30.0;
+        break;
+      case 4:
+        dPtotal = dPtotal / 10.0;
+        break;
+      case 5:
+        dPtotal = dPtotal / 10.0;
+        break;
+      default:
+        break;
+    }
+  }
+  if (is2MASS) {
+    switch (Nloop) {
+      case 0:
+        dPtotal = dPtotal / 1000.0;
+        break;
+      case 1:
+        dPtotal = dPtotal / 300.0;
+        break;
+      case 2:
+        dPtotal = dPtotal / 100.0;
+        break;
+      case 3:
+        dPtotal = dPtotal / 30.0;
+        break;
+      case 4:
+        dPtotal = dPtotal / 10.0;
+        break;
+      case 5:
+        dPtotal = dPtotal / 10.0;
+        break;
+      default:
+        break;
+    }
+  }
+  if (isTycho) {
+    switch (Nloop) {
+      case 0:
+      case 1:
+        dPtotal = dPtotal / 200.0;
+        break;
+      case 2:
+      case 3:
+        dPtotal = dPtotal / 100.0;
+        break;
+      case 4:
+      case 5:
+        dPtotal = dPtotal / 50.0;
+        break;
+      default:
+        break;
+    }
+  }
   return (dPtotal);
+}
 …
   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);

trunk/Ohana/src/relastro/src/ImageOps.c

-              r38062
+              r38986
 # include "relastro.h"
+int isGPC1chip (int photcode);
 # define USE_IMAGE_ID 1
 …
       raw[i].mask |= MARK_NAN_POS_ERROR;
+    }
+    // XXX A TEST: can we use only 2MASS measurements to fit the images?
+    // XXX Do NOT apply this for the real calibration
+    if (FALSE && !(measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS)) {
+      fprintf (stderr, "@");
+      raw[i].mask |= MARK_BIG_OFFSET;
+    }
     raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is
 …
     ref[i].R = catalog[c].average[n].R;
     ref[i].D = catalog[c].average[n].D;
+    // if we are applying the galaxy model, move the reference position...
+    if (USE_GALAXY_MODEL) {
+      // apply proper-motion from average position to measure epoch:
+      float dTime = (measure[0].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years
+      // 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
     if (USE_GALAXY_MODEL) {
+    if (0) {
       myAssert (!isnan(measure[0].RoffGAL), "oops");
       myAssert (!isnan(measure[0].DoffGAL), "oops");
 …
   if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
   free (R);
   free(dR);
   free(D);
   free(dD);
+  free (dR);
+  free (D);
+  free (dD);
+}
 …
   if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
   free (R);
   free(dR);
   free(D);
   free(dD);
   free(d2);
   free(index);
+  free (dR);
+  free (D);
+  free (dD);
+  free (d2);
+  free (index);
+}
 …
+  }
+  if (MinBadQF > 0.0) {
+  if ((MinBadQF > 0.0) && isGPC1chip(measure[0].photcode)) {
+    if (!isfinite(measure[0].psfQF)) return FALSE;
     if (measure[0].psfQF < MinBadQF) return FALSE;
+  }

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

-              r37261
+              r38986
 # 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
 …
 /* 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;
 …
   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

trunk/Ohana/src/relastro/src/Shutdown.c

r37038	r38986
21	21	va_end (argp);
22	22
23		~~SetProtect (TRUE);~~
24		~~gfits_db_close (db);~~
25	23	fprintf (stderr, "ERROR: relastro halted\n");
26	24	exit (1);

trunk/Ohana/src/relastro/src/StarMaps.c

r37807	r38986
163	163	}
164	164
165		if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts ~~: %f, %f\n", images[N].name, starmap[N].Npoints~~, dLmax, dMmax);
	165	if (VERBOSE2) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax);
166	166
167	167	if (dLmax > DPOS_MAX) {
168		~~if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints~~, dLmax, dMmax);
	168	if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax);
169	169	return (FALSE);
170	170	}
171	171	if (dMmax > DPOS_MAX) {
172		if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts ~~: %f, %f\n", images[N].name, starmap[N].Npoints~~, dLmax, dMmax);
	172	if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts (%d fitted) : %f, %f\n", images[N].name, starmap[N].Npoints, images[N].nFitAstrom, dLmax, dMmax);
173	173	return (FALSE);
174	174	}

trunk/Ohana/src/relastro/src/UpdateChips.c

-              r38062
+              r38986
 // update astrometry of all chips relative to the average positions
 // if NTHREADS is non-zero, call the threaded version of this function
 int UpdateChips (Catalog *catalog, int Ncatalog) {
+int UpdateChips (Catalog *catalog, int Ncatalog, int Nloop) {
   off_t Nskip, Nmosaic, NnewFit, NoldFit;
 …
   char *mode;
+  AstromErrorSetLoop (Nloop, TRUE);
   if (NTHREADS) {
     UpdateChips_threaded (catalog, Ncatalog);
 …
   // each chip is fitted independently, so we could do N at once in parallel
   for (i = 0; i < Nimage; i++) {
+    VERBOSE_IMAGE = !strcmp(image[i].name, "o5745g0516o.356887.cm.982631.smf[XY54]");
     // XXX looks like everything below is thread safe : we can unroll this into a set of
 …
+    }
+    if (VERBOSE_IMAGE) {
+      dump_stardata_pts (raw, Nraw, "testimage.raw.dat");
+      dump_stardata_pts (ref, Nref, "testimage.ref.dat");
+      fprintf (stderr, "dumped test image\n");
+    }
     // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
     assert (Nraw == Nref);
 …
     // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
     if (!FitChip (raw, ref, Nraw, &image[i])) {
       if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused %d of %d\n", image[i].name,  i,  Nraw, image[i].nFitAstrom, image[i].nstar);
       if (1) {
 …
     if (!checkStarMap (i)) {
       if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused: %d\n", image[i].name,  i,  Nraw, image[i].nFitAstrom);
       if (1) {
 …
     setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC);
     if (USE_GALAXY_MODEL) {
+      // XXX DEPRECATE?
       image[i].flags |= ID_IMAGE_ASTROM_GMM;
+    }
 …
     // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
     if (!FitChip (raw, ref, Nraw, &image[i])) {
       if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused %d of %d\n", image[i].name,  i,  Nraw, image[i].nFitAstrom, image[i].nstar);
       // restore status quo ante (replace truMap with tmpMap)
 …
     if (!checkStarMap (i)) {
       if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused: %d\n", image[i].name,  i,  Nraw, image[i].nFitAstrom);
       // restore status quo ante (replace truMap with tmpMap)

trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c

-              r38471
+              r38986
+  }
 # if (0)
+# if (1)
   UpdateObjectOffsets_parallel_table (table, sky);
 …
     char *command = NULL;
+    strextend (&command, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f",
+              group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR);
+    strextend (&command, "relastro_client -update-offsets");
+    strextend (&command, "-hostID %d", group->hosts[i][0].hostID);
+    strextend (&command, "-D CATDIR %s", CATDIR);
+    strextend (&command, "-hostdir %s", group->hosts[i][0].pathname);
+    strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    strextend (&command, "-statmode %s", STATMODE);
+    strextend (&command, "-minerror %f", MIN_ERROR);
     if (FIT_MODE == FIT_PM_ONLY)         strextend (&command, "-pm");
 …
     if (ExcludeBogus)    strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius);
+    if (USE_ALL_IMAGES)      strextend (&command, "-use-all-images");
     if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1");
 …
     if (PhotFlagPoor)       strextend (&command, "+photflagpoor %d", PhotFlagPoor);
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) {
 …
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) {
       strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG);

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

-              r37807
+              r38986
 # 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 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);
+}
+int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog);
 // 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?
 …
+  }
+  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, N_BOOTSTRAP_SAMPLES);
+  FitStats *fitStatsStack = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES);
+  AstromErrorSetLoop (Nloop, FALSE);
   int i;
 …
     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;
 …
       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);
+    }
+    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 ();
+  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);
+      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack);
+      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m);
+    }
+    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);
+    FitStatsSum (fitStatsChips, sumStatsChips);
+    FitStatsSum (fitStatsStack, sumStatsStack);
+  }
+  // DumpObjectsWith2MASS (catalog, Ncatalog);
+  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);
+}
 …
 // 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 setRefColor = areImagesMatched();
+int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog) {
+  int i;
+  for (i = 0; i < Ncatalog; i++) {
+    off_t j;
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      /* calculate the average value of R,D for a single star */
+      off_t m = catalog[i].average[j].measureOffset;
+      off_t k;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+        MeasureTiny *measure = &catalog[i].measureT[m+k];
+        if (measure->dbFlags & ID_MEAS_OBJECT_HAS_2MASS) {
+          fprintf (stderr, "0x%08x 0x%08x : %12.8f %12.8f %5d\n", catalog[i].average[j].objID, catalog[i].average[j].catID, measure->R, measure->D, measure->photcode);
+        }
+      }
+    }
+  }
+  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 cat, off_t measOff) {
+  int k;
   /* 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;
+  FitAstromResult fitPos, fitPM, fitPar;
+  FitAstromResultInit (&fitPos);
+  FitAstromResultInit (&fitPM);
+  FitAstromResultInit (&fitPar);
   // if we fail to fit the astrometry for some reason, we need to set/reset these
 …
   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
 …
   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
+  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, FALSE);
+  // if there are no exposure detections, use the stack position
+  if (fitStats->Npoints < 1) {
     if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
       average[0].R  = average[0].Rstk;
 …
       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;
+  FitAstromPoints_Project (fitStats, &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) ***
 …
     if (Trange < PM_DT_MIN) {
       mode = FIT_AVERAGE;
       goto skipPM;
+    }
     if (N <= PM_TOOFEW) {
+      goto justPosition;
+    }
+    if (fitStats->Npoints <= PM_TOOFEW) {
       mode = FIT_AVERAGE;
+      goto skipPM;
+    }
+    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);
+      goto justPosition;
+    }
+    if (fitStats->NfitAlloc == 1) {
+      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
+      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
+    } else {
+      fitStats->Nfit = 0;
+      for (k = 0; k < fitStats->NfitAlloc; k++) {
+        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
+        if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
+        fitStats->Nfit ++;
+      }
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
+    }
+    FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
     // 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);
+    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
     if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
-    // XXX : does this make sense at 0,360 boundary?
     fitPM.p  = fitPM.dp  = 0.0;
 …
+  }
-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);
     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);
+      goto justPosition;
+    }
+    if (fitStats->Npoints <= PAR_TOOFEW) {
+      mode = FIT_PM_ONLY;
+      goto justPosition;
+    }
+    if (fitStats->NfitAlloc == 1) {
+      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
+      FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
+    } else {
+      fitStats->Nfit = 0;
+      for (k = 0; k < fitStats->NfitAlloc; k++) {
+        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
+        FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
+        fitStats->Nfit ++;
+      }
+      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
+      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
+      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
+    }
+    FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
+    // project Ro, Do back to RA,DEC
+    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
+    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     average[0].flags |= ID_STAR_FIT_PAR;
     fitStats->Npar ++;
-    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     // XXX a hard-wired hack...
     if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
+    if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) {
       mode = FIT_PM_ONLY;
+    }
+  }
 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...
+    if (fitStats->NfitAlloc == 1) {
+      FitAstromResultSetPM (&fitPos, 1, average);
+      FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
+    } else {
+      fitStats->Nfit = 0;
+      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
+      for (k = 0; k < fitStats->NfitAlloc; k++) {
+        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
+        FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
+        fitStats->Nfit ++;
+      }
+      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+    }
+    FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
     // 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 (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
     average[0].flags |= ID_STAR_FIT_AVE;
     fitStats->Nave ++;
+  }
+  // update the bit flags of which points were used
+  for (k = 0; k < fitStats->Npoints; k++) {
+    int Nm = fitStats->points[k].measure;
+    myAssert (Nm >= 0, "oops");
+    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
+    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
+  }
+  // we can set the star reference-image color only if we have loaded the image data
+  int setRefColor = areImagesMatched();
   if (setRefColor) {
+    float *C_blue = NULL;
+    float *C_red = NULL;
+    ALLOCATE (C_blue, float, fitStats->Npoints);
+    ALLOCATE (C_red, float, fitStats->Npoints);
+    int NcBlue = 0;
+    int NcRed = 0;
+    for (k = 0; k < fitStats->Npoints; k++) {
+      int Nm = fitStats->points[k].measure;
+      float colorBlue = getColorBlue (measOff + Nm, cat);
+      if (!isnan(colorBlue)) {
+        C_blue[NcBlue] = colorBlue;
+        NcBlue++;
+      }
+      float colorRed = getColorRed (measOff + Nm, cat);
+      if (!isnan(colorRed)) {
+        C_red[NcRed] = colorRed;
+        NcRed++;
+      }
+    }
+    // need to reassign here if isfinite()
     float colorMedian;
     dsort (C_blue, NcBlue);
+    fsort (C_blue, NcBlue);
     colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
     average[0].refColorBlue = colorMedian;
     dsort (C_red, NcRed);
+    fsort (C_red, NcRed);
     colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
     average[0].refColorRed = colorMedian;
+    free (C_blue);
+    free (C_red);
+  }
 …
   // XXXX for now, just use the mode as the result:
   int result = mode;
+  FitAstromResult fit;
+  FitAstromResultInit (&fit);
   switch (result) {
     case FIT_AVERAGE:
       average[0].flags |= ID_STAR_USE_AVE;
       fit = fitAve;
+      fit = fitPos;
       break;
     case FIT_PM_ONLY:
 …
     case FIT_PM_AND_PAR:
       average[0].flags |= ID_STAR_USE_PAR;
       fit = fitPAR;
+      fit = fitPar;
       break;
+  }
 …
   // what is the offset relative to the mean fit position?
   coords.crval1 = average[0].R;
   coords.crval2 = average[0].D;
   if (isnan(coords.crval1)) {
+  fitStats->coords.crval1 = average[0].R;
+  fitStats->coords.crval2 = average[0].D;
+  if (isnan(fitStats->coords.crval1)) {
     return (FALSE);
+  }
   if (isnan(coords.crval2)) {
+  if (isnan(fitStats->coords.crval2)) {
     return (FALSE);
+  }
   double dXoff, dYoff;
   RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
+  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &fitStats->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);
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
+    }
     fitStats->Noffset ++;
 …
                       average[0].uR,
                       average[0].uD,
                       fitAve.chisq, fitPM.chisq, fitPAR.chisq);
+                      fitPos.chisq, fitPM.chisq, fitPar.chisq);
   average[0].R          = fit.Ro; // RA in degrees
 …
   average[0].dP         = fit.dp; // parallax error in arcsec
   average[0].ChiSqAve   = fitAve.chisq;
+  average[0].ChiSqAve   = fitPos.chisq;
   average[0].ChiSqPM    = fitPM.chisq;
   average[0].ChiSqPar   = fitPAR.chisq;
   average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
+  average[0].ChiSqPar   = fitPar.chisq;
+  average[0].Tmean      = (Tmean * 86400 * 365.25) + fitStats->T2000;
   average[0].Trange     = (Trange * 86400 * 365.25);
   average[0].Npos       = fit.Nfit;
 …
 // 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
 …
   /* calculate the average value of R,D for a single star */
+  PMFit fitAve;
+  memset (&fitAve, 0, sizeof(fitAve));
+  fitAve.chisq = NAN;
+  FitAstromResult fitPos;
+  FitAstromResultInit (&fitPos);
   if (average[0].Nmeasure == 0) return TRUE;
-  int N = 0;
   int XVERB = FALSE;
 …
   XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
+  // select the measurements to be used in this analysis
+  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
+  // too few measurements for average position (require 2 values)
+  if (fitStats->Npoints < 1) return FALSE; // XXX ??
+  double Tmean, Trange, parRange;
+  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
+  FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
+  FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
+  // project Ro, Do back to RA,DEC
+  XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
+  // 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,
+                      fitPos.Ro, fitPos.Do,
+*(average[0].R - fitPos.Ro),
+*(average[0].D - fitPos.Do));
+  // make sure that the fit succeeded
+  int status = TRUE;
+  status &= finite(fitPos.Ro);
+  status &= finite(fitPos.Do);
+  status &= finite(fitPos.dRo);
+  status &= finite(fitPos.dDo);
+  if (!status) {
+    fitStats->Nskip ++;
+    return FALSE;
+  }
+  // what is the offset relative to the mean fit position?
+  fitStats->coords.crval1 = average[0].R;
+  fitStats->coords.crval2 = average[0].D;
+  double dXoff, dYoff;
+  RD_to_XY (&dXoff, &dYoff, fitPos.Ro, fitPos.Do, &fitStats->coords);
+  float dPos = hypot (dXoff, dYoff);
+  if (dPos > MaxMeanOffset) {
+    if (fitStats->Noffset < 100) {
+      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", fitStats->coords.crval1, fitStats->coords.crval2, fitPos.Ro, fitPos.Do, dXoff, dYoff);
+    }
+    fitStats->Noffset ++;
+    return FALSE;
+  }
+  // set the stack position values
+  average[0].Rstk  = fitPos.Ro; // RA in degrees
+  average[0].Dstk  = fitPos.Do; // DEC in degrees
+  average[0].dRstk = fitPos.dRo; // RA scatter in arcsec
+  average[0].dDstk = fitPos.dDo; // DEC scatter in arcsec
+  return (TRUE);
+}
+int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack) {
+  // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries
+  int has2MASS = FALSE;
+  int Npoints = fit->Npoints = 0;
+  FitAstromPoint *points = fit->points;
+  int TESTPT2 = FALSE;
+  TESTPT2 |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC);
+  TESTPT2 |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST);
+  if (TESTPT2) {
+    fprintf (stderr, "got test det\n");
+  }
   // 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) {
+    if (0) {
       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);
+      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 everything except gpc1 stack data
+    if (!isGPC1stack(measure[k].photcode)) continue;
+    // SKIP gpc1 forced-warp data
+    if (isGPC1warp(measure[k].photcode)) continue;
+    // SKIP gpc1 stack data
+    if (isStack) {
+      if (!isGPC1stack(measure[k].photcode)) continue;
+    } else {
+      if ( isGPC1stack(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);
 …
+    }
+    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);
+    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) {
+      fit->coords.crval1 = average[0].R;
+      fit->coords.crval2 = average[0].D;
+      double Xi, Yi;
+      RD_to_XY (&Xi, &Yi, Ri, Di, &fit->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;
+    }
+    FitAstromPointInit (&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 - fit->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(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++;
+    if (isnan(points[Npoints].dX)) continue;
+    if (isnan(points[Npoints].dY)) continue;
+    points[Npoints].dT = measure[k].dt;
+    points[Npoints].measure = k;
+    Npoints++;
+    if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
+      has2MASS = TRUE;
+    }
+    myAssert (Npoints <= fit->NpointsAlloc, "oops");
   } // 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,
+*(average[0].R - fitAve.Ro),
+*(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
+  int TESTPT = FALSE;
+  TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (average[0].catID == CAT_ID_SRC) && (average[0].objID == OBJ_ID_SRC);
+  TESTPT |= CAT_ID_DST && OBJ_ID_DST && (average[0].catID == CAT_ID_DST) && (average[0].objID == OBJ_ID_DST);
+  if (TESTPT) {
+    fprintf (stderr, "got test det\n");
+  }
+  // XXX flag measurements from stars with 2MASS
+  for (k = 0; k < average[0].Nmeasure; k++) {
+    // reset the bit to note that a detection was used (or not)
+    if (has2MASS) {
+      measure[k].dbFlags |=  ID_MEAS_OBJECT_HAS_2MASS;
+    } else {
+      measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_2MASS;
+    }
+  }
+  fit->Npoints = Npoints;
+  return TRUE;
+}
+int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
+  int k;
+  int Npoints = fitStats->Npoints;
+  FitAstromPoint *points = fitStats->points;
+  // find Tmin & Tmax from the list of accepted measurements
+  double Tmin  = points[0].T;
+  double Tmax  = points[0].T;
+  double pRmin = +2.0;
+  double pRmax = -2.0;
+  double pDmin = +2.0;
+  double pDmax = -2.0;
+  *Tmean = 0.0;
+  double Tsum = 0.0;
+  double Wsum = 0.0;
+  for (k = 0; k < Npoints; k++) {
+    Tmin = MIN(Tmin, points[k].T);
+    Tmax = MAX(Tmax, points[k].T);
+    float wx = 1.0 / SQ(points[k].dX);
+    Tsum += points[k].T * wx;
+    Wsum += wx;
+    // at this point, T is in years since J2000
+    ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
+    pRmin = MIN (pRmin, points[k].pR);
+    pRmax = MAX (pRmax, points[k].pR);
+    pDmin = MIN (pDmin, points[k].pD);
+    pDmax = MAX (pDmax, points[k].pD);
+  }
+  *Trange = Tmax - Tmin;
+  // mean epoch
+  *Tmean = Tsum / Wsum;
+  // for HIGH_SPEED, just use the center of the range
+  if (RELASTRO_OP == OP_HIGH_SPEED) {
+    *Tmean = 0.5*(Tmax - Tmin);
+  }
+  *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
+  /* we need to do the fit in a locally linear space; choose a ref coordinate */
+  fitStats->coords.crval1 = points[0].R;
+  fitStats->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, &fitStats->coords);
+    points[k].T -= *Tmean;
+  }
+  return TRUE;
+}
+int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
+  int i, j;
+  int 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;
+}
+int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *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;
+}
+// calculate mean and sigma points for the 5 fit parameter
+int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
+  // generate a histogram for the selected element
+  double *values = NULL;
+  ALLOCATE (values, double, Nfit);
+  int i;
+  for (i = 0; i < Nfit; i++) {
+    switch (mode) {
+      case FIT_RESULT_RA:
+        values[i] = fit[i].Ro;
+        break;
+      case FIT_RESULT_DEC:
+        values[i] = fit[i].Do;
+        break;
+      case FIT_RESULT_uR:
+        values[i] = fit[i].uR;
+        break;
+      case FIT_RESULT_uD:
+        values[i] = fit[i].uD;
+        break;
+      case FIT_RESULT_PLX:
+        values[i] = fit[i].p;
+        break;
+      default:
+        myAbort ("invalid option");
+    }
+  }
+  dsort (values, Nfit);
+  double median;
+  if (Nfit % 2) {
+    int Ncenter = Nfit / 2;
+    median = values[Ncenter];
+  } else {
+    int Ncenter = Nfit / 2 - 1;
+    median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
+  }
+  double Slo = VectorFractionInterpolate (values, 0.158655, Nfit);
+  double Shi = VectorFractionInterpolate (values, 0.841345, Nfit);
+  double sigma = (Shi - Slo) / 2.0;
+  switch (mode) {
+    case FIT_RESULT_RA:
+      result->Ro = median;
+      result->dRo = sigma;
+      break;
+    case FIT_RESULT_DEC:
+      result->Do = median;
+      result->dDo = sigma;
+      break;
+    case FIT_RESULT_uR:
+      result->uR = median;
+      result->duR = sigma;
+      break;
+    case FIT_RESULT_uD:
+      result->uD = median;
+      result->duD = sigma;
+      break;
+    case FIT_RESULT_PLX:
+      result->p = median;
+      result->dp = sigma;
+      break;
+    default:
+      myAbort ("invalid option");
+  }
+  return TRUE;
+}
+double VectorFractionInterpolate (double *values, float fraction, int Npts) {
+  float F = fraction * Npts;
+  int   N = fraction * Npts;
+  if (N < 0        ) return NAN;
+  if (N >= Npts - 2) return NAN;
+  // interpolate between N,N+1
+  double S = (F - N) * (values[N+1] - values[N]) + values[N];
+  return S;
+}
+int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) {
+  int i;
+  // add up the chi square for the fit
+  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);
+  }
+  switch (mode) {
+    case FIT_AVERAGE:
+      fit->chisq = chisq / (2.0*Npoints - 2.0);
+      break;
+    case FIT_PM_ONLY:
+      fit->chisq = chisq / (2.0*Npoints - 4.0);
+      break;
+    case FIT_PM_AND_PAR:
+      fit->chisq = chisq / (2.0*Npoints - 5.0);
+      break;
+    default:
+      myAbort ("invalid mode");
+  }
+  fit->Nfit = Npoints;
   return (TRUE);
+}
+/* fitting proper-motion and parallax:
+   given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
+   let x,y be the coordinate in the linearized frame with y parallel to DEC lines
+   L,B are the ecliptic longitude and latitude of the object,
+   dL and dB are the offsets in the L and B directions
+   dL = sin(t - topp)
+   dB = cos(t - topp)*sin(B)
+   these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
+   for the source in the x,y frame is:
+   x = Ro + uR * (t - to) + p * pR
+   y = Do + uD * (t - to) + p * pD
+   the unknowns in these equations are Ro, uR, Do, uD, and p
+   XXX think through the concepts for the pole a bit better.  all objects near the pole
+   move the same way with the same phase.  choose a projection center and define dL,dB relative
+   to that center point coordinate system?
+*/
+int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average) {
+  int i;
+  if (USE_GALAXY_MODEL) {
+    for (i = 0; i < Nfit; i++) {
+      fit->uR = average->uRgal;
+      fit->uD = average->uDgal;
+    }
+  } else {
+    for (i = 0; i < Nfit; i++) {
+      fit->uR = 0.0;
+      fit->uD = 0.0;
+    }
+  }
+  return TRUE;
+}

trunk/Ohana/src/relastro/src/args.c

-              r38441
+              r38986
+  }
+  // If we are looking at the whole sky (or whole relevant sky), use the full image table
+  // -- this save substantial memory.  this could be automatic if the skyregion covers
+  // more than 2pi.
+  USE_ALL_IMAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-use-all-images"))) {
+    remove_argument (N, &argc, argv);
+    USE_ALL_IMAGES = TRUE;
+  }
   USE_BASIC_CHECK = FALSE;
   if ((N = get_argument (argc, argv, "-basic-image-search"))) {
 …
     remove_argument (N, &argc, argv);
     APPLY_OFFSETS = TRUE;
+  }
+  N_BOOTSTRAP_SAMPLES = 1;
+  if ((N = get_argument (argc, argv, "-bootstrap-samples"))) {
+    remove_argument (N, &argc, argv);
+    N_BOOTSTRAP_SAMPLES = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+    if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) {
+      fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n");
+      exit (2);
+    }
+  }
 …
+  }
+  // If we are looking at the whole sky (or whole relevant sky), use the full image table
+  // -- this save substantial memory.  this could be automatic if the skyregion covers
+  // more than 2pi.
+  USE_ALL_IMAGES = FALSE;
+  if ((N = get_argument (argc, argv, "-use-all-images"))) {
+    remove_argument (N, &argc, argv);
+    USE_ALL_IMAGES = TRUE;
+  }
   USE_BASIC_CHECK = FALSE;
   if ((N = get_argument (argc, argv, "-basic-image-search"))) {
 …
     remove_argument (N, &argc, argv);
     MaxDensityUse = TRUE;
+  }
+  N_BOOTSTRAP_SAMPLES = 1;
+  if ((N = get_argument (argc, argv, "-bootstrap-samples"))) {
+    remove_argument (N, &argc, argv);
+    N_BOOTSTRAP_SAMPLES = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+    if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) {
+      fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n");
+      exit (2);
+    }
+  }

trunk/Ohana/src/relastro/src/assign_images.c

r38441	r38986
26	26	char mapfile[DVO_MAX_PATH];
27	27	snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
28		AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE);
	28	AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE);
29	29
30	30	// assign images.coords.offsetMap -> table->map[i]

trunk/Ohana/src/relastro/src/bcatalog.c

-              r38441
+              r38986
     // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly
+    if (USE_GALAXY_MODEL && !catalog[0].average[i].Nstarpar) {
+      // we may want to add some starpar quality filters...
+      continue;
+    if (USE_GALAXY_MODEL) {
+      // if we use the galaxy model, we need to require the model pm to exist
+      if (!catalog[0].average[i].Nstarpar) continue;
+      if (!isfinite(catalog[0].average[i].uRgal)) continue;
+      if (!isfinite(catalog[0].average[i].uDgal)) continue;
+    }
+    int TESTPT = FALSE;
+    TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (catalog[0].average[i].catID == CAT_ID_SRC) && (catalog[0].average[i].objID == OBJ_ID_SRC);
+    TESTPT |= CAT_ID_DST && OBJ_ID_DST && (catalog[0].average[i].catID == CAT_ID_DST) && (catalog[0].average[i].objID == OBJ_ID_DST);
+    if (TESTPT) {
+      fprintf (stderr, "got test det\n");
+    }

trunk/Ohana/src/relastro/src/high_speed_catalogs.c

-              r38471
+              r38986
               PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    strextend (&command, "relastro_client -high-speed %s %s %f %s", PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, outputDir);
+    strextend (&command, "-hostID %d", table->hosts[i].hostID);
+    strextend (&command, "-D CATDIR %s", CATDIR);
+    strextend (&command, "-hostdir %s", table->hosts[i].pathname);
+    strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
     free (outputDir);
 …
     if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); }
+    if (USE_ALL_IMAGES)      { strextend (&command, "-use-all-images"); }
     if (USE_BASIC_CHECK)     { strextend (&command, "-basic-image-search"); }
     if (FlagOutlier)         { strextend (&command, "-clip %d", CLIP_THRESH); }
 …
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+    if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
+    if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
+    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
+    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
+    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     if (WHERE_A[0]) { strextend (&command, "-D WHERE_A \"%s\"", WHERE_A); }

trunk/Ohana/src/relastro/src/high_speed_objects.c

-              r38441
+              r38986
   fprintf (stderr, "found %d matches\n", Nmatch);
+  dvo_catalog_save (&catalogOut, VERBOSE2);
+  dvo_catalog_unlock (&catalogOut);
+  SetProtect (TRUE);
+  if (!dvo_catalog_save (&catalogOut, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalogOut.filename); exit (1); }
+  if (!dvo_catalog_unlock (&catalogOut)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalogOut.filename); exit (1); }
+  SetProtect (FALSE);
   dvo_catalog_free (&catalogOut);
   free (slowMoving);

trunk/Ohana/src/relastro/src/hpm_catalogs.c

-              r38471
+              r38986
     char *command = NULL;
+    strextend (&command, "relastro_client -hpm %f %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f",
+              RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    strextend (&command, "relastro_client -hpm %f %s", RADIUS, outputDir);
+    strextend (&command, " -hostID %d", table->hosts[i].hostID);
+    strextend (&command, " -D CATDIR %s", CATDIR);
+    strextend (&command, " -hostdir %s", table->hosts[i].pathname);
+    strextend (&command, " -region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
     free (outputDir);
 …
     if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); }
+    if (USE_ALL_IMAGES)      { strextend (&command, "-use-all-images"); }
     if (USE_BASIC_CHECK)     { strextend (&command, "-basic-image-search"); }
     if (FlagOutlier)         { strextend (&command, "-clip %d", CLIP_THRESH); }
 …
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+    if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
+    if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
+    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
+    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
+    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     if (TimeSelect) {

trunk/Ohana/src/relastro/src/hpm_objects.c

-              r38441
+              r38986
   if (catalog[0].Naverage == Nslow) {
     fprintf (stderr, "no possible fast objects, skipping this catalog\n");
+    dvo_catalog_save (&catalogOut, VERBOSE2);
+    dvo_catalog_unlock (&catalogOut);
+    SetProtect (TRUE);
+    if (!dvo_catalog_save (&catalogOut, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalogOut.filename); exit (1); }
+    if (!dvo_catalog_unlock (&catalogOut)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalogOut.filename); exit (1); }
+    SetProtect (FALSE);
     dvo_catalog_free (&catalogOut);
     free (slowMoving);

trunk/Ohana/src/relastro/src/launch_region_hosts.c

r38441	r38986
6	6	int i;
7	7
8		// do not remote the sync and fits files if we do a manual run -- user must clear if needed
	8	// do not remove the sync and fits files if we do a manual run -- user must clear if needed
9	9	if (!PARALLEL_REGIONS_MANUAL) {
10	10	// clear the I/O files

trunk/Ohana/src/relastro/src/load_catalogs.c

-              r38471
+              r38986
     char *command = NULL;
+    strextend (&command, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f -D RELASTRO_SIGMA_LIM %f",
+              table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM);
+    strextend (&command, "relastro_client -load-objects %s", table->hosts[i].results);
+    strextend (&command, " -hostID %d", table->hosts[i].hostID);
+    strextend (&command, " -D CATDIR %s", CATDIR);
+    strextend (&command, " -hostdir %s", table->hosts[i].pathname);
+    strextend (&command, " -region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    strextend (&command, " -statmode %s", STATMODE);
+    strextend (&command, " -minerror %f", MIN_ERROR);
+    strextend (&command, " -D RELASTRO_SIGMA_LIM %f", SIGMA_LIM);
     if (FIT_MODE == FIT_PM_ONLY)         strextend (&command, "-pm");
 …
+    }
+    if (USE_ALL_IMAGES)      strextend (&command, "-use-all-images");
     if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1");
     if (PHOTCODE_KEEP_LIST)  strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST);

trunk/Ohana/src/relastro/src/load_images.c

-              r38441
+              r38986
  */
 int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) {
+int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap, int UseAllImages) {
   Image     *image, *subset;
   off_t      Nimage, Nsubset;
   off_t     *LineNumber;
+  off_t     *LineNumber, i;
   INITTIME;
 …
   char mapfile[DVO_MAX_PATH];
   snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR);
   table = AstromOffsetMapLoad (mapfile, VERBOSE);
+  table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE);
   // assign images.coords.offsetMap -> table->map[i]
 …
   // select the images which overlap the selected sky regions
+  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap);
+  MARKTIME("  select images: %f sec\n", dtime);
+  if (UseAllImages) {
+    ALLOCATE (LineNumber, off_t, Nimage);
+    for (i = 0; i < Nimage; i++) LineNumber[i] = i;
+    subset  = image;
+    Nsubset = Nimage;
+  } else {
+    subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap);
+    MARKTIME("  select images: %f sec\n", dtime);
+  if (Nsubset == Nimage) {
+    free (subset);
+    subset = image;
+    if (Nsubset == Nimage) {
+      free (subset);
+      subset = image;
+    }
+  }
   initImages (subset, LineNumber, Nsubset);
   MARKTIME("  init images: %f sec\n", dtime);
   initMosaics (subset, Nsubset);
   MARKTIME("  init mosaics: %f sec\n", dtime);

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

r36630	r38986
4	4
5	5	/* get configuration info, args */
	6	SetSignals ();
6	7	initialize (argc, argv);
7	8

trunk/Ohana/src/relastro/src/relastro_client.c

-              r38441
+              r38986
   // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc)
+  SetSignals ();
   initialize_client (argc, argv);
 …
       /* load regions and images based on specified sky patch (default depth) */
       load_images (&db, skylist, FALSE);
+      load_images (&db, skylist, FALSE, USE_ALL_IMAGES);
       // I can free the database after I have loaded the images...

trunk/Ohana/src/relastro/src/relastro_images.c

-              r38441
+              r38986
   /* load regions and images based on specified sky patch (default depth) (require full overlap) */
   load_images (&db, skylist, TRUE);
+  load_images (&db, skylist, TRUE, USE_ALL_IMAGES);
   MARKTIME("load images: %f sec\n", dtime);
 …
       if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog);
       for (i = 0; i < NLOOP; i++) {
+        UpdateObjects (catalog, Ncatalog, (i > 0)); // calculate <R>,<D>; if (i > 0), apply Galaxy Motion Model (if desired)
+        if ((i > 1) || !USE_GALAXY_MODEL) {
+          // if GALAXY_MODEL is selected, we want to delay the frame correction until we have
+          // applied the galaxy model a couple of times.
+          FrameCorrectionSerial (catalog, Ncatalog);
+        }
+        UpdateChips (catalog, Ncatalog);   // measure.X,Y -> R,D, fit image.coords
+        UpdateObjects (catalog, Ncatalog, i); // calculate <R>,<D>; if (i > 0), apply Galaxy Motion Model (if desired)
+        UpdateChips (catalog, Ncatalog, i);   // measure.X,Y -> R,D, fit image.coords
         MARKTIME("update chips: %f sec\n", dtime);
+      }
 …
   freeImageBins (1);
+  // XXX why do I do this exactly?
+  reload_images (&db);
+  // If we did NOT use all images, then we applied the measured corrections to a subset of
+  // images.  we now need to set the images associated with db to have those values so
+  // they will be written out by the dvo_image_update() operations.  If we USE_ALL_IMAGES,
+  // then we do not need to do this, and we should call dvo_image_save (not dvo_image_update)
+  if (!USE_ALL_IMAGES) {
+    reload_images (&db);
+  }
   if (PARALLEL) {
 …
     // need to also save the image map table...
     save_astrom_table ();
+    dvo_image_update (&db, VERBOSE);
+    if (USE_ALL_IMAGES) {
+      dvo_image_save (&db, VERBOSE);
+    } else {
+      dvo_image_update (&db, VERBOSE);
+    }
     dvo_image_unlock (&db);
+  }
+  // NOTE: if we have parallel partitions, then we need to save the Images.dat and
+  // AstroMap.fits tables BEFORE the remote relastro_clients are launched (they load
+  // Images.dat and AstroMap.fits)
+  // if we do NOT have parallel partitions, we must NOT write them out yet: the act of
+  // writing out the files byte-swaps the data and makes the values invalide for the
+  // following calls to UpdateObjectOffsets (which attempt to apply the image values from
+  // the structure in memory)
   // iterate over catalogs to make detection coordinates consistant
+  UpdateObjectOffsets (skylist, 0, NULL);
+  if (APPLY_OFFSETS) {
+    UpdateObjectOffsets (skylist, 0, NULL);
+  }
   if (!PARALLEL) {
     // save the updated image parameters
     save_astrom_table ();
+    dvo_image_update (&db, VERBOSE);
+    if (USE_ALL_IMAGES) {
+      dvo_image_save (&db, VERBOSE);
+    } else {
+      dvo_image_update (&db, VERBOSE);
+    }
     dvo_image_unlock (&db);
+  }

trunk/Ohana/src/relastro/src/relastro_merge_source.c

-              r38441
+              r38986
     resort_catalog (&catalog_src);
+    dvo_catalog_save (&catalog_src, VERBOSE2);
+    dvo_catalog_unlock (&catalog_src);
+    SetProtect (TRUE);
+    if (!dvo_catalog_save (&catalog_src, VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog_src.filename); exit (1); }
+    if (!dvo_catalog_unlock (&catalog_src)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog_src.filename); exit (1); }
+    SetProtect (FALSE);
     dvo_catalog_free (&catalog_src);

trunk/Ohana/src/relastro/src/relastro_objects.c

-              r38471
+              r38986
     char *command = NULL;
+    strextend (&command, "relastro_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s",
+              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE);
+    strextend (&command, "relastro_client -update-objects");
+    strextend (&command, "-hostID %d", table->hosts[i].hostID);
+    strextend (&command, "-D CATDIR %s", CATDIR);
+    strextend (&command, "-hostdir %s", table->hosts[i].pathname);
+    strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    strextend (&command, "-statmode %s", STATMODE);
     if (FIT_MODE == FIT_PM_ONLY)         { strextend (&command, "-pm"); }
 …
     if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); }
     if (FlagOutlier)   { strextend (&command, "-clip %d", CLIP_THRESH); }
+    if (USE_ALL_IMAGES)      { strextend (&command, "-use-all-images"); }
     if (USE_FIXED_PIXCOORDS) { strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); }
     if (PHOTCODE_KEEP_LIST)  { strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); }
 …
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+    if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
+    if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
+    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
+    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
+    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     if (TimeSelect) {

trunk/Ohana/src/relastro/src/relastro_parallel_images.c

-              r38062
+              r38986
   char mapfile[DVO_MAX_PATH];
   snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.%d.fits", CATDIR, REGION_HOST_ID);
   AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE);
+  AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, 100000, VERBOSE);
   // assign images.coords.offsetMap -> table->map[i]
 …
+      }
       for (i = 0; i < NLOOP; i++) {
         UpdateObjects (catalog, Ncatalog, (i > 0));
+        UpdateObjects (catalog, Ncatalog, i);
         LOGRTIME("UpdateObjects loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime);
+        if ((i > 1) || !USE_GALAXY_MODEL) {
+          // if GALAXY_MODEL is selected, we want to delay the frame correction until we have
+          // applied the galaxy model a couple of times.
+          FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i);
+          LOGRTIME("FrameCorrection loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime);
+        }
         share_mean_pos (catalog, Ncatalog, regionHosts, i);
         LOGRTIME("share_mean_pos loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime);
 …
         LOGRTIME("slurp_mean_pos loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime);
         UpdateChips (catalog, Ncatalog);
+        UpdateChips (catalog, Ncatalog, i);
         LOGRTIME("UpdateChips loop %d on %s, host %d: %f sec\n", i, myHostName, REGION_HOST_ID, dtime);

trunk/Ohana/src/relastro/src/save_catalogs.c

-              r33652
+              r38986
     if (VERBOSE2) fprintf (stderr, "saving catalog %s\n", catalog[i].filename);
+    dvo_catalog_save (&catalog[i], VERBOSE2);
+    dvo_catalog_unlock (&catalog[i]);
+    SetProtect (TRUE);
+    if (!dvo_catalog_save (&catalog[i], VERBOSE2)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog[i].filename); exit (1); }
+    if (!dvo_catalog_unlock (&catalog[i])) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog[i].filename); exit (1); }
+    SetProtect (FALSE);
     dvo_catalog_free (&catalog[i]);
+  }

trunk/Ohana/src/relastro/src/share_images_pos.c

r38062	r38986
115	115	snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, regionHosts->hosts[i].hostID);
116	116
117		AstromOffsetTable *table = AstromOffsetMapLoad (mapname, VERBOSE);
	117	AstromOffsetTable *table = AstromOffsetMapLoad (mapname, 100000, VERBOSE);
118	118	LOGRTIME("image_maps_load host %d loop %d on %s, host %d: %f sec\n", i, nloop, myHostName, REGION_HOST_ID, dtime);
119	119

Context Navigation

Legend:

trunk/Ohana

trunk/Ohana/src/relastro/Makefile

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

trunk/Ohana/src/relastro/src/FitChip.c

trunk/Ohana/src/relastro/src/FitPM.c

trunk/Ohana/src/relastro/src/FitPMandPar.c

trunk/Ohana/src/relastro/src/FrameCorrection.c

trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c

trunk/Ohana/src/relastro/src/GetAstromError.c

trunk/Ohana/src/relastro/src/ImageOps.c

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

trunk/Ohana/src/relastro/src/Shutdown.c

trunk/Ohana/src/relastro/src/StarMaps.c

trunk/Ohana/src/relastro/src/UpdateChips.c

trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c

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

trunk/Ohana/src/relastro/src/args.c

trunk/Ohana/src/relastro/src/assign_images.c

trunk/Ohana/src/relastro/src/bcatalog.c

trunk/Ohana/src/relastro/src/high_speed_catalogs.c

trunk/Ohana/src/relastro/src/high_speed_objects.c

trunk/Ohana/src/relastro/src/hpm_catalogs.c

trunk/Ohana/src/relastro/src/hpm_objects.c

trunk/Ohana/src/relastro/src/launch_region_hosts.c

trunk/Ohana/src/relastro/src/load_catalogs.c

trunk/Ohana/src/relastro/src/load_images.c

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

trunk/Ohana/src/relastro/src/relastro_client.c

trunk/Ohana/src/relastro/src/relastro_images.c

trunk/Ohana/src/relastro/src/relastro_merge_source.c

trunk/Ohana/src/relastro/src/relastro_objects.c

trunk/Ohana/src/relastro/src/relastro_parallel_images.c

trunk/Ohana/src/relastro/src/save_catalogs.c

trunk/Ohana/src/relastro/src/share_images_pos.c

Download in other formats: