Index: trunk/ppStack/src/ppStackPrepare.c
===================================================================
--- trunk/ppStack/src/ppStackPrepare.c	(revision 35455)
+++ trunk/ppStack/src/ppStackPrepare.c	(revision 35529)
@@ -257,8 +257,10 @@
 
     bool mdok = false;
+    bool  simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
     float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
     float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
     float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below
     float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry
+    
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
     bool havePSFs = false;                                // Do we have any PSFs?
@@ -309,52 +311,104 @@
 
 #ifdef NEW_REJECTION
-    // Do GMM test with two modes to decide where to put the break.
-    double Punimodal;
-    int m = 2;
-    psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
-    psVector *means = psVectorAlloc(m,PS_TYPE_F32);
-    psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
-    psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
-    psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
-    float limit;
-    if (!psMM1DClass(options->inputSeeing,
-		     options->inputSeeing->n,
-		     modes,means,
-		     S,pi,P,
-		     2,
-		     &Punimodal)) {
-      // Handle error here
+    // Begin new rejection code.
+    float limit = INFINITY; // Default to a big value.
+    if (!isfinite(threshFWHM)) { // If this isn't set, initialize it to zero.  No inputs may be clipped if smaller than this
+      threshFWHM = 0.0;
+    }
+    if (!isfinite(maxFWHM)) { // No maxFWHM was set, accept all inputs
+      maxFWHM = INFINITY;
+    }
+    limit = maxFWHM; // Initialize the limit to be the maxFWHM
+    if (limit < threshFWHM) {
+      limit = threshFWHM;
+    }
+    if (!isfinite(clipFWHMnSig)) { // We cannot do a sigma clip with nSig, so use the maxFWHM to set the limit.
+      psLogMsg("ppStack",PS_LOG_INFO,
+	       "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
+	       limit,
+	       clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
+    }
+    if (simpleClip) { // Do a sigma clip like the old rejection code.
+      psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+      psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
+      psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
+      
+      options->clippedMean =  fwhmStats->clippedMean;
+      options->clippedStdev = fwhmStats->clippedStdev;
+      limit = options->clippedMean + clipFWHMnSig * options->clippedStdev;
+      if (limit < threshFWHM) {
+	limit = threshFWHM;
+      }
+      if (limit > maxFWHM) {
+	limit = maxFWHM;
+      }
+      psLogMsg("ppStack",PS_LOG_INFO,
+	       "PSF FWHM distribution: Used simple clip method.  Limit set to %f.",
+	       limit);
+      psFree(fwhmStats);
+    }
+    else { // Do the mixture model rejection.
+      if(!isfinite(asymmetryFWHM)) { // Or not, because the parameters aren't set.
+	psLogMsg("ppStack",PS_LOG_INFO,
+		 "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
+		 limit,
+		 clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
+      }
+      // Do GMM test with two modes to decide where to put the break.
+      double Punimodal;
+      int m = 2;
+      psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
+      psVector *means = psVectorAlloc(m,PS_TYPE_F32);
+      psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
+      psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
+      psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
+      
+      if (!psMM1DClass(options->inputSeeing,
+		       options->inputSeeing->n,
+		       modes,means,
+		       S,pi,P,
+		       2,
+		       &Punimodal)) {
+	// Handle error here
+	psFree(modes);
+	psFree(means);
+	psFree(S);
+	psFree(pi);
+	psFree(P);
+	return(false);
+      }
+      //    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
+      //    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
+      //    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
+      
+      // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
+      if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
+	  (num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
+	limit = maxFWHM;
+      }
+      else {                  // This is a bimodal distribution
+	if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
+	    (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
+	  limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
+	}
+	else {                                   // The smaller FWHM mode is more populated
+	  limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
+	}
+      }
+      if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
+      if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
+      psLogMsg("ppStack",PS_LOG_INFO,
+	       "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
+	       limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
+	       means->data.F32[1],S->data.F32[1],pi->data.F32[1],
+	       Punimodal);
+      psFree(means);
       psFree(modes);
-      psFree(means);
       psFree(S);
       psFree(pi);
       psFree(P);
-      return(false);
-    }
-    //    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
-    //    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
-    //    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
-
-    // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
-    if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
-	(num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
-      limit = maxFWHM;
-    }
-    else {                  // This is a bimodal distribution
-      if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
-	  (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
-	limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
-      }
-      else {                                   // The smaller FWHM mode is more populated
-	limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
-      }
-    }
-    if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
-    if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
-    psLogMsg("ppStack",PS_LOG_INFO,
-	     "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
-	     limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
-	     means->data.F32[1],S->data.F32[1],pi->data.F32[1],
-	     Punimodal);
+    } // End mixture model case
+
+    // Perform rejection using the limit set
     for (int i = 0; i < num; i++) {
       if (options->inputSeeing->data.F32[i] > limit) {
@@ -365,9 +419,4 @@
       }
     }
-    psFree(means);
-    psFree(modes);
-    psFree(S);
-    psFree(pi);
-    psFree(P);
     // End new rejection code
 #endif