Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c
===================================================================
--- trunk/psModules/src/imcombine/pmPSFEnvelope.c	(revision 15834)
+++ trunk/psModules/src/imcombine/pmPSFEnvelope.c	(revision 15836)
@@ -32,14 +32,12 @@
 
 #define TESTING                         // Enable test output
-#define PEAK_FLUX 1.0e3                 // peak flux for each source
+#define PEAK_FLUX 1.0e3                 // Peak flux for each source
+#define SKY_VALUE 0.0e0                 // Sky value for fake image
 #define WEIGHT_VAL 1.0e0                // Weighting for image
-#define PSF_STATS PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV // Statistics options for measuring PSF
+#define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF
 
 
 
 // XXX To do:
-//
-// * PSF representation doesn't normalise (x,y) positions in polynomials: need to find out how I can have a
-// star at (x1,y1) in the image that is "really" supposed to be at (x2,y2).
 //
 // * PSF variation when only a portion of the image is present (e.g., the edge of an FPA overlapping a
@@ -60,5 +58,6 @@
                      int instances, // Number of instances per dimension
                      int radius,        // Radius of each PSF
-                     const char *modelName // Name of PSF model to use
+                     const char *modelName,// Name of PSF model to use
+                     int xOrder, int yOrder // Order for PSF variation fit
                      )
 {
@@ -69,4 +68,6 @@
     PS_ASSERT_INT_POSITIVE(radius, NULL);
     PS_ASSERT_STRING_NON_EMPTY(modelName, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(xOrder, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(yOrder, NULL);
 
     float xOrigSpacing = (float)(numCols - 2 * radius) / (float)(instances - 1); // Spacing between instances
@@ -110,5 +111,5 @@
     // Generate fake images with each PSF, and take the envelope
     psImage *envelope = psImageAlloc(fakeSize, fakeSize, PS_TYPE_F32); // Image with envelope of PSFs
-    psImageInit(envelope, 0.0);
+    psImageInit(envelope, SKY_VALUE);
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
     for (int i = 0; i < inputs->n; i++) {
@@ -135,21 +136,5 @@
             float y = source->peak->yf + yOffset->data.S32[j]; // y coordinate of source
 
-#if 0
-            psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR,
-                                                                               fakeRO->image, NULL, NULL, 0,
-                                                                               0.0, 0.0, 1, 0, 0.0);
-            double flux;                 // Flux of peak
-            if (!psImageInterpolate(&flux, NULL, NULL, x, y, interp)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate source image.");
-                psFree(envelope);
-                psFree(yOffset);
-                psFree(xOffset);
-                psFree(fakes);
-                return NULL;
-            }
-#else
             double flux = fakeRO->image->data.F32[(int)y][(int)x];
-#endif
-
             float norm = PEAK_FLUX / flux; // Normalisation for source
             psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region
@@ -161,6 +146,4 @@
             psFree(subEnv);
         }
-
-
 
 #ifdef TESTING
@@ -214,4 +197,7 @@
         }
     }
+    psFree(xOffset);
+    psFree(yOffset);
+
     readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(1.0, PS_TYPE_F32));
     readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
@@ -227,8 +213,9 @@
 #endif
 
+    // Reset the sources to point to the new pixels, and measure the moments in preparation for PSF fitting
     for (int i = 0; i < numFakes; i++) {
         pmSource *source = fakes->data[i]; // Fake source
-        float x = source->peak->xf;// + xOffset->data.S32[i]; // x coordinates of source
-        float y = source->peak->yf;// + yOffset->data.S32[i]; // y coordinates of source
+        float x = source->peak->xf;     // x coordinates of source
+        float y = source->peak->yf;     // y coordinates of source
 
         psFree(source->pixels);
@@ -244,61 +231,26 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
             psFree(readout);
-            psFree(yOffset);
-            psFree(xOffset);
             psFree(fakes);
             return NULL;
         }
-
-#if 0
-        // Need to do some magic so that pmSourceMoments behaves --- the only offsets it understands are the
-        // col0,row0 in the psImage.
-        int oldCol0 = source->pixels->col0;
-        int oldRow0 = source->pixels->row0;
-        int col0 = oldCol0 - xOffset->data.S32[i];
-        int row0 = oldRow0 - yOffset->data.S32[i];
-
-        P_PSIMAGE_SET_COL0(source->pixels, col0);
-        P_PSIMAGE_SET_ROW0(source->pixels, row0);
-        P_PSIMAGE_SET_COL0(source->weight, col0);
-        P_PSIMAGE_SET_ROW0(source->weight, row0);
-        P_PSIMAGE_SET_COL0(source->maskObj, col0);
-        P_PSIMAGE_SET_ROW0(source->maskObj, row0);
-#endif
 
         if (!pmSourceMoments(source, radius)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source.");
             psFree(readout);
-            psFree(yOffset);
-            psFree(xOffset);
             psFree(fakes);
             return NULL;
         }
-
-#if 0
-        P_PSIMAGE_SET_COL0(source->pixels, oldCol0);
-        P_PSIMAGE_SET_ROW0(source->pixels, oldRow0);
-        P_PSIMAGE_SET_COL0(source->weight, oldCol0);
-        P_PSIMAGE_SET_ROW0(source->weight, oldRow0);
-        P_PSIMAGE_SET_COL0(source->maskObj, oldCol0);
-        P_PSIMAGE_SET_ROW0(source->maskObj, oldRow0);
-#endif
-    }
-
-    psMemCheckCorruption(stderr, true);
-
-
+    }
+
+    // Don't assume Poisson errors
     pmPSFOptions *options = pmPSFOptionsAlloc(); // Options for fitting a PSF
-    // Don't assume Poisson errors --- our fake system isn't Poisson.
     options->poissonErrorsPhotLMM = false;
     options->poissonErrorsPhotLin = false;
     options->poissonErrorsParams = false;
     options->stats = psStatsAlloc(PSF_STATS);
-
     options->radius = radius;
-
     options->psfTrendMode = PM_TREND_MAP;
-    options->psfTrendNx = 1;
-    options->psfTrendNy = 1;
-
+    options->psfTrendNx = xOrder;
+    options->psfTrendNy = yOrder;
     options->psfFieldNx = numCols;
     options->psfFieldNy = numRows;
@@ -313,6 +265,4 @@
         psError(PS_ERR_UNKNOWN, false, "Unable to fit PSF model to PSF envelope.");
         psFree(readout);
-        psFree(yOffset);
-        psFree(xOffset);
         psFree(fakes);
         return NULL;
@@ -334,5 +284,5 @@
 
         pmReadout *generated = pmReadoutAlloc(NULL); // Generated image
-        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius,
+        pmReadoutFakeFromSources(generated, numCols, numRows, fakes, NULL, NULL, psf, NAN, radius,
                                  false);
         {
@@ -341,5 +291,5 @@
             psFitsClose(fits);
         }
-        psBinaryOp(generated->image, generated->image, "-", envelope);
+        psBinaryOp(generated->image, generated->image, "-", readout->image);
         {
             psFits *fits = psFitsOpen("psf_field_resid.fits", "w");
@@ -352,7 +302,6 @@
 
     psFree(yOffset);
-    psFree(xOffset);
-    psFree(fakes);
     psFree(readout);
+
     return psf;
 }