Index: /trunk/psModules/src/camera/pmFPACopy.c
===================================================================
--- /trunk/psModules/src/camera/pmFPACopy.c	(revision 13033)
+++ /trunk/psModules/src/camera/pmFPACopy.c	(revision 13034)
@@ -211,4 +211,17 @@
         if (mdok && trimsec && !psRegionIsNaN(*trimsec)) {
             *trimsec = psImageBinningSetRuffRegion(binning, *trimsec);
+	    // force integer pixels : truncate x0, roundup x1:
+	    trimsec->x0 = (int)trimsec->x0;
+	    if (trimsec->x1 > (int)trimsec->x1) {
+		trimsec->x1 = (int)trimsec->x1 + 1;
+	    } else {
+		trimsec->x1 = (int)trimsec->x1;
+	    }		
+	    trimsec->y0 = (int)trimsec->y0;
+	    if (trimsec->y1 > (int)trimsec->y1) {
+		trimsec->y1 = (int)trimsec->y1 + 1;
+	    } else {
+		trimsec->y1 = (int)trimsec->y1;
+	    }		
         }
         psList *biassecs = psMetadataLookupPtr(&mdok, target->concepts, "CELL.BIASSEC"); // The bias sections
@@ -219,4 +232,17 @@
                 if (!psRegionIsNaN(*biassec)) {
                     *biassec = psImageBinningSetRuffRegion(binning, *biassec);
+		    // force integer pixels : truncate x0, roundup x1:
+		    biassec->x0 = (int)biassec->x0;
+		    if (biassec->x1 > (int)biassec->x1) {
+			biassec->x1 = (int)biassec->x1 + 1;
+		    } else {
+			biassec->x1 = (int)biassec->x1;
+		    }		
+		    biassec->y0 = (int)biassec->y0;
+		    if (biassec->y1 > (int)biassec->y1) {
+			biassec->y1 = (int)biassec->y1 + 1;
+		    } else {
+			biassec->y1 = (int)biassec->y1;
+		    }		
                 }
             }
Index: /trunk/psModules/src/camera/pmFPAfile.c
===================================================================
--- /trunk/psModules/src/camera/pmFPAfile.c	(revision 13033)
+++ /trunk/psModules/src/camera/pmFPAfile.c	(revision 13034)
@@ -17,4 +17,5 @@
 #include "pmFPAfile.h"
 #include "pmFPACopy.h"
+#include "pmConcepts.h"
 
 static void pmFPAfileFree(pmFPAfile *file)
@@ -326,6 +327,10 @@
     PS_ASSERT_PTR_NON_NULL(view, false);
 
+    // XXX this should be smarter (ie, only copy concepts from the current chips)
+    // but such a call is needed, so re-copy stuff rather than no copy
+    pmFPACopyConcepts (out, in);
+
     // pmFPAWrite takes care of all PHUs as needed
-    if ((view->chip == -1) || in->hdu) {
+    if (view->chip == -1) {
         status = pmFPACopyStructure (out, in, xBin, yBin);
         return status;
@@ -338,8 +343,8 @@
     pmChip *outChip = out->chips->data[view->chip];
 
-    if ((view->cell == -1) || inChip->hdu) {
+    if (view->cell == -1) {
         status = pmChipCopyStructure (outChip, inChip, xBin, yBin);
         return status;
-    }
+    } 
     if (view->cell >= inChip->cells->n) {
         psError(PS_ERR_IO, true, "Requested cell == %d>= inChip->cells->n == %ld",
@@ -350,10 +355,6 @@
     pmCell *outCell = outChip->cells->data[view->cell];
 
-    if ((view->readout == -1) || inCell->hdu) {
-        status = pmCellCopyStructure (outCell, inCell, xBin, yBin);
-        return status;
-    }
-    psError(PS_ERR_UNKNOWN, true, "Returning false");
-    return false;
+    status = pmCellCopyStructure (outCell, inCell, xBin, yBin);
+    return status;
 }
 
Index: /trunk/psModules/src/camera/pmFPAfileIO.c
===================================================================
--- /trunk/psModules/src/camera/pmFPAfileIO.c	(revision 13033)
+++ /trunk/psModules/src/camera/pmFPAfileIO.c	(revision 13034)
@@ -434,4 +434,5 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "closing %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         status = psFitsClose (file->fits);
@@ -447,5 +448,4 @@
       case PM_FPA_FILE_OBJ:
       case PM_FPA_FILE_CMP:
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -500,7 +500,7 @@
       case PM_FPA_FILE_CMP:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
         return true;
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -638,4 +638,5 @@
       case PM_FPA_FILE_FRINGE:
       case PM_FPA_FILE_CMF:
+      case PM_FPA_FILE_PSF:
         psTrace ("psModules.camera", 5, "opening %s (%s) (%d:%d:%d)\n", file->filename, file->name, view->chip, view->cell, view->readout);
         file->fits = psFitsOpen (file->filename, mode);
@@ -672,5 +673,4 @@
       case PM_FPA_FILE_CMP:
       case PM_FPA_FILE_RAW:
-      case PM_FPA_FILE_PSF:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -751,4 +751,5 @@
 
 // XXX this function is only called from pmFPAfileWrite
+// XXX for each data type, there should be a function which writes the PHU, if needed
 bool pmFPAfileWritePHU(pmFPAfile *file, const pmFPAview *view, const pmConfig *config)
 {
Index: /trunk/psModules/src/concepts/pmConcepts.c
===================================================================
--- /trunk/psModules/src/concepts/pmConcepts.c	(revision 13033)
+++ /trunk/psModules/src/concepts/pmConcepts.c	(revision 13034)
@@ -888,4 +888,8 @@
 
 
+// XXX EAM : Paul, please handle the biassec independently so the 
+// psList is copied without a verbose warning.  when this is fixed, you can
+// convert the trace back to a log (psMetadata.c:429)
+
 bool pmFPACopyConcepts(pmFPA *target, const pmFPA *source)
 {
Index: /trunk/psModules/src/concepts/pmConceptsUpdate.c
===================================================================
--- /trunk/psModules/src/concepts/pmConceptsUpdate.c	(revision 13033)
+++ /trunk/psModules/src/concepts/pmConceptsUpdate.c	(revision 13034)
@@ -44,4 +44,17 @@
 	    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
 	    *trimsec = psImageBinningSetRuffRegion (binning, *trimsec);
+	    // force integer pixels : truncate x0, roundup x1:
+	    trimsec->x0 = (int)trimsec->x0;
+	    if (trimsec->x1 > (int)trimsec->x1) {
+		trimsec->x1 = (int)trimsec->x1 + 1;
+	    } else {
+		trimsec->x1 = (int)trimsec->x1;
+	    }		
+	    trimsec->y0 = (int)trimsec->y0;
+	    if (trimsec->y1 > (int)trimsec->y1) {
+		trimsec->y1 = (int)trimsec->y1 + 1;
+	    } else {
+		trimsec->y1 = (int)trimsec->y1;
+	    }		
 	    psMetadataRemoveKey(cell->concepts, "CELL.TRIMSEC.UPDATE");
         }
@@ -54,4 +67,17 @@
 	    while ((biassec = psListGetAndIncrement(biassecsIter))) {
 		*biassec = psImageBinningSetRuffRegion (binning, *biassec);
+		// force integer pixels : truncate x0, roundup x1:
+		biassec->x0 = (int)biassec->x0;
+		if (biassec->x1 > (int)biassec->x1) {
+		    biassec->x1 = (int)biassec->x1 + 1;
+		} else {
+		    biassec->x1 = (int)biassec->x1;
+		}		
+		biassec->y0 = (int)biassec->y0;
+		if (biassec->y1 > (int)biassec->y1) {
+		    biassec->y1 = (int)biassec->y1 + 1;
+		} else {
+		    biassec->y1 = (int)biassec->y1;
+		}		
 	    }
 	    psFree(biassecsIter);
Index: /trunk/psModules/src/objects/pmModel.h
===================================================================
--- /trunk/psModules/src/objects/pmModel.h	(revision 13033)
+++ /trunk/psModules/src/objects/pmModel.h	(revision 13034)
@@ -5,6 +5,7 @@
  * @author EAM, IfA
  *
- * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-04-21 19:47:14 $
+ * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
+ *
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -26,4 +27,15 @@
     PM_MODEL_BADARGS   ///< model fit called with invalid args
 } pmModelStatus;
+
+typedef enum {
+    PM_MODEL_OP_NONE 	= 0x00,
+    PM_MODEL_OP_FUNC 	= 0x01,
+    PM_MODEL_OP_RES0 	= 0x02,
+    PM_MODEL_OP_RES1 	= 0x04,
+    PM_MODEL_OP_FULL   = 0x07,
+    PM_MODEL_OP_SKY  	= 0x08,
+    PM_MODEL_OP_CENTER = 0x10,
+    PM_MODEL_OP_NORM  	= 0x20,
+} pmModelOpMode;
 
 /** pmModel data structure
@@ -91,11 +103,9 @@
  */
 bool pmModelAdd(
-    psImage *image,   ///< The output image (float)
-    psImage *mask,   ///< The image pixel mask (valid == 0)
-    pmModel *model,   ///< The input pmModel
-    bool center,   ///< A boolean flag that determines whether pixels are centered
-    bool sky    ///< A boolean flag that determines if the sky is subtracted
+    psImage *image,			///< The output image (float)
+    psImage *mask,			///< The image pixel mask (valid == 0)
+    pmModel *model,			///< The input pmModel
+    pmModelOpMode mode		     ///< mode to control how the model is added into the image
 );
-
 
 /** pmModelSub()
@@ -110,9 +120,8 @@
  */
 bool pmModelSub(
-    psImage *image,   ///< The output image (float)
-    psImage *mask,   ///< The image pixel mask (valid == 0)
-    pmModel *model,   ///< The input pmModel
-    bool center,   ///< A boolean flag that determines whether pixels are centered
-    bool sky    ///< A boolean flag that determines if the sky is subtracted
+    psImage *image,			///< The output image (float)
+    psImage *mask,			///< The image pixel mask (valid == 0)
+    pmModel *model,			///< The input pmModel
+    pmModelOpMode mode		     ///< mode to control how the model is added into the image
 );
 
Index: /trunk/psModules/src/objects/pmPSF.c
===================================================================
--- /trunk/psModules/src/objects/pmPSF.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmPSF.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -424,5 +424,6 @@
     // place the reference object in the image center
     // no need to mask the source here
-    pmModelAdd (image, NULL, model, false, false);
+    // XXX should we measure this for the analytical model only or the full model?
+    pmModelAdd (image, NULL, model, PM_MODEL_OP_FULL);
 
     // loop over a range of source fluxes
@@ -436,4 +437,6 @@
         psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK);
         pmSourcePhotometryAper (&apMag, model, image, mask);
+
+	// XXX since we re-mask on each pass, this could be dropped.
         psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK));
 
Index: /trunk/psModules/src/objects/pmPSF_IO.c
===================================================================
--- /trunk/psModules/src/objects/pmPSF_IO.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmPSF_IO.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44,4 +44,488 @@
 #include "pmSourceIO.h"
 
+bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    pmFPA *fpa = file->fpa;
+
+    if (view->chip == -1) {
+        if (!pmFPAWritePSFmodel (fpa, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for fpa");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        if (!pmCellWritePSFmodel (cell, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for cell");
+            return false;
+        }
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
+	psError(PS_ERR_IO, false, "Failed to write PSF for readout");
+	return false;
+    }
+    return true;
+}
+
+// read in all chip-level PSFmodel files for this FPA
+bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+
+        pmChip *chip = fpa->chips->data[i];
+        if (!pmChipWritePSFmodel (chip, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// read in all cell-level PSFmodel files for this chip
+bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < chip->cells->n; i++) {
+
+        pmCell *cell = chip->cells->data[i];
+        if (!pmCellWritePSFmodel (cell, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth cell", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// read in all readout-level PSFmodel files for this cell
+bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < cell->readouts->n; i++) {
+
+        pmReadout *readout = cell->readouts->data[i];
+        if (!pmReadoutWritePSFmodel (readout, view, file, config)) {
+            psError(PS_ERR_IO, false, "Failed to write PSF for %dth readout", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
+// and header + image for the PSF residual images
+bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    bool status;
+
+    // select the psf of interest
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    if (!psf) {
+	psError(PS_ERR_UNKNOWN, true, "missing PSF for this readout");
+	return false;
+    }
+
+    // lookup the EXTNAME values used for table data and residual image segments
+    char *rule = NULL;
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+	psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+	return false;
+    }
+    // EXTNAME for table data
+    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
+    if (!rule) {
+	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+	return false;
+    }
+    char *tableName = pmFPAfileNameFromRule (rule, file, view);
+
+    // write the PSF model parameters in a FITS table
+
+    // we need to write a header for the table, 
+    psMetadata *header = psMetadataAlloc();
+
+    char *modelName = pmModelGetType (psf->type);
+    psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
+
+    psMetadataAddBool (header, PS_LIST_TAIL, "POISSON",  0, "Use Poisson errors in fits?", psf->poissonErrors);
+
+    int nPar = pmModelParameterCount (psf->type)    ;
+    psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+    // save the dimensions of each parameter
+    for (int i = 0; i < nPar; i++) {
+	char name[9];
+        psPolynomial2D *poly = psf->params_NEW->data[i];
+        if (poly == NULL) continue;
+	snprintf (name, 9, "PAR%02d_NX", i);
+        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);
+	snprintf (name, 9, "PAR%02d_NY", i);
+        psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY);
+    }
+
+    // other required information describing the PSF
+    psMetadataAddF32  (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
+    psMetadataAddF32  (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
+    psMetadataAddF32  (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
+    psMetadataAddS32  (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
+
+    // XXX can we drop this now?
+    psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
+
+    // build a FITS table of the PSF parameters
+    psArray *psfTable = psArrayAllocEmpty (100);
+    for (int i = 0; i < nPar; i++) {
+        psPolynomial2D *poly = psf->params_NEW->data[i];
+        if (poly == NULL) continue; // skip unset parameters (eg, XPOS)
+	for (int ix = 0; ix < poly->nX; ix++) {
+	    for (int iy = 0; iy < poly->nY; iy++) {
+
+		psMetadata *row = psMetadataAlloc ();
+		psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
+		psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+		psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+		psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+		psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+		psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
+
+		psArrayAdd (psfTable, 100, row);
+		psFree (row);
+	    }
+	}
+    }
+    
+    // write an empty FITS segment if we have no PSF information
+    if (psfTable->n == 0) {
+	// XXX this is probably an error (if we have a PSF, how do we have no data?)
+        psFitsWriteBlank (file->fits, header, tableName);
+    } else {
+	psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
+	if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
+	    psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
+	    psFree (tableName);
+	    psFree (psfTable);
+	    psFree (header);
+	    return false;
+	}
+    }
+    psFree (tableName);
+    psFree (psfTable);
+    psFree (header);
+
+    // EXTNAME for residual images
+    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+    if (!rule) {
+	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+	return false;
+    }
+    char *imageName = pmFPAfileNameFromRule (rule, file, view);
+
+    // write the residual images (3D) 
+    header = psMetadataAlloc ();
+    if (psf->residuals == NULL) {
+	// set some header keywords to make it clear there are no residuals?
+        psFitsWriteBlank (file->fits, header, imageName);
+	psFree (imageName);
+        psFree (header);
+        return true;
+    }
+
+    psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
+    psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
+
+    // write the residuals as three planes of the image
+    // this call creates an extension with NAXIS3 = 3
+    if (psf->residuals->Rx) {
+	// this call creates an extension with NAXIS3 = 3
+	psArray *images = psArrayAllocEmpty (3);
+	psArrayAdd (images, 1, psf->residuals->Ro);
+	psArrayAdd (images, 1, psf->residuals->Rx);
+	psArrayAdd (images, 1, psf->residuals->Ry);
+
+	psFitsWriteImageCube (file->fits, header, images, imageName);
+
+	// psFitsWriteImage(file->fits, header, psf->residuals->Ro, 3, imageName);
+	// psFitsUpdateImage(file->fits, psf->residuals->Rx, 0, 0, 1);
+	// psFitsUpdateImage(file->fits, psf->residuals->Ry, 0, 0, 2);
+    } else {
+	// this call creates an extension with NAXIS3 = 1
+	psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, imageName);
+    }
+  
+    psFree (imageName);
+    psFree (header);
+    return true;
+
+    // XXX save the growth curve
+    // XXX save the ApResid fit
+
+# if (0)
+    // build a FITS table of the fit to the Aperture Residuals
+    psArray *apresTable = psArrayAllocEmpty (100);
+    psPolynomial4D *poly = psf->ApTrend;
+    for (int ix = 0; ix < poly->nX; ix++) {
+	for (int iy = 0; iy < poly->nY; iy++) {
+
+	    row = psMetadataAlloc ();
+	    psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+	    psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+	    psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+	    psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+	    psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
+
+	    psArrayAdd (psfTable, 100, row);
+	    psFree (row);
+	}
+    }
+# endif
+}
+
+// XXX add in error handling
+bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    pmFPA *fpa = file->fpa;
+
+    if (view->chip == -1) {
+        pmFPAReadPSFmodel (fpa, view, file, config);
+        return true;
+    }
+
+    if (view->chip >= fpa->chips->n) {
+        return false;
+    }
+    pmChip *chip = fpa->chips->data[view->chip];
+
+    if (view->cell == -1) {
+        pmChipReadPSFmodel (chip, view, file, config);
+        return true;
+    }
+
+    if (view->cell >= chip->cells->n) {
+        return false;
+    }
+    pmCell *cell = chip->cells->data[view->cell];
+
+    if (view->readout == -1) {
+        pmCellReadPSFmodel (cell, view, file, config);
+        return true;
+    }
+
+    if (view->readout >= cell->readouts->n) {
+        return false;
+    }
+    pmReadout *readout = cell->readouts->data[view->readout];
+
+    pmReadoutReadPSFmodel (readout, view, file, config);
+    return true;
+}
+
+// read in all chip-level PSFmodel files for this FPA
+bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < fpa->chips->n; i++) {
+
+        pmChip *chip = fpa->chips->data[i];
+        pmChipReadPSFmodel (chip, view, file, config);
+    }
+    return true;
+}
+
+// read in all cell-level PSFmodel files for this chip
+bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < chip->cells->n; i++) {
+
+        pmCell *cell = chip->cells->data[i];
+        pmCellReadPSFmodel (cell, view, file, config);
+    }
+    return true;
+}
+
+// read in all readout-level PSFmodel files for this cell
+bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+
+    for (int i = 0; i < cell->readouts->n; i++) {
+
+        pmReadout *readout = cell->readouts->data[i];
+        pmReadoutReadPSFmodel (readout, view, file, config);
+    }
+    return true;
+}
+
+// for each Readout (ie, analysed image), we write out: header + table with PSF model parameters, 
+// and header + image for the PSF residual images
+bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+{
+    bool status; 
+    char *rule = NULL;
+    psMetadata *header = NULL;
+
+    // Menu of EXTNAME rules
+    psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+    if (!menu) {
+	psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+	return false;
+    }
+    // EXTNAME for table data
+    rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
+    if (!rule) {
+	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+	return false;
+    }
+    char *tableName = pmFPAfileNameFromRule (rule, file, view);
+    // EXTNAME for residual images
+    rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+    if (!rule) {
+	psError(PS_ERR_UNKNOWN, true, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+	return false;
+    }
+    char *imageName = pmFPAfileNameFromRule (rule, file, view);
+
+    // move fits pointer to table and read header
+    // advance to the table data extension
+    // since we have read the IMAGE header, the TABLE header should exist
+    if (!psFitsMoveExtName (file->fits, tableName)) {
+	psAbort("cannot find data extension %s in %s", tableName, file->filename);
+    }
+
+    // load the PSF model table header
+    header = psFitsReadHeader (NULL, file->fits);
+    if (!header) psAbort("cannot read table header");
+
+    // load the PSF model parameters from the FITS table
+    char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME");
+    pmModelType type = pmModelSetType (modelName);
+
+    bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON");
+
+    // we determine the PSF parameter polynomials from the MD-defined polynomials
+    pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL);
+
+    // check the number of expected parameters
+    int nPar = psMetadataLookupS32 (&status, header, "PSF_NPAR");
+    if (nPar != pmModelParameterCount (psf->type))
+        psAbort("mismatch model par count");
+
+    // load the dimensions of each parameter
+    for (int i = 0; i < nPar; i++) {
+	char name[9];
+	snprintf (name, 9, "PAR%02d_NX", i);
+        int nXorder = psMetadataLookupS32 (&status, header, name);
+	snprintf (name, 9, "PAR%02d_NY", i);
+        int nYorder = psMetadataLookupS32 (&status, header, name);
+	psf->params_NEW->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder);
+    }
+
+    // other required information describing the PSF
+    psf->ApResid   = psMetadataLookupF32 (&status, header, "AP_RESID");
+    psf->dApResid  = psMetadataLookupF32 (&status, header, "AP_ERROR");
+    psf->chisq     = psMetadataLookupF32 (&status, header, "CHISQ");
+    psf->nPSFstars = psMetadataLookupS32 (&status, header, "NSTARS");
+
+    // XXX can we drop this now?
+    psf->skyBias   = psMetadataLookupF32 (&status, header, "SKY_BIAS");
+
+    // read the raw table data
+    psArray *table = psFitsReadTable (file->fits);
+
+    // fill in the matching psf->params entries
+    for (int i = 0; i > table->n; i++) {
+	psMetadata *row = table->data[i];
+	int iPar = psMetadataLookupS32 (&status, row, "MODEL_TERM");
+	int xPow = psMetadataLookupS32 (&status, row, "X_POWER");
+	int yPow = psMetadataLookupS32 (&status, row, "Y_POWER");
+	// XXX sanity check here
+
+        psPolynomial2D *poly = psf->params_NEW->data[iPar];
+
+	poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+	poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
+	poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
+    }
+    psFree (header);
+    psFree (table);
+
+    // move fits pointer to residual image and read header
+    // advance to the table data extension
+    // since we have read the IMAGE header, the TABLE header should exist
+    if (!psFitsMoveExtName (file->fits, imageName)) {
+	psAbort("cannot find data extension %s in %s", imageName, file->filename);
+    }
+
+    header = psFitsReadHeader (NULL, file->fits);
+    int Naxis = psMetadataLookupS32 (&status, header, "NAXIS");
+    if (Naxis != 0) {
+
+	int Nx = psMetadataLookupS32 (&status, header, "NAXIS1");
+	int Ny = psMetadataLookupS32 (&status, header, "NAXIS2");
+	int Nz = psMetadataLookupS32 (&status, header, "NAXIS3");
+
+	int xBin  = psMetadataLookupS32 (&status, header, "XBIN");
+	int yBin  = psMetadataLookupS32 (&status, header, "YBIN");
+
+	int xSize = Nx / xBin;
+	int ySize = Ny / yBin;
+
+	psf->residuals = pmResidualsAlloc (xSize, ySize, xBin, yBin);
+
+	psf->residuals->xCenter = psMetadataLookupS32 (&status, header, "XCENTER");
+	psf->residuals->yCenter = psMetadataLookupS32 (&status, header, "YCENTER");
+
+	psRegion fullImage = {0, 0, 0, 0};
+	psFitsReadImageBuffer(psf->residuals->Ro, file->fits, fullImage, 0); // Desired pixels
+	if (Nz > 1) {
+	    assert (Nz == 3);
+	    psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1); // Desired pixels
+	    psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2); // Desired pixels
+	}
+    }
+
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
+    psFree (psf);
+
+    psFree (tableName);
+    psFree (imageName);
+
+    return true;
+}
+
+/************ old support functions, deprecate? **************/
+
 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf)
 {
@@ -124,81 +608,6 @@
 }
 
-bool pmFPAviewWritePSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    pmFPA *fpa = file->fpa;
-
-    if (view->chip == -1) {
-        pmFPAWritePSFmodel (fpa, view, file, config);
-        return true;
-    }
-
-    if (view->chip >= fpa->chips->n) {
-        return false;
-    }
-    pmChip *chip = fpa->chips->data[view->chip];
-
-    if (view->cell == -1) {
-        pmChipWritePSFmodel (chip, view, file, config);
-        return true;
-    }
-
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        pmCellWritePSFmodel (cell, view, file, config);
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    pmReadoutWritePSFmodel (readout, view, file, config);
-    return true;
-}
-
-// read in all chip-level PSFmodel files for this FPA
-bool pmFPAWritePSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < fpa->chips->n; i++) {
-
-        pmChip *chip = fpa->chips->data[i];
-        pmChipWritePSFmodel (chip, view, file, config);
-    }
-    return true;
-}
-
-// read in all cell-level PSFmodel files for this chip
-bool pmChipWritePSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < chip->cells->n; i++) {
-
-        pmCell *cell = chip->cells->data[i];
-        pmCellWritePSFmodel (cell, view, file, config);
-    }
-    return true;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellWritePSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-
-        pmReadout *readout = cell->readouts->data[i];
-        pmReadoutWritePSFmodel (readout, view, file, config);
-    }
-    return true;
-}
-
 // read in all readout-level Objects files for this cell
-bool pmReadoutWritePSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmReadoutWritePSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
     bool status;
@@ -228,83 +637,6 @@
 }
 
-
-
-bool pmFPAviewReadPSFmodel (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    pmFPA *fpa = file->fpa;
-
-    if (view->chip == -1) {
-        pmFPAReadPSFmodel (fpa, view, file, config);
-        return true;
-    }
-
-    if (view->chip >= fpa->chips->n) {
-        return false;
-    }
-    pmChip *chip = fpa->chips->data[view->chip];
-
-    if (view->cell == -1) {
-        pmChipReadPSFmodel (chip, view, file, config);
-        return true;
-    }
-
-    if (view->cell >= chip->cells->n) {
-        return false;
-    }
-    pmCell *cell = chip->cells->data[view->cell];
-
-    if (view->readout == -1) {
-        pmCellReadPSFmodel (cell, view, file, config);
-        return true;
-    }
-
-    if (view->readout >= cell->readouts->n) {
-        return false;
-    }
-    pmReadout *readout = cell->readouts->data[view->readout];
-
-    pmReadoutReadPSFmodel (readout, view, file, config);
-    return true;
-}
-
-// read in all chip-level PSFmodel files for this FPA
-bool pmFPAReadPSFmodel (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < fpa->chips->n; i++) {
-
-        pmChip *chip = fpa->chips->data[i];
-        pmChipReadPSFmodel (chip, view, file, config);
-    }
-    return true;
-}
-
-// read in all cell-level PSFmodel files for this chip
-bool pmChipReadPSFmodel (pmChip *chip, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < chip->cells->n; i++) {
-
-        pmCell *cell = chip->cells->data[i];
-        pmCellReadPSFmodel (cell, view, file, config);
-    }
-    return true;
-}
-
-// read in all readout-level PSFmodel files for this cell
-bool pmCellReadPSFmodel (pmCell *cell, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
-{
-
-    for (int i = 0; i < cell->readouts->n; i++) {
-
-        pmReadout *readout = cell->readouts->data[i];
-        pmReadoutReadPSFmodel (readout, view, file, config);
-    }
-    return true;
-}
-
 // read in all readout-level Objects files for this cell
-bool pmReadoutReadPSFmodel (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+bool pmReadoutReadPSFmodel_Config (pmReadout *readout, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
 
Index: /trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- /trunk/psModules/src/objects/pmPSFtry.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmPSFtry.c	(revision 13034)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -99,6 +99,4 @@
 {
     bool status;
-    float x;
-    float y;
     int Next = 0;
     int Npsf = 0;
@@ -120,13 +118,10 @@
             return NULL;
         }
-        x = source->peak->x;
-        y = source->peak->y;
-
-        // set temporary object mask and fit object
+
+        // set object mask to define valid pixels
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", PM_MASK_MARK);
+
         // fit model as EXT, not PSF
-
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
         status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT);
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
 
         // exclude the poor fits
@@ -158,11 +153,10 @@
         source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
         source->modelPSF->radiusFit = RADIUS;
-        x = source->peak->x;
-        y = source->peak->y;
-
-        // set the mask and fit the PSF model
-        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
+
+        // set object mask to define valid pixels
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", PM_MASK_MARK);
+
+        // fit the PSF model to the source
         status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF);
-        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
 
         // skip poor fits
Index: /trunk/psModules/src/objects/pmPeaks.c
===================================================================
--- /trunk/psModules/src/objects/pmPeaks.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmPeaks.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-12 18:56:52 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -134,4 +134,6 @@
     tmp->flux = 0;
     tmp->SN = 0;
+    tmp->xf = x;
+    tmp->yf = y;
 
     tmp->type = type;
Index: /trunk/psModules/src/objects/pmSource.c
===================================================================
--- /trunk/psModules/src/objects/pmSource.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSource.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -39,5 +39,7 @@
     psFree(tmp->pixels);
     psFree(tmp->weight);
-    psFree(tmp->mask);
+    psFree(tmp->maskObj);
+    psFree(tmp->maskView);
+    psFree(tmp->modelFlux);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -47,4 +49,5 @@
 }
 
+// free only the pixel data associated with this source
 void pmSourceFreePixels(pmSource *source)
 {
@@ -55,9 +58,13 @@
     psFree (source->pixels);
     psFree (source->weight);
-    psFree (source->mask);
+    psFree (source->maskObj);
+    psFree (source->maskView);
+    psFree (source->modelFlux);
 
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     return;
 }
@@ -76,5 +83,7 @@
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -114,24 +123,29 @@
     // the original values.
     pmSource *source = pmSourceAlloc ();
-    source->peak = psMemIncrRefCounter (in->peak);
-    source->pixels = psMemIncrRefCounter (in->pixels);
-    source->weight = psMemIncrRefCounter (in->weight);
-    source->mask = psMemIncrRefCounter (in->mask);
-    source->moments = psMemIncrRefCounter (in->moments);
-    source->blends = NULL;
-    source->modelPSF = NULL;
-    source->modelEXT = NULL;
+
+    // this is actually the same peak as the original, is this the correct way to handle this?
+    source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
+    source->peak->xf = in->peak->xf;
+    source->peak->yf = in->peak->yf;
+    source->peak->flux = in->peak->flux;
+    source->peak->SN = in->peak->SN;
+
+    // copy the values in the moments structure
+    source->moments  =  pmMomentsAlloc();
+    *source->moments = *in->moments;
+
+    // These images are all views to the parent.  
+    // We want a new view, but pointing at the same pixels.
+    source->pixels   = psImageCopyView (NULL, in->pixels);
+    source->weight   = psImageCopyView (NULL, in->weight);
+    source->maskView = psImageCopyView (NULL, in->maskView);
+
+    // the maskObj is a unique mask array; create a new mask image
+    if (in->maskObj) {
+	source->maskObj = psImageCopy (NULL, in->maskObj, PS_TYPE_MASK);
+    }
+
     source->type = in->type;
     source->mode = in->mode;
-    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
-
-    // default values are NAN
-    source->psfMag = NAN;
-    source->extMag = NAN;
-    source->errMag = NAN;
-    source->apMag  = NAN;
-    source->sky    = NAN;
-    source->skyErr = NAN;
-    source->pixWeight = NAN;
 
     return(source);
@@ -151,8 +165,12 @@
     srcRegion = psRegionForImage (readout->image, srcRegion);
 
-    mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
-    mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
-    mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
-    mySource->region = srcRegion;
+    // these images are subset images of the equivalent parents
+    mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
+    mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
+    mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
+    mySource->region   = srcRegion;
+
+    // the object mask is a copy, and used to define the source pixels
+    mySource->maskObj = psImageCopy (NULL, mySource->maskView, PS_TYPE_MASK);
 
     return true;
@@ -183,5 +201,6 @@
     extend |= (mySource->pixels == NULL);
     extend |= (mySource->weight == NULL);
-    extend |= (mySource->mask == NULL);
+    extend |= (mySource->maskObj == NULL);
+    extend |= (mySource->maskView == NULL);
 
     // extend = true;
@@ -190,12 +209,19 @@
         psFree (mySource->pixels);
         psFree (mySource->weight);
-        psFree (mySource->mask);
-
-        mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
-        mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
-        mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
-        mySource->region = newRegion;
-    }
-
+        psFree (mySource->maskView);
+
+        mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
+        mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
+        mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
+        mySource->region   = newRegion;
+
+	// re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
+	// pixels (eg, with psImageKeepCircle)
+	mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_MASK);
+	
+	// drop the old modelFlux pixels and force the user to re-create
+	psFree (mySource->modelFlux);
+	mySource->modelFlux = NULL;
+    }
     return extend;
 }
@@ -207,4 +233,6 @@
 
 // XXX EAM include a S/N cutoff in selecting the sources?
+// XXX this function should probably accept the values, not a recipe. wrap with a 
+// psphot-specific function which applies the recipe values
 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe)
 {
@@ -458,6 +486,6 @@
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
         inner = psRegionForSquare (source->peak->x, source->peak->y, 2);
-        inner = psRegionForImage (source->mask, inner);
-        int Nsatpix = psImageCountPixelMask (source->mask, inner, PM_MASK_SAT);
+        inner = psRegionForImage (source->maskView, inner);
+        int Nsatpix = psImageCountPixelMask (source->maskView, inner, PM_MASK_SAT);
 
         // saturated star (size consistent with PSF or larger)
@@ -572,5 +600,5 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
@@ -622,5 +650,5 @@
         psF32 *vPix = source->pixels->data.F32[row];
         psF32 *vWgt = source->weight->data.F32[row];
-        psU8  *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];
+        psU8  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.U8[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -719,12 +747,124 @@
 }
 
-pmModel *pmSourceSelectModel (pmSource *source)
-{
+// construct a realization of the source model 
+bool pmSourceCacheModel (pmSource *source) {
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+	
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // in some places (psphotEnsemble), we need a normalized version
+    // in others, we just want the model.  which is more commonly used?
+    // modelFlux always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM);
+    return true;
+}
+
+// should we call pmSourceCacheModel if it does not exist?
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add) {
+
+    bool status;
+
+    if (add) {
+	psTrace ("psphot", 3, "replacing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+	psTrace ("psphot", 3, "removing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    if (source->modelFlux) {
+	// add in the pixels from the modelFlux image
+	int dX = source->modelFlux->col0 - source->pixels->col0;
+	int dY = source->modelFlux->row0 - source->pixels->row0;
+	assert (dX >= 0);
+	assert (dY >= 0);
+	assert (dX + source->modelFlux->numCols <= source->pixels->numCols);
+	assert (dY + source->modelFlux->numRows <= source->pixels->numRows);
+	
+	// modelFlux has unity normalization
+	float Io = model->params->data.F32[PM_PAR_I0];
+	if (mode & PM_MODEL_OP_NORM) {
+	    Io = 1.0;
+	}
+
+	psU8 **mask = NULL;
+	if (source->maskObj) {
+	    mask = source->maskObj->data.U8;
+	}
+
+	// XXX need to respect the source and model masks?
+	for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+	    int oy = iy + dY;
+	    for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+		int ox = ix + dX;
+		if (mask && mask[iy][ix]) continue;
+		float value = Io*source->modelFlux->data.F32[iy][ix];
+		if (add) {
+		    source->pixels->data.F32[oy][ox] += value;
+		} else {
+		    source->pixels->data.F32[oy][ox] -= value;
+		}
+	    }
+	}
+	return true;
+    } 
+
+    if (add) {
+	status = pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    } else {
+	status = pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    }
+
+    return true;
+}
+
+bool pmSourceAdd (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, true);
+    return status;
+}
+
+bool pmSourceSub (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, false);
+    return status;
+}
+
+// given a source, which model is currently appropriate?
+// choose PSF or EXT based on source->type, but fall back on PSF
+// if the EXT model is NULL
+pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
+{
+
+    pmModel *model;
+
+    if (isPSF) {
+      *isPSF = false;
+    }
+
     switch (source->type) {
     case PM_SOURCE_TYPE_STAR:
-        return source->modelPSF;
+        model = source->modelPSF;
+        if (model == NULL)
+            return NULL;
+        if (isPSF) {
+            *isPSF = true;
+	}
+        return model;
 
     case PM_SOURCE_TYPE_EXTENDED:
-        return source->modelEXT;
+        model = source->modelEXT;
+        if (!model && source->modelPSF) {
+            if (isPSF) {
+                *isPSF = true;
+	    }
+	    return source->modelPSF;
+        }
+        return (model);
+        break;
 
     default:
Index: /trunk/psModules/src/objects/pmSource.h
===================================================================
--- /trunk/psModules/src/objects/pmSource.h	(revision 13033)
+++ /trunk/psModules/src/objects/pmSource.h	(revision 13034)
@@ -3,6 +3,6 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -13,21 +13,4 @@
 /// @addtogroup Objects Object Detection / Analysis Functions
 /// @{
-
-/**
- * In the object analysis process, we will use specific mask values to mark the
- * image pixels. The following structure defines the relevant mask values.
- *
- * XXX: This is probably a bad solution: we will want to set mask values
- * outside of the PSPHOT code.  Perhaps we can set up a registered set of mask
- * values with specific meanings that other functions can add to or define?
- 
- * XXX We will only use the PM_MASK_xxx mask values
-typedef enum {
-    PM_SOURCE_MASK_CLEAR     = 0x00,
-    PM_SOURCE_MASK_INVALID   = 0x01,
-    PM_SOURCE_MASK_SATURATED = 0x02,
-    PM_SOURCE_MASK_MARKED    = 0x08,
-} psphotMaskValues;
-*/
 
 /** pmSourceType enumeration
@@ -75,22 +58,24 @@
 typedef struct
 {
-    const int id;   ///< Unique ID for object
-    pmPeak *peak;   ///< Description of peak pixel.
-    psImage *pixels;   ///< Rectangular region including object pixels.
-    psImage *weight;   ///< Image variance.
-    psImage *mask;   ///< Mask which marks pixels associated with objects.
-    pmMoments *moments;   ///< Basic moments measure for the object.
-    pmModel *modelPSF;   ///< PSF Model fit (parameters and type)
-    pmModel *modelEXT;   ///< EXT (floating) Model fit (parameters and type).
-    pmSourceType type;   ///< Best identification of object.
-    pmSourceMode mode;   ///< Best identification of object.
+    const int id;			///< Unique ID for object
+    pmPeak *peak;			///< Description of peak pixel.
+    psImage *pixels;			///< Rectangular region including object pixels.
+    psImage *weight;			///< Image variance.
+    psImage *modelFlux;			///< cached copy of the model for this source
+    psImage *maskObj;			///< unique mask for this object which marks included pixels associated with objects.
+    psImage *maskView;			///< view into global image mask for this object region
+    pmMoments *moments;			///< Basic moments measure for the object.
+    pmModel *modelPSF;			///< PSF Model fit (parameters and type)
+    pmModel *modelEXT;			///< EXT (floating) Model fit (parameters and type).
+    pmSourceType type;			///< Best identification of object.
+    pmSourceMode mode;			///< Best identification of object.
     psArray *blends;
-    float psfMag;   ///< calculated from flux in modelPsf
-    float extMag;   ///< calculated from flux in modelEXT
-    float errMag;   ///< error in psfMag OR extMag (depending on type)
-    float apMag;   ///< apMag corresponding to psfMag or extMag (depending on type)
-    float pixWeight;   // model-weighted coverage of valid pixels
-    psRegion region;   // area on image covered by selected pixels
-    float sky, skyErr;   //?< The sky and its error at the center of the object
+    float psfMag;			///< calculated from flux in modelPsf
+    float extMag;			///< calculated from flux in modelEXT
+    float errMag;			///< error in psfMag OR extMag (depending on type)
+    float apMag;	       ///< apMag corresponding to psfMag or extMag (depending on type)
+    float pixWeight;			// model-weighted coverage of valid pixels
+    psRegion region;			// area on image covered by selected pixels
+    float sky, skyErr;			//?< The sky and its error at the center of the object
 }
 pmSource;
@@ -225,7 +210,11 @@
 );
 
-
-// select the model used for this source
-pmModel *pmSourceSelectModel (pmSource *source);
+pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
+
+bool pmSourceAdd (pmSource *source, pmModelOpMode mode);
+bool pmSourceSub (pmSource *source, pmModelOpMode mode);
+
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add);
+bool pmSourceCacheModel (pmSource *source);
 
 /// @}
Index: /trunk/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceFitModel.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceFitModel.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -57,5 +57,5 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
@@ -77,5 +77,5 @@
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
             // skip masked points
-            if (source->mask->data.U8[i][j]) {
+            if (source->maskObj->data.U8[i][j]) {
                 continue;
             }
@@ -334,5 +334,5 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->weight, false);
 
@@ -354,5 +354,5 @@
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
             // skip masked points
-            if (source->mask->data.U8[i][j]) {
+            if (source->maskObj->data.U8[i][j]) {
                 continue;
             }
Index: /trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -402,4 +402,5 @@
 
 // Given a FITS file pointer, read the table of object data
+// XXX add in error handling
 bool pmFPAviewReadObjects (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
Index: /trunk/psModules/src/objects/pmSourceIO_CMF.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_CMF.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_CMF.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -56,5 +56,7 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;
Index: /trunk/psModules/src/objects/pmSourceIO_CMP.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_CMP.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_CMP.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -107,5 +107,7 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;
Index: /trunk/psModules/src/objects/pmSourceIO_OBJ.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_OBJ.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_OBJ.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -61,5 +61,7 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;
Index: /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_0.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -63,5 +63,7 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
 
         if (model != NULL) {
Index: /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_SMPDATA.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -62,5 +62,7 @@
     for (i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model != NULL) {
 	    PAR = model->params->data.F32;
Index: /trunk/psModules/src/objects/pmSourceIO_SX.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceIO_SX.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceIO_SX.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -57,5 +57,7 @@
     for (int i = 0; i < sources->n; i++) {
         pmSource *source = (pmSource *) sources->data[i];
-        pmModel *model = pmSourceSelectModel (source);
+
+	// no difference between PSF and non-PSF model
+        pmModel *model = pmSourceGetModel (NULL, source);
         if (model == NULL)
             continue;
Index: /trunk/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- /trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 13034)
@@ -3,6 +3,6 @@
  *  @author EAM, IfA; GLG, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -106,6 +106,8 @@
     // replace source flux
     // XXX need to be certain which model and size of mask for prior subtraction
+    // XXX full model or just analytical?
+    // XXX use pmSourceAdd instead?
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelAdd (source->pixels, source->mask, model, false, false);
+        pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -180,14 +182,16 @@
     // set aperture mask circle to model radius
     // XXX use a different radius for apertures and fits...
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK);
+    // XXX can I remove this?  the source should have the mask defined when it is constructed or
+    // when the fit / aperture radius is changed...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "OR", PM_MASK_MARK);
 
     // measure the weight of included pixels
     // XXX is this supposed to use the weight or the flux?
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
-        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->mask);
+        pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj);
     }
 
     // measure object aperture photometry
-    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->mask);
+    status = pmSourcePhotometryAper  (&source->apMag, model, flux, source->maskObj);
 
     // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS
@@ -206,9 +210,11 @@
 
     // unmask aperture
-    psImageKeepCircle (source->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
+    // XXX can I remove this?  this will probably break things downstream...
+    psImageKeepCircle (source->maskObj, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
 
     // if source was originally subtracted, re-subtract object, leave local sky
+    // XXX replace with pmSourceSub...
     if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
-        pmModelSub (source->pixels, source->mask, model, false, false);
+        pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
     }
 
@@ -363,4 +369,5 @@
 }
 
+# if (0)
 double pmSourceCrossProduct (const pmSource *Mi,
                              const pmSource *Mj,
@@ -379,6 +386,6 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -435,6 +442,6 @@
     const psImage *Wi = Mi->weight;
 
-    const psImage *Ti = Mi->mask;
-    const psImage *Tj = Mj->mask;
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
 
     Xs = PS_MAX (Pi->col0, Pj->col0);
@@ -483,5 +490,5 @@
     const psImage *Pi = Mi->pixels;
     const psImage *Wi = Mi->weight;
-    const psImage *Ti = Mi->mask;
+    const psImage *Ti = Mi->maskObj;
 
     // note that this is addressing the same image pixels,
@@ -521,4 +528,5 @@
     return flux;
 }
+# endif
 
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight)
@@ -543,43 +551,158 @@
 }
 
-// given a source, which model is currently appropriate?
-// choose PSF or EXT based on source->type, but fall back on PSF
-// if the EXT model is NULL
-pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
-{
-
-    pmModel *model;
-
-    switch (source->type) {
-    case PM_SOURCE_TYPE_STAR:
-        model = source->modelPSF;
-        if (model == NULL)
-            return NULL;
-        if (isPSF)
-            *isPSF = true;
-        return model;
-
-    case PM_SOURCE_TYPE_EXTENDED:
-        model = source->modelEXT;
-        if (isPSF)
-            *isPSF = false;
-        if (model == NULL) {
-            model = source->modelPSF;
-            if (isPSF)
-                *isPSF = true;
-        }
-        if (model == NULL) {
-            if (isPSF)
-                *isPSF = FALSE;
-            return NULL;
-        }
-        return (model);
-        break;
-
-    default:
-        if (isPSF)
-            *isPSF = false;
-        return NULL;
-    }
-    return NULL;
-}
+
+double pmSourceModelWeight(const pmSource *Mi,
+                      int term,
+                      const bool unweighted_sum) // should the cross product be weighted?
+{
+    double flux = 0, wt = 0, factor = 0;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Wi = Mi->weight;
+    const psImage *Ti = Mi->maskObj;
+
+    for (int yi = 0; yi < Pi->numRows; yi++) {
+        for (int xi = 0; xi < Pi->numCols; xi++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (!unweighted_sum) {
+                wt = Wi->data.F32[yi][xi];
+                if (wt == 0)
+                    continue;
+            }
+
+            switch (term) {
+            case 0:
+                factor = 1;
+                break;
+            case 1:
+                factor = xi + Pi->col0;
+                break;
+            case 2:
+                factor = yi + Pi->row0;
+                break;
+            default:
+                psAbort("invalid term for pmSourceWeight");
+            }
+
+            if (unweighted_sum) {
+                flux += (factor * Pi->data.F32[yi][xi]);
+            } else {
+                flux += (factor * Pi->data.F32[yi][xi]) / wt;
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceModelDotModel (const pmSource *Mi,
+			      const pmSource *Mj,
+			      const bool unweighted_sum) // should the cross product be weighted?
+{
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->modelFlux;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+	    // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+
+double pmSourceDataDotModel (const pmSource *Mi,
+                             const pmSource *Mj,
+                             const bool unweighted_sum) // should the cross product be weighted?
+{
+
+    int Xs, Xe, Ys, Ye;
+    int xi, xj, yi, yj;
+    int xIs, xJs, yIs, yJs;
+    int xIe, yIe;
+    double flux, wt;
+
+    const psImage *Pi = Mi->pixels;
+    const psImage *Pj = Mj->modelFlux;
+
+    const psImage *Wi = Mi->weight;
+
+    const psImage *Ti = Mi->maskObj;
+    const psImage *Tj = Mj->maskObj;
+
+    Xs = PS_MAX (Pi->col0, Pj->col0);
+    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
+
+    Ys = PS_MAX (Pi->row0, Pj->row0);
+    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
+
+    xIs = Xs - Pi->col0;
+    xJs = Xs - Pj->col0;
+    yIs = Ys - Pi->row0;
+    yJs = Ys - Pj->row0;
+
+    xIe = Xe - Pi->col0;
+    yIe = Ye - Pi->row0;
+
+    // note that weight is addressing the same image pixels,
+    flux = 0;
+    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
+        for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
+            if (Ti->data.U8[yi][xi])
+                continue;
+            if (Tj->data.U8[yj][xj])
+                continue;
+
+	    // XXX skip the nonsense weight pixels?
+            if (unweighted_sum) {
+                flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]);
+            } else {
+                wt = Wi->data.F32[yi][xi];
+                if (wt > 0) {
+                    flux += (Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]) / wt;
+                }
+            }
+        }
+    }
+    return flux;
+}
+
Index: /trunk/psModules/src/objects/pmSourcePhotometry.h
===================================================================
--- /trunk/psModules/src/objects/pmSourcePhotometry.h	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourcePhotometry.h	(revision 13034)
@@ -4,6 +4,6 @@
  * @author EAM, IfA; GLG, MHPCC
  *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
+ * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-04-26 01:20:29 $
  * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -50,11 +50,16 @@
 bool pmSourceMagnitudesInit (psMetadata *config);
 bool pmSourceMagnitudes (pmSource *source, pmPSF *psf, pmSourcePhotometryMode mode);
-double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
-double pmSourceCrossWeight(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
 bool pmSourcePixelWeight (float *pixWeight, pmModel *model, psImage *image, psImage *mask);
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight);
-pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source);
 
-double pmSourceWeight(const pmSource *Mi, int term, const bool unweighted_sum);
+
+double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+double pmSourceModelDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum);
+
+// retire these:
+// double pmSourceCrossProduct(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+// double pmSourceCrossWeight(const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum);
+// double pmSourceWeight(const pmSource *Mi, int term, const bool unweighted_sum);
 
 /// @}
Index: /trunk/psModules/src/objects/pmSourceSky.c
===================================================================
--- /trunk/psModules/src/objects/pmSourceSky.c	(revision 13033)
+++ /trunk/psModules/src/objects/pmSourceSky.c	(revision 13034)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49,5 +49,5 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_IMAGE_NON_NULL(source->pixels, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_INT_POSITIVE(Radius, false);
@@ -61,5 +61,5 @@
 
     psImage *image = source->pixels;
-    psImage *mask  = source->mask;
+    psImage *mask  = source->maskObj;
     pmPeak *peak  = source->peak;
     psRegion srcRegion;
@@ -100,5 +100,5 @@
     PS_ASSERT_PTR_NON_NULL(source, false);
     PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
-    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->maskObj, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_INT_POSITIVE(Radius, false);
@@ -112,5 +112,5 @@
 
     psImage *image = source->weight;
-    psImage *mask  = source->mask;
+    psImage *mask  = source->maskObj;
     pmPeak *peak  = source->peak;
     psRegion srcRegion;
Index: /trunk/psModules/test/objects/tap_pmSourceFitModel.c
===================================================================
--- /trunk/psModules/test/objects/tap_pmSourceFitModel.c	(revision 13033)
+++ /trunk/psModules/test/objects/tap_pmSourceFitModel.c	(revision 13034)
@@ -136,5 +136,5 @@
 
     // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
-    pmModelAdd (source->pixels, source->mask, source->modelEXT, false, false);
+    pmModelAdd (source->pixels, source->mask, source->modelEXT, PM_MODEL_OP_FULL);
     int npix = 0;
     for (int j = 0; j < source->pixels->numRows; j++) {
Index: /trunk/psModules/test/objects/tap_pmSourceFitModel_Delta.c
===================================================================
--- /trunk/psModules/test/objects/tap_pmSourceFitModel_Delta.c	(revision 13033)
+++ /trunk/psModules/test/objects/tap_pmSourceFitModel_Delta.c	(revision 13034)
@@ -54,5 +54,5 @@
 
     // create an image with the model, and add noise: gain is 1, subtracted sky is 100, readnoise is 5
-    pmModelAdd (source->pixels, source->mask, source->modelEXT, false, false);
+    pmModelAdd (source->pixels, source->mask, source->modelEXT, PM_MODEL_OP_FULL);
     int npix = 0;
     for (int j = 0; j < source->pixels->numRows; j++) {