Index: /branches/pap_branch_20090108/ppStack/src/ppStackLoop.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackLoop.c	(revision 21105)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackLoop.c	(revision 21106)
@@ -348,4 +348,16 @@
         psFree(fileIter);
 
+        // Generate target PSF
+        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
+        psFree(psfs);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sourceLists);
+            psFree(view);
+            return false;
+        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
+                         "Target PSF for stack", targetPSF);
+
         // Zero point calibration
         sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
@@ -354,13 +366,4 @@
             psFree(sourceLists);
             psFree(targetPSF);
-            return false;
-        }
-
-        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
-        psFree(psfs);
-        if (!targetPSF) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
-            psFree(sourceLists);
-            psFree(view);
             return false;
         }
@@ -450,6 +453,11 @@
         psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
         psTimerStart("PPSTACK_MATCH");
+
+
+        // XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
+        // Passing common set of sources, rather than individual
+
         if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
-                          sourceLists->data[i], targetPSF, rng, config)) {
+                          sourceLists->data[0], targetPSF, rng, config)) {
             psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
             inputMask->data.U8[i] = PPSTACK_MASK_MATCH;
Index: /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21105)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackMatch.c	(revision 21106)
@@ -47,4 +47,20 @@
 #endif
 
+// Get coordinates from a source
+static void coordsFromSource(float *x, float *y, const pmSource *source)
+{
+    assert(x && y);
+    assert(source);
+
+    if (source->modelPSF) {
+        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+        *x = source->peak->xf;
+        *y = source->peak->yf;
+    }
+    return;
+}
+
 static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
                                    int exclusion // Exclusion zone, pixels
@@ -64,6 +80,5 @@
             continue;
         }
-        x->data.F32[numGood] = source->peak->xf;
-        y->data.F32[numGood] = source->peak->yf;
+        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
         numGood++;
     }
@@ -80,6 +95,9 @@
             continue;
         }
-        coords->data.F64[0] = source->peak->xf;
-        coords->data.F64[1] = source->peak->yf;
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        coords->data.F64[0] = xSource;
+        coords->data.F64[1] = ySource;
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
@@ -120,5 +138,8 @@
         }
 
-        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        int x = xSource + 0.5, y = ySource + 0.5; // Coordinates of source
 
         int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
@@ -319,7 +340,7 @@
                 pmSource *source = pmSourceAlloc();     // Source
                 sources->data[0] = source;
-                source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE);
+                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
                 source->psfMag = -13.0;
-                pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true);
+                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
                 float sum = 0.0;
                 for (int y = 0; y < test->image->numRows; y++) {
@@ -351,24 +372,30 @@
             }
 
+
+#if 0
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
+#endif
+
 
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "fake_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "real_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -410,12 +437,14 @@
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "conv_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "diff_%03d.fits", numInput);
@@ -423,5 +452,5 @@
                 psFree(name);
                 psBinaryOp(fake->image, output->image, "-", fake->image);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -582,4 +611,5 @@
 #define RADIUS 10                       // Radius of photometry
 #define MIN_ERR 0.05                    // Minimum photometric error, mag
+#define MAX_MAG -13                     // Maximum magnitude for source
 
     // Ensure the normalisation is correct
@@ -595,9 +625,10 @@
         pmSource *source = sources->data[i]; // Source of interest
         if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
-            source->errMag > MIN_ERR) {
+            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
             continue;
         }
 
-        float xSrc = source->peak->xf, ySrc = source->peak->yf; // Source coordinates
+        float xSrc, ySrc;              // Source coordinates
+        coordsFromSource(&xSrc, &ySrc, source);
         int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
         int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
@@ -635,9 +666,10 @@
 #ifdef TESTING
     {
+        pmHDU *hdu = pmHDUFromCell(readout->parent);
         psString name = NULL;
         psStringAppend(&name, "convolved_%03d.fits", numInput);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
-        psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
         psFitsClose(fits);
     }
Index: /branches/pap_branch_20090108/ppStack/src/ppStackSources.c
===================================================================
--- /branches/pap_branch_20090108/ppStack/src/ppStackSources.c	(revision 21105)
+++ /branches/pap_branch_20090108/ppStack/src/ppStackSources.c	(revision 21106)
@@ -9,4 +9,11 @@
 //#define TESTING                         // Enable debugging output
 
+// Size of fake image; set by hand because it's trouble to get it from other places
+#define FAKE_COLS 4861
+#define FAKE_ROWS 4913
+
+
+
+#define MIN_MAG -16                     // Minimum magnitude to consider
 
 float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
@@ -59,4 +66,31 @@
         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
         pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
+
+#if 1
+        // Excise bright sources that might be saturated
+        psArray *sources = sourceLists->data[i]; // Sources for image
+        for (int j = 0; j < sources->n; j++) {
+            pmSource *source = sources->data[j]; // Source of interest
+            if (source && source->psfMag < MIN_MAG) {
+                psFree(source);
+                sources->data[j] = NULL;
+            }
+        }
+#endif
+
+
+#if 1
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "start_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
 
         float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
@@ -105,7 +139,9 @@
         FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
         psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
+        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
         for (int i = 0; i < matches->n; i++) {
             pmSourceMatch *match = matches->data[i]; // Match of interest
             psVectorInit(mag, NAN);
+            psVectorInit(magErr, NAN);
             for (int j = 0; j < match->num; j++) {
                 if (match->mask->data.PS_TYPE_MASK_DATA[j]) {
@@ -114,11 +150,13 @@
                 int index = match->image->data.U32[j]; // Image index
                 mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
+                magErr->data.F32[index] = match->magErr->data.F32[j];
             }
             for (int j = 0; j < num; j++) {
-                fprintf(outMatches, "%f ", mag->data.F32[j]);
+                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
             }
             fprintf(outMatches, "\n");
         }
         psFree(mag);
+        psFree(magErr);
         fclose(outMatches);
     }
@@ -146,6 +184,31 @@
                 continue;
             }
+
+#if 1
+            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS],
+                y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; // Coordinates of source
+            if (x > 2187 && x < 2191 && y > 2280 && y < 2284) {
+                fprintf(stderr, "Image %d source: (%f,%f) %f --> %f\n",
+                        i, x, y, source->psfMag, source->psfMag + magCorr);
+            }
+#endif
+
             source->psfMag += magCorr;
         }
+
+#if 1
+// Size of fake image; set by hand because it's trouble to get it from other places
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sources, NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "corrected_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
     }
     psFree(trans);
@@ -162,4 +225,7 @@
         psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
                                                transThresh, starRej, starSys);
+        for (int i = 0; i < num; i++) {
+            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+        }
         psFree(trans);
     }