Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c
===================================================================
--- /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c	(revision 42509)
+++ /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c	(revision 42510)
@@ -95,5 +95,5 @@
     psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
 
-# define TESTRESID 1
+# define TESTRESID 0
 # if (TESTRESID)
     psArray *sourceRawCube = psArrayAllocEmpty (1);
@@ -101,4 +101,6 @@
     psArray *sourceVarCube = psArrayAllocEmpty (1);
     psArray *sourceMskCube = psArrayAllocEmpty (1);
+
+    psImage *refImage = NULL;
 # endif
 
@@ -121,6 +123,19 @@
 # if (TESTRESID)
 	psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32);
-	psArrayAdd (sourceRawCube, 1, tmpImage);
-	psFree (tmpImage);
+	if (refImage == NULL) {
+	  refImage = tmpImage;
+	  psArrayAdd (sourceRawCube, 1, tmpImage);
+	} else {
+	  if (tmpImage->numCols != refImage->numCols) {
+	    fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+	    continue;
+	  }
+	  if (tmpImage->numRows != refImage->numRows) {
+	    fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+	    continue;
+	  }
+	  psArrayAdd (sourceRawCube, 1, tmpImage);
+	  psFree (tmpImage);
+	}
 # endif
 
@@ -199,7 +214,4 @@
     psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK);
 
-    // statistic to use to determine baseline for clipping
-    psStats *fluxClip     = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psStats *fluxClipDef  = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     // statistic to use to determine output flux
     // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV?
@@ -207,4 +219,8 @@
     psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
 
+    // use this for the IRLS fitting below (defines niter)
+    psStats *statsIRLS = psStatsAlloc(PS_STAT_CLIPPED_MEAN); statsIRLS->clipIter = 10; // max number of iterations
+    psPolynomial2D *polyIRLS = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
+
     // Use psF64 to minimize overflow problems?
     psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
@@ -217,4 +233,11 @@
 
             int nGoodPixel = 0;              // pixel is off the image
+
+	    // skip pixels outside of the requested radius
+	    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;
+	    }
 
             // build the vector of data values for this output pixel
@@ -268,19 +291,21 @@
 
             // measure the robust median to determine a baseline reference value
-            *fluxClip = *fluxClipDef;
+	    psStats *fluxClip = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
             if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) {
 		psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
+		psFree (fluxClip);
 		return false;
 	    }
-	    if (isnan(fluxClip->robustMedian)) {
+	    if (isnan(fluxClip->sampleMedian)) {
                 resid->Ro->data.F32[oy][ox] = 0.0;
                 resid->Rx->data.F32[oy][ox] = 0.0;
                 resid->Ry->data.F32[oy][ox] = 0.0;
                 resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
+		psFree (fluxClip);
                 continue;
 	    }
 
-	    if (fabs(fluxClip->robustMedian) > 1.5) {
-	      fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->robustMedian);
+	    if (fabs(fluxClip->sampleMedian) > 1.5) {
+	      fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->sampleMedian);
 	    }
 
@@ -288,5 +313,5 @@
             int nKeep = 0;
             for (int i = 0; i < fluxes->n; i++) {
-                float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
+                float delta = fluxes->data.F32[i] - fluxClip->sampleMedian;
                 float sigma = sqrt (dfluxes->data.F32[i]);
                 float swing = fabs(delta) / sigma;
@@ -298,4 +323,5 @@
                 if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++;
             }
+	    psFree (fluxClip);
 
 	    if (nKeep < 5) {
@@ -311,18 +337,9 @@
 		}
 
-		// XXX this test should go at the top of the loop
-		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->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
 
                 resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
-                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
-
 		if (isnan(resid->Ro->data.F32[oy][ox])) {
 		    resid->Ro->data.F32[oy][ox] = 0.0;
-		    resid->Rx->data.F32[oy][ox] = 0.0;
-		    resid->Ry->data.F32[oy][ox] = 0.0;
 		    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
 		    continue;
@@ -335,45 +352,18 @@
                 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);
-                for (int i = 0; i < fluxes->n; i++) {
-                    if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
-                    B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i];
-                    B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i];
-                    B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
-
-                    A->data.F64[0][0] += 1.0/dfluxes->data.F32[i];
-                    A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i];
-                    A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i];
-
-                    A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i];
-                    A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i];
-                    A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i];
-                }
-
-                A->data.F64[0][1] = A->data.F64[1][0];
-                A->data.F64[0][2] = A->data.F64[2][0];
-                A->data.F64[2][1] = A->data.F64[1][2];
-
-                if (!psMatrixGJSolve(A, B)) {
-		    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;
-                }
-
-                resid->Ro->data.F32[oy][ox] = B->data.F64[0];
-                resid->Rx->data.F32[oy][ox] = B->data.F64[1];
-                resid->Ry->data.F32[oy][ox] = B->data.F64[2];
-
-                float dRo = sqrt(A->data.F32[0][0]);
-
-		if (fabs(B->data.F64[0]) > 1.5) {
-		  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, B->data.F64[0]);
+		// we have the following vectors to describe this pixel:
+		// fluxes, dfluxes, fmasks, xC, yC
+
+		psVectorIRLSFitPolynomial2D (polyIRLS, statsIRLS, fmasks, fmaskVal, fluxes, dfluxes, xC, yC);
+
+                resid->Ro->data.F32[oy][ox] = polyIRLS->coeff[0][0];
+                resid->Rx->data.F32[oy][ox] = polyIRLS->coeff[1][0];
+                resid->Ry->data.F32[oy][ox] = polyIRLS->coeff[0][1];
+
+		// is the error well-defined?
+                float dRo = polyIRLS->coeffErr[0][0];
+
+		if (fabs(resid->Ro->data.F32[oy][ox]) > 1.5) {
+		  fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, resid->Ro->data.F32[oy][ox]);
 		}
 
@@ -442,6 +432,4 @@
     psFree (fluxStats);
     psFree (fluxStatsDef);
-    psFree (fluxClip);
-    psFree (fluxClipDef);
 
     if (resid != NULL && psTraceGetLevel("psphot") > 5) {
@@ -453,4 +441,7 @@
     }
 
+    psFree (statsIRLS);
+    psFree (polyIRLS);
+
     psf->residuals = resid;
     return (resid != NULL) ? true : false;