Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h	(revision 33886)
@@ -36,4 +36,15 @@
 } Mosaic; 
 
+typedef enum {
+  STATS_NONE, 
+  STATS_MEAN, 
+  STATS_MEDIAN, 
+  STATS_WT_MEAN, 
+  STATS_INNER_MEAN, 
+  STATS_INNER_WTMEAN, 
+  STATS_CHI_INNER_MEAN, 
+  STATS_CHI_INNER_WTMEAN
+} ListStatsMode;
+
 typedef struct {
   double median;
@@ -48,4 +59,5 @@
   double total;
   int    Nmeas;
+  ListStatsMode statmode;
 } StatType;
 
@@ -245,6 +257,6 @@
 void          initialize          PROTO((int argc, char **argv));
 void          initialize_client   PROTO((int argc, char **argv));
-void          initstats           PROTO((char *mode));
-int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
+void          liststats_setmode   PROTO((StatType *stats, char *strmode));
+int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath));
 Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog));
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/GridOps.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/GridOps.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/GridOps.c	(revision 33886)
@@ -527,5 +527,7 @@
   double *list, *dlist;
   float Msys, Mrel, Mcal, Mmos;
+
   StatType stats;
+  liststats_setmode (&stats, STATMODE);
   
   if (!USE_GRID) return;
@@ -595,5 +597,5 @@
     }
 
-    liststats (list, dlist, N, &stats);
+    liststats (list, dlist, NULL, N, &stats);
     gridM[i] = stats.mean;
     gridS[i] = stats.sigma;
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c	(revision 33886)
@@ -399,5 +399,7 @@
   float Msys, Mrel, Mmos, Mgrid, Mflat;
   double *list, *dlist, *Mlist, *dMlist;
+
   StatType stats;
+  liststats_setmode (&stats, STATMODE);
 
   if (FREEZE_IMAGES) return;
@@ -518,5 +520,6 @@
     }
 
-    liststats (list, dlist, N, &stats);
+    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
+    liststats (list, dlist, NULL, N, &stats);
     image[i].Mcal  = stats.mean;
     image[i].dMcal = stats.error;
@@ -527,5 +530,5 @@
 
     // bright end scatter
-    liststats (Mlist, dMlist, Nbright, &stats);
+    liststats (Mlist, dMlist, NULL, Nbright, &stats);
     image[i].dMagSys = stats.sigma;
 
@@ -560,5 +563,4 @@
   double *mlist, *slist, *dlist;
   double MaxOffset, MaxScatter, MedOffset;
-  StatType stats;
 
   if (FREEZE_IMAGES) return;
@@ -578,9 +580,14 @@
     N++;
   }
-  initstats ("MEAN");
-  liststats (mlist, dlist, N, &stats);
+
+  // use a straight mean to find the global image statistics (no weighting)
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, dlist, NULL, N, &stats);
   MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
   MedOffset = stats.median;
-  liststats (slist, dlist, N, &stats);
+
+  liststats (slist, dlist, NULL, N, &stats);
   MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
@@ -603,5 +610,4 @@
 
   fprintf (stderr, "%d images marked poor\n", Nmark);
-  initstats (STATMODE);
   free (mlist);
   free (slist);
@@ -724,6 +730,6 @@
   double *list, *dlist;
   float Mcal, Mmos, Mgrid;
+
   StatType stats;
-
   bzero (&stats, sizeof (StatType));
   if (FREEZE_IMAGES) return (stats);
@@ -755,5 +761,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -783,5 +791,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -811,5 +821,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -839,5 +851,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c	(revision 33886)
@@ -764,5 +764,7 @@
   Image *imageReal;
   off_t j, NimageReal;
+
   StatType stats;
+  liststats_setmode (&stats, STATMODE);
 
   double *list   = info->list;
@@ -930,5 +932,6 @@
     }
   }
-  liststats (list, dlist, N, &stats);
+
+  liststats (list, dlist, NULL, N, &stats);
   if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
@@ -945,5 +948,5 @@
 
   // bright end scatter
-  liststats (Mlist, dMlist, Nbright, &stats);
+  liststats (Mlist, dMlist, NULL, Nbright, &stats);
   myMosaic[0].dMsys = stats.sigma;
 
@@ -1310,5 +1313,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1338,5 +1343,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1383,5 +1390,7 @@
   // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1411,5 +1420,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1423,5 +1434,4 @@
   double *mlist, *slist, *dlist;
   double MaxOffset, MedOffset, MaxScatter;
-  StatType stats;
 
   if (!MOSAIC_ZEROPT) return;
@@ -1442,9 +1452,13 @@
     N++;
   }
-  initstats ("MEAN");
-  liststats (mlist, dlist, N, &stats);
+
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, dlist, NULL, N, &stats);
   MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
   MedOffset = stats.median;
-  liststats (slist, dlist, N, &stats);
+
+  liststats (slist, dlist, NULL, N, &stats);
   MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
@@ -1476,5 +1490,5 @@
 
   fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
-  initstats (STATMODE);
+
   free (mlist);
   free (slist);
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c	(revision 33886)
@@ -13,4 +13,5 @@
   double *list;
   double *dlist;
+  double *wlist;
   double *aplist;
   double *daplist;
@@ -79,4 +80,5 @@
     ALLOCATE (results->list,  double, Nmax);
     ALLOCATE (results->dlist, double, Nmax);
+    ALLOCATE (results->wlist, double, Nmax);
   }
 }
@@ -85,4 +87,5 @@
   free (results->list);
   free (results->dlist);
+  free (results->wlist);
 }
 
@@ -159,5 +162,5 @@
   SetMrelInfo summary, results;
   SetMrelInfoInit (&summary, FALSE);
-  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist
+  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
   ALLOCATE (results.aplist, double, Nmax);
   ALLOCATE (results.daplist, double, Nmax);
@@ -303,12 +306,16 @@
   int N;
   float Msys, Mcal, Mmos, Mgrid;
+
   StatType stats, apstats;
+  liststats_setmode (&stats, STATMODE);
+  liststats_setmode (&apstats, STATMODE);
 
   double *list    = results->list;
   double *dlist   = results->dlist;
+  double *wlist   = results->wlist;
   double *aplist  = results->aplist;
   double *daplist = results->daplist;
 
-  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays
+  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
 
   int isSetMrelFinal = (pass >= 0);
@@ -390,6 +397,7 @@
 	  aplist[N] = Map - Mcal - Mmos - Mgrid;
 
-	  // count the extended detections
+	  // special options for PS1 data
 	  if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
+	    // count the extended detections
 	    if (!isnan(catalog[Nc].measure[m].Map)) {
 	      if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
@@ -427,6 +435,6 @@
 	}
 
-	// dlist gives the error, which is used as the weight in WT_MEAN.
-	// we can modify the resulting weight in a few ways:
+	// dlist gives the error per measurement, wlist gives the weight
+	// we can modify the error and weight in a few ways:
 	// 1) MIN_ERROR guarantees a floor
 	// 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
@@ -434,8 +442,9 @@
 	// 4) some reference photcode of some kind can be specified as fixed and have a high weight
 	dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
+	wlist[N] = 1.0;
 
 	// up-weight the ubercal values (or convergence can take a long time...)
 	if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-	  dlist[N] = MAX (0.1*catalog[Nc].measureT[m].dM, MIN_ERROR);
+	  wlist[N] = 10.0;
 	}
 
@@ -445,7 +454,5 @@
 	if (refPhotcode) {
 	  if (code->code == refPhotcode->code) {
-	    // tiny error -> large weight
-	    // dlist[N] = MAX (0.01*catalog[Nc].measureT[m].dM, MIN_ERROR);
-	    dlist[N] = 0.0001;
+	    wlist[N] = 100.0;
 	  }
 	}
@@ -469,5 +476,5 @@
       }	
 
-      liststats (list, dlist, N, &stats);
+      liststats (list, dlist, wlist, N, &stats);
 
       catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
@@ -481,5 +488,4 @@
 
       if (isSetMrelFinal) {
-	liststats (aplist, daplist, N, &apstats);
 	catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE;
 
@@ -490,4 +496,9 @@
 	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
 	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
+
+	// NOTE : use the modified weight for apmags as well as psf mags
+	liststats (aplist, daplist, wlist, N, &apstats);
+
+	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean; 
 
 	// NOTE: for 2MASS measurements, Next should be 1, as should N
@@ -567,197 +578,4 @@
 }
 
-# if (0)
-
-/* grab Nsec for named photcode */
-# define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME)	\
-  short MY_NSEC = -1;				\
-  {						\
-  PhotCode *code = GetPhotcodebyName (NAME);	\
-  if (code) {					\
-    MY_NSEC = GetPhotcodeNsec (code->equiv);	\
-  } }
-
-// For each average object, set the average mags based on existing equiv photometry.
-// NOTE: this function operates on the real Measure & Average structures, not the
-// MeasureTiny & AverageTiny structures
-// NOTE: this function is called on the remote machine -- make sure any (global) options 
-// are passed to the relphot_client program
-int setMave (Catalog *catalog, int Ncatalog) {
-
-  off_t j, k, m, Nmax;
-  int i, Ns, Nsecfilt, N, Nc;
-  float Msys;
-  double *list, *dlist;
-  StatType stats;
-  DVOAverageFlags flagBits;
-
-  flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
-
-  // pre-allocate a list for stats purposes
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
-    }
-  }
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
-
-  Nsecfilt = GetPhotcodeNsecfilt ();
-
-  // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data
-  NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J");
-  NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H");
-  NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K");
-
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      // update average photometry for each of the average filters
-
-      // XXX Note that this would be faster if we had an array of results and accumulated
-      // them in a single pass
-
-      int Next = 0;
-
-      for (Ns = 0; Ns < Nsecfilt; Ns++) {
-
-	PhotCode *code = GetPhotcodebyNsec (Ns);
-	Nc = code[0].code;
-	
-	N = 0;
-	m = catalog[i].average[j].measureOffset;
-	for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-	  if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-	  if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
-
-	  Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-	  if (isnan(Msys)) continue;
-
-	  // reject POOR detections (PSF_QF < 0.85) and count extended detections
-	  // XXX only apply this filter for psphot data from GPC1 for now...
-	  if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
-	    if (catalog[i].measure[m].psfQual < 0.85) continue;
-	    if (!isnan(catalog[i].measure[m].Map)) {
-	      if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
-		Next ++;
-	      }
-	    }
-	  }
-
-	  // XXX make it optional to apply Mcal (do I need an extra field for advisory Mcal?)
-	  list[N] = Msys - catalog[i].measure[m].Mcal;
-	  dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-	  N++;
-	}
-	if (N < 1) continue;
-
-	liststats (list, dlist, N, &stats);
-
-	/* use sigma or error in dM for output? */
-	catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
-	catalog[i].secfilt[Nsecfilt*j+Ns].dM = stats.error;
-	catalog[i].secfilt[Nsecfilt*j+Ns].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
-	catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-	catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
-	catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
-
-	catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80;
-	catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20;
-
-	if ((Next > 0) && (Next > 0.5*N)) {
-	  catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT;
-	}
-      }
-
-      // update average flags based on the detection stats.  
-
-      int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
-      int goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
-      int good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
-      int nEXT = 0;
-      int nPSF = 0;     // number of PS1 PSF vs EXT measurements
-      int have2MASS = FALSE;
-      
-      // count, flag good and extended detections
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-
-	// PS1 data :
-	if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
-	  if (catalog[i].measure[m].psfQual > 0.85) {
-	    goodPS1 = TRUE;
-	    if (!isnan(catalog[i].measure[m].Map)) {
-	      if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
-		nEXT ++;
-	      } else {
-		nPSF ++;
-	      }
-	    }
-	  }
-	}
-	  
-	// 2MASS data J-band flags
-	if (catalog[i].measure[m].photcode == 2011) {
-	  // only need to do this once (always have JHK triplet; galaxy flag is same for all 3)
-	  have2MASS = TRUE;
-	  if (catalog[i].measure[m].photFlags & 0x00c00000) {
-	    Galaxy2MASS = TRUE;
-	  }
-	  if (catalog[i].measure[m].photFlags & 0x00000007) {
-	    good2MASS = TRUE; 
-	    if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0;
-	  } else {
-	    if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1;
-	  }
-	}  
-	// 2MASS data H-band flags
-	if (catalog[i].measure[m].photcode == 2012) {
-	  if (catalog[i].measure[m].photFlags & 0x00000007) {
-	    good2MASS = TRUE; 
-	    if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0;
-	  } else {
-	    if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1;
-	  }
-	}  
-	// 2MASS data K-band flags
-	if (catalog[i].measure[m].photcode == 2013) {
-	  if (catalog[i].measure[m].photFlags & 0x00000007) {
-	    good2MASS = TRUE; 
-	    if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0;
-	  } else {
-	    if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1;
-	  }
-	}  
-      }
-
-      // we attempt to set a few flags here; reset those bits before trying:
-      catalog[i].average[j].flags &= ~flagBits;
-
-      // XXX set the secfilt bits?
-      if (nEXT && (nEXT > nPSF)) {
-	catalog[i].average[j].flags |= ID_OBJ_EXT;
-      }
-      if (goodPS1) {
-	catalog[i].average[j].flags |= ID_OBJ_GOOD;
-      }
-      if (Galaxy2MASS) {
-	catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
-	if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT;
-	if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT;
-	if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT;
-      }
-      if (good2MASS) {
-	catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
-      }
-    }
-  }
-
-  free (list);
-  free (dlist);
-  return (TRUE);
-}
-# endif
-
 /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
 int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
@@ -804,5 +622,7 @@
   double Chisq, MaxScatter, MaxChisq;
   double *xlist, *slist, *dlist;
+
   StatType stats;
+  liststats_setmode (&stats, "MEAN");
 
   if (VERBOSE) fprintf (stderr, "marking poor stars\n");
@@ -818,11 +638,5 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  // XX int oldPLOTSTUFF = PLOTSTUFF;
-  // XX PLOTSTUFF = TRUE;
-  // XX plot_chisq (catalog, Ncatalog);
-  // XX PLOTSTUFF = oldPLOTSTUFF;
-
   // eliminate bad stars using the stats for a single secfilt at a time
-  // XXX DEP replace average.flags with secfilt flags
   for (Ns = 0; Ns < Nphotcodes; Ns ++) {
     
@@ -843,8 +657,8 @@
     }
   
-    initstats ("MEAN");
-    liststats (xlist, dlist, Ntot, &stats);
+    liststats (xlist, dlist, NULL, Ntot, &stats);
     MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
-    liststats (slist, dlist, Ntot, &stats);
+
+    liststats (slist, dlist, NULL, Ntot, &stats);
     MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
     fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
@@ -870,5 +684,4 @@
     }
     fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
-    initstats (STATMODE);
   }
   free (xlist);
@@ -876,4 +689,9 @@
   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
 
 # define NSIGMA_CLIP 3.0
@@ -886,5 +704,4 @@
   double *tlist, *list, *dlist;
   float Msys, Mcal, Mmos, Mgrid;
-  StatType stats;
   int Ncal, Nmos, Ngrid, Nfew;
 
@@ -906,4 +723,9 @@
   /* 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;
@@ -953,11 +775,10 @@
 
 	// calculated mean of inner 50%
-	initstats ("INNER_MEAN");
-	liststats (list, dlist, N, &stats);
-	stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
+	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] - stats.median) < NSIGMA_CLIP*stats.sigma) {
+	  if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
 	    list[m] = list[k];
 	    m++;
@@ -965,6 +786,5 @@
 	}
 	// recalculate the mean & sigma of the accepted measurements
-	initstats ("MEAN");
-	liststats (list, dlist, m, &stats);
+	liststats (list, dlist, NULL, m, &stats);
 	stats.sigma = MAX (MIN_ERROR, stats.sigma);
 
@@ -1015,12 +835,10 @@
     }
   }
-  initstats (STATMODE);
   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);
-
-
 }
 
@@ -1072,6 +890,7 @@
   }
 
-  // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1112,5 +931,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -1151,5 +972,7 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/global_stats.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/global_stats.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/global_stats.c	(revision 33886)
@@ -12,6 +12,4 @@
 
   // struct timeval startTimer, stopTimer;
-
-  initstats ("MEAN");
 
   fprintf (stderr, "\n");
@@ -60,8 +58,4 @@
   fprintf (stderr, "dMcal mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
   fprintf (stderr, "chisq mosaic:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);
-  
-
-  initstats (STATMODE);
-
 }
 
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c	(revision 33886)
@@ -33,6 +33,4 @@
     exit (1);
   }
-
-  initstats (STATMODE);
 
   IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
Index: /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/liststats.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/liststats.c	(revision 33885)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/liststats.c	(revision 33886)
@@ -1,30 +1,26 @@
 # include "relphot.h"
 
-enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
-      M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
+void liststats_setmode (StatType *stats, char *strmode) {
 
-static int statmode;
+  stats->statmode = -1;
+  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
+  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
+  if (!strcmp (strmode, "WT_MEAN"))          { stats->statmode = STATS_WT_MEAN; return; }
+  if (!strcmp (strmode, "INNER_MEAN"))       { stats->statmode = STATS_INNER_MEAN; return; }
+  if (!strcmp (strmode, "INNER_WTMEAN"))     { stats->statmode = STATS_INNER_WTMEAN; return; }
+  if (!strcmp (strmode, "CHI_INNER_MEAN"))   { stats->statmode = STATS_CHI_INNER_MEAN; return; }
+  if (!strcmp (strmode, "CHI_INNER_WTMEAN")) { stats->statmode = STATS_CHI_INNER_WTMEAN; return; }
 
-void initstats (char *mode) {
-
-  statmode = -1;
-  if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
-  if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
-  if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN;
-  if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
-  if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
-  if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
-  if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
-
-  if (statmode == -1) {
-    fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
-    exit (1);
-  }
+  fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode);
+  exit (1);
 }
 
-int liststats (double *value, double *dvalue, int N, StatType *stats) {
+int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) {
   
-  int i, ks, ke, Nm;
-  double Mo, dMo, M, dM, X2, dS, *chi;
+  int i, ks, ke;
+  double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;
+
+  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
+  myAssert (weight, "work out logic for NULL weight");
 
   ke = ks = dMo = 0;
@@ -37,66 +33,116 @@
   if (N < 1) return (FALSE);
 
-  dsortpair (value, dvalue, N);
-  stats[0].median = value[(int)(0.5*N)];
-  stats[0].min    = value[0];
-  stats[0].max    = value[N-1];
   int N80 = MIN (N-1, 0.8*N);
   int N20 = MAX (0, 0.2*N);
+
+  if (weight) {
+    dsortthree (value, dvalue, weight, N);
+  } else {
+    dsortpair (value, dvalue, N);
+  }
+
+  // these values do not depend on the errors or weighting scheme
+  stats[0].median  = value[(int)(0.5*N)];
+  stats[0].min     = value[0];
+  stats[0].max     = value[N-1];
   stats[0].Upper80 = value[N80];
   stats[0].Lower20 = value[N20];
 
-  switch (statmode) {
-  case M_MEDIAN:
-    ks = 0;
-    ke = N;
-    Mo = stats[0].median;
-    Nm = N;
-    goto chisq;
-    break;
-  case M_MEAN:
-  case M_WT_MEAN:
-    ks = 0;
-    ke = N;
-    break;
-  case M_INNER_MEAN:
-  case M_INNER_WTMEAN:
-  case M_CHI_INNER_MEAN:
-  case M_CHI_INNER_WTMEAN:
-    ks = 0.25*N + 0.50;
-    ke = 0.75*N + 0.25;
-    if (N <= 3) {
+  switch (stats->statmode) {
+    case STATS_MEDIAN:
       ks = 0;
       ke = N;
-    }
-    break;
+      Mo = stats[0].median;
+      Nm = N;
+      goto chisq;
+      break;
+    case STATS_MEAN:
+    case STATS_WT_MEAN:
+      ks = 0;
+      ke = N;
+      break;
+    case STATS_INNER_MEAN:
+    case STATS_INNER_WTMEAN:
+    case STATS_CHI_INNER_MEAN:
+    case STATS_CHI_INNER_WTMEAN:
+      ks = 0.25*N + 0.50;
+      ke = 0.75*N + 0.25;
+      if (N <= 3) {
+	ks = 0;
+	ke = N;
+      }
+      break;
+    case STATS_NONE:
+      myAbort ("undefined stats");
   }    
 
-  if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
+  // for these two modes, I need a vector of the chi-square contribution
+  // I'm actually just using chisq to get the correct sorting order
+  if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
     ALLOCATE (chi, double, N);
     for (i = 0; i < N; i++) {
       chi[i] = (value[i] - stats[0].median) / dvalue[i];
     }
-    dsortthree (chi, value, dvalue, N);
+    if (weight) {
+      dsortfour (chi, value, dvalue, weight, N);
+    } else {
+      dsortthree (chi, value, dvalue, N);
+    }
     free (chi);
   }
 
-  /* calculating the per-star offset based on the weighted average */
-  M = dM = Nm = 0;
-  if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
-    for (i = ks; i < ke; i++) {
-      M   += value[i] / SQ (dvalue[i]);
-      dM  += 1.0 / SQ (dvalue[i]);
-      Nm  ++;  
-    }	
-    Mo = M / dM;
-    dMo = sqrt (1.0 / dM);
+  int WeightedMean = FALSE;
+  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
+  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
+  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
+
+  /* calculating the per-star offset based on the desired weighting scheme */
+  M = dM = Nm = W = R = 0;
+  if (weight) {
+    // the weight value is multiplied by whichever nominal weighting scheme is provided 
+    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
+    // and 0.01 for under-weight
+    if (WeightedMean) {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i] / SQ(dvalue[i]);
+	W  +=            weight[i] / SQ(dvalue[i]);
+	dM += SQ (weight[i] / dvalue[i]);
+	R  += weight[i] / SQ(dvalue[i]);
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM / R);
+    } else {
+      for (i = ks; i < ke; i++) {
+	M  += value[i] * weight[i];
+	W  +=            weight[i];
+	dM += SQ (weight[i] * dvalue[i]);
+	R  += weight[i];
+	Nm += 1.0;  
+      }	
+      Mo  = M / W;
+      dMo = sqrt (dM / R);
+    }
   } else {
-    for (i = ks; i < ke; i++) {
-      M   += value[i];
-      dM  += SQ (dvalue[i]);
-      Nm  ++;  
-    }	
-    Mo = M / (double) Nm;
-    dMo = sqrt (dM / (double) Nm);
+    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
+    if (WeightedMean) {
+      // weighted by inverse-variance
+      for (i = ks; i < ke; i++) {
+	M   += value[i] / SQ (dvalue[i]);
+	dM  += 1.0 / SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / dM;
+      dMo = sqrt (1.0 / dM);
+    } else {
+      // pure un-weighted
+      for (i = ks; i < ke; i++) {
+	M   += value[i];
+	dM  += SQ (dvalue[i]);
+	Nm  += 1.0;  
+      }	
+      Mo = M / Nm;
+      dMo = sqrt (dM) / Nm;
+    }
   }
 
@@ -111,5 +157,5 @@
   }
   X2 = X2 / (Nm - 1);
-  dS = sqrt (dS / Nm);
+  dS = sqrt (dS / (Nm - 1));
 
   stats[0].mean  = Mo;
@@ -122,2 +168,11 @@
 }
 
+// we could define the weight to be the only scale factor:
+// \mu      = \sum (value_i * weight_i) / \sum (weight_i)
+// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
+// R = \sum (weight_i^2)
+
+// or, we could define the weight to be a scale factor times the inverse error:
+// \mu      = \sum (value_i * weight_i / sigma_i) / \sum (weight_i)
+// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2) 
+// R = \sum (weight_i^2)