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);
 }