Index: /trunk/psphot/src/psphotSourceSize.c
===================================================================
--- /trunk/psphot/src/psphotSourceSize.c	(revision 20877)
+++ /trunk/psphot/src/psphotSourceSize.c	(revision 20878)
@@ -100,30 +100,26 @@
         // Compare the central pixel with those on either side, for the four possible lines through it.
 
+        // Soften weights (add systematic error)
+        float softening = soft * PS_SQR(source->peak->flux); // Softening for weights
+
         // Across the middle: y = 0
         float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
         float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1];
-        float nX = cX / sqrtf(dcX);
+        float nX = cX / sqrtf(dcX + softening);
 
         // Up the centre: x = 0
         float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
         float dcY = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak+0] + weight[yPeak+1][xPeak+0];
-        float nY = cY / sqrtf(dcY);
+        float nY = cY / sqrtf(dcY + softening);
 
         // Diagonal: x = y
         float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
         float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1];
-        float nL = cL / sqrtf(dcL);
+        float nL = cL / sqrtf(dcL + softening);
 
         // Diagonal: x = - y
         float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
         float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1];
-        float nR = cR / sqrtf(dcR);
-
-        // Soften weights (add systematic error)
-        float softening = 6 * soft * source->peak->flux; // Softening for weights
-        dcX += softening;
-        dcY += softening;
-        dcL += softening;
-        dcR += softening;
+        float nR = cR / sqrtf(dcR + softening);
 
         // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
