Index: trunk/psastro/src/psastro.h
===================================================================
--- trunk/psastro/src/psastro.h	(revision 25299)
+++ trunk/psastro/src/psastro.h	(revision 25906)
@@ -146,6 +146,9 @@
 
 bool 		  psastroZeroPoint (pmConfig *config);
-bool 		  psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime);
+psVector         *psastroZeroPointReadoutAccum(psVector *dMag, pmReadout *readout, float exptime);
+bool              psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean);
 bool 		  psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe);
+
+psStats          *psastroStatsPercentile (psVector *myVector, float flimit);
 
 // masking functions
Index: trunk/psastro/src/psastroZeroPoint.c
===================================================================
--- trunk/psastro/src/psastroZeroPoint.c	(revision 25299)
+++ trunk/psastro/src/psastroZeroPoint.c	(revision 25906)
@@ -15,4 +15,5 @@
 bool psastroZeroPoint (pmConfig *config) {
 
+    bool status;
     float zeropt, exptime;
     pmChip *chip = NULL;
@@ -27,4 +28,9 @@
     }
 
+    // recipe options
+    bool byExposure = psMetadataLookupBool (&status, recipe, "ZERO.POINT.BY.EXPOSURE");
+    bool useMean    = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
+    float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
+
     // select the input data sources
     pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
@@ -55,5 +61,9 @@
     }
 
+    // if we measure the zero point by exposure, accumulate the dMag values here:
+    psVector *dMag = NULL;
+
     // this loop selects the matched stars for all chips
+    // XXX optionally measure zero point for entire exposure in a single statistic
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
         psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
@@ -71,8 +81,22 @@
 
 		// calculate dMag for the matched stars
-		psastroZeroPointReadout (readout, zeropt, exptime);
-
+		dMag = psastroZeroPointReadoutAccum (dMag, readout, exptime);
+
+		if (!byExposure) {
+		    // calculate dMag for the matched stars just for this readout (well, chip)
+		    psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
+		    psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
+		    psFree (dMag);
+		    dMag = NULL;
+		}
 	    }
 	}
+    }
+
+    if (byExposure) {
+	psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
+	psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
+	psFree (dMag);
+	dMag = NULL;
     }
 
@@ -82,7 +106,7 @@
 
 /**
- * we measure <dMag> and \sigma_dMag and write them to the header
+ * accumulate the dMag values from this readout
  */
-bool psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime) {
+psVector *psastroZeroPointReadoutAccum(psVector *dMag, pmReadout *readout, float exptime) {
 
     bool status;
@@ -90,16 +114,17 @@
     // select the raw objects for this readout
     psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
-    if (rawstars == NULL) return false;
+    if (rawstars == NULL) return dMag;
 
     // select the raw objects for this readout
     psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
-    if (refstars == NULL) return false;
+    if (refstars == NULL) return dMag;
 
     psArray *matches = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.MATCH");
-    if (matches == NULL) return false;
-
-    psVector *dMag  = psVectorAllocEmpty (100, PS_TYPE_F32);
-
-    int Npts = 0;
+    if (matches == NULL) return dMag;
+
+    if (!dMag) {
+	dMag  = psVectorAllocEmpty (100, PS_TYPE_F32);
+    }
+
     for (int i = 0; i < matches->n; i++) {
 
@@ -110,27 +135,42 @@
       pmAstromObj *ref = refstars->data[match->ref];
 
-      dMag->data.F32[Npts] = ref->Mag - raw->Mag - 2.5*log10(exptime);
-      psVectorExtend (dMag, 100, 1);
-      Npts++;
-    }
-
-    psTrace ("psModules.astrom", 4, "Npts: %d\n", Npts);
-
-    if (Npts < 3) {
+      float value = ref->Mag - raw->Mag - 2.5*log10(exptime);
+      psVectorAppend (dMag, value);
+    }
+    return dMag;
+}
+
+bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean) {
+
+    // XXX make this depend on the mode?
+    if (dMag->n < 3) {
       fprintf (stderr, "zero point NaN +/- NaN\n");
-      psFree (dMag);
       return false;
     }
 
-    // stats structure and mask for use in measuring the clipping statistic
-    // this analysis has too few data points to use the robust median method
-    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    if (!psVectorStats (stats, dMag, NULL, NULL, 0)) {
-	psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-	return false;
-    }
-    fprintf (stderr, "zero point %f +/- %f using %d stars; transparency %f\n", stats->clippedMean, stats->clippedStdev, Npts, zeropt - stats->clippedMean);
-
-    psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
+    // the zero point analysis depends on the type of desired statistic.  For comparisons
+    // against a high-quality reference catalog with a good match to the actual filter used, it
+    // is best to use a standard clipped mean or global mean statistic.  If the reference
+    // catalog has some unmodeled extra parameter, as is the case for the synthetic grizy
+    // database vs PS1, then it is best to use some consistent feature in the color
+    // distribution.
+
+    psStats *stats = NULL;
+    if (useMean) {
+	// stats structure and mask for use in measuring the clipping statistic
+	// this analysis has too few data points to use the robust median method
+	stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	if (!psVectorStats (stats, dMag, NULL, NULL, 0)) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by mean");
+	    return false;
+	}
+    } else {
+	stats = psastroStatsPercentile (dMag, edgeFraction);
+	if (!stats) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge");
+	    return false;
+	}
+    }
+    fprintf (stderr, "zero point %f +/- %f using %ld stars; transparency %f\n", stats->clippedMean, stats->clippedStdev, dMag->n, zeropt - stats->clippedMean);
 
     psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_OBS", PS_META_REPLACE, "measured zero point", stats->clippedMean);
@@ -139,8 +179,157 @@
     psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_OFF", PS_META_REPLACE, "measured zero point", zeropt - stats->clippedMean);
 
-    psFree (dMag);
     psFree (stats);
     return true;
 }
+
+#define MAG_RESOLUTION 0.001
+
+// set the bin closest to the corresponding value.  if USE_END is +/- 1,
+// out-of-range saturates on lower/upper bin REGARDLESS of actual value
+#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \
+        psVectorBinaryDisectResult result; \
+        psScalar tmpScalar; \
+        tmpScalar.type.type = PS_TYPE_F32; \
+        tmpScalar.data.F32 = (VALUE); \
+        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
+        switch (result) { \
+          case PS_BINARY_DISECT_PASS: \
+            break; \
+          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
+            psTrace("psastro", 6, "selected bin outside range"); \
+            if (USE_END == -1) { RESULT = 0; } \
+            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
+            break; \
+          case PS_BINARY_DISECT_INVALID_INPUT: \
+          case PS_BINARY_DISECT_INVALID_TYPE: \
+            psAbort ("programming error"); \
+            break; \
+        } }
+
+# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
+        float dX, dY, Xo, Yo, Xt; \
+        if (BIN == BOUNDS->n - 1) { \
+            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
+            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \
+            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
+            Yo = VECTOR->data.F32[BIN]; \
+        } else { \
+            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
+            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \
+            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
+            Yo = VECTOR->data.F32[BIN]; \
+        } \
+        if (dY != 0) { \
+            Xt = (VALUE - Yo)*dX/dY + Xo; \
+        } else { \
+            Xt = Xo; \
+        } \
+        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
+        psTrace("psastro", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
+                Xo, Yo, dX, dY, Xt, VALUE); \
+        RESULT = Xt; }
+
+// measure the edge of the sample at flimit
+// return results in stats->sampleMean, sampleStdev
+// XXX this is a misuse of psStats -- make our own structure?
+psStats *psastroStatsPercentile (psVector *myVector, float flimit) {
+
+    // search for the 'blue' edge of the dMag distribution:
+    // the distribution is not a normal population, but instead has a broad range with fairly hard edges.
+    // construct a histogram and look for the 
+
+    // stats is first used to find the data range
+    psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+    psHistogram *histogram = NULL;      // Histogram of the data
+    psHistogram *cumulative = NULL;     // Cumulative histogram of the data
+    float min = NAN, max = NAN;         // Mimimum and maximum values
+
+    // Get the minimum and maximum values
+    if (!psVectorStats(stats, myVector, NULL, NULL, 0)) goto escape;
+    min = stats->min;
+    max = stats->max;
+    if (isnan(min) || isnan(max)) goto escape;
+    psTrace("psastro", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
+
+    // If all data points have the same value, then we set the appropriate members of stats and return.
+    if (fabs(max - min) <= FLT_EPSILON) {
+	stats->clippedMean = min;
+	stats->clippedStdev = NAN;
+	stats->results |= PS_STAT_CLIPPED_MEAN;
+	stats->results |= PS_STAT_CLIPPED_STDEV;
+	return stats;
+    }
+    
+    // Define the histogram bin size.  
+    float binSize = MAG_RESOLUTION;
+    long numBins = (max - min) / binSize; // Number of bins
+    psTrace("psastro", 6, "Numbins is %ld\n", numBins);
+    psTrace("psastro", 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
+
+    // Generate the histogram
+    histogram = psHistogramAlloc(min, max, numBins);
+    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
+      // if psVectorHistogram returns false, we have a programming error
+      psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for psastroZeroPointAnalysis\n");
+      psFree(histogram);
+      psFree(stats);
+      return NULL;
+    }
+    if (psTraceGetLevel("psastro") >= 8) {
+      PS_VECTOR_PRINT_F32(histogram->bounds);
+      PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    // Convert the specific histogram to a cumulative histogram
+    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
+    cumulative = psHistogramAlloc(min, max, numBins);
+    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
+    for (long i = 1; i < histogram->nums->n; i++) {
+      cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
+      cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
+    }
+    if (psTraceGetLevel("psastro") >= 8) {
+      PS_VECTOR_PRINT_F32(cumulative->bounds);
+      PS_VECTOR_PRINT_F32(cumulative->nums);
+    }
+
+    // Find the bin which contains the first data point above the limit
+    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
+    psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints);
+
+    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
+    long bin;
+    PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0);
+    psTrace("psastro", 6, "The bin is %ld (%.2f to %.2f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]);
+
+    // Linear interpolation to the limit value in bin units
+    float value;
+    PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit);
+    psTrace("psastro", 6, "limit value is %f\n", value);
+
+    stats->clippedMean = value;
+    stats->clippedStdev = 0.0; // XXX derive correct error value
+    stats->results |= PS_STAT_CLIPPED_MEAN;
+    stats->results |= PS_STAT_CLIPPED_STDEV;
+    stats->options = PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV;
+
+    // Clean up
+    psFree(histogram);
+    psFree(cumulative);
+    return stats;
+
+escape:
+    stats->clippedMean = NAN;
+    stats->clippedStdev = NAN;
+    stats->results |= PS_STAT_CLIPPED_MEAN;
+    stats->results |= PS_STAT_CLIPPED_STDEV;
+    stats->options = PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV;
+
+    psFree(histogram);
+    psFree(cumulative);
+
+    return stats;
+}
+
 
 # define ESCAPE(MSG) { \