Index: /branches/eam_branches/ppStack.20211015/src/Makefile.am
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/Makefile.am	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/Makefile.am	(revision 41842)
@@ -29,4 +29,6 @@
 	ppStackPrepare.c	\
 	ppStackConvolve.c	\
+	ppStackLoopByPercent.c  \
+	ppStackCombinePercent.c \
 	ppStackCombinePrepare.c	\
 	ppStackCombineInitial.c	\
Index: /branches/eam_branches/ppStack.20211015/src/ppStack.h
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStack.h	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStack.h	(revision 41842)
@@ -115,4 +115,16 @@
 bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
     );
+
+// Thread entry point for ppStackReadoutPercent
+bool ppStackReadoutPercentThread(psThreadJob *job);
+
+bool ppStackReadoutPercent(const pmConfig *config,
+			   pmReadout *outRO,
+			   const psArray *readouts,
+			   const psVector *mask,
+			   const psVector *weightings,
+			   const psVector *exposures,
+			   const psVector *addVariance);
+
 // Perform median stacking for background
 bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config);
Index: /branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c	(revision 41842)
@@ -1,6 +1,3 @@
 #include "ppStack.h"
-
-// This is the doomsday switch.
-// #define TESTING                         // Enable test output
 
 bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config)
@@ -40,15 +37,4 @@
     }
 
-    // Harvest the jobs, gathering the inspection lists
-    options->inspect = psArrayAlloc(options->num);
-    options->rejected = psArrayAlloc(options->num);
-    for (int i = 0; i < options->num; i++) {
-        if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
-            continue;
-        }
-        options->inspect->data[i] = psArrayAllocEmpty(numChunk);
-        options->rejected->data[i] = psArrayAllocEmpty(numChunk);
-    }
-
     psThreadJob *job;               // Completed job
     while ((job = psThreadJobGetDone())) {
Index: /branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c	(revision 41842)
@@ -5,5 +5,6 @@
     assert(config);
 
-    // Convolve inputs (this must be called even if the convolution is turned off)
+    // Convolve inputs (this must be called even if the convolution is turned off because
+    // it performs the normalization of the inputs).  
     psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
     if (!ppStackConvolve(options, config)) {
@@ -15,6 +16,4 @@
     if (options->quality) return true; // Can't do anything else
 
-    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
-
     // Define threading elements
     ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true);
@@ -24,5 +23,5 @@
     }
 
-    // Prepare for combination
+    // Prepare for combination (define outputs, generate output cells)
     if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", options->doBackground ? "PPSTACK.OUTPUT.BKGMODEL" : NULL, PPSTACK_FILES_STACK, stack, options, config)) {
         psError(psErrorCodeLast(), false, "Unable to prepare for combination.");
Index: /branches/eam_branches/ppStack.20211015/src/ppStackOptions.c
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackOptions.c	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackOptions.c	(revision 41842)
@@ -58,6 +58,4 @@
 
     options->clipPercent = false;
-    options->clipMinRange = 0.16; // 1-sigma default
-    options->clipMaxRange = 0.84; // 1-sigma default
 
     options->psf = NULL;
Index: /branches/eam_branches/ppStack.20211015/src/ppStackOptions.h
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackOptions.h	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackOptions.h	(revision 41842)
@@ -18,6 +18,4 @@
 
     bool clipPercent;                   // use percentile range to clip?
-    psF32 clipMinRange;                 // min percentile
-    psF32 clipMaxRange;                 // max percentile
 
     // Prepare
Index: /branches/eam_branches/ppStack.20211015/src/ppStackReadout.c
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackReadout.c	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackReadout.c	(revision 41842)
@@ -68,6 +68,5 @@
     psVector *addVariance = options->matchChi2; // Additional variance when rejecting
 
-    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
-                                         weightings, exposures, addVariance);
+    job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask, weightings, exposures, addVariance);
     thread->busy = false;
 
@@ -340,5 +339,7 @@
 }
 
-psArray *ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+bool pmStackCombineByPercentile(pmReadout *outRO, psArray *stack, float minRange, float maxRange, psImageMaskType maskBad, psImageMaskType maskSuspect, psImageMaskType maskBlank);
+
+bool ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
 			       const psVector *mask, const psVector *weightings, const psVector *exposures,
 			       const psVector *addVariance)
@@ -357,30 +358,22 @@
 
     bool mdok;                          // Status of MD lookup
-    float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations
-    float minRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MIN.RANGE"); // Combination threshold
-    float maxRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MAX.RANGE"); // Combination threshold
-
-    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
-    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
-
-    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask for bad
-    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
-
-    psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    float minRange = psMetadataLookupF32(&mdok, recipe, "COMBINE.MIN.RANGE"); // Combination threshold
+    float maxRange = psMetadataLookupF32(&mdok, recipe, "COMBINE.MAX.RANGE"); // Combination threshold
+
+    char defaultBlankStr[16] = "BLANK";
+
+    psString maskBadStr     = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask for bad
+    psString maskSuspectStr = psMetadataLookupStr(&mdok, recipe, "MASK.SUSPECT"); // Name of suspect mask bits
+    psString maskBlankStr   = psMetadataLookupStr(&mdok, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels
+    if (!maskBlankStr) {
+      maskBlankStr = psMetadataLookupStr(&mdok, recipe, "MASK.BAD"); // Old name for MASK.BLANK
+    }
+    if (!maskBlankStr) {
+      maskBlankStr = defaultBlankStr; // this is statically allocated above
+    }
+
+    psImageMaskType maskBad     = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
     psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits
-
-    bool status = false;
-    psImageMaskType maskBlank;
-    psString maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels
-    if (maskBlankStr) {
-      maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
-    } else {
-      maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BAD"); // Old name for MASK.BLANK
-      if (maskBlankStr) {
-	maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
-      } else {
-	maskBlank = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels
-      }
-    }
+    psImageMaskType maskBlank   = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels
 
     int num = readouts->n;              // Number of inputs
@@ -412,10 +405,6 @@
     sectionNum++;
 
-    psArray *results = psArrayAlloc(2); // Array of results
-    results->data[0] = inspect;
-    results->data[1] = reject;
-
-    return results;
-}
-
-
+    return true;
+}
+
+
Index: /branches/eam_branches/ppStack.20211015/src/ppStackSetup.c
===================================================================
--- /branches/eam_branches/ppStack.20211015/src/ppStackSetup.c	(revision 41841)
+++ /branches/eam_branches/ppStack.20211015/src/ppStackSetup.c	(revision 41842)
@@ -20,7 +20,4 @@
 
     options->clipPercent  = psMetadataLookupBool(NULL, recipe, "COMBINE.PERCENT"); // use percentile range to clip?
-    options->clipMinRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MIN.RANGE"); // min percentile
-    options->clipMaxRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MAX.RANGE"); // max percentile
-
 
     if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {