Changeset 36833 for trunk/Ohana/src/relastro

Timestamp:

Jun 7, 2014, 6:33:11 AM (12 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20140423: add dcr-color option, ref_color, bad detection clipping to relastro; fix RA,DEC for mosaics in imextract; allow photdbc to restrict images to output region; change PHOTCODE_SKIP_LIST to PHOTCODE_KEEP_LIST in photdbc; add -reset-to-photcode option to relastro; ensure relastro child processes get the complete set of options; allow relastro region hosts to continue after load-object clients are done, but before the data is loaded; invalidate extreme proper motion; limit relastro catalog density on a grid; require more detections for higher order fits; option to delete duplicate images based on exp time/date and photcode; have kapa use XSync instead of XFlush to avoid flashing; handle some cases to avoid crashing kapa; extend kapa connection timeout; allow kapa to handle > 2GB images; ensure dvomerge works even if incoming db has no images; add dvo fuctions coordmosaic & psastro_model

Location:

trunk/Ohana/src/relastro

Files:

• Makefile (modified) (1 diff)
• include/relastro.h (modified) (7 diffs)
• src/BrightCatalog.c (modified) (8 diffs)
• src/FitChip.c (modified) (7 diffs)
• src/GetAstromError.c (modified) (2 diffs)
• src/ImageOps.c (modified) (12 diffs)
• src/UpdateObjectOffsets.c (modified) (1 diff)
• src/UpdateObjects.c (modified) (14 diffs)
• src/args.c (modified) (6 diffs)
• src/bcatalog.c (modified) (10 diffs)
• src/dvo_astrom_ops.c (modified) (2 diffs)
• src/initialize.c (modified) (2 diffs)
• src/launch_region_hosts.c (modified) (3 diffs)
• src/load_catalogs.c (modified) (8 diffs)
• src/relastro_client.c (modified) (1 diff)
• src/relastro_images.c (modified) (2 diffs)
• src/relastro_objects.c (modified) (1 diff)
• src/relastro_parallel_images.c (modified) (1 diff)

trunk/Ohana/src/relastro/Makefile

@@ -127 +127 @@
 $(SRC)/plotstuff.$(ARCH).o           \
 $(SRC)/relastroVisual.$(ARCH).o \
+$(SRC)/syncfile.$(ARCH).o \
 $(SRC)/BrightCatalog.$(ARCH).o 
 

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

@@ -83 +83 @@
   double L, M;  /* Focal Plane   - pixels  */
   double X, Y;  /* Chip Coords   - pixels  */
-  double Mag;
-  double dMag;
-  double dPos;
+  float Mag;
+  float Color;
+  float dMag;
+  float dPos;
   int mask;
   int Nmeas;
@@ -207 +208 @@
 double MaxDensityValue;
 
-char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
-int           NphotcodesKeep,      NphotcodesSkip;
-PhotCode     **photcodesKeep,     **photcodesSkip;
+char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST, *PHOTCODE_RESET_LIST;
+int           NphotcodesKeep,      NphotcodesSkip,      NphotcodesReset;
+PhotCode     **photcodesKeep,     **photcodesSkip,     **photcodesReset;
 
 char          *PHOTCODE_A_LIST,  *PHOTCODE_B_LIST;
@@ -217 +218 @@
 SkyRegionSelection SELECTION;
 
+char         *DCR_COLOR_POS,    *DCR_COLOR_NEG;
+PhotCode     *DCR_PHOTCODE_POS, *DCR_PHOTCODE_NEG; 
+int           DCR_NSEC_POS,      DCR_NSEC_NEG; 
+
 int ImagSelect;
 double ImagMin, ImagMax;
@@ -233 +238 @@
 int    CLIP_THRESH;
 int USE_BASIC_CHECK;
+
+int ExcludeBogus;
+double ExcludeBogusRadius;
 
 FitMode FIT_MODE;
@@ -307 +315 @@
 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
 int           liststats_pos       PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB));
-Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile));
 int           load_images         PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap));
 Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap));
@@ -417 +425 @@
 int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
 int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt);
+double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanD_Big (Measure *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);
@@ -521 +531 @@
 int ImageTableSave (char *filename, Image *images, off_t Nimages);
 int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage);
+
+float getColor (off_t meas, int cat);
+
+int strextend (char *input, char *format,...);
+
+int areImagesLoaded ();
+int areImagesMatched ();

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

@@ -167 +167 @@
     GET_COLUMN(measureOffset,  "OFF_MEASURE", int);
     GET_COLUMN(missingOffset,  "OFF_MISSING", int);
-    GET_COLUMN(extendOffset,   "OFF_EXTEND",  int);
+    GET_COLUMN(refColor,   "REF_COLOR",  float);
     GET_COLUMN(flags,          "FLAGS",       int);
     GET_COLUMN(photFlagsUpper, "PHOTFLAGS_U", int);
@@ -202 +202 @@
       average[i].measureOffset   = measureOffset[i]   ; 
       average[i].missingOffset   = missingOffset[i]   ; 
-      average[i].extendOffset    = extendOffset[i]    ;  
+      average[i].refColor    = refColor[i]    ;  
       average[i].flags           = flags[i]           ;
       average[i].photFlagsUpper  = photFlagsUpper[i]  ;
@@ -234 +234 @@
     free (measureOffset);
     free (missingOffset);
-    free (extendOffset);
+    free (refColor);
     free (flags);
     free (photFlagsUpper);
@@ -476 +476 @@
     gfits_define_bintable_column (&theader, "J", "OFF_MEASURE", "offset to first psf measurement",                        "", 1.0, 1.0*0x8000);
     gfits_define_bintable_column (&theader, "J", "OFF_MISSING", "offset to first missing obs",                            "", 1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "OFF_EXTEND",  "offset to first extended measurement",                   "", 1.0, 1.0*0x8000);
+    gfits_define_bintable_column (&theader, "E", "REF_COLOR",   "reference color",                                        "", 1.0, 0.0);
     gfits_define_bintable_column (&theader, "J", "FLAGS",       "average object flags (star; ghost; etc)",                "", 1.0, 1.0*0x8000);
     gfits_define_bintable_column (&theader, "J", "PHOTFLAGS_U", "upper bit of 2 bit summary of per-measure photflags",    "", 1.0, 1.0*0x8000);
@@ -510 +510 @@
     int      *measureOffset ; ALLOCATE (measureOffset , int     , catalog->Naverage);
     int      *missingOffset ; ALLOCATE (missingOffset , int     , catalog->Naverage);
-    int      *extendOffset  ; ALLOCATE (extendOffset  , int     , catalog->Naverage);
+    float    *refColor      ; ALLOCATE (refColor  , float     , catalog->Naverage);
     int      *flags         ; ALLOCATE (flags         , int     , catalog->Naverage);
     int      *photFlagsUpper; ALLOCATE (photFlagsUpper, int     , catalog->Naverage);
@@ -543 +543 @@
       measureOffset[i]   = average[i].measureOffset   ; 
       missingOffset[i]   = average[i].missingOffset   ; 
-      extendOffset[i]    = average[i].extendOffset    ;  
+      refColor[i]    = average[i].refColor    ;  
       flags[i]           = average[i].flags           ;
       photFlagsUpper[i]  = average[i].photFlagsUpper  ;
@@ -575 +575 @@
     gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset,   catalog->Naverage);
     gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset,   catalog->Naverage);
-    gfits_set_bintable_column (&theader, &ftable, "OFF_EXTEND",  extendOffset,    catalog->Naverage);
+    gfits_set_bintable_column (&theader, &ftable, "REF_COLOR",   refColor,        catalog->Naverage);
     gfits_set_bintable_column (&theader, &ftable, "FLAGS",       flags,           catalog->Naverage);
     gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_U", photFlagsUpper,  catalog->Naverage);
@@ -605 +605 @@
     free (measureOffset);
     free (missingOffset);
-    free (extendOffset);
+    free (refColor);
     free (flags);
     free (photFlagsUpper);

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

@@ -4 +4 @@
 int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
 
-  int i, NstatFull, Nstat, Niter;
+  int i, NstatFull, Nstat, Niter, Nkeep;
   float dLsig, dMsig, dRsig;
   float dLsigFull, dMsigFull, dRsigFull;
@@ -24 +24 @@
 
     // measure the scatter distribution (use only the bright end detections)
+    Nkeep = 0;
     for (i = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
@@ -42 +43 @@
         continue;
       }
-    }
+      Nkeep ++;
+    }
+
+    // I'm rejecting some points from the fit above; I count the remainders and actually
+    // use that count to set order_use below
 
     // figures to assess the fitting process:
@@ -50 +55 @@
     relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
 
-    // fit the requested order polynomial
-    if (CHIPORDER > 0) {
-      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;
-    }
+    // set the maximum order for the polynomial (based on number of stars kept above)
+    int order_use = 0;
+    if (Nkeep >  5) order_use = 1; // 4 stars per polynomial term (per dimension)
+    if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension)
+    if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension)
+    if (order_use < 1) {
+      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep);
+      image[0].flags |= ID_IMAGE_ASTROM_FEW;
+      return FALSE;
+    }
+    image[0].coords.Npolyterms = order_use;
 
     if (fit) fit_free (fit);
@@ -71 +75 @@
       fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
     }
-
-# if (0)
-    // check if the fit has enough data points for the polynomial order
-    skip = FALSE;
-    switch (image[0].coords.Npolyterms) {
-      case 0:
-      case 1:
-        skip = (fit[0].Npts < 8);
-        break;
-      case 2:
-        skip = (fit[0].Npts < 11);
-        break;
-      case 3:
-        skip = (fit[0].Npts < 15);
-        break;
-      default:
-        fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms);
-        skip = TRUE;
-    }
-    if (skip) {
-      if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
-      fit_free (fit);
-      image[0].flags |= ID_IMAGE_ASTROM_FEW;
-      return FALSE;
-    }
-# endif
 
     // measure the fit, update the coords & object coordinates
@@ -134 +112 @@
   }
 
+  int Ncolor = 0;
+  float *colorList = NULL;
+  ALLOCATE (colorList, float, Nmatch);
+
+  for (i = 0; i < Nmatch; i++) {
+    if (!raw[i].mask) continue;
+    if (!isnan(ref[i].Color)) {
+      colorList[Ncolor] = ref[i].Color;
+      Ncolor ++;
+    }
+  }
+  fsort (colorList, Ncolor);
+  float colorMedian = (Ncolor > 0) ? colorList[(int)(0.5*Ncolor)] : NAN;
+  image[0].refColor = colorMedian;
+  free (colorList);
+
   GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
   GetScatterRawRef(&dLsig,     &dMsig,     &dRsig,     &Nstat,     raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
@@ -151 +145 @@
   image[0].dYpixSys = dMsig;
   image[0].nFitAstrom = fit[0].Npts;
+
+  // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom);
 
   if (fit) fit_free (fit);

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

@@ -30 +30 @@
 
   code  = GetPhotcodebyCode (measure[0].photcode);
+  if (!code) return NAN;
 
   // do not raise an exception, just send back the result
@@ -72 +73 @@
 
   code  = GetPhotcodebyCode (measure[0].photcode);
+  if (!code) return NAN;
 
   // do not raise an exception, just send back the result

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

@@ -39 +39 @@
 // N_onImage was 'Nlist'
 // N_ONIMAGE was 'NLIST'
+
+int areImagesLoaded () {
+
+  if (image) return TRUE;
+  return FALSE;
+}
+
+int areImagesMatched () {
+
+  if (MeasureToImage) return TRUE;
+  return FALSE;
+}
 
 Image *getimages (off_t *N, off_t **line_number) {
@@ -333 +345 @@
 */
 
+float getColor (off_t meas, int cat) {
+
+  off_t i;
+
+  i = MeasureToImage[cat][meas];
+  if (i == -1) return (NAN);
+  return (image[i].refColor);
+}
+
 void plot_images () {
 
@@ -468 +489 @@
     }
 
+    X = measure[0].Xccd;
+    Y = measure[0].Yccd;
     if (USE_FIXED_PIXCOORDS) {
-      X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
-      Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
-    } else {
-      X = measure[0].Xccd;
-      Y = measure[0].Yccd;
+      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
+        float dX = measure[0].Xfix - measure[0].Xccd;
+        float dY = measure[0].Yfix - measure[0].Yccd;
+        if (hypot(dX,dY) < 2.0) {
+          X = measure[0].Xfix;
+          Y = measure[0].Yfix;
+        } 
+      }
     }
     n = measure[0].averef;
@@ -560 +586 @@
 
   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);
+  if (VERBOSE2 && (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);
 
@@ -609 +635 @@
     Measure *measure = &catalog[c].measure[m];
 
+    X = measure[0].Xccd;
+    Y = measure[0].Yccd;
     if (USE_FIXED_PIXCOORDS) {
-      X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
-      Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
-    } else {
-      X = measure[0].Xccd;
-      Y = measure[0].Yccd;
+      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
+        float dX = measure[0].Xfix - measure[0].Xccd;
+        float dY = measure[0].Yfix - measure[0].Yccd;
+        if (hypot(dX,dY) < 2.0) {
+          X = measure[0].Xfix;
+          Y = measure[0].Yfix;
+        } 
+      }
     }
     n = measure[0].averef;
@@ -679 +710 @@
 
     /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    raw[i].X = measure[0].Xccd;
+    raw[i].Y = measure[0].Yccd;
     if (USE_FIXED_PIXCOORDS) {
-      raw[i].X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd;
-      raw[i].Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd;
-    } else {
-      raw[i].X = measure[0].Xccd;
-      raw[i].Y = measure[0].Yccd;
-    }
-
+      if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) {
+        float dX = measure[0].Xfix - measure[0].Xccd;
+        float dY = measure[0].Yfix - measure[0].Yccd;
+        if (hypot(dX,dY) < 2.0) {
+          raw[i].X = measure[0].Xfix;
+          raw[i].Y = measure[0].Yfix;
+        } 
+      }
+    }
     raw[i].Mag  = measure[0].M;
     raw[i].dMag = measure[0].dM;
@@ -738 +773 @@
   StarData *ref;
 
+  int Nsecfilt = GetPhotcodeNsecfilt();
+
   mosaic = NULL;
   moscoords = NULL;
@@ -765 +802 @@
     ref[i].dMag = measure[0].dM;
     ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS);
+
+    if ((DCR_NSEC_POS >= 0) && (DCR_NSEC_NEG >= -1)) {
+      ref[i].Color = catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_POS].M - catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_NEG].M;
+    }
 
     ref[i].mask = FALSE;
@@ -859 +900 @@
     N = 0;
     for (k = 0; k < catalog[0].average[j].Nmeasure; k++) {
-      //reset flag on each invocation
+      // reset flag on each invocation
       measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
 
@@ -898 +939 @@
   double *R, *D, *dR, *dD, *d2;
   StatType statsR, statsD;
+
+  // XXX we are not going to use this for now
+  return;
 
   Nsecfilt = GetPhotcodeNsecfilt();
@@ -1063 +1107 @@
     } else {
       code = GetPhotcodebyCode (measure[0].photcode);
+      if (!code) return FALSE;
       mask = code[0].astromBadMask;
     }
@@ -1129 +1174 @@
     } else {
       code = GetPhotcodebyCode (measure[0].photcode);
+      if (!code) return FALSE;
       mask = code[0].astromBadMask;
     }

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

@@ -126 +126 @@
               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 (UPDATE)        { snprintf (tmpline, 1024, "%s -update",         command);                    strcpy (command, tmpline); }
-
-    if (RESET_BAD_IMAGES) { snprintf (tmpline, 1024, "%s -reset-bad-images", 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 (FIT_MODE == FIT_PM_ONLY)         strextend (command, "-pm");
+    if (FIT_MODE == FIT_PAR_ONLY)        strextend (command, "-par");
+    if (FIT_MODE == FIT_PM_AND_PAR)      strextend (command, "-pmpar");
+
+    if (VERBOSE)       strextend (command, "-v");
+    if (VERBOSE2)      strextend (command, "-vv");
+    if (RESET)         strextend (command, "-reset");
+    if (UPDATE)        strextend (command, "-update");
+
+    if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images");
+
+    if (ImagSelect)    strextend (command, "-instmag %f %f", ImagMin, ImagMax);
+    if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue);
     
-    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 (USE_BASIC_CHECK) strextend (command, "-basic-image-search");
+    if (FlagOutlier)     strextend (command, "-clip %d", CLIP_THRESH);
+    if (ExcludeBogus)    strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     
-    if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command);     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); }
+    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); 
+
+    if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
+    if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
+    if (PhotFlagSelect)     strextend (command, "+photflags"); 
+    if (PhotFlagBad)        strextend (command, "+photflagbad %d", PhotFlagBad);
+    if (PhotFlagPoor)       strextend (command, "+photflagpoor %d", PhotFlagPoor);
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
 
-    if (MinBadQF > 0.0)    { snprintf (tmpline, 1024, "%s -min-bad-psfqf %f",        command, MinBadQF);      strcpy (command, tmpline); }
-    if (MaxMeanOffset != 10.0) { snprintf (tmpline, 1024, "%s -max-mean-offset  %f", command, MaxMeanOffset); strcpy (command, tmpline); }
+    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
+      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); 
+    }
+
+    if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF); 
+    if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
 
     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); 
+      strextend (command, "-time %s %s", tstart, tstop); 
       free (tstart);
       free (tstop);
-      strcpy (command, tmpline); 
     }
     fprintf (stderr, "command: %s\n", command);

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

@@ -10 +10 @@
 static double *T;
 static double *dT;
+static double *C;
 
 void initObjectData (Catalog *catalog, int Ncatalog) {
@@ -36 +37 @@
   ALLOCATE (pX, double, MAX (1, Nmax));
   ALLOCATE (pY, double, MAX (1, Nmax));
+
+  ALLOCATE (C, double, MAX (1, Nmax));
 }  
 
@@ -54 +57 @@
   free (pX);
   free (pY);
+
+  free (C);
 }  
 
@@ -61 +66 @@
 
   off_t j, k, m;
-  int i, N, Nsecfilt, mode, result, status, XVERB;
+  int i, N, Nc, Nsecfilt, mode, result, status, XVERB;
   StatType statsR, statsD;
   Coords coords;
@@ -75 +80 @@
   memset (&fitPAR, 0, sizeof(fitPAR));
   initObjectData (catalog, Ncatalog);
+
+  int setRefColor = areImagesMatched();
 
   /* project coordinates to a plane centered on the object with units of arcsec */
@@ -112 +119 @@
       fitPAR.chisq = NAN;
 
+      // if we fail to fit the astrometry for some reason, we need to set/reset these
+      catalog[i].average[j].flags |= ID_STAR_NO_ASTROM;
+      catalog[i].average[j].ChiSqAve  = NAN;
+      catalog[i].average[j].ChiSqPM   = NAN;
+      catalog[i].average[j].ChiSqPar  = NAN;
+      catalog[i].average[j].Npos = 0;
+
       if (catalog[i].average[j].Nmeasure == 0) {
-          continue;
-      }
-
+        continue;
+      }
+
+      Nc = 0;
       N = 0;
       m = catalog[i].average[j].measureOffset;
@@ -138 +153 @@
         }
 
+        // reset the bit to note that a detection was used (or not)
+        measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+        if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
+
         // does the measurement pass the supplied filtering constraints?
         // MeasFilterTestTiny does not test psfQF
+        // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
         int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
         if (!keepMeasure) {
-          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
-          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
 
+        // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
+        if (ExcludeBogus) {
+          double Ri = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+          double Di = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+          coords.crval1 = catalog[i].average[j].R;
+          coords.crval2 = catalog[i].average[j].D;
+          double Xi, Yi;
+          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
+          double radius = hypot(Xi, Yi);
+          if (radius > ExcludeBogusRadius) {
+              measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
+              if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
+              continue;
+          }
+          measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
+        }
+
         // outlier rejection
-        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; }
+        if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
           continue;
         }
-
-        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]);
@@ -170 +201 @@
         // 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;
-          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
-          continue;
-        }
-        if (isnan(dY[N])) {
-          measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
-          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
-          continue;
-        }
+        if (isnan(dX[N])) continue;
+        if (isnan(dY[N])) continue;
 
         // add systematic error in quadrature, if desired
@@ -201 +224 @@
         dD[N] = dY[N] / 3600.0;
 
+        if (setRefColor) {
+          float color = getColor (m+k, i);
+          if (!isnan(color)) {
+            C[Nc] = color;
+            Nc++;
+          }
+        }
+
+        measure[k].dbFlags |= ID_MEAS_USED_OBJ;
+        if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }
+
         N++;
-      }
+      } // loop over measurements : catalog[i].average[j].Nmeasure 
 
       // if we have too few good detections for the desired fit, or too limited a
       // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
       // for even the average position, consider including the lower-quality detections?
-
-      catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM;
 
       // find Tmin & Tmax from the list of accepted measurements
@@ -231 +263 @@
       // too few measurements for average position (require 2 values)
       if (N < SRC_MEAS_TOOFEW) {
-        // XXX need to define PHOTOM and ASTROM object flags
-        // XXX reset the average value fields?
-        catalog[i].average[j].flags |= ID_STAR_NO_ASTROM;
-        catalog[i].average[j].ChiSqAve  = NAN;
-        catalog[i].average[j].ChiSqPM   = NAN;
-        catalog[i].average[j].ChiSqPar  = NAN;
         if (N < 2) continue;
       }
@@ -268 +294 @@
         catalog[i].average[j].flags |= ID_STAR_FIT_PM;
         Npm ++;
+
+        // XXX a hard-wired hack...
+        if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
+          mode = FIT_AVERAGE;
+          catalog[i].average[j].flags |= ID_STAR_BAD_PM;
+        }
       }
 
       if (mode == FIT_PM_AND_PAR) {
-        // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__);
-
         float pXmin = +2.0;
         float pXmax = -2.0;
@@ -296 +326 @@
           catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
           Npar ++;
+
+          // XXX a hard-wired hack...
+          if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
+            mode = FIT_AVERAGE;
+            catalog[i].average[j].flags |= ID_STAR_BAD_PM;
+          }
         } else {
           // need to set mode = FIT_PM_ONLY if we do not fit for parallax
@@ -321 +357 @@
         catalog[i].average[j].flags |= ID_STAR_FIT_AVE;
         Nave ++;
+      }
+
+      if (setRefColor) {
+        dsort (C, Nc);
+        float colorMedian = (Nc > 0) ? C[(int)(0.5*Nc)] : NAN;
+        catalog[i].average[j].refColor = colorMedian;
       }
 
@@ -411 +453 @@
       catalog[i].average[j].Trange = (Trange * 86400 * 365.25);
       catalog[i].average[j].Npos = fit.Nfit;
+
+      // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
+      catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM;
+
       if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
                           catalog[i].average[j].R,

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

@@ -179 +179 @@
   FlagOutlier = FALSE;
   if ((N = get_argument (argc, argv, "-clip"))) {
+    fprintf (stderr, "-clip is currently disabled\n");
     remove_argument (N, &argc, argv);
     CLIP_THRESH = atof (argv[N]);
     remove_argument (N, &argc, argv);
     FlagOutlier = TRUE;
+  }
+
+  ExcludeBogus = FALSE;
+  ExcludeBogusRadius = 0.0;
+  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
+    remove_argument (N, &argc, argv);
+    ExcludeBogusRadius = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ExcludeBogus = TRUE;
   }
 
@@ -284 +294 @@
   }
 
+  DCR_COLOR_POS = NULL;
+  DCR_COLOR_NEG = NULL;
+  if ((N = get_argument (argc, argv, "-dcr-color"))) {
+    remove_argument (N, &argc, argv);
+    DCR_COLOR_POS = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    DCR_COLOR_NEG = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PHOTCODE_RESET_LIST = NULL;
+  if ((N = get_argument (argc, argv, "-reset-to-photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_RESET_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   PHOTCODE_KEEP_LIST = NULL;
   if ((N = get_argument (argc, argv, "+photcode"))) {
@@ -291 +318 @@
   }
 
-  PHOTCODE_SKIP_LIST = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9");
+  char *SuperCOSMOS_SKIP = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9");
+
+  PHOTCODE_SKIP_LIST = NULL;
   if ((N = get_argument (argc, argv, "-photcode"))) {
     remove_argument (N, &argc, argv);
-    char *tmp1 = strcreate(argv[N]);
-
-    int Ntotal = strlen(tmp1) + strlen(PHOTCODE_SKIP_LIST) + 5;
-
-    char *tmp2 = NULL;
-    ALLOCATE (tmp2, char, Ntotal);
-    snprintf (tmp2, Ntotal, "%s,%s", PHOTCODE_SKIP_LIST, tmp1);
-
-    free (tmp1);
-    free (PHOTCODE_SKIP_LIST);
-
-    PHOTCODE_SKIP_LIST = tmp2;
+    char *RawSkip = strcreate(argv[N]);
+
+    char *GotSkip = strstr (RawSkip, SuperCOSMOS_SKIP);
+    if (!GotSkip) {
+      int Ntotal = strlen(RawSkip) + strlen(SuperCOSMOS_SKIP) + 5;
+      ALLOCATE (PHOTCODE_SKIP_LIST, char, Ntotal);
+      snprintf (PHOTCODE_SKIP_LIST, Ntotal, "%s,%s", SuperCOSMOS_SKIP, RawSkip);
+      free (RawSkip);
+      free (SuperCOSMOS_SKIP);
+    } else {
+      PHOTCODE_SKIP_LIST = RawSkip;
+    }
     remove_argument (N, &argc, argv);
   }
@@ -582 +611 @@
   FlagOutlier = FALSE;
   if ((N = get_argument (argc, argv, "-clip"))) {
+    fprintf (stderr, "-clip is currently disabled\n");
     remove_argument (N, &argc, argv);
     CLIP_THRESH = atof (argv[N]);
     remove_argument (N, &argc, argv);
     FlagOutlier = TRUE;
+  }
+
+  ExcludeBogus = FALSE;
+  ExcludeBogusRadius = 0.0;
+  if ((N = get_argument (argc, argv, "-exclude-bogus"))) {
+    remove_argument (N, &argc, argv);
+    ExcludeBogusRadius = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ExcludeBogus = TRUE;
   }
 
@@ -655 +694 @@
     remove_argument (N, &argc, argv);
     TimeSelect = TRUE;
+  }
+
+  DCR_COLOR_POS = NULL;
+  DCR_COLOR_NEG = NULL;
+  if ((N = get_argument (argc, argv, "-dcr-color"))) {
+    remove_argument (N, &argc, argv);
+    DCR_COLOR_POS = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    DCR_COLOR_NEG = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PHOTCODE_RESET_LIST = NULL;
+  if ((N = get_argument (argc, argv, "-reset-to-photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_RESET_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -861 +917 @@
   exit (2);
 }
+
+int strextend (char *input, char *format,...) {
+
+  char tmpextra[1024], tmpline[1024];
+  va_list argp;
+
+  va_start (argp, format);
+  vsnprintf (tmpextra, 1024, format, argp);
+  snprintf (tmpline, 1024, "%s %s", input, tmpextra);
+  strcpy (input, tmpline);
+
+  return TRUE;
+}

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

@@ -1 +1 @@
 # include "relastro.h"
 
+static int Nkeep1 = 0;
+static int Nkeep2 = 0;
 static int Nskip1 = 0;
 static int Nskip2 = 0;
-static unsigned int Tref = 1323353985;
-static short Cref = 10355;
+
+FILE *fbogus = NULL;
+static int NskipBogus = 0;
+
+// test image: 2013/06/15,13:25:51, GPC1.r.XY50
+static int CHECK_TEST_IMAGE = FALSE;
+// static unsigned int Tref = 1378812312; 
+// static short Cref = 10001;
+
+static unsigned int Tref = 1379570672;
+static short Cref = 10341;
+
+int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog);
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
@@ -11 +24 @@
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   int Nsecfilt;
+  Coords coords;
+
+  /* for outlier rejection, 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, "DEC--SIN");
 
   // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
@@ -34 +59 @@
     ID_STAR_NO_ASTROM ; 
 
+  if (VERBOSE2 && ExcludeBogus && (fbogus == NULL)) {
+    char name[1024];
+    snprintf (name, 1024, "%s/bogus.%02d.dat", CATDIR, HOST_ID);
+    fbogus = fopen (name, "w");
+    if (!fbogus) {
+      fprintf (stderr, "trouble opening bogus detection dump : %s\n", name);
+    }
+  }
+
   /* exclude stars not in range or with too few measurements */
   for (i = 0; i < catalog[0].Naverage; i++) {
@@ -60 +94 @@
       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)) {
+        if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
           Nskip1 ++;
         }
         continue;
       }
-
-      // filter out outliers
-      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
+      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+        Nkeep1 ++;
+      }
+
+      // filter out outliers - these are detections inconsistent with the offset distribution
+      // XXX disable this for now
+      if (FALSE && 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)) {
@@ -76 +114 @@
       catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
 
+      // exclude bogus 
+      if (ExcludeBogus) {
+          double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
+          double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
+          coords.crval1 = catalog[0].average[i].R;
+          coords.crval2 = catalog[0].average[i].D;
+          double Xi, Yi;
+          RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
+          double radius = hypot(Xi, Yi);
+          if (radius > ExcludeBogusRadius) {
+              NskipBogus ++;
+              if (VERBOSE2) {
+                FILE *foutput = fbogus ? fbogus : stderr;
+                char *date = ohana_sec_to_date(catalog[0].measure[offset].t);
+                fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date);
+                free (date);
+              }
+              continue;
+          }
+      }
+
       // re-assess on each run of relastro if a measurement should be used
 
@@ -87 +146 @@
       // but before the final average properties are calculated, these measurements may be
       // allowed.
+
+      if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) {
+        Nkeep2 ++;
+      }
 
       CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
@@ -125 +188 @@
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
-    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
+    LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog);
   } else {
     if (VERBOSE2) {
@@ -134 +197 @@
     }
   }
+  if (CHECK_TEST_IMAGE && (Nkeep1 + Nkeep2 + Nskip1 + Nskip2 > 0)) {
+    fprintf (stderr, "kept %d %d, skipped %d %d\n", Nkeep1, Nkeep2, Nskip1, Nskip2);
+  }
+
   return (TRUE);
 }
 
 void bcatalog_show_skips () {
-  fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
+  if (ExcludeBogus) {
+    if (fbogus) fclose (fbogus);
+    fprintf (stderr, "NskipBogus: %d\n", NskipBogus);
+    // fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2);
+    // fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2);
+  }
 }
 
@@ -155 +227 @@
 # undef COMPARE
 
+}
+
+/* generate a grid in a locally-projected coordinate system, try to select average entries 
+   from each grid cell in decending Nmeasure order.
+ */ 
+int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) {
+
+  off_t i, j;
+  int ix, iy;
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  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);
+
+  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 (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // generate a grid in locally projected space
+  double Rc = 0.5*(Rmin + Rmax);
+  double Dc = 0.5*(Dmin + Dmax);
+
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  Coords coords;
+  coords.crval1 = Rc;
+  coords.crval2 = Dc;
+  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, "DEC--SIN");
+
+  // convert all average R,D values to X,Y:
+  double *X, *Y;
+  ALLOCATE (X, double, Naverage);
+  ALLOCATE (Y, double, Naverage);
+  float Xmin = +10000.0, Ymin = +10000.0;
+  float Xmax = -10000.0, Ymax = -10000.0;
+  for (i = 0; i < Naverage; i++) {
+    X[i] = NAN;
+    Y[i] = NAN;
+    // skip any stars which are outside of nominal catalog range
+    if (subcatalog[0].average[i].R < Rmin) continue;
+    if (subcatalog[0].average[i].R > Rmax) continue;
+    if (subcatalog[0].average[i].D < Dmin) continue;
+    if (subcatalog[0].average[i].D > Dmax) continue;
+    RD_to_XY (&X[i], &Y[i], subcatalog[0].average[i].R, subcatalog[0].average[i].D, &coords);
+    Xmin = MIN (Xmin, X[i]);
+    Xmax = MAX (Xmax, X[i]);
+    Ymin = MIN (Ymin, Y[i]);
+    Ymax = MAX (Ymax, Y[i]);
+  }
+
+  // how many grid cells? what is the grid spacing? 
+  float dX = Xmax - Xmin;
+  float dY = Ymax - Ymin;
+
+  // *** XXX for the moment, I'm using a hard-wired cell size (200 arcsec ~ 3.3 arcmin)
+  int NX = (int)(dX / 200) + 1;
+  int NY = (int)(dY / 200) + 1;
+  // fprintf (stderr, "Density Grid: %d x %d\n", NX, NY);
+  // XXX check that NX,NY are sensible (5 degrees / 200 arcsec seems like the absolute max)
+  if (NX > 1000) { 
+    fprintf (stderr, "serious problem with %s: NX = %d\n", subcatalog[0].filename, NX); 
+    exit (3); 
+  }
+  if (NY > 1000) { 
+    fprintf (stderr, "serious problem with %s: NY = %d\n", subcatalog[0].filename, NY); 
+    exit (3); 
+  }
+  
+  // kind of ugly : generate a grid of index, Nmeasure arrays
+  // to be filled below (I also need NN and Nn to track the number of 
+  // entries in each).
+  int    **NN_grid;
+  int    **Nn_grid;
+  int   ***Nm_grid;
+  off_t ***idxgrid;
+  ALLOCATE (NN_grid, int *, NX);
+  ALLOCATE (Nn_grid, int *, NX);
+  ALLOCATE (Nm_grid, int **, NX);
+  ALLOCATE (idxgrid, off_t **, NX);
+
+  for (ix = 0; ix < NX; ix++) {
+    ALLOCATE (NN_grid[ix], int, NY);
+    ALLOCATE (Nn_grid[ix], int, NY);
+    ALLOCATE (Nm_grid[ix], int *, NY);
+    ALLOCATE (idxgrid[ix], off_t *, NY);
+    for (iy = 0; iy < NY; iy++) {
+      Nn_grid[ix][iy] = 0;
+      NN_grid[ix][iy] = 100;
+      ALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
+      ALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
+    }
+  }
+
+  // assign all of the average entries to a grid cell
+  for (i = 0; i < Naverage; i++) {
+    if (isnan(X[i])) continue;
+    if (isnan(Y[i])) continue;
+    ix = (X[i] - Xmin) / 200.0;
+    iy = (Y[i] - Ymin) / 200.0;
+    int Nn = Nn_grid[ix][iy];
+    Nm_grid[ix][iy][Nn] = subcatalog[0].average[i].Nmeasure;
+    
+    // if we are resetting to a given photcode, we need to have that photcode...
+    if (NphotcodesReset) {
+      int k;
+      int foundReset = FALSE;
+      int m = subcatalog[0].average[i].measureOffset;
+      MeasureTiny *measure = &subcatalog[0].measureT[m];
+      for (j = 0; (j < subcatalog[0].average[i].Nmeasure) && !foundReset; j++) {
+        if (CHECK_TEST_IMAGE && (abs(measure[j].t - Tref) < 10) && (measure[j].photcode == Cref)) {
+          fprintf (stderr, ".");
+        }
+        for (k = 0; (k < NphotcodesReset) && !foundReset; k++) {
+          if (photcodesReset[k][0].code == measure[j].photcode) foundReset = TRUE;
+        }
+      }
+      if (!foundReset) {
+        Nm_grid[ix][iy][Nn] = 0;
+      }
+    }
+
+    idxgrid[ix][iy][Nn] = i;
+    Nn_grid[ix][iy] ++;
+    if (Nn_grid[ix][iy] >= NN_grid[ix][iy]) {
+      NN_grid[ix][iy] += 100;
+      REALLOCATE (Nm_grid[ix][iy], int,   NN_grid[ix][iy]);
+      REALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]);
+    }
+  }
+    
+  // sort all of the grid cells
+  for (ix = 0; ix < NX; ix++) {
+    for (iy = 0; iy < NY; iy++) {
+      sort_by_Nmeasure (Nm_grid[ix][iy], idxgrid[ix][iy], Nn_grid[ix][iy]);
+      NN_grid[ix][iy] = 0; // I'm going to use this array to track which element I've already selected
+    }
+  }
+
+  // cycle over the grid until we ready Nmax
+  off_t *keepidx = NULL;
+  ALLOCATE (keepidx, off_t, Naverage);
+  memset (keepidx, 0, Naverage*sizeof(off_t));
+  int Nkeep = 0;
+
+  for (i = 0; (i < 20) && (Nkeep < Nmax); i++) {
+    for (ix = 0; (ix < NX) && (Nkeep < Nmax); ix++) {
+      for (iy = 0; (iy < NY) && (Nkeep < Nmax); iy++) {
+        if (NN_grid[ix][iy] >= Nn_grid[ix][iy]) continue; // all used up!
+        int Nn = NN_grid[ix][iy];
+        keepidx[Nkeep] = idxgrid[ix][iy][Nn];
+        Nkeep ++;
+        NN_grid[ix][iy] ++;
+      }
+    }
+  }
+
+  // count the number of measurements this selection will yield
+  off_t ave, NMEASURE = 0;
+  for (i = 0; i < Nkeep; i++) {
+    ave = keepidx[i];
+    NMEASURE += subcatalog[0].average[ave].Nmeasure;
+  }
+
+  // test catID : 37262 37261 37257 37258
+  int dumpit = FALSE;
+  dumpit |= (oldcatalog[0].catID == 37007);
+  // dumpit |= (oldcatalog[0].catID == 37261);
+  // dumpit |= (oldcatalog[0].catID == 37257);
+  // dumpit |= (oldcatalog[0].catID == 37258);
+  if (dumpit) {
+    char name[64];
+    snprintf (name, 64, "cat.%05d.dump.dat", oldcatalog[0].catID);
+    FILE *fdump = fopen (name, "w");
+    for (i = 0; i < Nkeep; i++) {
+      ave = keepidx[i];
+      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
+    }
+    fclose (fdump);
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.average,  Average,     Nkeep);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt,     Nkeep * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nkeep selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nkeep; i++) {
+    ave = keepidx[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) {
+    char *basename = filebasename (oldcatalog[0].filename);
+    fprintf (stderr, "limited to %d ("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", 
+             Nkeep, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
+    free (basename);
+  }
+
+  free (X);
+  free (Y);
+
+  for (ix = 0; ix < NX; ix++) {
+    for (iy = 0; iy < NY; iy++) {
+      free (Nm_grid[ix][iy]);
+      free (idxgrid[ix][iy]);
+    }
+    free (NN_grid[ix]);
+    free (Nn_grid[ix]);
+    free (Nm_grid[ix]);
+    free (idxgrid[ix]);
+  }
+  free (NN_grid);
+  free (Nn_grid);
+  free (Nm_grid);
+  free (idxgrid);
+
+  free (keepidx);
+
+  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 = Nkeep;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
+
+  return (TRUE);
 }
 
@@ -205 +543 @@
     NMEASURE += subcatalog[0].average[ave].Nmeasure;
   }
+
+# if (0)
+  if (oldcatalog[0].catID == 59962) {
+    FILE *fdump = fopen ("cat.dump.dat", "w");
+    for (i = 0; i < Nmax; i++) {
+      ave = index[i];
+      fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure);
+    }
+    fclose (fdump);
+  }
+# endif
 
   // allocate the output data 

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

@@ -91 +91 @@
 }
 
+double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double ra;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  ra = average[0].R - measure[0].dR / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (ra);
+}
+
+double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double dec;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  dec = average[0].D - measure[0].dD / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (dec);
+}
+
 int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
 
@@ -144 +184 @@
   return (TRUE);
 }
+

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

@@ -8 +8 @@
   if (RELASTRO_OP == OP_MERGE_SOURCE) return;
 
-  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);
+  if (DCR_COLOR_POS)  fprintf (stderr, "DCR_COLOR_POS:  %s\n", DCR_COLOR_POS);
+  if (DCR_COLOR_NEG)  fprintf (stderr, "DCR_COLOR_NEG:  %s\n", DCR_COLOR_NEG);
 
-  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);
+  if (PHOTCODE_KEEP_LIST)  fprintf (stderr, "PHOTCODE_KEEP_LIST:  %s\n", PHOTCODE_KEEP_LIST);
+  if (PHOTCODE_SKIP_LIST)  fprintf (stderr, "PHOTCODE_SKIP_LIST:  %s\n", PHOTCODE_SKIP_LIST);
+  if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST);
+  if (PHOTCODE_A_LIST)     fprintf (stderr, "PHOTCODE_A_LIST:     %s\n", PHOTCODE_A_LIST);
+  if (PHOTCODE_B_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);
+  photcodesReset  = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset,  FALSE);
+  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST,     &NphotcodesGroupA, TRUE);
+  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST,     &NphotcodesGroupB, TRUE);
+
+  DCR_NSEC_POS = DCR_NSEC_NEG = -1;
+  if (DCR_COLOR_POS) {
+    DCR_PHOTCODE_POS = GetPhotcodebyName (DCR_COLOR_POS);
+    if (!DCR_PHOTCODE_POS) {
+      fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_POS);
+      exit (1);
+    }
+    DCR_NSEC_POS = GetPhotcodeNsec (DCR_PHOTCODE_POS[0].code);
+  }
+  if (DCR_COLOR_NEG) {
+    DCR_PHOTCODE_NEG = GetPhotcodebyName (DCR_COLOR_NEG);
+    if (!DCR_PHOTCODE_NEG) {
+      fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_NEG);
+      exit (1);
+    }
+    DCR_NSEC_NEG = GetPhotcodeNsec (DCR_PHOTCODE_NEG[0].code);
+  }
 
   initstats (STATMODE);
@@ -42 +65 @@
   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);
+  if (PHOTCODE_KEEP_LIST)  fprintf (stderr, "PHOTCODE_KEEP_LIST:  %s\n", PHOTCODE_KEEP_LIST);
+  if (PHOTCODE_SKIP_LIST)  fprintf (stderr, "PHOTCODE_SKIP_LIST:  %s\n", PHOTCODE_SKIP_LIST);
+  if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST);
+  if (PHOTCODE_A_LIST)     fprintf (stderr, "PHOTCODE_A_LIST:     %s\n", PHOTCODE_A_LIST);
+  if (PHOTCODE_B_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);
+  photcodesKeep   = ParsePhotcodeList (PHOTCODE_KEEP_LIST,  &NphotcodesKeep,   FALSE);
+  photcodesSkip   = ParsePhotcodeList (PHOTCODE_SKIP_LIST,  &NphotcodesSkip,   FALSE);
+  photcodesReset  = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset,  FALSE);
+  photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST,     &NphotcodesGroupA, TRUE);
+  photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST,     &NphotcodesGroupB, TRUE);
 
   initstats (STATMODE);

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

@@ -1 +1 @@
 # include "relastro.h"
 # define DEBUG 0
-
-int strextend (char *input, char *format,...) {
-
-  char tmpextra[1024], tmpline[1024];
-  va_list argp;
-
-  va_start (argp, format);
-  vsnprintf (tmpextra, 1024, format, argp);
-  snprintf (tmpline, 1024, "%s %s", input, tmpextra);
-  strcpy (input, tmpline);
-
-  return TRUE;
-}
 
 int launch_region_hosts (RegionHostTable *regionHosts) {
@@ -88 +75 @@
     if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); 
     if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
+    if (PHOTCODE_RESET_LIST) strextend (command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST);
 
     if (MaxDensityUse)      strextend (command, "-max-density %f", MaxDensityValue);
     if (ImagSelect)         strextend (command, "-instmag %f %f", ImagMin, ImagMax);
+    if (ExcludeBogus)       strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
+
+    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
+      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); 
+    }
+
     if (PhotFlagSelect)     strextend (command, "+photflags"); 
     if (PhotFlagBad)        strextend (command, "+photflagbad %d", PhotFlagBad); 
@@ -97 +91 @@
     if (MinBadQF > 0.0)        strextend (command, "-min-bad-psfqf %f", MinBadQF);
     if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset  %f", MaxMeanOffset);
+
+    strextend (command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW);
+    strextend (command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM);
+    strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX);
+    strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS);
+
+    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
 
     if (TimeSelect) { 

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

@@ -1 +1 @@
 # include "relastro.h"
 
-Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog);
+Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile);
 void bcatalog_show_skips ();
 
-Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) {
-
-  int i, j, k, m, Nstar;
+Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile) {
+
+  int i, j, Nstar;
+  // int k, m;
   Catalog *catalog, *pcatalog, tcatalog;
 
   // XXX need to decide how to determine PARALLEL mode...
   if (PARALLEL && !hostID) {
-    catalog = load_catalogs_parallel (skylist, Ncatalog);
+    catalog = load_catalogs_parallel (skylist, Ncatalog, syncfile);
     return catalog;
   }
@@ -73 +74 @@
       if (RESET) {
         for (j = 0; j < catalog[Ncat].Naverage; j++) {
+# if (0)          
           catalog[Ncat].average[j].flags = 0;
           m = catalog[Ncat].average[j].measureOffset;
@@ -78 +80 @@
             catalog[Ncat].measure[m+k].dbFlags = 0;
           }
+# endif
         }
       }
@@ -84 +87 @@
   }
 
-  // XXX TEST : bcatalog_show_skips();
+  bcatalog_show_skips();
 
   Nstar = 0;
@@ -113 +116 @@
 // CATDIR is supplied globally
 # define DEBUG 1
-Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
+Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile) {
 
   char uniquer[12];
@@ -150 +153 @@
               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 (FIT_MODE == FIT_PM_ONLY)         strextend (command, "-pm");
+    if (FIT_MODE == FIT_PAR_ONLY)        strextend (command, "-par");
+    if (FIT_MODE == FIT_PM_AND_PAR)      strextend (command, "-pmpar");
+
+    if (VERBOSE)       strextend (command, "-v");
+    if (VERBOSE2)      strextend (command, "-vv");
+    if (RESET)         strextend (command, "-reset");
+    if (ImagSelect)    strextend (command, "-instmag %f %f", ImagMin, ImagMax);
+    if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue);
+    if (FlagOutlier)   strextend (command, "-clip %d", CLIP_THRESH);
+    if (ExcludeBogus)  strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     
-    if (FlagOutlier)   { snprintf (tmpline, 1024, "%s -clip %d",        command, CLIP_THRESH);       strcpy (command, tmpline); }
-    
-    if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command);     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); }
+    if (DCR_COLOR_POS && DCR_COLOR_NEG) {
+      strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); 
+    }
+
+    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
+    if (PHOTCODE_KEEP_LIST)  strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST);
+    if (PHOTCODE_SKIP_LIST)  strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST);
+    if (PhotFlagSelect)      strextend (command, "+photflags");
+    if (PhotFlagBad)         strextend (command, "+photflagbad %d", PhotFlagBad); 
+    if (PhotFlagPoor)        strextend (command, "+photflagpoor %d", PhotFlagPoor);
     // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
 
@@ -175 +180 @@
       char *tstart = ohana_sec_to_date (TSTART);
       char *tstop  = ohana_sec_to_date (TSTOP);
-      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); 
+      strextend (command, "-time %s %s", tstart, tstop); 
       free (tstart);
       free (tstop);
-      strcpy (command, tmpline); 
     }
 
@@ -211 +215 @@
   }
 
+  // update syncfile here (save lots of I/O time)
+
+  // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
+  // relastro_client requests.  I need to have the master mediate this.  the master
+  // will not launch the next remote job until this one says it is done
+  if (syncfile) {
+    update_sync_file (syncfile, 1);
+  }
+
   // 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

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

@@ -33 +33 @@
       // USAGE: relastro_client -load-objects
       int Ncatalog;
-      Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR);
+      Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR, NULL);
       if (!catalog) {
         fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR);

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

@@ -28 +28 @@
   // XXX pass in the image table
   // XXX who carries the image grid?
-  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL);
+
+  // photcodesKeep is used here to allow measurements from the images being calibrated
+  // note if -reset-to-photcode is selected, photocodesKeep is replaced with below with photcodesReset
+  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, NULL);
   MARKTIME("load catalog data: %f sec\n", dtime);
+
+  if (photcodesReset) {
+    photcodesKeep  = photcodesReset;
+    NphotcodesKeep = NphotcodesReset;
+  }
 
   if (Ncatalog == 0) {
@@ -50 +58 @@
   // set test points based on the starmap
   createStarMap (catalog, Ncatalog);
+
+  // XXX NOTE : for 2mass reset, photcodesKeep should now limit to 2MASS measurements
 
   /* major modes */

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

@@ -40 +40 @@
     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.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

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

@@ -45 +45 @@
 
   /* load catalog data from region files (hostID is 0 since we are not a client */
-  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL);
+  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
+  catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, syncfile);
   MARKTIME("-- load catalog data: %f sec\n", dtime);
+  free (syncfile);
   
-  // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
-  // relastro_client requests.  I need to have the master mediate this.  the master
-  // will not launch the next remote job until this one says it is done
-  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
-  update_sync_file (syncfile, 1);
-  
+  if (photcodesReset) {
+    photcodesKeep  = photcodesReset;
+    NphotcodesKeep = NphotcodesReset;
+  }
+
   // generate tables go from catID,objID -> catSeq,objSeq
   indexCatalogs (catalog, Ncatalog);