Changeset 27581

Timestamp:

Apr 2, 2010, 3:55:04 PM (16 years ago)

Author:

eugene

Message:

various improvements to relastro to calculate the significance of the astrometry model (merged dev work from branches/eam_brances/relastro.20100326)

Location:

trunk/Ohana/src/relastro

Files:

• . (modified) (1 prop)
• Makefile (modified) (1 diff)
• include/relastro.h (modified) (10 diffs)
• src/CoordOps.c (copied) (copied from branches/eam_branches/relastro.20100326/src/CoordOps.c )
• src/FitChip.c (modified) (4 diffs)
• src/FixProblemImages.c (copied) (copied from branches/eam_branches/relastro.20100326/src/FixProblemImages.c )
• src/GetAstromError.c (modified) (1 diff)
• src/ImageOps.c (modified) (19 diffs)
• src/MosaicOps.c (modified) (3 diffs)
• src/ParFactor.c (modified) (2 diffs)
• src/UpdateChips.c (modified) (2 diffs)
• src/UpdateObjectOffsets.c (modified) (2 diffs)
• src/UpdateObjects.c (modified) (11 diffs)
• src/UpdateSimple.c (modified) (1 diff)
• src/args.c (modified) (6 diffs)
• src/fitpoly.c (modified) (6 diffs)
• src/load_images.c (modified) (1 diff)
• src/mkpolyterm.c (modified) (1 diff)
• src/relastro.c (modified) (5 diffs)
• src/relastro_objects.c (modified) (1 diff)

trunk/Ohana/src/relastro/Makefile

@@ -56 +56 @@
 $(SRC)/save_catalogs.$(ARCH).o       \
 $(SRC)/write_coords.$(ARCH).o        \
+$(SRC)/CoordOps.$(ARCH).o        \
+$(SRC)/FixProblemImages.$(ARCH).o        \
 $(SRC)/relastroVisual.$(ARCH).o
 

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

@@ -91 +91 @@
 
 int    VERBOSE;
+int    VERBOSE2;
 
 int    RESET;
+int    NLOOP;
 int    UPDATE;
 int    PLOTSTUFF;
@@ -103 +105 @@
 int    PLOTDELAY;
 int    CHIPORDER;
+
+int UserCatalog;
+double UserCatalogRA, UserCatalogDEC;
 
 char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
@@ -151 +156 @@
 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 +171 @@
 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 +183 @@
 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 +204 @@
 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 +255 @@
 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 +261 @@
 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 +278 @@
 
 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 +308 @@
                                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);

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

@@ -17 +17 @@
 // 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 +84 @@
       default:
         fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
-        abort ();
+        skip = TRUE;
     }
     if (skip) {
@@ -90 +90 @@
       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 +115 @@
 
   free (values);
-  return;
+  return TRUE;
 }
 

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

@@ -8 +8 @@
   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:

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

@@ -11 +11 @@
 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 +36 @@
 }
 
+int *getCatlist (int *N, off_t im) {
+
+  *N = Ncatlist[im];
+  return (catlist[im]);
+}
+
 void initImages (Image *input, off_t N) {
 
@@ -86 +96 @@
 }
 
-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 +113 @@
 
   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 +149 @@
 
 /* match measurements to images */
-void findImages (Catalog *catalog, int Ncatalog) {
+void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) {
 
   off_t i, j;
@@ -136 +159 @@
       // 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 +172 @@
 # 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 +183 @@
   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 +202 @@
     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 +224 @@
 # 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 +257 @@
     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 +271 @@
   return (&image[i].coords);
 }
+*/
 
 void plot_images () {
@@ -280 +323 @@
 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 +371 @@
     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 +412 @@
     }
   }
+
+  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 +495 @@
   Coords *moscoords;
   StarData *raw;
-
-  ALLOCATE (raw, StarData, Nlist[im]);
 
   mosaic = NULL;
@@ -363 +502 @@
     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 +571 @@
   StarData *ref;
 
-  ALLOCATE (ref, StarData, Nlist[im]);
-
   mosaic = NULL;
   moscoords = NULL;
@@ -438 +577 @@
     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++) {

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

@@ -146 +146 @@
 
     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 +188 @@
     // 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 +224 @@
     // 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

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

@@ -42 +42 @@
 
 /* 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 +49 @@
   /* 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);
 

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

@@ -7 +7 @@
   Image *image;
   StarData *raw, *ref;
+  Coords *oldCoords;
 
   image = getimages (&Nimage);
@@ -17 +18 @@
     /* 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);

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

@@ -1 +1 @@
 # 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 +34 @@
 
     // 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

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

@@ -40 +40 @@
 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 +63 @@
   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 +72 @@
   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 +83 @@
       // 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 +114 @@
         // 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 +125 @@
         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 +153 @@
 
       // 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 +160 @@
       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 +171 @@
       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 +296 @@
       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 +320 @@
 /* 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?
 
 */

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

@@ -18 +18 @@
     /* 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);

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

@@ -50 +50 @@
 
   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 +77 @@
     remove_argument (N, &argc, argv);
   } else {
-    usage ();
+    if (!UserCatalog) {
+      usage ();
+    }
   }
 
@@ -111 +123 @@
   }
 
-  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 +258 @@
   }
 
+  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 +271 @@
 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 +287 @@
   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");

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

@@ -103 +103 @@
 
 /* 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 +146 @@
     // 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 +168 @@
   array_free (matrix, fit[0].Nelems);
   array_free (vector, fit[0].Nelems);
+  return (TRUE);
 }
 
@@ -271 +274 @@
 /* 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 +284 @@
   // 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 +335 @@
   /* keep the order and type from initial values */
   
-  // XXX if desired in the future, return modfit (and free above)
-  return (NULL);
-}
-
+  return (TRUE);
+}
+

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

@@ -24 +24 @@
   
   // 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);
 

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

@@ -73 +73 @@
     if (dPos > dPosRef) {
       fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n");
+      return FALSE;
     }
     array_free (alpha, 2);

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

@@ -9 +9 @@
 int main (int argc, char **argv) {
 
-  int status, Ncatalog;
+  int i, status, Ncatalog;
   Catalog *catalog;
   FITS_DB db;
@@ -40 +40 @@
   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 +47 @@
 
   /* 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 +61 @@
   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 +95 @@
   UpdateObjectOffsets (skylist);
 
+  // iterate over catalogs to make detection coordinates consistant
+  FixProblemImages (skylist);
+
   // save the updated image parameters
   dvo_image_update (&db, VERBOSE);

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

@@ -15 +15 @@
   
   // 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