Index: /branches/eam_branches/20091113/psastro/src/psastroZeroPoint.c
===================================================================
--- /branches/eam_branches/20091113/psastro/src/psastroZeroPoint.c	(revision 26223)
+++ /branches/eam_branches/20091113/psastro/src/psastroZeroPoint.c	(revision 26224)
@@ -30,6 +30,4 @@
     // 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
@@ -86,5 +84,5 @@
 		    // 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);
+		    psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
 		    psFree (dMag);
 		    dMag = NULL;
@@ -96,5 +94,5 @@
     if (byExposure) {
 	psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
-	psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
+	psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
 	psFree (dMag);
 	dMag = NULL;
@@ -141,5 +139,5 @@
 }
 
-bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean) {
+bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe) {
 
     // XXX make this depend on the mode?
@@ -148,4 +146,7 @@
       return false;
     }
+
+    bool status;
+    bool useMean = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
 
     // the zero point analysis depends on the type of desired statistic.  For comparisons
@@ -166,5 +167,5 @@
 	}
     } else {
-	stats = psastroStatsPercentile (dMag, edgeFraction);
+	stats = psastroStatsPercentile (dMag, recipe);
 	if (!stats) {
 	    psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge");
@@ -233,9 +234,14 @@
 // return results in stats->sampleMean, sampleStdev
 // XXX this is a misuse of psStats -- make our own structure?
-psStats *psastroStatsPercentile (psVector *myVector, float flimit) {
+psStats *psastroStatsPercentile (psVector *myVector, psMetadata *recipe) {
 
     // 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 
+
+    bool status;
+    float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
+    int edgeSample = psMetadataLookupS32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE");
+    float edgeSampleFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE.FRACTION");
 
     // stats is first used to find the data range
@@ -245,4 +251,6 @@
     float min = NAN, max = NAN;         // Mimimum and maximum values
 
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);  
+
     // Get the minimum and maximum values
     if (!psVectorStats(stats, myVector, NULL, NULL, 0)) goto escape;
@@ -250,5 +258,5 @@
     max = stats->max;
     if (isnan(min) || isnan(max)) goto escape;
-    psTrace("psastro", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
+    psTrace("psastro", 5, "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.
@@ -264,50 +272,39 @@
     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
+    psTrace("psastro", 5, "Numbins is %ld\n", numBins);
+    psTrace("psastro", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
+
+    // allocate the histogram containers
     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
+
+    // find the mean value:
+    stats->clippedMean = psastroStatsPercentileValue (histogram, cumulative, myVector, edgeFraction);
+
+    int nSubset = myVector->n * edgeSampleFraction;
+    psVector *subset = psVectorAlloc (nSubset, PS_TYPE_F32);
+
+    float Sum = 0.0;
+    float S2 = 0.0;
+    for (int i = 0; i < edgeSample; i++) {
+	
+	// generate the subset vector
+	for (long i = 0; i < nSubset; i++) {
+	    double frnd = psRandomUniform(rng);
+	    int entry = PS_MIN(myVector->n - 1, PS_MAX(0, myVector->n * frnd));
+
+	    subset->data.F32[i] = myVector->data.F32[entry];
+	}
+
+	float value = psastroStatsPercentileValue (histogram, cumulative, subset, edgeFraction);
+
+	Sum += value;
+	S2 += value*value;
+    }
+    psTrace("psastro", 6, "subset stats: Sum: %f, S2: %f, Npts: %d\n", Sum, S2, edgeSample);
+
+    stats->clippedStdev = sqrt(S2 / edgeSample - PS_SQR(Sum/edgeSample));
+    psTrace("psastro", 5, "percentile stats %f +/- %f\n", stats->clippedMean, stats->clippedStdev);
+
     stats->results |= PS_STAT_CLIPPED_MEAN;
     stats->results |= PS_STAT_CLIPPED_STDEV;
@@ -317,4 +314,6 @@
     psFree(histogram);
     psFree(cumulative);
+    psFree(subset);
+    psFree(rng);
     return stats;
 
@@ -332,4 +331,50 @@
 }
 
+
+// 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?
+float psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit) {
+
+    // need to initialize the histogram on each pass
+    psVectorInit (histogram->nums, 0);
+    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
+	// if psVectorHistogram returns false, we have a programming error
+	psAbort ("Unable to generate histogram for psastroZeroPointAnalysis");
+    }
+    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->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 numBins = cumulative->nums->n;
+    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 (%.4f to %.4f)\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);
+
+    return value;
+}
 
 # define ESCAPE(MSG) { \