Index: trunk/ppSub/src/ppSubReadout.c
===================================================================
--- trunk/ppSub/src/ppSubReadout.c	(revision 14316)
+++ trunk/ppSub/src/ppSubReadout.c	(revision 14328)
@@ -11,4 +11,5 @@
 
 #define TESTING
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
 
 bool ppSubReadout(pmConfig *config, const pmFPAview *view)
@@ -121,4 +122,60 @@
         }
     }
+    psFree(stamps);
+    stamps = NULL;
+
+#ifdef TESTING
+    {
+        // Generate image with convolution kernels
+        int fullSize = 2 * size + 1 + 1;    // Full size of kernel
+        psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
+        psImageInit(convKernels, NAN);
+        for (int j = -2; j <= 2; j++) {
+            for (int i = -2; i <= 2; i++) {
+                psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
+                                                           (float)j / 2.0); // Image of the kernel
+                if (!kernel) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+                    psFree(convKernels);
+                    goto ERROR;
+                }
+
+                if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
+                                          (j + 2) * fullSize, "=") == 0) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
+                    psFree(kernel);
+                    psFree(convKernels);
+                    goto ERROR;
+                }
+                psFree(kernel);
+            }
+        }
+
+        // XXX What do we do with this image?
+
+        psFits *kernelFile = psFitsOpen("kernel.fits", "w");
+        (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL);
+        psFitsClose(kernelFile);
+
+        psFree(convKernels);
+    }
+
+    {
+        // Generate images of the kernel components
+        psMetadata *header = psMetadataAlloc(); // Header
+        for (int i = 0; i < solution->n; i++) {
+            psString name = NULL;       // Header keyword
+            psStringAppend(&name, "SOLN%04d", i);
+            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
+            psFree(name);
+        }
+        psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
+        psFits *kernelFile = psFitsOpen("kernels.fits", "w");
+        (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
+        psFitsClose(kernelFile);
+        psFree(kernelImages);
+        psFree(header);
+    }
+#endif
 
     psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images
@@ -130,4 +187,7 @@
     psFree(subMask);
     subMask = NULL;
+
+    psFree(kernels);
+    psFree(solution);
 
     // Do the subtraction
@@ -165,41 +225,4 @@
     }
 
-#ifdef TESTING
-    // Generate image with convolution kernels
-    int fullSize = 2 * size + 1 + 1;    // Full size of kernel
-    psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
-    psImageInit(convKernels, NAN);
-    for (int j = -2; j <= 2; j++) {
-        for (int i = -2; i <= 2; i++) {
-            psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
-                                                       (float)j / 2.0); // Image of the kernel
-            if (!kernel) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-                psFree(convKernels);
-                goto ERROR;
-            }
-            if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
-                                      (j + 2) * fullSize, "=") == 0) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
-                psFree(kernel);
-                psFree(convKernels);
-                goto ERROR;
-            }
-        }
-    }
-
-    // XXX What do we do with this image?
-
-    psFits *kernelFile = psFitsOpen("kernel.fits", "w");
-    (void)psFitsWriteImage(kernelFile, NULL, convKernels, 0, NULL);
-    psFitsClose(kernelFile);
-
-    psFree(convKernels);
-#endif
-
-    psFree(kernels);
-    psFree(stamps);
-    psFree(solution);
-
     if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) {
         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
@@ -210,4 +233,29 @@
         pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
     }
+
+    // Copy astrometry over
+    pmFPA *refFPA = refRO->parent->parent->parent; // Reference FPA
+    pmHDU *refHDU = refFPA->hdu;        // Reference HDU
+    pmFPA *outFPA = outCell->parent->parent; // Output FPA
+    pmHDU *outHDU = outFPA->hdu; // Output HDU
+    if (!outHDU || !refHDU) {
+        psWarning("Unable to find HDU at FPA level to copy astrometry.");
+    } else {
+        if (!pmAstromReadWCS(outFPA, outCell->parent, refHDU->header, 1.0)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from reference FPA.");
+            psFree(outRO);
+            return false;
+        }
+        if (!outHDU->header) {
+            outHDU->header = psMetadataAlloc();
+        }
+        if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA.");
+            psFree(outRO);
+            return false;
+        }
+    }
+
+    psFree(outRO);
 
     return true;
@@ -218,4 +266,5 @@
     psFree(stamps);
     psFree(solution);
+    psFree(outRO);
     return false;
 }