Index: /tags/pap_tags/pap_root_080207/ppStack/.cvsignore
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/.cvsignore	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/.cvsignore	(revision 22291)
@@ -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: /tags/pap_tags/pap_root_080207/ppStack/Makefile.am
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/Makefile.am	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/Makefile.am	(revision 22291)
@@ -0,0 +1,3 @@
+SUBDIRS = src
+
+CLEANFILES = *~ core core.*
Index: /tags/pap_tags/pap_root_080207/ppStack/README.txt
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/README.txt	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/README.txt	(revision 22291)
@@ -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: /tags/pap_tags/pap_root_080207/ppStack/autogen.sh
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/autogen.sh	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/autogen.sh	(revision 22291)
@@ -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: /tags/pap_tags/pap_root_080207/ppStack/configure.ac
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/configure.ac	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/configure.ac	(revision 22291)
@@ -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: /tags/pap_tags/pap_root_080207/ppStack/src/.cvsignore
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/.cvsignore	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/.cvsignore	(revision 22291)
@@ -0,0 +1,10 @@
+.deps
+Makefile
+Makefile.in
+.libs
+*.lo
+*.la
+ppStack
+config.h
+config.h.in
+stamp-h1
Index: /tags/pap_tags/pap_root_080207/ppStack/src/Makefile.am
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/Makefile.am	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/Makefile.am	(revision 22291)
@@ -0,0 +1,24 @@
+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		\
+	ppStackReadout.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: /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.c	(revision 22291)
@@ -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: /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.h
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.h	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStack.h	(revision 22291)
@@ -0,0 +1,44 @@
+#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
+    );
+
+// Perform stacking on a readout
+bool ppStackReadout(pmConfig *config,   // Configuration
+                    const pmFPAview *view // View for 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 *output,    ///< Convolved readout
+                  const pmReadout *input, ///< Readout to be convolved
+                  const pmReadout *sourcesRO, ///< Readout with sources
+                  const pmPSF *psf,     ///< Target PSF
+                  const pmConfig *config ///< Configuration
+    );
+
+#endif
Index: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackArguments.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackArguments.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackArguments.c	(revision 22291)
@@ -0,0 +1,183 @@
+#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);
+    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("-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: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackCamera.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackCamera.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackCamera.c	(revision 22291)
@@ -0,0 +1,291 @@
+#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"
+
+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;
+            }
+        }
+
+        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: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackLoop.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackLoop.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackLoop.c	(revision 22291)
@@ -0,0 +1,128 @@
+#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"
+
+bool ppStackLoop(pmConfig *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 *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file
+    if (!input) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n");
+        return false;
+    }
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file
+    if (!output) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n");
+        return false;
+    }
+
+    pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
+    pmHDU *lastHDU = NULL;              // Last HDU that was updated
+
+    // Iterate over the FPA hierarchy
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        return false;
+    }
+
+    pmChip *chip;                       // Chip of interest
+    while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) {
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            return false;
+        }
+
+        pmCell *cell;                // Cell of interest
+        while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                return false;
+            }
+
+            // Put version information into the header
+            pmHDU *hdu = pmHDUFromCell(cell);
+            if (hdu && hdu != lastHDU) {
+                if (!hdu->header) {
+                    hdu->header = psMetadataAlloc();
+                }
+                ppStackVersionMetadata(hdu->header);
+                lastHDU = hdu;
+            }
+
+            pmReadout *readout;         // Readout of interest
+            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) {
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    return false;
+                }
+
+                // Perform the analysis
+                if (!ppStackReadout(config, view)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
+                    return false;
+                }
+
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                    return false;
+                }
+            }
+
+            // Perform statistics on the cell
+            if (stats) {
+                ppStatsFPA(stats, output->fpa, view, maskBlank, config);
+            }
+
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                return false;
+            }
+        }
+
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            return false;
+        }
+    }
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        return false;
+    }
+
+    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: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackMatch.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackMatch.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackMatch.c	(revision 22291)
@@ -0,0 +1,100 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+//#define TESTING
+
+bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO,
+                  const pmPSF *psf, const pmConfig *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, input->image->numCols, input->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
+    if (!pmSubtractionMatch(output, NULL, input, 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);
+        return false;
+    }
+    psFree(fake);
+    psFree(optWidths);
+
+    return true;
+}
Index: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackReadout.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackReadout.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackReadout.c	(revision 22291)
@@ -0,0 +1,416 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+#define ARRAY_BUFFER 16                 // Number to add to array at a time
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+
+//#define NO_CONVOLUTION                  // Don't perform convolutions?
+//#define CONVOLUTION_FILES               // Write convolutions?
+//#define REJECTION_FILES                 // Write rejection mask?
+//#define INSPECTION_FILES                // Write inspection mask?
+
+bool ppStackReadout(pmConfig *config, const pmFPAview *view)
+{
+    // Get the recipe values
+    bool mdok;                          // Status of MD lookup
+    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
+    bool photometry = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // Perform photometry?
+    int psfInstances = psMetadataLookupS32(&mdok, 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(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF
+
+    // Get the output target
+    pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+    pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout
+    pmFPA *outFPA = outCell->parent->parent; // Output FPA
+
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    assert(num > 0);
+
+    // Get the input sources
+    psArray *stack = psArrayAllocEmpty(ARRAY_BUFFER); // The stack of inputs
+    pmReadout *sources = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Sources for stamps
+    psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
+                                                           "^PPSTACK.INPUT$"); // Iterator over input files
+    psMetadataItem *fileItem;           // Item from iteration
+    int fileNum = 0;                    // Number of file
+    psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope
+    pmPSF *outPSF = NULL;               // Ouptut PSF
+    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
+        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+        pmCell *cell = ro->parent;      // The cell
+        pmChip *chip = cell->parent;    // The chip: holds the PSF
+
+        bool mdok;                      // Status of MD lookup
+        pmPSF *psf = psMetadataLookupPtr(&mdok, chip->analysis, "PSPHOT.PSF");
+        if (mdok && psf) {
+            psListAdd(psfList, PS_LIST_TAIL, psf);
+            if (numCols == 0 && numRows == 0) {
+                numCols = ro->image->numCols;
+                numRows = ro->image->numRows;
+            }
+        }
+    }
+
+    bool convolve = false;              // Convolve the input images?
+    if (psfList->n > 0) {
+        psArray *psfArray = psListToArray(psfList); // Array of PSFs
+        outPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel,
+                               psfOrder, psfOrder);
+        psFree(psfArray);
+        if (!outPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            // XXX Free stuff
+            return false;
+        }
+        convolve = true;
+        PS_ASSERT_PTR_NON_NULL(sources, false);
+    }
+
+    // Iterate through again to get the convolved images (or not)
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics
+    psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator
+    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
+
+    psMetadataIteratorSet(fileIter, PS_LIST_HEAD);
+    fileNum = 0;
+    while ((fileItem = psMetadataGetAndIncrement(fileIter))) {
+        assert(fileItem->type == PS_DATA_UNKNOWN);
+        pmFPAfile *inputFile = fileItem->data.V; // An input file
+        pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout
+
+        bool mdok;                      // Status of MD lookup
+        float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis,
+                                              "PPSTACK.WEIGHTING"); // Relative weighting
+        if (!mdok || !isfinite(weighting)) {
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", fileNum);
+            weighting = 1.0;
+        }
+        float scale = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
+        if (!mdok || !isfinite(scale)) {
+            psWarning("No SCALE supplied for image %d --- set to unity.", fileNum);
+            scale = 1.0;
+        }
+
+        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
+#if 0
+        if (!isfinite(exposure)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "CELL.EXPOSURE is not set for input file %ld", stack->n);
+            psFree(stats);
+            psFree(rng);
+            psFree(fileIter);
+            psFree(fpaList);
+            psFree(cellList);
+            psFree(outRO);
+            psFree(stack);
+            return false;
+        }
+#endif
+        totExposure += exposure;        // Total exposure time
+
+        // Generate convolved version of input
+        pmReadout *convolved;
+        if (convolve) {
+            convolved = pmReadoutAlloc(NULL); // Convolved version of input readout
+            // Background subtraction, scaling and normalisation is performed automatically by the image
+            // matching
+            if (!ppStackMatch(convolved, ro, sources, outPSF, config)) {
+                psWarning("Unable to match image %d --- ignoring.", fileNum);
+                psErrorClear();
+                psFree(convolved);
+                // XXX Free the bad image so it's not taking up good memory!
+                continue;
+            }
+
+#ifdef CONVOLUTION_FILES
+            if (convolved->image) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_image.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->image, 0, NULL);
+                psFitsClose(fits);
+            }
+            if (convolved->mask) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_mask.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL);
+                psFitsClose(fits);
+            }
+            if (convolved->weight) {
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "convolved%03d_weight.fits", fileNum);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL);
+                psFitsClose(fits);
+            }
+#endif // CONVOLUTION_FILES
+
+        } else {
+            // No convolution --- just use the unconvolved images as the "convolved" images, with some tweaks
+            convolved = psMemIncrRefCounter(ro);
+            sources = NULL;
+
+            // Brain-dead background subtraction
+            if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskBad, rng)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n);
+                psFree(stats);
+                psFree(rng);
+                psFree(fileIter);
+                psFree(fpaList);
+                psFree(cellList);
+                psFree(stack);
+                psFree(outRO);
+                psFree(convolved);
+                return false;
+            }
+            psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian);
+            (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32));
+
+            // Apply scaling
+            (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, PS_TYPE_F32));
+
+            // Normalise each input by the exposure time
+            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(stats);
+                psFree(rng);
+                psFree(fileIter);
+                psFree(fpaList);
+                psFree(cellList);
+                psFree(outRO);
+                psFree(convolved);
+                psFree(stack);
+                return false;
+            }
+
+            (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32));
+            if (ro->weight) {
+                (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32));
+            }
+        }
+
+        if (fileNum == 0) {
+            // Copy astrometry over
+            pmFPA *fpa = ro->parent->parent->parent; // Template FPA
+            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, outCell->parent, 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, outCell->parent, WCS_TOLERANCE)) {
+                        psErrorClear();
+                        psWarning("Unable to write WCS astrometry to output FPA.");
+                    }
+                }
+            }
+        }
+
+        // Don't need the original any more!
+        psFree(ro);
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!convolved->mask) {
+            convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows,
+                                           PS_TYPE_MASK);
+            psImageInit(convolved->mask, 0);
+        }
+
+        psListAdd(fpaList, PS_LIST_TAIL, inputFile->fpa);
+        psListAdd(cellList, PS_LIST_TAIL, ro->parent);
+
+        pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack
+        psFree(convolved);
+        psArrayAdd(stack, ARRAY_BUFFER, data);
+        psFree(data);                   // Drop reference
+
+        fileNum++;
+    }
+    psFree(fileIter);
+    psFree(stats);
+    psFree(rng);
+
+    if (fileNum == 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not enough good files to combine.");
+        psFree(fpaList);
+        psFree(cellList);
+        psFree(stack);
+        psFree(outRO);
+        return false;
+    }
+
+    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);
+        psFree(outRO);
+        return false;
+    }
+
+#ifdef INSPECTION_FILES
+    {
+        psFits *fits = psFitsOpen("combined.fits", "w");
+        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,
+                                                        0, data->readout->image->numRows),
+                                            maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "inspect%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, inspected, 0, NULL);
+        psFree(inspected);
+        psFitsClose(fits);
+    }
+#endif
+
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        pmReadout *readout = data->readout; // Readout for this image
+
+        // Extract the regions and solutions used in the image matching
+        psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
+        {
+            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
+                                                               "^SUBTRACTION.REGION$"); // Iterator
+            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);
+        }
+        psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions
+        {
+            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
+                                                               "^SUBTRACTION.SOLUTION$"); // Iterator
+            psMetadataItem *item = NULL;// Item from iteration
+            while ((item = psMetadataGetAndIncrement(iter))) {
+                assert(item->type == PS_DATA_VECTOR);
+                solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V);
+            }
+            psFree(iter);
+        }
+        assert(regions->n == solutions->n);
+
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
+                                                            "SUBTRACTION.KERNEL"); // Kernels
+
+        psPixels *reject = pmStackReject(data->pixels, threshold, regions,
+                                         solutions, kernels); // List of pixels to reject
+        psFree(data->pixels);
+        data->pixels = reject;
+
+        psFree(solutions);
+        psFree(regions);
+    }
+
+#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,
+                                                       0, data->readout->image->numRows),
+                                           maskBlank);
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "reject%03d.fits", 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);
+        psFree(outRO);
+        return false;
+    }
+
+    if (!convolve) {
+        // Restore image to counts using the total exposure time
+        (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
+        if (outRO->weight) {
+            (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32));
+        }
+    }
+
+    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);
+
+    if (photometry) {
+        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
+        pmFPACopy(photFile->fpa, outRO->parent->parent->parent);
+
+        if (!psphotReadout(config, view)) {
+            psWarning("Unable to perform photometry on stacked image.");
+            psErrorStackPrint(stderr, "Error stack from photometry:");
+            psErrorClear();
+        }
+
+        pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
+    }
+
+    psFree(stack);
+    psFree(outRO);                      // Drop reference
+
+    return true;
+}
Index: /tags/pap_tags/pap_root_080207/ppStack/src/ppStackVersion.c
===================================================================
--- /tags/pap_tags/pap_root_080207/ppStack/src/ppStackVersion.c	(revision 22291)
+++ /tags/pap_tags/pap_root_080207/ppStack/src/ppStackVersion.c	(revision 22291)
@@ -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;
+}