Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 27551)
+++ trunk/Ohana/src/relastro/Makefile	(revision 27581)
@@ -56,4 +56,6 @@
 $(SRC)/save_catalogs.$(ARCH).o       \
 $(SRC)/write_coords.$(ARCH).o        \
+$(SRC)/CoordOps.$(ARCH).o        \
+$(SRC)/FixProblemImages.$(ARCH).o        \
 $(SRC)/relastroVisual.$(ARCH).o
 
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 27551)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 27581)
@@ -91,6 +91,8 @@
 
 int    VERBOSE;
+int    VERBOSE2;
 
 int    RESET;
+int    NLOOP;
 int    UPDATE;
 int    PLOTSTUFF;
@@ -103,4 +105,7 @@
 int    PLOTDELAY;
 int    CHIPORDER;
+
+int UserCatalog;
+double UserCatalogRA, UserCatalogDEC;
 
 char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
@@ -151,5 +156,5 @@
 void          dump_grid           PROTO((void));
 int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
-void          findImages          PROTO((Catalog *catalog, int Ncatalog));
+void          findImages          PROTO((Catalog *catalog, int Ncatalog, int MATCHCAT));
 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
 Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
@@ -166,5 +171,5 @@
 void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
 int           gcatalog            PROTO((Catalog *catalog, int FINAL));
-Coords       *getCoords           PROTO((off_t meas, int cat));
+// Coords       *getCoords           PROTO((off_t meas, int cat));
 float         getMcal             PROTO((off_t meas, int cat));
 float         getMgrid            PROTO((off_t meas, int cat));
@@ -178,5 +183,5 @@
 void          initGrid            PROTO((int dX, int dY));
 void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
-void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
+void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int FULLINIT));
 void          initImages          PROTO((Image *input, off_t N));
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
@@ -199,5 +204,5 @@
 int           main                PROTO((int argc, char **argv));
 void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
-void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
+void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat, int MATCHCAT));
 void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
 GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
@@ -250,5 +255,5 @@
 void fit_free (CoordFit *fit);
 void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
-void fit_eval (CoordFit *fit);
+int 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);
@@ -256,8 +261,8 @@
 double **poly2d_copy (double **poly, int Nx, int Ny);
 double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
-CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords);
+int 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);
-void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
 void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
 void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
@@ -273,5 +278,5 @@
 
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
-int ParFactor (double *pR, double *pD, double R, double D, time_t T);
+int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean);
 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
@@ -303,2 +308,15 @@
 			       StatType statsR, StatType statsD, double thresh);
 
+
+
+int FixProblemImages (SkyList *skylist);
+int *getCatlist (int *N, off_t im);
+
+void initCoords (void);
+void getOffsets (double *dPos, off_t *nPos, off_t N);
+void saveOffsets (double dPos, off_t nPos, off_t N);
+void setBadCoords (off_t N);
+int badCoords (off_t N);
+Coords *getCoords (off_t N);
+int saveCoords (Coords *coords, off_t N);
+void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);
Index: trunk/Ohana/src/relastro/src/CoordOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/CoordOps.c	(revision 27581)
+++ trunk/Ohana/src/relastro/src/CoordOps.c	(revision 27581)
@@ -0,0 +1,83 @@
+# include "relastro.h"
+
+static double *dPosSum    = NULL; // sum of dPos^2 for all measures on each image
+static off_t  *nPosSum    = NULL; // sum of measures on each image (used for dPosSum)
+static int    *isBadCoord = NULL; // keep or reject each image?
+
+static Coords *oldCoords;   // list of available images
+static off_t  NoldCoords;   // number of available images
+
+void initCoords (void) {
+
+  off_t N;
+  Image *images;
+
+  images = getimages (&N);
+
+  NoldCoords = N;
+  ALLOCATE (oldCoords, Coords, NoldCoords);
+  ALLOCATE (dPosSum, double, NoldCoords);
+  ALLOCATE (nPosSum, off_t,  NoldCoords);
+  ALLOCATE (isBadCoord, int, NoldCoords);
+  memset (oldCoords,  0, N*sizeof(Coords));
+  memset (dPosSum,    0, N*sizeof(double));
+  memset (nPosSum,    0, N*sizeof(off_t));
+  memset (isBadCoord, 0, N*sizeof(int));
+}
+
+int saveCoords (Coords *coords, off_t N) {
+
+  if (N < 0) return FALSE;
+  if (N >= NoldCoords) return FALSE;
+
+  memcpy (&oldCoords[N], coords, sizeof(Coords));
+  return TRUE;
+}
+
+Coords *getCoords (off_t N) {
+
+  if (N < 0) return NULL;
+  if (N >= NoldCoords) return NULL;
+
+  return (&oldCoords[N]);
+}
+
+int badCoords (off_t N) {
+
+  if (N < 0) return FALSE;
+  if (N >= NoldCoords) return FALSE;
+
+  return (isBadCoord[N]);
+}
+  
+void setBadCoords (off_t N) {
+
+  if (N < 0) return;
+  if (N >= NoldCoords) return;
+
+  isBadCoord[N] = TRUE;
+  return;
+}
+  
+void saveOffsets (double dPos, off_t nPos, off_t N) {
+
+  if (N < 0) return;
+  if (N >= NoldCoords) return;
+
+  dPosSum[N] += dPos;
+  nPosSum[N] += nPos;
+  
+  return;
+}
+  
+void getOffsets (double *dPos, off_t *nPos, off_t N) {
+
+  if (N < 0) return;
+  if (N >= NoldCoords) return;
+
+  *dPos = dPosSum[N];
+  *nPos = nPosSum[N];
+  
+  return;
+}
+  
Index: trunk/Ohana/src/relastro/src/FitChip.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitChip.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/FitChip.c	(revision 27581)
@@ -17,5 +17,5 @@
 // XXX save measurements of the fit quality (scatter, chisq) in the image table
 
-void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
 
   int i, Nscatter, Niter, skip;
@@ -84,5 +84,5 @@
       default:
         fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
-        abort ();
+	skip = TRUE;
     }
     if (skip) {
@@ -90,12 +90,20 @@
       fit_free (fit);
       free (values);
-      return;
+      return FALSE;
     }
 
-    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
+    // fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
 
     // measure the fit, update the coords & object coordinates
-    fit_eval (fit);
-    fit_apply_coords (fit, coords);
+    if (!fit_eval (fit)) {
+      fprintf (stderr, "failed to fit new model\n");
+      return FALSE;
+    }
+
+    if (!fit_apply_coords (fit, coords)) {
+      fprintf (stderr, "failed to fit new model\n");
+      return FALSE;
+    }
+
     fit_free (fit);
 
@@ -107,5 +115,5 @@
 
   free (values);
-  return;
+  return TRUE;
 }
 
Index: trunk/Ohana/src/relastro/src/FixProblemImages.c
===================================================================
--- trunk/Ohana/src/relastro/src/FixProblemImages.c	(revision 27581)
+++ trunk/Ohana/src/relastro/src/FixProblemImages.c	(revision 27581)
@@ -0,0 +1,83 @@
+# include "relastro.h"
+
+// loop over all images.  for any images which have a bad coordinate solution, replace the
+// original coordinates and recalculate the positions
+
+int FixProblemImages (SkyList *skylist) {
+
+  int Nbad;
+  off_t i, Nimage;
+  Image *image;
+  SkyList sublist;
+
+  // allocate so we can reallocate below
+  ALLOCATE (sublist.regions, SkyRegion *, 1);
+  ALLOCATE (sublist.filename, char *, 1);
+
+  image = getimages (&Nimage);
+
+  Nbad = 0;
+  // first check on the dPos reported for each image
+  for (i = 0; i < Nimage; i++) {
+    double dPosSum, dPos;
+    off_t nPos;
+
+    // check if this image should be fixed
+    if (badCoords(i)) {
+      Nbad ++;
+      continue;
+    }
+
+    getOffsets (&dPosSum, &nPos, i);
+    dPos = sqrt(dPosSum / nPos);
+    if (dPos > 4.0) {
+      setBadCoords (i);
+      Nbad ++;
+    }
+  }
+
+  fprintf (stderr, "fixing %d images\n", Nbad);
+
+  for (i = 0; i < Nimage; i++) {
+    int j, cat, Ncat, *catlist, Ncatlist;
+    Catalog *catalog;
+
+    // check if this image should be fixed
+    if (!badCoords(i)) continue;
+
+    fprintf (stderr, "fixing %s\n", image[i].name);
+
+    // I need a list of the catalogs for this image
+    catlist = getCatlist(&Ncatlist, i);
+
+    // allocate Ncatlist skylist regions
+    REALLOCATE (sublist.regions, SkyRegion *, Ncatlist);
+    REALLOCATE (sublist.filename, char *, Ncatlist);
+    sublist.Nregions = Ncatlist;
+    sublist.ownElements = FALSE; // this list is only holding a view to the elements
+
+    // copy the desired catalogs from skylist to skylistSubset
+    for (j = 0; j < Ncatlist; j++) {
+      cat = catlist[j];
+      sublist.filename[j] = skylist[0].filename[cat];
+      sublist.regions[j] = skylist[0].regions[cat];
+    }
+
+    catalog = load_catalogs (&sublist, &Ncat, FALSE);
+    assert (Ncat == Ncatlist);
+
+    // match measurements with images
+    initImageBins (catalog, Ncat, FALSE);
+    findImages (catalog, Ncat, FALSE);
+
+    // update the detection coordinates using the new image parameters
+    resetImageRaw (catalog, Ncat, i);
+
+    freeImageBins (Ncat);
+
+    // write the updated detections to disk
+    save_catalogs (catalog, Ncat);
+  }
+  
+  return (TRUE);
+}
Index: trunk/Ohana/src/relastro/src/GetAstromError.c
===================================================================
--- trunk/Ohana/src/relastro/src/GetAstromError.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/GetAstromError.c	(revision 27581)
@@ -8,11 +8,11 @@
   switch (mode) {
     case ERROR_MODE_RA:
-      dPobs = measure[0].dXccd;  // need to redefine this as RAerr
+      dPobs = measure[0].dXccd / 100.0;  // need to redefine this as RAerr
       break;
     case ERROR_MODE_DEC:
-      dPobs = measure[0].dYccd;  // need to redefine this as RAerr
+      dPobs = measure[0].dYccd / 100.0;  // 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
+      dPobs = hypot (measure[0].dXccd, measure[0].dYccd) / 100.0;  // need to redefine this as RAerr
       break;
     default:
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 27581)
@@ -11,4 +11,8 @@
 static Image        *image;   // list of available images
 static off_t        Nimage;   // number of available images
+
+static int         *Ncatlist;  // catalogs associated with each image
+static int         *NCATLIST;  // catalogs associated with each image
+static int         **catlist;  // catalogs associated with each image
 
 // if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
@@ -32,4 +36,10 @@
 }
 
+int *getCatlist (int *N, off_t im) {
+
+  *N = Ncatlist[im];
+  return (catlist[im]);
+}
+
 void initImages (Image *input, off_t N) {
 
@@ -86,5 +96,6 @@
 }
 
-void initImageBins (Catalog *catalog, int Ncatalog) {
+// these are really image & catalog indexes
+void initImageBins (Catalog *catalog, int Ncatalog, int FULLINIT) {
 
   off_t i, j;
@@ -102,8 +113,20 @@
 
   for (i = 0; i < Nimage; i++) {
-    Nlist[i] = 0;
+    Nlist[i] =   0;
     NLIST[i] = 100;
     ALLOCATE (clist[i], int, NLIST[i]);
     ALLOCATE (mlist[i], off_t, NLIST[i]);
+  }
+
+  if (FULLINIT) {
+    ALLOCATE (Ncatlist, int,  Nimage);
+    ALLOCATE (NCATLIST, int,  Nimage);
+    ALLOCATE (catlist, int *, Nimage);
+
+    for (i = 0; i < Nimage; i++) {
+      Ncatlist[i] =  0;
+      NCATLIST[i] = 32;
+      ALLOCATE (catlist[i], int, NCATLIST[i]);
+    }
   }
 }
@@ -126,5 +149,5 @@
 
 /* match measurements to images */
-void findImages (Catalog *catalog, int Ncatalog) {
+void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) {
 
   off_t i, j;
@@ -136,9 +159,9 @@
       // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
       // if (photcode[0].code != ecode) continue;
-      matchImage (catalog, j, i);
-    }
-  }
-
-  for (i = 0; VERBOSE && (i < Nimage); i++) {
+      matchImage (catalog, j, i, MATCHCAT);
+    }
+  }
+
+  for (i = 0; VERBOSE2 && (i < Nimage); i++) {
     name = GetPhotcodeNamebyCode (image[i].photcode);
     fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
@@ -149,8 +172,9 @@
 # if USE_IMAGE_ID
 // this is the imageID-based match
-void matchImage (Catalog *catalog, off_t meas, int cat) {
+void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {
 
   off_t idx, ID;
   Measure *measure;
+  int i, found;
 
   measure = &catalog[cat].measure[meas];
@@ -159,6 +183,6 @@
   idx = getImageByID (ID);
   if (idx == -1) {
-    fprintf (stderr, "can't match detection to image?\n");
-    abort();
+    if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
+    return;
   }
 
@@ -178,4 +202,21 @@
     REALLOCATE (mlist[idx], off_t, NLIST[idx]);
   }
+
+  if (MATCHCAT) {
+    // index for image -> catalog list
+    found = FALSE;
+    for (i = 0; !found && (i < Ncatlist[idx]); i++) {
+      if (catlist[idx][i] == cat) found = TRUE;
+    }
+    if (!found) {
+      catlist[idx][Ncatlist[idx]] = cat;
+      Ncatlist[idx] ++;
+      if (Ncatlist[idx] == NCATLIST[idx]) {
+	NCATLIST[idx] += 32;
+	REALLOCATE (catlist[idx], int, NCATLIST[idx]);
+      }
+    }
+  }
+
   return;
 }
@@ -183,5 +224,5 @@
 # else
 // this is the time-based match
-void matchImage (Catalog *catalog, off_t meas, int cat) {
+void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {
 
   off_t i;
@@ -216,9 +257,10 @@
     return;
   }
-  fprintf (stderr, "can't match detection to image?\n");
-  abort();
+  if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
+  return;
 }
 # endif
 
+/*
 Coords *getCoords (off_t meas, int cat) {
 
@@ -229,4 +271,5 @@
   return (&image[i].coords);
 }
+*/
 
 void plot_images () {
@@ -280,16 +323,28 @@
 void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
 
-  off_t i, m, c, n;
+  off_t i, m, c, n, nPos;
   double X, Y, L, M, P, Q, R, D, dR, dD;
+  double dPos;
 
   Mosaic *mosaic;
   Coords *moscoords, *imcoords;
 
+  // check if this image is bad and should be skipped
+  if (badCoords(im)) return;
+
+  // WRP images need to have an associated mosaic
   moscoords = NULL;
-  mosaic = getMosaicForImage (im);
-  if (mosaic != NULL) {
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) return;  // if we cannot find the associated image, skip it
     moscoords = &mosaic[0].coords;
   }
   imcoords = &image[im].coords;
+
+  // accumulate the rms position offsets.  if this value, or any specific entry, is too
+  // large, we will reset the image to the original coords at the end of the analysis
+
+  dPos = 0.0;
+  nPos = 0;
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -316,14 +371,30 @@
     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????
+    // complain if the new location is far from the average location
+    if (fabs(dR) > 2.0) {
       fprintf (stderr, "!");
-      // abort ();
-    }
-    if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
+      setBadCoords (im); // report a failure for this image
+      return;
+    }
+    if (fabs(dD) > 2.0) {
       fprintf (stderr, "*");
-      // abort ();
-    }
+      setBadCoords (im); // report a failure for this image
+      return;
+    }
+
+    // complain if the new location is far from the old location
+    if (fabs(catalog[c].measure[m].dR - dR) > 2.0) {
+      fprintf (stderr, "@");
+      setBadCoords (im); // report a failure for this image
+      return;
+    }
+    if (fabs(catalog[c].measure[m].dD - dD) > 2.0) {
+      fprintf (stderr, "#");
+      setBadCoords (im); // report a failure for this image
+      return;
+    }
+
+    dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD);
+    nPos ++;
 
     catalog[c].measure[m].dR = dR;
@@ -341,4 +412,74 @@
     }
   }
+
+  saveOffsets (dPos, nPos, im);
+
+  return;
+}
+
+// return StarData values for detections in the specified image, converting coordinates from the
+// chip positions: X,Y -> L,M -> P,Q -> R,D
+void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
+
+  off_t i, m, c, n;
+  double X, Y, L, M, P, Q, R, D, dR, dD;
+
+  Mosaic *mosaic;
+  Coords *moscoords, *imcoords, *oldcoords;
+
+  // check if this image is bad and should be skipped
+  if (!badCoords(im)) {
+    fprintf (stderr, "ERROR: inconsistent result?");
+    abort();
+  }
+
+  // replace the current coords with the old coords:
+  oldcoords = getCoords (im);
+  memcpy (&image[im].coords, oldcoords, sizeof(Coords));
+
+  // WRP images need to have an associated mosaic
+  moscoords = NULL;
+  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) return;  // if we cannot find the associated image, skip it
+    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;
+
+    dR = dD = 0.0;
+    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);
+    }
+
+    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;
 }
@@ -354,6 +495,4 @@
   Coords *moscoords;
   StarData *raw;
-
-  ALLOCATE (raw, StarData, Nlist[im]);
 
   mosaic = NULL;
@@ -363,8 +502,10 @@
     if (mosaic == NULL) {
       fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
-      exit (1);
+      return NULL;
     }
     moscoords = &mosaic[0].coords;
   }
+
+  ALLOCATE (raw, StarData, Nlist[im]);
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -430,6 +571,4 @@
   StarData *ref;
 
-  ALLOCATE (ref, StarData, Nlist[im]);
-
   mosaic = NULL;
   moscoords = NULL;
@@ -438,8 +577,10 @@
     if (mosaic == NULL) {
       fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
-      exit (1);
+      return NULL;
     }
     moscoords = &mosaic[0].coords;
   }
+
+  ALLOCATE (ref, StarData, Nlist[im]);
 
   for (i = 0; i < Nlist[im]; i++) {
Index: trunk/Ohana/src/relastro/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 27581)
@@ -146,5 +146,8 @@
 
     Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
-    if (Nmos == -1) continue;
+    if (Nmos == -1) {
+      fprintf (stderr, "cannot match mosaic for %s\n", image[i].name);
+      continue;
+    }
 
     // mosaic corresponding to this image
@@ -185,4 +188,8 @@
     // this function does the reverse-lookup for the mosaic corresponding to this image
     new = getImageRaw (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    if (!new) {
+      fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (1)\n");
+      abort();
+    }
     
     // merge new and raw
@@ -217,4 +224,8 @@
     // this function does the reverse-lookup for the mosaic corresponding to this image
     new = getImageRef (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    if (!new) {
+      fprintf (stderr, "inconsistent: missing mosaic for image already associated with a mosaic? (2)\n");
+      abort();
+    }
     
     // merge new and ref
Index: trunk/Ohana/src/relastro/src/ParFactor.c
===================================================================
--- trunk/Ohana/src/relastro/src/ParFactor.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/ParFactor.c	(revision 27581)
@@ -42,5 +42,5 @@
 
 /* given RA, DEC, Time, calculate the parallax factor */
-int ParFactor (double *pR, double *pD, double R, double D, time_t T) {
+int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean) {
 
   double jd;
@@ -49,5 +49,5 @@
   /* given a time T in UNIX seconds, determine the solar longitude S */
 
-  jd = ohana_sec_to_jd (T);
+  jd = ohana_sec_to_jd (365.25*86400.0*(T + Tmean));
   sun_ecliptic (jd, &L, &B, &E);
 
Index: trunk/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 27581)
@@ -7,4 +7,5 @@
   Image *image;
   StarData *raw, *ref;
+  Coords *oldCoords;
 
   image = getimages (&Nimage);
@@ -17,14 +18,22 @@
     /* convert measure coordinates to raw entries */
     raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC);
+    if (!raw) continue;
 
     /* convert average coordinates to ref entries */
     ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC);
+    if (!ref) continue;
 
     // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
     assert (Nraw == Nref);
 
+    saveCoords (&image[i].coords, i);
+
     // FitChip does iterative, clipped fitting
-    fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw);
-    FitChip (raw, ref, Nraw, &image[i].coords);
+    // fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw);
+    if (!FitChip (raw, ref, Nraw, &image[i].coords)) {
+      fprintf (stderr, "reject fit for image %s (%lld) : Nstars: %lld\n", image[i].name, (long long) i, (long long) Nraw);
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+    }
 
     free (raw);
Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c	(revision 27581)
@@ -1,3 +1,9 @@
 # include "relastro.h"
+
+// We run through each DVO catalog, updating the measures that come from the modified images
+// We need to watch for failures:
+// * in UpdateMeasures, in fixImageRaw, we track the cumulative offset for each image
+// * after all updates are done, we can check for any bad images and reset them to the
+//   original coordinates
 
 int UpdateObjectOffsets (SkyList *skylist) {
@@ -28,6 +34,6 @@
 
     // match measurements with images
-    initImageBins (&catalog, 1);
-    findImages (&catalog, 1);
+    initImageBins (&catalog, 1, FALSE);
+    findImages (&catalog, 1, FALSE);
 
     // update the detection coordinates using the new image parameters
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 27581)
@@ -40,11 +40,13 @@
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
-  int i, j, k, m, N, Nsecfilt;
+  off_t j, k, m;
+  int i, N, Nsecfilt, mode, result, status, XVERB;
   StatType statsR, statsD;
   Coords coords;
-  PMFit fit;
-  time_t To;
-  int mode, Nave, Npm, Npar, Nskip;
-  double Tmin, Tmax;
+  PMFit fitAve, fitPM, fitPAR, fit;
+  time_t T2000;
+  off_t Nave, Npm, Npar, Nskip;
+  off_t NaveSum, NpmSum, NparSum, NskipSum;
+  double Tmin, Tmax, Tmean, Trange;
 
   initObjectData (catalog, Ncatalog);
@@ -61,7 +63,8 @@
   strcpy (coords.ctype, "RA---SIN");
 
+  XVERB = FALSE;
+
   // use J2000 as a reference time
-  To = ohana_date_to_sec ("2000/01/01");
-  Nave = Npar = Npm = 0;
+  T2000 = ohana_date_to_sec ("2000/01/01");
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -69,9 +72,10 @@
   assert (catalog[0].Nsecfilt == Nsecfilt);
 
+  NaveSum = NparSum = NpmSum = NskipSum = 0;
   for (i = 0; i < Ncatalog; i++) {
 
     if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure);
 
-    Nskip = 0;
+    Nave = Npar = Npm = Nskip = 0;
     for (j = 0; j < catalog[i].Naverage; j++) {
       /* calculate the average value of R,D for a single star */
@@ -79,17 +83,19 @@
       // skip objects which are known to be problematic
       // XXX include this code or not?
-      # if (0)
+# if (0)
       if (catalog[i].average[j].code & STAR_BAD) {
 	Nskip ++;
 	continue;  
       }
-      # endif
+# endif
 
       N = 0;
       m = catalog[i].average[j].measureOffset;
 
-      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
+      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
       mode = FIT_MODE;
 
+      // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
+      Tmean = 0;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
 
@@ -108,10 +114,10 @@
 	// exclude measurements by previous outlier detection
 	// XXX include this code or not?
-	# if (0)
+# if (0)
 	if (catalog[i].measure[m].dbFlags & MEAS_BAD) { 
 	  catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
 	  continue;
 	}
-	# endif
+# endif
 
 	catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
@@ -119,12 +125,21 @@
 	R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
 	D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-	T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years
+	T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
 
 	Tmin = MIN(Tmin, T[N]);
 	Tmax = MAX(Tmax, T[N]);
-
-	dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
-	dD[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+	Tmean += T[N];
+
+	// dX, dY : error in arcsec -- 
+	dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
+	dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
 	dT[N] = catalog[i].measure[m].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[k] = dX[k] / 3600.0;
+	dD[k] = dY[k] / 3600.0;
 
 	N++;
@@ -138,5 +153,6 @@
 
       // XXX add the parallax factor range as a criterion as well
-      if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE;
+      Trange = Tmax - Tmin;
+      if (Trange < PM_DT_MIN) mode = FIT_AVERAGE;
       if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
 
@@ -144,5 +160,9 @@
       if (N < SRC_MEAS_TOOFEW) {
 	// XXX need to define PHOTOM and ASTROM object flags
+	// XXX reset the average value fields?
 	catalog[i].average[j].flags |= ID_STAR_FEW;
+	catalog[i].average[j].ChiSqAve  = NAN;
+	catalog[i].average[j].ChiSqPM   = NAN;
+	catalog[i].average[j].ChiSqPar  = NAN;
 	if (N < 2) continue;
       }
@@ -151,103 +171,112 @@
       coords.crval1 = R[0];
       coords.crval2 = D[0];
+      Tmean /= (float) N;
       
-      /* project all of the R,D coordinates to a plane centered on this coordinate */
-      for (k = 0; k < N; k++) {
-	RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
-	dX[k] =  dR[k];
-	dY[k] =  dD[k];
-	// fprintf (stderr, "%d %f %f %f  %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]);
+      XVERB = FALSE && (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
+
+      // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
+
+      if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
+	// project all of the R,D coordinates to a plane centered on this coordinate set
+	// the times to be relative to Tmean (this is required for parallax as well)
+	for (k = 0; k < N; k++) {
+	  RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
+	  T[k] -= Tmean;
+	  if (XVERB) {
+	    fprintf (stderr, "%lld %f %f %f  %f %f +/- %f %f\n", (long long) k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
+	  }
+	}	  
+
+	FitPM (&fitPM, X, dX, Y, dY, T, N);
+	if (XVERB) fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
+
+	// project Ro, Do back to RA,DEC
+	XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
+	if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
+
+	fitPM.p  = fitPM.dp  = 0.0;
+	catalog[i].average[j].flags |= ID_STAR_FIT_PM;
+	Npm ++;
+      }
+
+      if (mode == FIT_PM_AND_PAR) {
+	fprintf (stderr, "parallax fitting is still untested (%s, %d)", __FILE__, __LINE__);
+	exit (2);
+
+	for (k = 0; k < N; k++) {
+	  ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
+	}
+	FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N);
+	XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+	catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
+	Npar ++;
       }	  
 
-      /* fit the model components as needed */
-      switch (mode) {
+      // fit the average model
+      if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
+	liststats (R, dR, N, &statsR); // WARNING: this function modifies R (do not use after here)
+	liststats (D, dD, N, &statsD); // WARNING: this function modifies D (do not use after here)
+
+	fitAve.Ro = statsR.mean;
+	fitAve.dRo = 3600.0*statsR.sigma;
+
+	fitAve.Do = statsD.mean;
+	fitAve.dDo = 3600.0*statsD.sigma;
+
+	fitAve.chisq = 0.5*(statsR.chisq + statsD.chisq);
+	fitAve.Nfit = N;
+
+	fitAve.uR = fitAve.duR = 0.0;
+	fitAve.uD = fitAve.duD = 0.0;
+	fitAve.p  = fitAve.dp  = 0.0;
+	catalog[i].average[j].flags |= ID_STAR_FIT_AVE;
+	Nave ++;
+      }
+
+      /* choose the result based on the chisq values */
+      // XXXX for now, just use the mode as the result:
+      result = mode;
+
+      switch (result) {
 	case FIT_AVERAGE:
-	  liststats (R, dR, N, &statsR);
-	  liststats (D, dD, N, &statsD);
-
-	  fit.Ro = statsR.mean;
-	  fit.dRo = 3600.0*statsR.sigma;
-
-	  fit.Do = statsD.mean;
-	  fit.dDo = 3600.0*statsD.sigma;
-
-	  fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
-	  fit.Nfit = N;
-
-	  fit.uR = fit.duR = 0.0;
-	  fit.uD = fit.duD = 0.0;
-	  fit.p  = fit.dp  = 0.0;
-
-	  Nave ++;
+	  catalog[i].average[j].flags |= ID_STAR_USE_AVE;
+	  fit = fitAve;
 	  break;
-
 	case FIT_PM_ONLY:
-	  FitPM (&fit, X, dX, Y, dY, T, N);
-	  // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
-	  // project Ro, Do back to RA,DEC
-	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-	  // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
-	  // continue;
-
-	  fit.p  = fit.dp  = 0.0;
-
-	  Npm ++;
+	  catalog[i].average[j].flags |= ID_STAR_USE_PM;
+	  fit = fitPM;
 	  break;
-
-	case FIT_PAR_ONLY:
-	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
-	  exit (2);
-
-	  for (k = 0; k < N; k++) {
-	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
-	  }
-	  FitPar (&fit, X, dX, Y, dY, pX, pY, N);
-
-	  // project Ro, Do back to RA,DEC
-	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-
-	  fit.uR = fit.duR = 0.0;
-	  fit.uD = fit.duD = 0.0;
-
-	  Npar ++;
+	case FIT_PM_AND_PAR:
+	  catalog[i].average[j].flags |= ID_STAR_USE_PAR;
+	  fit = fitPAR;
 	  break;
-
-	case FIT_PM_AND_PAR:
-	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
-	  exit (2);
-
-	  for (k = 0; k < N; k++) {
-	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
-	  }
-	  FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
-	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
-	  Npar ++;
-	  break;
-
-	default:
-	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
-	  exit (2);
-      }	  
-
-      if (0 && (j < 100)) {
-	fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
-		 catalog[i].average[j].R, 
-		 catalog[i].average[j].D, 
-		 fit.Ro, fit.Do, 
-		 3600*(catalog[i].average[j].R - fit.Ro), 
-		 3600*(catalog[i].average[j].D - fit.Do));
-      }
+      }
+
+      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
+			  catalog[i].average[j].R, 
+			  catalog[i].average[j].D, 
+			  fit.Ro, fit.Do, 
+			  3600*(catalog[i].average[j].R - fit.Ro), 
+			  3600*(catalog[i].average[j].D - fit.Do));
 
       //make sure that the fit succeeded
-      assert(finite(fit.Ro) && finite(fit.Do) && 
-	     finite(fit.dRo) && finite(fit.dDo) &&
-	     finite(fit.uR) && finite(fit.uD) &&
-	     finite(fit.duR) && finite(fit.duD) &&
-	     finite(fit.p) && finite(fit.dp));
+      status  = finite(fit.Ro);
+      status &= finite(fit.Do);
+      status &= finite(fit.dRo);
+      status &= finite(fit.dDo);
+      status &= finite(fit.uR);
+      status &= finite(fit.uD);
+      status &= finite(fit.duR);
+      status &= finite(fit.duD);
+      status &= finite(fit.p);
+      status &= finite(fit.dp);
+      if (!status) {
+	Nskip ++;
+	continue;
+      }
 
       // the measure fields must be updated before the average fields
       m = catalog[i].average[j].measureOffset;
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	// XXX why was this here?? if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
 	setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
 	setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
@@ -267,11 +296,23 @@
       catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec
 
-      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
+      // Xp is supposed to be the position scatter, not the chisq : fix this:
+      // catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
+      catalog[i].average[j].ChiSqAve  = fitAve.chisq;
+      catalog[i].average[j].ChiSqPM   = fitPM.chisq;
+      catalog[i].average[j].ChiSqPar  = fitPAR.chisq;
+      catalog[i].average[j].Xp        = 0.0;
+      catalog[i].average[j].Tmean = (Tmean * 86400 * 365.26) + T2000;
+      catalog[i].average[j].Trange = (Trange * 86400 * 365.26);
+      catalog[i].average[j].Npos = fit.Nfit;
     }
 
-    if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip);
+    NaveSum += Nave;
+    NpmSum += Npm;
+    NparSum += Npar;
+    NskipSum += Nskip;
+    if (VERBOSE) fprintf (stderr, "catalog %lld : %lld ave, %lld pm, %lld par : Nskip %lld\n", (long long) i, (long long) Nave, (long long) Npm, (long long) Npar, (long long) Nskip);
   }
 
-  if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar);
+  fprintf (stderr, "fitted %lld objects (%lld ave, %lld pm, %lld par), skipped %lld\n", (long long) (NaveSum + NpmSum + NparSum), (long long) NaveSum, (long long) NpmSum, (long long) NparSum, (long long) NskipSum);
   return (TRUE);
 }
@@ -279,25 +320,25 @@
 /* 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?
+   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?
 
 */
Index: trunk/Ohana/src/relastro/src/UpdateSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 27581)
@@ -18,7 +18,9 @@
     /* convert measure coordinates to raw entries */
     raw = getImageRaw (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE);
+    if (!raw) continue;
 
     /* convert average coordinates to ref entries */
     ref = getImageRef (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE);
+    if (!ref) continue;
 
     FitSimple (raw, ref, Nstars, &image[i].coords);
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 27581)
@@ -50,4 +50,14 @@
 
   if (FIT_TARGET == TARGET_NONE) usage();
+
+  UserCatalog = FALSE;
+  if ((N = get_argument (argc, argv, "-catalog"))) {
+    remove_argument (N, &argc, argv);
+    UserCatalogRA = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    UserCatalogDEC = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    UserCatalog = TRUE;
+  }
 
   /* specify portion of the sky : allow default of all sky? */
@@ -67,5 +77,7 @@
     remove_argument (N, &argc, argv);
   } else {
-    usage ();
+    if (!UserCatalog) {
+      usage ();
+    }
   }
 
@@ -111,7 +123,11 @@
   }
 
-  VERBOSE = FALSE;
+  VERBOSE = VERBOSE2 = FALSE;
   if ((N = get_argument (argc, argv, "-v"))) {
     VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-vv"))) {
+    VERBOSE = VERBOSE2 = TRUE;
     remove_argument (N, &argc, argv);
   }
@@ -242,4 +258,11 @@
   }
 
+  NLOOP = 4;
+  if ((N = get_argument (argc, argv, "-nloop"))) {
+    remove_argument (N, &argc, argv);
+    NLOOP = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 1) usage ();
   return TRUE;
@@ -248,4 +271,5 @@
 void usage () {
   fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
+  fprintf (stderr, "       OR:    relastro -catalog (ra) (dec)\n");
   fprintf (stderr, "  working options: \n");
   fprintf (stderr, "  -update-objects\n");
@@ -263,4 +287,5 @@
   fprintf (stderr, "  -statmode (mode)\n");
   fprintf (stderr, "  -reset");
+  fprintf (stderr, "  -nloop (N) : number of image-fit iterations");
   fprintf (stderr, "  -update : apply new fit to database\n");
   fprintf (stderr, "  -params\n");
Index: trunk/Ohana/src/relastro/src/fitpoly.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpoly.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/fitpoly.c	(revision 27581)
@@ -103,5 +103,5 @@
 
 /* convert the xsum,ysum,sum terms into vector,matrix and solve */
-void fit_eval (CoordFit *fit) {
+int fit_eval (CoordFit *fit) {
 
   int i, j, ix, iy, jx, jy;
@@ -146,6 +146,8 @@
     // ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }	
-
-  dgaussjordan (matrix, vector, fit[0].Nelems, 2); 
+  
+  if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) {
+    return (FALSE);
+  }
 
   for (i = 0; i < fit[0].Nelems; i++) {
@@ -166,4 +168,5 @@
   array_free (matrix, fit[0].Nelems);
   array_free (vector, fit[0].Nelems);
+  return (TRUE);
 }
 
@@ -271,5 +274,5 @@
 /* this should only apply to the polynomial, not the projection terms */
 /* compare with psastro supporting code */
-CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) {
+int fit_apply_coords (CoordFit *fit, Coords *coords) {
 
   double Xo, Yo, R1, R2;
@@ -281,5 +284,8 @@
   // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
 
-  CoordsGetCenter (fit, 0.001, &Xo, &Yo);
+  if (!CoordsGetCenter (fit, 0.001, &Xo, &Yo)) {
+    fprintf (stderr, "failed to modify model\n");
+    return (FALSE);
+  }
   coords[0].crpix1 = Xo;
   coords[0].crpix2 = Yo;
@@ -329,6 +335,5 @@
   /* keep the order and type from initial values */
   
-  // XXX if desired in the future, return modfit (and free above)
-  return (NULL);
-}
-
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/load_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_images.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/load_images.c	(revision 27581)
@@ -24,5 +24,9 @@
   
   // determine the populated SkyRegions overlapping the requested area
-  skylist = SkyListByPatch (sky, -1, region);
+  if (UserCatalog) {
+    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
+  } else {
+    skylist = SkyListByPatch (sky, -1, region);
+  }
   MARKTIME("  setup sky: %f sec\n", dtime);
 
Index: trunk/Ohana/src/relastro/src/mkpolyterm.c
===================================================================
--- trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 27581)
@@ -73,4 +73,5 @@
     if (dPos > dPosRef) {
       fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n");
+      return FALSE;
     }
     array_free (alpha, 2);
Index: trunk/Ohana/src/relastro/src/relastro.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/relastro.c	(revision 27581)
@@ -9,5 +9,5 @@
 int main (int argc, char **argv) {
 
-  int status, Ncatalog;
+  int i, status, Ncatalog;
   Catalog *catalog;
   FITS_DB db;
@@ -40,4 +40,6 @@
   MARKTIME("load images: %f sec\n", dtime);
 
+  initCoords();
+
   /* load catalog data from region files : subselect high-quality measurements */
   catalog = load_catalogs (skylist, &Ncatalog, TRUE);
@@ -45,8 +47,8 @@
 
   /* match measurements with images */
-  initImageBins (catalog, Ncatalog);
+  initImageBins (catalog, Ncatalog, TRUE);
   MARKTIME("make image bins: %f sec\n", dtime);
 
-  findImages (catalog, Ncatalog);
+  findImages (catalog, Ncatalog, TRUE);
   MARKTIME("set up image indexes: %f sec\n", dtime);
 
@@ -59,13 +61,23 @@
   switch (FIT_TARGET) {
     case TARGET_SIMPLE:
-      UpdateSimple (catalog, Ncatalog);
+      for (i = 0; i < NLOOP; i++) {
+	UpdateObjects (catalog, Ncatalog);
+	UpdateSimple (catalog, Ncatalog);
+      }
       break;
 
     case TARGET_CHIPS:
-      UpdateChips (catalog, Ncatalog);
+      for (i = 0; i < NLOOP; i++) {
+	UpdateObjects (catalog, Ncatalog);
+	UpdateChips (catalog, Ncatalog);
+	MARKTIME("update chips: %f sec\n", dtime);
+      }
       break;
 
     case TARGET_MOSAICS:
-      UpdateMosaic (catalog, Ncatalog);
+      for (i = 0; i < NLOOP; i++) {
+	UpdateObjects (catalog, Ncatalog);
+	UpdateMosaic (catalog, Ncatalog);
+      }
       break;
 
@@ -83,4 +95,7 @@
   UpdateObjectOffsets (skylist);
 
+  // iterate over catalogs to make detection coordinates consistant
+  FixProblemImages (skylist);
+
   // save the updated image parameters
   dvo_image_update (&db, VERBOSE);
Index: trunk/Ohana/src/relastro/src/relastro_objects.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro_objects.c	(revision 27551)
+++ trunk/Ohana/src/relastro/src/relastro_objects.c	(revision 27581)
@@ -15,5 +15,9 @@
   
   // determine the populated SkyRegions overlapping the requested area (default depth)
-  skylist = SkyListByPatch (sky, -1, &UserPatch);
+  if (UserCatalog) {
+    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
+  } else {
+    skylist = SkyListByPatch (sky, -1, &UserPatch);
+  }
 
   // load data from each region file, only use bright stars