Index: trunk/Ohana/src/relphot/Makefile
===================================================================
--- trunk/Ohana/src/relphot/Makefile	(revision 39635)
+++ trunk/Ohana/src/relphot/Makefile	(revision 39636)
@@ -20,7 +20,10 @@
 relphot_client: $(BIN)/relphot_client.$(ARCH)
 
+test_liststats: $(BIN)/test_liststats.$(ARCH)
+
 install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client
 
 RELPHOT = \
+$(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
 $(SRC)/GridOps.$(ARCH).o	 \
@@ -79,4 +82,5 @@
 
 RELPHOT_CLIENT = \
+$(SRC)/StatDataSetOps.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
 $(SRC)/GridOps.$(ARCH).o	 \
@@ -120,2 +124,11 @@
 $(RELPHOT_CLIENT): $(INC)/relphot.h
 $(BIN)/relphot_client.$(ARCH): $(RELPHOT_CLIENT)
+
+TEST_LISTSTATS = \
+$(SRC)/StatDataSetOps.$(ARCH).o	 \
+$(SRC)/liststats.$(ARCH).o	 \
+$(SRC)/test_liststats.$(ARCH).o
+
+$(TEST_LISTSTATS): $(INC)/relphot.h
+$(BIN)/test_liststats.$(ARCH): $(TEST_LISTSTATS)
+
Index: trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- trunk/Ohana/src/relphot/include/relphot.h	(revision 39635)
+++ trunk/Ohana/src/relphot/include/relphot.h	(revision 39636)
@@ -119,4 +119,5 @@
   int    *ranking;	      // weights to use for mean mags
   int    *measSeq;	      // weights to use for mean mags
+  int    *msklist;	      // mask modifications
   int     Nlist;
 } StatDataSet;
@@ -428,4 +429,6 @@
 int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
 int           liststats_init      PROTO((StatType *stats));
+int           liststats_irls      PROTO((StatDataSet *dataset, int Npoints, StatType *stats));
+
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
 Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile));
Index: trunk/Ohana/src/relphot/src/StarOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/StarOps.c	(revision 39635)
+++ trunk/Ohana/src/relphot/src/StarOps.c	(revision 39636)
@@ -52,34 +52,4 @@
   value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
   return (value);
-}
-
-StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax) {
-
-  int i;
-
-  StatDataSet *dataset = NULL;
-  ALLOCATE (dataset, StatDataSet, Nsecfilt);
-  for (i = 0; i < Nsecfilt; i++) {
-    ALLOCATE (dataset[i].flxlist, double, Nmax);
-    ALLOCATE (dataset[i].wgtlist, double, Nmax);
-    ALLOCATE (dataset[i].errlist, double, Nmax);
-    ALLOCATE (dataset[i].ranking,    int, Nmax);
-    ALLOCATE (dataset[i].measSeq,    int, Nmax);
-  }  
-  return dataset;
-}
-
-void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) {
-
-  int i;
-
-  for (i = 0; i < Nsecfilt; i++) {
-    FREE (dataset[i].flxlist);
-    FREE (dataset[i].wgtlist);
-    FREE (dataset[i].errlist);
-    FREE (dataset[i].ranking);
-    FREE (dataset[i].measSeq);
-  }  
-  FREE (dataset);
 }
 
Index: trunk/Ohana/src/relphot/src/StatDataSetOps.c
===================================================================
--- trunk/Ohana/src/relphot/src/StatDataSetOps.c	(revision 39636)
+++ trunk/Ohana/src/relphot/src/StatDataSetOps.c	(revision 39636)
@@ -0,0 +1,34 @@
+# include "relphot.h"
+
+StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax) {
+
+  int i;
+
+  StatDataSet *dataset = NULL;
+  ALLOCATE (dataset, StatDataSet, Nsecfilt);
+  for (i = 0; i < Nsecfilt; i++) {
+    ALLOCATE (dataset[i].flxlist, double, Nmax);
+    ALLOCATE (dataset[i].wgtlist, double, Nmax);
+    ALLOCATE (dataset[i].errlist, double, Nmax);
+    ALLOCATE (dataset[i].ranking,    int, Nmax);
+    ALLOCATE (dataset[i].measSeq,    int, Nmax);
+    ALLOCATE (dataset[i].msklist,    int, Nmax);
+  }  
+  return dataset;
+}
+
+void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) {
+
+  int i;
+
+  for (i = 0; i < Nsecfilt; i++) {
+    FREE (dataset[i].flxlist);
+    FREE (dataset[i].wgtlist);
+    FREE (dataset[i].errlist);
+    FREE (dataset[i].ranking);
+    FREE (dataset[i].measSeq);
+    FREE (dataset[i].msklist);
+  }  
+  FREE (dataset);
+}
+
Index: trunk/Ohana/src/relphot/src/liststats.c
===================================================================
--- trunk/Ohana/src/relphot/src/liststats.c	(revision 39635)
+++ trunk/Ohana/src/relphot/src/liststats.c	(revision 39636)
@@ -188,4 +188,176 @@
 }
 
+// These should probably be tunable:
+# define MAX_ITERATIONS 10
+# define FIT_TOLERANCE 1e-4
+# define FLT_TOLERANCE 1e-6
+# define WEIGHT_THRESHOLD 0.3
+# define NBOOTSTRAP 100
+
+int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts);
+double VectorFractionInterpolate (double *values, float fraction, int Npts);
+double weight_cauchy (double x);
+
+// this is a zero-order fit (constant value only)
+int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats) {
+
+  double value;
+
+  liststats_init (stats);
+
+  // OLS
+  if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
+  
+  // XXX add to dataset elements?
+  ALLOCATE_PTR (wt, double, Npoints);
+
+  int converged = FALSE;
+  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
+    
+    for (int i = 0; i < Npoints; i++) {
+      // we are only including the formal error, not the weight in the definition of wt[]
+      wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+    }
+    
+    double oldValue = value;
+    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wt, Npoints)) {
+      value = oldValue;
+      break;
+    }
+    
+    converged = TRUE;
+    if ((fabs(value - oldValue) > FIT_TOLERANCE * fabs(value)) && 
+	(fabs(value - oldValue) > FLT_TOLERANCE)) {
+      converged = FALSE;
+    }
+  }
+  stats->mean = value;
+      
+  // calculate the weight thresholds to mask the bad points:
+  double Sum_W = 0.0;
+  for (int i = 0; i < Npoints; i++) {
+    wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
+    Sum_W += wt[i];
+  }
+  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
+
+  // generate the unmasked subset
+  // XXX add these to the dataset elements?
+  ALLOCATE_PTR ( ykeep, double, Npoints);
+  ALLOCATE_PTR (dykeep, double, Npoints);
+  ALLOCATE_PTR (wtkeep, double, Npoints);
+
+  // save unmasked points
+  int Nkeep = 0;
+  stats->min = NAN;
+  stats->max = NAN;
+  double dChi = 0.0;
+  double dSig = 0.0;
+  for (int i = 0; i < Npoints; i++) {
+    if ((wt[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
+      dataset->msklist[i] = TRUE; // mark the masked points
+      continue;
+    }
+     ykeep[Nkeep] = dataset->flxlist[i];
+    dykeep[Nkeep] = dataset->errlist[i];
+    wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
+    Nkeep ++;
+    
+    // record mean, error, chisq, min, max, sigma, Nmeas (unmasked points)
+    stats->min = isfinite(stats->min) ? MIN(dataset->flxlist[i], stats->min) : dataset->flxlist[i];
+    stats->max = isfinite(stats->max) ? MAX(dataset->flxlist[i], stats->max) : dataset->flxlist[i];
+    double dValue2 = SQ(dataset->flxlist[i] - value);
+    dChi += dValue2 / SQ (dataset->errlist[i]);
+    dSig += dValue2;
+  }
+  stats->Nmeas = Nkeep;
+  stats->chisq = dChi / (Nkeep - 1);
+  stats->sigma = sqrt (dSig / (Nkeep - 1));
+
+  // bootstrap resampling to generate the errorbars
+  // XXX add these to the dataset elements?
+  ALLOCATE_PTR (ysample,  double, Nkeep);
+  ALLOCATE_PTR (dysample, double, Nkeep);
+  ALLOCATE_PTR (wtsample, double, Nkeep);
+  ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
+
+  int Nboot = 0;
+  for (int iboot = 0; iboot < NBOOTSTRAP; iboot++) {
+    
+    // resample
+    for (int i = 0; i < Nkeep; i++) {
+      // I need to draw Npoints random entries from 'points' with replacement:
+      int N = Nkeep * drand48();
+       ysample[i] =  ykeep[N];
+      dysample[i] = dykeep[N];
+      wtsample[i] = wtkeep[N];
+    }
+
+    if (!fit_least_squares (&value, ysample, dysample, wtsample, NULL, Nkeep)) continue;
+
+    bvalue[Nboot] = value;
+    Nboot ++;
+  }
+
+  dsort (bvalue, Nboot);
+  
+  double Slo = VectorFractionInterpolate (bvalue, 0.158655, Nboot);
+  double Shi = VectorFractionInterpolate (bvalue, 0.841345, Nboot);
+  stats->error = (Shi - Slo) / 2.0;
+
+  free (bvalue);
+  free ( ysample);
+  free (dysample);
+  free (wtsample);
+  free ( ykeep);
+  free (dykeep);
+  free (wtkeep);
+  free (wt);
+
+  return TRUE;
+}
+
+// wgt is externally-supplied weight, wt is optional
+int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts) { 
+
+  int i;
+
+  double S0 = 0;
+  double S1 = 0;
+
+  /* perform linear fit */
+  for (i = 0; i < Npts; i++, y++, dy++, wgt++) {
+    if (!finite(*y)) continue;
+
+    // wt is optional
+    double dY = wt ? wt[i] * (*wgt) / SQ(*dy) : (*wgt) / SQ(*dy);
+
+    S0 +=    dY;
+    S1 += *y*dY;
+  }
+  if (S0 == 0.0) return FALSE;
+  *fit = S1 / S0;
+  return TRUE;
+}
+
+double weight_cauchy (double x) {
+  double r = x / 2.385;
+  return (1.0 / (1.0 + SQ(r)));
+}
+
+double VectorFractionInterpolate (double *values, float fraction, int Npts) {
+
+  float F = fraction * Npts;
+  int   N = fraction * Npts;
+
+  if (N < 0        ) return NAN;
+  if (N >= Npts - 2) return NAN;
+
+  // interpolate between N,N+1
+    
+  double S = (F - N) * (values[N+1] - values[N]) + values[N];
+  return S;
+}
+
 // we could define the weight to be the only scale factor:
 // \mu      = \sum (value_i * weight_i) / \sum (weight_i)
Index: trunk/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 39635)
+++ trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 39636)
@@ -10,8 +10,12 @@
 
 int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats);
-int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
+int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats);
+int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag);
 void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
 void sort_entry_by_offset (double *offset, int *entry, int N);
 void sort_StatDataSet (StatDataSet *dataset);
+int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats);
+
+# define MAG_STATS_BY_RANKING magStatsByRankingIRLS
 
 # define UBERCAL_WEIGHT 100.0
@@ -255,4 +259,5 @@
 	results->aperData[Nsec].ranking[Nap] = measureRank[k];
 	results->aperData[Nsec].measSeq[Nap] = k;
+	results->aperData[Nsec].msklist[Nap] = 0;
 	results->aperData[Nsec].Nlist ++;
       }
@@ -270,4 +275,5 @@
 	results->kronData[Nsec].ranking[Nkron] = measureRank[k];
 	results->kronData[Nsec].measSeq[Nkron] = k;
+	results->kronData[Nsec].msklist[Nkron] = 0;
 	results->kronData[Nsec].Nlist ++;
       }
@@ -283,4 +289,5 @@
       results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
       results->psfData[Nsec].measSeq[Npsf] = k;
+      results->psfData[Nsec].msklist[Npsf] = 0;
       results->psfData[Nsec].Nlist ++;
     }
@@ -362,5 +369,5 @@
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *psfstats = &results->psfstats;
-    int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
+    int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     if (Nranking < Nminmeas) { 
       secfilt[Nsec].flags |= ID_OBJ_FEW;
@@ -370,5 +377,5 @@
       secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
     }
-    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
+    int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
 
     // when running -averages, we have no information about the images, so we cannot set this
@@ -383,5 +390,5 @@
 
       // mark the measurements matching this ranking 
-      markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_PHOTOM_PSF);
+      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
 
       if (Nranking) {
@@ -391,5 +398,4 @@
 	secfilt[Nsec].Mmin   = psfstats->min;
       }
-
       secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
       secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
@@ -399,5 +405,5 @@
 
       StatType *apstats = &results->apstats;
-      Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
+      Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
       if (Nranking) {
 	secfilt[Nsec].Map     = apstats->mean; 
@@ -406,8 +412,9 @@
 	secfilt[Nsec].NusedAp = Nranking;
       }
-      markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER);
+      int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
 
       StatType *kronstats = &results->kronstats;
-      Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
+      Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
       if (Nranking) {
 	secfilt[Nsec].Mkron     = kronstats->mean; 
@@ -416,5 +423,6 @@
 	secfilt[Nsec].NusedKron = Nranking;
       }
-      markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON);
+      int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
+      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
 
       // does this object appear extended in > 50% of measurements?
@@ -423,5 +431,5 @@
       }
 
-      switch (minRank) {
+      switch (minRankPSF) {
 	case 0:
 	  secfilt[Nsec].flags |= ID_SECF_RANK_0;
@@ -853,4 +861,5 @@
       results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
       results->psfData[Nsec].measSeq[Npsf] = k;
+      results->psfData[Nsec].msklist[Npsf] = 0;
       results->psfData[Nsec].Nlist ++;
     }     
@@ -871,4 +880,5 @@
       results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0;
       results->aperData[Nsec].measSeq[Naper] = k;
+      results->aperData[Nsec].msklist[Naper] = 0;
       results->aperData[Nsec].Nlist ++;
     }
@@ -884,4 +894,5 @@
       results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
       results->kronData[Nsec].measSeq[Nkron] = k;
+      results->kronData[Nsec].msklist[Nkron] = 0;
       results->kronData[Nsec].Nlist ++;
     }
@@ -905,7 +916,9 @@
     }
 
+    int Nranking, minRank;
+
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *psfstats = &results->psfstats;
-    int Nranking = magStatsByRankingClipped (&results->psfData[Nsec], psfstats);
+    Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
     if (Nranking) {
       secfilt[Nsec].FpsfWrp  = psfstats->mean;
@@ -915,10 +928,10 @@
       secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
     }
-    int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
-    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED);
+    minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
+    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
 
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *apstats = &results->apstats;
-    Nranking = magStatsByRankingClipped (&results->aperData[Nsec], apstats);
+    Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
     if (Nranking) {
       secfilt[Nsec].FapWrp     = apstats->mean;
@@ -928,8 +941,10 @@
       secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
     }
+    minRank = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
+    markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
 
     // if too few valid measurements meet the minimum criteria, go to the next entry
     StatType *kronstats = &results->kronstats;
-    Nranking = magStatsByRankingClipped (&results->kronData[Nsec], kronstats);
+    Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
     if (Nranking) {
       secfilt[Nsec].FkronWrp     = kronstats->mean;
@@ -939,4 +954,6 @@
       secfilt[Nsec].MkronWrp     = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
     }
+    minRank = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
+    markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
 
     secfilt[Nsec].Nwarp     = results->Nmeas[Nsec];
@@ -967,5 +984,4 @@
 }
 
-# if (0)
 int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats) {
 
@@ -985,8 +1001,7 @@
   for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);
 
-  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
+  liststats_irls (dataset, Nranking, stats);
   return (Nranking);
 }
-# endif
 
 // outlier warp measurements are driving bad mean values.  
@@ -1061,5 +1076,5 @@
 }
 
-int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) {
+int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag) {
   int i;
 
@@ -1067,8 +1082,9 @@
     int k = dataset->measSeq[i];
     if (dataset->ranking[i] > minrank) {
-      measure[k].dbFlags &= ~flags;
+      measure[k].dbFlags &= ~keepflag;
     } else {
-      measure[k].dbFlags |= flags;
-    }
+      measure[k].dbFlags |= keepflag;
+    }
+    if (dataset->msklist[i]) measure[k].dbFlags |= maskflag;
   }
   return TRUE;
Index: trunk/Ohana/src/relphot/src/test_liststats.c
===================================================================
--- trunk/Ohana/src/relphot/src/test_liststats.c	(revision 39636)
+++ trunk/Ohana/src/relphot/src/test_liststats.c	(revision 39636)
@@ -0,0 +1,155 @@
+# include "relphot.h"
+
+int mkstar (StatDataSet *dataset, double flux, int Npoints, int Nbad, int GaussianOutliers);
+
+static int Nrand = 0;
+double drand48_cnt () {
+  double value = drand48();
+  Nrand ++;
+  return value;
+}
+
+typedef enum {
+  FIT_FLUX_NONE,
+  FIT_FLUX_MEAN,
+  FIT_FLUX_IRLS,
+} TestModeType;
+
+int main (int argc, char **argv) {
+  
+  int N_STARS     = 3000;
+  int N_POINTS    =  100;
+  int N_OUTLIERS  =   10;
+  // int N_BOOTSTRAP =  100;
+
+  int N;
+  if ((N = get_argument (argc, argv, "-Nstars"))) {
+    remove_argument (N, &argc, argv);
+    N_STARS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Npoints"))) {
+    remove_argument (N, &argc, argv);
+    N_POINTS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-Noutliers"))) {
+    remove_argument (N, &argc, argv);
+    N_OUTLIERS = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int GaussianOutliers = TRUE;
+  if ((N = get_argument (argc, argv, "-flat-outliers"))) {
+    remove_argument (N, &argc, argv);
+    GaussianOutliers = FALSE;
+  }
+
+  if (argc != 2) {
+    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
+    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -flat-outliers\n");
+    exit (2);
+  }
+
+  TestModeType mode = FIT_FLUX_NONE;
+  if (!strcasecmp (argv[1], "irls")) mode = FIT_FLUX_IRLS;
+  if (!strcasecmp (argv[1], "mean")) mode = FIT_FLUX_MEAN;
+  if (!mode) {
+    fprintf (stderr, " mode options: mean, irls\n");
+    exit (2);
+  }
+
+  // plan_tests (14);
+
+  // diag ("relastro fitpm tests");
+  
+  // init the random seed
+  { 
+    struct timeval now;
+    gettimeofday (&now, NULL);
+    long A = now.tv_sec + now.tv_usec * 1000000;
+    srand48(A);
+    // srand48(1);
+  }
+
+  ohana_gaussdev_init ();
+
+  // allow for 5 filters and 1-to-1 outliers
+  StatDataSet *dataset = StatDataSetAlloc (1, N_POINTS + N_OUTLIERS);
+
+  StatType stats;
+  liststats_init (&stats);
+
+  int Nstars = N_STARS;
+  for (int i = 0; i < Nstars; i++) {
+
+    double lflux = 2.0 + 3.0*drand48_cnt(); // uniform in log-flux
+    double flux = pow(10.0, lflux);
+
+    int Npoints   = N_POINTS;
+    int Noutliers = N_OUTLIERS;
+
+    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
+    mkstar (dataset, flux, Npoints, Noutliers, GaussianOutliers);
+
+    switch (mode) {
+      case FIT_FLUX_MEAN:
+	liststats_setmode (&stats, "MEAN");
+	liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Npoints + Noutliers, &stats);
+	break;
+
+      case FIT_FLUX_IRLS:
+	liststats_irls (dataset, Npoints + Noutliers, &stats);
+	break;
+
+      default: myAbort ("programming error");
+    }
+
+    fprintf (stdout, "%8.2f %8.2f : %6.2f %6.2f : %6.3f : %4d %4d %4d\n",
+	     flux, stats.mean, stats.sigma, stats.error, stats.chisq, Npoints, Noutliers, stats.Nmeas);
+  }
+  // return exit_status();
+  ohana_gaussdev_free();
+
+  StatDataSetFree (dataset, 1);
+
+  ohana_memcheck (TRUE);
+  ohana_memdump (TRUE);
+
+  exit (0);
+}
+
+# define READ_NOISE_SQR 25.0
+int mkstar (StatDataSet *dataset, double flux, int Npoints, int Nbad, int GaussianOutliers) {
+  
+  double FluxError = sqrt(flux);
+
+  for (int i = 0; i < Npoints; i++) {
+    double dF = ohana_gaussdev_rnd(0.0, FluxError);
+
+    double fluxObs = flux + dF;
+    double dFluxObs = sqrt(fluxObs + READ_NOISE_SQR); 
+
+    dataset->flxlist[i] = fluxObs;
+    dataset->errlist[i] = dFluxObs;
+    dataset->wgtlist[i] = 1.0;
+    dataset->ranking[i] = 1;
+    dataset->measSeq[i] = i;
+  }
+
+  for (int i = 0; i < Nbad; i++) {
+
+    // double dF = GaussianOutliers ? ohana_gaussdev_rnd(0.0, 10.0*FluxError) : 10.0*FluxError*(drand48_cnt() - 0.5);
+    double dF = GaussianOutliers ? ohana_gaussdev_rnd(10.0*FluxError, 10.0*FluxError) : 10.0*FluxError*drand48_cnt();
+
+    double fluxObs = flux + dF;
+    double dFluxObs = sqrt(fluxObs + READ_NOISE_SQR); 
+
+    dataset->flxlist[Npoints + i] = fluxObs;
+    dataset->errlist[Npoints + i] = dFluxObs;
+    dataset->wgtlist[Npoints + i] = 1.0;
+    dataset->ranking[Npoints + i] = 1;
+    dataset->measSeq[Npoints + i] = i;
+  }
+
+  return TRUE;
+}
Index: trunk/Ohana/src/relphot/test/mana.sh
===================================================================
--- trunk/Ohana/src/relphot/test/mana.sh	(revision 39636)
+++ trunk/Ohana/src/relphot/test/mana.sh	(revision 39636)
@@ -0,0 +1,63 @@
+
+macro go
+  if ($0 != 2)
+    echo "USAGE: go (datafile)"
+    break
+  end
+
+  data $1
+  read Fo 1 Fx 2 dFo 4 dFx 5 chisq 7 Nstar 9 Nbad 10 Nfit 11
+  set dF = Fo - Fx
+  set n = ramp(Fo)
+
+  # predicted errorbar compared to measured errorbar
+  set dfo = dFo / sqrt(Nstar)
+  lim -n 1 dfo dFx; clear; box; plot dfo dFx
+  line -c red 0 0 to 50 50 
+
+  # number fitted
+  set fFit = Nfit / Nstar
+  lim -n 2 n fFit; clear; box; plot n fFit
+
+  # fitted mean vs truth scaled by sigma
+  set dNsigma = dF / dFx
+  lim -n 0 n dNsigma; clear; box; plot n dNsigma
+
+  histogram dNsigma NdNsigma -5.0 5.0 0.1 -range dx
+  $C0 = 0
+  $C1 = 1.0
+  $C2 = 100
+  $C3 = 0
+  vgauss -q dx NdNsigma con NdNsigmaF
+  vstat -q dNsigma
+
+  lim -n 3 dx NdNsigma; clear; box; plot -x 1 dx NdNsigma; plot -x 0 dx NdNsigmaF -c red
+
+  echo $C0 $C1 $MEAN $SIGMA
+end
+
+macro testset
+  if ($0 != 2)
+    echo "USAGE: testset (mode)"
+    break
+  end
+
+  foreach Nout 0 1 2 4 8 16
+    echo --- $Nout ---
+    exec bin/test_liststats.lin64 $1 -Noutliers $Nout  > test.$1.dat; go test.$1.dat
+    echo -------------
+  end
+end
+
+macro testset.flat
+  if ($0 != 2)
+    echo "USAGE: testset (mode)"
+    break
+  end
+
+  foreach Nout 0 1 2 4 8 16
+    echo --- $Nout ---
+    exec bin/test_liststats.lin64 $1 -Noutliers $Nout -flat-outliers > test.$1.dat; go test.$1.dat
+    echo -------------
+  end
+end