Index: /trunk/Ohana/src/relphot/include/relphot.h
===================================================================
--- /trunk/Ohana/src/relphot/include/relphot.h	(revision 41670)
+++ /trunk/Ohana/src/relphot/include/relphot.h	(revision 41671)
@@ -163,12 +163,16 @@
 
 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
+  // these modes are primary and mutually exclusive:
+  STATS_NONE               = 0x0000, 
+  STATS_MEAN               = 0x0001, 
+  STATS_MEDIAN             = 0x0002, 
+  STATS_WT_MEAN            = 0x0003, 
+  STATS_INNER_MEAN         = 0x0004, 
+  STATS_INNER_WTMEAN       = 0x0005, 
+  STATS_CHI_INNER_MEAN     = 0x0006, 
+  STATS_CHI_INNER_WTMEAN   = 0x0007, 
+  // these modes are additional options
+  STATS_PRIMARY            = 0x0007, // use this to mask out the optional bits
+  STATS_VARSTATS           = 0x0010, 
 } ListStatsMode;
 
@@ -185,4 +189,8 @@
   double Lower20;
   double Lower10;
+  double Upper90Nsig;
+  double Upper80Nsig;
+  double Lower20Nsig;
+  double Lower10Nsig;
   double total;
   int    Nmeas;
@@ -419,4 +427,5 @@
 double IMFIT_SYS_SIGMA_LIM;
 double CLOUD_TOLERANCE;
+int    VARIABILITY_STATS;
 
 int    PARALLEL;
Index: /trunk/Ohana/src/relphot/src/args.c
===================================================================
--- /trunk/Ohana/src/relphot/src/args.c	(revision 41670)
+++ /trunk/Ohana/src/relphot/src/args.c	(revision 41671)
@@ -434,4 +434,10 @@
     remove_argument (N, &argc, argv);
     KEEP_UBERCAL = FALSE;
+  }
+
+  VARIABILITY_STATS = FALSE;
+  if ((N = get_argument (argc, argv, "-varstats"))) {
+    remove_argument (N, &argc, argv);
+    VARIABILITY_STATS = TRUE;
   }
 
@@ -865,4 +871,10 @@
 
 
+  VARIABILITY_STATS = FALSE;
+  if ((N = get_argument (argc, argv, "-varstats"))) {
+    remove_argument (N, &argc, argv);
+    VARIABILITY_STATS = TRUE;
+  }
+
   MIN_ERROR = 0.001;
   if ((N = get_argument (argc, argv, "-minerror"))) {
Index: /trunk/Ohana/src/relphot/src/liststats.c
===================================================================
--- /trunk/Ohana/src/relphot/src/liststats.c	(revision 41670)
+++ /trunk/Ohana/src/relphot/src/liststats.c	(revision 41671)
@@ -3,5 +3,5 @@
 void liststats_setmode (StatType *stats, char *strmode) {
 
-  stats->statmode = -1;
+  stats->statmode = STATS_NONE;
   if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
   if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
@@ -17,17 +17,21 @@
 
 int liststats_init (StatType *stats) {
-  stats->median  = NAN;
-  stats->mean    = NAN;
-  stats->sigma   = NAN;
-  stats->error   = NAN;
-  stats->chisq   = NAN;
-  stats->min     = NAN;
-  stats->max     = NAN;
-  stats->Upper80 = NAN;
-  stats->Lower20 = NAN;
-  stats->Upper90 = NAN;
-  stats->Lower10 = NAN;
-  stats->total   = NAN;
-  stats->Nmeas = 0;
+  stats->median      = NAN;
+  stats->mean        = NAN;
+  stats->sigma       = NAN;
+  stats->error       = NAN;
+  stats->chisq       = NAN;
+  stats->min         = NAN;
+  stats->max         = NAN;
+  stats->Upper80     = NAN;
+  stats->Lower20     = NAN;
+  stats->Upper90     = NAN;
+  stats->Lower10     = NAN;
+  stats->Upper80Nsig = NAN;
+  stats->Lower20Nsig = NAN;
+  stats->Upper90Nsig = NAN;
+  stats->Lower10Nsig = NAN;
+  stats->total       = NAN;
+  stats->Nmeas       = 0;
   return TRUE;
 }
@@ -39,4 +43,5 @@
 
   myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
+  ListStatsMode myMode = stats->statmode & STATS_PRIMARY; // exclude the option bits
 
   liststats_init (stats);
@@ -75,5 +80,5 @@
   stats[0].Lower10 = value[N10];
 
-  switch (stats->statmode) {
+  switch (myMode) {
     case STATS_MEDIAN:
       ks = 0;
@@ -100,4 +105,5 @@
       break;
     case STATS_NONE:
+    default:
       myAbort ("undefined stats");
   }    
@@ -105,5 +111,5 @@
   // 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)) {
+  if ((myMode == STATS_CHI_INNER_MEAN) || (myMode == STATS_CHI_INNER_WTMEAN)) {
     ALLOCATE (chi, double, N);
     for (i = 0; i < N; i++) {
@@ -119,7 +125,7 @@
 
   int WeightedMean = FALSE;
-  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
-  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
-  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
+  WeightedMean |= (myMode == STATS_WT_MEAN);
+  WeightedMean |= (myMode == STATS_INNER_WTMEAN);
+  WeightedMean |= (myMode == STATS_CHI_INNER_WTMEAN);
   if (!dvalue) WeightedMean = FALSE;  // warn the user?
 
@@ -310,13 +316,10 @@
       
   // calculate the weight thresholds to mask the bad points:
-  // double Sum_W = 0.0;
   for (int i = 0; i < Npoints; i++) {
     dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     dataset->wtlist[i] = dataset->wtvals[i];
-    // Sum_W += dataset->wtvals[i];
   }
   dsort (dataset->wtlist, Npoints);
   double WtMedian = (Npoints % 2) ? dataset->wtlist[midpt] : 0.5*(dataset->wtlist[midpt] + dataset->wtlist[midpt-1]);
-  // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
   double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
 
@@ -348,4 +351,37 @@
   stats->sigma = sqrt (dSig / (Nkeep - 1));
 
+  // if percentile ranges are desired, calculate them
+  if (stats->statmode & STATS_VARSTATS) {
+    // save ykeep values in a vector to be sorted
+    for (int i = 0; i < Nkeep; i++) {
+      dataset->wtvals[i] = dataset->ykeep[i];
+    }
+    dsort (dataset->wtvals, Nkeep);
+
+    int N90 = MIN (Nkeep-1, 0.9*Nkeep);
+    int N80 = MIN (Nkeep-1, 0.8*Nkeep);
+    int N20 = MAX (0, 0.2*Nkeep);
+    int N10 = MAX (0, 0.1*Nkeep);
+
+    stats->Upper90 = dataset->wtvals[N90];
+    stats->Upper80 = dataset->wtvals[N80];
+    stats->Lower20 = dataset->wtvals[N20];
+    stats->Lower10 = dataset->wtvals[N10];
+    
+    int NkeepMid = 0.5 * Nkeep;
+    double ClipMedian = (Nkeep % 2) ? dataset->wtlist[NkeepMid] : 0.5*(dataset->wtlist[NkeepMid] + dataset->wtlist[NkeepMid-1]);
+
+    // save ykeep values in a vector to be sorted
+    for (int i = 0; i < Nkeep; i++) {
+      dataset->wtvals[i] = (dataset->ykeep[i] - ClipMedian) / dataset->dykeep[i];
+    }
+    dsort (dataset->wtvals, Nkeep);
+
+    stats->Upper90Nsig = dataset->wtvals[N90];
+    stats->Upper80Nsig = dataset->wtvals[N80];
+    stats->Lower20Nsig = dataset->wtvals[N20];
+    stats->Lower10Nsig = dataset->wtvals[N10];
+  }
+
   int Nboot = 0;
   for (int iboot = 0; iboot < NBOOTSTRAP; iboot++) {
Index: /trunk/Ohana/src/relphot/src/reload_catalogs.c
===================================================================
--- /trunk/Ohana/src/relphot/src/reload_catalogs.c	(revision 41670)
+++ /trunk/Ohana/src/relphot/src/reload_catalogs.c	(revision 41671)
@@ -237,4 +237,5 @@
     if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
     // if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
+    if (VARIABILITY_STATS) { strextend (&command, "-varstats"); }
     if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
 
Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c
===================================================================
--- /trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41670)
+++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c	(revision 41671)
@@ -1,3 +1,13 @@
 # include "relphot.h"
+
+# define UPPER90_VAL MpsfStk
+# define UPPER80_VAL FpsfStk
+# define LOWER20_VAL MkronStk
+# define LOWER10_VAL FkronStk
+# define UPPER90_SIG MpsfWrp 
+# define UPPER80_SIG FpsfWrp 
+# define LOWER20_SIG MkronWrp
+# define LOWER10_SIG FkronWrp
+
 
 # if (0)
@@ -435,4 +445,7 @@
     StatType *psfstats = &results->psfstats;
 
+    // if isSetMrelFinal and varstats are requested, request percentile ranges
+    if (isSetMrelFinal && VARIABILITY_STATS) psfstats->statmode |= STATS_VARSTATS;
+
     int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
     if (NrankingPSF < Nminmeas) { 
@@ -469,4 +482,15 @@
 	secfilt[Nsec].Mmax     = psfstats->max;
 	secfilt[Nsec].Mmin     = psfstats->min;
+      }
+      if (VARIABILITY_STATS) {
+	// XXX NOTE: these are overloaded on warp and stack mags and fluxes
+	secfilt[Nsec].UPPER90_VAL = psfstats->Upper90;
+	secfilt[Nsec].UPPER80_VAL = psfstats->Upper80;
+	secfilt[Nsec].LOWER20_VAL = psfstats->Lower20;
+	secfilt[Nsec].LOWER10_VAL = psfstats->Lower10;
+	secfilt[Nsec].UPPER90_SIG = psfstats->Upper90Nsig;
+	secfilt[Nsec].UPPER80_SIG = psfstats->Upper80Nsig;
+	secfilt[Nsec].LOWER20_SIG = psfstats->Lower20Nsig;
+	secfilt[Nsec].LOWER10_SIG = psfstats->Lower10Nsig;
       }
       secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];