Index: branches/simmosaic_branches/psphot/src/psphotMakeResiduals.c
===================================================================
--- branches/simmosaic_branches/psphot/src/psphotMakeResiduals.c	(revision 24860)
+++ branches/simmosaic_branches/psphot/src/psphotMakeResiduals.c	(revision 27839)
@@ -1,3 +1,5 @@
 # include "psphotInternal.h"
+
+# define RESIDUAL_SOFTENING 0.005 
 
 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal) {
@@ -31,4 +33,7 @@
 
     float pixelSN = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.PIX.SN");
+    PS_ASSERT (status, false);
+
+    float radiusMax = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.RADIUS");
     PS_ASSERT (status, false);
 
@@ -171,5 +176,4 @@
                 bool offImage = false;
                 if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
-                    // fprintf (stderr, "off image: %f %f : %f %f\n", ix, iy, flux, dflux);
                     // This pixel is off the image
                     offImage = true;
@@ -179,6 +183,9 @@
                 }
                 fluxes->data.F32[i] = flux;
-                dfluxes->data.F32[i] = dflux;
+                dfluxes->data.F32[i] = hypot(dflux, RESIDUAL_SOFTENING);
                 if (isnan(flux)) {
+                    fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
+                }
+                if (isnan(dflux)) {
                     fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i] = badMask;
                 }
@@ -234,4 +241,10 @@
 		}
 
+		float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+		if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+		  continue;
+                }
+
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
                 resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
@@ -248,9 +261,13 @@
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
                 }
-
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n", ox, oy, resid->Ro->data.F32[oy][ox], fluxStats->sampleStdev, fluxStats->sampleStdev/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
-
             } else {
                 assert (SPATIAL_ORDER == 1);
+
+		float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows));
+		if (radius > radiusMax) {
+                  resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+		  continue;
+                }
+
                 psImageInit(A, 0.0);
                 psVectorInit(B, 0.0);
@@ -275,8 +292,7 @@
 
                 if (!psMatrixGJSolve(A, B)) {
-                    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (y,x) = (%d,%d)'s residuals",
-                            oy, ox);
-                    psFree(resid); resid = NULL;
-                    break;
+		    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
+                    psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox);
+		    continue;
                 }
 
@@ -286,11 +302,11 @@
 
                 float dRo = sqrt(A->data.F32[0][0]);
-                // fprintf (stderr, "res: %2d %2d : %6.4f  %6.4f  %6.4f   %3d  %1d\n",
-                // ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]);
 
                 if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) {
                   resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1;
-                }
-                //resid->variance->data.F32[oy][ox] = XXX;
+		  resid->Ro->data.F32[oy][ox] = 0.0;
+		  resid->Rx->data.F32[oy][ox] = 0.0;
+		  resid->Ry->data.F32[oy][ox] = 0.0;
+                }
             }
         }