Index: /trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- /trunk/psphot/src/psphotMakeResiduals.c	(revision 13805)
+++ /trunk/psphot/src/psphotMakeResiduals.c	(revision 13806)
@@ -14,7 +14,7 @@
     PS_ASSERT (status, false);
     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;
+        psError(PSPHOT_ERR_CONFIG, true, "PSF.RESIDUALS.SPATIAL_ORDER must be 0 or 1 (not %d)",
+                SPATIAL_ORDER);
+        return false;
     }
 
@@ -35,6 +35,6 @@
     psImageInterpolateMode mode = psImageInterpolateModeFromString (modeString);
     if (mode == PS_INTERPOLATE_NONE) {
-	psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
-	return false;
+        psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
+        return false;
     }
 
@@ -44,6 +44,6 @@
     psStatsOptions statOption = psStatsOptionFromString (statString);
     if (!statOption) {
-	psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
-	return false;
+        psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
+        return false;
     }
 
@@ -56,5 +56,5 @@
     // - construct a residual image, renormalized
     // - construct a renormalized weight image
-    // - construct a new mask image 
+    // - construct a new mask image
 
     // construct the output residual table (Nx*DX,Ny*DY)
@@ -66,9 +66,9 @@
     // - set output pixel, weight, and mask
 
-    const int badMask = 1;		// mask bits
-    const int poorMask = 2;		//       from psImageInterpolate
-    const int clippedMask = 4;		// mask bit set for clipped values
-
-    bool offImage = false;		// pixel is off the image
+    const int badMask = 1;              // mask bits
+    const int poorMask = 2;             //       from psImageInterpolate
+    const int clippedMask = 4;          // mask bit set for clipped values
+
+    bool offImage = false;              // pixel is off the image
 
     // determine the maximum image size from the input sources
@@ -87,7 +87,7 @@
         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
 
-	// which model to use?
-	pmModel *model = pmSourceGetModel (&isPSF, source);
-	if (model == NULL) continue;  // model must be defined
+        // which model to use?
+        pmModel *model = pmSourceGetModel (&isPSF, source);
+        if (model == NULL) continue;  // model must be defined
 
         psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
@@ -95,35 +95,35 @@
         psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_U8);
         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC);
-	
-	// re-normalize image and weight
-	float Io = model->params->data.F32[PM_PAR_I0];
-	psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
-	psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
-
-	// we will interpolate the image and weight - include the mask or not?
-	// XXX consider better values for the mask bits
-	psImageInterpolateOptions *interp =
-	    psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0);
-	psArrayAdd (input,  100, interp);
-
-	// save the X,Y position for future reference 
-	xC->data.F32[xC->n] = model->params->data.F32[PM_PAR_XPOS];
-	yC->data.F32[yC->n] = model->params->data.F32[PM_PAR_YPOS];
-	psVectorExtend (xC, 100, 1);
-	psVectorExtend (yC, 100, 1);
-
-	xSize = PS_MAX (xSize, image->numCols);
-	ySize = PS_MAX (ySize, image->numRows);
-
-	// free up the excess references 
-	psFree (mask);
-	psFree (image);
-	psFree (weight);
-	psFree (interp);
+
+        // re-normalize image and weight
+        float Io = model->params->data.F32[PM_PAR_I0];
+        psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32));
+        psBinaryOp (weight, weight, "/", psScalarAlloc(Io*Io, PS_TYPE_F32));
+
+        // we will interpolate the image and weight - include the mask or not?
+        // XXX consider better values for the mask bits
+        psImageInterpolateOptions *interp =
+            psImageInterpolateOptionsAlloc(mode, image, weight, NULL, 0xff, 0.0, 0.0, badMask, poorMask, 0.0);
+        psArrayAdd (input,  100, interp);
+
+        // save the X,Y position for future reference
+        xC->data.F32[xC->n] = model->params->data.F32[PM_PAR_XPOS];
+        yC->data.F32[yC->n] = model->params->data.F32[PM_PAR_YPOS];
+        psVectorExtend (xC, 100, 1);
+        psVectorExtend (yC, 100, 1);
+
+        xSize = PS_MAX (xSize, image->numCols);
+        ySize = PS_MAX (ySize, image->numRows);
+
+        // free up the excess references
+        psFree (mask);
+        psFree (image);
+        psFree (weight);
+        psFree (interp);
     }
     pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin);
 
     // x(resid) = (x(image) - Xo)*xBin + xCenter
-    
+
     psVector *fluxes  = psVectorAlloc (input->n, PS_TYPE_F32);
     psVector *dfluxes = psVectorAlloc (input->n, PS_TYPE_F32);
@@ -144,88 +144,89 @@
     // (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++) {
-	    
-	    // build the vector of data values for this output pixel
-	    for (int i = 0; i < input->n; i++) {
-
-		psImageInterpolateOptions *interp = input->data[i];
-		
-		// 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;
-		offImage = false;
-		if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
-		    // This pixel is off the image
-		    offImage = true;
-		}
-		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]);
-	    }
-
-	    // 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] = clippedMask;
-		}
-	    }		    
-
-	    if (offImage || SPATIAL_ORDER == 0) {
-		// 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->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
-		//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 (y,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];
-		//resid->weight->data.F32[oy][ox] = XXX;
-	    }
-	}
+        for (int ox = 0; ox < resid->Ro->numCols; ox++) {
+
+            // build the vector of data values for this output pixel
+            for (int i = 0; i < input->n; i++) {
+
+                psImageInterpolateOptions *interp = input->data[i];
+
+                // 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;
+                offImage = false;
+                if (psImageInterpolate (&flux, &dflux, &mflux, ix, iy, interp) == PS_INTERPOLATE_STATUS_OFF) {
+                    // This pixel is off the image
+                    offImage = true;
+                }
+                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]);
+            }
+
+            // 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] = clippedMask;
+                }
+            }
+
+            if (offImage || SPATIAL_ORDER == 0) {
+                // measure the desired statistic on the unclipped pixels
+                *fluxStats = *fluxStatsDef;
+                psVectorStats (fluxStats, fluxes, NULL, fmasks, 0xff);
+                psErrorClear();         // clear (ignore) any outstanding errors
+
+                resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
+                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
+                //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 (y,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];
+                //resid->weight->data.F32[oy][ox] = XXX;
+            }
+        }
     }