Index: /trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- /trunk/psphot/src/psphotMakeResiduals.c	(revision 13406)
+++ /trunk/psphot/src/psphotMakeResiduals.c	(revision 13407)
@@ -1,4 +1,3 @@
 # include "psphotInternal.h"
-# define ZERO_ORDER 0
 
 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
@@ -12,9 +11,25 @@
     if (!psMetadataLookupBool(&status, recipe, "PSF.RESIDUALS")) return true;
 
+    int SPATIAL_ORDER = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.SPATIAL_ORDER");
+#if 0					// True when PSF.RESIDUALS.SPATIAL_ORDER's in default recipe
+    PS_ASSERT (status, false);
+#else
+    if (!status) {
+	SPATIAL_ORDER = 1;		// i.e. linear
+    }
+#endif
+    if (SPATIAL_ORDER != 0 && SPATIAL_ORDER != 1) {
+	psError(PSPHOT_ERR_CONFIG, true, "PSF.RESIDUALS.SPATIAL_ORDER must be 0 or 1 (not %d)",
+		SPATIAL_ORDER);
+	return false;
+    }
+
     int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN");
     PS_ASSERT (status, false);
+    psErrorClear(); PS_ASSERT (xBin != 0, false);
 
     int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN");
     PS_ASSERT (status, false);
+    psErrorClear(); PS_ASSERT (yBin != 0, false);
 
     float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
@@ -122,8 +137,10 @@
     psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
 
-// this section builds just the 0th order term Ro
-# if (ZERO_ORDER)
-    // build Ro = DATA - MODEL (rebinned image) pixel-by-pixel
-    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
+    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
+    psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
+
+    // Solve MODEL = R + x R_x + y R_y in pixel-by-pixel a least-squares sense
+    // (If SPATIAL_ORDER == 0, just solve MODEL = R)
+    for (int oy = 0; resid != NULL && oy < resid->Ro->numRows; oy++) {
 	for (int ox = 0; ox < resid->Ro->numCols; ox++) {
 	    
@@ -148,4 +165,5 @@
 	    *fluxClip = *fluxClipDef;
 	    psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
+	    psErrorClear();		// clear (ignore) any outstanding errors 
 
 	    // mark input pixels which are more than N sigma from the median
@@ -161,86 +179,45 @@
 	    }		    
 
-	    // measure the desired statistic on the unclipped pixels
-	    *fluxStats = *fluxStatsDef;
-	    psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
-
-	    resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
-	    resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
-
-	    // clear (ignore) any outstanding errors 
-	    psErrorClear(); 
-	}
-    }
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generated 0-th order residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
-
-# else
-
-    psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix
-    psVector *B = psVectorAlloc(3, PS_TYPE_F64); // Least-squares vector
-
-    // build (x,y)*(DATA - MODEL - Ro) pixel-by-pixel
-    for (int oy = 0; oy < resid->Ro->numRows; oy++) {
-	for (int ox = 0; ox < resid->Ro->numCols; ox++) {
-	    
-	    // build the vector of data values for this output pixel
-	    // XXX this is identical to the pass above: we could cache the results for speed
-	    for (int i = 0; i < input->n; i++) {
-
-		psImageInterpolateOptions *interp = input->data[i];
+	    if (SPATIAL_ORDER == 0) {
+		// measure the desired statistic on the unclipped pixels
+		*fluxStats = *fluxStatsDef;
+		psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
 		
-		// fractional image position
-		float ix = (ox + 0.5 - resid->xCenter) / (float) xBin + xC->data.F32[i] - interp->image->col0;
-		float iy = (oy + 0.5 - resid->yCenter) / (float) yBin + yC->data.F32[i] - interp->image->row0;
-
-		mflux = 0;
-		psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp);
-		fluxes->data.F32[i] = flux;
-		dfluxes->data.F32[i] = dflux;
-		fmasks->data.U8[i] = mflux;
-		// fprintf (stderr, "%f %f : %f %f (%d)\n", ix, iy, flux, dflux, fmasks->data.U8[i]);
+		resid->Ro->data.F32[oy][ox] = psStatsGetValue (fluxStats, statOption);
+		//resid->weight->data.F32[oy][ox] = fluxStats->sampleStdev;
+	    } else {
+		assert (SPATIAL_ORDER == 1);
+		psImageInit(A, 0.0);
+		psVectorInit(B, 0.0);
+		for (int i = 0; i < fluxes->n; i++) {
+		    if (fmasks->data.U8[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)) {
+		    psError(PSPHOT_ERR_PSF, false, "Singular matrix solving for (r,x) = (%d,%d)'s residuals",
+			    oy, ox);
+		    psFree(resid); resid = NULL;
+		    break;		    
+		}
+
+		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];
 	    }
-
-	    // measure the robust median to determine a baseline reference value
-	    *fluxClip = *fluxClipDef;
-	    psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
-	    psErrorClear();		// clear (ignore) any outstanding errors 
-
-	    // mark input pixels which are more than N sigma from the median
-	    for (int i = 0; i < fluxes->n; i++) {
-		float delta = fluxes->data.F32[i] - fluxClip->robustMedian;
-		float sigma = sqrt (dfluxes->data.F32[i]);
-		float swing = fabs(delta) / sigma;
-
-		// make this a user option
-		if (swing > nSigma) {
-		    fmasks->data.U8[i] = 1;
-		}
-	    }		    
-
-	    psImageInit(A, 0.0);
-	    psVectorInit(B, 0.0);
-	    for (int i = 0; i < fluxes->n; i++) {
-		if (fmasks->data.U8[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];
-	    psMatrixGJSolve(A, B);
-
-	    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];
 	}
     }
@@ -248,6 +225,4 @@
     psFree (A);
     psFree (B);
-
-# endif
 
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
@@ -266,5 +241,5 @@
     psFree (fluxClipDef);
 
-    if (psTraceGetLevel("psphot") > 5) {
+    if (resid != NULL && psTraceGetLevel("psphot") > 5) {
       psphotSaveImage (NULL, resid->Ro,     "resid.ro.fits");
       psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
@@ -275,4 +250,4 @@
 
     psf->residuals = resid;
-    return true;
+    return (resid != NULL) ? true : false;
 }