Index: /branches/pap_branches/pap_branch_080207/ppStack/.cvsignore
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/.cvsignore	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/.cvsignore	(revision 22229)
@@ -0,0 +1,17 @@
+aclocal.m4
+compile
+config.guess
+config.log
+config.status
+config.sub
+configure
+depcomp
+install-sh
+libtool
+ltmain.sh
+Makefile
+Makefile.in
+missing
+ppStac.pc
+test
+autom4te.cache
Index: /branches/pap_branches/pap_branch_080207/ppStack/Makefile.am
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/Makefile.am	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/Makefile.am	(revision 22229)
@@ -0,0 +1,3 @@
+SUBDIRS = src
+
+CLEANFILES = *~ core core.*
Index: /branches/pap_branches/pap_branch_080207/ppStack/README.txt
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/README.txt	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/README.txt	(revision 22229)
@@ -0,0 +1,25 @@
+ppStack needs to be optimised, principally to reduce the memory
+footprint (10 or so images will get it to over 4 GB).  The proposed
+strategy follows, more or less, the same as for ppMerge.
+
+* Convolution steps:
+  - For each file:
+    + read entire image
+    + convolve to reference PSF
+    + write convolved image
+
+* Combination step:
+  - For each chunk:
+    + For each file, read chunk
+    + combine chunk
+  - Probably have to repeat, to do the multi-pass rejection
+
+
+To do:
+* Need file PPSTACK.OUTPUT.CONV (image, mask, weight versions)
+* Need file PPSTACK.INPUT.CONV (image, mask, weight versions)
+* Check: does pmFPAfileChecks do piece-by-piece read?
+* Can we reset a file, so that we can read it twice through?
+* pmStack functions probably need work so that they operate on a piece
+  of the image at a time.  Will there be issues with the masks in the
+  second pass?
Index: /branches/pap_branches/pap_branch_080207/ppStack/autogen.sh
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/autogen.sh	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/autogen.sh	(revision 22229)
@@ -0,0 +1,103 @@
+#!/bin/sh
+# Run this to generate all the initial makefiles, etc.
+
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+ORIGDIR=`pwd`
+cd $srcdir
+
+PROJECT=ppStack
+TEST_TYPE=-f
+# change this to be a unique filename in the top level dir
+FILE=autogen.sh
+
+DIE=0
+
+LIBTOOLIZE=libtoolize
+ACLOCAL="aclocal $ACLOCAL_FLAGS"
+AUTOHEADER=autoheader
+AUTOMAKE=automake
+AUTOCONF=autoconf
+
+($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $LIBTOOlIZE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/"
+        DIE=1
+}
+
+($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $ACLOCAL installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOHEADER installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOMAKE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOCONF installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+test $TEST_TYPE $FILE || {
+        echo "You must run this script in the top-level $PROJECT directory"
+        exit 1
+}
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed"
+$ACLOCAL || echo "$ACLOCAL failed"
+$AUTOHEADER || echo "$AUTOHEADER failed"
+$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed"
+$AUTOCONF || echo "$AUTOCONF failed"
+
+cd $ORIGDIR
+
+run_configure=true
+for arg in $*; do
+    case $arg in
+        --no-configure)
+            run_configure=false
+            ;;
+        *)
+            ;;
+    esac
+done
+
+if $run_configure; then
+    $srcdir/configure --enable-maintainer-mode "$@"
+    echo
+    echo "Now type 'make' to compile $PROJECT."
+else
+    echo
+    echo "Now run 'configure' and 'make' to compile $PROJECT."
+fi
Index: /branches/pap_branches/pap_branch_080207/ppStack/configure.ac
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/configure.ac	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/configure.ac	(revision 22229)
@@ -0,0 +1,34 @@
+AC_PREREQ(2.61)
+
+AC_INIT([ppStack], [0.1.1], [ipp-support@ifa.hawaii.edu])
+AC_CONFIG_SRCDIR([src])
+
+AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2])
+AM_CONFIG_HEADER([src/config.h])
+AM_MAINTAINER_MODE
+
+IPP_STDCFLAGS
+
+AC_LANG(C)
+AC_GNU_SOURCE
+AC_PROG_CC_C99
+AC_PROG_INSTALL
+AC_PROG_LIBTOOL
+AC_SYS_LARGEFILE 
+
+PKG_CHECK_MODULES([PSLIB],    [pslib >= 1.0.0])
+PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0])
+PKG_CHECK_MODULES([PSPHOT],   [psphot >= 0.8.0]) 
+PKG_CHECK_MODULES([PPSTATS],  [ppStats >= 1.0.0]) 
+
+IPP_STDOPTS
+CFLAGS="${CFLAGS=} -Wall -Werror"
+
+AC_SUBST([PPSTACK_CFLAGS])
+AC_SUBST([PPSTACK_LIBS])
+
+AC_CONFIG_FILES([
+  Makefile
+  src/Makefile
+])
+AC_OUTPUT
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/.cvsignore
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/.cvsignore	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/.cvsignore	(revision 22229)
@@ -0,0 +1,10 @@
+.deps
+Makefile
+Makefile.in
+.libs
+*.lo
+*.la
+ppStack
+config.h
+config.h.in
+stamp-h1
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/Makefile.am
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/Makefile.am	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/Makefile.am	(revision 22229)
@@ -0,0 +1,26 @@
+bin_PROGRAMS = ppStack
+
+ppStack_CFLAGS 	= $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS)
+ppStack_LDFLAGS = $(PSLIB_LIBS)   $(PSMODULE_LIBS)   $(PSPHOT_LIBS)   $(PPSTATS_LIBS)   $(PPSTACK_LIBS)
+
+ppStack_SOURCES =		\
+	ppStack.c		\
+	ppStackArguments.c	\
+	ppStackCamera.c		\
+	ppStackLoop.c		\
+	ppStackPSF.c		\
+	ppStackReadout.c	\
+	ppStackPhotometry.c	\
+	ppStackVersion.c	\
+	ppStackMatch.c
+
+noinst_HEADERS = 		\
+	ppStack.h
+
+CLEANFILES = *~
+
+clean-local:
+	-rm -f TAGS
+# Tags for emacs
+tags:
+	etags `find . -name \*.[ch] -print`
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.c	(revision 22229)
@@ -0,0 +1,63 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+#define TIMER_NAME "PPSTACK"            // Name of timer
+
+int main(int argc, char *argv[])
+{
+    psExit exitValue = PS_EXIT_SUCCESS; // Exit value
+    psTimerStart(TIMER_NAME);
+    psLibInit(NULL);
+
+    pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration
+    if (!config) {
+        psErrorStackPrint(stderr, "Error reading configuration.");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppStackArguments(argc, argv, config)) {
+        psErrorStackPrint(stderr, "Error reading arguments.\n");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!pmModelClassInit()) {
+        psErrorStackPrint(stderr, "Error initialising model classes.\n");
+        exitValue = PS_EXIT_PROG_ERROR;
+        goto die;
+    }
+
+    if (!ppStackCamera(config)) {
+        psErrorStackPrint(stderr, "Error reading configuration.\n");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppStackLoop(config)) {
+        psErrorStackPrint(stderr, "Error performing combination.\n");
+        exitValue = PS_EXIT_DATA_ERROR;
+        goto die;
+    }
+
+
+     // Common code for the death.
+die:
+    psTrace("ppStack", 1, "Finished at %f sec\n", psTimerMark(TIMER_NAME));
+    psTimerStop();
+
+    psFree(config);
+    pmModelClassCleanup();
+    pmConfigDone();
+    psLibFinalize();
+
+    exit(exitValue);
+}
+
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.h
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.h	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStack.h	(revision 22229)
@@ -0,0 +1,60 @@
+#ifndef PP_STACK_H
+#define PP_STACK_H
+
+#define PPSTACK_RECIPE "PPSTACK"
+
+#include <psmodules.h>
+
+// Parse command-line arguments
+bool ppStackArguments(int argc, char *argv[], // Command-line arguments
+                      pmConfig *config  // Configuration
+    );
+
+// Parse cameras
+bool ppStackCamera(pmConfig *config     // Configuration
+    );
+
+// Loop over the inputs, doing the combination
+bool ppStackLoop(pmConfig *config       // Configuration
+    );
+
+// Determine target PSF for input images
+pmPSF *ppStackPSF(const pmConfig *config, // Configuration
+                  int numCols, int numRows, // Size of image
+                  const psList *list    // List of input PSFs
+    );
+
+// Perform stacking on a readout
+bool ppStackReadout(const pmConfig *config,   // Configuration
+                    pmReadout *outRO,   // Output readout
+                    const psArray *readouts, // Input readouts
+                    const psArray *regions, // Array with array of regions used in each PSF matching
+                    const psArray *kernels // Array with array of kernels used in each PSF matching
+    );
+
+// Perform photometry on stack
+bool ppStackPhotometry(pmConfig *config, // Configuration
+                       const pmReadout *readout, // Readout to be photometered
+                       const pmFPAview *view // View to readout
+    );
+
+// Return software version
+psString ppStackVersion(void);
+
+// Return long description of software version
+psString ppStackVersionLong(void);
+
+// Supplement metadata with software version
+void ppStackVersionMetadata(psMetadata *metadata // Metadata to supplement
+    );
+
+/// Convolve image to match specified seeing
+bool ppStackMatch(pmReadout *readout, ///< Readout to be convolved; replaced with output
+                  psArray **regions, // Array of regions used in each PSF matching, returned
+                  psArray **kernels, // Array of kernels used in each PSF matching, returned
+                  const pmReadout *sourcesRO, ///< Readout with sources
+                  const pmPSF *psf,     ///< Target PSF
+                  const pmConfig *config ///< Configuration
+    );
+
+#endif
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackArguments.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackArguments.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackArguments.c	(revision 22229)
@@ -0,0 +1,185 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+// XXX add in the version info as in ppImage
+
+// Print usage information and die
+static void usage(const char *program,  // Name of the program
+                  psMetadata *arguments, // Command-line arguments
+                  pmConfig *config      // Configuration
+    )
+{
+    fprintf(stderr, "\nPan-STARRS Image combination\n\n");
+    fprintf(stderr,
+            "Usage: %s INPUTS.mdc OUTPUT_ROOT [-sources STAMPS.cmf | -stamps STAMPS.dat]\n"
+            "where INPUTS.mdc contains various METADATAs, each with:\n"
+            "\tIMAGE(STR):     Image filename\n"
+            "\tMASK(STR):      Mask filename\n"
+            "\tWEIGHT(STR)     Weight map filename\n"
+            "\tPSF(STR)        PSF filename\n"
+            "\tWEIGHTING(F32): Relative weighting to be applied\n"
+            "\tSCALE(F32):     Relative scaling to be applied\n",
+            program);
+    fprintf(stderr, "\n");
+    psArgumentHelp(arguments);
+    psFree(arguments);
+    psFree(config);
+    pmConfigDone();
+    psLibFinalize();
+    exit(PS_EXIT_CONFIG_ERROR);
+}
+
+// Get a float-point value from the command-line or recipe, and add it to the arguments
+#define VALUE_ARG_RECIPE_FLOAT(ARGNAME, RECIPENAME, TYPE) { \
+    ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \
+    if (isnan(value)) { \
+        bool mdok; \
+        value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
+        if (!mdok) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+                RECIPENAME, PPSTACK_RECIPE); \
+            goto ERROR; \
+        } \
+    } \
+    psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
+}
+
+// Get an integer value from the command-line or recipe, and add it to the arguments
+#define VALUE_ARG_RECIPE_INT(ARGNAME, RECIPENAME, TYPE, UNSET) { \
+    ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \
+    if (value == UNSET) { \
+        bool mdok; \
+        value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \
+        if (!mdok) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+                RECIPENAME, PPSTACK_RECIPE); \
+            goto ERROR; \
+        } \
+    } \
+    psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
+}
+
+// Get a mask value from the command-line or recipe, and add it to the arguments
+#define VALUE_ARG_RECIPE_MASK(ARGNAME, RECIPENAME) { \
+    bool mdok; \
+    const char *name = psMetadataLookupStr(&mdok, arguments, ARGNAME); \
+    if (!mdok || !name || strlen(name) == 0) { \
+        name = psMetadataLookupStr(NULL, recipe, RECIPENAME); \
+        if (!name) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \
+                RECIPENAME, PPSTACK_RECIPE); \
+            goto ERROR; \
+        } \
+    } \
+    psMaskType value = pmConfigMask(name, config); \
+    psMetadataAddU8(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \
+}
+
+// Get a string value from the command-line and add it to the target
+static bool valueArgStr(pmConfig *config,      // Configuration
+                        psMetadata *arguments, // Command-line arguments
+                        const char *argName, // Argument name in the command-line arguments
+                        const char *mdName, // Name for value in the metadata
+                        psMetadata *target // Target metadata to which to add value
+                        )
+{
+    psString value = psMetadataLookupStr(NULL, arguments, argName); // Value of interest
+    if (value && strlen(value) > 0) {
+        return psMetadataAddStr(target, PS_LIST_TAIL, mdName, 0, NULL, value);
+    }
+    return false;
+}
+
+bool ppStackArguments(int argc, char *argv[], pmConfig *config)
+{
+    assert(config);
+    // bool mdok;                          // Status of MD lookup
+
+    pmConfigFileSetsMD(config->arguments, &argc, argv, "PPSTACK.SOURCES", "-sources", NULL);
+
+    psMetadata *arguments = psMetadataAlloc(); // Command-line arguments
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0, "Combination rejection thresold (sigma)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-convolve-rej", 0, "Convolution rejection thresold (sigma)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-extent", 0, "Extent of convolution (sigma)", NAN);
+    psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0);
+    psMetadataAddU8(arguments,  PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-rows", 0, "Rows to read at once", 128);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Do photometry on stacked image?", false);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-instances", 0, "Number of instances for PSF generation", 5);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-psf-radius", 0, "Radius for PSF generation", 20.0);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-psf-model", 0, "Model name for PSF generation", "PS_MODEL_RGAUSS");
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-order", 0, "Spatial order for PSF generation", 3);
+
+    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) {
+        usage(argv[0], arguments, config);
+    }
+
+    const char *stampsName = psMetadataLookupStr(NULL, arguments, "-stamps"); // Name of stamps file
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STAMPS", 0, "Stamps file", stampsName);
+    // if (!stampsName && !psMetadataLookupPtr(&mdok, config->arguments, "PPSTACK.SOURCES")) {
+    // psError(PS_ERR_BAD_PARAMETER_VALUE, true, "One of -sources or -stamps must be specified.");
+    //goto ERROR;
+    // }
+
+    unsigned int numBad = 0;                     // Number of bad lines
+    psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[1], false); // Information about inputs
+    if (!inputs || numBad > 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs.");
+        goto ERROR;
+    }
+    psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0,
+                          "Metadata with input details", inputs);
+    psFree(inputs);
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0,
+                     "Root name of the output image list", argv[2]);
+
+    valueArgStr(config, arguments, "-stats", "STATS", config->arguments);
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppSim
+    if (!recipe) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
+        goto ERROR;
+    }
+
+    VALUE_ARG_RECIPE_INT("-iter",             "ITER",           S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-combine-rej",    "COMBINE.REJ",    F32);
+    VALUE_ARG_RECIPE_FLOAT("-convolve-rej",   "CONVOLVE.REJ",   F32);
+    VALUE_ARG_RECIPE_FLOAT("-extent",         "EXTENT",         F32);
+    VALUE_ARG_RECIPE_MASK("-mask-bad",        "MASK.BAD");
+    VALUE_ARG_RECIPE_MASK("-mask-blank",      "MASK.BLANK");
+    VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
+    VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
+
+    VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-psf-radius",  "PSF.RADIUS",    F32);
+    VALUE_ARG_RECIPE_INT("-psf-order",     "PSF.ORDER",     S32, 0);
+    valueArgStr(config, arguments, "-psf-model", "PSF.MODEL", config->arguments);
+
+    if (psMetadataLookupBool(NULL, arguments, "-photometry") ||
+        psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
+        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PHOTOMETRY", 0,
+                          "Do photometry on stacked image?", true);
+    }
+
+    psTrace("ppStack", 1, "Done reading command-line arguments\n");
+    psFree(arguments);
+    return true;
+
+ERROR:
+    psFree(arguments);
+    return false;
+}
+
+
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackCamera.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackCamera.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackCamera.c	(revision 22229)
@@ -0,0 +1,374 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+
+#if 0
+// Define an output convolved image file
+static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name
+                                        pmFPA *fpa, // FPA to bind
+                                        const pmConfig *config, // Configuration
+                                        pmFPAfileType type // Expected type
+    )
+{
+    pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, name);
+    if (!file) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to define output convolved file %s", name);
+        return NULL;
+    }
+    if (file->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTCONV is not of type %s", pmFPAfileStringFromType(type));
+        return NULL;
+    }
+
+    return file;
+}
+
+// Define an input convolved image file, using the output as a basis
+static pmFPAfile *defineInputConvolved(const char *inputName, // Input FPA file name
+                                       pmFPAfile *outFile, // Corresponding output FPA file
+                                       pmConfig *config, // Configuration
+                                       pmFPAfileType type // Expected type
+    )
+{
+    pmFPAview *view = pmFPAviewAlloc(0); // View into sky cells
+    view->chip = view->cell = view->readout = 0;
+
+    psString imageName = pmFPANameFromRule(outFile->filerule, outFile->fpa, view);
+    psArray *imageNames = psArrayAlloc(1);
+    imageNames->data[0] = imageName;
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INCONV.FILENAMES", PS_META_REPLACE,
+                       "Filenames of input convolved image files", imageNames);
+    psFree(imageNames);
+    bool found = false;                 // Found the file?
+    pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INCONV",
+                                                   "INCONV.FILENAMES");
+    psMetadataRemoveKey(config->arguments, "INCONV.FILENAMES");
+    if (!imageFile || !found) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to define %s file", inputName);
+        return NULL;
+    }
+    if (imageFile->type != type) {
+        psError(PS_ERR_IO, true, "PPSTACK.INCONV is not of type %s",
+                pmFPAfileStringFromType(type));
+        return NULL;
+    }
+
+    return imageFile;
+}
+#endif
+
+
+
+bool ppStackCamera(pmConfig *config)
+{
+    bool haveWeights = false;           // Do we have weight maps?
+    bool havePSFs = false;              // Do we have PSFs?
+
+    psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info
+    psMetadataIterator *iter = psMetadataIteratorAlloc(inputs, PS_LIST_HEAD, NULL); // Iterator
+    psMetadataItem *item;               // Item from iteration
+    int i = 0;                          // Counter
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        if (item->type != PS_DATA_METADATA) {
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                    "Component %s of the input metadata is not of type METADATA", item->name);
+            psFree(iter);
+            return false;
+        }
+
+        psMetadata *input = item->data.md; // The input metadata of interest
+
+        psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image
+        if (!image || strlen(image) == 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name);
+            psFree(iter);
+            return false;
+        }
+
+        bool mdok;
+        psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask
+        psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map
+        psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF
+
+        float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting
+        if (!mdok || isnan(weighting) || weighting == 0.0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks WEIGHTING of type F32", item->name);
+            psFree(iter);
+            return false;
+        }
+
+        float scale = psMetadataLookupF32(&mdok, input, "SCALE"); // Relative scale
+        if (!mdok || isnan(scale) || scale == 0.0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks SCALE of type F32", item->name);
+            psFree(iter);
+            return false;
+        }
+
+        // Add the image file
+        psArray *imageFiles = psArrayAlloc(1); // Array of filenames for this FPA
+        imageFiles->data[0] = psMemIncrRefCounter(image);
+        psMetadataAddArray(config->arguments, PS_LIST_TAIL, "IMAGE.FILENAMES", PS_META_REPLACE,
+                           "Filenames of image files", imageFiles);
+        psFree(imageFiles);
+
+        bool found = false;             // Found the file?
+        pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INPUT", "IMAGE.FILENAMES");
+        if (!imageFile || !found) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
+            return false;
+        }
+        if (imageFile->type != PM_FPA_FILE_IMAGE) {
+            psError(PS_ERR_IO, true, "PPSTACK.INPUT is not of type IMAGE");
+            return false;
+        }
+
+        // Optionally add the mask file
+        if (mask && strlen(mask) > 0) {
+            psArray *maskFiles = psArrayAlloc(1); // Array of filenames for this FPA
+            maskFiles->data[0] = psMemIncrRefCounter(mask);
+            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "MASK.FILENAMES", PS_META_REPLACE,
+                               "Filenames of mask files", maskFiles);
+            psFree(maskFiles);
+
+            bool status;
+            pmFPAfile *maskFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.MASK",
+                                                        "MASK.FILENAMES");
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
+                return false;
+            }
+            if (maskFile->type != PM_FPA_FILE_MASK) {
+                psError(PS_ERR_IO, true, "PPSTACK.INPUT.MASK is not of type MASK");
+                return false;
+            }
+        }
+
+        // Optionally add the weight file
+        if (weight && strlen(weight) > 0) {
+            haveWeights = true;
+            psArray *weightFiles = psArrayAlloc(1); // Array of filenames for this FPA
+            weightFiles->data[0] = psMemIncrRefCounter(weight);
+            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "WEIGHT.FILENAMES", PS_META_REPLACE,
+                               "Filenames of weight files", weightFiles);
+            psFree(weightFiles);
+
+            bool status;
+            pmFPAfile *weightFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.WEIGHT",
+                                                          "WEIGHT.FILENAMES");
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from weight %d (%s)", i, weight);
+                return false;
+            }
+            if (weightFile->type != PM_FPA_FILE_WEIGHT) {
+                psError(PS_ERR_IO, true, "PPSTACK.INPUT.WEIGHT is not of type WEIGHT");
+                return false;
+            }
+        }
+
+        // Optionally add the psf file
+        if (psf && strlen(psf) > 0) {
+            havePSFs = true;
+            psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA
+            psfFiles->data[0] = psMemIncrRefCounter(psf);
+            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE,
+                               "Filenames of PSF files", psfFiles);
+            psFree(psfFiles);
+
+            bool status;
+            pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PSPHOT.PSF.LOAD",
+                                                       "PSF.FILENAMES");
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf);
+                return false;
+            }
+            if (psfFile->type != PM_FPA_FILE_PSF) {
+                psError(PS_ERR_IO, true, "PSPHOT.PSF.LOAD is not of type PSF");
+                return false;
+            }
+        }
+
+#if 0
+        // Output convolved files
+        pmFPAfile *outconvImage  = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
+                                                         PM_FPA_FILE_IMAGE);
+        pmFPAfile *outconvMask   = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config,
+                                                         PM_FPA_FILE_MASK);
+        pmFPAfile *outconvWeight = defineOutputConvolved("PPSTACK.OUTCONV.WEIGHT", imageFile->fpa, config,
+                                                         PM_FPA_FILE_WEIGHT);
+        if (!outconvImage || !outconvMask || !outconvWeight) {
+            return false;
+        }
+
+        // Input convolved files
+        pmFPAfile *inconvImage  = defineInputConvolved("PPSTACK.INCONV", outconvImage, config,
+                                                       PM_FPA_FILE_IMAGE);
+        pmFPAfile *inconvMask   = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config,
+                                                       PM_FPA_FILE_MASK);
+        pmFPAfile *inconvWeight = defineInputConvolved("PPSTACK.INCONV.WEIGHT", outconvWeight, config,
+                                                       PM_FPA_FILE_WEIGHT);
+        if (!inconvImage || !inconvMask || !inconvWeight) {
+            return false;
+        }
+#endif
+
+        psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
+                         "Relative weighting for image", weighting);
+        psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.SCALE", 0,
+                         "Relative scale for image", scale);
+
+        i++;
+    }
+    psFree(iter);
+    psMetadataRemoveKey(config->arguments, "IMAGE.FILENAMES");
+    if (psMetadataLookup(config->arguments, "MASK.FILENAMES")) {
+        psMetadataRemoveKey(config->arguments, "MASK.FILENAMES");
+    }
+    if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) {
+        psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES");
+    }
+    if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) {
+        psMetadataRemoveKey(config->arguments, "PSF.FILENAMES");
+    }
+
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i);
+
+    // Output image
+    pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the output
+    if (!fpa) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
+        return false;
+    }
+    pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSTACK.OUTPUT");
+    psFree(fpa);                        // Drop reference
+    if (!output) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT"));
+        return false;
+    }
+    if (output->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT is not of type IMAGE");
+        return false;
+    }
+    output->save = true;
+
+    // Generate the correct structure
+    pmFPALevel phuLevel = pmFPAPHULevel(output->format); // Level at which PHU goes
+    pmFPAview *view = pmFPAviewAlloc(0);// View for current level
+    if (phuLevel == PM_FPA_LEVEL_FPA) {
+        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+            psFree(fpa);
+            psFree(view);
+            return NULL;
+        }
+    } else {
+        pmChip *chip;                       // Chip from FPA
+        while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
+            if (phuLevel == PM_FPA_LEVEL_CHIP) {
+                if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+                    psFree(fpa);
+                    psFree(view);
+                    return NULL;
+                }
+            } else {
+                pmCell *cell;                   // Cell from chip
+                while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
+                    if (phuLevel == PM_FPA_LEVEL_CELL) {
+                        if (!pmFPAAddSourceFromView(fpa, "Stack", view, output->format)) {
+                            psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
+                            psFree(fpa);
+                            psFree(view);
+                            return NULL;
+                        }
+                    }
+                }
+            }
+        }
+    }
+    psFree(view);
+
+    // Output mask
+    pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.MASK");
+    if (!outMask) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK"));
+        return false;
+    }
+    if (outMask->type != PM_FPA_FILE_MASK) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.MASK is not of type MASK");
+        return false;
+    }
+    outMask->save = true;
+
+    // Output weight
+    if (haveWeights) {
+        pmFPAfile *outWeight = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.WEIGHT");
+        if (!outWeight) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.WEIGHT"));
+            return false;
+        }
+        if (outWeight->type != PM_FPA_FILE_WEIGHT) {
+            psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.WEIGHT is not of type WEIGHT");
+            return false;
+        }
+        outWeight->save = true;
+    }
+
+
+    // Output PSF
+    if (havePSFs) {
+        pmFPAfile *outPSF = pmFPAfileDefineOutput(config, output->fpa, "PSPHOT.PSF.SAVE");
+        if (!outPSF) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.PSF.SAVE"));
+            return false;
+        }
+        if (outPSF->type != PM_FPA_FILE_PSF) {
+            psError(PS_ERR_IO, true, "PSPHOT.PSF.SAVE is not of type PSF");
+            return false;
+        }
+        outPSF->save = true;
+    }
+
+    // Sources for use as stamps
+    bool status = false;                // Found the file?
+    pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSTACK.SOURCES", "PPSTACK.SOURCES");
+    if (!status) {
+        psError(PS_ERR_IO, false, "Failed to load file definition PPSTACK.SOURCES");
+        return false;
+    }
+    if (sources && sources->type != PM_FPA_FILE_CMF) {
+        psError(PS_ERR_IO, true, "PPSTACK.SOURCES is not of type CMF");
+        return false;
+    }
+
+    // For photometry, we operate on the chip-mosaicked image
+    // we create a copy of the mosaicked image for psphot so we can write out a clean image
+    bool mdok = false;
+    bool doPhotom = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // perform photometry
+    if (doPhotom) {
+        // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by
+        // psphotDefineFiles
+        pmFPAfile *psphotInput = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT");
+        if (!psphotInput) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.INPUT"));
+            return false;
+        }
+
+        // Define associated psphot input/output files
+        if (!psphotDefineFiles(config, psphotInput)) {
+            psError(PSPHOT_ERR_CONFIG, false,
+                    "Trouble defining the additional input/output files for psphot");
+            return false;
+        }
+    }
+
+    return true;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackLoop.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackLoop.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackLoop.c	(revision 22229)
@@ -0,0 +1,483 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <ppStats.h>
+
+#include "ppStack.h"
+
+// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
+
+#if 0
+// All files in the system
+static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
+                            "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                            "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES",
+                            "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV",
+                            "PSPHOT.BACKGND", "PSPHOT.BACKSUB", "SOURCE.PLOT.MOMENTS",
+                            "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF",
+                            0 };
+#endif
+
+// Files required in preparation for convolution
+static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
+
+// Files required for the convolution
+static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 };
+
+// Output files for the combination
+static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 };
+
+// Files for photometry
+static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
+                             "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
+                             "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
+                             "PSPHOT.INPUT.CMF", 0 };
+
+//#define CONVOLVED_ALREADY               // Already have the convolution products --- testing
+
+
+// Activate/deactivate a list of files
+static void fileActivation(pmConfig *config, // Configuration
+                           char **files, // Files to turn on/off
+                           bool state   // Activation state
+    )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfileActivate(config->files, state, files[i]);
+    }
+    return;
+}
+
+// Activate/deactivate a single element for a list; return array of files
+static psArray *fileActivationSingle(pmConfig *config, // Configuration
+                                     char **files, // Files to turn on/off
+                                     bool state,   // Activation state
+                                     int num // Number of file in sequence
+                                     )
+{
+    assert(config);
+    assert(files);
+
+    psList *list = psListAlloc(NULL);   // List of files
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
+        psListAdd(list, PS_LIST_TAIL, file);
+    }
+
+    psArray *array = psListToArray(list);
+    psFree(list);
+
+    return array;
+}
+
+#if 0
+// Set the data level for a list of files
+static void fileSetDataLevel(pmConfig *config, // Configuration
+                             char **files, // Files for which to set level
+                             pmFPALevel level // Level to set
+                             )
+{
+    assert(config);
+    assert(files);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest
+        for (int j = 0; j < selected->n; j++) {
+            pmFPAfile *file = selected->data[j];
+            assert(file);
+            file->dataLevel = level;
+        }
+        psFree(selected);
+    }
+    return;
+}
+#endif
+
+// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
+static pmFPAview *filesIterateDown(pmConfig *config // Configuration
+                                  )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->chip = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->cell = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return NULL;
+    }
+    return view;
+}
+
+// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells
+static bool filesIterateUp(pmConfig *config // Configuration
+                           )
+{
+    assert(config);
+
+    pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy
+    view->chip = view->cell = view->readout = 0;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->readout = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->cell = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    view->chip = -1;
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+    return true;
+}
+
+#ifndef CONVOLVED_ALREADY
+// Write an image to a FITS file
+static bool writeImage(const char *name, // Name of image
+                       psMetadata *header, // Header
+                       const psImage *image // Image
+                       )
+{
+    assert(name);
+    assert(image);
+
+    psFits *fits = psFitsOpen(name, "w");
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open FITS file to write image.");
+        return false;
+    }
+    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
+        psError(PS_ERR_IO, false, "Unable to write FITS image.");
+        return false;
+    }
+    psFitsClose(fits);
+    return true;
+}
+#endif
+
+
+bool ppStackLoop(pmConfig *config)
+{
+    assert(config);
+
+    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
+
+    bool mdok;                          // Status of MD lookup
+    const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics
+    psMetadata *stats = NULL;           // Container for statistics
+    FILE *statsFile = NULL;             // File stream for statistics
+    if (statsName && strlen(statsName) > 0) {
+        psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename
+        statsFile = fopen(resolved, "w");
+        if (!statsFile) {
+            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
+            psFree(resolved);
+            return false;
+        } else {
+            stats = psMetadataAlloc();
+        }
+        psFree(resolved);
+    }
+
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
+    if (!output) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!");
+        return false;
+    }
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    int numScans = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of scans to read at once
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int overlap = 2 * psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
+
+    // Preparation iteration: Load the sources, and get a target PSF model
+    pmReadout *sources = NULL;          // Readout with sources to use for PSF matching
+    pmPSF *targetPSF = NULL;            // Target PSF
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, prepareFiles, true);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            return false;
+        }
+
+        // We want to hang on to the 'sources' even when its host FPA is blown away
+        sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"));
+        if (!sources) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find sources.");
+            psFree(view);
+            return false;
+        }
+
+        // Generate list of PSFs
+        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
+                                                               "^PPSTACK.INPUT$"); // Iterator
+        psMetadataItem *fileItem; // Item from iteration
+        psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+        int numCols = 0, numRows = 0;   // Size of image
+        while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
+            assert(fileItem->type == PS_DATA_UNKNOWN);
+            pmFPAfile *inputFile = fileItem->data.V; // An input file
+            pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF
+            pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF
+            if (!psf) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to find PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
+                return false;
+            }
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+
+            pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest
+            pmHDU *hdu = pmHDUFromCell(cell);
+            assert(hdu && hdu->header);
+            int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns
+            int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows
+            if (naxis1 <= 0 || naxis2 <= 0) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF.");
+                psFree(view);
+                psFree(sources);
+                psFree(fileIter);
+                psFree(psfList);
+                return false;
+            }
+            if (numCols == 0 && numRows == 0) {
+                numCols = naxis1;
+                numRows = naxis2;
+            }
+        }
+        psFree(fileIter);
+        psFree(view);
+
+        targetPSF = ppStackPSF(config, numCols, numRows, psfList);
+        psFree(psfList);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sources);
+            return false;
+        }
+
+        filesIterateUp(config);
+    }
+
+    const char *suffix = "conv.fits";   // Suffix for convolved images; ultimately this will be from recipe
+    const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+    assert(outName);
+    psArray *imageNames = psArrayAlloc(num);
+    psArray *maskNames = psArrayAlloc(num);
+    psArray *weightNames = psArrayAlloc(num);
+    for (int i = 0; i < num; i++) {
+        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
+        psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix);
+        psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix);
+        psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix);
+        imageNames->data[i] = imageName;
+        maskNames->data[i] = maskName;
+        weightNames->data[i] = weightName;
+    }
+
+    // Generate convolutions and write them to disk
+    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
+    psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
+    psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking
+    for (int i = 0; i < num; i++) {
+        pmFPAfileActivate(config->files, false, NULL);
+        psArray *files = fileActivationSingle(config, convolveFiles, true, i);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout
+        psFree(view);
+
+#ifndef CONVOLVED_ALREADY
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
+        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
+        if (!ppStackMatch(readout, &regions, &kernels, sources, targetPSF, config)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i);
+            psFree(sources);
+            psFree(targetPSF);
+            return false;
+        }
+
+        subRegions->data[i] = regions;
+        subKernels->data[i] = kernels;
+
+        // Write the temporary convolved files
+        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
+        assert(hdu);
+        writeImage(imageNames->data[i],  hdu->header, readout->image);
+        writeImage(maskNames->data[i],   hdu->header, readout->mask);
+        writeImage(weightNames->data[i], hdu->header, readout->weight);
+#endif
+
+        cells->data[i] = psMemIncrRefCounter(readout->parent);
+        filesIterateUp(config);
+    }
+    psFree(sources);
+    psFree(targetPSF);
+
+    // Stack the convolved files
+    {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, combineFiles, true);
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) {
+            fileActivation(config, photFiles, true);
+        }
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(cells);
+            psFree(subKernels);
+            psFree(subRegions);
+            return false;
+        }
+        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+        pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+
+        psArray *readouts = psArrayAlloc(num); // Readouts for convolved images
+        for (int i = 0; i < num; i++) {
+            readouts->data[i] = pmReadoutAlloc(cells->data[i]); // Readout into which to read
+        }
+        psFree(cells);
+
+        // FITS files
+        psArray *imageFits  = psArrayAlloc(num);
+        psArray *maskFits   = psArrayAlloc(num);
+        psArray *weightFits = psArrayAlloc(num);
+        for (int i = 0; i < num; i++) {
+            imageFits->data[i] = psFitsOpen(imageNames->data[i], "r");
+            maskFits->data[i] = psFitsOpen(maskNames->data[i], "r");
+            weightFits->data[i] = psFitsOpen(weightNames->data[i], "r");
+        }
+
+        // Read convolutions by chunks
+        bool more = true;               // More to read?
+        for (int numChunk = 0; more; numChunk++) {
+            for (int i = 0; i < num; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+
+                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) ||
+                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) {
+                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
+                    psFree(readouts);
+                    psFree(subKernels);
+                    psFree(subRegions);
+                    psFree(outRO);
+                    psFree(view);
+                    return false;
+                }
+            }
+
+#ifdef TESTING
+            {
+                pmReadout *ro = readouts->data[0];
+                psTrace("ppStack", 1, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
+                        ro->row0, ro->row0 + ro->image->numRows);
+            }
+#endif
+
+            if (!ppStackReadout(config, outRO, readouts, subRegions, subKernels)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                psFree(readouts);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(outRO);
+                psFree(view);
+                return false;
+            }
+
+            for (int i = 0; i < num && more; i++) {
+                pmReadout *readout = readouts->data[i];
+                assert(readout);
+                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans);
+                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans);
+                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans);
+            }
+        }
+
+        psFree(readouts);
+        psFree(subKernels);
+        psFree(subRegions);
+        for (int i = 0; i < num; i++) {
+            psFitsClose(imageFits->data[i]);
+            psFitsClose(maskFits->data[i]);
+            psFitsClose(weightFits->data[i]);
+        }
+
+        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") &&
+            !ppStackPhotometry(config, outRO, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
+            psFree(outRO);
+            psFree(view);
+            return false;
+        }
+
+        // Statistics on output
+        if (stats) {
+            ppStatsFPA(stats, outCell->parent->parent, view, maskBlank, config);
+        }
+
+        // Put version information into the header
+        pmHDU *hdu = pmHDUFromCell(outRO->parent);
+        if (!hdu) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
+            return false;
+        }
+        if (!hdu->header) {
+            hdu->header = psMetadataAlloc();
+        }
+        ppStackVersionMetadata(hdu->header);
+
+        // Write out the output files
+        fileActivation(config, combineFiles, true);
+        filesIterateUp(config);
+
+        psFree(outRO);
+        psFree(view);
+    }
+
+
+    // Write out summary statistics
+    if (stats) {
+        const char *statsMDC = psMetadataConfigFormat(stats);
+        if (!statsMDC || strlen(statsMDC) == 0) {
+            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
+        } else {
+            fprintf(statsFile, "%s", statsMDC);
+        }
+        psFree((void *)statsMDC);
+        fclose(statsFile);
+
+        psFree(stats);
+    }
+
+    return true;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackMatch.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackMatch.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackMatch.c	(revision 22229)
@@ -0,0 +1,159 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+#define ARRAY_BUFFER 16                 // Number to add to array at a time
+
+
+//#define TESTING
+
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
+                  const pmReadout *sourcesRO, const pmPSF *psf, const pmConfig *config)
+{
+    assert(readout);
+    assert(regions && !*regions);
+    assert(kernels && !*kernels);
+    assert(sourcesRO);
+    assert(psf);
+    assert(config);
+
+    // Look up appropriate values from the ppSub recipe
+    bool mdok;                          // Status of MD lookup
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size
+    int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
+    float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
+    float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
+    int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
+    float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
+    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
+    pmSubtractionKernelsType type =
+        pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type
+    psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
+    psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
+    int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
+    int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
+    int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
+    psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BAD"),
+                                      config); // Value to mask
+    psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"),
+                                        config); // Mask for blank reg.
+    float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
+    bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
+    float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
+    float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
+    float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
+    float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
+    int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
+
+    // These values are specified specifically for stacking
+    const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
+
+    psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
+    if (optimum) {
+        optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+    }
+
+    psArray *sources = NULL;            // Sources in image
+    if (sourcesRO) {
+        sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
+    }
+
+    // Generate a fake image to match to
+    float maxMag = -INFINITY;           // Maximum magnitude of sources
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (source->psfMag > maxMag) {
+            maxMag = source->psfMag;
+        }
+    }
+
+    pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
+
+    if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL,
+                                  psf, powf(10.0, -0.4 * maxMag), 0, false)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+        psFree(fake);
+        psFree(optWidths);
+        return false;
+    }
+
+#ifdef TESTING
+    {
+        psFits *fits = psFitsOpen("fake.fits", "w");
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    // Do the image matching
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
+    if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
+                            sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
+                            binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
+                            maskBlank, badFrac, PM_SUBTRACTION_MODE_1)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+        psFree(fake);
+        psFree(optWidths);
+        psFree(output);
+        return false;
+    }
+    psFree(fake);
+    psFree(optWidths);
+
+    // Replace original images with convolved
+    psFree(readout->image);
+    psFree(readout->mask);
+    psFree(readout->weight);
+    readout->image  = psMemIncrRefCounter(output->image);
+    readout->mask   = psMemIncrRefCounter(output->mask);
+    readout->weight = psMemIncrRefCounter(output->weight);
+
+    // Extract the regions and solutions used in the image matching
+    // This stops them from being freed when we iterate back up the FPA
+    *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_REGION);
+            *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
+        }
+        psFree(iter);
+    }
+    *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
+    {
+        psString regex = NULL;          // Regular expression
+        psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
+        psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
+        psFree(regex);
+        psMetadataItem *item = NULL;// Item from iteration
+        while ((item = psMetadataGetAndIncrement(iter))) {
+            assert(item->type == PS_DATA_UNKNOWN);
+            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
+            // images by scans.
+            pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
+            kernel->numCols = readout->image->numCols;
+            kernel->numRows = readout->image->numRows;
+
+            *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
+        }
+        psFree(iter);
+    }
+    assert((*regions)->n == (*kernels)->n);
+
+
+    psFree(output);
+
+    return true;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPSF.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPSF.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPSF.c	(revision 22229)
@@ -0,0 +1,30 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+pmPSF *ppStackPSF(const pmConfig *config, int numCols, int numRows, const psList *list)
+{
+    // Get the recipe values
+    int psfInstances = psMetadataLookupS32(NULL, config->arguments, "PSF.INSTANCES"); // Number of instances for PSF
+    float psfRadius = psMetadataLookupF32(NULL, config->arguments, "PSF.RADIUS"); // Radius for PSF
+    const char *psfModel = psMetadataLookupStr(NULL, config->arguments, "PSF.MODEL"); // Model for PSF
+    int psfOrder = psMetadataLookupS32(NULL, config->arguments, "PSF.ORDER"); // Spatial order for PSF
+
+    // Solve for the target PSF
+    psArray *array = psListToArray(list); // Array of PSFs
+    pmPSF *psf = pmPSFEnvelope(numCols, numRows, array, psfInstances, psfRadius, psfModel,
+                               psfOrder, psfOrder);
+    psFree(array);
+    if (!psf) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+        return NULL;
+    }
+
+    return psf;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c	(revision 22229)
@@ -0,0 +1,27 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)
+{
+    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
+    pmFPACopy(photFile->fpa, readout->parent->parent->parent);
+
+    if (!psphotReadout(config, view)) {
+        // Clear the error, so that the output files are written.
+        psWarning("Unable to perform photometry on stacked image.");
+        psErrorStackPrint(stderr, "Error stack from photometry:");
+        psErrorClear();
+    }
+
+    pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
+
+    return true;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackReadout.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackReadout.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackReadout.c	(revision 22229)
@@ -0,0 +1,205 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
+//#define REJECTION_FILES                 // Write rejection mask?
+//#define INSPECTION_FILES                // Write inspection mask?
+
+static int sectionNum = 0;              // Section number; for debugging outputs
+
+
+bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                    const psArray *regions, const psArray *kernels)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(regions);
+    assert(kernels);
+    assert(readouts->n == regions->n);
+    assert(regions->n == kernels->n);
+
+
+    // Get the recipe values
+    int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations
+    float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold
+    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
+    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
+    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    pmCell *outCell = outRO->parent;    // Output cell
+    pmChip *outChip = outCell->parent;  // Output chip
+    pmFPA *outFPA = outChip->parent;    // Output FPA
+
+    float totExposure = 0.0;            // Total exposure time
+    psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
+    psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+        pmFPA *fpa = ro->parent->parent->parent; // Parent FPA
+
+        bool mdok;                      // Status of MD lookup
+        float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight
+        if (!mdok || !isfinite(weighting)) {
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
+            weighting = 1.0;
+        }
+
+        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
+        if (!isfinite(exposure)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "CELL.EXPOSURE is not set for input file %ld", stack->n);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(stack);
+            return false;
+        }
+        totExposure += exposure;        // Total exposure time
+
+#if 0
+        if (i == 0) {
+            // Copy astrometry over
+            pmHDU *hdu = fpa->hdu; // Template HDU
+            pmHDU *outHDU = outFPA->hdu; // Output HDU
+            if (!outHDU || !hdu) {
+                psWarning("Unable to find HDU at FPA level to copy astrometry.");
+            } else {
+                if (!pmAstromReadWCS(outFPA, outChip, hdu->header, 1.0)) {
+                    psErrorClear();
+                    psWarning("Unable to read WCS astrometry from input FPA.");
+                } else {
+                    if (!outHDU->header) {
+                        outHDU->header = psMetadataAlloc();
+                    }
+                    if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) {
+                        psErrorClear();
+                        psWarning("Unable to write WCS astrometry to output FPA.");
+                    }
+                }
+            }
+        }
+#endif
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        psListAdd(fpaList, PS_LIST_TAIL, fpa);
+        psListAdd(cellList, PS_LIST_TAIL, ro->parent);
+
+        stack->data[i] = pmStackDataAlloc(ro, weighting);
+    }
+
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, iter, combineRej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(fpaList);
+        psFree(cellList);
+        psFree(stack);
+        return false;
+    }
+
+#ifdef INSPECTION_FILES
+    {
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
+        psFitsClose(fits);
+    }
+
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        psImage *inspected = psPixelsToMask(NULL, data->pixels,
+                                            psRegionSet(0, data->readout->image->numCols - 1,
+                                                        0, data->readout->image->numRows - 1),
+                                            maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
+        psFree(inspected);
+        psFitsClose(fits);
+    }
+#endif
+
+    // Reject pixels
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        pmReadout *readout = data->readout; // Readout of interest
+        int col0 = readout->col0, row0 = readout->row0; // Offset for readout
+        int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
+
+        psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej
+        psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i],
+                                         kernels->data[i]); // Pixels to reject
+        psFree(valid);
+        psFree(data->pixels);
+        data->pixels = reject;
+    }
+
+#ifdef REJECTION_FILES
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        if (!data) {
+            continue;
+        }
+        psImage *rejected = psPixelsToMask(NULL, data->pixels,
+                                           psRegionSet(0, data->readout->image->numCols - 1,
+                                                       0, data->readout->image->numRows - 1),
+                                           maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, rejected, 0, NULL);
+        psFree(rejected);
+        psFitsClose(fits);
+    }
+#endif
+
+    if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, combineRej)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
+        psFree(fpaList);
+        psFree(cellList);
+        psFree(stack);
+        return false;
+    }
+
+    psMetadataAddF32(outCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
+                     "Summed exposure time (sec)", totExposure);
+    psMetadataAddF32(outCell->parent->parent->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE,
+                     "Summed exposure time (sec)", totExposure);
+
+    outRO->data_exists = true;
+    outCell->data_exists = true;
+    outCell->parent->data_exists = true;
+
+    // Copy other concepts
+    pmConceptsAverageFPAs(outFPA, fpaList);
+    pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
+    psFree(fpaList);
+    psFree(cellList);
+
+    psFree(stack);
+
+    sectionNum++;
+
+    return true;
+}
Index: /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackVersion.c
===================================================================
--- /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackVersion.c	(revision 22229)
+++ /branches/pap_branches/pap_branch_080207/ppStack/src/ppStackVersion.c	(revision 22229)
@@ -0,0 +1,60 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <ppStats.h>
+
+#include "ppStack.h"
+
+static const char *cvsTag = "$Name: not supported by cvs2svn $";// CVS tag name
+
+psString ppStackVersion(void)
+{
+    psString version = NULL;            // Version, to return
+    psStringAppend(&version, "%s-%s",PACKAGE_NAME,PACKAGE_VERSION);
+    return version;
+}
+
+psString ppStackVersionLong(void)
+{
+    psString version = ppStackVersion(); // Version, to return
+    psString tag = psStringStripCVS(cvsTag, "Name"); // CVS tag
+    psStringAppend(&version, " (cvs tag %s) %s, %s", tag, __DATE__, __TIME__);
+    psFree(tag);
+    return version;
+}
+
+
+void ppStackVersionMetadata(psMetadata *metadata)
+{
+    PS_ASSERT_METADATA_NON_NULL(metadata,);
+
+    psString pslib = psLibVersionLong();// psLib version
+    psString psmodules = psModulesVersionLong(); // psModules version
+    psString ppStats = ppStatsVersionLong(); // ppStats version
+    psString ppStack = ppStackVersionLong(); // ppStack version
+
+    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
+    psString timeString = psTimeToISO(time); // The time in an ISO string
+    psFree(time);
+    psString head = NULL;               // Head string
+    psStringAppend(&head, "ppStack processing at %s. Component information:", timeString);
+    psFree(timeString);
+
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, head, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, pslib, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, psmodules, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppStats, "");
+    psMetadataAddStr(metadata, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, ppStack, "");
+
+    psFree(head);
+    psFree(pslib);
+    psFree(psmodules);
+    psFree(ppStats);
+    psFree(ppStack);
+
+    return;
+}