Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 16140)
+++ trunk/Ohana/src/relastro/Makefile	(revision 16810)
@@ -37,4 +37,6 @@
 $(SRC)/UpdateObjects.$(ARCH).o   \
 $(SRC)/UpdateSimple.$(ARCH).o    \
+$(SRC)/UpdateMeasures.$(ARCH).o  \
+$(SRC)/GetAstromError.$(ARCH).o  \
 $(SRC)/args.$(ARCH).o		 \
 $(SRC)/bcatalog.$(ARCH).o	 \
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 16140)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 16810)
@@ -10,4 +10,6 @@
 } CoordMode;
 
+typedef enum {ERROR_MODE_RA, ERROR_MODE_DEC, ERROR_MODE_POS} ErrorMode;
+
 typedef struct {
   double R, D;  /* Sky Coords    - degrees */
@@ -15,5 +17,7 @@
   double L, M;  /* Focal Plane   - pixels  */
   double X, Y;  /* Chip Coords   - pixels  */
-  double Mag, dMag;
+  double Mag;
+  double dMag;
+  double dPos;
   int mask;
 } StarData;
@@ -263,9 +267,10 @@
 void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
 void fit_eval (CoordFit *fit);
+void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1);
 double **poly2d_dx (double **poly, int Nx, int Ny);
 double **poly2d_dy (double **poly, int Nx, int Ny);
 double **poly2d_copy (double **poly, int Nx, int Ny);
 double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
-void fit_apply_coords (CoordFit *fit, Coords *coords);
+CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords);
 int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
 CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
@@ -278,4 +283,7 @@
 int UpdateChips (Catalog *catalog, int Ncatalog);
 int UpdateMosaic (Catalog *catalog, int Ncatalog);
+int UpdateMeasures (Catalog *catalog, int Ncatalog);
+void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
+
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
 int ParFactor (double *pR, double *pD, double R, double D, time_t T);
Index: trunk/Ohana/src/relastro/src/FitChip.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitChip.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/FitChip.c	(revision 16810)
@@ -1,24 +1,67 @@
 # include "relastro.h"
+# define SCATTER_MAX_ERROR 0.05
 
 void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
 
-  int i;
+  int i, Nscatter, Niter;
   CoordFit *fit;
+  double dL, dM, dR, dRmax, *values;
 
-  fit = fit_init (coords[0].Npolyterms);
-  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, 1.0);
-  }
-  fit_eval (fit);
-  fit_apply_coords (fit, coords);
-  fit_free (fit);
+  ALLOCATE (values, double, Nmatch);
 
-  // apply new coords to raw (X,Y -> L,M)
-  for (i = 0; i < Nmatch; i++) {
-    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
+  for (Niter = 0; Niter < 3; Niter ++) {
+    // measure the scatter distribution (use only the bright end detections)
+    for (i = Nscatter = 0; i < Nmatch; i++) {
+      if (raw[i].mask) continue;
+      if (raw[i].dMag > SCATTER_MAX_ERROR) continue;
+
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+    
+      values[Nscatter] = dR;
+      Nscatter++;
+    }
+
+    // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
+    dsort (values, Nscatter);
+    dRmax = 3.0*values[(int)(0.40*Nscatter)];
+
+    fit = fit_init (coords[0].Npolyterms);
+
+    for (i = 0; i < Nmatch; i++) {
+      if (raw[i].mask) continue;
+
+      // require radius of XXX arcsec
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+
+      if (dR > dRmax) {
+	continue;
+	raw[i].mask = TRUE;
+      }
+    
+      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
+    }
+    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n",
+	     dRmax, Nscatter, fit[0].Npts, Nmatch);
+
+    if (fit[0].Npts < 25) {
+      fit_free (fit);
+      free (values);
+      return;
+    }
+    fit_eval (fit);
+    fit_apply_coords (fit, coords);
+    fit_free (fit);
+
+    for (i = 0; i < Nmatch; i++) {
+      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
+    }
   }
 
-  fit_free (fit);
+  free (values);
+  return;
 }
 
@@ -34,3 +77,28 @@
 */
 
-/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+/* example using fit_apply() :
+
+  f = fopen ("test3.dat", "w");
+
+  // apply new coords to raw (X,Y -> L,M)
+  for (i = 0; i < Nmatch; i++) {
+    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
+    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
+    fprintf (f, "%f %f\n", L1, M1);
+  }
+  fclose (f);
+*/
+
+/*** XXX this function can be improved in 4 important ways:
+
+1) the initial value of the clipping radius could be set based on the distribution of the
+   value (eg, inner 50%) (better to use dL, dM rather than dP, dQ)
+
+2) multiple clipping passes could be performed
+
+3) the per-star astrometric errors could be included in the fit (non-systematic terms at
+   least, or scaled values of the mag error, or simply the mag error itself)
+
+4) we could ignore input detections on the basis of the photFlags
+
+*/
Index: trunk/Ohana/src/relastro/src/FitMosaic.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitMosaic.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/FitMosaic.c	(revision 16810)
@@ -5,9 +5,23 @@
   int i;
   CoordFit *fit;
+  double dP, dQ, dR;
 
   fit = fit_init (coords[0].Npolyterms);
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
+
+    // require radius of XXX arcsec
+    dP = raw[i].P - ref[i].P;
+    dQ = raw[i].Q - ref[i].Q;
+    dR = 3600.0 * hypot (dP, dQ);
+
+    // XXX the value needs to be set in a more intelligent way
+    if (dR > 0.15) continue;
+    
     fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0);
+  }
+  if (fit[0].Npts == 0) {
+    fit_free (fit);
+    return;
   }
   fit_eval (fit);
Index: trunk/Ohana/src/relastro/src/FitSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitSimple.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/FitSimple.c	(revision 16810)
@@ -5,9 +5,23 @@
   int i;
   CoordFit *fit;
+  double dP, dQ, dR;
 
   fit = fit_init (coords[0].Npolyterms);
   for (i = 0; i < Nmatch; i++) {
     if (raw[i].mask) continue;
+
+    // require radius of XXX arcsec
+    dP = raw[i].P - ref[i].P;
+    dQ = raw[i].Q - ref[i].Q;
+    dR = 3600.0 * hypot (dP, dQ);
+
+    // XXX the value needs to be set in a more intelligent way
+    if (dR > 0.15) continue;
+    
     fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);
+  }
+  if (fit[0].Npts == 0) {
+    fit_free (fit);
+    return;
   }
   fit_eval (fit);
@@ -21,6 +35,4 @@
     raw[i].Q = raw[i].M;
   }
-
-  fit_free (fit);
 }
 
@@ -37,2 +49,4 @@
 
 /* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+/* XXX See notes in FitChips.c */
Index: trunk/Ohana/src/relastro/src/GetAstromError.c
===================================================================
--- trunk/Ohana/src/relastro/src/GetAstromError.c	(revision 16810)
+++ trunk/Ohana/src/relastro/src/GetAstromError.c	(revision 16810)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+float GetAstromError (Measure *measure, int mode) {
+
+  PhotCode *code;
+  float dPobs, dPsys, dPtotal, dM, AS, MS;
+
+  switch (mode) {
+    case ERROR_MODE_RA:
+      dPobs = measure[0].dXccd;  // need to redefine this as RAerr
+      break;
+    case ERROR_MODE_DEC:
+      dPobs = measure[0].dYccd;  // need to redefine this as RAerr
+      break;
+    case ERROR_MODE_POS:
+      dPobs = hypot (measure[0].dXccd, measure[0].dYccd);  // need to redefine this as RAerr
+      break;
+    default:
+      abort();
+  }
+
+  /* the astrometric errors are not being carried yet (but should be!) */
+  /* we use the photometric mag error as a weighting term */
+
+  code 	= GetPhotcodebyCode (measure[0].photcode);
+  AS   	= code[0].astromErrScale;
+  MS   	= code[0].astromErrMagScale;
+  dPsys = code[0].astromErrSys;
+  dM    = measure[0].dM;
+  
+  dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
+  dPtotal = MAX (dPtotal, MIN_ERROR);
+
+  return (dPtotal);
+}
+
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 16810)
@@ -33,6 +33,6 @@
 
   for (i = 0; i < Nimage; i++) {
-    start[i] = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
-    stop[i]  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
+    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
   }
 }
@@ -81,4 +81,5 @@
 
   int i, j;
+  char *name;
 
   for (i = 0; i < Ncatalog; i++) {
@@ -87,4 +88,10 @@
     }
   }
+
+  for (i = 0; VERBOSE && (i < Nimage); i++) {
+    name = GetPhotcodeNamebyCode (image[i].photcode);
+    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i], 
+	     ohana_sec_to_date(image[i].tzero), name);
+  } 
 }
 
@@ -182,7 +189,75 @@
 // return StarData values for detections in the specified image, converting coordinates from the
 // chip positions: X,Y -> L,M -> P,Q -> R,D
+void fixImageRaw (Catalog *catalog, int Ncatalog, int im) {
+
+  int i, m, c, n;
+  double X, Y, L, M, P, Q, R, D, dR, dD;
+  
+  Mosaic *mosaic;
+  Coords *moscoords, *imcoords;
+  
+  moscoords = NULL;
+  mosaic = getMosaicForImage (im);
+  if (mosaic != NULL) {
+      moscoords = &mosaic[0].coords;
+  }
+  imcoords = &image[im].coords;
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+
+    X = catalog[c].measure[m].Xccd;
+    Y = catalog[c].measure[m].Yccd;
+    n = catalog[c].measure[m].averef;
+
+    if (moscoords == NULL) {
+      // this is a Simple image (not a mosaic)
+      // note that for a Simple image, L,M = P,Q
+      XY_to_LM (&L, &M, X, Y, imcoords);
+      LM_to_RD (&R, &D, L, M, imcoords);
+    } else {
+      XY_to_LM (&L, &M, X, Y, imcoords);
+      XY_to_LM (&P, &Q, L, M, moscoords);
+      LM_to_RD (&R, &D, P, Q, moscoords);
+    }
+
+    // new dR, dD : test
+    dR = 3600.0*(catalog[c].average[n].R - R);
+    dD = 3600.0*(catalog[c].average[n].D - D);
+    
+    if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
+      // XXXXX running into this still for last megacam exposure: wrong mosaic?
+      // ???? inconsistently hitting this????
+      fprintf (stderr, "!");
+      abort ();
+    }
+    if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
+      fprintf (stderr, "*");
+      abort ();
+    }
+
+    catalog[c].measure[m].dR = dR;
+    catalog[c].measure[m].dD = dD;
+
+    if (catalog[c].measure[m].dR > +180.0*3600.0) {
+      // average on high end of boundary, move star up
+      R += 360.0;
+      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+    if (catalog[c].measure[m].dR < -180.0*3600.0) {
+      // average on low end of boundary, move star down
+      R -= 360.0;
+      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+  }  
+  return;
+}
+
+// return StarData values for detections in the specified image, converting coordinates from the
+// chip positions: X,Y -> L,M -> P,Q -> R,D
 StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
 
-  int i, m, c;
+  int i, m, c, n;
   
   Mosaic *mosaic;
@@ -213,6 +288,12 @@
     raw[i].Mag  = catalog[c].measure[m].M;
     raw[i].dMag = catalog[c].measure[m].dM;
+    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
+
+    n = catalog[c].measure[m].averef;
 
     raw[i].mask = FALSE;
+    if (catalog[c].average[n].Nmeasure < 2) {
+      raw[i].mask = TRUE;
+    }
 
     switch (mode) {
@@ -273,4 +354,5 @@
     ref[i].Mag  = catalog[c].measure[m].M;
     ref[i].dMag = catalog[c].measure[m].dM;
+    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
 
     ref[i].mask = FALSE;
@@ -286,5 +368,5 @@
       case MODE_MOSAIC:
       RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
-      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords);
+      LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
       LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
       break;
Index: trunk/Ohana/src/relastro/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 16810)
@@ -38,6 +38,6 @@
 
     /* set image time range */
-    start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
-    stop  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
+    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
 
     /* a new mosaic, define ranges */
@@ -77,6 +77,6 @@
 
     /* set image time range */
-    start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
-    stop  = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1);
+    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
 
     /* find existing mosaic with this time range */
Index: trunk/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 16810)
@@ -21,4 +21,5 @@
     ref = getImageRef (catalog, Ncatalog, i, &Nstars, MODE_MOSAIC);
 
+    fprintf (stderr, "image %d : Nstars: %d\n", i, Nstars);
     FitChip (raw, ref, Nstars, &image[i].coords);
 
Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateMeasures.c	(revision 16810)
+++ trunk/Ohana/src/relastro/src/UpdateMeasures.c	(revision 16810)
@@ -0,0 +1,21 @@
+# include "relastro.h"
+
+int UpdateMeasures (Catalog *catalog, int Ncatalog) {
+
+  int i, Nimage;
+  Image *image;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip DIS images (since they have no associated detections) */
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* convert measure coordinates to raw entries */
+    fixImageRaw (catalog, Ncatalog, i);
+  }
+
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 16810)
@@ -18,5 +18,5 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nm);
+      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
     }
   }
@@ -48,5 +48,4 @@
   double Tmin, Tmax;
   float errorScale;
-  PhotCode *code;
 
   initObjectData (catalog, Ncatalog);
@@ -81,10 +80,10 @@
 
       N = 0;
-      m = catalog[i].average[j].offset;
+      m = catalog[i].average[j].measureOffset;
 
       Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
       mode = FIT_MODE;
 
-      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
 	
@@ -96,11 +95,6 @@
 	Tmax = MAX(Tmax, T[N]);
 
-	/* the astrometric errors are not being carried yet (but should be!) */
-	/* we use the photometric mag error as a weighting term */
-
-	code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
-	errorScale = code[0].astromScale;
-	dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
-	dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
+	dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
+	dD[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
 	dT[N] = catalog[i].measure[m].dt;
 
@@ -193,6 +187,6 @@
 
       // the measure fields must be updated before the average fields
-      m = catalog[i].average[j].offset;
-      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+      m = catalog[i].average[j].measureOffset;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
 	setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 16810)
@@ -202,16 +202,27 @@
   fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
   fprintf (stderr, "  working options: \n");
-  fprintf (stderr, "  -time (start) (stop)\n");
-  fprintf (stderr, "  -v\n");
+  fprintf (stderr, "  -update-objects\n");
+  fprintf (stderr, "    -pm\n");
+  fprintf (stderr, "    -par\n");
+  fprintf (stderr, "    -pmpar\n");
+  fprintf (stderr, "  -update-simple\n");
+  fprintf (stderr, "  -update-chips\n");
+  fprintf (stderr, "  -update-mosaics\n");
+  fprintf (stderr, "  -time (start)(stop)\n");
+  fprintf (stderr, "  -photcode (code)[,code,code...]\n");
   fprintf (stderr, "  -plot\n");
   fprintf (stderr, "  -plotdelay (seconds)\n");
-  fprintf (stderr, "  -update\n");
+  fprintf (stderr, "  -statmode (mode)\n");
+  fprintf (stderr, "  -reset");
+  fprintf (stderr, "  -update : apply new fit to database\n");
   fprintf (stderr, "  -params\n");
-  fprintf (stderr, "  -reset\n");
+  fprintf (stderr, "  -plrange\n");
+  fprintf (stderr, "  -minerror\n");
+  fprintf (stderr, "  -area\n");
   fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
   fprintf (stderr, "  -instmag min max\n\n");
-  fprintf (stderr, "  planned options: \n");
-  fprintf (stderr, "  -photcode code[,code,code,..]\n");
+  fprintf (stderr, "  -v\n");
   fprintf (stderr, "  \n");
   exit (2);
 } 
+
Index: trunk/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relastro/src/bcatalog.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/bcatalog.c	(revision 16810)
@@ -22,9 +22,9 @@
   /* exclude stars not in range or with too few measurements */
   for (i = 0; i < catalog[0].Naverage; i++) {
-    if (catalog[0].average[i].Nm < 2) continue; 
+    if (catalog[0].average[i].Nmeasure < 2) continue; 
 
     /* start with all stars good */
     subcatalog[0].average[Naverage] = catalog[0].average[i];
-    subcatalog[0].average[Naverage].offset = Nmeasure;
+    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
     for (j = 0; j < PhotNsec; j++) {
       subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j];
@@ -38,7 +38,7 @@
 
     Nm = 0;
-    for (j = 0; j < catalog[0].average[i].Nm; j++) {
+    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
 
-      offset = catalog[0].average[i].offset + j;
+      offset = catalog[0].average[i].measureOffset + j;
 
       /* select measurements by photcode, if specified */
@@ -95,5 +95,5 @@
       continue; 
     }
-    subcatalog[0].average[Naverage].Nm = Nm;
+    subcatalog[0].average[Naverage].Nmeasure = Nm;
     Naverage ++;
     if (Naverage == NAVERAGE) {
Index: trunk/Ohana/src/relastro/src/fitpoly.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpoly.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/fitpoly.c	(revision 16810)
@@ -128,4 +128,14 @@
       matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
     }
+
+    // mask the terms not represented by the Coords terms
+    if (ix + iy > fit[0].Norder) {
+      for (j = 0; j < fit[0].Nelems; j++) {
+	matrix[i][j] = 0.0;
+      }
+      vector[i][0] = 0.0;
+      vector[i][1] = 0.0;
+      matrix[i][i] = 1.0;
+    }      
   }
 
@@ -133,6 +143,6 @@
     ix = i % fit[0].Nterms;
     iy = i / fit[0].Nterms;
-    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
-	     ix, iy, vector[i][0], ix, iy, vector[i][1]);
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }	
 
@@ -142,6 +152,6 @@
     ix = i % fit[0].Nterms;
     iy = i / fit[0].Nterms;
-    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
-	     ix, iy, vector[i][0], ix, iy, vector[i][1]);
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }	
 
@@ -158,4 +168,28 @@
 }
 
+void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) {
+  
+  int ix, iy;
+  double xterm, yterm, term;
+  double x, y;
+
+  x = 0.0;
+  y = 0.0;
+
+  xterm = 1;
+  for (ix = 0; ix < fit[0].Nterms; ix++) {
+    yterm = 1;
+    for (iy = 0; iy < fit[0].Nterms; iy++) {
+      term = xterm*yterm;
+      x += fit[0].xfit[ix][iy]*term;
+      y += fit[0].yfit[ix][iy]*term;
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  *x2 = x;
+  *y2 = y;
+}
+
 // Nx, Ny is the number of terms in x and in y
 double **poly2d_dx (double **poly, int Nx, int Ny) {
@@ -209,5 +243,5 @@
   for (i = 0; i < Nx; i++) {
     for (j = 0; j < Ny; j++) {
-      out[i][j] = poly[i][j+1];
+      out[i][j] = poly[i][j];
     }
   }
@@ -237,9 +271,9 @@
 /* this should only apply to the polynomial, not the projection terms */
 /* compare with psastro supporting code */
-void fit_apply_coords (CoordFit *fit, Coords *coords) {
+CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) {
 
   double c11, c12;
   double c21, c22;
-  double Xo, Yo, R;
+  double Xo, Yo, R1, R2;
   CoordFit *modfit;
 
@@ -256,48 +290,41 @@
   modfit = CoordsSetCenter (fit, Xo, Yo);
 
+  /* we do not modify crval1,2: these are kept at the default values */
+
+  // set cdelt1, cdelt2
+  coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]);
+  coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]);
+  R1 = 1 / coords[0].cdelt1;
+  R2 = 1 / coords[0].cdelt2;
+
   // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
-  coords[0].pc1_1 = modfit[0].xfit[1][0];
-  coords[0].pc1_2 = modfit[0].xfit[0][1];
-  coords[0].pc2_1 = modfit[0].yfit[1][0];
-  coords[0].pc2_2 = modfit[0].yfit[0][1];
+  coords[0].pc1_1 = modfit[0].xfit[1][0] * R1;
+  coords[0].pc1_2 = modfit[0].xfit[0][1] * R2;
+  coords[0].pc2_1 = modfit[0].yfit[1][0] * R1;
+  coords[0].pc2_2 = modfit[0].yfit[0][1] * R2;
 
   // set the polyterm elements 
   if (coords->Npolyterms > 1) {
-    coords[0].polyterms[0][0] = modfit[0].xfit[2][0];
-    coords[0].polyterms[1][0] = modfit[0].xfit[1][1];
-    coords[0].polyterms[2][0] = modfit[0].xfit[0][2];
-
-    coords[0].polyterms[0][1] = modfit[0].yfit[2][0];
-    coords[0].polyterms[1][1] = modfit[0].yfit[1][1];
-    coords[0].polyterms[2][1] = modfit[0].yfit[0][2];
+    coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1;
+    coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2;
+    coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2;
+
+    coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1;
+    coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2;
+    coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2;
   }
 
   // I need to validate Norder
   if (coords->Npolyterms > 2) {
-    coords[0].polyterms[3][0] = modfit[0].xfit[3][0];
-    coords[0].polyterms[4][0] = modfit[0].xfit[2][1];
-    coords[0].polyterms[5][0] = modfit[0].xfit[1][2];
-    coords[0].polyterms[6][0] = modfit[0].xfit[0][3];
-
-    coords[0].polyterms[3][1] = modfit[0].yfit[3][0];
-    coords[0].polyterms[4][1] = modfit[0].yfit[2][1];
-    coords[0].polyterms[5][1] = modfit[0].yfit[1][2];
-    coords[0].polyterms[6][1] = modfit[0].yfit[0][3];
-  }
-
-  /* we do not modify crval1,2: these are kept at the default values */
-
-  /* normalize pc11,etc */
-  c11 = coords[0].pc1_1;
-  c21 = coords[0].pc1_1;
-  c12 = coords[0].pc1_1;
-  c22 = coords[0].pc1_1;
-  coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
-  R = 1 / coords[0].cdelt1;
-
-  coords[0].pc1_1  = c11*R;
-  coords[0].pc2_1  = c21*R;
-  coords[0].pc1_2  = c12*R;
-  coords[0].pc2_2  = c22*R;
+    coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2;
+
+    coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2;
+  }
 
   fit_free (modfit);
Index: trunk/Ohana/src/relastro/src/load_catalogs.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_catalogs.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/load_catalogs.c	(revision 16810)
@@ -39,6 +39,6 @@
 	for (j = 0; j < catalog[i].Naverage; j++) {
 	  catalog[i].average[j].code = 0;
-	  m = catalog[i].average[j].offset;
-	  for (k = 0; k < catalog[i].average[j].Nm; k++) {
+	  m = catalog[i].average[j].measureOffset;
+	  for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
 	    catalog[i].measure[m+k].dbFlags = 0;
 	  }
Index: trunk/Ohana/src/relastro/src/mkpolyterm.c
===================================================================
--- trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 16810)
@@ -20,5 +20,5 @@
 
   int i, Nx, Ny;
-  double R, Xo, Yo, dPos;
+  double R, Xo, Yo, dPos, dPosRef;
   double **XdX, **XdY, **YdX, **YdY, **alpha, **beta;
   double **xfit, **yfit;
@@ -51,4 +51,5 @@
      */
     dPos = tol + 1;
+    dPosRef = dPos;
     for (i = 0; (dPos > tol) && (i < 20); i++) {
       // NOTE: order for alpha is: [y][x]
@@ -66,4 +67,10 @@
       Yo -= beta[1][0];
       dPos = hypot(beta[0][0], beta[1][0]);
+      if (i == 0) {
+	dPosRef = dPos;
+      }
+    }
+    if (dPos > dPosRef) {
+      fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n");
     }
     array_free (alpha, 2);
Index: trunk/Ohana/src/relastro/src/plot_scatter.c
===================================================================
--- trunk/Ohana/src/relastro/src/plot_scatter.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/plot_scatter.c	(revision 16810)
@@ -12,5 +12,5 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Ntot += catalog[i].average[j].Nm;
+      Ntot += catalog[i].average[j].Nmeasure;
     }
   }
@@ -25,7 +25,7 @@
       /* calculate the average value for a single star */
       if (catalog[i].average[j].code & STAR_BAD) continue;  
-      m = catalog[i].average[j].offset;
+      m = catalog[i].average[j].measureOffset;
 
-      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
 
Index: trunk/Ohana/src/relastro/src/relastro.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/relastro.c	(revision 16810)
@@ -64,4 +64,20 @@
     save_catalogs (catalog, Ncatalog);
   } else {
+    // XXX for now, reload all catalogs at once; as memory use gets large, reload one-at-a-time
+
+    // load catalog data from region files : do not subselect
+    catalog = load_catalogs (skylist, &Ncatalog, FALSE);
+    
+    // match measurements with images
+    initImageBins (catalog, Ncatalog);
+    findImages (catalog, Ncatalog);
+
+    // update the detection coordinates using the new image parameters
+    UpdateMeasures (catalog, Ncatalog);
+
+    // write the updated detections to disk
+    save_catalogs (catalog, Ncatalog);
+
+    // save the updated image parameters
     dvo_image_update (&db, VERBOSE);
     dvo_image_unlock (&db); 
Index: trunk/Ohana/src/relastro/src/setExclusions.c
===================================================================
--- trunk/Ohana/src/relastro/src/setExclusions.c	(revision 16140)
+++ trunk/Ohana/src/relastro/src/setExclusions.c	(revision 16810)
@@ -11,6 +11,6 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      m = catalog[i].average[j].offset;
-      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+      m = catalog[i].average[j].measureOffset;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 
 	/* select measurements by photcode */