Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 15426)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 15427)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.70 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-10-31 21:09:33 $
+ *  @version $Revision: 1.71 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-11-01 01:11:03 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -688,11 +688,24 @@
         psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
 
+#ifdef TESTING
+        for (int j = 0; j < numKernels; j++) {
+            psKernel *conv = convolutions->data[j]; // Convolution of interest
+            psString filename = NULL; // Filename for output
+            psStringAppend(&filename, "conv_%03d_%03d.fits", i, j);
+            psFits *fits = psFitsOpen(filename, "w"); // FITS file pointer
+            psFree(filename);
+            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+            psFitsClose(fits);
+        }
+#endif
+
         // Generate least-squares vector and matrix
-        bool bad = false;           // Are there bad values?
-        for (int j = 0; j < numKernels; j++) {
+        bool bad = false;               // Are there bad values?
+        psString badString = NULL;      // Details on bad values
+        for (int j = 0; j < numKernels && !bad; j++) {
             psKernel *jConv = convolutions->data[j]; // Convolution for j-th element
 
             // Generate matrix
-            for (int k = j; k < numKernels; k++) {
+            for (int k = j; k < numKernels && !bad; k++) {
                 psKernel *kConv = convolutions->data[k]; // Convolution for k-th element
                 double sumCC = 0.0; // Sum of the convolution products
@@ -703,5 +716,7 @@
                 }
                 if (!isfinite(sumCC)) {
+                    psStringAppend(&badString, "\nBad sumCC at %d, %d", j, k);
                     bad = true;
+                    continue;
                 }
                 // Generate the pseudo-convolutions from the spatial polynomial terms
@@ -732,7 +747,15 @@
                 }
             }
-            if (!isfinite(sumC) || !isfinite(sumIC)) {
+            if (!isfinite(sumC)) {
+                psStringAppend(&badString, "\nBad sumC at %d", j);
                 bad = true;
-            }
+                continue;
+            }
+            if (!isfinite(sumIC)) {
+                psStringAppend(&badString, "\nBad sumIC detected at %d", j);
+                bad = true;
+                continue;
+            }
+
             for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
                 for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
@@ -759,7 +782,12 @@
             }
         }
-        if (!isfinite(sum1) || !isfinite(sumI)) {
+        if (!isfinite(sum1)) {
+            psStringAppend(&badString, "\nBad sum1 detected");
             bad = true;
-        }
+        } else if (!isfinite(sumI)) {
+            psStringAppend(&badString, "\nBad sumI detected");
+            bad = true;
+        }
+
         stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
         stampVector->data.F64[bgIndex] = sumI;
@@ -767,6 +795,7 @@
         if (bad) {
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
-                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+            psWarning("Rejecting stamp %d (%d,%d) because of bad equation:%s",
+                      i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), badString);
+            psFree(badString);
         } else {
             stamp->status = PM_SUBTRACTION_STAMP_USED;
@@ -854,5 +883,5 @@
     if (psTraceGetLevel("psModules.imcombine") >= 7) {
         for (int i = 0; i < solution->n; i++) {
-            printf("       Solution %d: %f\n", i, solution->data.F64[i]);
+            psTrace("psModules.imcombine", 7, "Solution %d: %f\n", i, solution->data.F64[i]);
         }
     }