Index: trunk/psphot/src/Makefile.am
===================================================================
--- trunk/psphot/src/Makefile.am	(revision 12923)
+++ trunk/psphot/src/Makefile.am	(revision 12950)
@@ -56,4 +56,6 @@
 	psphotDeblendSatstars.c	 \
 	psphotMosaicSubimage.c	 \
+	psphotMakeResiduals.c	 \
+	psphotTestSourceOutput.c \
 	psphotAddNoise.c
 
Index: trunk/psphot/src/psphot.c
===================================================================
--- trunk/psphot/src/psphot.c	(revision 12923)
+++ trunk/psphot/src/psphot.c	(revision 12950)
@@ -18,5 +18,5 @@
     // load input data (config and images (signal, noise, mask)
     if (!psphotParseCamera (config)) {
-        psErrorStackPrint(stderr, "Error setting up the camera");
+        psErrorStackPrint(stderr, "Error setting up the camera\n");
         exit (psphotGetExitStatus());
     }
@@ -24,5 +24,5 @@
     // call psphot for each readout
     if (!psphotImageLoop (config)) {
-        psErrorStackPrint(stderr, "Error in the psphot image loop");
+        psErrorStackPrint(stderr, "Error in the psphot image loop\n");
         exit (psphotGetExitStatus());
     }
Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 12923)
+++ trunk/psphot/src/psphot.h	(revision 12950)
@@ -21,8 +21,9 @@
 bool            psphotDefineFiles (pmConfig *config, pmFPAfile *input);
 
+
 // XXX test functions
 bool            psphotTestPSF (pmReadout *readout, psArray *sources, psMetadata *recipe);
 bool            pmPSFtestModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask);
-
+bool psphotTestSourceOutput (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf);
 psArray        *psphotFakeSources ();
 
@@ -104,3 +105,5 @@
 bool 		psphotDeblendSatstars (psArray *sources, psMetadata *recipe);
 
+bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf);
+
 #endif
Index: trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- trunk/psphot/src/psphotChoosePSF.c	(revision 12923)
+++ trunk/psphot/src/psphotChoosePSF.c	(revision 12950)
@@ -221,4 +221,10 @@
     }
 
+    // build a PSF residual image
+    if (!psphotMakeResiduals (try->sources, recipe, try->psf)) {
+	psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF");
+	return NULL;
+    }
+
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) { 
Index: trunk/psphot/src/psphotGuessModels.c
===================================================================
--- trunk/psphot/src/psphotGuessModels.c	(revision 12923)
+++ trunk/psphot/src/psphotGuessModels.c	(revision 12950)
@@ -55,4 +55,5 @@
     // set the source PSF model
     source->modelPSF = modelPSF;
+    source->modelPSF->residuals = psf->residuals;
   }
   psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot"));
Index: trunk/psphot/src/psphotMakeResiduals.c
===================================================================
--- trunk/psphot/src/psphotMakeResiduals.c	(revision 12950)
+++ trunk/psphot/src/psphotMakeResiduals.c	(revision 12950)
@@ -0,0 +1,278 @@
+# include "psphotInternal.h"
+# define ZERO_ORDER 0
+
+bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf) {
+
+    bool status, isPSF;
+    double flux, dflux;
+    psU8 mflux;
+
+    psTimerStart ("residuals");
+
+    if (!psMetadataLookupBool(&status, recipe, "PSF.RESIDUALS")) return true;
+
+    int xBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.XBIN");
+    PS_ASSERT (status, false);
+
+    int yBin = psMetadataLookupS32(&status, recipe, "PSF.RESIDUALS.YBIN");
+    PS_ASSERT (status, false);
+
+    float nSigma = psMetadataLookupF32(&status, recipe, "PSF.RESIDUALS.NSIGMA");
+    PS_ASSERT (status, false);
+
+    char *modeString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.INTERPOLATION");
+    PS_ASSERT (status, false);
+
+    psImageInterpolateMode mode = psImageInterpolateModeFromString (modeString);
+    if (mode == PS_INTERPOLATE_NONE) {
+	psError(PSPHOT_ERR_CONFIG, false, "invalid interpolation in psphot.config");
+	return false;
+    }
+
+    char *statString = psMetadataLookupStr(&status, recipe, "PSF.RESIDUALS.STATISTIC");
+    PS_ASSERT (status, false);
+
+    psStatsOptions statOption = psStatsOptionFromString (statString);
+    if (!statOption) {
+	psError(PSPHOT_ERR_CONFIG, false, "invalid residual statistic in psphot.config");
+	return false;
+    }
+
+    // user parameters:
+    // size of aperture (determine from source images?)
+    // binning factor
+
+    // select the subset of sources which are the PSFSTARs
+    // for each input source:
+    // - construct a residual image, renormalized
+    // - construct a renormalized weight image
+    // - construct a new mask image 
+
+    // construct the output residual table (Nx*DX,Ny*DY)
+    // for each output pixel:
+    // - construct a histogram of the values & weights (interpolate to the common pixel coordinate)
+    // - measure the robust median & sigma
+    // - reject (mask) input pixels which are outliers
+    // - re-measure the robust median & sigma
+    // - set output pixel, weight, and mask
+
+    // determine the maximum image size from the input sources
+    int xSize = 0;
+    int ySize = 0;
+
+    psVector *xC = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    // build (DATA - MODEL) [an image] for each psf star
+    psArray *input = psArrayAllocEmpty (100);
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        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
+
+        psImage *image  = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+        psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32);
+        psImage *mask   = psImageCopy (NULL, source->mask,   PS_TYPE_U8);
+        pmModelSub (image, mask, model, false, false);
+	
+	// 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, 1, 2, 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);
+    psVector *fmasks  = psVectorAlloc (input->n, PS_TYPE_U8);
+
+    // 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?
+    psStats *fluxStats    = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV);
+    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++) {
+	fprintf (stderr, ".");
+	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;
+		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]);
+	    }
+
+	    // measure the robust median to determine a baseline reference value
+	    *fluxClip = *fluxClipDef;
+	    psVectorStats (fluxClip, fluxes, NULL, fmasks, 0xff);
+
+	    // 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;
+		}
+	    }		    
+
+	    // 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++) {
+	fprintf (stderr, ".");
+	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];
+		
+		// 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]);
+	    }
+
+	    // 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];
+	}
+    }
+
+    psFree (A);
+    psFree (B);
+
+# endif
+
+    psFree (xC);
+    psFree (yC);
+    psFree (input);
+
+    psFree (fluxes);
+    psFree (dfluxes);
+    psFree (fmasks);
+
+    psFree (fluxStats);
+    psFree (fluxStatsDef);
+    psFree (fluxClip);
+    psFree (fluxClipDef);
+
+    psphotSaveImage (NULL, resid->Ro,     "resid.ro.fits");
+    psphotSaveImage (NULL, resid->Rx,     "resid.rx.fits");
+    psphotSaveImage (NULL, resid->Ry,     "resid.ry.fits");
+    psphotSaveImage (NULL, resid->weight, "resid.wt.fits");
+    psphotSaveImage (NULL, resid->mask,   "resid.mk.fits");
+
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals"));
+
+    psf->residuals = resid;
+    return true;
+}
Index: trunk/psphot/src/psphotReadout.c
===================================================================
--- trunk/psphot/src/psphotReadout.c	(revision 12923)
+++ trunk/psphot/src/psphotReadout.c	(revision 12950)
@@ -64,4 +64,7 @@
     }
 
+    // psArray *bounds = psphotGetBounds ();
+    // psArray *sources = psphotBoundsToSources ();
+
     // construct sources and measure basic stats
     psArray *sources = psphotSourceStats (readout, recipe, peaks);
@@ -106,4 +109,7 @@
     psphotGuessModels (readout, sources, recipe, psf);
 
+    // XXX test output of models 
+    // psphotTestSourceOutput (readout, sources, recipe, psf);
+
     if (dump) psphotSaveImage (NULL, readout->image,  "image.v0.fits");
 
Index: trunk/psphot/src/psphotTestSourceOutput.c
===================================================================
--- trunk/psphot/src/psphotTestSourceOutput.c	(revision 12950)
+++ trunk/psphot/src/psphotTestSourceOutput.c	(revision 12950)
@@ -0,0 +1,158 @@
+# include "psphotInternal.h"
+
+enum {
+    PSPHOT_ADD_NONE = 0,
+    PSPHOT_ADD_MODEL = 1,
+    PSPHOT_ADD_R0 = 2,
+    PSPHOT_ADD_R1 = 4,
+};
+
+bool psphotAddModel(psImage *image,
+		    pmModel *model,
+		    int mode
+    )
+{
+    psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+
+    PS_ASSERT_PTR_NON_NULL(model, false);
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
+
+    psVector *x = psVectorAlloc(2, PS_TYPE_F32);
+    psVector *params = model->params;
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
+    psS32 imageCol;
+    psS32 imageRow;
+    psF32 skyValue = params->data.F32[0];
+    psF32 pixelValue;
+    
+    float xCenter = model->params->data.F32[PM_PAR_XPOS];
+    float yCenter = model->params->data.F32[PM_PAR_YPOS];
+    float Io = model->params->data.F32[PM_PAR_I0];
+
+    int xBin = 1;
+    int yBin = 1;
+    float xResidCenter = 0.0;
+    float yResidCenter = 0.0;
+
+    psImageInterpolateOptions *Ro = NULL;
+    psImageInterpolateOptions *Rx = NULL;
+    psImageInterpolateOptions *Ry = NULL;
+    if (model->residuals && (mode & (PSPHOT_ADD_R0 | PSPHOT_ADD_R1))) {
+	Ro = psImageInterpolateOptionsAlloc(
+	    PS_INTERPOLATE_BILINEAR,
+	    model->residuals->Ro, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+	Rx = psImageInterpolateOptionsAlloc(
+	    PS_INTERPOLATE_BILINEAR,
+	    model->residuals->Rx, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+	Ry = psImageInterpolateOptionsAlloc(
+	    PS_INTERPOLATE_BILINEAR,
+	    model->residuals->Ry, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+
+	xBin = model->residuals->xBin;
+	yBin = model->residuals->yBin;
+	xResidCenter = model->residuals->xCenter;
+	yResidCenter = model->residuals->yCenter;
+    }
+
+    for (psS32 iy = 0; iy < image->numRows; iy++) {
+        for (psS32 ix = 0; ix < image->numCols; ix++) {
+
+            // Convert i/j to image coord space:
+	    imageCol = ix + image->col0;
+	    imageRow = iy + image->row0;
+
+            x->data.F32[0] = (float) imageCol;
+            x->data.F32[1] = (float) imageRow;
+
+            // set the appropriate pixel value for this coordinate
+	    if (mode & PSPHOT_ADD_MODEL) {
+		pixelValue = modelFunc (NULL, params, x) - skyValue;
+	    } else {
+		pixelValue = 0.0;
+	    }
+
+	    // get the contribution from the residual model
+	    // XXX for a test, do this for all sources and all pixels
+	    if (Ro) {
+		// fractional image position
+		// this is wrong for the 'center' case
+		float ox = xBin*(ix + 0.5 + image->col0 - xCenter) + xResidCenter;
+		float oy = yBin*(iy + 0.5 + image->row0 - yCenter) + yResidCenter;
+
+		psU8 mflux = 0;
+		double Fo = 0.0;
+		double Fx = 0.0;
+		double Fy = 0.0;
+		psImageInterpolate (&Fo, NULL, &mflux, ox, oy, Ro);
+		psImageInterpolate (&Fx, NULL, &mflux, ox, oy, Rx);
+		psImageInterpolate (&Fy, NULL, &mflux, ox, oy, Ry);
+
+		if (!mflux && isfinite(Fo) && isfinite(Fx) && isfinite(Fy)) {
+		    if (mode & PSPHOT_ADD_R0) {
+			pixelValue += Io*Fo;
+		    }
+		    if (mode & PSPHOT_ADD_R1) {
+			pixelValue += Io*(xCenter*Fx + yCenter*Fy);
+		    }
+		}
+	    }
+	    image->data.F32[iy][ix] += pixelValue;
+        }
+    }
+    psFree(x);
+    psFree(Ro);
+    psFree(Rx);
+    psFree(Ry);
+    psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
+    return(true);
+}
+
+// construct an initial PSF model for each object 
+bool psphotTestSourceOutput (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+
+    psImage *imMo = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+    psImage *imR0 = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+    psImage *imR1 = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
+    
+    // create template model
+    pmModel *modelRef = pmModelAlloc(psf->type);
+    modelRef->params->data.F32[PM_PAR_SKY] = 0;
+    modelRef->params->data.F32[PM_PAR_I0] = 1000;
+
+    int dx = 25;
+    int dy = 25;
+
+    // generate a grid of fake sources with amplitude 1000
+    for (int iy = 50; iy < imMo->numRows; iy += 100) {
+	for (int ix = 50; ix < imMo->numCols; ix += 100) {
+	    
+	    // assign the x and y coords to the image center
+	    modelRef->params->data.F32[PM_PAR_XPOS] = ix;
+	    modelRef->params->data.F32[PM_PAR_YPOS] = iy;
+	    
+	    // create modelPSF from this model
+	    pmModel *model = pmModelFromPSF (modelRef, psf);
+	    model->residuals = psf->residuals;
+
+	    // generate working image for this source
+	    psRegion region = {ix - dx, ix + dx, iy - dy, iy + dy};
+
+	    psImage *vM = psImageSubset (imMo, region);
+	    psImage *v0 = psImageSubset (imR0, region);
+	    psImage *v1 = psImageSubset (imR1, region);
+
+	    // we want to make one image o
+	    psphotAddModel (vM, model, PSPHOT_ADD_MODEL);
+	    psphotAddModel (v0, model, PSPHOT_ADD_R0);
+	    psphotAddModel (v1, model, PSPHOT_ADD_R1);
+	}
+    }
+
+    psphotSaveImage (NULL, imMo, "grid.Mo.fits");
+    psphotSaveImage (NULL, imR0, "grid.R0.fits");
+    psphotSaveImage (NULL, imR1, "grid.R1.fits");
+
+    exit (0);
+}
+