Changeset 27840 for branches/simtest_nebulous_branches/Ohana/src/relastro

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/relastro (modified) (1 prop)
• Ohana/src/relastro/Makefile (modified) (1 diff)
• Ohana/src/relastro/doc/notes.txt (modified) (1 diff)
• Ohana/src/relastro/include/relastro.h (modified) (12 diffs)
• Ohana/src/relastro/src/ConfigInit.c (modified) (1 diff)
• Ohana/src/relastro/src/CoordOps.c (copied) (copied from trunk/Ohana/src/relastro/src/CoordOps.c )
• Ohana/src/relastro/src/FitChip.c (modified) (4 diffs)
• Ohana/src/relastro/src/FitSimple.c (modified) (1 diff)
• Ohana/src/relastro/src/FixProblemImages.c (copied) (copied from trunk/Ohana/src/relastro/src/FixProblemImages.c )
• Ohana/src/relastro/src/GetAstromError.c (modified) (1 diff)
• Ohana/src/relastro/src/ImageOps.c (modified) (19 diffs)
• Ohana/src/relastro/src/MosaicOps.c (modified) (10 diffs)
• Ohana/src/relastro/src/ParFactor.c (modified) (2 diffs)
• Ohana/src/relastro/src/UpdateChips.c (modified) (2 diffs)
• Ohana/src/relastro/src/UpdateMeasures.c (modified) (1 diff)
• Ohana/src/relastro/src/UpdateMosaic.c (modified) (1 diff)
• Ohana/src/relastro/src/UpdateObjectOffsets.c (modified) (2 diffs)
• Ohana/src/relastro/src/UpdateObjects.c (modified) (13 diffs)
• Ohana/src/relastro/src/UpdateSimple.c (modified) (2 diffs)
• Ohana/src/relastro/src/args.c (modified) (6 diffs)
• Ohana/src/relastro/src/bcatalog.c (modified) (2 diffs)
• Ohana/src/relastro/src/fitpoly.c (modified) (6 diffs)
• Ohana/src/relastro/src/initialize.c (modified) (1 diff)
• Ohana/src/relastro/src/load_images.c (modified) (1 diff)
• Ohana/src/relastro/src/mkpolyterm.c (modified) (1 diff)
• Ohana/src/relastro/src/plotstuff.c (modified) (1 diff)
• Ohana/src/relastro/src/relastro.c (modified) (4 diffs)
• Ohana/src/relastro/src/relastro_objects.c (modified) (1 diff)
• Ohana/src/relastro/src/select_images.c (modified) (11 diffs)

branches/simtest_nebulous_branches/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
 

branches/simtest_nebulous_branches/Ohana/src/relastro/doc/notes.txt

@@ +1 @@
+
+2010.03.24
+
+  I have finished only part of the optimizations listed below: I have
+  made the detection->image and image->mosaic lookups use bracket
+  searches.  I have NOT converted the image->mosaic lookups to use the
+  parentID, so these matches are still based on time ranges.  I have
+  also added a bit of positional constraint to the image match
+  function in 'select_images' (comparing Rmax of the sky regions to
+  Rmin of the images).  These have had a huge improvement.  Here is
+  the execution time for the SAS data (~450 images) with the old vs
+  the new code on pikake:
+
+  old : 39min
+  new :  7min
+
+  The bulk of the time for the new code is now disk I/O.
+
+2010.03.23
+
+  I've been working on optimization.  I've updated ImageOps to make
+  the measure->image relationship generation use the image IDs along
+  with a bracket search -- this goes much faster as a result.  More
+  work is needed:
+
+  * select_images needs to apply astrometry to all image corners: 
+    -> cache the image center and radius and use this to narrow down
+       the searches.
+
+  * MosaicOps uses the time to match mosaics.  I've added sorting to
+    speed this up, but this should be done based on the image and
+    parent index.
+
+    -> make a tool to create the parent ID for existing DBs.
+    -> make sure addstar is populating the parent IDs
+    -> use the parent IDs to make the link.
+
+load image data: 0.006347 sec
+  setup sky: 0.091475 sec
+  convert image table: 0.091498 sec
+  select images: 2.925978 sec
+  init images: 2.926957 sec
+  init mosaics: 2.927319 sec
+load images: 2.933762 sec
+load catalog data: 9.772259 sec
+make image bins: 9.780618 sec
+set up image indexes: 10.196519 sec
+
+load image data: 0.006746 sec
+  setup sky: 0.091440 sec
+  convert image table: 0.091467 sec
+  select images: 0.374193 sec
+  init images: 0.375355 sec
+  init mosaics: 0.375753 sec
+load images: 0.382606 sec
+load catalog data: 7.260094 sec
+make image bins: 7.268564 sec
+set up image indexes: 7.687485 sec
+
 
 2008.03.01

branches/simtest_nebulous_branches/Ohana/src/relastro/include/relastro.h

@@ -91 +91 @@
 
 int    VERBOSE;
+int    VERBOSE2;
 
 int    RESET;
+int    NLOOP;
 int    UPDATE;
 int    PLOTSTUFF;
@@ -104 +106 @@
 int    CHIPORDER;
 
+int UserCatalog;
+double UserCatalogRA, UserCatalogDEC;
+
 char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
 int           NphotcodesKeep,      NphotcodesSkip;
@@ -140 +145 @@
 char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
 void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
-int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
+off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
 int           args                PROTO((int argc, char **argv));
 int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
@@ -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, int Nregion, int *Nimage, int **LineNum));
+Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
 void set_db (FITS_DB *in);
-int Shutdown (char *format, ...);
+int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2);
 void TrapSignal (int sig);
 void SetProtect (int mode);
 int SetSignals (void);
 
-GSCRegion    *find_regions        PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion));
+GSCRegion    *find_regions        PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion));
 void          freeGridBins        PROTO((int Ncatalog));
 void          freeImageBins       PROTO((int Ncatalog));
@@ -166 +171 @@
 void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
 int           gcatalog            PROTO((Catalog *catalog, int FINAL));
-Coords       *getCoords           PROTO((int meas, int cat));
-float         getMcal             PROTO((int meas, int cat));
-float         getMgrid            PROTO((int meas, int cat));
-float         getMmos             PROTO((int meas, int cat));
-float         getMrel             PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions));
-void          getfullregion       PROTO((Image *image, int Nimage, GSCRegion *fullregion));
-Image        *getimage            PROTO((int N));
-Image        *getimages           PROTO((int *N));
+// Coords       *getCoords           PROTO((off_t meas, int cat));
+float         getMcal             PROTO((off_t meas, int cat));
+float         getMgrid            PROTO((off_t meas, int cat));
+float         getMmos             PROTO((off_t meas, int cat));
+float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
+GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions));
+void          getfullregion       PROTO((Image *image, off_t Nimage, GSCRegion *fullregion));
+Image        *getimage            PROTO((off_t N));
+Image        *getimages           PROTO((off_t *N));
 void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
 void          initGrid            PROTO((int dX, int dY));
 void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
-void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
-void          initImages          PROTO((Image *input, int N));
+void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int FULLINIT));
+void          initImages          PROTO((Image *input, off_t N));
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
-void          initMosaicGrid      PROTO((Image *image, int Nimage));
-void          initMosaics         PROTO((Image *image, int Nimage));
+void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
+void          initMosaics         PROTO((Image *image, off_t Nimage));
 void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
 void          initialize          PROTO((int argc, char **argv));
@@ -189 +194 @@
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
 SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
-Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
 
 void check_permissions (char *basefile);
@@ -198 +203 @@
 
 int           main                PROTO((int argc, char **argv));
-void          mark_images         PROTO((Image *image, int Nimage, Image *timage, int Ntimage));
-void          matchImage          PROTO((Catalog *catalog, int meas, int cat));
-void          matchMosaics        PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *name_region         PROTO((char *name, int *Nregions));
+void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
+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));
 double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
 void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
@@ -224 +229 @@
 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
 void          set_ZP              PROTO((double ZERO));
-int           setrefcode          PROTO((Image *image, int Nimage));
+int           setrefcode          PROTO((Image *image, off_t Nimage));
 void          skip_measurements   PROTO((Catalog *catalog, int pass));
 void          sortA               PROTO((double *X, int N));
@@ -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);
@@ -268 +273 @@
 int UpdateMosaic (Catalog *catalog, int Ncatalog);
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
-void fixImageRaw (Catalog *catalog, int Ncatalog, int im);
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
 void FlagOutliers(Catalog *catalog);
 int MeasFilterTest(Measure *measure);
 
 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);
 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
 
-Mosaic *getMosaicForImage (int N);
-
-StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
-StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode);
-
-Mosaic *getmosaics (int *N);
-void initMosaics (Image *image, int Nimage);
-StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
-StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars);
-Mosaic *getMosaicForImage (int im);
+Mosaic *getMosaicForImage (off_t N);
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+
+Mosaic *getmosaics (off_t *N);
+void initMosaics (Image *image, off_t Nimage);
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+Mosaic *getMosaicForImage (off_t im);
 
 double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
@@ -297 +302 @@
 int relastro_objects (void);
 int UpdateObjectOffsets (SkyList *skylist);
+
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+int relastroVisualPlotScatter(double values[], double thresh, int npts);
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
+                               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);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/ConfigInit.c

@@ -60 +60 @@
   /* update master photcode table if not defined */
   sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR);
-  if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) {
+  if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile, TRUE)) {
     fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile);
     exit (1);

branches/simtest_nebulous_branches/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;
 }
 

branches/simtest_nebulous_branches/Ohana/src/relastro/src/FitSimple.c

@@ -17 +17 @@
 
     // XXX the value needs to be set in a more intelligent way
-    if (dR > 0.15) continue;
+    if (dR > 0.50) continue;
     
     fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);

branches/simtest_nebulous_branches/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:

branches/simtest_nebulous_branches/Ohana/src/relastro/src/ImageOps.c

@@ -1 +1 @@
 # include "relastro.h"
-#define TESTING
-
+
+# define USE_IMAGE_ID 1
+
+static off_t       **bin;     // link from catalog,measure to image
+static int         **clist;   // catalog which supplied measurement on image
+static off_t       **mlist;   // measure reference for measurement on image
+static off_t        *Nlist;   // number of measurements on image
+static off_t        *NLIST;   // allocated number of measurements on image   
+
+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
+// index
+
+# if USE_IMAGE_ID
+static off_t        *imageIDs; // list of all image IDs
+static off_t        *imageIdx; // list of index for image IDs 
+# else
 static unsigned int *start;
 static unsigned int *stop;
-static int         **bin;
-
-static int         **clist;
-static int         **mlist;
-static int          *Nlist;
-static int          *NLIST;
-
-static Image        *image;
-static int          Nimage;
-
-Image *getimages (int *N) {
+# endif
+
+Image *getimages (off_t *N) {
   *N = Nimage;
   return (image);
 }
 
-Image *getimage (int N) {
+Image *getimage (off_t N) {
   return (&image[N]);
 }
 
-void initImages (Image *input, int N) {
-
-  int i;
+int *getCatlist (int *N, off_t im) {
+
+  *N = Ncatlist[im];
+  return (catlist[im]);
+}
+
+void initImages (Image *input, off_t N) {
+
+  off_t i;
 
   image = input;
   Nimage = N;
 
-  ALLOCATE (start,   unsigned, Nimage);
-  ALLOCATE (stop,    unsigned, Nimage);
+# if USE_IMAGE_ID
+  ALLOCATE (imageIDs, off_t, Nimage);
+  ALLOCATE (imageIdx, off_t, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    imageIdx[i] = i;
+    imageIDs[i] = image[i].imageID;
+  }
+  llsortpair (imageIDs, imageIdx, Nimage);
+# else
+  ALLOCATE (start, unsigned, Nimage);
+  ALLOCATE (stop, unsigned, Nimage);
 
   for (i = 0; i < Nimage; i++) {
@@ -37 +66 @@
     stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
   }
-}
-
-void initImageBins (Catalog *catalog, int Ncatalog) {
-
-  int i, j;
-
-  ALLOCATE (bin, int *, Ncatalog);
+# endif
+}
+
+off_t getImageByID (off_t ID) {
+
+  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
+  // use bisection to find the specified image ID
+
+# if USE_IMAGE_ID
+  off_t Nlo, Nhi, N;
+
+  Nlo = 0; Nhi = Nimage;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (imageIDs[N] < ID) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nimage);
+    }
+  }
+
+  for (N = Nlo; N < Nhi; N++) {
+    if (imageIDs[N] == ID)
+      return (imageIdx[N]);
+  }
+# endif
+
+  return (-1);
+}
+
+// these are really image & catalog indexes
+void initImageBins (Catalog *catalog, int Ncatalog, int FULLINIT) {
+
+  off_t i, j;
+
+  ALLOCATE (bin, off_t *, Ncatalog);
   for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1));
+    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
     for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
   }
 
-  ALLOCATE (Nlist, int, Nimage);
-  ALLOCATE (NLIST, int, Nimage);
-  ALLOCATE (clist, int *, Nimage);
-  ALLOCATE (mlist, int *, Nimage);
+  ALLOCATE (Nlist, off_t,   Nimage);
+  ALLOCATE (NLIST, off_t,   Nimage);
+  ALLOCATE (clist, int *,   Nimage);
+  ALLOCATE (mlist, off_t *, Nimage);
 
   for (i = 0; i < Nimage; i++) {
-    Nlist[i] = 0;
+    Nlist[i] =   0;
     NLIST[i] = 100;
     ALLOCATE (clist[i], int, NLIST[i]);
-    ALLOCATE (mlist[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]);
+    }
   }
 }
@@ -64 +134 @@
 void freeImageBins (int Ncatalog) {
 
-  int i;
+  off_t i;
 
   for (i = 0; i < Ncatalog; i++) {
@@ -79 +149 @@
 
 /* match measurements to images */
-void findImages (Catalog *catalog, int Ncatalog) {
-
-  int i, j;
+void findImages (Catalog *catalog, int Ncatalog, int MATCHCAT) {
+
+  off_t i, j;
   char *name;
 
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
-      matchImage (catalog, j, i);
-    }
-  }
-
-  for (i = 0; VERBOSE && (i < Nimage); i++) {
+      // XXX : update to optionally restrict by photcode equiv a la relphot
+      // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
+      // if (photcode[0].code != ecode) continue;
+      matchImage (catalog, j, i, MATCHCAT);
+    }
+  }
+
+  for (i = 0; VERBOSE2 && (i < Nimage); i++) {
     name = GetPhotcodeNamebyCode (image[i].photcode);
-    fprintf (stderr, "image %d has %d measures (%s, %s)\n", i, Nlist[i],
+    fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
              ohana_sec_to_date(image[i].tzero), name);
   }
 }
 
-/* modify this function to use the measure->imageID field */
-void matchImage (Catalog *catalog, int meas, int cat) {
-
-  int i;
+# if USE_IMAGE_ID
+// this is the imageID-based match
+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];
+
+  ID = measure[0].imageID;
+  idx = getImageByID (ID);
+  if (idx == -1) {
+    if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
+    return;
+  }
+
+  // index for (catalog, measure) -> image
+  bin[cat][meas] = idx;
+
+  // index for image, Nentry -> catalog
+  clist[idx][Nlist[idx]] = cat;
+
+  // index for image, Nentry -> measure
+  mlist[idx][Nlist[idx]] = meas;
+  Nlist[idx] ++;
+
+  if (Nlist[idx] == NLIST[idx]) {
+    NLIST[idx] += 100;
+    REALLOCATE (clist[idx], int,   NLIST[idx]);
+    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;
+}
+
+# else
+// this is the time-based match
+void matchImage (Catalog *catalog, off_t meas, int cat, int MATCHCAT) {
+
+  off_t i;
   Measure *measure;
 
@@ -107 +233 @@
   /* find the image that supplied this measurement */
   for (i = 0; i < Nimage; i++) {
+    // let's try the very slow method first before adding a bisection search
+    // if (image[i].imageID != measure[0].imageID) continue;
     if (image[0].photcode == -1) continue;
     if (measure[0].photcode != image[i].photcode) continue;
@@ -124 +252 @@
     if (Nlist[i] == NLIST[i]) {
       NLIST[i] += 100;
-      REALLOCATE (clist[i], int, NLIST[i]);
-      REALLOCATE (mlist[i], int, NLIST[i]);
+      REALLOCATE (clist[i], int,   NLIST[i]);
+      REALLOCATE (mlist[i], off_t, NLIST[i]);
     }
     return;
   }
-}
-
-Coords *getCoords (int meas, int cat) {
-
-  int i;
+  if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n");
+  return;
+}
+# endif
+
+/*
+Coords *getCoords (off_t meas, int cat) {
+
+  off_t i;
 
   i = bin[cat][meas];
@@ -139 +271 @@
   return (&image[i].coords);
 }
+*/
 
 void plot_images () {
 
-  int i, bin;
+  off_t i, bin;
   double *xlist, *Mlist, *dlist;
   Graphdata graphdata;
@@ -188 +321 @@
 // 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;
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
+
+  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) {
-      moscoords = &mosaic[0].coords;
+  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++) {
@@ -226 +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;
@@ -251 +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;
 }
@@ -257 +488 @@
 // chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
 // which the detections have already been filtered when they were loaded (bcatalog)
-StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
-
-  int i, m, c, n;
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
 
   Mosaic *mosaic;
@@ -265 +496 @@
   StarData *raw;
 
-  ALLOCATE (raw, StarData, Nlist[im]);
-
   mosaic = NULL;
   moscoords = NULL;
   if (mode == MODE_MOSAIC) {
-      mosaic = getMosaicForImage (im);
-      if (mosaic == NULL) {
-        fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
-        exit (1);
-      }
-      moscoords = &mosaic[0].coords;
-  }
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) {
+      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
+      return NULL;
+    }
+    moscoords = &mosaic[0].coords;
+  }
+
+  ALLOCATE (raw, StarData, Nlist[im]);
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -319 +550 @@
         LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
         break;
-    default:
-      fprintf (stderr, "error: invalid mode in getImageRaw");
-      abort ();
+      default:
+        fprintf (stderr, "error: invalid mode in getImageRaw");
+        abort ();
     }
   }
@@ -332 +563 @@
 // the sky positions: R,D -> P,Q -> L,M -> X,Y
 
-StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, CoordMode mode) {
-
-  int i, m, c, n;
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
 
   Mosaic *mosaic;
   Coords *moscoords;
   StarData *ref;
-
-  ALLOCATE (ref, StarData, Nlist[im]);
 
   mosaic = NULL;
@@ -348 +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++) {
@@ -371 +602 @@
     switch (mode) {
       case MODE_SIMPLE:
-      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
-      ref[i].L = ref[i].P;
-      ref[i].M = ref[i].Q;
-      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
-      break;
+        RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
+        ref[i].L = ref[i].P;
+        ref[i].M = ref[i].Q;
+        LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+        break;
       case MODE_MOSAIC:
-      RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, 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;
+        RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, 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;
       default:
         fprintf (stderr, "invalid case");
@@ -399 +630 @@
   return;
 
-  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
   double Ns, theta, x, y;
   double *R, *D, *dR, *dD;
@@ -490 +722 @@
 void FlagOutliers2D (Catalog *catalog) {
 
-  int i, j, k, m, N, Ndel, Nave, Nmax, TOOFEW, Nsecfilt;
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
   double *index;
   double Ns, theta, x, y;
@@ -569 +802 @@
     // recalculate image center, sigma based on closest 50% of points
     for(k = 0;  k < N; k++) {
-      int ind = (int) index[k];
+      off_t ind = (off_t) index[k];
       R[k] = catalog[0].measure[ind].dR;
       D[k] = catalog[0].measure[ind].dD;

branches/simtest_nebulous_branches/Ohana/src/relastro/src/MosaicOps.c

@@ -2 +2 @@
 
 // array of mosaic definition structures
-static int    Nmosaic;
+static off_t   Nmosaic;
 static Mosaic *mosaic;
 
 // list of all images associated with a mosaic
-static int   *Nmosaic_own_images; // number of images for this mosaic
-static int   *Amosaic_own_images; // size of allocated array
-static int   **mosaic_own_images; // array of arrays: mosaic -> images
+static off_t   *Nmosaic_own_images; // number of images for this mosaic
+static off_t   *Amosaic_own_images; // size of allocated array
+static off_t   **mosaic_own_images; // array of arrays: mosaic -> images
 
 // list of mosaic associated with each image  
-static int    Nmosaic_for_images; // number of images (for internal checks)
-static int    *mosaic_for_images; // array of: image -> mosaic
-
-Mosaic *getmosaics (int *N) {
+static off_t    Nmosaic_for_images; // number of images (for off_ternal checks)
+static off_t    *mosaic_for_images; // array of: image -> mosaic
+
+Mosaic *getmosaics (off_t *N) {
   *N = Nmosaic;
   return (mosaic);
 }
 
+off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nmosaic;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (startMos[N] < start) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nmosaic);
+    }
+  }
+
+  // check for the matched mosaic starting from Nlo 
+  // we may have to go much beyond Nlo since stop is not sorted
+  // can we use a sorted version of stop to check when we are beyond the valid range??
+  for (N = Nlo; N < Nmosaic; N++) { 
+    if (stop  < stopMos[N]) continue;
+    if (start > startMos[N])  continue;
+    if (stop  < startMos[N]) return (-1);
+    return (indexMos[N]);
+  }
+
+  return (-1);
+}
+
+// sort two times vectors and an index by first time vector
+void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \
+  tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \
+  tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \
+  tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \
+}
+# define COMPARE(A,B)(S[A] < S[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+// first, let's continue to use the time to make the match, but use bracketing to make it faster:
+
 // find mosaic frames (unique time periods & photcode name matches mosaic) 
-void initMosaics (Image *image, int Nimage) {
-
-  int i, j, found, NMOSAIC;
+void initMosaics (Image *image, off_t Nimage) {
+
+  off_t i, Nmos, NMOSAIC;
   unsigned int start, stop;
+
+  unsigned int *startMos, *stopMos;
+  off_t *indexMos;
 
   Nmosaic = 0;
@@ -29 +81 @@
   ALLOCATE (mosaic, Mosaic, NMOSAIC);
 
-  ALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-  ALLOCATE (Amosaic_own_images, int, NMOSAIC);
-  ALLOCATE (mosaic_own_images, int *, NMOSAIC);
+  ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+  ALLOCATE (startMos, unsigned int, NMOSAIC);
+  ALLOCATE (stopMos, unsigned int, NMOSAIC);
+  ALLOCATE (indexMos, off_t, NMOSAIC);
 
   /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
@@ -54 +109 @@
     Nmosaic_own_images[Nmosaic] = 0;
     Amosaic_own_images[Nmosaic] = 10;
-    ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]);
+    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
+
+    startMos[Nmosaic] = start;
+    stopMos[Nmosaic] = stop;
+    indexMos[Nmosaic] = Nmosaic;
 
     Nmosaic ++;
@@ -60 +119 @@
       NMOSAIC += 10;
       REALLOCATE (mosaic, Mosaic, NMOSAIC);
-      REALLOCATE (mosaic_own_images, int *, NMOSAIC);
-      REALLOCATE (Nmosaic_own_images, int, NMOSAIC);
-      REALLOCATE (Amosaic_own_images, int, NMOSAIC);
-    }
-  }
-
+      REALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+      REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (startMos, unsigned int, NMOSAIC);
+      REALLOCATE (stopMos, unsigned int, NMOSAIC);
+      REALLOCATE (indexMos, off_t, NMOSAIC);
+    }
+  }
+
+  // sort the index, start, and stop by the start times:
+  sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic);
+  
   // array to store image->mosaic index
   Nmosaic_for_images = Nimage;
-  ALLOCATE (mosaic_for_images, int, Nmosaic_for_images);
+  ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images);
 
   /* now assign the WRP images to these mosaics */
@@ -80 +145 @@
     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
 
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-      break;
-    }
-    /* if no matching mosaic exists, skip this image */
-    if (!found) continue;
+    Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
+    if (Nmos == -1) {
+      fprintf (stderr, "cannot match mosaic for %s\n", image[i].name);
+      continue;
+    }
 
     // mosaic corresponding to this image
-    mosaic_for_images[i] = j;
+    mosaic_for_images[i] = Nmos;
 
     // add image to mosaic_own_image list 
-    mosaic_own_images[j][Nmosaic_own_images[j]] = i;
-    Nmosaic_own_images[j] ++;
-    if (Nmosaic_own_images[j] == Amosaic_own_images[j]) {
-      Amosaic_own_images[j] += 10;
-      REALLOCATE (mosaic_own_images[j], int, Amosaic_own_images[j]);
-    }
-    assert (Nmosaic_own_images[j] < Amosaic_own_images[j]);
-  }
-
+    mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i;
+    Nmosaic_own_images[Nmos] ++;
+    if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) {
+      Amosaic_own_images[Nmos] += 10;
+      REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]);
+    }
+    assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]);
+  }
+
+  free (startMos);
+  free (stopMos);
+  free (indexMos);
   return;
 }
@@ -109 +172 @@
 // return StarData values for detections in the specified image, converting coordinates from the
 // chip positions: X,Y -> L,M -> P,Q -> R,D
-StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nraw, Nnew;
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nraw, Nnew;
   StarData *raw, *new;
 
@@ -125 +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
@@ -142 +209 @@
 // return StarData values for averages positions in the specified image, converting coordinates from
 // the sky positions: R,D -> P,Q -> L,M -> X,Y
-StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {
-
-  int i, j, im, Nref, Nnew;
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nref, Nnew;
   StarData *ref, *new;
   
@@ -157 +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
@@ -172 +243 @@
 }
 
-Mosaic *getMosaicForImage (int im) {
-
-  int mos;
+Mosaic *getMosaicForImage (off_t im) {
+
+  off_t mos;
 
   if (im < 0) abort();

branches/simtest_nebulous_branches/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);
 

branches/simtest_nebulous_branches/Ohana/src/relastro/src/UpdateChips.c

@@ -4 +4 @@
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  int i, Nimage, Nraw, Nref;
+  off_t i, Nimage, Nraw, Nref;
   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 %d : Nstars: %d\n", i, 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);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/UpdateMeasures.c

@@ -3 +3 @@
 int UpdateMeasures (Catalog *catalog, int Ncatalog) {
 
-  int i, Nimage;
+  off_t i, Nimage;
   Image *image;
 

branches/simtest_nebulous_branches/Ohana/src/relastro/src/UpdateMosaic.c

@@ -4 +4 @@
 
   /* we can measure new image parameters for each mosaic independently */
-  int i, Nmosaic, Nstars;
+  off_t i, Nmosaic, Nstars;
   Mosaic *mosaic;
   StarData *raw, *ref;

branches/simtest_nebulous_branches/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

branches/simtest_nebulous_branches/Ohana/src/relastro/src/UpdateObjects.c

@@ -1 +1 @@
 # include "relastro.h"
 
-static int Nmax;
+static off_t   Nmax;
 static double *X, *dX;
 static double *Y, *dY;
@@ -13 +13 @@
 void initObjectData (Catalog *catalog, int Ncatalog) {
 
-  int i, j;
+  off_t i, j;
   
   Nmax = 0;
@@ -40 +40 @@
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
-  int i, j, k, m, N, Nsecfilt, found, kp;
+  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;
-  float mag;
-  int mask;
-  PhotCode *code;
-
+  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;
+
+  memset (&fit,    0, sizeof(fit));
+  memset (&fitAve, 0, sizeof(fitAve));
+  memset (&fitPM,  0, sizeof(fitPM));
+  memset (&fitPAR, 0, sizeof(fitPAR));
+  
   initObjectData (catalog, Ncatalog);
 
@@ -64 +68 @@
   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?
@@ -72 +77 @@
   assert (catalog[0].Nsecfilt == Nsecfilt);
 
+  NaveSum = NparSum = NpmSum = NskipSum = 0;
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure);
-
-    Nskip = 0;
+    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure);
+
+    Nave = Npar = Npm = Nskip = 0;
     for (j = 0; j < catalog[i].Naverage; j++) {
       /* calculate the average value of R,D for a single star */
@@ -82 +88 @@
       // 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++) {
 
@@ -111 +119 @@
         // 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;
@@ -122 +130 @@
         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++;
@@ -141 +158 @@
 
       // 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;
 
@@ -147 +165 @@
       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;
       }
@@ -154 +176 @@
       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]);
@@ -270 +301 @@
       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);
 }
@@ -282 +325 @@
 /* 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?
 
 */

branches/simtest_nebulous_branches/Ohana/src/relastro/src/UpdateSimple.c

@@ -4 +4 @@
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  int i, Nimage, Nstars;
+  off_t i, Nimage, Nstars;
   Image *image;
   StarData *raw, *ref;
@@ -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);

branches/simtest_nebulous_branches/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");

branches/simtest_nebulous_branches/Ohana/src/relastro/src/bcatalog.c

@@ -3 +3 @@
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
 
-  int i, j, k, offset, found;
-  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
-  float mag;
-  int mask;
-  PhotCode *code;
+  off_t i, j, offset;
+  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
+  int Nsecfilt;
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -106 +104 @@
 
   if (VERBOSE) {
-    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i,
-             subcatalog[0].Naverage, subcatalog[0].Nmeasure);
+    fprintf (stderr, "%lld: using %lld stars (%lld measures) for catalog\n", (long long) i, (long long) subcatalog[0].Naverage, (long long) subcatalog[0].Nmeasure);
    }
   return (TRUE);

branches/simtest_nebulous_branches/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);
+}
+

branches/simtest_nebulous_branches/Ohana/src/relastro/src/initialize.c

@@ -5 +5 @@
   int NPHOTCODES;
   char *codename, *ptr, *list;
+
+  ptr = NULL;
 
   ConfigInit (&argc, argv);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/load_images.c

@@ -1 +1 @@
 # include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 SkyList *load_images (FITS_DB *db, SkyRegion *region) {
 
   Image     *image, *subset;
-  int        Nimage, Nsubset;
-  int       *LineNumber;
+  off_t      Nimage, Nsubset;
+  off_t     *LineNumber;
+  struct timeval start, stop;
 
   SkyTable *sky = NULL;
   SkyList *skylist = NULL;
 
+  gettimeofday (&start, (void *) NULL);
+
   // load the current sky table (layout of all SkyRegions) 
-  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE);
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
   SkyTableSetFilenames (sky, CATDIR, "cpt");
   
   // 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);
 
   // convert database table to internal structure
   image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+  MARKTIME("  convert image table: %f sec\n", dtime);
 
   // select the images which overlap the selected sky regions
   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+  MARKTIME("  select images: %f sec\n", dtime);
 
   gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
 
   initImages (subset, Nsubset);
+  MARKTIME("  init images: %f sec\n", dtime);
+
   initMosaics (subset, Nsubset);
+  MARKTIME("  init mosaics: %f sec\n", dtime);
   
   /* unlock, if we can (else, unlocked below) */

branches/simtest_nebulous_branches/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);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/plotstuff.c

@@ -125 +125 @@
   if (PLOTDELAY > 500000) {
     fprintf (stdout, "press return\n"); 
-    fscanf (stdin, "%*c");
+    if (fscanf (stdin, "%*c") != 0) fprintf (stderr, "\n");
   } else {
     usleep (PLOTDELAY);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/relastro.c

@@ -1 +1 @@
 # include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
 
 int main (int argc, char **argv) {
 
-  int status, Ncatalog;
+  int i, status, Ncatalog;
   Catalog *catalog;
   FITS_DB db;
+  struct timeval start, stop;
+  SkyList *skylist = NULL;
 
-  SkyList *skylist = NULL;
+  gettimeofday (&start, (void *) NULL);
 
   /* get configuration info, args */
@@ -26 +34 @@
   if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
   if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+  MARKTIME("load image data: %f sec\n", dtime);
 
   /* load regions and images based on specified sky patch (default depth) */
   skylist = load_images (&db, &UserPatch);
+  MARKTIME("load images: %f sec\n", dtime);
+
+  initCoords();
 
   /* load catalog data from region files : subselect high-quality measurements */
   catalog = load_catalogs (skylist, &Ncatalog, TRUE);
+  MARKTIME("load catalog data: %f sec\n", dtime);
 
   /* match measurements with images */
-  initImageBins (catalog, Ncatalog);
-  findImages (catalog, Ncatalog);
+  initImageBins (catalog, Ncatalog, TRUE);
+  MARKTIME("make image bins: %f sec\n", dtime);
+
+  findImages (catalog, Ncatalog, TRUE);
+  MARKTIME("set up image indexes: %f sec\n", dtime);
 
   if (PLOTSTUFF) {
@@ -45 +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;
 
@@ -69 +95 @@
   UpdateObjectOffsets (skylist);
 
+  // iterate over catalogs to make detection coordinates consistant
+  FixProblemImages (skylist);
+
   // save the updated image parameters
   dvo_image_update (&db, VERBOSE);

branches/simtest_nebulous_branches/Ohana/src/relastro/src/relastro_objects.c

@@ -11 +11 @@
 
   // load the current sky table (layout of all SkyRegions) 
-  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE);
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
   SkyTableSetFilenames (sky, CATDIR, "cpt");
   
   // 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

branches/simtest_nebulous_branches/Ohana/src/relastro/src/select_images.c

@@ -13 +13 @@
 } SkyRegionCoords;
 
-Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) {
+void dsortindex (double *X, off_t *Y, int N);
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
+Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
   
   Image *image;
-  int i, j, k, m, found, nimage, NIMAGE;
-  int InRange;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t *line_number;
+  int InRange, found;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
-  int *line_number;
   Coords tcoords;
   SkyRegionCoords *skycoords;
   
+  double *RmaxSky;
+  off_t *index;
+
   if (skylist[0].Nregions < 1) {
     *Nimage = 0;
@@ -38 +44 @@
 
   ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
+
+  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
+  ALLOCATE (index, off_t, skylist[0].Nregions);
 
   /* compare with each region file */
@@ -56 +65 @@
     skycoords[i].Yc[4] = skycoords[i].Yc[0];    
 
+    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
+    index[i] = i;
+
     dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
     dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
@@ -65 +77 @@
   }
 
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
+
   if (VERBOSE) fprintf (stderr, "finding images\n");
   BuildChipMatch (timage, Ntimage);
@@ -71 +85 @@
   NIMAGE = 100;
   ALLOCATE (image, Image, NIMAGE);
-  ALLOCATE (line_number, int, NIMAGE);
+  ALLOCATE (line_number, off_t, NIMAGE);
   
   // go through the complete list of images, selecting ones which overlap any region
@@ -100 +114 @@
     }
     
-    if (!FindMosaicForImage (timage, Ntimage, i)) continue;
+    if (!FindMosaicForImage (timage, Ntimage, i)) {
+      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
+      continue;
+    }
 
     /* define image corners */
@@ -110 +127 @@
     found = FALSE;
 
-    /* transform to ra,dec */
+    /* transform corners to ra,dec */
+    double RminImage = 360.0;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
-    }
+      RminImage = MIN(RminImage, Ri[j]);
+    }
+
+    // RA(nStart) is guaranteed to be < RminImage:
+    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
 
     /* compare with each region file */
-    for (m = 0; (m < skylist[0].Nregions) && !found; m++) { 
+    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+
+      m = index[iSky];
 
       /* we make positional comparisons in the projection of catalog */
@@ -122 +146 @@
       tcoords.crval2 = skycoords[m].Dc;
 
-      /* transform to ra,dec */
+      /* transform corner coords to X,Y in this catalog system */
       InRange = TRUE;
       for (j = 0; (j < 5) && InRange; j++) {
@@ -151 +175 @@
       }
       if (RESET) {
-        // XXX do we need this in relastro?
+        // XXX do we need / want to do this in relastro?
         assignMcal (&image[nimage], (double *) NULL, -1);
         image[nimage].dMcal = NAN;
@@ -161 +185 @@
         NIMAGE += 100;
         REALLOCATE (image, Image, NIMAGE);
-        REALLOCATE (line_number, int, NIMAGE);
+        REALLOCATE (line_number, off_t, NIMAGE);
       }
     }
   }
       
-  if (VERBOSE) fprintf (stderr, "found %d images\n", nimage);
+  if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage);
 
   REALLOCATE (image, Image, MAX (nimage, 1));
-  REALLOCATE (line_number, int, MAX (nimage, 1));
+  REALLOCATE (line_number, off_t, MAX (nimage, 1));
   free (skycoords);
 
@@ -253 +277 @@
 }
 
+void dsortindex (double *X, off_t *Y, int N) {
+
+# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
+  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
+  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nregions;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (Rref[N] < R) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nregions);
+    }
+  }
+  return (Nlo);
+}