Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 29543)
+++ trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 29598)
@@ -25,4 +25,5 @@
 # define PENALTY false
 # define MOMENTS (!PENALTY)
+# define MOMENTS_PENALTY_SCALE 2e-5 // XXX this value is not completely arbitrary, but I don't understand why it needs to be this value...
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -325,28 +326,15 @@
 		    // add in second moments
 		    if (MOMENTS) {
-			matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA;
-			matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA;
-
-			matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA;
-			matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA;
-
-			matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB;
-			matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB;
-
-			matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB;
-			matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB;
-
-			// XXX this uses the single moments, first try
-			// matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
-			// matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
-			// 
-			// matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
-			// matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
-			// 
-			// matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
-			// matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
-			// 
-			// matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
-			// matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
+			matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
+			matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
+
+			matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
+			matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
+
+			matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
+			matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
+
+			matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
+			matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
 		    }
 		    matrix->data.F64[iIndex][jIndex + numParams] = ab;
@@ -465,9 +453,9 @@
 	    // add in second moments
 	    if (MOMENTS) {
-		vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1;
-		vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1;
-
-		vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2;
-		vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2;
+		vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1 * MOMENTS_PENALTY_SCALE;
+		vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1 * MOMENTS_PENALTY_SCALE;
+
+		vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2 * MOMENTS_PENALTY_SCALE;
+		vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2 * MOMENTS_PENALTY_SCALE;
 	    }
         }
@@ -793,4 +781,5 @@
         }
     }
+
     stamp->norm = normI2 / normI1;
     stamp->normI1 = normI1;