Changeset 33652

Timestamp:

Apr 1, 2012, 3:00:19 PM (14 years ago)

Author:

eugene

Message:

updates to relastro to support distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relastro does not reset photom bits); skip objects with NAN error bars (can be used to exclude some photcodes from astrometry); enum for mask bits in FitChip and related; choose max order using same rules as psastro; plug leaks; reduce verbosity; big re-org of main (for parallel code); only use accepted measurements for Tmin,Tmax,Trange,Tmean; require a min parallax factor to fit parallaxes; use MeasureTiny where possible; use new ParsePhotcodeList function; testparallax function; threaded version of UpdateChips; bcatalog limit density selecting objects by number of measurements; high-speed now copies support tables to output database; high-speed may write to a parallel database; test for unsorted database (& fail if true)

Location:

trunk/Ohana/src/relastro

Files:

• Makefile (modified) (4 diffs)
• doc/parallel.txt (copied) (copied from branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/parallel.txt )
• doc/plx_factor.pro (copied) (copied from branches/eam_branches/ipp-20111122/Ohana/src/relastro/doc/plx_factor.pro )
• include/relastro.h (modified) (15 diffs)
• src/BrightCatalog.c (copied) (copied from branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/BrightCatalog.c )
• src/ConfigInit.c (modified) (2 diffs)
• src/FitChip.c (modified) (10 diffs)
• src/FixProblemImages.c (modified) (2 diffs)
• src/GetAstromError.c (modified) (2 diffs)
• src/ImageOps.c (modified) (31 diffs)
• src/ParFactor.c (modified) (5 diffs)
• src/StarMaps.c (modified) (2 diffs)
• src/UpdateChips.c (modified) (7 diffs)
• src/UpdateMeasures.c (modified) (1 diff)
• src/UpdateObjectOffsets.c (modified) (4 diffs)
• src/UpdateObjects.c (modified) (9 diffs)
• src/args.c (modified) (11 diffs)
• src/bcatalog.c (modified) (14 diffs)
• src/dvo_astrom_ops.c (modified) (5 diffs)
• src/fitpoly.c (modified) (1 diff)
• src/high_speed_catalogs.c (modified) (4 diffs)
• src/high_speed_objects.c (modified) (11 diffs)
• src/high_speed_utils.c (modified) (8 diffs)
• src/initialize.c (modified) (3 diffs)
• src/load_catalogs.c (modified) (4 diffs)
• src/load_images.c (modified) (5 diffs)
• src/relastro.c (modified) (1 diff)
• src/relastro_client.c (copied) (copied from branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_client.c )
• src/relastro_images.c (copied) (copied from branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_images.c )
• src/relastro_merge_source.c (modified) (4 diffs)
• src/relastro_objects.c (modified) (5 diffs)
• src/resort_catalog.c (modified) (1 diff)
• src/save_catalogs.c (modified) (1 diff)
• src/select_images.c (modified) (8 diffs)
• src/testparallax.c (modified) (1 diff)

trunk/Ohana/src/relastro/Makefile

@@ -1 @@
-default: relastro
+default: relastro relastro_client
 help:
         @echo "make options: relastro (default)"
@@ -18 +18 @@
 
 relastro: $(BIN)/relastro.$(ARCH)
-install: $(DESTBIN)/relastro
+relastro_client: $(BIN)/relastro_client.$(ARCH)
+
+install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client $(DESTBIN)/testparallax
 
 RELASTRO = \
@@ -53 +55 @@
 $(SRC)/select_images.$(ARCH).o       \
 $(SRC)/relastro.$(ARCH).o            \
+$(SRC)/relastro_images.$(ARCH).o    \
 $(SRC)/relastro_objects.$(ARCH).o    \
 $(SRC)/save_catalogs.$(ARCH).o       \
-$(SRC)/write_coords.$(ARCH).o        \
 $(SRC)/CoordOps.$(ARCH).o            \
 $(SRC)/FixProblemImages.$(ARCH).o    \
@@ -64 +66 @@
 $(SRC)/relastro_merge_source.$(ARCH).o  \
 $(SRC)/resort_catalog.$(ARCH).o  \
+$(SRC)/BrightCatalog.$(ARCH).o  \
 $(SRC)/relastroVisual.$(ARCH).o
 
-$(RELASTRO): $(INC)/relastro.h $(KAPA_INCS)
-$(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS)
+#$(SRC)/write_coords.$(ARCH).o        \
+
+$(RELASTRO): $(INC)/relastro.h
+$(BIN)/relastro.$(ARCH): $(RELASTRO)
+
+RELASTRO_CLIENT = \
+$(SRC)/ConfigInit.$(ARCH).o          \
+$(SRC)/FitPM.$(ARCH).o               \
+$(SRC)/FitPar.$(ARCH).o              \
+$(SRC)/FitPMandPar.$(ARCH).o         \
+$(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/ImageOps.$(ARCH).o            \
+$(SRC)/MosaicOps.$(ARCH).o           \
+$(SRC)/ParFactor.$(ARCH).o           \
+$(SRC)/SetSignals.$(ARCH).o          \
+$(SRC)/Shutdown.$(ARCH).o            \
+$(SRC)/UpdateObjects.$(ARCH).o       \
+$(SRC)/UpdateObjectOffsets.$(ARCH).o \
+$(SRC)/UpdateMeasures.$(ARCH).o      \
+$(SRC)/GetAstromError.$(ARCH).o      \
+$(SRC)/args.$(ARCH).o                \
+$(SRC)/bcatalog.$(ARCH).o            \
+$(SRC)/dvo_astrom_ops.$(ARCH).o      \
+$(SRC)/fitpoly.$(ARCH).o             \
+$(SRC)/initialize.$(ARCH).o          \
+$(SRC)/liststats.$(ARCH).o           \
+$(SRC)/load_catalogs.$(ARCH).o       \
+$(SRC)/load_images.$(ARCH).o         \
+$(SRC)/mkpolyterm.$(ARCH).o          \
+$(SRC)/select_images.$(ARCH).o       \
+$(SRC)/relastro_client.$(ARCH).o        \
+$(SRC)/relastro_objects.$(ARCH).o    \
+$(SRC)/save_catalogs.$(ARCH).o       \
+$(SRC)/CoordOps.$(ARCH).o            \
+$(SRC)/high_speed_catalogs.$(ARCH).o  \
+$(SRC)/high_speed_objects.$(ARCH).o  \
+$(SRC)/high_speed_utils.$(ARCH).o  \
+$(SRC)/relastro_merge_source.$(ARCH).o  \
+$(SRC)/resort_catalog.$(ARCH).o  \
+$(SRC)/plot_scatter.$(ARCH).o        \
+$(SRC)/plotstuff.$(ARCH).o           \
+$(SRC)/relastroVisual.$(ARCH).o \
+$(SRC)/BrightCatalog.$(ARCH).o 
+
+#$(SRC)/FitChip.$(ARCH).o             \
+#$(SRC)/FitMosaic.$(ARCH).o           \
+#$(SRC)/UpdateChips.$(ARCH).o         \
+#$(SRC)/UpdateMosaic.$(ARCH).o        \
+#$(SRC)/UpdateSimple.$(ARCH).o        \
+#$(SRC)/write_coords.$(ARCH).o        \
+#$(SRC)/FixProblemImages.$(ARCH).o    \
+#$(SRC)/StarMaps.$(ARCH).o    \
+
+$(RELASTRO_CLIENT): $(INC)/relastro.h
+$(BIN)/relastro_client.$(ARCH): $(RELASTRO_CLIENT)
+
+testparallax: $(BIN)/testparallax.$(ARCH)
+
+TESTPAR = \
+$(SRC)/FitPMandPar.$(ARCH).o         \
+$(SRC)/ParFactor.$(ARCH).o           \
+$(SRC)/fitpoly.$(ARCH).o             \
+$(SRC)/mkpolyterm.$(ARCH).o            \
+$(SRC)/testparallax.$(ARCH).o
+
+$(TESTPAR): $(INC)/relastro.h
+$(BIN)/testparallax.$(ARCH): $(TESTPAR)

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

@@ -4 +4 @@
 # include <signal.h>
 # include <assert.h>
+# include <pthread.h>
+
+// choose off_t or int depending on full-scale relphot analysis resources
+// # define IDX_T off_t
+# define IDX_T int 
 
 typedef enum {
@@ -15 +20 @@
 typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;
 
-typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE} FitTarget;
+typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS} FitTarget;
+
+typedef enum {
+  MARK_MEAS_DEFAULT  = 0x0000,
+  MARK_TOO_FEW_MEAS  = 0x0001,
+  MARK_NAN_POS_ERROR = 0x0002,
+  MARK_NAN_MAG_ERROR = 0x0004,
+  MARK_BIG_MAG_ERROR = 0x0008,
+  MARK_BIG_OFFSET    = 0x0010,
+} MeasurementMask;
+
+typedef struct {
+  Average     *average;       // array of (minimal) average data
+  MeasureTiny *measure;       // array of (minimal) measure data 
+  SecFilt     *secfilt;       // array of secfilt data (matched to average by Nsecfilt)
+  off_t       Naverage;
+  off_t       Nmeasure;
+} BrightCatalog;
+
+typedef struct {
+  Catalog *catalog;           // array of catalogs generated
+  int NCATALOG;               // number of catalogs allocated
+  int Ncatalog;               // number of catalogs generated
+  int Nsecfilt;               // number of catalogs generated
+  off_t *NAVERAGE;            // allocated Averages per catalog
+  off_t *NMEASURE;            // allocated Measures per catalog
+  int   *index;               // lookup table catID -> catalog[i]
+  int   *catIDs;              // lookup table catID <- catalog[i]
+  int    maxID;               // max catID value to date
+} CatalogSplitter;
 
 typedef struct {
@@ -86 +120 @@
 
 /* global variables set in parameter file */
-char   ImageCat[256];
-char   GSCFILE[256];
-char   CATDIR[256];
+# define MAX_PATH_LENGTH 1024
+char   ImageCat[MAX_PATH_LENGTH];
+char   GSCFILE[MAX_PATH_LENGTH];
+char   CATDIR[MAX_PATH_LENGTH];
 char   *HIGH_SPEED_DIR;
 char   CATMODE[16];    /* raw, mef, split, mysql */
 char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
-char   SKY_TABLE[256];
+char   SKY_TABLE[MAX_PATH_LENGTH];
 int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
+
+int          HOST_ID;
+char        *HOSTDIR;
+char        *BCATALOG;
 
 unsigned int OBJ_ID_SRC;
@@ -110 +149 @@
 double RADIUS; // match radius for high-speed objects
 
+int    PARALLEL;
+int    PARALLEL_MANUAL;
+int    PARALLEL_SERIAL;
+
+int    PARALLEL_OUTPUT;
+
 int    VERBOSE;
 int    VERBOSE2;
@@ -115 +160 @@
 int    RESET;
 int    NLOOP;
+int    NTHREADS;
 int    UPDATE;
 int    PLOTSTUFF;
@@ -123 +169 @@
 int    PM_TOOFEW;
 double PM_DT_MIN;
+double PAR_FACTOR_MIN;
 int    PLOTDELAY;
 int    CHIPORDER;
@@ -128 +175 @@
 int MaxDensityUse;
 double MaxDensityValue;
-
-int UserCatalog;
-double UserCatalogRA, UserCatalogDEC;
 
 char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
@@ -142 +186 @@
 SkyRegionSelection SELECTION;
 
-int AreaSelect;
-double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
-
-int ImagSelect, ImagMin, ImagMax;
+int ImagSelect;
+double ImagMin, ImagMax;
 
 double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
 
-int DophotSelect, DophotValue;
 int PhotFlagSelect, PhotFlagPoor, PhotFlagBad;
 
@@ -181 +222 @@
 char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
 void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
-off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
+off_t        *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
 int           args                PROTO((int argc, char **argv));
+int           args_client         PROTO((int argc, char **argv));
 int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
 void          clean_images        PROTO((void));
@@ -226 +268 @@
 void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
 void          initialize          PROTO((int argc, char **argv));
+void          initialize_client   PROTO((int argc, char **argv));
 void          initstats           PROTO((char *mode));
 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
-Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
-SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath));
+int           load_images         PROTO((FITS_DB *db, SkyList *skylist));
 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
 
@@ -314 +357 @@
 void FlagOutliers(Catalog *catalog);
 int MeasFilterTest(Measure *measure, int applySigmaLim);
+int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim);
 
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius);
@@ -333 +377 @@
 Mosaic *getMosaicForImage (off_t im);
 
-double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
-double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);
-int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
-int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanR (MeasureTiny *measure, Average *average, SecFilt *secfilt);
+double getMeanD (MeasureTiny *measure, Average *average, SecFilt *secfilt);
+int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
+int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
+int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
+int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);
 
 float GetAstromError (Measure *measure, int mode);
-int relastro_objects (void);
-int UpdateObjectOffsets (SkyList *skylist);
+float GetAstromErrorTiny (MeasureTiny *measure, int mode);
+int relastro_objects (SkyList *skylist, int hostID, char *hostpath);
+int relastro_images (SkyList *skylist);
+int UpdateObjectOffsets (SkyList *skylist, int hostID, char *hostpath);
 
 int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj);
@@ -361 +409 @@
 void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);
 
-int high_speed_catalogs ();
+int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath);
 int high_speed_objects (SkyRegion *region, Catalog *catalog);
 int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset);
@@ -369 +417 @@
 int createStarMapPoints();
 int checkStarMap(int N);
+int createStarMap (Catalog *catalog, int Ncatalog);
 
 int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog);
+int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog);
 
 int initializeConstraints();
@@ -378 +428 @@
 void setupAreaSelection(SkyRegion *region);
 
-int relastro_merge_source ();
+int relastro_merge_source (SkyTable *sky);
 void resort_catalog (Catalog *catalog);
+
+BrightCatalog *BrightCatalogLoad(char *filename);
+int BrightCatalogSave(char *filename, BrightCatalog *catalog);
+BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog);
+CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt);
+int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog);
+int BrightCatalogSplitFree (CatalogSplitter *catalogs);
+
+PhotCode **ParsePhotcodeList (char *rawlist, int *nphotcodes, int needAve);

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

@@ -26 +26 @@
   if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
 
-  // XXX these are used in relphot to identify poor stars / images -- define
-  // an equivalent concept here?
-  // GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
-  // GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
-  // GetConfig (config, "IMAGE_OFFSET",           "%lf", 0, &IMAGE_OFFSET);
-  // GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
-
   GetConfig (config, "PM_DT_MIN",              "%lf", 0, &PM_DT_MIN);
   GetConfig (config, "PM_TOOFEW",              "%d",  0, &PM_TOOFEW);
   GetConfig (config, "POS_TOOFEW",             "%d",  0, &POS_TOOFEW);
+
+  GetConfig (config, "PAR_FACTOR_MIN",         "%lf", 0, &PAR_FACTOR_MIN);
 
   GetConfig (config, "RELASTRO_MAP_NX",        "%d",  0, &NX_MAP);
@@ -42 +37 @@
   GetConfig (config, "ADDSTAR_RADIUS",         "%lf", 0, &ADDSTAR_RADIUS);
 
+  // force CATDIR to be absolute (so parallel mode will work)
+  GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
+  char *tmpcatdir = abspath (CATDIR, MAX_PATH_LENGTH);
+  strcpy (CATDIR, tmpcatdir);
+  free (tmpcatdir);
+
   GetConfig (config, "GSCFILE",                "%s",  0, GSCFILE);
-  GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
   ScanConfig(config, "CATMODE",                "%s",  0, CATMODE);
   ScanConfig(config, "CATFORMAT",              "%s",  0, CATFORMAT);

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

@@ -4 +4 @@
 int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
 
-  int i, NstatFull, Nstat, Niter, skip;
+  int i, NstatFull, Nstat, Niter;
   float dLsig, dMsig, dRsig;
   float dLsigFull, dMsigFull, dRsigFull;
@@ -27 +27 @@
       if (raw[i].mask) continue;
       if (isnan(raw[i].dMag)) { 
-        raw[i].mask |= 0x0004; 
+        raw[i].mask |= MARK_NAN_MAG_ERROR; 
         continue; 
       }
       if (raw[i].dMag > SIGMA_LIM) { 
-        raw[i].mask |= 0x0008; 
+        raw[i].mask |= MARK_BIG_MAG_ERROR; 
         continue; 
       } // is this redundant with ImageOps.c:915?
@@ -39 +39 @@
       dR = hypot (dL, dM);
       if (dR > dRmax) {
-        raw[i].mask |= 0x0010;
+        raw[i].mask |= MARK_BIG_OFFSET;
         continue;
       }
@@ -52 +52 @@
     // fit the requested order polynomial
     if (CHIPORDER > 0) {
-      image[0].coords.Npolyterms = CHIPORDER;
-    }
+      int Ndof_min = 3;
+      int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3);
+      int order_use = MIN (CHIPORDER, order_max);
+      if (order_use < 1) {
+        if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch);
+        image[0].flags |= ID_IMAGE_ASTROM_FEW;
+        return FALSE;
+      }
+      image[0].coords.Npolyterms = order_use;
+    }
+
     if (fit) fit_free (fit);
     fit = fit_init (image[0].coords.Npolyterms);
@@ -63 +72 @@
     }
 
+# if (0)
     // check if the fit has enough data points for the polynomial order
     skip = FALSE;
@@ -86 +96 @@
       return FALSE;
     }
+# endif
 
     // measure the fit, update the coords & object coordinates
@@ -91 +102 @@
       fprintf (stderr, "failed to fit new model\n");
       image[0].flags |= ID_IMAGE_ASTROM_FAIL;
+      if (fit) fit_free (fit);
       return FALSE;
     }
@@ -97 +109 @@
       fprintf (stderr, "failed to fit new model\n");
       image[0].flags |= ID_IMAGE_ASTROM_FAIL;
+      if (fit) fit_free (fit);
       return FALSE;
     }
@@ -114 +127 @@
   for (i = 0; i < Nmatch; i++) {
     if (!raw[i].mask) continue;
-    if (raw[i].mask & 0x01) nMask1 ++;
-    if (raw[i].mask & 0x02) nMask2 ++;
-    if (raw[i].mask & 0x04) nMask3 ++;
-    if (raw[i].mask & 0x08) nMask4 ++;
-    if (raw[i].mask & 0x10) nMask5 ++;
+    if (raw[i].mask & MARK_TOO_FEW_MEAS ) nMask1 ++;
+    if (raw[i].mask & MARK_NAN_POS_ERROR) nMask2 ++;
+    if (raw[i].mask & MARK_NAN_MAG_ERROR) nMask3 ++;
+    if (raw[i].mask & MARK_BIG_MAG_ERROR) nMask4 ++;
+    if (raw[i].mask & MARK_BIG_OFFSET   ) nMask5 ++;
   }
 
   GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
-  GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
+  GetScatterRawRef(&dLsig,     &dMsig,     &dRsig,     &Nstat,     raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
 
   int Nm = 0;
@@ -133 +146 @@
   image[0].nLinkAstrom = (Nm / Ns);
 
-  if (VERBOSE) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);
+  if (VERBOSE2) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);
 
   image[0].dXpixSys = dLsig;

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

@@ -5 +5 @@
 
 int FixProblemImages (SkyList *skylist) {
+
+    fprintf (stderr, "this function is currently not used : check on change from measure -> measureT in load_catalogs\n");
+    abort();
 
   int Nbad;
@@ -64 +67 @@
     }
 
-    catalog = load_catalogs (&sublist, &Ncat, FALSE);
+    // XXX use a different function here
+    // catalog = load_catalogs (&sublist, &Ncat, FALSE);
     assert (Ncat == Ncatlist);
 

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

@@ -1 +1 @@
 # include "relastro.h"
 # define WEIGHTED_ERRORS 1
+
+float GetAstromErrorTiny (MeasureTiny *measure, int mode) {
+
+  PhotCode *code;
+  float dPobs, dPsys, dPtotal, dM, AS, MS, dX, dY;
+
+  if (!WEIGHTED_ERRORS) {
+    // if we don't understand the errors at all, this at least lets us get things roughly
+    // right:
+    return 0.1;
+  }
+
+  switch (mode) {
+    case ERROR_MODE_RA:
+      dPobs = FromShortPixels(measure[0].dXccd);  // dXccd is a value in pixels
+      break;
+    case ERROR_MODE_DEC:
+      dPobs = FromShortPixels(measure[0].dYccd);  // dYccd is a value in pixels
+      break;
+    case ERROR_MODE_POS:
+      dX = FromShortPixels(measure[0].dXccd);  // dXccd is a value in pixels
+      dY = FromShortPixels(measure[0].dYccd);  // dYccd is a value in pixels
+      dPobs = hypot (dX, dY);
+      break;
+    default:
+      abort();
+  }
+
+  code  = GetPhotcodebyCode (measure[0].photcode);
+
+  // do not raise an exception, just send back the result
+  if (isnan(code[0].astromErrSys)) return NAN;
+
+  AS    = code[0].astromErrScale;
+  MS    = code[0].astromErrMagScale;
+  dPsys = code[0].astromErrSys;
+  dM    = measure[0].dM;
+  dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM));
+
+  dPtotal = MAX (dPtotal, MIN_ERROR);
+  return (dPtotal);
+}
 
 float GetAstromError (Measure *measure, int mode) {
@@ -30 +72 @@
 
   code  = GetPhotcodebyCode (measure[0].photcode);
+
+  // do not raise an exception, just send back the result
+  if (isnan(code[0].astromErrSys)) return NAN;
+
   AS    = code[0].astromErrScale;
   MS    = code[0].astromErrMagScale;

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

@@ -3 +3 @@
 # 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   
-
+// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage)
 static Image        *image;   // list of available images
 static off_t        Nimage;   // number of available images
+
+// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row
+// each image belongs
 static off_t       *LineNumber; // match of subset to full image table
 
-static int         *Ncatlist;  // catalogs associated with each image
-static int         *NCATLIST;  // catalogs associated with each image
+static off_t        *N_onImage;   // number of measurements on image
+static off_t        *N_ONIMAGE;   // allocated number of measurements on image   
+
+static int          *Ncatlist;  // catalogs associated with each image
+static int          *NCATLIST;  // catalogs associated with each image
 static int         **catlist;  // catalogs associated with each image
+
+static IDX_T       **MeasureToImage;     // link from catalog,measure to image
+static IDX_T       **ImageToCatalog;   // catalog which supplied measurement on image
+static IDX_T       **ImageToMeasure;   // measure reference for measurement on image
 
 // if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
@@ -28 +33 @@
 # endif
 
+// MeasureToImage was 'bin'
+// ImageToCatalog was 'clist'
+// ImageToMeasure was 'mlist'
+
+// N_onImage was 'Nlist'
+// N_ONIMAGE was 'NLIST'
+
 Image *getimages (off_t *N, off_t **line_number) {
 
@@ -105 +117 @@
   off_t i, j;
 
-  ALLOCATE (bin, off_t *, Ncatalog);
+  ALLOCATE (MeasureToImage, IDX_T *, Ncatalog);
   for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
-    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
-  }
-
-  ALLOCATE (Nlist, off_t,   Nimage);
-  ALLOCATE (NLIST, off_t,   Nimage);
-  ALLOCATE (clist, int *,   Nimage);
-  ALLOCATE (mlist, off_t *, Nimage);
+    ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1));
+    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1;
+  }
+
+  ALLOCATE (N_onImage, off_t,   Nimage);
+  ALLOCATE (N_ONIMAGE, off_t,   Nimage);
+  ALLOCATE (ImageToCatalog, IDX_T *,   Nimage);
+  ALLOCATE (ImageToMeasure, IDX_T *, Nimage);
 
   for (i = 0; i < Nimage; i++) {
-    Nlist[i] =   0;
-    NLIST[i] = 100;
-    clist[i] = NULL;  // we allocate these iff they are needed in matchImage
-    mlist[i] = NULL;  // we allocate these iff they are needed in matchImage
+    N_onImage[i] =   0;
+    N_ONIMAGE[i] =  30;
+    ImageToCatalog[i] = NULL;  // we allocate these iff they are needed in matchImage
+    ImageToMeasure[i] = NULL;  // we allocate these iff they are needed in matchImage
   }
 
@@ -141 +153 @@
 
   for (i = 0; i < Ncatalog; i++) {
-    free (bin[i]);
-  }
-  free (bin);
+    free (MeasureToImage[i]);
+  }
+  free (MeasureToImage);
   for (i = 0; i < Nimage; i++) {
-    if (clist[i]) { free (clist[i]); }
-    if (mlist[i]) { free (mlist[i]); }
-  }
-  free (clist);
-  free (mlist);
-  free (Nlist);
-  free (NLIST);
+    if (ImageToCatalog[i]) { free (ImageToCatalog[i]); }
+    if (ImageToMeasure[i]) { free (ImageToMeasure[i]); }
+  }
+  free (ImageToCatalog);
+  free (ImageToMeasure);
+  free (N_onImage);
+  free (N_ONIMAGE);
 }
 
@@ -169 +181 @@
   }
 
-  for (i = 0; VERBOSE2 && (i < Nimage); i++) {
-    name = GetPhotcodeNamebyCode (image[i].photcode);
-    fprintf (stderr, "image "OFF_T_FMT" has "OFF_T_FMT" measures (%s, %s)\n",  i,  Nlist[i],
-             ohana_sec_to_date(image[i].tzero), name);
-  }
+  // watch for chips with few stars
+  int Nfew = 0;
+  int Nbad = 0;
+  for (i = 0; i < Nimage; i++) {
+    // ignore the PHU entries
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+    if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) {
+      name = GetPhotcodeNamebyCode (image[i].photcode);
+      fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)",  i,  image[i].imageID, image[i].name, N_onImage[i], image[i].nstar,
+               ohana_sec_to_date(image[i].tzero), name);
+      if (N_onImage[i] < 20) {
+        fprintf (stderr, "*");
+        Nfew ++;
+      } 
+      if (N_onImage[i] < 15) {
+        fprintf (stderr, "*** warning ***");
+        Nbad ++;
+      } 
+      fprintf (stderr, "\n");
+    }
+  }
+  fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad);
 }
 
@@ -181 +210 @@
 
   off_t idx, ID;
-  Measure *measure;
+  MeasureTiny *measure;
   int i, found;
 
-  measure = &catalog[cat].measure[meas];
+  measure = &catalog[cat].measureT[meas];
 
   ID = measure[0].imageID;
@@ -194 +223 @@
 
   // index for (catalog, measure) -> image
-  bin[cat][meas] = idx;
+  MeasureToImage[cat][meas] = idx;
 
   // if we need to allocate an image index table, do so here
-  if (!clist[idx]) {
-      ALLOCATE (clist[idx], int, NLIST[idx]);
-  }
-  if (!mlist[idx]) {
-      ALLOCATE (mlist[idx], off_t, NLIST[idx]);
+  if (!ImageToCatalog[idx]) {
+      ALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]);
+  }
+  if (!ImageToMeasure[idx]) {
+      ALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]);
   }
 
   // index for image, Nentry -> catalog
-  clist[idx][Nlist[idx]] = cat;
+  ImageToCatalog[idx][N_onImage[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]);
+  ImageToMeasure[idx][N_onImage[idx]] = meas;
+  N_onImage[idx] ++;
+
+  if (N_onImage[idx] == N_ONIMAGE[idx]) {
+    N_ONIMAGE[idx] += 30;
+    REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]);
+    REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]);
   }
 
@@ -241 +270 @@
 
   off_t i;
-  Measure *measure;
-
-  measure = &catalog[cat].measure[meas];
+  MeasureTiny *measure;
+
+  measure = &catalog[cat].measureT[meas];
 
   /* find the image that supplied this measurement */
@@ -255 +284 @@
 
     // index for (catalog, measure) -> image
-    bin[cat][meas] = i;
+    MeasureToImage[cat][meas] = i;
 
     // index for image, Nentry -> catalog
-    clist[i][Nlist[i]] = cat;
+    ImageToCatalog[i][N_onImage[i]] = cat;
 
     // index for image, Nentry -> measure
-    mlist[i][Nlist[i]] = meas;
-    Nlist[i] ++;
-
-    if (Nlist[i] == NLIST[i]) {
-      NLIST[i] += 100;
-      REALLOCATE (clist[i], int,   NLIST[i]);
-      REALLOCATE (mlist[i], off_t, NLIST[i]);
+    ImageToMeasure[i][N_onImage[i]] = meas;
+    N_onImage[i] ++;
+
+    if (N_onImage[i] == N_ONIMAGE[i]) {
+      N_ONIMAGE[i] += 30;
+      REALLOCATE (ImageToCatalog[i], IDX_T, N_ONIMAGE[i]);
+      REALLOCATE (ImageToMeasure[i], IDX_T, N_ONIMAGE[i]);
     }
     return;
@@ -281 +310 @@
   off_t i;
 
-  i = bin[cat][meas];
+  i = MeasureToImage[cat][meas];
   if (i == -1) return (NULL);
   return (&image[i].coords);
@@ -333 +362 @@
 }
 
-void dump_measures(Average *average, Measure *measure) {
+void dump_measures (Average *average, Measure *measure) {
 
   off_t j, off;
@@ -402 +431 @@
   }      
 
-  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;
+  int NoffRAave = 0;  int NoffRAori = 0; 
+  int NoffDECave = 0; int NoffDECori = 0;
+
+  for (i = 0; i < N_onImage[im]; i++) {
+    m = ImageToMeasure[im][i];
+    c = ImageToCatalog[im][i];
+
+    Measure *measure = &catalog[c].measure[m];
+
+    X = measure[0].Xccd;
+    Y = measure[0].Yccd;
+    n = measure[0].averef;
 
     if (moscoords == NULL) {
@@ -428 +462 @@
     // NOTE: This should never happen, or our StarMap tests are not working
     if (fabs(dR) > 3.0*ADDSTAR_RADIUS) {
-      fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
-      dump_measures (&catalog[c].average[n], catalog[c].measure);
+      NoffRAave ++;
+      if (VERBOSE2) {
+        fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
+        dump_measures (&catalog[c].average[n], catalog[c].measure);
+      }
       // abort ();
     }
     if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
-      fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
-      dump_measures (&catalog[c].average[n], catalog[c].measure);
+      NoffDECave ++;
+      if (VERBOSE2) {
+        fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
+        dump_measures (&catalog[c].average[n], catalog[c].measure);
+      }
       // abort ();
     }
 
     // complain if the new location is far from the old location
-    if (fabs(catalog[c].measure[m].dR - dR) > DPOS_MAX_ASEC) {
-      fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
-      dump_measures (&catalog[c].average[n], catalog[c].measure);
+    if (fabs(catalog[c].measureT[m].dR - dR) > DPOS_MAX_ASEC) {
+      NoffRAori ++;
+      if (VERBOSE2) {
+        fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
+        dump_measures (&catalog[c].average[n], catalog[c].measure);
+      }
       // abort();
     }
-    if (fabs(catalog[c].measure[m].dD - dD) > DPOS_MAX_ASEC) {
-      fprintf (stderr, "measurement is far from original location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
-      dump_measures (&catalog[c].average[n], catalog[c].measure);
-      // abort();
-    }
-
-    dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD);
+    if (fabs(measure[0].dD - dD) > DPOS_MAX_ASEC) {
+      NoffDECori ++;
+      if (VERBOSE2) {
+        fprintf (stderr, "measurement is far from original location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD);
+        dump_measures (&catalog[c].average[n], catalog[c].measure);
+      }      // abort();
+    }
+
+    dPos += SQ(measure[0].dR - dR) + SQ(measure[0].dD - dD);
     nPos ++;
 
-    catalog[c].measure[m].dR = dR;
-    catalog[c].measure[m].dD = dD;
+    measure[0].dR = dR;
+    measure[0].dD = dD;
     
-    if (catalog[c].measure[m].dR > +180.0*3600.0) {
+    if (measure[0].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) {
+      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+    if (measure[0].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);
+      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
     }
 
     // set the systematic error for this image:
-    catalog[c].measure[m].dRsys = ToShortPixels(dPosSys);
-  }
-
+    measure[0].dRsys = ToShortPixels(dPosSys);
+  }
+
+  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
+  if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);
   saveOffsets (dPos, nPos, im);
 
@@ -487 +534 @@
   Mosaic *mosaic;
   Coords *moscoords, *imcoords, *oldcoords;
+
+  fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n");
+  abort();
 
   // check if this image is bad and should be skipped
@@ -507 +557 @@
   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;
+  for (i = 0; i < N_onImage[im]; i++) {
+    m = ImageToMeasure[im][i];
+    c = ImageToCatalog[im][i];
+
+    // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated
+    Measure *measure = &catalog[c].measure[m];
+
+    X = measure[0].Xccd;
+    Y = measure[0].Yccd;
+    n = measure[0].averef;
 
     dR = dD = 0.0;
@@ -527 +580 @@
     }
 
-    catalog[c].measure[m].dR = dR;
-    catalog[c].measure[m].dD = dD;
-
-    if (catalog[c].measure[m].dR > +180.0*3600.0) {
+    measure[0].dR = dR;
+    measure[0].dD = dD;
+
+    if (measure[0].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) {
+      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+    if (measure[0].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);
+      measure[0].dR = 3600.0*(catalog[c].average[n].R - R);
     }
   }
@@ -567 +620 @@
   }
 
-  ALLOCATE (raw, StarData, Nlist[im]);
-
-  for (i = 0; i < Nlist[im]; i++) {
-    m = mlist[im][i];
-    c = clist[im][i];
+  ALLOCATE (raw, StarData, N_onImage[im]);
+
+  for (i = 0; i < N_onImage[im]; i++) {
+    m = ImageToMeasure[im][i];
+    c = ImageToCatalog[im][i];
+
+    MeasureTiny *measure = &catalog[c].measureT[m];
 
     /* apply the current image transformation or use the current value of R+dR, D+dD? */
-    raw[i].X = catalog[c].measure[m].Xccd;
-    raw[i].Y = catalog[c].measure[m].Yccd;
-
-    raw[i].Mag  = catalog[c].measure[m].M;
-    raw[i].dMag = catalog[c].measure[m].dM;
-    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
-
-    n = catalog[c].measure[m].averef;
+    raw[i].X = measure[0].Xccd;
+    raw[i].Y = measure[0].Yccd;
+
+    raw[i].Mag  = measure[0].M;
+    raw[i].dMag = measure[0].dM;
+    raw[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS);
+
+    n = measure[0].averef;
 
     // an object with only one detection provides no information about the image calibration
     // XXX this is already taken care of in bcatalog
-    raw[i].mask = 0x0000;
+    raw[i].mask = MARK_MEAS_DEFAULT;
     if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
-      raw[i].mask |= 0x0001;
-    }
-    if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
-      raw[i].mask |= 0x0002;
+      raw[i].mask |= MARK_TOO_FEW_MEAS;
+    }
+    if (!finite(measure[0].dR) || !finite(measure[0].dD)) {
+      raw[i].mask |= MARK_NAN_POS_ERROR;
     }
     raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is
@@ -613 +668 @@
   }
 
-  *Nstars = Nlist[im];
+  *Nstars = N_onImage[im];
   return (raw);
 }
@@ -639 +694 @@
   }
 
-  ALLOCATE (ref, StarData, Nlist[im]);
-
-  for (i = 0; i < Nlist[im]; i++) {
-    m = mlist[im][i];
-    c = clist[im][i];
-    n = catalog[c].measure[m].averef;
+  ALLOCATE (ref, StarData, N_onImage[im]);
+
+  for (i = 0; i < N_onImage[im]; i++) {
+    m = ImageToMeasure[im][i];
+    c = ImageToCatalog[im][i];
+
+    MeasureTiny *measure = &catalog[c].measureT[m];
+    n = measure[0].averef;
 
     /* apply the current image transformation or use the current value of R+dR, D+dD? */
@@ -650 +707 @@
     ref[i].D = catalog[c].average[n].D;
 
-    ref[i].Mag  = catalog[c].measure[m].M;
-    ref[i].dMag = catalog[c].measure[m].dM;
-    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
+    ref[i].Mag  = measure[0].M;
+    ref[i].dMag = measure[0].dM;
+    ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS);
 
     ref[i].mask = FALSE;
@@ -675 +732 @@
   }
 
-  *Nstars = Nlist[im];
+  *Nstars = N_onImage[im];
   return (ref);
 }
@@ -682 +739 @@
 void FlagOutliers2D(Catalog *catalog);
 
+// operates on Full values (not tiny)
 void FlagOutliers (Catalog *catalog) {
-  // XXX FlagOutliers is just using FlagOutliers2D
+
+  // XXX FlagOutliers is now just using FlagOutliers2D
   FlagOutliers2D(catalog);
   return;
@@ -696 +755 @@
   assert(catalog[0].Nsecfilt == Nsecfilt);
 
-  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+  if (VERBOSE2) fprintf (stderr, "marking poor measures\n");
   Nmax = 0;
   for (i = 0; i < catalog[0].Naverage; i++) {
@@ -716 +775 @@
   for (j = 0; j < catalog[0].Naverage; j++) {
     
+    // pointer to this set of measurements
+    m = catalog[0].average[j].measureOffset;
+    Measure *measure = &catalog[0].measure[m];
+
     /* accumulate list of valid measurements */
-    m = catalog[0].average[j].measureOffset;
     N = 0;
-    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
-      
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
-      R[N] = catalog[0].measure[m].dR;
-      D[N] = catalog[0].measure[m].dD;
-      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
-      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (!MeasFilterTest(&measure[k], FALSE)) continue;
+      R[N] = measure[k].dR;
+      D[N] = measure[k].dD;
+      dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA);
+      dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC);
       if (isnan(R[N]) || isnan(D[N])) continue;
       N++;
@@ -741 +802 @@
     
     /* compare per-object distance to this standard deviation, and flag outliers*/
-    m = catalog[0].average[j].measureOffset;
     N = 0;
-    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
       //reset flag on each invocation
-      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
 
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
+      if (!MeasFilterTest(&measure[k], FALSE)) continue;
       
-      x = catalog[0].measure[m].dR - statsR.median;
-      y = catalog[0].measure[m].dD - statsD.median;
+      x = measure[k].dR - statsR.median;
+      y = measure[k].dD - statsD.median;
       theta = atan2(y,x);
       if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + 
                          SQR(statsD.sigma * Ns * sin(theta)))) {   
-        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
         Ndel++;
       }
@@ -787 +847 @@
   assert(catalog[0].Nsecfilt == Nsecfilt);
 
-  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+  if (VERBOSE2) fprintf (stderr, "marking poor measures\n");
   Nmax = 0;
   for (i = 0; i < catalog[0].Naverage; i++) {
@@ -809 +869 @@
   for (j = 0; j < catalog[0].Naverage; j++) {
     
+    // pointer to this set of measurements
+    m = catalog[0].average[j].measureOffset;
+    Measure *measure = &catalog[0].measure[m];
+
     /* accumulate list of valid measurements */
-    m = catalog[0].average[j].measureOffset;
     N = 0;
-    
-    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
-      
-      //reset flag on each invocation
-      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
+
+      // reset flag on each invocation
+      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
       
       // skip measurements based on user selected criteria
-      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;
-      R[N] = catalog[0].measure[m].dR;
-      D[N] = catalog[0].measure[m].dD;
-      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
-      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (!MeasFilterTest(&measure[k], FALSE)) continue;
+      R[N] = measure[k].dR;
+      D[N] = measure[k].dD;
+      dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA);
+      dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC);
       if (isnan(R[N]) || isnan(D[N])) continue;
       N++;
@@ -837 +899 @@
     
     /* calculate deviations of all points*/
-    m = catalog[0].average[j].measureOffset;
     N = 0;
-    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
-      //skip bad measurements
-      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;  
-      x = catalog[0].measure[m].dR - statsR.median;
-      y = catalog[0].measure[m].dD - statsD.median;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
+      // skip bad measurements
+      if (!MeasFilterTest(&measure[k], FALSE)) continue;  
+      x = measure[k].dR - statsR.median;
+      y = measure[k].dD - statsD.median;
       theta = atan2(y,x);
       d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
                              SQR(statsD.sigma * Ns * sin(theta)));      
-      index[N] = m;
+      index[N] = k;
       N++;
     }
@@ -856 +917 @@
 
     // recalculate image center, sigma based on closest 50% of points
-    for(k = 0;  k < N; k++) {
+    for (k = 0;  k < N; k++) {
       off_t ind = (off_t) index[k];
-      R[k] = catalog[0].measure[ind].dR;
-      D[k] = catalog[0].measure[ind].dD;
-      dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
-      dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
+      R[k] = measure[ind].dR;
+      D[k] = measure[ind].dD;
+      dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA);
+      dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC);
     }
     liststats (R, dR, N, &statsR);
@@ -869 +930 @@
     
     // use these new statistics to flag outliers 
-    m = catalog[0].average[j].measureOffset;
     N = 0;
-    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
       //skip bad measurements
-      if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue;  
-      x = catalog[0].measure[m].dR - statsR.median;
-      y = catalog[0].measure[m].dD - statsD.median;
+      if (!MeasFilterTest(&measure[k], FALSE)) continue;  
+      x = measure[k].dR - statsR.median;
+      y = measure[k].dD - statsD.median;
       theta = atan2(y,x);
       d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
                              SQR(statsD.sigma * Ns * sin(theta)));      
       if ((d2[N]) > 1) {
-        catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+        measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
         Ndel ++;
       }
       N++;
       Nave++;
-    }  //done rejecting outliers
+    }  // done rejecting outliers
 
     // examine results
     // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
     
-  } //done looping over objects
+  } // done looping over objects
   
   if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
@@ -901 +961 @@
 }
 
-
+/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
+// we only optionally apply the sigma limit: for object averages, this should not be used (should it?)
+int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) {
+  int found, k;
+  long mask;
+  PhotCode *code;
+  float mag;
+
+  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
+  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
+  if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks?
+  
+  /* select measurements by photcode, or equiv photcode, if specified */
+  if (NphotcodesKeep > 0) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
+      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (!found) return FALSE;
+  }
+  
+  if (NphotcodesSkip > 0) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
+      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (found) return FALSE;
+  }  
+  
+  /* select measurements by time */
+  if (TimeSelect) {
+    if (measure[0].t < TSTART) return FALSE;
+    if (measure[0].t > TSTOP) return FALSE;
+  }
+  
+  /* select measurements by quality */
+  if (PhotFlagSelect) {
+    if (PhotFlagBad) {
+      mask = PhotFlagBad;
+    } else {
+      code = GetPhotcodebyCode (measure[0].photcode);
+      mask = code[0].astromBadMask;
+    }
+    if (mask & measure[0].photFlags) return FALSE;
+  }
+
+  /* select measurements by measurement error */
+  if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
+    return FALSE;
+  }
+  
+  /* select measurements by mag limit */
+  if (ImagSelect) {
+    mag = PhotInstTiny (measure);
+    if (mag < ImagMin || mag > ImagMax) return FALSE;
+  }
+  
+  return TRUE;
+}
 
 /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 

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

@@ -1 +1 @@
 # include "relastro.h"
+# define J2000 2451545.       /* Julian date at standard epoch */
 
 # if (0)
@@ -31 +32 @@
   double n, L, g;
 
-# define J2000 2451545.       /* Julian date at standard epoch */
 
   n = jd - J2000;
@@ -46 +46 @@
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) {
 
-# define J2000 2451545.       /* Julian date at standard epoch */
-
   double n = jd - J2000;              // day number
   double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration)
@@ -60 +58 @@
 
 /* given RA, DEC, Time, calculate the parallax factor */
+// Time is relative to Tmean, Tmean is years relative to J2000
 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) {
 
@@ -66 +65 @@
   /* given a time T in UNIX seconds, determine the solar longitude S */
 
-  jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
+  // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean));
+  jd = 365.25*(Time + Tmean) + J2000;
+  // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd);
+
   sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius);
 

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

@@ -68 +68 @@
   for (i = 0; i < catalog[0].Nmeasure; i++) {
     
-    N = getImageByID(catalog[0].measure[i].imageID);
+    MeasureTiny *measure = &catalog[0].measureT[i];
+
+    N = getImageByID(measure[0].imageID);
     if (N < 0) continue;
 
-    xbin = catalog[0].measure[i].Xccd / starmap[N].Nx;
-    ybin = catalog[0].measure[i].Yccd / starmap[N].Ny;
+    xbin = measure[0].Xccd / starmap[N].Nx;
+    ybin = measure[0].Yccd / starmap[N].Ny;
 
     xbin = MAX(0, MIN(NX_MAP-1, xbin));
@@ -82 +84 @@
 
   return (TRUE);
+}
+
+int createStarMap (Catalog *catalog, int Ncatalog) {
+
+  int i;
+
+  initStarMaps();
+
+  for (i = 0; i < Ncatalog; i++) {
+    // check coverage per chip (operates on (Average, MeasureTiny, Secfilt)
+    updateStarMaps (&catalog[i]);
+  }
+
+  createStarMapPoints();
+  return TRUE;
 }
 

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

@@ -1 +1 @@
 # include "relastro.h"
+
 int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
 int saveCenter (Image *image, double *Ro, double *Do, int im);
-
+off_t getNextImageForThread ();
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog);
+void *UpdateChips_worker (void *data);
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  double *Ro;
+  double *Do;
+  char *mode;
+  off_t Nskip;
+  off_t Nmosaic;
+  off_t NnewFit;
+  off_t NoldFit;
+  Catalog *catalog;
+  int Ncatalog;
+} ThreadInfo;
+
+static Image *image = NULL;
+static off_t Nimage = 0;
+static off_t nextImage = 0;
+
+// update astrometry of all chips relative to the average positions
+// if NTHREADS is non-zero, call the threaded version of this function
 int UpdateChips (Catalog *catalog, int Ncatalog) {
 
-  int Nskip, Nmosaic, NnewFit, NoldFit;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref, nFitAstr;
-  Image *image;
+  off_t i, Nraw, Nref, nFitAstr;
   StarData *raw, *ref;
   Coords *oldCoords;
@@ -16 +41 @@
   char *mode;
 
+  if (NTHREADS) {
+    UpdateChips_threaded (catalog, Ncatalog);
+    return TRUE;
+  }
+
   Nskip = Nmosaic = NnewFit = NoldFit = 0;
 
@@ -25 +55 @@
   ALLOCATE (mode, char, Nimage);
 
+  // XXX for faster processing in the future, this can be easily run in parallel
+  // each chip is fitted independently, so we could do N at once in parallel
   for (i = 0; i < Nimage; i++) {
+
+    // XXX looks like everything below is thread safe : we can unroll this into a set of
+    // helper functions that grab the next available chip....
 
     /* skip all except WRP images */
@@ -47 +82 @@
       Nskip ++;
       mode[i] = 0;
+      free (raw);
       continue;
     }
@@ -62 +98 @@
     // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
     if (!FitChip (raw, ref, Nraw, &image[i])) {
-      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
 
       // restore status quo ante
@@ -105 +141 @@
 
   plotChipFits (Ro, Do, mode, Nimage);
-
-  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+           (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
   return (TRUE);
 }
 
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog) {
+
+  int i;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
+  double *Ro, *Do;
+  char *mode;
+
+  image = getimages (&Nimage, NULL);
+  nextImage = 0;
+
+  // save fit results for summary plot
+  ALLOCATE (Ro, double, Nimage);
+  ALLOCATE (Do, double, Nimage);
+  ALLOCATE (mode, char, Nimage);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].Ro = Ro;
+    threadinfo[i].Do = Do;
+    threadinfo[i].mode = mode;
+    threadinfo[i].Nskip = 0;
+    threadinfo[i].Nmosaic = 0;
+    threadinfo[i].NnewFit = 0;
+    threadinfo[i].NoldFit = 0;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].Ncatalog = Ncatalog;
+    pthread_create (&threads[i], NULL, UpdateChips_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  plotChipFits (Ro, Do, mode, Nimage);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  // report stats & summary from the threads
+  Nskip = Nmosaic = NnewFit = NoldFit = 0;
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "UpdateChips thread %d : %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+             i, (int) threadinfo[i].NnewFit, (int) threadinfo[i].NoldFit, (int) threadinfo[i].Nskip, (int) threadinfo[i].Nmosaic);
+
+    NnewFit += threadinfo[i].NnewFit;
+    NoldFit += threadinfo[i].NoldFit;
+    Nskip   += threadinfo[i].Nskip;
+    Nmosaic += threadinfo[i].Nmosaic;
+  }
+  free (threadinfo);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+           (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
+  return (TRUE);
+}
+
+void *UpdateChips_worker (void *data) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  off_t Nraw, Nref, nFitAstr;
+  StarData *raw, *ref;
+  Coords *oldCoords;
+  float dXpixSys, dYpixSys;
+
+  ThreadInfo *threadinfo = data;
+
+  while (1) {
+
+    off_t i = getNextImageForThread();
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
+      threadinfo->Nmosaic ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, &Nraw, MODE_MOSAIC);
+    if (!raw) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (threadinfo->catalog, threadinfo->Ncatalog, i, &Nref, MODE_MOSAIC);
+    if (!ref) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      free (raw);
+      continue;
+    }
+
+    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
+    assert (Nraw == Nref);
+
+    // save these in case of failure
+    saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
+
+    // FitChip does iterative, clipped fitting
+    // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
+    if (!FitChip (raw, ref, Nraw, &image[i])) {
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
+
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 1;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    }
+
+    if (!checkStarMap (i)) {
+      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 2;
+      image[i].flags |= ID_IMAGE_ASTROM_POOR;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    } 
+
+    saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+    threadinfo->mode[i] = 3;
+    threadinfo->NnewFit ++;
+    free (raw);
+    free (ref);
+  }
+  
+  // we should never reach here...
+  return NULL;
+}
+
+// mutex to lock UpdateChips_worker operations 
+static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+// we have an array of chips (image, Nimage).  we need to hand out images one at a time to
+// the worker threads as they need
+off_t getNextImageForThread () {
+
+  pthread_mutex_lock (&UpdateChips_mutex);
+  if (nextImage >= Nimage) {
+    pthread_mutex_unlock (&UpdateChips_mutex);
+    return (-1);
+  }
+  off_t thisImage = nextImage;
+  nextImage ++;
+
+  pthread_mutex_unlock (&UpdateChips_mutex);
+  return (thisImage);
+}
+
+// This function uses the mosaic array in MosaicOps.c to associate PHU + CHIP this is
+// thread safe : the lookup (getMosaicForImage) is selecting an element of a previously
+// allocated and assigned static array.  This function does NOT use the thread-unsafe (and
+// somewhat slower) functions in libdvo/src/mosaic_astrom.c to associate CHIP and PHU
 int saveCenter (Image *image, double *Ro, double *Do, int im) {
 
@@ -136 +379 @@
     // 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);
+    LM_to_RD (&R, &D, L, M, imcoords); // because of the block above, ctype is not -WRP 
   } 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);
-  }
-
-  double Rmid;
-  if (UserCatalog) {
-      Rmid = UserCatalogRA;
-  } else {
-      Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
-  }
+    LM_to_RD (&R, &D, P, Q, moscoords); // moscoords.ctype is -DIS
+  }
+
+  double Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
 
   R = ohana_normalize_angle_to_midpoint (R, Rmid);

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

@@ -1 +1 @@
 # include "relastro.h"
 
+// this function operates on Measure, not MeasureTiny
 int UpdateMeasures (Catalog *catalog, int Ncatalog) {
 

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

@@ -7 +7 @@
 //   original coordinates
 
-int UpdateObjectOffsets (SkyList *skylist) {
+int UpdateObjectOffsets_parallel (SkyList *sky);
+
+int UpdateObjectOffsets (SkyList *skylist, int hostID, char *hostpath) {
 
   int i;
   Catalog catalog;
 
+  // XXX need to decide how to determine PARALLEL mode...
+  if (PARALLEL && !hostID) {
+    UpdateObjectOffsets_parallel (skylist);
+    return TRUE;
+  }
+
   // load data from each region file, only use bright stars
   for (i = 0; i < skylist[0].Nregions; i++) {
 
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
+
+    // define the catalog file name
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
+    catalog.filename = hostID ? hostfile : skylist[0].filename[i];
+
     // set up the basic catalog info
-    catalog.filename  = skylist[0].filename[i];
     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
@@ -22 +37 @@
     catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
 
-    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "w")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
       exit (1);
     }
     if (!catalog.Naves_disk) {
-      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
       dvo_catalog_unlock (&catalog);
       dvo_catalog_free (&catalog);
       continue;
     }
+
+    // set the values in MeasureTiny needed by UpdateObjects
+    populate_tiny_values(&catalog, DVO_TV_MEASURE);
 
     // match measurements with images
@@ -40 +58 @@
     UpdateMeasures (&catalog, 1);
 
+    // uses MeasureTiny values, but will also update Measure values when populated
     UpdateObjects (&catalog, 1);
 
@@ -50 +69 @@
   return (TRUE);
 }
+
+// CATDIR is supplied globally
+# define DEBUG 1
+int UpdateObjectOffsets_parallel (SkyList *sky) {
+
+  // launch the setphot_client jobs to the parallel hosts
+
+  // load the list of hosts
+  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
+  if (!table) {
+    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
+    exit (1);
+  }    
+
+  int i;
+  for (i = 0; i < table->Nhosts; i++) {
+
+    // ensure that the paths are absolute path names
+    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
+    free (table->hosts[i].pathname);
+    table->hosts[i].pathname = tmppath;
+
+    // options / arguments that can affect relastro_client -load:
+    // VERBOSE, VERBOSE2
+    // RESET (-reset)
+    // TimeSelect -time
+    // (note that psfQual is applied rigidly at 0.85, as is the galaxy test)
+    // ImagSelect, ImagMin, ImagMax
+    // MaxDensityUse, MaxDensityValue
+
+    char command[1024];
+    snprintf (command, 1024, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", 
+              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR);
+
+    char tmpline[1024];
+    if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);           strcpy (command, tmpline); }
+
+    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
+    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
+    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
+    if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
+    if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
+    
+    if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command);              strcpy (command, tmpline); }
+    if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
+    
+    if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
+    if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
+    if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                     strcpy (command, tmpline); }
+    if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);     strcpy (command, tmpline); }
+    if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);   strcpy (command, tmpline); }
+    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+
+    if (TimeSelect) { 
+      char *tstart = ohana_sec_to_date (TSTART);
+      char *tstop  = ohana_sec_to_date (TSTOP);
+      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); 
+      free (tstart);
+      free (tstop);
+      strcpy (command, tmpline); 
+    }
+    fprintf (stderr, "command: %s\n", command);
+
+    if (PARALLEL_MANUAL) continue;
+
+    if (PARALLEL_SERIAL) {
+      int status = system (command);
+      if (status) {
+        fprintf (stderr, "ERROR running relastro_client\n");
+        exit (2);
+      }
+    } else {
+      // launch the job on the remote machine (no handshake)
+      int errorInfo = 0;
+      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
+      if (!pid) {
+        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
+        exit (1);
+      }
+      table->hosts[i].pid = pid; // save for future reference
+    }
+  }
+
+  if (PARALLEL_MANUAL) {
+    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
+    getchar();
+  }
+  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
+    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
+  }
+
+  return TRUE;
+}      

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

@@ -38 +38 @@
 }  
 
+// This function operates on both Measure and MeasureTiny.  In the big stages, this should
+// be called with just MeasureTiny set and Measure == NULL
 int UpdateObjects (Catalog *catalog, int Ncatalog) {
 
@@ -72 +74 @@
   T2000 = ohana_date_to_sec ("2000/01/01");
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+
   Nsecfilt = GetPhotcodeNsecfilt ();
-  assert (catalog[0].Nsecfilt == Nsecfilt);
+  if (Ncatalog) {
+    assert (catalog[0].Nsecfilt == Nsecfilt);
+  }
 
   NaveSum = NparSum = NpmSum = NskipSum = 0;
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
+    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);
 
     Nave = Npar = Npm = Nskip = 0;
@@ -86 +91 @@
       XVERB = FALSE;
 
-      // skip objects which are known to be problematic
-      // XXX include this code or not?
-# if (0)
-      if (catalog[i].average[j].code & STAR_BAD) {
-        Nskip ++;
-        continue;  
-      }
-# endif
-
       if (catalog[i].average[j].Nmeasure == 0) {
           continue;
@@ -101 +97 @@
       N = 0;
       m = catalog[i].average[j].measureOffset;
-      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
+      MeasureTiny *measure = &catalog[i].measureT[m];
+      Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL;
+      // when we update the output measure values, we need to do it here
+
       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++) {
-
-        //does the measurement pass the supplied filtering constraints?
-        if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) {
-          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+
+        // does the measurement pass the supplied filtering constraints?
+        if (!MeasFilterTestTiny(&measure[k], FALSE)) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
 
         //outlier rejection
-        if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
-          catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+        if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
 
-        // exclude measurements by previous outlier detection
-        // XXX include this code or not?
-# if (0)
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) { 
-          catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
+        R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+
+        // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects
+        // to get Time in years relative to UNIX Tzero (1970/01/01)
+        T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years
+
+        // dX, dY : error in arcsec -- 
+        dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
+        dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);
+
+        // allow a given photcode or measurement to be
+        // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
+        if (isnan(dX[N])) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
           continue;
         }
-# endif
-
-        catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
-
-        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 - T2000) / (86400*365.25) ; // time relative to J2000 in years
-
-        Tmin = MIN(Tmin, T[N]);
-        Tmax = MAX(Tmax, T[N]);
-        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);
+        if (isnan(dY[N])) {
+          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          continue;
+        }
 
         // add systematic error in quadrature, if desired
         // only do this after the fit has converged (or you will never improve the poor images)
         // if (INCLUDE_SYS_ERR) {
-        // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
+        // float dRsys = FromShortPixels(measure[k].dRsys);
         // dX[N] = hypot(dX[N], dRsys);
         // dY[N] = hypot(dY[N], dRsys);
@@ -154 +156 @@
         // dY[N] = 0.1;
 
-        dT[N] = catalog[i].measure[m].dt;
+        dT[N] = measure[k].dt;
 
         // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
@@ -172 +174 @@
       catalog[i].average[j].flags &= ~ID_STAR_FEW;
 
+      // find Tmin & Tmax from the list of accepted measurements
+      Tmean = 0;
+      Tmin = Tmax = T[0];
+      for (k = 0; k < N; k++) {
+        Tmin = MIN(Tmin, T[k]);
+        Tmax = MAX(Tmax, T[k]);
+        Tmean += T[k];
+      }
       // XXX add the parallax factor range as a criterion as well
       Trange = Tmax - Tmin;
       if (Trange < PM_DT_MIN) mode = FIT_AVERAGE;
       if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
+
+      if (FIT_TARGET == TARGET_HIGH_SPEED) {
+          Tmean = 0.5*(Tmax - Tmin);
+      } else {
+          Tmean /= (float) N;
+      }
 
       // too few measurements for average position (require 2 values)
@@ -192 +208 @@
       coords.crval2 = D[0];
 
-      if (FIT_TARGET == TARGET_HIGH_SPEED) {
-          Tmean = 0.5*(Tmax - Tmin);
-      } else {
-          Tmean /= (float) N;
-      }
-      
       // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
       XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
-      // XVERB = (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
@@ -232 +241 @@
         // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__);
 
+        float pXmin = +2.0;
+        float pXmax = -2.0;
+        float pYmin = +2.0;
+        float pYmax = -2.0;
         for (k = 0; k < N; k++) {
           ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
-        }
-        FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
-        XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
-        catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
-        Npar ++;
+          pXmin = MIN (pXmin, pX[k]);
+          pXmax = MAX (pXmax, pX[k]);
+          pYmin = MIN (pYmin, pY[k]);
+          pYmax = MAX (pYmax, pY[k]);
+        }
+        float dXRange = pXmax - pXmin;
+        float dYRange = pYmax - pYmin;
+        float parRange = hypot (dXRange, dYRange);
+        
+        if (parRange >= PAR_FACTOR_MIN) {
+          FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
+          XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
+          catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
+          Npar ++;
+        }
       }   
 
@@ -307 +330 @@
 
       // 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++) {
-        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]);
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+        setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        if (measureBig) {
+          setMeanR_Big (fit.Ro, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+          setMeanD_Big (fit.Do, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+        }
       }      
 

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

@@ -1 +1 @@
 # include "relastro.h"
 void usage (void);
+void usage_client (void);
 void usage_merge_source (void);
 void usage_merge_source_id (char *name);
@@ -42 +43 @@
   if ((N = get_argument (argc, argv, "-update-objects"))) {
     remove_argument (N, &argc, argv);
-    FIT_TARGET = TARGET_OBJECTS;
+    FIT_TARGET = TARGET_UPDATE_OBJECTS;
+  }
+
+  if ((N = get_argument (argc, argv, "-update-offsets"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_UPDATE_OFFSETS;
   }
 
@@ -92 +98 @@
     RADIUS = atof(argv[N]);
     remove_argument (N, &argc, argv);
-    HIGH_SPEED_DIR = strcreate(argv[N]);
-    remove_argument (N, &argc, argv);
+    HIGH_SPEED_DIR = abspath(argv[N], MAX_PATH_LENGTH);
+    remove_argument (N, &argc, argv);
+  }
+
+  PARALLEL_OUTPUT = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel-output"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL_OUTPUT = TRUE;
+    if (FIT_TARGET != TARGET_HIGH_SPEED) {
+      fprintf (stderr, "-parallel-output only valid for -high-speed mode\n");
+      exit (1);
+    }
   }
 
@@ -118 +134 @@
 
   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? */
@@ -144 +150 @@
     UserPatch.Dmax = atof (argv[N]);
     remove_argument (N, &argc, argv);
-  } else {
-    if (!UserCatalog) {
-      usage ();
+  } 
+  if ((N = get_argument (argc, argv, "-catalog"))) {
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmin = atof (argv[N]);
+    UserPatch.Rmax = UserPatch.Rmin + 0.001;
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmin = atof (argv[N]);
+    UserPatch.Dmax = UserPatch.Dmin + 0.001;
+    remove_argument (N, &argc, argv);
+  }
+
+  // XXX for the moment, make this selection manual.  it needs to be automatic 
+  // based on the state of the SkyTable
+  PARALLEL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel"))) {
+    PARALLEL = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  // this is a test mode : rather than launching the remote jobs and waiting for completion,
+  // relastro will simply list the remote command and wait for the user to signal completion
+  PARALLEL_MANUAL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel-manual"))) {
+    PARALLEL = TRUE; // -parallel-manual implies -parallel
+    PARALLEL_MANUAL = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  // this is a test mode : rather than launching the relastro_client jobs remotely, they are 
+  // run in serial via 'system'
+  PARALLEL_SERIAL = FALSE;
+  if ((N = get_argument (argc, argv, "-parallel-serial"))) {
+    if (PARALLEL_MANUAL) {
+      fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n");
+      exit (1);
     }
+    PARALLEL = TRUE; // -parallel-serial implies -parallel
+    PARALLEL_SERIAL = TRUE;
+    remove_argument (N, &argc, argv);
   }
 
@@ -290 +329 @@
   }  
 
-  AreaSelect = FALSE;
-  if ((N = get_argument (argc, argv, "-area"))) {
-    remove_argument (N, &argc, argv);
-    AreaXmin = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    AreaXmax = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    AreaYmin = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    AreaYmax = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    AreaSelect = TRUE;
-  }
-
   ImagSelect = FALSE;
   if ((N = get_argument (argc, argv, "-instmag"))) {
@@ -336 +361 @@
   }
 
-  /* XXX drop this? */
-  DophotSelect = FALSE;
-  if ((N = get_argument (argc, argv, "-dophot"))) {
-    remove_argument (N, &argc, argv);
-    DophotValue = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    DophotSelect = TRUE;
-  }
-
   NLOOP = 4;
   if ((N = get_argument (argc, argv, "-nloop"))) {
     remove_argument (N, &argc, argv);
     NLOOP = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  NTHREADS = 0;
+  if ((N = get_argument (argc, argv, "-threads"))) {
+    remove_argument (N, &argc, argv);
+    NTHREADS = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -356 +379 @@
 }
 
+int args_client (int argc, char **argv) {
+
+  int N;
+  double trange;
+  char *endptr;
+
+  /* possible operations */
+  FIT_TARGET = TARGET_NONE;
+  FIT_MODE = FIT_AVERAGE;
+
+  // by definition, the client is not parallel 
+  PARALLEL = FALSE;
+  PARALLEL_MANUAL = FALSE;
+  PARALLEL_SERIAL = FALSE;
+
+  BCATALOG = NULL;
+
+  HOST_ID = 0;
+  if ((N = get_argument (argc, argv, "-hostID"))) {
+    remove_argument (N, &argc, argv);
+    HOST_ID = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if (!HOST_ID) usage_client();
+
+  HOSTDIR = NULL;
+  if ((N = get_argument (argc, argv, "-hostdir"))) {
+    remove_argument (N, &argc, argv);
+    HOSTDIR = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if (!HOSTDIR) usage_client();
+
+  if ((N = get_argument (argc, argv, "-load-objects"))) {
+    remove_argument (N, &argc, argv);
+    BCATALOG = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_LOAD_OBJECTS;
+  }
+
+  if ((N = get_argument (argc, argv, "-update-objects"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_UPDATE_OBJECTS;
+  }
+
+  if ((N = get_argument (argc, argv, "-update-offsets"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_UPDATE_OFFSETS;
+  }
+
+  // check for object fitting modes
+  if ((N = get_argument (argc, argv, "-pm"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PM_ONLY;
+  }
+  if ((N = get_argument (argc, argv, "-par"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PAR_ONLY;
+  }
+  if ((N = get_argument (argc, argv, "-pmpar"))) {
+    remove_argument (N, &argc, argv);
+    FIT_MODE = FIT_PM_AND_PAR;
+  }
+
+  if ((N = get_argument (argc, argv, "-high-speed"))) {
+    // XXX include a parallax / no-parallax option
+    if (N >= argc - 3) usage_client();
+    FIT_TARGET = TARGET_HIGH_SPEED;
+    remove_argument (N, &argc, argv);
+    PHOTCODE_A_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    PHOTCODE_B_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    RADIUS = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    HIGH_SPEED_DIR = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj1"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj2"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
+  }
+
+  FlagOutlier = FALSE;
+  if ((N = get_argument (argc, argv, "-clip"))) {
+    remove_argument (N, &argc, argv);
+    CLIP_THRESH = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    FlagOutlier = TRUE;
+  }
+
+  if (FIT_TARGET == TARGET_NONE) usage_client();
+
+  /* specify portion of the sky : allow default of all sky? */
+  UserPatch.Rmin = 0;
+  UserPatch.Rmax = 360;
+  UserPatch.Dmin = -90;
+  UserPatch.Dmax = +90;
+  if ((N = get_argument (argc, argv, "-region"))) {
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  } 
+  if ((N = get_argument (argc, argv, "-catalog"))) {
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmin = atof (argv[N]);
+    UserPatch.Rmax = UserPatch.Rmin + 0.001;
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmin = atof (argv[N]);
+    UserPatch.Dmax = UserPatch.Dmin + 0.001;
+    remove_argument (N, &argc, argv);
+  }
+
+  USE_BASIC_CHECK = FALSE;
+  if ((N = get_argument (argc, argv, "-basic-image-search"))) {
+    remove_argument (N, &argc, argv);
+    USE_BASIC_CHECK = TRUE;
+  }
+
+  MaxDensityUse = FALSE;
+  if ((N = get_argument (argc, argv, "-max-density"))) {
+    remove_argument (N, &argc, argv);
+    MaxDensityValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    MaxDensityUse = TRUE;
+  }
+
+  /* define time */
+  TimeSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-time"))) {
+    remove_argument (N, &argc, argv);
+    if (!ohana_str_to_time (argv[N], &TSTART)) { 
+      fprintf (stderr, "ERROR: syntax error\n");
+      return (FALSE);
+    }
+    remove_argument (N, &argc, argv);
+    if (!ohana_str_to_dtime (argv[N], &trange)) { 
+      if (!ohana_str_to_time (argv[N], &TSTOP)) { 
+        fprintf (stderr, "ERROR: syntax error\n");
+        return (FALSE);
+      }
+    } else {
+      if (trange < 0) {
+        trange = fabs (trange);
+        TSTOP = TSTART;
+        TSTART -= trange;
+      } else {
+        TSTOP = TSTART + trange;
+      }
+    }
+    remove_argument (N, &argc, argv);
+    TimeSelect = TRUE;
+  }
+
+  PHOTCODE_KEEP_LIST = NULL;
+  if ((N = get_argument (argc, argv, "+photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_KEEP_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PHOTCODE_SKIP_LIST = NULL;
+  if ((N = get_argument (argc, argv, "-photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_SKIP_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  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);
+  }
+
+  // by default, require > 10pts to clip
+  strcpy (STATMODE, "CHI_INNER_80_WTMEAN");
+  if ((N = get_argument (argc, argv, "-statmode"))) {
+    remove_argument (N, &argc, argv);
+    strcpy (STATMODE, argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  RESET = FALSE;
+  if ((N = get_argument (argc, argv, "-reset"))) {
+    remove_argument (N, &argc, argv);
+    RESET = TRUE;
+  }
+
+  UPDATE = FALSE;
+  if ((N = get_argument (argc, argv, "-update"))) {
+    remove_argument (N, &argc, argv);
+    UPDATE = TRUE;
+  }
+
+  // do not hold the lock on the Image table for relastro_client -update-offsets
+  if (FIT_TARGET == TARGET_UPDATE_OFFSETS) {
+    UPDATE = FALSE;
+  }
+
+  /* XXX update this */
+  MIN_ERROR = 0.001;
+  if ((N = get_argument (argc, argv, "-minerror"))) {
+    remove_argument (N, &argc, argv);
+    MIN_ERROR = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    /* require MIN_ERROR > 0 */
+  }  
+
+  ImagSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-instmag"))) {
+    remove_argument (N, &argc, argv);
+    ImagMin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ImagMax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ImagSelect = TRUE;
+  }
+  
+  // for now, make the default to ignore the photflags
+  // XXX make it true by default instead?
+  PhotFlagSelect = FALSE;
+  if ((N = get_argument (argc, argv, "+photflags"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagSelect = TRUE;
+  }
+
+  PhotFlagBad = 0;
+  if ((N = get_argument (argc, argv, "-photflagbad"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagBad = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PhotFlagPoor = 0;
+  if ((N = get_argument (argc, argv, "-photflagpoor"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagPoor = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 1) usage_client ();
+  return TRUE;
+}
+
 void usage () {
-  fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
-  fprintf (stderr, "       OR:    relastro -catalog (ra) (dec)\n\n");
-  fprintf (stderr, "       OR:    relastro -high-speed [options]\n\n");
+  fprintf (stderr, "ERROR: USAGE: relastro -update-simple [options]\n");
+  fprintf (stderr, "       OR:    relastro -update-chips [options]\n");
+  fprintf (stderr, "       OR:    relastro -update-mosaic [options]\n");
+  fprintf (stderr, "       OR:    relastro -update-objects [options]\n");
+  fprintf (stderr, "       OR:    relastro -high-speed [options]\n");
   fprintf (stderr, "       OR:    relastro -merge-source [options]\n\n");
+
   fprintf (stderr, "  specify one of the following modes: \n");
   fprintf (stderr, "  -update-objects\n");
@@ -369 +667 @@
   fprintf (stderr, "  -update-chips\n");
   fprintf (stderr, "  -update-mosaics\n");
-  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n\n");
-  fprintf (stderr, "  additional options: \n");
+  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n");
+  fprintf (stderr, "  -merge-source (objID) (catID) into (objID) (catID)\n\n");
+
+  fprintf (stderr, " additional options: \n");
+  fprintf (stderr, "  -region RA RA DEC DEC\n");
+  fprintf (stderr, "  -catalog (ra) (dec)\n\n");
   fprintf (stderr, "  -time (start)(stop)\n");
   fprintf (stderr, "  +photcode (code)[,code,code...]\n");
@@ -383 +685 @@
   fprintf (stderr, "  -plrange\n");
   fprintf (stderr, "  -minerror\n");
-  fprintf (stderr, "  -area\n");
-  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
   fprintf (stderr, "  -instmag min max\n\n");
   fprintf (stderr, "  +photflags\n");
@@ -394 +694 @@
 } 
 
+void usage_client () {
+  fprintf (stderr, "ERROR: USAGE: relastro_client -load\n");
+  fprintf (stderr, "       OR:    relastro_client -update-offsets\n");
+  fprintf (stderr, "       OR:    relastro_client -update-objects\n");
+  fprintf (stderr, "       OR:    relastro_client -high-speed\n\n");
+
+  fprintf (stderr, "  specify one of the following modes: \n");
+  fprintf (stderr, "  -load : load the bright source detections for analysis\n");
+  fprintf (stderr, "  -update-offsets : apply the updated image parameters\n");
+  fprintf (stderr, "  -update-objects : calculate average astrometric properties of objects\n");
+  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n");
+  fprintf (stderr, "    support for the special high-speed mode\n");
+
+  fprintf (stderr, " additional options: \n");
+  fprintf (stderr, "  -pm : calculate proper motions\n");
+  fprintf (stderr, "  -par : calculate parallaxes\n");
+  fprintf (stderr, "  -pmpar : calculate motions and parallaxes\n");
+  fprintf (stderr, "  -region Rmin Rmax Dmin Dmax");
+  fprintf (stderr, "  -catalog RA DEC");
+  fprintf (stderr, "  -time (start)(stop)\n");
+  fprintf (stderr, "  +photcode (code)[,code,code...]\n");
+  fprintf (stderr, "  -photcode (code)[,code,code...]\n");
+  fprintf (stderr, "  -statmode (mode)\n");
+  fprintf (stderr, "  -reset\n");
+  fprintf (stderr, "  -minerror\n");
+  fprintf (stderr, "  -instmag min max\n\n");
+  fprintf (stderr, "  +photflags\n");
+  fprintf (stderr, "  -photflagbad\n");
+  fprintf (stderr, "  -photflagpoor\n");
+  fprintf (stderr, "  -v\n");
+  fprintf (stderr, "  \n");
+  exit (2);
+} 
+
 void usage_merge_source_id (char *name) {
 

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

@@ -1 +1 @@
 # include "relastro.h"
+
+static int Nskip1 = 0;
+static int Nskip2 = 0;
+static unsigned int Tref = 1323353985;
+static short Cref = 10355;
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -12 +17 @@
 
   /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
+  memset(subcatalog, 0, sizeof(Catalog));
   NAVERAGE = 50;
   NMEASURE = 1000;
-  ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
-  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
-  ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+  ALLOCATE (subcatalog[0].average,  Average,     NAVERAGE);
+  ALLOCATE (subcatalog[0].secfilt,  SecFilt,     NAVERAGE*Nsecfilt);
+  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
   Nmeasure = Naverage = 0;
+
+  DVOAverageFlags averageBits = 
+    ID_STAR_FEW     |
+    ID_STAR_POOR    |
+    ID_STAR_FIT_AVE |
+    ID_STAR_FIT_PM  |
+    ID_STAR_FIT_PAR |
+    ID_STAR_USE_AVE |
+    ID_STAR_USE_PM  |
+    ID_STAR_USE_PAR; 
 
   /* exclude stars not in range or with too few measurements */
@@ -23 +39 @@
     if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
       continue;
+      // XXX can't (easily) count this for a given chip
     }
     
@@ -33 +50 @@
 
     if (RESET) {
-      // XXX reset the ra,dec coords?
-      // XXX define astrometry vs photometry average.dbFlags
-      // XXX put photometry flags in secfilt (ie, per average mag?)
-      subcatalog[0].average[Naverage].flags &= ~ID_STAR_FEW;
-      subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR;
+      subcatalog[0].average[Naverage].flags &= ~averageBits;
     }
 
@@ -48 +61 @@
       if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+          Nskip1 ++;
+        }
         continue;
       }
@@ -54 +70 @@
       if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
         catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+        if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+          Nskip2 ++;
+        }
         continue;
       }
@@ -70 +89 @@
       // allowed.
 
-      // CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
-
-      subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
-      subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
-      subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
-      subcatalog[0].measure[Nmeasure].averef   = Naverage;
+      CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
+      // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset];
+      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
+      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
+      subcatalog[0].measureT[Nmeasure].averef   = Naverage;
       if (RESET) {
-        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
-        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
       }
 
@@ -85 +103 @@
       if (Nmeasure == NMEASURE) {
         NMEASURE += 1000;
-        REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
       }
     }
@@ -96 +114 @@
     }
   }
-  REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
-  REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
-  REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt);
+  REALLOCATE (subcatalog[0].average,  Average,     MAX (Naverage, 1));
+  REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1));
+  REALLOCATE (subcatalog[0].secfilt,  SecFilt,     MAX (Naverage, 1)*Nsecfilt);
   subcatalog[0].Naverage = Naverage;
   subcatalog[0].Nmeasure = Nmeasure;
   subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
-  subcatalog[0].Nsecf_mem = Naverage * Nsecfilt;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].catID    = catalog[0].catID;
   assert (Nsecfilt == catalog[0].Nsecfilt);
 
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog (subcatalog, catalog);
-  }
-
-  if (VERBOSE) {
-    fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n",  i,  subcatalog[0].Naverage,  subcatalog[0].Nmeasure);
-   }
+    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
+  } else {
+    if (VERBOSE2) {
+      char *basename = filebasename (catalog[0].filename);
+      fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", 
+               subcatalog[0].Naverage, catalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Nmeasure, basename);
+      free (basename);
+    }
+  }
   return (TRUE);
 }
 
+void bcatalog_show_skips () {
+  fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
+}
+
+// sort by decreasing Nmeasure (X)
+void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {
+
+# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
+  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
+  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
+}
+# define COMPARE(A,B)(X[A] > X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
 /* this version does NOT use AverageTiny, MeasureTiny */ 
-int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
 
   Catalog tmpcatalog;
@@ -125 +167 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
-  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
-  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
-  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
 
   if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
@@ -147 +189 @@
   // select a random subset of Nmax stars from subcatalog using Fisher-Yates
 
-  // we are going to select Nmax entries by generating a random-sorted index list
-  off_t *index, tmp, i, j, ave;
+  // we are going to select Nmax entries by choosing the brightest objects
+  int *value;
+  off_t *index, i, j, ave;
   ALLOCATE (index, off_t, Naverage);
+  ALLOCATE (value, int, Naverage);
   for (i = 0; i < Naverage; i++) {
     index[i] = i;
-  }
-  for (i = 0; i < Naverage; i++) {
-    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
-    tmp = index[j];
-    index[j] = index[i];
-    index[i] = tmp;
-  }
+    value[i] = subcatalog[0].average[i].Nmeasure;
+  }
+  sort_by_Nmeasure (value, index, Naverage);
 
   // count the number of measurements this selection will yield
@@ -168 +208 @@
 
   // allocate the output data 
-  ALLOCATE (tmpcatalog.average, Average, Nmax);
-  ALLOCATE (tmpcatalog.measure, Measure, NMEASURE);
-  ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt);
+  ALLOCATE (tmpcatalog.average,  Average,     Nmax);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nmax * Nsecfilt);
 
   off_t Nmeasure = 0;
@@ -181 +221 @@
     for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
       off_t offset = subcatalog[0].average[ave].measureOffset + j;
-      tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset];
-      tmpcatalog.measure[Nmeasure].averef = i;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
       Nmeasure ++;
     }
-  }
-
-  if (VERBOSE) {
+    for (j = 0; j < Nsecfilt; j++) {
+      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
+    }
+  }
+
+  if (VERBOSE2) {
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", 
+             Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
+  }
+
+  free (index);
+  free (value);
+  free (subcatalog[0].average);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].average = tmpcatalog.average;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+
+/* this version does NOT use AverageTiny, MeasureTiny */ 
+int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+
+  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE2) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
+
+  // we are going to select Nmax entries by generating a random-sorted index list
+  off_t *index, tmp, i, j, ave;
+  ALLOCATE (index, off_t, Naverage);
+  for (i = 0; i < Naverage; i++) {
+    index[i] = i;
+  }
+  for (i = 0; i < Naverage; i++) {
+    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
+    tmp = index[j];
+    index[j] = index[i];
+    index[i] = tmp;
+  }
+
+  // count the number of measurements this selection will yield
+  off_t NMEASURE = 0;
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    NMEASURE += subcatalog[0].average[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average,  Average,     Nmax);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nmax * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    tmpcatalog.average[i] = subcatalog[0].average[ave];
+    tmpcatalog.average[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].average[ave].measureOffset + j;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+    for (j = 0; j < Nsecfilt; j++) {
+      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
+    }
+  }
+
+  if (VERBOSE2) {
     fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
              Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
@@ -193 +334 @@
 
   free (subcatalog[0].average);
-  free (subcatalog[0].measure);
+  free (subcatalog[0].measureT);
   free (subcatalog[0].secfilt);
 
   subcatalog[0].average = tmpcatalog.average;
-  subcatalog[0].measure = tmpcatalog.measure;
+  subcatalog[0].measureT = tmpcatalog.measureT;
   subcatalog[0].secfilt = tmpcatalog.secfilt;
   subcatalog[0].Naverage = Nmax;

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

@@ -1 +1 @@
 # include "relastro.h"
 
-double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) {
+double getMeanR (MeasureTiny *measure, Average *average, SecFilt *secfilt) {
 
   double ra;
@@ -21 +21 @@
 }
 
-double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) {
+double getMeanD (MeasureTiny *measure, Average *average, SecFilt *secfilt) {
 
   double dec;
@@ -41 +41 @@
 }
 
-int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt) {
 
   /* math to get from new fitted position to new measure offset
@@ -66 +66 @@
 }
 
-int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt) {
 
   /* math to get from new fitted position to new measure offset
@@ -90 +90 @@
   return (TRUE);
 }
+
+int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  if (!measure) return TRUE;
+
+  /* math to get from new fitted position to new measure offset
+     ra_obs = average[0].R - measure[0].dR / 3600.0;
+     measure[0].dR = (ra_fit - ra_obs) * 3600.0;
+     measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0
+     measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dR += (ra_fit - average[0].R) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}
+
+int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  if (!measure) return TRUE;
+
+  /* math to get from new fitted position to new measure offset
+     dec_obs = average[0].D - measure[0].dD / 3600.0;
+     measure[0].dD = (dec_fit - dec_obs) * 3600.0;
+     measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0
+     measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dD += (dec_fit - average[0].D) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}

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

@@ -149 +149 @@
   
   if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) {
+    array_free (matrix, fit[0].Nelems);
+    array_free (vector, fit[0].Nelems);
     return (FALSE);
   }

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

@@ -1 +1 @@
 # include "relastro.h"
 
-int high_speed_catalogs () {
+int high_speed_catalogs_parallel (SkyList *sky);
+
+int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath) {
 
   int i;
-
-  SkyTable *sky = NULL;
-  SkyList *skylist = NULL;
   Catalog catalog;
 
@@ -12 +11 @@
   dvo_set_catdir(CATDIR);
 
-  // load the current sky table (layout of all SkyRegions) 
-  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE);
-  SkyTableSetFilenames (sky, CATDIR, "cpt");
-  
-  // determine the populated SkyRegions overlapping the requested area (default depth)
-  if (UserCatalog) {
-    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
-    setupAreaSelection(skylist->regions[0]);
-  } else {
-    skylist = SkyListByPatch (sky, -1, &UserPatch);
-    setupAreaSelection(&UserPatch);
-  }
-
   initializeConstraints();
+
+  // XXX need to decide how to determine PARALLEL mode...
+  if (PARALLEL && !hostID) {
+    high_speed_catalogs_parallel (skylist);
+    goto finish;
+  }
 
   // load data from each region file, only use bright stars
   for (i = 0; i < skylist[0].Nregions; i++) {
 
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
+
     // set up the basic catalog info
-    catalog.filename  = skylist[0].filename[i];
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
+    catalog.filename  = hostID ? hostfile : skylist[0].filename[i];
     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
@@ -37 +34 @@
     catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
 
-    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) {
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "r")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
       exit (1);
     }
     if (!catalog.Naves_disk) {
-      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
       dvo_catalog_unlock (&catalog);
       dvo_catalog_free (&catalog);
@@ -54 +51 @@
   }
   
+finish:
+
+  // in relastro (NOT relastr_client), save the Images.dat, Photcodes.dat, SkyTable.fits files
+  if (!hostID) {
+    char photcodeFile[1024];
+    snprintf (photcodeFile, 1024, "%s/Photcodes.dat", HIGH_SPEED_DIR);
+    if (!SavePhotcodesFITS (photcodeFile)) {
+      fprintf (stderr, "error saving photcode table %s\n", photcodeFile);
+      exit (1);
+    }
+    
+    // need to copy across the Images, SkyTable, and Photcode tables:
+    char *skyfile = SkyTableFilename (HIGH_SPEED_DIR);
+    SkyTableSave (sky, skyfile);
+    
+    char line[2048];
+    snprintf (line, 2048, "cp %s/Images.dat %s/Images.dat", CATDIR, HIGH_SPEED_DIR);
+    system (line);
+  }
+
   return (TRUE);
 }
+
+// CATDIR is supplied globally
+# define DEBUG 1
+int high_speed_catalogs_parallel (SkyList *skylist) {
+
+  // launch the setphot_client jobs to the parallel hosts
+
+  // load the list of hosts
+  HostTable *table = HostTableLoad (CATDIR, skylist->hosts);
+  if (!table) {
+    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", skylist->hosts, CATDIR);
+    exit (1);
+  }    
+
+  // If we want to parallelize the output, the remote client needs to load the HostTable
+  // for the output catdir.  Require the output HostTable to match the input (ie, same
+  // hostIDs must exist).  Preferred, but not required, that they match in the host
+  // location.
+  HostTable *tableOut = NULL;
+  if (PARALLEL_OUTPUT) {
+    // load the list of hosts
+    tableOut = HostTableLoad (HIGH_SPEED_DIR, skylist->hosts);
+    if (!tableOut) {
+      fprintf (stderr, "ERROR: failure reading Host Table %s for output database %s\n", skylist->hosts, HIGH_SPEED_DIR);
+      exit (1);
+    }    
+    if (table->Nhosts != tableOut->Nhosts) {
+      fprintf (stderr, "ERROR: output HostTable must have matching hosts (%d in, %d out)\n", table->Nhosts, tableOut->Nhosts);
+      exit (1);
+    }
+  }
+
+  int i;
+  for (i = 0; i < table->Nhosts; i++) {
+
+    // ensure that the paths are absolute path names
+    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
+    free (table->hosts[i].pathname);
+    table->hosts[i].pathname = tmppath;
+
+    // use this directory to save the output catalogs (distributed if parallel)
+    char *outputDir = strcreate(HIGH_SPEED_DIR);
+    if (PARALLEL_OUTPUT) {
+      // find the matching host entry in the output HostTable (match by ID, not index)
+      int hostID = table->hosts[i].hostID;
+      int index = table->index[hostID];
+      if (index == -1) {
+        fprintf (stderr, "ERROR: output HostTable must have matching hosts (host ID %d not found)\n", hostID);
+        exit (1);
+      }
+
+      free (outputDir);
+      outputDir = abspath (tableOut->hosts[index].pathname, MAX_PATH_LENGTH);
+    }
+
+    char command[1024];
+    snprintf (command, 1024, "relastro_client -high-speed %s %s %f %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", 
+              PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+
+    free (outputDir);
+
+    // options / arguments that can affect relastro_client -high-speed
+    char tmpline[1024];
+    if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);              strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);              strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);              strcpy (command, tmpline); }
+
+    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                       strcpy (command, tmpline); }
+    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                       strcpy (command, tmpline); }
+    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                       strcpy (command, tmpline); }
+    if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);     strcpy (command, tmpline); }
+    if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);      strcpy (command, tmpline); }
+    
+    if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command);                 strcpy (command, tmpline); }
+
+    if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);          strcpy (command, tmpline); }
+    
+    if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
+    if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
+    if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                      strcpy (command, tmpline); }
+    if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);      strcpy (command, tmpline); }
+    if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);    strcpy (command, tmpline); }
+    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+
+    if (WHERE_A[0]) { snprintf (tmpline, 1024, "%s -D WHERE_A \"%s\"", command, WHERE_A);               strcpy (command, tmpline); }
+    if (WHERE_B[0]) { snprintf (tmpline, 1024, "%s -D WHERE_B \"%s\"", command, WHERE_B);               strcpy (command, tmpline); }
+
+    if (TimeSelect) { 
+      char *tstart = ohana_sec_to_date (TSTART);
+      char *tstop  = ohana_sec_to_date (TSTOP);
+      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); 
+      free (tstart);
+      free (tstop);
+      strcpy (command, tmpline); 
+    }
+
+    fprintf (stderr, "command: %s\n", command);
+
+    if (PARALLEL_MANUAL) continue;
+
+    if (PARALLEL_SERIAL) {
+      int status = system (command);
+      if (status) {
+        fprintf (stderr, "ERROR running relastro_client\n");
+        exit (2);
+      }
+    } else {
+      // launch the job on the remote machine (no handshake)
+      int errorInfo = 0;
+      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
+      if (!pid) {
+        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
+        exit (1);
+      }
+      table->hosts[i].pid = pid; // save for future reference
+    }
+  }
+
+  if (PARALLEL_MANUAL) {
+    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
+    getchar();
+  }
+  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
+    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
+  }
+
+  return TRUE;
+}      

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

@@ -6 +6 @@
 int high_speed_objects (SkyRegion *region, Catalog *catalog) {
 
-  // XXX seems silly to require region here just to find the catalog bounds / center
-
-  off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE;
-  int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj;
-  int foundA, foundB;
+  off_t i, j, m, J, ni, nj, *N1;
+  off_t Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
+  off_t l, i1, NAVERAGE, NMEASURE;
+  int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas, Nepoch, nv[2], Nmatchmeasobj;
+  int foundA, foundB, XVERB;
   double *X1, *Y1;
   double dX, dY, dR, RADIUS2;
   Coords tcoords;
-  Catalog catalog1;
+  Catalog catalogOut;
 
   int Nsecfilt;
   char filename[1024];
-  Noff = strlen(CATDIR);
-  sprintf (filename, "%s%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]);
-  fprintf(stderr, "%s\n",filename);
-  dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/
-  catalog1.filename = strcreate(filename);
+
+  snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name);
+  fprintf (stderr, "%s\n",filename);
+
+  dvo_catalog_init (&catalogOut, TRUE); /*initialise new catalogue*/
+  catalogOut.filename = strcreate(filename);
 
   SIGMA_LIM = 0.0;
 
   Nsecfilt = GetPhotcodeNsecfilt();
-  Naverage=catalog[0].Naverage;
-  Nmeasure=catalog[0].Nmeasure;
-  catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion
-  catalog1.Nsecfilt  = Nsecfilt;
-  // always load all of the data (if any exists)
-  catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+  // off_t Naverage = catalog[0].Naverage;
+  // off_t Nmeasure = catalog[0].Nmeasure;
+  catalogOut.filename = filename; // based on the input name, need to keep everything below the catdir portion
+  catalogOut.Nsecfilt = Nsecfilt;
+  catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data
   
-  catalog1.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
-  catalog1.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
-
-  if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) {
+  catalogOut.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
+  catalogOut.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
+
+  if (!dvo_catalog_open (&catalogOut, region, VERBOSE2,"w")) {
     fprintf (stderr, "ERROR: failure to open catalog file %s\n",
              filename);
@@ -45 +45 @@
   NAVERAGE = 1000;
   NMEASURE = 10000;
-  REALLOCATE (catalog1.average, Average, NAVERAGE);
-  REALLOCATE (catalog1.measure, Measure, NMEASURE);
-  REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
+  REALLOCATE (catalogOut.average, Average, NAVERAGE);
+  REALLOCATE (catalogOut.measure, Measure, NMEASURE);
+  REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
 
   // high-speed between different surveys (easier case):
@@ -77 +77 @@
   // mark (exclude) objects with both sets of target photcodes
   for (i = 0; i < catalog[0].Naverage; i++) {
+
+    XVERB  = (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC);
+    XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST);
+    if (XVERB) {
+      fprintf (stderr, "test object\n");
+    }
+
     // do any of the measures for this object match group A?
     m = catalog[0].average[i].measureOffset;
@@ -114 +121 @@
             // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n",  i);
         }
-
-    } else if (foundB && !foundA) {
-
+        continue;
+    } 
+    if (foundB && !foundA) {
         if (applyConstraintsB(catalog, i)) {
             groupB[i] = TRUE;
@@ -124 +131 @@
             NgroupBbad++;
         }
+        continue;
     }
     // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow
@@ -173 +181 @@
     nj = N1[j];
 
+    XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
+    XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
+    if (XVERB) {
+      fprintf (stderr, "test object %d or %d\n", (int) ni, (int) nj);
+    }
+
     if (slowMoving[ni]) NEXT_I;
     if (slowMoving[nj]) NEXT_J;
@@ -185 +199 @@
     // look for pairs that are within the maximum separation
     dX = X1[i] - X1[j];
+    if (XVERB) {
+      fprintf (stderr, "%d %d : %f\n", (int) ni, (int) nj, dX);
+    }
 
     if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small
@@ -197 +214 @@
       nj = N1[J];
       if (!groupB[nj]) continue;
+
+      XVERB  = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC);
+      XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST);
+      if (XVERB) {
+        fprintf (stderr, "test object pt2 %d or %d : %f\n", (int) ni, (int) nj, dX);
+        fprintf (stderr, ".\n");
+      }
 
       dY = Y1[i] - Y1[J];
@@ -209 +233 @@
       nv[1]=nj;
 
-      catalog1.average[Nmatch]=catalog[0].average[nv[0]];
+      catalogOut.average[Nmatch]=catalog[0].average[nv[0]];
       /*Loop over index values and set measurements in new catalogue*/
       Nmatchmeasobj=0;
-      catalog1.average[Nmatch].measureOffset=Nmatchmeas;
+      catalogOut.average[Nmatch].measureOffset=Nmatchmeas;
       for(l=0;l<Nepoch;l++) {
           m = catalog[0].average[nv[l]].measureOffset;
           for(i1=0;i1<catalog[0].average[nv[l]].Nmeasure;i1++) {
-              catalog1.measure[Nmatchmeas]=catalog[0].measure[m+i1];
+              catalogOut.measure[Nmatchmeas]=catalog[0].measure[m+i1];
               /*Set offset RA and Dec wrt correct average value*/
-              catalog1.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
-              catalog1.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
-              catalog1.measure[Nmatchmeas].averef = Nmatch;
+              catalogOut.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
+              catalogOut.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
+              catalogOut.measure[Nmatchmeas].averef = Nmatch;
               Nmatchmeasobj++;
               Nmatchmeas++;
@@ -226 +250 @@
               if (Nmatchmeas == NMEASURE - 1) {
                 NMEASURE += 10000;
-                REALLOCATE (catalog1.measure, Measure, NMEASURE);
+                REALLOCATE (catalogOut.measure, Measure, NMEASURE);
               }
           }
       }
-      catalog1.average[Nmatch].Nmeasure=Nmatchmeasobj;
+      catalogOut.average[Nmatch].Nmeasure = Nmatchmeasobj;
       Nmatch ++;
 
       if (Nmatch == NAVERAGE - 1) {
         NAVERAGE += 1000;
-        REALLOCATE (catalog1.average, Average, NAVERAGE);
-        REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
+        REALLOCATE (catalogOut.average, Average, NAVERAGE);
+        REALLOCATE(catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt);
       }
     }
     i++;
   }
-  catalog1.Naverage=Nmatch;
-  catalog1.Nmeasure=Nmatchmeas;
-  catalog1.Nsecfilt=Nsecfilt;
-  catalog1.Nsecf_mem=Nmatch*Nsecfilt;
-
-  UpdateObjects(&catalog1, 1);
+  catalogOut.Naverage=Nmatch;
+  catalogOut.Nmeasure=Nmatchmeas;
+  catalogOut.Nsecfilt=Nsecfilt;
+  catalogOut.Nsecf_mem=Nmatch*Nsecfilt;
+
+  populate_tiny_values (&catalogOut, DVO_TV_MEASURE);
+  UpdateObjects (&catalogOut, 1);
 
   fprintf (stderr, "found %d matches\n", Nmatch);
-  dvo_catalog_save (&catalog1, VERBOSE);
-  dvo_catalog_unlock (&catalog1);
-  dvo_catalog_free (&catalog1);
+  dvo_catalog_save (&catalogOut, VERBOSE2);
+  dvo_catalog_unlock (&catalogOut);
+  dvo_catalog_free (&catalogOut);
   free (slowMoving);
   free (groupA);
@@ -262 +287 @@
 }
 
-// return TRUE if any measure->photcodes match the photcodeSet
+// return TRUE if measure->photcode match any in the photcodeSet
+// (but, return FALSE if the measurement is bad)
 int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) {
 
   int found, k;
 
-  assert (Nset > 0);
-
-  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
-  if (!finite(measure[0].M)) return FALSE;
+  if (!Nset) return TRUE;
+  // XXX require a set or not?  assert (Nset > 0);
+
+  if (!finite(measure[0].dR)) return FALSE;
+  if (!finite(measure[0].dD)) return FALSE;
+  if (!finite(measure[0].M))  return FALSE;
   
+  float dX = GetAstromError (measure, ERROR_MODE_RA);
+  if (isnan(dX)) return FALSE;
+
+  float dY = GetAstromError (measure, ERROR_MODE_DEC);
+  if (isnan(dY)) return FALSE;
+
   /* select measurements by photcode, or equiv photcode, if specified */
-  if (Nset > 0) {
-    found = FALSE;
-    for (k = 0; (k < Nset) && !found; k++) {
-      if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
-      if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
-    }
-    if (!found) return FALSE;
-  }
-  
-  return TRUE;
+  found = FALSE;
+  for (k = 0; (k < Nset) && !found; k++) {
+    if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
+    if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+  }
+
+  return found;
 }

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

@@ -17 +17 @@
 static dbStack *stackA = NULL;
 static int NstackA = 0;
+
 static dbField *fieldsB = NULL;
 static int NfieldsB = 0;
 static dbStack *stackB = NULL;
 static int NstackB = 0;
+
 static dbValue *valuesA;
 // static int NvaluesA;
@@ -41 +43 @@
     ALLOCATE(words, char *, NwordsAllocated);
 
-    char *copy_of_str = strdup(string);
+    char *copy_of_str = strcreate(string);
 
     char *saveptr = NULL;
@@ -60 +62 @@
 // Create the data strucutres required to apply avextract like cuts
 int setupConstraints(char *whereString, dbField **pFields, int *pNfields, dbStack **pStack, int *pNstack) {
+
+  // an empty where is allowed (no further filtering)
+  if (!whereString[0]) return TRUE;
 
     // split the string into words
@@ -89 +94 @@
 
 int initializeConstraints() {
+
   // if constaints haven't been set by the user use 2mass
-    if (!WHERE_A[0]) {
-        strcpy(WHERE_A, DEFAULT_WHERE_A);
-    }
-    if (!WHERE_B[0]) {
-        strcpy(WHERE_B, DEFAULT_WHERE_B);
-    }
+  // XXX  if (!WHERE_A[0]) {
+  // XXX      strcpy(WHERE_A, DEFAULT_WHERE_A);
+  // XXX  }
+  // XXX  if (!WHERE_B[0]) {
+  // XXX      strcpy(WHERE_B, DEFAULT_WHERE_B);
+  // XXX  }
+
+    fprintf (stderr, "where A: %s\n", WHERE_A);
+    fprintf (stderr, "where B: %s\n", WHERE_B);
+
+    // pull TIMEFORMAT & TIMEREF from the Config system
+    dbExtractAveragesInit();
 
     if (! setupConstraints(WHERE_A, &fieldsA, &NfieldsA, &stackA, &NstackA)) {
@@ -102 +114 @@
         exit(1);
     }
-    ALLOCATE(valuesA, dbValue, NfieldsA);
+    if (NfieldsA) {
+      ALLOCATE(valuesA, dbValue, NfieldsA);
+    }
+
     if (! setupConstraints(WHERE_B, &fieldsB, &NfieldsB, &stackB, &NstackB)) {
         print_error();
@@ -108 +123 @@
         exit(1);
     }
-    ALLOCATE(valuesB, dbValue, NfieldsB);
+    if (NfieldsB) {
+      ALLOCATE(valuesB, dbValue, NfieldsB);
+    }
+
     return TRUE;
 }
@@ -116 +134 @@
 int applyConstraintsA(Catalog *catalog, off_t i) {
     off_t m = catalog[0].average[i].measureOffset;
-    off_t n;
     int Nsecfilt = GetPhotcodeNsecfilt();
     
-    dbExtractAveragesInit();
+    // an empty WHERE matches all objects
+    if (!NfieldsA) return TRUE;
+
+    off_t n;
     for (n = 0; n < NfieldsA; n++) {
         valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]);
@@ -130 +150 @@
     int Nsecfilt = GetPhotcodeNsecfilt();
     
-    dbExtractAveragesInit();
+    // an empty WHERE matches all objects
+    if (!NfieldsB) return TRUE;
 
     off_t n;

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

@@ -2 +2 @@
 
 void initialize (int argc, char **argv) {
-
-  int NPHOTCODES;
-  char *codename, *ptr, *list;
-
-  ptr = NULL;
 
   ConfigInit (&argc, argv);
@@ -13 +8 @@
   if (FIT_TARGET == TARGET_MERGE_SOURCE) return;
 
-  /* build a list of accepted photcodes. these will be used by bcatalog to accept or
-     reject loaded data */
+  fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST);
+  fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST);
+  fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
+  fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
 
-  NphotcodesKeep = 0;
-  photcodesKeep = NULL;
-  if (PHOTCODE_KEEP_LIST != NULL) {
-    NPHOTCODES = 10;
-    ALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES);
-
-    /* parse the comma-separated list of photcodesKeep */
-    list = PHOTCODE_KEEP_LIST;
-    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
-      list = NULL; // pass NULL on successive strtok_r calls
-      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_KEEP_LIST);
-      fprintf (stderr, "codename: %s\n", codename);
-      if ((photcodesKeep[NphotcodesKeep] = GetPhotcodebyName (codename)) == NULL) {
-        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
-        exit (1);
-      }
-      NphotcodesKeep ++;
-      CHECK_REALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES, NphotcodesKeep, 10);
-    }
-  }
-
-  NphotcodesSkip = 0;
-  photcodesSkip = NULL;
-  if (PHOTCODE_SKIP_LIST != NULL) {
-    NPHOTCODES = 10;
-    ALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES);
-
-    /* parse the comma-separated list of photcodesSkip */
-    list = PHOTCODE_SKIP_LIST;
-    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
-      list = NULL; // pass NULL on successive strtok_r calls
-      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_SKIP_LIST);
-      fprintf (stderr, "codename: %s\n", codename);
-      if ((photcodesSkip[NphotcodesSkip] = GetPhotcodebyName (codename)) == NULL) {
-        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
-        exit (1);
-      }
-      NphotcodesSkip ++;
-      CHECK_REALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES, NphotcodesSkip, 10);
-    }
-  }
-
-  NphotcodesGroupA = 0;
-  photcodesGroupA = NULL;
-  if (PHOTCODE_A_LIST != NULL) {
-    NPHOTCODES = 10;
-    ALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES);
-
-    /* parse the comma-separated list of photcodesGroupA */
-    fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
-    list = PHOTCODE_A_LIST;
-    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
-      list = NULL; // pass NULL on successive strtok_r calls
-      fprintf (stderr, "codename: %s\n", codename);
-      if ((photcodesGroupA[NphotcodesGroupA] = GetPhotcodebyName (codename)) == NULL) {
-        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
-        exit (1);
-      }
-      NphotcodesGroupA ++;
-      CHECK_REALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES, NphotcodesGroupA, 10);
-    }
-  }
-
-  NphotcodesGroupB = 0;
-  photcodesGroupB = NULL;
-  if (PHOTCODE_B_LIST != NULL) {
-    NPHOTCODES = 10;
-    ALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES);
-
-    fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
-    /* parse the comma-separated list of photcodesGroupB */
-    list = PHOTCODE_B_LIST;
-    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
-      list = NULL; // pass NULL on successive strtok_r calls
-      fprintf (stderr, "codename: %s\n", codename);
-      if ((photcodesGroupB[NphotcodesGroupB] = GetPhotcodebyName (codename)) == NULL) {
-        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
-        exit (1);
-      }
-      NphotcodesGroupB ++;
-      CHECK_REALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES, NphotcodesGroupB, 10);
-    }
-  }
+  photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
+  photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
+  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE);
+  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE);
 
   initstats (STATMODE);
@@ -123 +41 @@
 }
 
+void initialize_client (int argc, char **argv) {
+
+  ConfigInit (&argc, argv);
+  args_client (argc, argv);
+
+  fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST);
+  fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST);
+  fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
+  fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
+
+  photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE);
+  photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE);
+  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE);
+  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE);
+
+  initstats (STATMODE);
+
+  // IMAGE_BAD = ID_IMAGE_ASTROM_POOR | ID_IMAGE_ASTROM_FEW | ID_IMAGE_ASTROM_SKIP;
+  // STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
+  // MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_ASTROM | ID_MEAS_SKIP_ASTROM | ID_MEAS_AREA;
+}
+

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

@@ -1 +1 @@
 # include "relastro.h"
 
-Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
+Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog);
+void bcatalog_show_skips ();
+
+Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) {
 
   int i, j, k, m, Nstar;
   Catalog *catalog, *pcatalog, tcatalog;
 
+  // XXX need to decide how to determine PARALLEL mode...
+  if (PARALLEL && !hostID) {
+    catalog = load_catalogs_parallel (skylist, Ncatalog);
+    return catalog;
+  }
+
   if (VERBOSE) fprintf (stderr, "loading catalog data\n");
 
   ALLOCATE (catalog, Catalog, skylist[0].Nregions);
 
-  initStarMaps();
-
   // load data from each region file, only use bright stars
+  int Ncat = 0;
   for (i = 0; i < skylist[0].Nregions; i++) {
+
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
 
     // we only allow output if we do not use a subset.  in this case,
@@ -19 +30 @@
     pcatalog = subselect ? &tcatalog : &catalog[i];
 
+    // define the catalog file name
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
+    pcatalog->filename = hostID ? hostfile : skylist[0].filename[i];
+
     // set up the basic catalog info
-    pcatalog[0].filename  = skylist[0].filename[i];
-    pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
-    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
-    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
-    pcatalog[0].Nsecfilt  = GetPhotcodeNsecfilt ();
-
+    pcatalog->catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    pcatalog->catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    pcatalog->catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
+    pcatalog->Nsecfilt  = GetPhotcodeNsecfilt ();
+
+    // loads Average, Measure, SecFilt
     if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE2, "w")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
       exit (1);
     }
-    if (VERBOSE2 && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
-
-    //outlier rejection
+    if (!pcatalog[0].Naves_disk) {
+      if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
+      dvo_catalog_unlock (pcatalog);
+      dvo_catalog_free (pcatalog);
+      continue;
+    }
+
+    if (!pcatalog->sorted) {
+      fprintf (stderr, "this database is not sorted.  please sort using addstar -resort\n");
+      exit (1);
+    }
+
+    // outlier rejection (operates on full tables (Average, Measure, Secfilt))
     if (FlagOutlier) {
       FlagOutliers(&tcatalog);
     }
 
-    updateStarMaps (&tcatalog);
-
     // select only the brighter stars
     if (subselect) {
-      bcatalog (&catalog[i], &tcatalog);
+      // results are in Average, Secfilt, and MeasureTiny
+      // Ncat tracks the actually used catalogs
+      bcatalog (&catalog[Ncat], &tcatalog);
       dvo_catalog_unlock (&tcatalog);
       dvo_catalog_free (&tcatalog);
     } else {
       if (RESET) {
-        for (j = 0; j < catalog[i].Naverage; j++) {
-          catalog[i].average[j].flags = 0;
-          m = catalog[i].average[j].measureOffset;
-          for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
-            catalog[i].measure[m+k].dbFlags = 0;
+        for (j = 0; j < catalog[Ncat].Naverage; j++) {
+          catalog[Ncat].average[j].flags = 0;
+          m = catalog[Ncat].average[j].measureOffset;
+          for (k = 0; k < catalog[Ncat].average[j].Nmeasure; k++) {
+            catalog[Ncat].measure[m+k].dbFlags = 0;
           }
         }
       }
     }
-  }
-
-  // XXX keep this test?
+    Ncat ++;
+  }
+
+  // XXX TEST : bcatalog_show_skips();
+
   Nstar = 0;
-  for (i = 0; i < skylist[0].Nregions; i++) {
+  for (i = 0; i < Ncat; i++) {
     Nstar += catalog[i].Naverage;
+    if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) {
+      fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n");
+      fprintf (stderr, "  If you need to do this, and you can afford the RAM, rebuild with IDX_T set to off_t (see relastro.h, ImagesOps.c)\n");
+      exit (3);
+    }
   }
   if (Nstar < 2) { 
@@ -66 +99 @@
   }
 
-  // XXX consider only returning the populated catalogs
-  *Ncatalog = skylist[0].Nregions;
+  // only return the populated catalogs
+  REALLOCATE (catalog, Catalog, Ncat);
+  *Ncatalog = Ncat;
   return (catalog);
 }
@@ -76 +110 @@
    need to use an XCLD lock here.  
 */
+
+// CATDIR is supplied globally
+# define DEBUG 1
+Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
+
+  // load the list of hosts
+  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
+  if (!table) {
+    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
+    exit (1);
+  }    
+
+  int i;
+  for (i = 0; i < table->Nhosts; i++) {
+
+    // ensure that the paths are absolute path names
+    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
+    free (table->hosts[i].pathname);
+    table->hosts[i].pathname = tmppath;
+
+    ALLOCATE (table->hosts[i].results, char, 1024);
+    snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname);
+
+    // options / arguments that can affect relastro_client -load:
+    // VERBOSE, VERBOSE2
+    // RESET (-reset)
+    // TimeSelect -time
+    // (note that psfQual is applied rigidly at 0.85, as is the galaxy test)
+    // ImagSelect, ImagMin, ImagMax
+    // MaxDensityUse, MaxDensityValue
+
+    char command[1024];
+    snprintf (command, 1024, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f -D RELASTRO_SIGMA_LIM %f", 
+              table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM);
+
+    char tmpline[1024];
+    if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);           strcpy (command, tmpline); }
+
+    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
+    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
+    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
+    if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
+    if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
+    
+    if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
+    
+    if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
+    if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
+    if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                     strcpy (command, tmpline); }
+    if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);     strcpy (command, tmpline); }
+    if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);   strcpy (command, tmpline); }
+    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+
+    if (TimeSelect) { 
+      char *tstart = ohana_sec_to_date (TSTART);
+      char *tstop  = ohana_sec_to_date (TSTOP);
+      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); 
+      free (tstart);
+      free (tstop);
+      strcpy (command, tmpline); 
+    }
+
+    fprintf (stderr, "command: %s\n", command);
+
+    if (PARALLEL_MANUAL) continue;
+
+    if (PARALLEL_SERIAL) {
+      int status = system (command);
+      if (status) {
+        fprintf (stderr, "ERROR running relastro_client\n");
+        exit (2);
+      }
+    } else {
+      // launch the job on the remote machine (no handshake)
+      int errorInfo = 0;
+      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
+      if (!pid) {
+        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
+        exit (1);
+      }
+      table->hosts[i].pid = pid; // save for future reference
+    }
+  }
+
+  if (PARALLEL_MANUAL) {
+    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
+    getchar();
+  }
+  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
+    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
+  }
+
+  // each host generates a BrightCatalog structure, with the measure, average, etc value
+  // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt).  I need to
+  // split out the per-catalog measurements into separate catalog entries.
+
+  // set up an initial array of catalogs
+  int Nsecfilt  = GetPhotcodeNsecfilt (); 
+  CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt);
+
+  for (i = 0; i < table->Nhosts; i++) {
+
+    BrightCatalog *bcatalog = BrightCatalogLoad (table->hosts[i].results);
+    assert (bcatalog);
+    
+    BrightCatalogSplit (catalogs, bcatalog);
+
+    free (bcatalog->average);
+    free (bcatalog->measure);
+    free (bcatalog->secfilt);
+    free (bcatalog);
+  }
+
+  Catalog *catalog = catalogs->catalog;
+  *Ncatalog = catalogs->Ncatalog;
+  BrightCatalogSplitFree (catalogs);
+
+  int Nmeasure = 0;
+  int Naverage = 0;
+  for (i = 0; i < catalogs->Ncatalog; i++) {
+    Nmeasure += catalogs->catalog[i].Nmeasure;
+    Naverage += catalogs->catalog[i].Naverage;
+  }
+
+  fprintf (stderr, "loaded %d catalogs, using a total of %d stars (%d measures)\n", catalogs->Ncatalog, Naverage, Nmeasure);
+
+  return (catalog);
+}      

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

@@ -7 +7 @@
   fprintf (stderr, MSG, __VA_ARGS__); }
 
-SkyList *load_images (FITS_DB *db, SkyRegion *region) {
+int load_images (FITS_DB *db, SkyList *skylist) {
 
   Image     *image, *subset;
@@ -14 +14 @@
   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, TRUE, SKY_DEPTH, VERBOSE);
-  SkyTableSetFilenames (sky, CATDIR, "cpt");
-  
-  // determine the populated SkyRegions overlapping the requested area
-  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
@@ -37 +22 @@
       exit (2);
   }
-
   MARKTIME("  convert image table: %f sec\n", dtime);
 
@@ -53 +37 @@
   if (!UPDATE) dvo_image_unlock (db); 
 
-  return (skylist);
+  return TRUE;
 }
 
@@ -71 +55 @@
     memcpy (vtable[0].buffer[i], &image[i], Nx);
   }
-  return (TRUE);
+  return TRUE;
 }

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

@@ -3 +3 @@
 int main (int argc, char **argv) {
 
-  int i, status, Ncatalog;
-  Catalog *catalog;
-  FITS_DB db;
-  struct timeval start, stop;
-  SkyList *skylist = NULL;
-
-  gettimeofday (&start, (void *) NULL);
-
   /* get configuration info, args */
   initialize (argc, argv);
 
-  /* the object analysis is a separate process iterating over catalogs */
-  if (FIT_TARGET == TARGET_OBJECTS) {
-    relastro_objects ();
-    exit (0);
-  }
+  SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch);
 
-  /* the object analysis is a separate process iterating over catalogs */
-  if (FIT_TARGET == TARGET_HIGH_SPEED) {
-    high_speed_catalogs ();
-    exit (0);
-  }
+  switch (FIT_TARGET) {
+    case TARGET_UPDATE_OBJECTS:
+      /* the object analysis is a separate process iterating over catalogs */
+      relastro_objects (skylist, 0, NULL);
+      exit (0);
 
-  /* the object analysis is a separate process iterating over catalogs */
-  if (FIT_TARGET == TARGET_MERGE_SOURCE) {
-    relastro_merge_source ();
-    exit (0);
-  }
+    case TARGET_HIGH_SPEED:
+      /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */
+      high_speed_catalogs (sky, skylist, 0, NULL);
+      exit (0);
 
-  /* register database handle with shutdown procedure */
-  set_db (&db);
+    case TARGET_MERGE_SOURCE:
+      /* a special method to manually merge unlinked detections of sources togther (not parallel) */
+      relastro_merge_source (sky);
+      exit (0);
 
-  /* lock and load the image db table */
-  status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
-  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
-  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);
+    case TARGET_SIMPLE:
+    case TARGET_CHIPS:
+    case TARGET_MOSAICS:
+      relastro_images (skylist);
+      exit (0);
 
-  /* 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 */
-  // XXX pass in the image table
-  // XXX who carries the image grid?
-  catalog = load_catalogs (skylist, &Ncatalog, TRUE);
-  MARKTIME("load catalog data: %f sec\n", dtime);
-
-  /* match measurements with images */
-  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) {
-    // plot_star_coords (catalog, Ncatalog);
-    // plot_mosaic_fields (catalog);
-  }
-
-  // set test points based on the starmap
-  createStarMapPoints();
-
-  /* major modes */
-  switch (FIT_TARGET) {
-    case TARGET_SIMPLE:
-      for (i = 0; i < NLOOP; i++) {
-        UpdateObjects (catalog, Ncatalog);
-        UpdateSimple (catalog, Ncatalog);
-      }
-      break;
-
-    case TARGET_CHIPS:
-      for (i = 0; i < NLOOP; i++) {
-        UpdateObjects (catalog, Ncatalog);
-        UpdateChips (catalog, Ncatalog);
-        MARKTIME("update chips: %f sec\n", dtime);
-      }
-      // create summary plots of the process
-      // relastroVisualSummaryChips();
-      break;
-
-    case TARGET_MOSAICS:
-      for (i = 0; i < NLOOP; i++) {
-        UpdateObjects (catalog, Ncatalog);
-        UpdateMosaic (catalog, Ncatalog);
-      }
-      break;
+    case TARGET_UPDATE_OFFSETS:
+      // iterate over catalogs to make detection coordinates consistant
+      UpdateObjectOffsets (skylist, 0, NULL);
 
     default:
-      fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__);
-      exit (2);
+      fprintf (stderr, "impossible!\n");
+      abort();
   }
-
-  if (!UPDATE) exit (0);
-
-  // free the image / measurement pointers
-  freeImageBins (1);
-
-  reload_images (&db);
-    
-  // iterate over catalogs to make detection coordinates consistant
-  UpdateObjectOffsets (skylist);
-
-  // iterate over catalogs to make detection coordinates consistant
-  // FixProblemImages (skylist);
-
-  // save the updated image parameters
-  dvo_image_update (&db, VERBOSE);
-  dvo_image_unlock (&db); 
-
-  exit (0);
+  exit (1);
 }

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

@@ -1 +1 @@
 # include "relastro.h"
 
-int relastro_merge_source () {
+int relastro_merge_source (SkyTable *sky) {
 
   int i, m;
-  SkyTable *sky = NULL;
   Catalog catalog_src, catalog_dst;
 
@@ -13 +12 @@
   dvo_set_catdir(CATDIR);
 
-  // load the current sky table (layout of all SkyRegions) 
-  sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, SKY_DEPTH, VERBOSE);
-  SkyTableSetFilenames (sky, CATDIR, "cpt");
-  
   // we are merging the detections of the src into the dst : OBJ_ID_SRC,CAT_ID_SRC -> OBJ_ID_DST,CAT_ID_DST
 
@@ -52 +47 @@
   catalog_src.Nsecfilt  = GetPhotcodeNsecfilt ();
 
-  if (!dvo_catalog_open (&catalog_src, region_src, VERBOSE, "w")) {
+  if (!dvo_catalog_open (&catalog_src, region_src, VERBOSE2, "w")) {
     fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog_src.filename);
     exit (1);
@@ -127 +122 @@
     resort_catalog (&catalog_src);
 
-    dvo_catalog_save (&catalog_src, VERBOSE);
+    dvo_catalog_save (&catalog_src, VERBOSE2);
     dvo_catalog_unlock (&catalog_src);
     dvo_catalog_free (&catalog_src);

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

@@ -1 +1 @@
 # include "relastro.h"
 
-int relastro_objects () {
+int relastro_objects_parallel (SkyList *sky);
+
+int relastro_objects (SkyList *skylist, int hostID, char *hostpath) {
 
   int i, j, k, m;
 
-  SkyTable *sky = NULL;
-  SkyList *skylist = NULL;
   Catalog catalog;
 
+  DVOMeasureFlags measureBits = 
+    ID_MEAS_POOR_ASTROM |
+    ID_MEAS_SKIP_ASTROM | 
+    ID_MEAS_USED_OBJ    |
+    ID_MEAS_USED_CHIP;
 
-  // load the current sky table (layout of all SkyRegions) 
-  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
-  SkyTableSetFilenames (sky, CATDIR, "cpt");
-  
-  // determine the populated SkyRegions overlapping the requested area (default depth)
-  if (UserCatalog) {
-    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
-  } else {
-    skylist = SkyListByPatch (sky, -1, &UserPatch);
+  DVOAverageFlags averageBits = 
+    ID_STAR_FEW     |
+    ID_STAR_POOR    |
+    ID_STAR_FIT_AVE |
+    ID_STAR_FIT_PM  |
+    ID_STAR_FIT_PAR |
+    ID_STAR_USE_AVE |
+    ID_STAR_USE_PM  |
+    ID_STAR_USE_PAR; 
+
+  // XXX need to decide how to determine PARALLEL mode...
+  if (PARALLEL && !hostID) {
+    relastro_objects_parallel (skylist);
+    return TRUE;
   }
 
@@ -24 +34 @@
   for (i = 0; i < skylist[0].Nregions; i++) {
 
+    // does this host ID match the desired location for the table?
+    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
+
     // set up the basic catalog info
-    catalog.filename  = skylist[0].filename[i];
+    char hostfile[1024];
+    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
+    catalog.filename  = hostID ? hostfile : skylist[0].filename[i];
     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
@@ -31 +46 @@
     catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
 
-    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+    // loads Average, Measure, SecFilt
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "w")) {
       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
       exit (1);
     }
     if (!catalog.Naves_disk) {
-      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
       dvo_catalog_unlock (&catalog);
       dvo_catalog_free (&catalog);
@@ -42 +58 @@
     }
 
+    // operates on full tables (Average, Measure, Secfilt)
     if (FlagOutlier) {
       FlagOutliers(&catalog);
     }
 
-    // XXX consider what gets reset (only ASTROM flags)
+    // reset only the astrometry-related average and measure bits
     if (RESET) {
       for (j = 0; j < catalog.Naverage; j++) {
-        catalog.average[j].flags = 0;
+        catalog.average[j].flags &= ~averageBits;
         m = catalog.average[j].measureOffset;
         for (k = 0; k < catalog.average[j].Nmeasure; k++) {
-          catalog.measure[m+k].dbFlags = 0;
+          catalog.measure[m+k].dbFlags &= ~measureBits;
         }
       }
     }
 
+    populate_tiny_values(&catalog, DVO_TV_MEASURE);
+
     UpdateObjects (&catalog, 1);
+
+    free_tiny_values(&catalog);
 
     if (!UPDATE) {
@@ -70 +91 @@
   return (TRUE);
 }
+
+// CATDIR is supplied globally
+# define DEBUG 1
+int relastro_objects_parallel (SkyList *sky) {
+
+  // launch the setphot_client jobs to the parallel hosts
+
+  // load the list of hosts
+  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
+  if (!table) {
+    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
+    exit (1);
+  }    
+
+  int i;
+  for (i = 0; i < table->Nhosts; i++) {
+
+    // ensure that the paths are absolute path names
+    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
+    free (table->hosts[i].pathname);
+    table->hosts[i].pathname = tmppath;
+
+    // options / arguments that can affect relastro_client -update-objects:
+    // VERBOSE, VERBOSE2
+    // RESET (-reset)
+    // TimeSelect -time
+    // (note that psfQual is applied rigidly at 0.85, as is the galaxy test)
+    // ImagSelect, ImagMin, ImagMax
+    // MaxDensityUse, MaxDensityValue
+
+    // FIT_MODE
+    // PM_TOOFEW
+    // SRC_MEAS_TOOFEW
+    // FIT_TARGET
+
+    char command[1024];
+    snprintf (command, 1024, "relastro_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", 
+              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+
+    char tmpline[1024];
+    if (FIT_MODE == FIT_PM_ONLY)         { snprintf (tmpline, 1024, "%s -pm",    command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PAR_ONLY)        { snprintf (tmpline, 1024, "%s -par",   command);           strcpy (command, tmpline); }
+    if (FIT_MODE == FIT_PM_AND_PAR)      { snprintf (tmpline, 1024, "%s -pmpar", command);           strcpy (command, tmpline); }
+
+    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
+    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
+    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
+    if (UPDATE)        { snprintf (tmpline, 1024, "%s -update",         command);                    strcpy (command, tmpline); }
+
+    if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
+    if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
+    
+    if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
+    
+    if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); }
+    if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); }
+    if (PhotFlagSelect)    { snprintf (tmpline, 1024, "%s +photflags",   command);                     strcpy (command, tmpline); }
+    if (PhotFlagBad)       { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad);     strcpy (command, tmpline); }
+    if (PhotFlagPoor)      { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor);   strcpy (command, tmpline); }
+    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
+
+    if (TimeSelect) { 
+      char *tstart = ohana_sec_to_date (TSTART);
+      char *tstop  = ohana_sec_to_date (TSTOP);
+      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); 
+      free (tstart);
+      free (tstop);
+      strcpy (command, tmpline); 
+    }
+
+    fprintf (stderr, "command: %s\n", command);
+
+    if (PARALLEL_MANUAL) continue;
+
+    if (PARALLEL_SERIAL) {
+      int status = system (command);
+      if (status) {
+        fprintf (stderr, "ERROR running relastro_client\n");
+        exit (2);
+      }
+    } else {
+      // launch the job on the remote machine (no handshake)
+      int errorInfo = 0;
+      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
+      if (!pid) {
+        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
+        exit (1);
+      }
+      table->hosts[i].pid = pid; // save for future reference
+    }
+  }
+
+  if (PARALLEL_MANUAL) {
+    fprintf (stderr, "run the relastro_client commands above.  when these are done, hit return\n");
+    getchar();
+  }
+  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
+    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
+  }
+
+  return TRUE;
+}      

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

@@ -1 +1 @@
 # include "relastro.h"
-
-# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); }
-# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } }
 
 // sort the measure Sequence based on the average Sequence entries

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

@@ -8 +8 @@
   for (i = 0; i < Ncatalog; i++) {
 
-    if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename);
-    dvo_catalog_save (&catalog[i], VERBOSE); 
+    if (VERBOSE2) fprintf (stderr, "saving catalog %s\n", catalog[i].filename);
+    dvo_catalog_save (&catalog[i], VERBOSE2); 
     dvo_catalog_unlock (&catalog[i]);
     dvo_catalog_free (&catalog[i]);

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

@@ -25 +25 @@
   
   Image *image;
-  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE, D_NIMAGE;
   off_t *line_number;
   int InRange, found;
@@ -68 +68 @@
   DminSkyRegion = +90.0;
   DmaxSkyRegion = -90.0;
+
+  D_NIMAGE = 6000;
 
   /* compare with each region file */
@@ -113 +115 @@
 
   nimage = 0;
-  NIMAGE = 100;
+  NIMAGE = D_NIMAGE;
   ALLOCATE (image, Image, NIMAGE);
   ALLOCATE (line_number, off_t, NIMAGE);
@@ -167 +169 @@
 
     /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */
-    double RminImage = 360.0;
-    double RmaxImage =   0.0;
+    double RminImage = RmidSkyRegion + 180.0;
+    double RmaxImage = RmidSkyRegion - 180.0;
     double DminImage = +90.0;
     double DmaxImage = -90.0;
@@ -194 +196 @@
     if (RmaxImage < RminSkyRegion) continue;
 
+    if (!strncmp(timage[i].name, "o5903g0638o", 10) && (timage[i].photcode == 10355)) {
+      fprintf (stderr, "test image\n");
+    }
+
+    // require that the full image be inside the region of interest (only for non PHU
+    // images) XXX : if we calibrate the mosaic, require that the full mosaic be inside
+    // the region and all of its chips..
+# define FULL_OVERLAP 1
+    if (FULL_OVERLAP && strcmp(&timage[i].coords.ctype[4], "-DIS")) {
+      if (RmaxImage > UserPatch.Rmax) continue;
+      if (RminImage < UserPatch.Rmin) continue;
+      if (DmaxImage > UserPatch.Dmax) continue;
+      if (DminImage < UserPatch.Dmin) continue;
+    }
+
     // image overlaps region, keep it
     if (USE_BASIC_CHECK) goto found_it;
@@ -243 +260 @@
     }
     if (RESET) {
-      // XXX do we need / want to do this in relastro?
-      // assignMcal (&image[nimage], (double *) NULL, -1);
-      // image[nimage].Mcal = NAN;
-      // image[nimage].dMcal = NAN;
+      // this only resets the astrometry image flags, not the photometry ones
       image[nimage].flags &= ~badImage;
     }
@@ -252 +266 @@
     nimage ++;
     if (nimage == NIMAGE) {
-      NIMAGE += 100;
+      NIMAGE += D_NIMAGE;
+      D_NIMAGE = MAX (100000, D_NIMAGE * 1.5);
       REALLOCATE (image, Image, NIMAGE);
       REALLOCATE (line_number, off_t, NIMAGE);
@@ -264 +279 @@
   REALLOCATE (line_number, off_t, MAX (nimage, 1));
   free (skycoords);
+  free (RmaxSky);
+  free (index);
 
   *Nimage  = nimage;

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

@@ -1 +1 @@
 # include "relastro.h"
+# define NPTS 64
+# define J2000 2451545.       /* Julian date at standard epoch */
+# define T1970 2440587.500000 /* JD at UNIX ref time */
 
 int main (int argc, char **argv) {
+   
+  int i, Npts;
+  char line[1024];
+  double R[NPTS], D[NPTS], Time[NPTS], JD[NPTS];
+  double X[NPTS], Y[NPTS], dX[NPTS], dY[NPTS], pX[NPTS], pY[NPTS];
+  double Tref[NPTS], Tjyrs[NPTS], TrefS, TjyrsS, TrefMean, TjyrsMean, Ro, Do;
+  Coords coords;
+  PMFit fitPAR;
 
-    if (argc != N) {
-        fprintf (stderr, "USAGE: testparallax (file.dat)\n");
-        exit (1);
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: testparallax (file.dat)\n");
+    exit (1);
+  }
+
+  // test parallax program
+  FILE *f = fopen (argv[1], "r");
+  if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); }
+
+  // XXX uncomment to skip first line scan_line(f, line);
+    
+  Npts = 0;
+  for (i = 0; scan_line(f, line) != EOF; i++) {
+    if (Npts == NPTS) {
+      fprintf (stderr, "too many point: use dynamic alloc\n");
+      exit (2);
     }
+    if (line[0] == '#') continue;
+    dparse (&R[Npts],    1, line);
+    dparse (&D[Npts],    2, line);
+    dparse (&JD[Npts],   3, line);
+    dparse (&Time[Npts], 3, line);
 
-    // test parallax program
-    FILE *f = fopen (argv[1], "r");
-    if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); }
+    dparse (&dX[Npts],   4, line);
+    dparse (&dY[Npts],   4, line);
 
-    char buffer[1024];
-    scan_line(f, line);
-    
-    while (scan_line(f, line) != EOF) {
-        dparse (&jd,   1, line);
-        dparse (&time, 2, line);
-        dparse (&ra,   3, line);
-        dparse (&dec,  4, line);
+    // trent's file
+    // dparse (&JD[Npts],   1, line);
+    // dparse (&Time[Npts], 2, line);
+    // dparse (&R[Npts],    3, line);
+    // dparse (&D[Npts],    4, line);
+    Npts ++;
+  }
 
-    }
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  coords.crval1 = 0;
+  coords.crval2 = 0;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "RA---SIN");
 
-    // project to local coords
+  // use one point as a local reference
+  coords.crval1 = R[0];
+  coords.crval2 = D[0];
 
-    // run fitter
-    FitPMandPar
+  TrefS = TjyrsS = 0.0;
 
+  // project to local coords
+  for (i = 0; i < Npts; i++) {
+    RD_to_XY (&X[i], &Y[i], R[i], D[i], &coords);
 
+    Tjyrs[i] = Time[i] / 365.25;
+    Tref[i] = Time[i] / 365.25;
+    // Tunix[i] = (JD[i] - T1970) / 365.25; // time relative to T1970 in years
+    // Tjyrs[i] = (JD[i] - J2000) / 365.25; // time relative to J2000 in years
+
+    TrefS   += Tref[i];
+    TjyrsS   += Tjyrs[i];
+
+    // nominal for PS1
+    // dX[i] = 0.020;
+    // dY[i] = 0.020;
+  }
+  TrefMean = TrefS / Npts;
+  TjyrsMean = TjyrsS / Npts;
+
+  for (i = 0; i < Npts; i++) {
+    Tref[i] -= TrefMean;
+    ParFactor (&pX[i], &pY[i], R[i], D[i], Tjyrs[i], 0.0);
+    fprintf (stderr, "%f %f : %f %f : %f %f : %f %f\n", R[i], D[i], X[i], Y[i], Tref[i], Tjyrs[i], pX[i], pY[i]);
+  }
+
+  // run fitter
+  // FitPMandPar (&fitPAR, X, dX, Y, dY, Tref, pX, pY, Npts, TRUE);
+  FitPMandPar (&fitPAR, X, dX, Y, dY, Tjyrs, pX, pY, Npts, TRUE);
+  XY_to_RD (&Ro, &Do, fitPAR.Ro, fitPAR.Do, &coords);
+
+  fprintf (stderr, "Rx   : %f\n", fitPAR.Ro);
+  fprintf (stderr, "Dx   : %f\n", fitPAR.Do);
+  fprintf (stderr, "Ro   : %f\n", Ro);
+  fprintf (stderr, "Do   : %f\n", Do);
+  fprintf (stderr, "dRo  : %f\n", fitPAR.dRo);
+  fprintf (stderr, "dDo  : %f\n", fitPAR.dDo);
+  fprintf (stderr, "uR   : %f\n", fitPAR.uR);
+  fprintf (stderr, "uD   : %f\n", fitPAR.uD);
+  fprintf (stderr, "duR  : %f\n", fitPAR.duR);
+  fprintf (stderr, "duD  : %f\n", fitPAR.duD);
+  fprintf (stderr, "p    : %f\n", fitPAR.p);
+  fprintf (stderr, "dp   : %f\n", fitPAR.dp);
+  fprintf (stderr, "uTot : %f\n", hypot(fitPAR.uR,fitPAR.uD));
+  fprintf (stderr, "PA   : %f\n", DEG_RAD*atan2(fitPAR.uR,fitPAR.uD) + 360.0);
+  
+  exit (0);
 }
+