Index: /trunk/psphot/src/psphotModelWithPSF.c
===================================================================
--- /trunk/psphot/src/psphotModelWithPSF.c	(revision 14346)
+++ /trunk/psphot/src/psphotModelWithPSF.c	(revision 14347)
@@ -185,23 +185,34 @@
     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
 
+    psImageInit (pcm->model, 0.0);
+    for (int n = 0; n < params->n; n++) {
+	if (!pcm->dmodels->data[n]) continue;
+	psImageInit (pcm->dmodels->data[n], 0.0);
+    }
+
     // fill in the coordinate and value entries
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
+
+	    // XXX can we skip some of the data points where the model
+	    // is not going to be fitted??
+
             // skip masked points
 	    // XXX probably should not skipped masked points: 
 	    // XXX skip if convolution of unmasked pixels will not see this pixel
-            if (source->maskObj->data.U8[i][j]) {
-                continue;
-            }
+            // if (source->maskObj->data.U8[i][j]) {
+	    // continue;
+	    // }
+
             // skip zero-weight points
 	    // XXX why is this not masked?
-            if (source->weight->data.F32[i][j] == 0) {
-                continue;
-            }
+            // if (source->weight->data.F32[i][j] == 0) {
+	    // continue;
+	    // }
             // skip nan value points
 	    // XXX why is this not masked?
-            if (!isfinite(source->pixels->data.F32[i][j])) {
-                continue;
-            }
+            // if (!isfinite(source->pixels->data.F32[i][j])) {
+	    // continue;
+	    // }
 
             // Convert i/j to image space:
@@ -221,5 +232,4 @@
     psFree(deriv);
 
-    psphotSaveImage (NULL, pcm->model, "model1.fits");
 
     // convolve model and dmodel arrays with PSF
@@ -233,4 +243,6 @@
 
     // XXX TEST : SAVE IMAGES
+# if (1) 
+    psphotSaveImage (NULL, psf->image, "psf.fits");
     psphotSaveImage (NULL, pcm->model, "model.fits");
     psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits");
@@ -238,5 +250,5 @@
     psphotSaveImage (NULL, source->maskObj, "mask.fits");
     psphotSaveImage (NULL, source->weight, "weight.fits");
-    exit (0);
+# endif
 
     // zero alpha and beta for summing below
@@ -330,17 +342,17 @@
 
     // Allocate storage images for raw model and derivative images
-    pcm->model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    pcm->model = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     pcm->dmodels = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      pcm->dmodels->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+      pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     }
 
     // Allocate storage images for convolved model and derivative images
-    pcm->modelConv = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    pcm->modelConv = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     pcm->dmodelsConv = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
       if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      pcm->dmodelsConv->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+      pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     }
 
@@ -350,5 +362,6 @@
 psImage *pmPCMDataSaveImage (pmPCMData *pcm) {
 
-    psImage *model = psImageCopy (NULL, pcm->model, PS_TYPE_F32);
+    psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32);
+   
     return model;
 }