Index: /tags/ipp-20130712/ippconfig/gpc1/ppStack.config
===================================================================
--- /tags/ipp-20130712/ippconfig/gpc1/ppStack.config	(revision 35940)
+++ /tags/ipp-20130712/ippconfig/gpc1/ppStack.config	(revision 35941)
@@ -63,4 +63,5 @@
     PSF.INPUT.THRESH        F32   10.0
     PSF.INPUT.ASYMMETRY     F32   0.2
+    PSF.TARGET.AS.MAX	    BOOL  TRUE
    MATCH.REJ        F32    4.0
    SAFE		 BOOL	F
Index: /tags/ipp-20130712/ippconfig/recipes/ppStack.config
===================================================================
--- /tags/ipp-20130712/ippconfig/recipes/ppStack.config	(revision 35940)
+++ /tags/ipp-20130712/ippconfig/recipes/ppStack.config	(revision 35941)
@@ -95,4 +95,7 @@
 PSF.INPUT.THRESH        F32  NAN	# Set minimum limit below which we do not exclude an input (defaults to 0.0)
 PSF.INPUT.ASYMMETRY     F32  NAN	# Set difference in mixture model populations to consider equal.
+PSF.TARGET.AS.MAX       BOOL F          # Set the target PSF FWHM as the maximum of accepted input FWHM values.
+PSF.TARGET.AS.MAX.EPSILON F32 0.1       # Amount to set the target PSF FWHM larger than the maximum input. Target = eps + max(input)
+
 
 TEMP.IMAGE	STR	conv.im.fits	# Suffix for temporary convolved images
Index: /tags/ipp-20130712/ppStack/src/ppStackPrepare.c
===================================================================
--- /tags/ipp-20130712/ppStack/src/ppStackPrepare.c	(revision 35940)
+++ /tags/ipp-20130712/ppStack/src/ppStackPrepare.c	(revision 35941)
@@ -257,9 +257,10 @@
 
     bool mdok = false;
-    bool  simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
+    bool  simpleClip = psMetadataLookupBool(&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
@@ -477,6 +478,21 @@
                          "Target PSF for stack", options->psf);
         options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target
-        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);
-
+
+
+	bool psfTargetAsMax = psMetadataLookupBool(&mdok, recipe, "PSF.TARGET.AS.MAX");
+	if (psfTargetAsMax) { // Should we use the largest input as the target?
+	  float psfTargetEpsilon = psMetadataLookupF32(&mdok,recipe,"PSF.TARGET.AS.MAX.EPSILON");
+	  if (!mdok) { psfTargetEpsilon = 0.0; }
+	  options->targetSeeing = 0.0;
+	  for (int i = 0; i < num; i++) {
+            if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+	    options->targetSeeing = PS_MAX(options->targetSeeing,options->inputSeeing->data.F32[i]);
+	  }
+	  psLogMsg("ppStack", PS_LOG_INFO, "Using MAX accepted input FWHM as target (max: %f epsilon: %f target %f)\n",
+		   options->targetSeeing,psfTargetEpsilon,options->targetSeeing + psfTargetEpsilon);
+	  options->targetSeeing = options->targetSeeing + psfTargetEpsilon;
+	}
+        psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing);	
+	
         pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip
         psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
Index: /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionMatch.c	(revision 35940)
+++ /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionMatch.c	(revision 35941)
@@ -540,5 +540,5 @@
     // Bail here if we're doing the simple matching
     if (type == PM_SUBTRACTION_KERNEL_SIMPLE) {
-      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor)) {
+      if (!pmSubtractionSimpleMatch(conv1,conv2,ro1,ro2,sources,size,maskVal,maskBad,maskPoor,optThreshold)) {
 	return false;
       }
Index: /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.c
===================================================================
--- /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35940)
+++ /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.c	(revision 35941)
@@ -63,4 +63,20 @@
 }  
 
+bool simple_apply_mask(psImage *image, psImage *weight, psImage *mask,
+		       psImageMaskType maskVal) {
+  for (int y = 0; y < mask->numRows; y++) {
+    for (int x = 0; x < mask->numCols; x++) {
+      if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+	image->data.F32[y][x] = NAN;
+	if (weight) {
+	  weight->data.F32[y][x] = NAN;
+	}
+      }
+    }
+  }
+  return(true);
+}
+		       
+
 // Copied from pmSubtraction
 static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated
@@ -90,5 +106,6 @@
 			      psImageMaskType maskVal,
 			      psImageMaskType maskBad,
-			      psImageMaskType maskPoor
+			      psImageMaskType maskPoor,
+			      float deconvolveThreshold
 			      ) {
   //
@@ -115,4 +132,6 @@
   psImage *varC2 = NULL;
 
+  psImage *maskTemp = NULL;
+  
   // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called.  
   int numCols = ro1->image->numCols;
@@ -190,13 +209,17 @@
   if (!conv1) {
     if (convolution_direction == 1) {
-      chisq = 100;
-    }
-    convolution_direction = 2;
+      if (sigma1 - sigma2 > deconvolveThreshold) {
+	chisq = 100;
+      }
+    }
+    //    convolution_direction = 2;
   }
   if (!conv2) {
     if (convolution_direction == 2) {
-      chisq = 100;
-    }
-    convolution_direction = 1;
+      if (sigma2 - sigma1 > deconvolveThreshold) {
+	chisq = 100;
+      }
+    }
+    //    convolution_direction = 1;
   }
   
@@ -205,5 +228,5 @@
   int maskBox = (int) ceil(sigmaKern * 1.1774); // diameter is 1/2 FWHM
   int maskBlank = 8;  // I should be able to get this from a reference, right?
-
+  int maskPoorVal = 16384; // Another value that should be found elsewhere.
   //
   // Make a fake pmSubtractionKernels element so we can add it appropriately.
@@ -249,9 +272,19 @@
   // Do convolutions
   if (convolution_direction == 1) {
-    psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
-    psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
-    maskC1 = psImageConvolveMask(maskC1,mask1,maskVal,maskBad,
-				 -maskBox,maskBox,-maskBox,maskBox);
-    conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance);
+    if (conv1) {
+      psImageSmoothMask_Threaded(imageC1,image1,mask1,maskVal,sigmaKern,6,1e-6);
+      psImageSmoothMask_Threaded(varC1,var1,mask1,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+
+      maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+      maskTemp = psImageConvolveMask(maskTemp,mask1,maskVal,maskBad,      // Mask bad values
+				   -maskBox,maskBox,-maskBox,maskBox);
+      maskC1 = psImageConvolveMask(maskC1,maskTemp,maskPoorVal,maskPoor,  // Mask poor values
+				   -maskBox,maskBox,-maskBox,maskBox);
+      psFree(maskTemp);
+
+      conv1->covariance = psImageCovarianceCalculate(kernel,ro1->covariance);
+      pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
+      simple_apply_mask(imageC1,varC1,maskC1,maskBad);
+    }
     if (conv2) {
       imageC2 = psImageCopy(imageC2,image2,PS_TYPE_F32);
@@ -260,13 +293,21 @@
       conv2->covariance = psMemIncrRefCounter(ro2->covariance);
     }
-    pmSubtractionBorder(imageC1,varC1,maskC1,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC1,varC1,maskC1,PM_SUBTRACTION_MODE_1);
   }
   else if (convolution_direction == 2) {
-    psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
-    psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
-    maskC2 = psImageConvolveMask(maskC2,mask2,maskVal,maskBad,
-				 -maskBox,maskBox,-maskBox,maskBox);
-    conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance);
+    if (conv2) {
+      psImageSmoothMask_Threaded(imageC2,image2,mask2,maskVal,sigmaKern,6,1e-6);
+      psImageSmoothMask_Threaded(varC2,var2,mask2,maskVal,sigmaKern * M_SQRT1_2,6,1e-6);
+
+      maskTemp = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+      maskTemp = psImageConvolveMask(maskTemp,mask2,maskVal,maskBad,      // Mask bad values
+				   -maskBox,maskBox,-maskBox,maskBox);
+      maskC2 = psImageConvolveMask(maskC2,maskTemp,maskPoorVal,maskPoor,  // Mask poor values
+				   -maskBox,maskBox,-maskBox,maskBox);
+      psFree(maskTemp);
+
+      conv2->covariance = psImageCovarianceCalculate(kernel,ro2->covariance);
+      pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
+      simple_apply_mask(imageC2,varC2,maskC2,maskBad);
+    }
     if (conv1) {
       imageC1 = psImageCopy(imageC1,image1,PS_TYPE_F32);
@@ -275,6 +316,4 @@
       conv1->covariance = psMemIncrRefCounter(ro1->covariance);
     }
-    pmSubtractionBorder(imageC2,varC2,maskC2,maskBox,maskBlank);
-    pmSubtractionMaskApply(imageC2,varC2,maskC2,PM_SUBTRACTION_MODE_2);
   }    
 
@@ -294,22 +333,18 @@
     float flux1,flux2;
 
-    if (convolution_direction == 1) {
+    if (conv1) {
       simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
-      if (conv2) {
-	simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
-      }
-      else {
-	simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
-      }
-    }
-    else if (convolution_direction == 2) {
+    }
+    else {
+      simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
+    }
+
+    if (conv2) {
       simple_do_boxphot(&nPix2,&flux2,source,imageC2,maskC2,maskBad,photRadius);
-      if (conv1) {
-	simple_do_boxphot(&nPix1,&flux1,source,imageC1,maskC1,maskBad,photRadius);
-      }
-      else {
-	simple_do_boxphot(&nPix1,&flux1,source,image1,mask1,maskBad,photRadius);
-      }
-    }
+    }
+    else {
+      simple_do_boxphot(&nPix2,&flux2,source,image2,mask2,maskBad,photRadius);
+    }
+
     logFluxDifferences->data.F32[i] = flux2 - flux1;
     fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
@@ -321,5 +356,4 @@
     //    fprintf(stderr,"SOURCES: %d %g %g %g -> %d %d %g %g %d %g\n",i,source->peak->xf,source->peak->yf,source->psfMag,
     //	    nPix1,nPix2,flux1,flux2,fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i],logFluxDifferences->data.F32[i]);
-    
   }
 
Index: /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.h
===================================================================
--- /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.h	(revision 35940)
+++ /tags/ipp-20130712/psModules/src/imcombine/pmSubtractionSimple.h	(revision 35941)
@@ -16,5 +16,6 @@
 			      psImageMaskType maskVal,
 			      psImageMaskType maskBad,
-			      psImageMaskType maskPoor
+			      psImageMaskType maskPoor,
+			      float deconvolveThreshold
 			      );