Index: /branches/eam_branches/relphot.20210521/include/relphot.h
===================================================================
--- /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/include/relphot.h	(revision 41625)
@@ -8,4 +8,11 @@
 
 # define NBOOTSTRAP 100
+
+// MEAS_BAD is used in ImageOps.c, MosaicOps.c, TGroupOps.c to skip specific measurements.
+// ID_MEAS_AREA is raised for measurements outside the user-specified region of the chips
+// ID_MEAS_NOCAL is raised for measurements outside the valid time range or for inactive photcodes
+// Note: these bits are only raised on the temporary MeasureTiny elements and are not saved
+// to the database files.
+# define MEAS_BAD (ID_MEAS_AREA | ID_MEAS_NOCAL)
 
 # ifndef MAX_INT
@@ -453,11 +460,8 @@
 char   MOSAICNAME[256];
 char   STATMODE[32];
-int    STAR_BAD;
-int    MEAS_BAD;
 int    STAR_TOOFEW;
 int    GRID_TOOFEW;
 int    IMAGE_TOOFEW;
 double IMAGE_GOOD_FRACTION;
-int    IMAGE_BAD;
 int    CALIBRATE_STACKS_AND_WARPS;
 
@@ -848,5 +852,5 @@
 int GetZptIteration (void);
 
-void SetZeroPointModes (void);
+void SetZeroPointModes (Catalog *catalog, int Ncatalog);
 int UseStandardOLS (ZptFitModeType mode);
 int GetActivePhotcodeIndex (int photcode);
Index: /branches/eam_branches/relphot.20210521/src/GridOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/GridOps.c	(revision 41625)
@@ -192,6 +192,5 @@
       for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) {
 
-	// XXX I need to skip bad measurements, but I am not certain this
-	// flag is set correctly.  see note below 
+	// skip measurements marked by AREA or TIME
 	if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue;
 
Index: /branches/eam_branches/relphot.20210521/src/ImageOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/ImageOps.c	(revision 41625)
@@ -385,5 +385,7 @@
   if (i == -1) return NAN;
 
-  if (image[i].flags & IMAGE_BAD) return NAN;  
+  // allow context dependence?: in relphot_images, we want to exclude detections from poor images
+  // but in the final pass, we want to allow even bad detections 
+  if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) return NAN;  
 
   switch (class) {
@@ -409,5 +411,4 @@
   if (i == -1) return (1000);
 
-  if (image[i].flags & IMAGE_BAD) return (1000);  
   distance = image[i].ubercalDist; // was dummy3 in structure
   return (distance);
@@ -689,6 +690,9 @@
     if (mark) {
       image[i].flags |= ID_IMAGE_PHOTOM_FEW;
+      image[i].McalPSF    = 0.0;
+      image[i].McalAPER   = 0.0;
+      image[i].McalChiSq  = NAN;
       image[i].dMcal      = NAN;
-      image[i].McalChiSq  = NAN;
+      image[i].dMagSys    = NAN;
       image[i].nFitPhotom = 0;
       Nfew ++;
@@ -894,5 +898,5 @@
   // measure stats for Mcal and dMcal
   for (i = N = 0; i < Nimage; i++) {
-    if (image[i].flags & IMAGE_BAD) continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -922,6 +926,15 @@
 
     mark = FALSE;
-    image[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
-    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset);
+
+    if (image[i].flags & ID_IMAGE_PHOTOM_FEW) {
+      mark = TRUE;
+    }
+    if (image[i].dMcal > MaxScatter) {
+      mark = TRUE;
+    }
+    if (fabs(image[i].McalPSF - MedOffset) > MaxOffset) {
+      mark = TRUE;
+    }
+
     if (mark) { 
       Nmark ++;
@@ -1089,5 +1102,5 @@
   n = 0;
   for (i = 0; i < Nimage; i++) {
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1138,5 +1151,5 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1171,5 +1184,5 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1204,5 +1217,5 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
Index: /branches/eam_branches/relphot.20210521/src/MagResidSave.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/MagResidSave.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/MagResidSave.c	(revision 41625)
@@ -215,5 +215,4 @@
     off_t c = ImageToCatalog[j];
       
-    // XXX drop this?
     if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
Index: /branches/eam_branches/relphot.20210521/src/MosaicOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41625)
@@ -780,5 +780,7 @@
   if (i == -1) return (0.0);
 
-  if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
+  // if the mosaic cannot be calibrated (too few measurements, skip it)
+  // XXX: if the mosaic is bad (few) but the chips are OK, I can used them.
+  // XXX: if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) return (NAN);  
   value = mosaic[i].McalPSF;
   return (value);
@@ -1330,205 +1332,4 @@
 }
 
-// When we rationalize the images/mosaics, we are driving the negative cloud images back
-// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
-// driven by the coupling of common stars. 
-// XXX this function is no longer used because we force significantly negative clouds to 0.0
-int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
-
-  double *imageOffset, **starOffset;
-  int **starNcount, *seclist;
-  int **Slist, *NSlist, *NSLIST;
-  int i, j, k, m, nMos, Ns, found;
-
-  off_t Nimage;
-  Image *image;
-
-  // set a test value for now
-  float CLOUD_TOLERANCE = 0.025;
-
-  if (!MOSAIC_ZEROPT) return (FALSE);
-  if (FREEZE_MOSAICS) return (FALSE);
-
-  image = getimages (&Nimage, NULL);
-
-  ALLOCATE (imageOffset, double, Nmosaic);
-  ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
-  ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
-  ALLOCATE (NSLIST, int,   Nmosaic); // number of Slist entries allocated
-  memset (Slist, 0, Nmosaic*sizeof(int *));
-
-  // find the images / mosaics with negative clouds and save their offset
-  for (i = 0; i < Nmosaic; i++) {
-  
-    NSlist[i] =   0;
-    NSLIST[i] = 100;
-    ALLOCATE (Slist[i], int, NSLIST[i]);
-
-    imageOffset[i] = 0.0;
-
-    if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-    if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) {
-      imageOffset[i] = -mosaic[i].McalPSF; 
-      // NOTE the negative sign: down below, we are going to add in the negative of McalPSF
-      // to this image, and the propagated mean values for other images
-      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-    if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-  }
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-  ALLOCATE (seclist, int, Nphotcodes);
-  for (Ns = 0; Ns < Nphotcodes; Ns ++) { 
-    int thisCode = photcodes[Ns][0].code;
-    seclist[Ns] = GetPhotcodeNsec(thisCode);
-  }
-
-  // allocate an array for star offsets
-  int Nstars = 0;
-  int NSTARS = 1000;
-  ALLOCATE (starOffset, double *, NSTARS);
-  ALLOCATE (starNcount, int *,    NSTARS);
-  memset (starOffset, 0, NSTARS*sizeof(double *));
-  memset (starNcount, 0, NSTARS*sizeof(int *));
-
-  // find the mean offset for each star
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      ALLOCATE (starOffset[Nstars], double, Nphotcodes);
-      ALLOCATE (starNcount[Nstars], int,    Nphotcodes);
-      memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double));
-      memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int));
-      
-      m = catalog[i].averageT[j].measureOffset;
-
-      // determine the mosaic for each measurement
-      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-
-	// skip unused measurements
-	if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-
-	// skip unused measurements
-
-	// XXX replace with:
-	// Ns = GetActivePhotcodeIndex (photcode);
-	// if (Ns < 0) continue;
-
-	int Nsec;
-	int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-	found = FALSE;
-	for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
-	  if (ecode == photcodes[Ns][0].code) {
-	    found = TRUE;
-	    break;
-	  }
-	}
-	if (!found) continue;
-	Nsec = seclist[Ns];
-
-	// bad stars for this secfilt
-	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-
-	// skip REF measurements (not tied to an image)
-	if (getImageEntry (m, i) < 0) continue;
-
-	// find the source of this measurement (skip unassigned measurements)
-	nMos = MeasureToMosaic[i][m];
-	if (nMos == -1) continue;
-
-	// XXX how can this not be true?  
-	if (mosaic[nMos].photcode != ecode) {
-	  fprintf (stderr, "*");
-	}
-
-	assert (Ns >= 0);
-	assert (Ns < Nphotcodes);
-
-	// accumulate the offsets from the negative cloud images (others have 0.0 value)
-	starOffset[Nstars][Ns] += imageOffset[nMos];
-	starNcount[Nstars][Ns] ++;
-	
-	// record the mosaic->star reference
-	Slist[nMos][NSlist[nMos]] = Nstars;
-	NSlist[nMos] ++;
-	if (NSlist[nMos] == NSLIST[nMos]) {
-	  NSLIST[nMos] += 100;
-	  REALLOCATE (Slist[nMos], int, NSLIST[nMos]);
-	}	  
-      }
-      Nstars ++;
-      if (Nstars == NSTARS) {
-	NSTARS += 1000;
-	REALLOCATE (starOffset, double *, NSTARS);
-	REALLOCATE (starNcount, int *,    NSTARS);
-	memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *));
-	memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *));
-      }
-    }
-  }
-
-  // find the mean offset of the images without negative clouds
-  for (i = 0; i < Nmosaic; i++) {
-
-    found = FALSE;
-    for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
-      if (mosaic[i].photcode == photcodes[Ns][0].code) {
-	found = TRUE;
-	break;
-      }
-    }
-    if (!found) {
-      fprintf (stderr, "invalid photcode for mosaic?\n");
-      abort();
-    }
-
-    // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF)
-    if (imageOffset[i] > 0.0) continue;
-
-    // we need to actually have cross-references to count
-    if (NSlist[i] < 2) continue;
-
-    float dM = 0.0;
-    for (j = 0; j < NSlist[i]; j++) {
-      k = Slist[i][j];
-      if (starNcount[k][Ns] > 1) {
-	dM += (starOffset[k][Ns] / starNcount[k][Ns]);
-      }
-    }
-    imageOffset[i] = dM / NSlist[i];
-    // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]);
-  }
-
-  // for (i = 0; i < Nmosaic; i++) {
-  //   fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]);
-  // } 
-
-  // apply all offset values to the mosaics
-  // find the images / mosaics with negative clouds and save their offset
-  for (i = 0; i < Nmosaic; i++) {
-    mosaic[i].McalPSF += imageOffset[i];
-  }
-
-  for (i = 0; i < Nstars; i++) {
-    free (starOffset[i]);
-    free (starNcount[i]);
-  }
-  free (starOffset);
-  free (starNcount);
-
-  free (seclist);
-  free (NSlist);
-  free (NSLIST);
-  for (i = 0; i < Nmosaic; i++){
-    free (Slist[i]);
-  }
-  free (Slist);
-  free (imageOffset);
-  return (TRUE);
-}
-
 StatType statsMosaicM (Catalog *catalog) {
   OHANA_UNUSED_PARAM(catalog);
@@ -1547,5 +1348,5 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1580,5 +1381,5 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1613,5 +1414,5 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD)  continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW))  continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1663,5 +1464,5 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1795,5 +1596,5 @@
       c = MosaicToCatalog[i][j];
       
-      if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
       // ave = catalog[c].measureT[m].averef;
Index: /branches/eam_branches/relphot.20210521/src/ResetOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ResetOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/ResetOps.c	(revision 41625)
@@ -75,4 +75,7 @@
     }
   }
+  if (RESET_FLATCORR) {
+    measure->Mflat = 0.0;
+  }
   measure->dbFlags &= ~photomBitsMeasure;
 }
@@ -88,5 +91,5 @@
   for (off_t i = 0; i < Nsubset; i++) {
     // reset these bits regardless (we will re-determine)
-    subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW);
+    subset[i].flags &= ~(ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW);
     if (RESET) {
       if (RESET_ZEROPTS) {
@@ -117,10 +120,12 @@
 
   for (off_t i = 0; i < catalog->Naverage; i++) {
-    catalog[0].average[i].psfQF     = NAN;	// force recalculation in setMrelCatalog
-    catalog[0].average[i].psfQFperf = NAN;	// force recalculation in setMrelCatalog
-    catalog[0].average[i].stargal   = NAN;	// force recalculation in setMrelCatalog
-    catalog[0].average[i].photFlagsUpper = 0;	// reset (will be re-calculated)
-    catalog[0].average[i].photFlagsLower = 0;	// reset (will be re-calculated)
-    catalog[0].average[i].NwarpOK        = 0;	// reset (will be re-calculated)
+    catalog->average[i].psfQF     = NAN;	// force recalculation in setMrelCatalog
+    catalog->average[i].psfQFperf = NAN;	// force recalculation in setMrelCatalog
+    catalog->average[i].stargal   = NAN;	// force recalculation in setMrelCatalog
+    catalog->average[i].photFlagsUpper = 0;	// reset (will be re-calculated)
+    catalog->average[i].photFlagsLower = 0;	// reset (will be re-calculated)
+    catalog->average[i].NwarpOK        = 0;	// reset (will be re-calculated)
+
+    // RESET_FLATCORR independent of RESET?
 
     if (!RESET) continue;
@@ -128,38 +133,61 @@
     catalog->average[i].flags    &= ~photomBitsAverage; // reset all photometry bits (but not astrom)
 
-    for (int Ns = 0; Ns < Nsecfilt; Ns++) {
+    // Reset secfilt values only for the active photcodes
+    for (int Ns = 0; Ns < Nphotcodes; Ns++) {
+      
+      int thisCode = photcodes[Ns][0].code;
+      int Nsec = GetPhotcodeNsec(thisCode);
+      off_t N = Nsecfilt*i+Nsec;
+      
+      dvo_secfilt_init (&catalog->secfilt[N], SECFILT_RESET_ALL);
+    }
+
+    // Reset measure values only for the active photcodes
+    
+    off_t m = catalog->average[i].measureOffset;
+    
+    for (off_t j = 0; j < catalog->average[i].Nmeasure; j++, m++) {
+      
+      /* select measurements for active photcodes */
+      PhotCode *code = GetPhotcodebyCode (catalog->measure[m].photcode);
+      if (!code) continue;
+      // photcode exists, but is it one of the active ones?
+      int found = FALSE;
+      for (int Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+	if (code->equiv != photcodes[Ns][0].code) continue;
+	found = TRUE;
+      }
+      if (!found) continue; 
+
+      // reset the Mflat value if requested
+      if (RESET_FLATCORR) {
+	catalog->measure[m].Mflat = 0.0;
+      }
+      
+      /* select measurements by time */
+      if (TimeSelect) {
+	if (catalog->measure[m].t < TSTART) continue;
+	if (catalog->measure[m].t > TSTOP) continue;
+      }
 	
-      off_t N = Nsecfilt*i+Ns;
-      dvo_secfilt_init (&catalog->secfilt[N], SECFILT_RESET_ALL);
+      // only reset Mcal for measures with a matching image
+      // do not reset Mcal for ubercal images unless explicitly requested
 
-      off_t m = catalog->average[i].measureOffset;
-
-      for (off_t j = 0; j < catalog->average[i].Nmeasure; j++, m++) {
-	
-	/* select measurements by time */
-	if (TimeSelect) {
-	  if (catalog->measure[m].t < TSTART) continue;
-	  if (catalog->measure[m].t > TSTOP) continue;
+      if (catalog->measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
+	if (!KEEP_UBERCAL) {
+	  catalog->measure[m].McalPSF  = 0.0;
+	  catalog->measure[m].McalAPER = 0.0;
+	  catalog->measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
+	} 
+      } else {
+	if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
+	  catalog->measure[m].McalPSF  = 0.0;
+	  catalog->measure[m].McalAPER = 0.0;
 	}
-	
-	// only reset Mcal for measures with a matching image
-	// do not reset Mcal for ubercal images unless explicitly requested
-
-	if (catalog->measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-	  if (!KEEP_UBERCAL) {
-	    catalog->measure[m].McalPSF  = 0.0;
-	    catalog->measure[m].McalAPER = 0.0;
-	    catalog->measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
-	  } 
-	} else {
-	  if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
-	    catalog->measure[m].McalPSF  = 0.0;
-	    catalog->measure[m].McalAPER = 0.0;
-	  }
-	}
-	catalog->measure[m].dbFlags &= ~photomBitsMeasure;
       }
+      catalog->measure[m].dbFlags &= ~photomBitsMeasure;
     }
   }
 }
 
+
Index: /branches/eam_branches/relphot.20210521/src/StarOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/StarOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/StarOps.c	(revision 41625)
@@ -1,3 +1,5 @@
 # include "relphot.h"
+
+# define STAR_BAD (ID_SECF_STAR_FEW | ID_SECF_STAR_POOR)
 
 enum {THREAD_RUN, THREAD_DONE};
@@ -290,11 +292,14 @@
 
 // setMrel and setMrelOutput are extremely similar, but have slightly different implications:
+// * both call setMrelCatalog, but setMrel set isSetMrelFinal to FALSE
 // * setMrel only uses the internal Tiny structures
-// * setMrel skips stars for which there are too few good measurements (setMrelOutput assigns an average mag with only 1)
+// * setMrel skips stars for which there are too few good measurements
+// -> (setMrelOutput assigns an average mag for all objects, even those with only 1 measurement)
 // * setMrelOutput sets psfQF, psfQFperf, extNsigma, rank info, etc
 // * setMrelOutput sets average Map
 // * setMrelOutput updates 2MASS average flags
 // * setMrelOutput updates average EXT flags (PS1 and 2MASS)
-// * setMrelOutput sets average stack & warp photometry (if requested)
+// * setMrelOutput sets average stack & warp photometry (if requested), setMrel always skips those
+// * setMrelOutput is NOT threaded it reads/writes catalogs one at a time
 
 int setMrel (Catalog *catalog, int Ncatalog) {
@@ -325,5 +330,4 @@
 }
 
-// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
 int setMrelOutput (Catalog *catalog, int Ncatalog) {
 
@@ -443,13 +447,4 @@
 }
 
-# define SKIP_THIS_MEAS(REASON) {				\
-    catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;	\
-    if (catalog[Nc].measure) {					\
-      catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;	\
-    }								\
-    results->REASON ++;						\
-    continue; }
-
-
 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
 
@@ -471,5 +466,5 @@
 }
 
-/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
+/* set measure.Mcal for all measures except ID_MEAS_NOCAL */
 int setMcalOutput (Catalog *catalog, int Ncatalog) {
 
@@ -477,7 +472,4 @@
   off_t j, k, m;
 
-  MEAS_BAD = ID_MEAS_NOCAL;
-  IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
@@ -485,14 +477,18 @@
       m = catalog[i].averageT[j].measureOffset;
       for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-	if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
+	if (catalog[i].measureT[m].dbFlags & ID_MEAS_NOCAL) continue;
 	float McalPSF  = getMcal  (m, i, MAG_CLASS_PSF);
 	float McalAPER = getMcal  (m, i, MAG_CLASS_APER);
 	if (isnan(McalPSF)) continue;
-	float Mmos  = getMmos  (m, i);
-	if (isnan(Mmos)) continue;
-	float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
-	if (isnan(Mgrp)) continue;
-	float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
-	if (isnan(Mgrid)) continue;
+
+	// Mmos, Mgrp have already been transfered to Mcal and should be zero here
+	// float Mmos  = getMmos  (m, i);
+	// if (isnan(Mmos)) continue;
+	// float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
+	// if (isnan(Mgrp)) continue;
+
+	// Mgrid is applied to Mflat and should be ignored here
+	// float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
+	// if (isnan(Mgrid)) continue;
 
 	// Note that this operation is setting measure->McalAPER to image->McalAPER only
@@ -507,8 +503,6 @@
 
 	// set the output calibration
-	// float Moff = Mmos + Mgrp + Mgrid;
-	float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
-	catalog[i].measure[m].McalPSF  = McalPSF + Moff;
-	catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
+	catalog[i].measure[m].McalPSF  = McalPSF;
+	catalog[i].measure[m].McalAPER = useStackAper ? McalAPER : McalPSF;
 
 	if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
@@ -612,5 +606,5 @@
     for (i = Ntot = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
-	if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
+	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
 	float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
 	if (isnan(Mchisq)) continue;
@@ -621,10 +615,14 @@
       }
     }
-  
+
+    // XXX the limits for MaxChisq and MaxScatter should be
+    // user-defined.  This uses 
     liststats (xlist, dlist, NULL, Ntot, &stats);
-    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
+    float ChiSqUpper90 = stats.Upper90;
+    MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90);
 
     liststats (slist, dlist, NULL, Ntot, &stats);
-    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
+    float ScatterUpper90 = stats.Upper90;
+    MaxScatter = MAX (STAR_SCATTER, ScatterUpper90);
     fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
 
@@ -636,5 +634,5 @@
 	mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
 	if (mark) {
-	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_OBJ_POOR;
+	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR;
 	  Ndel ++;
 	  if (dM > MaxScatter)   { Nscat ++; }
@@ -642,5 +640,5 @@
 	  if (Mchisq > MaxChisq) { Nchi ++; }
 	} else {
-	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_POOR;
+	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
 	}
 	Nave ++;
@@ -652,181 +650,4 @@
   free (slist);
   free (dlist);
-}
-
-// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
-// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
-// sigma of the INNER 50% mean.  it then flags any measurements which are more than
-// NSIGMA_REJECT (5) sigma of the mean
-
-// this function only operations on the PSF magnitudes
-
-// XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers
-
-# define NSIGMA_CLIP 3.0
-# define NSIGMA_REJECT 5.0
-void clean_measures (Catalog *catalog, int Ncatalog, int final) {
-
-  off_t j, k, m, Nmax, Ndel, Nave;
-  int i, N, image_bad, TOOFEW;
-  off_t *ilist;
-  double *tlist, *list, *dlist;
-  float Msys, Mcal, Mmos, Mgrp, Mgrid;
-  int Ncal, Nmos, Ngrp, Ngrid, Nfew;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
-  /* Nmeasure is now different, need to reallocate */
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
-    }
-  }
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
-  ALLOCATE (ilist, off_t, Nmax);
-  ALLOCATE (tlist, double, Nmax);
-  
-  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
-  TOOFEW = MAX (5, STAR_TOOFEW);
-
-  // stats structures for inner and full stats
-  StatType instats, stats;
-  liststats_setmode (&instats, "INNER_MEAN");
-  liststats_setmode (&stats, "MEAN");
-
-  Ndel = Nave = 0;
-  Ncal = Nmos = Ngrp = Ngrid = Nfew = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-	
-	int thisCode = photcodes[Ns][0].code;
-	int Nsec = GetPhotcodeNsec(thisCode);
-	
-	/* on final processing, skip stars already measured */
-	if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue;  
-
-	/* skip bad stars to prevent them from becoming good (on inner sample) */
-	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-
-	/* accumulate list of valid measurements */
-	m = catalog[i].averageT[j].measureOffset;
-	N = 0;
-	for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-	  // skip measurements that do not match the current photcode
-	  int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-	  if (ecode != thisCode) { continue; }
-
-	  // SKIP gpc1 stack data
-	  if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
-	  
-	  // SKIP gpc1 forced-warp data
-	  if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
-	  
-	  // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-	  Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
-	  if (isnan(Mcal)) { Ncal ++; continue; }
-	  Mmos  = getMmos  (m, i);
-	  if (isnan(Mmos)) { Nmos ++; continue; }
-	  Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
-	  if (isnan(Mgrp)) { Ngrp ++; continue; }
-	  Mgrid = getMgridTiny (&catalog[i].measureT[m]);
-	  if (isnan(Mgrid)) { Ngrid ++; continue; }
-
-	  Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-	  list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
-	  dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-	  N++;
-	}
-	if (N <= TOOFEW) { 
-	  Nfew ++; 
-	  continue; 
-	}
-
-	/* 3-sigma clip based on stats of inner 50% */
-
-	// calculated mean of inner 50%
-	liststats (list, dlist, NULL, N, &instats);
-	instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
-
-	// ignore entries > 3sigma from inner mean
-	for (k = m = 0; k < N; k++) {
-	  if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
-	    list[m] = list[k];
-	    m++;
-	  }
-	}
-	// recalculate the mean & sigma of the accepted measurements
-	liststats (list, dlist, NULL, m, &stats);
-	stats.sigma = MAX (MIN_ERROR, stats.sigma);
-
-	/* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
-	image_bad = IMAGE_BAD;
-	IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-	m = catalog[i].averageT[j].measureOffset;
-	N = 0;
-	for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-	  // skip measurements that do not match the current photcode
-	  int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-	  if (ecode != thisCode) { continue; }
-
-	  // SKIP gpc1 stack data
-	  if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
-	  
-	  // SKIP gpc1 forced-warp data
-	  if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
-	  
-	  // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-	  Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
-	  if (isnan(Mcal)) continue;
-	  Mmos  = getMmos  (m, i);
-	  if (isnan(Mmos)) continue;
-	  Mgrid = getMgridTiny (&catalog[i].measureT[m]);
-	  if (isnan(Mgrid)) continue;
-
-	  Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-	  list[N] = Msys - Mcal - Mmos - Mgrid;
-	  dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-	  ilist[N] = m;
-	  N++;
-	  Nave ++;
-	}
-	if (N < TOOFEW) continue;
-
-	/* mark bad measures (> 5 sigma deviant) */
-	for (k = 0; k < N; k++) {
-	  // treat the scatter of the star as a systematic term.  this is a bit of an
-	  // over-estimage (a perfect Gauss distribution with perfect errors would have
-	  // mySigma = sqrt(2) too large)
-	  float mySigma = hypot (stats.sigma, dlist[k]);
-	  if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
-	    catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-	    if (final) {
-	      // for the final pass, we have a duplicate set of values in measure and measureT
-	      catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-	    }
-	    Ndel ++;
-	  } else {
-	    catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-	    if (final) {
-	      // for the final pass, we have a duplicate set of values in measure and measureT
-	      catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-	    }
-	  }
-	}
-	IMAGE_BAD = image_bad;
-      }
-    }
-  }
-  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
-
-  free (list);
-  free (dlist);
-  free (ilist);
-  free (tlist);
 }
 
Index: /branches/eam_branches/relphot.20210521/src/TGroupOps.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/TGroupOps.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/TGroupOps.c	(revision 41625)
@@ -1172,11 +1172,9 @@
   ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
 
-  int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW;
-
   int n = 0;
   for (int i = 0; i < NtgroupTimes; i++) {
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & NIGHT_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].McalPSF;
       dlist[n] = 1;
@@ -1209,5 +1207,5 @@
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & IMAGE_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].dMcal;
       dlist[n] = 1;
@@ -1240,5 +1238,5 @@
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & IMAGE_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].McalChiSq;
       dlist[n] = 1;
@@ -1396,5 +1394,5 @@
 	c = tgroup[i].catalog[j];
       
-	if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
+	if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
 	// ave = catalog[c].measureT[m].averef;
Index: /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/ZeroPointModes.c	(revision 41625)
@@ -22,5 +22,5 @@
 }
 
-void SetZeroPointModes (void) {
+void SetZeroPointModes (Catalog *catalog, int Ncatalog) {
 
   if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
@@ -59,4 +59,5 @@
       clean_mosaics(); // do this on every pass 
       clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
     }
     if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) {
@@ -71,4 +72,5 @@
       clean_mosaics(); // do this on every pass 
       clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
     }
   }
Index: /branches/eam_branches/relphot.20210521/src/initialize.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/initialize.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/initialize.c	(revision 41625)
@@ -27,8 +27,4 @@
   }    
     
-  IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-  STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-  MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-
   if (SHOW_PARAMS) {
     int Ns;
Index: /branches/eam_branches/relphot.20210521/src/plot_scatter.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/plot_scatter.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/plot_scatter.c	(revision 41625)
@@ -28,32 +28,32 @@
     N = 0;
     for (i = 0; i < Ncatalog; i++) {
-	for (j = 0; j < catalog[i].Naverage; j++) {
+      for (j = 0; j < catalog[i].Naverage; j++) {
 
-	    /* calculate the average value for a single star */
-	    if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-	    m = catalog[i].averageT[j].measureOffset;
+	/* calculate the average value for a single star */
+	// if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
+	m = catalog[i].averageT[j].measureOffset;
 
-	    for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-		// skip measurements that do not match the current photcode
-		int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-		if (ecode != thisCode) { continue; }
+	for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
+	  // skip measurements that do not match the current photcode
+	  int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
+	  if (ecode != thisCode) { continue; }
 
-		if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-		Mcal = getMcal  (m, i, MAG_CLASS_PSF);
-		if (isnan(Mcal)) continue;
-		Mmos = getMmos  (m, i);
-		if (isnan(Mmos)) continue;
-		Mgrid = getMgridTiny (&catalog[i].measureT[m]);
-		if (isnan(Mgrid)) continue;
+	  if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
+	  Mcal = getMcal  (m, i, MAG_CLASS_PSF);
+	  if (isnan(Mcal)) continue;
+	  Mmos = getMmos  (m, i);
+	  if (isnan(Mmos)) continue;
+	  Mgrid = getMgridTiny (&catalog[i].measureT[m]);
+	  if (isnan(Mgrid)) continue;
 
-		Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
-		if (isnan(Mrel)) continue;
+	  Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
+	  if (isnan(Mrel)) continue;
 
-		xlist[N] = Mrel;
-		ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
-		ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
-		N++;
-	    }
+	  xlist[N] = Mrel;
+	  ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
+	  ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
+	  N++;
 	}
+      }
     }
 
Index: /branches/eam_branches/relphot.20210521/src/relphot_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/relphot_images.c	(revision 41625)
@@ -83,4 +83,5 @@
     SAVEPLOT = FALSE;
 
+    // mark ID_MEAS_AREA (user-selected chip region) and ID_MEAS_NOCAL (time range & active photcode)
     setExclusions (catalog, Ncatalog, TRUE);
 
@@ -114,5 +115,5 @@
 
       global_stats (catalog, Ncatalog, i);
-      SetZeroPointModes ();
+      SetZeroPointModes (catalog, Ncatalog);
 
       MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
Index: /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/relphot_parallel_images.c	(revision 41625)
@@ -8,4 +8,6 @@
  * update the unowned detections for owned objects from neighbor regions
  * update the unowned objects for owned detections 
+
+ * TGROUP and GRID corrections are not calculated.
  */
 
@@ -105,5 +107,4 @@
     setMcal  (catalog);
     setMmos  (catalog);
-    // setMgrp  (catalog, FALSE); XXX think this through: this may not make sense for tgroups
     MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
     
@@ -117,5 +118,5 @@
     global_stats (catalog, Ncatalog, i);
 
-    SetZeroPointModes ();
+    SetZeroPointModes (catalog, Ncatalog);
     MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
   }
@@ -138,5 +139,4 @@
 
   setMcalFromMosaics (); // copy per-mosaic calibrations to the images
-  // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images
 
   share_image_mags (regionHosts, -1);
@@ -152,5 +152,4 @@
   freeImageBins(Ncatalog, TRUE);
   freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins ();
   freeImages((char *)image);
   free (image);
Index: /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/setMrelCatalog.c	(revision 41625)
@@ -225,4 +225,18 @@
     }
 
+    /* some things to note here:
+
+       Mcal, Mmos, Mgrid are only relevant (non-zero) during the relphot_images and
+       relphot_parallel_images analysis steps. For the final output steps (setMrelFinal),
+       these are turned off
+
+       During the image analysis, Mflat is the prior stored value while Mgrid is the new
+       measurement.  In the final output step, Mgrid is transferred to Mflat and zero
+
+       On the final calculation of Mrel, the value of Mgrid has been applied to Mflat, but
+       it has not been removed from the GridOps.c structures.  In order to avoid
+       double-counting, we need to skip Mgrid on the final calculation.
+    */
+
     // ** Choose the calibration (depends on the mode : do I have an image reference or not?) 
     float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
@@ -235,18 +249,21 @@
       Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value
     } else {
-      // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat
+      // getMcal returns image[].Mcal
       Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
       if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
-      Mmos  = getMmos  (meas, cat);
-      if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
-      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
-      if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
-      Mgrid = isSetMrelFinal ? 0.0 : getMgridTiny (&measureT[k]); 
-      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+
+      // the flat-field correction is stored in measure.Mflat
       Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
-    }
-    // on the final calculation of Mrel the value of Mgrid has been applied to Mflat, but
-    // it has not been removed from the GridOps.c structures.  In order to avoid double-counting,
-    // we need to skip Mgrid on the final calculation.
+
+      // see note above re: final output vs image analysis
+      if (!isSetMrelFinal) {
+	Mmos  = getMmos  (meas, cat);
+	if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+	Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
+	if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
+	Mgrid = getMgridTiny (&measureT[k]); 
+	if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+      }
+    }
 
     int myUbercalDist = getUbercalDist(meas, cat);
@@ -367,6 +384,7 @@
     }
 
+    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
+    // NOTE: This test only applies if this bit was set by an earlier analysis
     if (PRESERVE_PS1 && (secfilt[Nsec].flags & ID_SECF_HAS_PS1)) continue;
-    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
 
     // XXX hardwired grizy = (01234) JHK = (567) w = (8)
@@ -419,9 +437,10 @@
     int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
     if (NrankingPSF < Nminmeas) { 
-      secfilt[Nsec].flags |= ID_OBJ_FEW;
+      secfilt[Nsec].flags |= ID_SECF_STAR_FEW;
     } else {
       secfilt[Nsec].MpsfChp  = psfstats->mean;
       secfilt[Nsec].dMpsfChp = psfstats->error;
       secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
+      secfilt[Nsec].flags &= ~ID_SECF_STAR_FEW;
     }
 
@@ -629,6 +648,4 @@
       haveStackObject = TRUE;
 
-      // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
-
       int isPrimary = FALSE;
       
Index: /branches/eam_branches/relphot.20210521/src/setMrelFinal.c
===================================================================
--- /branches/eam_branches/relphot.20210521/src/setMrelFinal.c	(revision 41624)
+++ /branches/eam_branches/relphot.20210521/src/setMrelFinal.c	(revision 41625)
@@ -13,7 +13,7 @@
   ResetAverageAndMeasure (catalog);
   
-  // this sets flags in the measureT element, not the measure element
+  // This sets ID_MEAS flags in the measureT element, not the measure element
+  // These bits are thus ephemeral and not saved.  
   setExclusions (catalog, 1, VERBOSE);  /* mark by area */
-  // XXXXX THis probably does not make sense
 
   /* set catalog[0].found[i] = FALSE */
@@ -31,37 +31,5 @@
   setMrelOutput (catalog, 1); // sets the values secfilt.MpsfChp = <measure.M - image.Mcal - measure.Mflat>
   setMcalOutput (catalog, 1); // sets measure.Mcal = image.Mcal
-
-  // clear ID_OBJ_POOR, ID_OBJ_FEW, ID_MEAS_NOCAL values before writing 
-  // These are internal bits, so they should be cleared
-  // XXX are these still used?
-  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
-  for (i = 0; i < catalog[0].Naverage; i++) {
-    catalog[0].average[i].flags &= ~ID_OBJ_FEW;
-    catalog[0].average[i].flags &= ~ID_OBJ_POOR;
-    off_t j;
-    for (j = 0; j < Nsecfilt; j++) {
-	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_FEW; 
-	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_POOR; 
-    }
-    off_t m = catalog[0].average[i].measureOffset;
-    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
-      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
-    }
-  }
 }
-
-/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
-// For each star & photcode, we are selecting the measurements to use in the average
-// photometry value.  This function is called several times, lowering the bar to
-// acceptance on each pass. The rejections at the different passes are:
-
-// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
-// 1 : accept measurements thought to be POOR based on photflags
-// 2 : accept the measurements marked as outliers; accept images marked as outliers
-// 3 : accept measurements thought to be BAD based on photflags
-// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)
-
-// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
-// non-outliers with Imag range
 
 // int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
@@ -84,129 +52,11 @@
 // }
 
-int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
-
-// XXX this function is no longer used:
-void skip_measurements (Catalog *catalog, int pass) {
-
-  off_t i, k, m;
-  off_t Ntot, Ntry, Nkeep, Nskip;
-  float mag;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  Ntot = Ntry = Nskip = Nkeep = 0;
-
-  /* allow measures from images marked POOR and FEW */
-  if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-  
-  /* allow measures marked as outliers (POOR) and off image region (AREA) */
-  if (pass >= 3) {
-    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
-  } else {
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-  }
-
-  /* mark measures which should be ignored on this pass */
-  for (i = 0; i < catalog[0].Naverage; i++) {
-    Ntot += catalog[0].average[i].Nmeasure;
-
-    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
-      fprintf (stderr, "test obj\n");
-      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
-    }
-
-    // mark measurements for each secfilt separately
-    int Ns;
-    for (Ns = 0; Ns < Nphotcodes; Ns++) {
-      
-      int thisCode = photcodes[Ns][0].code;
-      int Nsec = GetPhotcodeNsec(thisCode);
-      
-      /* star/photcodes already calibrated */
-      if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue;  
-
-      m = catalog[0].average[i].measureOffset;
-      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
-	Ntry++;
-
-	// skip measurements not related to this photcode
-	PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
-	if (!code) continue;
-	if (code->equiv != thisCode) continue;
-
-	// skip measurements by time range (mark as skipped for this secfilt)
-	// XXX note that this is a bit dangerous : some objects may never get calibrated 
-	if (TimeSelect) {
-	  if (catalog[0].measure[m].t < TSTART) goto skip;
-	  if (catalog[0].measure[m].t > TSTOP) goto skip;
-	}
-	
-	// we now have a measurement relavant to this photcode (and time range)
-
-	// clear SKIP for all measures at first 
-	catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-	catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-	// skip measurements from BAD images and mosaics (not REF mags)
-	// do NOT skip measurements without a matching image (REF mags)
-	off_t Nim = getImageEntry (m, 0);
-	if (Nim > -1) {
-	  if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip;
-	  if (isnan(getMmos (m, 0))) goto skip;
-	}
-
-	// PASS 4 : skip measurements by inst mag limit (not REF mags)
-	if ((pass < 4) && ImagSelect) {
-	  if (Nim > -1) {
-	    mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
-	    if (mag < ImagMin) goto skip;
-	    if (mag > ImagMax) goto skip;
-	  }
-	}
-
-	// PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers
-	if (pass < 3) {
-	  if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
-	  if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
-	    if (catalog[0].measure[m].psfQF < 0.85) goto skip;
-	  }
-	}
-	
-	// PASS 2 : (no additional cut)
-
-	// PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
-	if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;
-
-	// PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)
-
-	Nkeep ++;
-	
-	continue;
-
-      skip:
-	catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
-	catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
-	Nskip ++;
-      }
-    }
-  }
-  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
-			pass, Ntot, Ntry, Nskip, Nkeep);
-}
-
-// XXX review this, especially the good/bad images
-// by this point, we have propagated the mosaic and tgroup flags to each image
+// This function is only called for the final output step.  By this point, we have
+// propagated the mosaic and tgroup flags to each image.
 void setMeasureRank (Catalog *catalog) {
 
   int i;
 
-  // make these global
-  int IMAGE_BAD  = ID_IMAGE_PHOTOM_NOCAL;
-  int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-  int MEAS_BAD   = ID_MEAS_NOCAL        | ID_MEAS_SKIP_PHOTOM;
-  int MEAS_POOR  = ID_MEAS_POOR_PHOTOM  | ID_MEAS_AREA;
-
   Measure     *measure     = catalog[0].measure;
-  MeasureTiny *measureT    = catalog[0].measureT;
   char        *measureRank = catalog[0].measureRank;
 
@@ -215,9 +65,5 @@
     measureRank[i] = 11; // start at a low rank
 
-    // clear SKIP for all measures at first 
-    measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-    measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-    // skip measurements without a valid photcode
+    // measurements without a valid photcode have lowest rank (should not be used anyway)
     PhotCode *code = GetPhotcodebyCode (measure[i].photcode);
     if (!code) continue;
@@ -225,49 +71,50 @@
     // measurements outside time range have poor rank
     if (TimeSelect) {
-      if (measure[i].t < TSTART) 		     { measureRank[i] = 10; continue; }
-      if (measure[i].t > TSTOP)  		     { measureRank[i] = 10; continue; }
+      if (measure[i].t < TSTART) 		                        { measureRank[i] = 10; continue; }
+      if (measure[i].t > TSTOP)  		                        { measureRank[i] = 10; continue; }
     }
     
-    // measurements from BAD images and mosaics (not REF mags)
-    off_t Nim = getImageEntry (i, 0);
+    // Nim < 0 for REF mags, imageFlags have bits for IMAGE, MOSAIC, NIGHT
+    off_t Nim      = getImageEntry (i, 0);
     int imageFlags = getImageFlags (i, 0);
-    int mosaicFlags = getMosaicFlags (i, 0);
 
     if (Nim > -1) {
-      if (imageFlags  & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
-
-      if (mosaicFlags & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
-    
       // measurements ranked by inst mag limit (not REF mags)
       if (ImagSelect) {
 	float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
-	if (mag < ImagMin) 			     { measureRank[i] = 8; continue; }
-	if (mag > ImagMax) 			     { measureRank[i] = 8; continue; }
+	if (mag < ImagMin) 			                        { measureRank[i] = 9; continue; }
+	if (mag > ImagMax) 			                        { measureRank[i] = 9; continue; }
       }
     }
 
-    // RANK 1 : BAD photFlags
-    if (measure[i].dbFlags & MEAS_BAD)               { measureRank[i] = 7; continue; }
-
-    // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers
-    if (measure[i].photFlags & code->photomBadMask)  { measureRank[i] = 6; continue; }
-
+    // RANK 8 : Poor image
     if (Nim > -1) {
-      if (imageFlags  & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
-      if (mosaicFlags & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
+      if (imageFlags & ID_IMAGE_PHOTOM_POOR)                            { measureRank[i] = 8; continue; }
     }
 
-    // RANK 2 : psfQF value
-    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; }
+    // RANK 7 : BAD photFlags (eg, SAT, CR), internal outliers
+    if (measure[i].photFlags & code->photomBadMask)                     { measureRank[i] = 7; continue; }
+
+    // RANK 6 : bad psfQF value
+    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85)         { measureRank[i] = 6; continue; }
 	
-    // RANK 1 : POOR photFlags
-    if (measure[i].dbFlags & MEAS_POOR)              { measureRank[i] = 3; continue; }
+    // RANK 5 : not in valid chip region 
+    if (measure[i].dbFlags & ID_MEAS_AREA)                              { measureRank[i] = 5; continue; }
 
-    // RANK 1 : POOR photFlags
-    if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; }
+    // RANK 4 : POOR photFlags
+    if (measure[i].photFlags & code->photomPoorMask)                    { measureRank[i] = 4; continue; }
     
-    // RANK 1 : psfQFperfect value
-    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; }
+    // RANK 3 : bad psfQFperfect value
+    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 3; continue; }
 	
+    // RANK 2 : Poor mosaic
+    if (Nim > -1) {
+      if (imageFlags & ID_IMAGE_MOSAIC_POOR)                            { measureRank[i] = 2; continue; }
+    }
+
+    // RANK 1 : Poor night
+    if (Nim > -1) {
+      if (imageFlags & ID_IMAGE_NIGHT_POOR)                             { measureRank[i] = 1; continue; }
+    }
     // RANK 0 : perfect measurement:
     measureRank[i] = 0;