Index: /branches/eam_02_branch/psModules/src/objects/pmSource.c
===================================================================
--- /branches/eam_02_branch/psModules/src/objects/pmSource.c	(revision 12953)
+++ /branches/eam_02_branch/psModules/src/objects/pmSource.c	(revision 12954)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.25.4.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 00:03:31 $
+ *  @version $Revision: 1.25.4.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-22 22:23:20 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -733,2 +733,103 @@
     return NULL;
 }
+
+// construct a realization of the source model 
+bool pmSourceCacheModel (pmSource *source) {
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) continue;  // model must be defined
+	
+    // if we already have a cached image, re-use that memory
+    if (source->modelFlux) { 
+	psImageRecycle (source->modelFlux);
+    } else {
+	source->modelFlux = psImageCopy (source->pixels);
+    }
+    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?
+    pmModelAdd (source->modelFlux, source->mask, model, PM_MODEL_ADD_FULL | PM_MODEL_ADD_NORM);
+    return true;
+}
+
+// XXX should this API take an operator (+/-)?
+// XXX should we specify norm, etc?
+bool pmSourceAdd (pmSource *source) {
+
+    if (source->modelFlux) {
+	// add in the pixels from the modelFlux image
+	psImageSubsetAdd (source->pixels, source->modelFlux);
+    } else {
+	pmModel *model = pmSourceGetModel (NULL, source);
+	if (model == NULL) continue;  // model must be defined
+	pmModelAdd (source->pixels, source->mask, model, PM_MODEL_ADD_FULL);
+    }	
+
+    // we would need to save X,Y somehow (or use the peak?)
+    // psTrace ("psphot", 3, "replacing object at %f,%f\n", 
+    // model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    return true;
+}
+
+bool pmSourceSub (pmSource *source) {
+
+    if (source->modelFlux) {
+	// add in the pixels from the modelFlux image
+	psImageSubsetSub (source->pixels, source->modelFlux);
+    } else {
+	pmModel *model = pmSourceGetModel (NULL, source);
+	if (model == NULL) continue;  // model must be defined
+	pmModelSub (source->pixels, source->mask, model, PM_MODEL_ADD_FULL);
+    }	
+    return true;
+}
+
+// XXX merge this function with psImageOverlay, which is probably mid-defined
+bool psImageSubsetAdd (psImage *input, psImage *subset) {
+
+    // col0,row0 in both images refer to parent coordinates
+
+    int dX = subset->col0 - input->col0;
+    int dY = subset->row0 - input->row0;
+    assert (dX >= 0);
+    assert (dY >= 0);
+    assert (dX + subset->numCols < input->numCols);
+    assert (dY + subset->numRows < input->numRows);
+    assert (input->type.type == PS_TYPE_F32);
+    assert (subset->type.type == PS_TYPE_F32);
+
+    for (int iy = 0; iy < subset->numRows; iy++) {
+	oy = iy + dY;
+	for (int ix = 0; ix < subset->numCols; ix++) {
+	    ox = ix + dX;
+	    input->data.F32[oy][ox] += subset->data.F32[iy][ix];
+	}
+    }
+    return true;
+}
+
+// XXX merge this function with psImageOverlay, which is probably mid-defined
+bool psImageSubsetSub (psImage *input, psImage *subset) {
+
+    // col0,row0 in both images refer to parent coordinates
+
+    int dX = subset->col0 - input->col0;
+    int dY = subset->row0 - input->row0;
+    assert (dX >= 0);
+    assert (dY >= 0);
+    assert (dX + subset->numCols < input->numCols);
+    assert (dY + subset->numRows < input->numRows);
+    assert (input->type.type == PS_TYPE_F32);
+    assert (subset->type.type == PS_TYPE_F32);
+
+    for (int iy = 0; iy < subset->numRows; iy++) {
+	oy = iy + dY;
+	for (int ix = 0; ix < subset->numCols; ix++) {
+	    ox = ix + dX;
+	    input->data.F32[oy][ox] -= subset->data.F32[iy][ix];
+	}
+    }
+    return true;
+}
