Index: trunk/psModules/src/objects/pmModel.c
===================================================================
--- trunk/psModules/src/objects/pmModel.c	(revision 9810)
+++ trunk/psModules/src/objects/pmModel.c	(revision 12949)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-31 19:38:44 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-21 19:47:14 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21,4 +21,5 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmResiduals.h"
 #include "pmModel.h"
 
@@ -49,5 +50,7 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
+    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
 
     tmp->type = type;
@@ -57,4 +60,5 @@
     tmp->radiusFit = 0;
     tmp->status = PM_MODEL_UNTRIED;
+    tmp->residuals = NULL;		// XXX should the model own this memory?
 
     psS32 Nparams = pmModelParameterCount(type);
@@ -72,5 +76,4 @@
     }
 
-    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
     return(tmp);
@@ -88,4 +91,7 @@
     new->status = model->status;
     new->radiusFit = model->radiusFit;
+
+    // XXX note that model->residuals is just a reference
+    new->residuals = model->residuals;
 
     for (int i = 0; i < new->params->n; i++) {
@@ -147,8 +153,37 @@
     psF32 skyValue = params->data.F32[0];
     psF32 pixelValue;
-
-    for (psS32 i = 0; i < image->numRows; i++) {
-        for (psS32 j = 0; j < image->numCols; j++) {
-            if ((mask != NULL) && mask->data.U8[i][j])
+    
+    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) {
+	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++) {
+            if ((mask != NULL) && mask->data.U8[iy][ix])
                 continue;
 
@@ -156,9 +191,9 @@
             // 'center' option changes meaning of i,j:
             if (center) {
-                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
-                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
+                imageCol = ix - 0.5*image->numCols + xCenter;
+                imageRow = iy - 0.5*image->numRows + yCenter;
             } else {
-                imageCol = j + image->col0;
-                imageRow = i + image->row0;
+                imageCol = ix + image->col0;
+                imageRow = iy + image->row0;
             }
 
@@ -173,16 +208,38 @@
             }
 
+	    // 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)) {
+		    double flux = Fo + xCenter*Fx + yCenter*Fy;
+		    pixelValue += Io*flux;
+		}
+	    }
 
             // add or subtract the value
-            if (add
-               ) {
-                image->data.F32[i][j] += pixelValue;
-            }
-            else {
-                image->data.F32[i][j] -= pixelValue;
+            if (add) {
+                image->data.F32[iy][ix] += pixelValue;
+            } else {
+                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);