Index: /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.c	(revision 18350)
@@ -480,5 +480,5 @@
 #endif
 
-bool pmFPAviewFitsWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) {
+bool pmFPAviewFitsWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config) {
 
     bool status = false;
@@ -498,14 +498,14 @@
     switch (file->fileLevel) {
       case PM_FPA_LEVEL_FPA:
-        status = pmFPAWrite(fpa, file->fits, NULL, true, false);
+        status = pmFPAWrite(fpa, file->fits, config, true, false);
         break;
       case PM_FPA_LEVEL_CHIP: {
           pmChip *chip = pmFPAviewThisChip(view, fpa);
-          status = pmChipWrite(chip, file->fits, NULL, true, false);
+          status = pmChipWrite(chip, file->fits, config, true, false);
           break;
       }
       case PM_FPA_LEVEL_CELL: {
           pmCell *cell = pmFPAviewThisCell(view, fpa);
-          status = pmCellWrite(cell, file->fits, NULL, true);
+          status = pmCellWrite(cell, file->fits, config, true);
           break;
       }
Index: /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileFitsIO.h	(revision 18350)
@@ -5,6 +5,6 @@
  * @author PAP, IfA
  *
- * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-06-17 22:16:38 $
+ * @version $Revision: 1.12.4.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-06-27 23:56:41 $
  * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
  */
@@ -70,5 +70,5 @@
 
 /// Write a PHU for a fits image if needed
-bool pmFPAviewFitsWritePHU (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmFPAviewFitsWritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
 
 /// Free the data for the specified view
Index: /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.c	(revision 18350)
@@ -433,5 +433,8 @@
         pmChip *chip = pmFPAviewThisChip(view, file->fpa); // Chip of interest, or NULL
         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest, or NULL
-        pmFPAUpdateNames(fpa, chip, cell);
+        if (!pmFPAUpdateNames(fpa, chip, cell)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to update names in header.");
+            return false;
+        }
 
         pmConceptSource sources = PM_CONCEPT_SOURCE_HEADER | PM_CONCEPT_SOURCE_CELLS |
@@ -918,5 +921,5 @@
 // XXX this function is only called from pmFPAfileWrite
 // XXX for each data type, there should be a function which writes the PHU, if needed
-bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, const pmConfig *config)
+bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, pmConfig *config)
 {
     PS_ASSERT_PTR_NON_NULL(file, false);
Index: /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/camera/pmFPAfileIO.h	(revision 18350)
@@ -5,6 +5,6 @@
  * @author PAP, IfA
  *
- * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-02-07 00:09:40 $
+ * @version $Revision: 1.7.16.1 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-06-27 23:56:41 $
  * Copyright 2004-2005 Institute for Astronomy, University of Hawaii
  */
@@ -49,5 +49,5 @@
 bool pmFPAfileIOChecks (pmConfig *config, const pmFPAview *view, pmFPAfilePlace place);
 
-bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, const pmConfig *config);
+bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, pmConfig *config);
 bool pmFPAfileReadPHU (pmFPAfile *file, const pmFPAview *view, pmConfig *config);
 
Index: /branches/eam_branch_20080620/psModules/src/camera/pmReadoutFake.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/camera/pmReadoutFake.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/camera/pmReadoutFake.c	(revision 18350)
@@ -26,4 +26,6 @@
 #define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
 #define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
+#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
 
 
@@ -92,4 +94,7 @@
     for (int i = 0; i < numSources; i++) {
         pmSource *source = sources->data[i]; // Source of interest
+        if (source->mode & SOURCE_MASK) {
+            continue;
+        }
         if (!isfinite(source->psfMag)) {
             continue;
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtraction.c	(revision 18350)
@@ -3,7 +3,4 @@
  *  @author Paul Price, IfA
  *  @author GLG, MHPCC
- *
- *  @version $Revision: 1.95 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-06-19 03:03:53 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -512,5 +509,9 @@
     PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
 
-    double totalSquareDev = 0.0;        // Total square deviation from zero
+    // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
+    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
+    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
+
+    double mean = 0.0;                  // Mean deviation
     int numStamps = 0;                  // Number of used stamps
     for (int i = 0; i < stamps->num; i++) {
@@ -519,9 +520,18 @@
             continue;
         }
-        totalSquareDev += PS_SQR(deviations->data.F32[i]);
+        mean += deviations->data.F32[i];
         numStamps++;
     }
-
-    float rms = sqrt(totalSquareDev / (double)numStamps); // Convert to RMS
+    mean /= numStamps;
+
+    double rms = 0.0;                   // Standard deviation
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            continue;
+        }
+        rms += PS_SQR(deviations->data.F32[i] - mean);
+    }
+    rms = sqrt(rms / (numStamps - 1));
 
     if (rmsPtr) {
@@ -549,9 +559,15 @@
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
-    double newSquareDev = 0.0;          // New square deviation
+    double newMean = 0.0;               // New mean
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
         if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
-            if (deviations->data.F32[i] > limit) {
+            // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
+            // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
+            // they're anomalously low, compared to the rest of the population which (we hope) is indicative
+            // of normality?  Besides, the standard deviation is going to be blown up by stars that didn't
+            // subtract well, in which case very few (if any) stars will be legitimately rejected for being
+            // low.
+            if (deviations->data.F32[i] - mean > limit) {
                 // Mask out the stamp in the image so you it's not found again
                 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
@@ -587,18 +603,18 @@
             } else {
                 numGood++;
-                newSquareDev += PS_SQR(deviations->data.F32[i]);
+                newMean += deviations->data.F32[i];
             }
         }
     }
+    newMean /= numGood;
 
     if (numRejected > 0) {
         psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
-                 numGood, numRejected, rms,
-                 sqrt(newSquareDev / (double)numGood));
+                 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
+                 numGood, numRejected, mean, newMean);
     } else {
         psLogMsg("psModules.imcombine", PS_LOG_INFO,
-                 "%d good stamps; 0 rejected.\nRMS deviation: %f\n",
-                 numGood, rms);
+                 "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
+                 numGood, mean);
     }
 
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.c	(revision 18350)
@@ -30,5 +30,5 @@
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
-            sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x];
+            sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x];
         }
     }
@@ -194,5 +194,5 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    sumC += conv->kernel[y][x] / weight->kernel[y][x];
+                    sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
                 }
             }
@@ -217,5 +217,5 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                double invNoise2 = 1.0 / weight->kernel[y][x];
+                double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x];
                 double value = input->kernel[y][x] * invNoise2;
                 sumI += value;
@@ -276,5 +276,5 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                    sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x];
+                sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
             }
         }
@@ -296,5 +296,5 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                float value = target->kernel[y][x] / weight->kernel[y][x];
+                float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x];
                 sumIT += value * input->kernel[y][x];
                 sumT += value;
@@ -365,5 +365,5 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                sumC += conv->kernel[y][x] / weight->kernel[y][x];
+                sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
             }
         }
@@ -382,7 +382,32 @@
 }
 
+
+// Add in penalty term to least-squares vector
+static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term
+                             const pmSubtractionKernels *kernels // Kernel parameters
+    )
+{
+    if (kernels->penalty == 0.0) {
+        return true;
+    }
+
+    psVector *penalties = kernels->penalties; // Penalties for each kernel component
+    int spatialOrder = kernels->spatialOrder; // Order of spatial variations
+    int numKernels = kernels->num; // Number of kernel components
+    for (int i = 0; i < numKernels; i++) {
+        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
+            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
+                vector->data.F64[index] -= penalties->data.F32[i];
+            }
+        }
+    }
+
+    return true;
+}
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Semi-public functions
-// XXX EAM these cannot be inline :: talk to Josh
+// XXX We might like to define these functions as "extern inline" but gcc currently doesn't handle this in c99
+// mode.  See http://gcc.gnu.org/ml/gcc/2006-11/msg00006.html
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -675,4 +700,5 @@
         psVectorInit(sumVector, 0.0);
         psImageInit(sumMatrix, 0.0);
+        int numStamps = 0;              // Number of good stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -680,6 +706,8 @@
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
-            }
-        }
+                numStamps++;
+            }
+        }
+        calculatePenalty(sumVector, kernels);
 
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
@@ -716,4 +744,5 @@
         psVectorInit(sumVector2, 0.0);
 
+        int numStamps = 0;              // Number of good stamps
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -724,6 +753,9 @@
                 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
                 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
-            }
-        }
+                numStamps++;
+            }
+        }
+        calculatePenalty(sumVector1, kernels);
+        calculatePenalty(sumVector2, kernels);
 
 #if 0
@@ -778,5 +810,5 @@
         psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A);
         assert(F->numRows == numParams && F->numCols == numParams);
-        float det = NAN;
+        float det = 0.0;
         psImage *Fi = psMatrixInvert(NULL, F, &det);
         assert(Fi->numRows == numParams && Fi->numCols == numParams);
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionEquation.h	(revision 18350)
@@ -21,25 +21,25 @@
 
 /// Calculate the value of a polynomial, specified by coefficients and polynomial values
-inline double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients
-                                                 const psImage *polyValues, ///< Polynomial values
-                                                 int order, ///< Order of polynomials
-                                                 int index, ///< Index at which to begin
-                                                 int step ///< Step between subsequent indices
+double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients
+                                          const psImage *polyValues, ///< Polynomial values
+                                          int order, ///< Order of polynomials
+                                          int index, ///< Index at which to begin
+                                          int step ///< Step between subsequent indices
     );
 
 /// Return the specified coefficient in the solution
-inline double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                           const psImage *polyValues, ///< Polynomial values
-                                           int index, ///< Coefficient index to calculate
-                                           bool wantDual ///< Calculate the coefficient for the dual solution?
+double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                    const psImage *polyValues, ///< Polynomial values
+                                    int index, ///< Coefficient index to calculate
+                                    bool wantDual ///< Calculate the coefficient for the dual solution?
     );
 
 /// Return the normalisation in the solution
-inline double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters
+double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters
     );
 
 /// Return the background (difference) in the solution
-inline double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                                const psImage *polyValues ///< Polynomial values
+double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                         const psImage *polyValues ///< Polynomial values
     );
 
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.c	(revision 18350)
@@ -24,4 +24,5 @@
     psFree(kernels->vStop);
     psFree(kernels->preCalc);
+    psFree(kernels->penalties);
     psFree(kernels->solution1);
     psFree(kernels->solution2);
@@ -60,4 +61,5 @@
     kernels->v = psVectorRealloc(kernels->v, start + numNew);
     kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew);
+    kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     kernels->inner = start;
 
@@ -74,4 +76,5 @@
             kernels->v->data.S32[index] = v;
             kernels->preCalc->data[index] = NULL;
+            kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v));
 
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
@@ -84,5 +87,5 @@
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
                                                     const psVector *fwhms, const psVector *orders,
-                                                    pmSubtractionMode mode)
+                                                    float penalty, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
@@ -104,7 +107,7 @@
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS,
-                                                              size, spatialOrder, mode); // The kernels
-    psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
+                                                              spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements",
@@ -122,8 +125,10 @@
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
                 double sum = 0.0;       // Normalisation
+                double moment = 0.0;    // Moment, for penalty
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
                         sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
+                        moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v));
                     }
                 }
@@ -146,7 +151,8 @@
                 }
                 kernels->preCalc->data[index] = preCalc;
-
-                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index,
-                        fwhms->data.F32[i], uOrder, vOrder);
+                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
+                        fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
             }
         }
@@ -161,5 +167,6 @@
 
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
-                                                int size, int spatialOrder, pmSubtractionMode mode)
+                                                int size, int spatialOrder, float penalty,
+                                                pmSubtractionMode mode)
 {
     pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
@@ -173,4 +180,6 @@
     kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     kernels->preCalc = psArrayAlloc(numBasisFunctions);
+    kernels->penalty = penalty;
+    kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32);
     kernels->uStop = NULL;
     kernels->vStop = NULL;
@@ -188,5 +197,6 @@
 }
 
-pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode)
+pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
+                                               pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -195,7 +205,7 @@
     int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
-                                                              size, spatialOrder, mode); // The kernels
-    psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size,
+                                                              spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
              size, spatialOrder, num);
@@ -211,8 +221,8 @@
 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
                                                const psVector *fwhms, const psVector *orders,
-                                               pmSubtractionMode mode)
-{
-    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders, mode); // Kernels
+                                               float penalty, pmSubtractionMode mode)
+{
+    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
+                                                                  penalty, mode); // Kernels
     if (!kernels) {
         return NULL;
@@ -242,5 +252,5 @@
 
 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
-                                               pmSubtractionMode mode)
+                                               float penalty, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -262,7 +272,8 @@
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
-                                                              size, spatialOrder, mode); // The kernels
-    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
+                                                              spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
+                   penalty);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements",
@@ -324,9 +335,13 @@
     psFree(widths);
 
+    psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels.");
+    psVectorInit(kernels->penalties, 0.0);
+
     return kernels;
 }
 
 
-pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode)
+pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
+                                                pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -354,7 +369,7 @@
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
-                                                              size, spatialOrder, mode); // The kernels
-    psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
+                                                              spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements",
@@ -414,4 +429,7 @@
     psFree(stop);
 
+    psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels.");
+    psVectorInit(kernels->penalties, 0.0);
+
     return kernels;
 }
@@ -419,5 +437,6 @@
 // Grid United with Normal Kernel
 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
-                                               const psVector *orders, int inner, pmSubtractionMode mode)
+                                               const psVector *orders, int inner, float penalty,
+                                               pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -431,6 +450,8 @@
     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
 
-    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders, mode); // Kernels
+    // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed
+
+    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
+                                                                  penalty, mode); // Kernels
     psStringPrepend(&kernels->description, "GUNK=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
@@ -447,5 +468,5 @@
 // RINGS --- just what it says
 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
-                                                pmSubtractionMode mode)
+                                                float penalty, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -477,7 +498,8 @@
     int num = numRings * numPoly; // Total number of basis functions
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS,
-                                                              size, spatialOrder, mode); // The kernels
-    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder);
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
+                                                              spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
+                   penalty);
 
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements",
@@ -520,4 +542,5 @@
                 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords
                 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial
+                double moment = 0.0;    // Moment, for penalty
 
                 if (i == 0) {
@@ -527,4 +550,5 @@
                     uCoords->n = vCoords->n = poly->n = 1;
                     radiusLast = 0;
+                    moment = 0.0;
                 } else {
                     int j = 0;          // Index for data
@@ -546,4 +570,5 @@
                                     poly->data.F32[j] = polyVal;
                                     norm += polyVal;
+                                    moment += polyVal * (PS_SQR(u) + PS_SQR(v));
 
                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
@@ -571,4 +596,5 @@
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                     }
+//                    moment /= norm;
                 }
 
@@ -578,4 +604,5 @@
                 kernels->u->data.S32[index] = uOrder;
                 kernels->v->data.S32[index] = vOrder;
+                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
 
                 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,
@@ -590,19 +617,20 @@
 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
                                                    const psVector *fwhms, const psVector *orders, int inner,
-                                                   int binning, int ringsOrder, pmSubtractionMode mode)
+                                                   int binning, int ringsOrder, float penalty,
+                                                   pmSubtractionMode mode)
 {
     switch (type) {
       case PM_SUBTRACTION_KERNEL_POIS:
-        return pmSubtractionKernelsPOIS(size, spatialOrder, mode);
+        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
       case PM_SUBTRACTION_KERNEL_ISIS:
-        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode);
+        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
       case PM_SUBTRACTION_KERNEL_SPAM:
-        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode);
+        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
       case PM_SUBTRACTION_KERNEL_FRIES:
-        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode);
+        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
       case PM_SUBTRACTION_KERNEL_GUNK:
-        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode);
+        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
       case PM_SUBTRACTION_KERNEL_RINGS:
-        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode);
+        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
@@ -657,4 +685,5 @@
     int binning = 0;                    // Binning to use
     int ringsOrder = 0;                 // Polynomial order for rings
+    float penalty = 0.0;                // Penalty for wideness
 
     if (strncmp(description, "ISIS", 4) == 0) {
@@ -684,5 +713,6 @@
 
             ptr++;                      // Eat ','
-            spatialOrder = parseStringInt(ptr);
+            PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
+            penalty = parseStringFloat(ptr);
         }
     } else if (strncmp(description, "RINGS", 5) == 0) {
@@ -692,5 +722,6 @@
         PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt);
-        PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt);
+        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
+        PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
     } else {
         psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported.");
@@ -699,5 +730,5 @@
 
     return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders,
-                                        inner, binning, ringsOrder, mode);
+                                        inner, binning, ringsOrder, penalty, mode);
 }
 
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionKernels.h	(revision 18350)
@@ -33,4 +33,6 @@
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
+    float penalty;                      ///< Penalty for wideness
+    psVector *penalties;                ///< Penalty for each kernel component
     int size;                           ///< The half-size of the kernel
     int inner;                          ///< The size of an inner region
@@ -107,4 +109,5 @@
                                                 int size, ///< Half-size of kernel
                                                 int spatialOrder, ///< Order of spatial variations
+                                                float penalty, ///< Penalty for wideness
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -113,4 +116,5 @@
 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
                                                int spatialOrder, ///< Order of spatial variations
+                                               float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -121,4 +125,5 @@
                                                     const psVector *fwhms, ///< Gaussian FWHMs
                                                     const psVector *orders, ///< Polynomial order of gaussians
+                                                    float penalty, ///< Penalty for wideness
                                                     pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -129,4 +134,5 @@
                                                const psVector *fwhms, ///< Gaussian FWHMs
                                                const psVector *orders, ///< Polynomial order of gaussians
+                                               float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
                                                );
@@ -137,4 +143,5 @@
                                                int inner, ///< Inner radius to preserve unbinned
                                                int binning, ///< Kernel binning factor
+                                               float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -144,4 +151,5 @@
                                                 int spatialOrder, ///< Order of spatial variations
                                                 int inner, ///< Inner radius to preserve unbinned
+                                                float penalty, ///< Penalty for wideness
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -153,4 +161,5 @@
                                                const psVector *orders, ///< Polynomial order of gaussians
                                                int inner, ///< Inner radius containing grid of delta functions
+                                               float penalty, ///< Penalty for wideness
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -161,4 +170,5 @@
                                                 int inner, ///< Inner radius to preserve unbinned
                                                 int ringsOrder, ///< Polynomial order
+                                                float penalty, ///< Penalty for wideness
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -174,4 +184,5 @@
                                                    int binning, ///< Kernel binning factor
                                                    int ringsOrder, ///< Polynomial order for RINGS
+                                                   float penalty, ///< Penalty for wideness
                                                    pmSubtractionMode mode ///< Mode for subtraction
     );
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.c	(revision 18350)
@@ -117,7 +117,8 @@
                         pmSubtractionKernelsType type, int size, int spatialOrder,
                         const psVector *isisWidths, const psVector *isisOrders,
-                        int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
-                        int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
-                        psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
+                        int inner, int ringsOrder, int binning, float penalty,
+                        bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold,
+                        int iter, float rej, psMaskType maskBad, psMaskType maskBlank,
+                        float badFrac, pmSubtractionMode mode)
 {
     if (mode != PM_SUBTRACTION_MODE_2) {
@@ -320,5 +321,5 @@
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold, mode);
+                                                          stamps, footprint, optThreshold, penalty, mode);
                 if (!kernels) {
                     psErrorClear();
@@ -329,5 +330,5 @@
                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder, mode);
+                                                       inner, binning, ringsOrder, penalty, mode);
             }
 
@@ -450,5 +451,9 @@
                 // Generate image with convolution kernels
                 int fullSize = 2 * size + 1 + 1;    // Full size of kernel
-                psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
+                int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
+                psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
+                                                    imageSize - 1 +
+                                                    (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
+                                                    imageSize - 1, PS_TYPE_F32);
                 psImageInit(convKernels, NAN);
                 for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
@@ -471,4 +476,27 @@
                         }
                         psFree(kernel);
+
+                        if (mode == PM_SUBTRACTION_MODE_DUAL) {
+                            kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
+                                                              (float)j / (float)KERNEL_MOSAIC,
+                                                              true); // Image of the kernel
+                            if (!kernel) {
+                                psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                                psFree(convKernels);
+                                goto MATCH_ERROR;
+                            }
+
+                            if (psImageOverlaySection(convKernels, kernel,
+                                                      (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize +
+                                                      4,
+                                                      (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
+                                psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                                psFree(kernel);
+                                psFree(convKernels);
+                                goto MATCH_ERROR;
+                            }
+                            psFree(kernel);
+                        }
+
                     }
                 }
@@ -509,4 +537,26 @@
                 goto MATCH_ERROR;
             }
+
+            // Set the variance factors
+            switch (mode) {
+              case PM_SUBTRACTION_MODE_1: {
+                  recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
+                                       region);
+                  break;
+              }
+              case PM_SUBTRACTION_MODE_2:
+                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
+                                     region);
+                break;
+              case PM_SUBTRACTION_MODE_DUAL:
+                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false),
+                                     region);
+                recordVarianceFactor(conv2, pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true),
+                                     region);
+                break;
+              default:
+                psAbort("Should never reach here.");
+            }
+
             psFree(kernels);
             kernels = NULL;
@@ -526,28 +576,4 @@
                     conv2->parent->parent->data_exists = true;
                 }
-            }
-
-
-            // Set the variance factors
-            float xMid = 0.5 * (region ? (region->x0 + region->x1) : numCols); // Average x position
-            float yMid = 0.5 * (region ? (region->y0 + region->y1) : numRows); // Average y position
-            switch (mode) {
-              case PM_SUBTRACTION_MODE_1: {
-                  recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, xMid, yMid, false),
-                                       region);
-                  break;
-              }
-              case PM_SUBTRACTION_MODE_2:
-                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, xMid, yMid, false),
-                                     region);
-                break;
-              case PM_SUBTRACTION_MODE_DUAL:
-                recordVarianceFactor(conv1, pmSubtractionVarianceFactor(kernels, xMid, yMid, false),
-                                     region);
-                recordVarianceFactor(conv2, pmSubtractionVarianceFactor(kernels, xMid, yMid, true),
-                                     region);
-                break;
-              default:
-                psAbort("Should never reach here.");
             }
         }
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionMatch.h	(revision 18350)
@@ -37,4 +37,5 @@
                         int ringsOrder, ///< RINGS polynomial order
                         int binning,    ///< SPAM kernel binning
+                        float penalty,  ///< Penalty for wideness
                         bool optimum,   ///< Search for optimum ISIS kernel?
                         const psVector *optFWHMs, ///< FWHMs for optimum search
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.c
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.c	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.c	(revision 18350)
@@ -204,5 +204,5 @@
                                                       int spatialOrder, const psVector *fwhms, int maxOrder,
                                                       const pmSubtractionStampList *stamps, int footprint,
-                                                      float tolerance, pmSubtractionMode mode)
+                                                      float tolerance, float penalty, pmSubtractionMode mode)
 {
     if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
@@ -231,6 +231,6 @@
     psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders
     psVectorInit(orders, maxOrder);
-    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders, mode); // Kernels
+    pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
+                                                                  penalty, mode); // Kernels
     psFree(orders);
     psFree(kernels->description);
Index: /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.h
===================================================================
--- /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.h	(revision 18349)
+++ /branches/eam_branch_20080620/psModules/src/imcombine/pmSubtractionParams.h	(revision 18350)
@@ -16,4 +16,5 @@
                                                       int footprint, ///< Convolution footprint for stamps
                                                       float tolerance, ///< Maximum difference in chi^2
+                                                      float penalty, ///< Penalty for wideness
                                                       pmSubtractionMode mode // Mode for subtraction
     );