Index: /tags/pap_merge_20090108/ppStack/.cvsignore
===================================================================
--- /tags/pap_merge_20090108/ppStack/.cvsignore	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/.cvsignore	(revision 21200)
@@ -0,0 +1,18 @@
+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
+Doxyfile
Index: /tags/pap_merge_20090108/ppStack/Doxyfile.in
===================================================================
--- /tags/pap_merge_20090108/ppStack/Doxyfile.in	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/Doxyfile.in	(revision 21200)
@@ -0,0 +1,1310 @@
+# Doxyfile 1.5.4
+
+# This file describes the settings to be used by the documentation system
+# doxygen (www.doxygen.org) for a project
+#
+# All text after a hash (#) is considered a comment and will be ignored
+# The format is:
+#       TAG = value [value, ...]
+# For lists items can also be appended using:
+#       TAG += value [value, ...]
+# Values that contain spaces should be placed between quotes (" ")
+
+#---------------------------------------------------------------------------
+# Project related configuration options
+#---------------------------------------------------------------------------
+
+# This tag specifies the encoding used for all characters in the config file that 
+# follow. The default is UTF-8 which is also the encoding used for all text before 
+# the first occurrence of this tag. Doxygen uses libiconv (or the iconv built into 
+# libc) for the transcoding. See http://www.gnu.org/software/libiconv for the list of 
+# possible encodings.
+
+DOXYFILE_ENCODING      = UTF-8
+
+# The PROJECT_NAME tag is a single word (or a sequence of words surrounded 
+# by quotes) that should identify the project.
+
+PROJECT_NAME           = ppStack
+
+# The PROJECT_NUMBER tag can be used to enter a project or revision number. 
+# This could be handy for archiving the generated documentation or 
+# if some version control system is used.
+
+PROJECT_NUMBER         = ipp-2.7.dev
+
+# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) 
+# base path where the generated documentation will be put. 
+# If a relative path is entered, it will be relative to the location 
+# where doxygen was started. If left blank the current directory will be used.
+
+OUTPUT_DIRECTORY       = ./docs
+
+# If the CREATE_SUBDIRS tag is set to YES, then doxygen will create 
+# 4096 sub-directories (in 2 levels) under the output directory of each output 
+# format and will distribute the generated files over these directories. 
+# Enabling this option can be useful when feeding doxygen a huge amount of 
+# source files, where putting all generated files in the same directory would 
+# otherwise cause performance problems for the file system.
+
+CREATE_SUBDIRS         = NO
+
+# The OUTPUT_LANGUAGE tag is used to specify the language in which all 
+# documentation generated by doxygen is written. Doxygen will use this 
+# information to generate all constant output in the proper language. 
+# The default language is English, other supported languages are: 
+# Afrikaans, Arabic, Brazilian, Catalan, Chinese, Chinese-Traditional, 
+# Croatian, Czech, Danish, Dutch, Finnish, French, German, Greek, Hungarian, 
+# Italian, Japanese, Japanese-en (Japanese with English messages), Korean, 
+# Korean-en, Lithuanian, Norwegian, Polish, Portuguese, Romanian, Russian, 
+# Serbian, Slovak, Slovene, Spanish, Swedish, and Ukrainian.
+
+OUTPUT_LANGUAGE        = English
+
+# If the BRIEF_MEMBER_DESC tag is set to YES (the default) Doxygen will 
+# include brief member descriptions after the members that are listed in 
+# the file and class documentation (similar to JavaDoc). 
+# Set to NO to disable this.
+
+BRIEF_MEMBER_DESC      = YES
+
+# If the REPEAT_BRIEF tag is set to YES (the default) Doxygen will prepend 
+# the brief description of a member or function before the detailed description. 
+# Note: if both HIDE_UNDOC_MEMBERS and BRIEF_MEMBER_DESC are set to NO, the 
+# brief descriptions will be completely suppressed.
+
+REPEAT_BRIEF           = YES
+
+# This tag implements a quasi-intelligent brief description abbreviator 
+# that is used to form the text in various listings. Each string 
+# in this list, if found as the leading text of the brief description, will be 
+# stripped from the text and the result after processing the whole list, is 
+# used as the annotated text. Otherwise, the brief description is used as-is. 
+# If left blank, the following values are used ("$name" is automatically 
+# replaced with the name of the entity): "The $name class" "The $name widget" 
+# "The $name file" "is" "provides" "specifies" "contains" 
+# "represents" "a" "an" "the"
+
+ABBREVIATE_BRIEF       = 
+
+# If the ALWAYS_DETAILED_SEC and REPEAT_BRIEF tags are both set to YES then 
+# Doxygen will generate a detailed section even if there is only a brief 
+# description.
+
+ALWAYS_DETAILED_SEC    = NO
+
+# If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all 
+# inherited members of a class in the documentation of that class as if those 
+# members were ordinary class members. Constructors, destructors and assignment 
+# operators of the base classes will not be shown.
+
+INLINE_INHERITED_MEMB  = NO
+
+# If the FULL_PATH_NAMES tag is set to YES then Doxygen will prepend the full 
+# path before files name in the file list and in the header files. If set 
+# to NO the shortest path that makes the file name unique will be used.
+
+FULL_PATH_NAMES        = NO
+
+# If the FULL_PATH_NAMES tag is set to YES then the STRIP_FROM_PATH tag 
+# can be used to strip a user-defined part of the path. Stripping is 
+# only done if one of the specified strings matches the left-hand part of 
+# the path. The tag can be used to show relative paths in the file list. 
+# If left blank the directory from which doxygen is run is used as the 
+# path to strip.
+
+STRIP_FROM_PATH        = 
+
+# The STRIP_FROM_INC_PATH tag can be used to strip a user-defined part of 
+# the path mentioned in the documentation of a class, which tells 
+# the reader which header file to include in order to use a class. 
+# If left blank only the name of the header file containing the class 
+# definition is used. Otherwise one should specify the include paths that 
+# are normally passed to the compiler using the -I flag.
+
+STRIP_FROM_INC_PATH    = 
+
+# If the SHORT_NAMES tag is set to YES, doxygen will generate much shorter 
+# (but less readable) file names. This can be useful is your file systems 
+# doesn't support long names like on DOS, Mac, or CD-ROM.
+
+SHORT_NAMES            = NO
+
+# If the JAVADOC_AUTOBRIEF tag is set to YES then Doxygen 
+# will interpret the first line (until the first dot) of a JavaDoc-style 
+# comment as the brief description. If set to NO, the JavaDoc 
+# comments will behave just like regular Qt-style comments 
+# (thus requiring an explicit @brief command for a brief description.)
+
+JAVADOC_AUTOBRIEF      = NO
+
+# If the QT_AUTOBRIEF tag is set to YES then Doxygen will 
+# interpret the first line (until the first dot) of a Qt-style 
+# comment as the brief description. If set to NO, the comments 
+# will behave just like regular Qt-style comments (thus requiring 
+# an explicit \brief command for a brief description.)
+
+QT_AUTOBRIEF           = NO
+
+# The MULTILINE_CPP_IS_BRIEF tag can be set to YES to make Doxygen 
+# treat a multi-line C++ special comment block (i.e. a block of //! or /// 
+# comments) as a brief description. This used to be the default behaviour. 
+# The new default is to treat a multi-line C++ comment block as a detailed 
+# description. Set this tag to YES if you prefer the old behaviour instead.
+
+MULTILINE_CPP_IS_BRIEF = NO
+
+# If the DETAILS_AT_TOP tag is set to YES then Doxygen 
+# will output the detailed description near the top, like JavaDoc.
+# If set to NO, the detailed description appears after the member 
+# documentation.
+
+DETAILS_AT_TOP         = NO
+
+# If the INHERIT_DOCS tag is set to YES (the default) then an undocumented 
+# member inherits the documentation from any documented member that it 
+# re-implements.
+
+INHERIT_DOCS           = YES
+
+# If the SEPARATE_MEMBER_PAGES tag is set to YES, then doxygen will produce 
+# a new page for each member. If set to NO, the documentation of a member will 
+# be part of the file/class/namespace that contains it.
+
+SEPARATE_MEMBER_PAGES  = NO
+
+# The TAB_SIZE tag can be used to set the number of spaces in a tab. 
+# Doxygen uses this value to replace tabs by spaces in code fragments.
+
+TAB_SIZE               = 4
+
+# This tag can be used to specify a number of aliases that acts 
+# as commands in the documentation. An alias has the form "name=value". 
+# For example adding "sideeffect=\par Side Effects:\n" will allow you to 
+# put the command \sideeffect (or @sideeffect) in the documentation, which 
+# will result in a user-defined paragraph with heading "Side Effects:". 
+# You can put \n's in the value part of an alias to insert newlines.
+
+ALIASES                = 
+
+# Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C 
+# sources only. Doxygen will then generate output that is more tailored for C. 
+# For instance, some of the names that are used will be different. The list 
+# of all members will be omitted, etc.
+
+OPTIMIZE_OUTPUT_FOR_C  = NO
+
+# Set the OPTIMIZE_OUTPUT_JAVA tag to YES if your project consists of Java 
+# sources only. Doxygen will then generate output that is more tailored for Java. 
+# For instance, namespaces will be presented as packages, qualified scopes 
+# will look different, etc.
+
+OPTIMIZE_OUTPUT_JAVA   = NO
+
+# If you use STL classes (i.e. std::string, std::vector, etc.) but do not want to 
+# include (a tag file for) the STL sources as input, then you should 
+# set this tag to YES in order to let doxygen match functions declarations and 
+# definitions whose arguments contain STL classes (e.g. func(std::string); v.s. 
+# func(std::string) {}). This also make the inheritance and collaboration 
+# diagrams that involve STL classes more complete and accurate.
+
+BUILTIN_STL_SUPPORT    = NO
+
+# If you use Microsoft's C++/CLI language, you should set this option to YES to
+# enable parsing support.
+
+CPP_CLI_SUPPORT        = NO
+
+# Set the SIP_SUPPORT tag to YES if your project consists of sip sources only. 
+# Doxygen will parse them like normal C++ but will assume all classes use public 
+# instead of private inheritance when no explicit protection keyword is present.
+
+SIP_SUPPORT            = NO
+
+# If member grouping is used in the documentation and the DISTRIBUTE_GROUP_DOC 
+# tag is set to YES, then doxygen will reuse the documentation of the first 
+# member in the group (if any) for the other members of the group. By default 
+# all members of a group must be documented explicitly.
+
+DISTRIBUTE_GROUP_DOC   = NO
+
+# Set the SUBGROUPING tag to YES (the default) to allow class member groups of 
+# the same type (for instance a group of public functions) to be put as a 
+# subgroup of that type (e.g. under the Public Functions section). Set it to 
+# NO to prevent subgrouping. Alternatively, this can be done per class using 
+# the \nosubgrouping command.
+
+SUBGROUPING            = YES
+
+# When TYPEDEF_HIDES_STRUCT is enabled, a typedef of a struct (or union) is 
+# documented as struct with the name of the typedef. So 
+# typedef struct TypeS {} TypeT, will appear in the documentation as a struct 
+# with name TypeT. When disabled the typedef will appear as a member of a file, 
+# namespace, or class. And the struct will be named TypeS. This can typically 
+# be useful for C code where the coding convention is that all structs are 
+# typedef'ed and only the typedef is referenced never the struct's name.
+
+TYPEDEF_HIDES_STRUCT   = NO
+
+#---------------------------------------------------------------------------
+# Build related configuration options
+#---------------------------------------------------------------------------
+
+# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in 
+# documentation are documented, even if no documentation was available. 
+# Private class members and static file members will be hidden unless 
+# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES
+
+EXTRACT_ALL            = NO
+
+# If the EXTRACT_PRIVATE tag is set to YES all private members of a class 
+# will be included in the documentation.
+
+EXTRACT_PRIVATE        = NO
+
+# If the EXTRACT_STATIC tag is set to YES all static members of a file 
+# will be included in the documentation.
+
+EXTRACT_STATIC         = NO
+
+# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs) 
+# defined locally in source files will be included in the documentation. 
+# If set to NO only classes defined in header files are included.
+
+EXTRACT_LOCAL_CLASSES  = YES
+
+# This flag is only useful for Objective-C code. When set to YES local 
+# methods, which are defined in the implementation section but not in 
+# the interface are included in the documentation. 
+# If set to NO (the default) only methods in the interface are included.
+
+EXTRACT_LOCAL_METHODS  = NO
+
+# If this flag is set to YES, the members of anonymous namespaces will be extracted 
+# and appear in the documentation as a namespace called 'anonymous_namespace{file}', 
+# where file will be replaced with the base name of the file that contains the anonymous 
+# namespace. By default anonymous namespace are hidden.
+
+EXTRACT_ANON_NSPACES   = NO
+
+# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all 
+# undocumented members of documented classes, files or namespaces. 
+# If set to NO (the default) these members will be included in the 
+# various overviews, but no documentation section is generated. 
+# This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_MEMBERS     = NO
+
+# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all 
+# undocumented classes that are normally visible in the class hierarchy. 
+# If set to NO (the default) these classes will be included in the various 
+# overviews. This option has no effect if EXTRACT_ALL is enabled.
+
+HIDE_UNDOC_CLASSES     = NO
+
+# If the HIDE_FRIEND_COMPOUNDS tag is set to YES, Doxygen will hide all 
+# friend (class|struct|union) declarations. 
+# If set to NO (the default) these declarations will be included in the 
+# documentation.
+
+HIDE_FRIEND_COMPOUNDS  = NO
+
+# If the HIDE_IN_BODY_DOCS tag is set to YES, Doxygen will hide any 
+# documentation blocks found inside the body of a function. 
+# If set to NO (the default) these blocks will be appended to the 
+# function's detailed documentation block.
+
+HIDE_IN_BODY_DOCS      = NO
+
+# The INTERNAL_DOCS tag determines if documentation 
+# that is typed after a \internal command is included. If the tag is set 
+# to NO (the default) then the documentation will be excluded. 
+# Set it to YES to include the internal documentation.
+
+INTERNAL_DOCS          = NO
+
+# If the CASE_SENSE_NAMES tag is set to NO then Doxygen will only generate 
+# file names in lower-case letters. If set to YES upper-case letters are also 
+# allowed. This is useful if you have classes or files whose names only differ 
+# in case and if your file system supports case sensitive file names. Windows 
+# and Mac users are advised to set this option to NO.
+
+CASE_SENSE_NAMES       = YES
+
+# If the HIDE_SCOPE_NAMES tag is set to NO (the default) then Doxygen 
+# will show members with their full class and namespace scopes in the 
+# documentation. If set to YES the scope will be hidden.
+
+HIDE_SCOPE_NAMES       = NO
+
+# If the SHOW_INCLUDE_FILES tag is set to YES (the default) then Doxygen 
+# will put a list of the files that are included by a file in the documentation 
+# of that file.
+
+SHOW_INCLUDE_FILES     = YES
+
+# If the INLINE_INFO tag is set to YES (the default) then a tag [inline] 
+# is inserted in the documentation for inline members.
+
+INLINE_INFO            = YES
+
+# If the SORT_MEMBER_DOCS tag is set to YES (the default) then doxygen 
+# will sort the (detailed) documentation of file and class members 
+# alphabetically by member name. If set to NO the members will appear in 
+# declaration order.
+
+SORT_MEMBER_DOCS       = YES
+
+# If the SORT_BRIEF_DOCS tag is set to YES then doxygen will sort the 
+# brief documentation of file, namespace and class members alphabetically 
+# by member name. If set to NO (the default) the members will appear in 
+# declaration order.
+
+SORT_BRIEF_DOCS        = NO
+
+# If the SORT_BY_SCOPE_NAME tag is set to YES, the class list will be 
+# sorted by fully-qualified names, including namespaces. If set to 
+# NO (the default), the class list will be sorted only by class name, 
+# not including the namespace part. 
+# Note: This option is not very useful if HIDE_SCOPE_NAMES is set to YES.
+# Note: This option applies only to the class list, not to the 
+# alphabetical list.
+
+SORT_BY_SCOPE_NAME     = NO
+
+# The GENERATE_TODOLIST tag can be used to enable (YES) or 
+# disable (NO) the todo list. This list is created by putting \todo 
+# commands in the documentation.
+
+GENERATE_TODOLIST      = YES
+
+# The GENERATE_TESTLIST tag can be used to enable (YES) or 
+# disable (NO) the test list. This list is created by putting \test 
+# commands in the documentation.
+
+GENERATE_TESTLIST      = YES
+
+# The GENERATE_BUGLIST tag can be used to enable (YES) or 
+# disable (NO) the bug list. This list is created by putting \bug 
+# commands in the documentation.
+
+GENERATE_BUGLIST       = YES
+
+# The GENERATE_DEPRECATEDLIST tag can be used to enable (YES) or 
+# disable (NO) the deprecated list. This list is created by putting 
+# \deprecated commands in the documentation.
+
+GENERATE_DEPRECATEDLIST= YES
+
+# The ENABLED_SECTIONS tag can be used to enable conditional 
+# documentation sections, marked by \if sectionname ... \endif.
+
+ENABLED_SECTIONS       = 
+
+# The MAX_INITIALIZER_LINES tag determines the maximum number of lines 
+# the initial value of a variable or define consists of for it to appear in 
+# the documentation. If the initializer consists of more lines than specified 
+# here it will be hidden. Use a value of 0 to hide initializers completely. 
+# The appearance of the initializer of individual variables and defines in the 
+# documentation can be controlled using \showinitializer or \hideinitializer 
+# command in the documentation regardless of this setting.
+
+MAX_INITIALIZER_LINES  = 30
+
+# Set the SHOW_USED_FILES tag to NO to disable the list of files generated 
+# at the bottom of the documentation of classes and structs. If set to YES the 
+# list will mention the files that were used to generate the documentation.
+
+SHOW_USED_FILES        = YES
+
+# If the sources in your project are distributed over multiple directories 
+# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy 
+# in the documentation. The default is NO.
+
+SHOW_DIRECTORIES       = YES
+
+# The FILE_VERSION_FILTER tag can be used to specify a program or script that 
+# doxygen should invoke to get the current version for each file (typically from the 
+# version control system). Doxygen will invoke the program by executing (via 
+# popen()) the command <command> <input-file>, where <command> is the value of 
+# the FILE_VERSION_FILTER tag, and <input-file> is the name of an input file 
+# provided by doxygen. Whatever the program writes to standard output 
+# is used as the file version. See the manual for examples.
+
+FILE_VERSION_FILTER    = 
+
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+
+# The QUIET tag can be used to turn on/off the messages that are generated 
+# by doxygen. Possible values are YES and NO. If left blank NO is used.
+
+QUIET                  = NO
+
+# The WARNINGS tag can be used to turn on/off the warning messages that are 
+# generated by doxygen. Possible values are YES and NO. If left blank 
+# NO is used.
+
+WARNINGS               = YES
+
+# If WARN_IF_UNDOCUMENTED is set to YES, then doxygen will generate warnings 
+# for undocumented members. If EXTRACT_ALL is set to YES then this flag will 
+# automatically be disabled.
+
+WARN_IF_UNDOCUMENTED   = YES
+
+# If WARN_IF_DOC_ERROR is set to YES, doxygen will generate warnings for 
+# potential errors in the documentation, such as not documenting some 
+# parameters in a documented function, or documenting parameters that 
+# don't exist or using markup commands wrongly.
+
+WARN_IF_DOC_ERROR      = YES
+
+# This WARN_NO_PARAMDOC option can be abled to get warnings for 
+# functions that are documented, but have no documentation for their parameters 
+# or return value. If set to NO (the default) doxygen will only warn about 
+# wrong or incomplete parameter documentation, but not about the absence of 
+# documentation.
+
+WARN_NO_PARAMDOC       = NO
+
+# The WARN_FORMAT tag determines the format of the warning messages that 
+# doxygen can produce. The string should contain the $file, $line, and $text 
+# tags, which will be replaced by the file and line number from which the 
+# warning originated and the warning text. Optionally the format may contain 
+# $version, which will be replaced by the version of the file (if it could 
+# be obtained via FILE_VERSION_FILTER)
+
+WARN_FORMAT            = "$file:$line: $text "
+
+# The WARN_LOGFILE tag can be used to specify a file to which warning 
+# and error messages should be written. If left blank the output is written 
+# to stderr.
+
+WARN_LOGFILE           = 
+
+#---------------------------------------------------------------------------
+# configuration options related to the input files
+#---------------------------------------------------------------------------
+
+# The INPUT tag can be used to specify the files and/or directories that contain 
+# documented source files. You may enter file names like "myfile.cpp" or 
+# directories like "/usr/src/myproject". Separate the files or directories 
+# with spaces.
+
+INPUT                  = ./src
+
+# This tag can be used to specify the character encoding of the source files that 
+# doxygen parses. Internally doxygen uses the UTF-8 encoding, which is also the default 
+# input encoding. Doxygen uses libiconv (or the iconv built into libc) for the transcoding. 
+# See http://www.gnu.org/software/libiconv for the list of possible encodings.
+
+INPUT_ENCODING         = UTF-8
+
+# If the value of the INPUT tag contains directories, you can use the 
+# FILE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
+# and *.h) to filter out the source-files in the directories. If left 
+# blank the following patterns are tested: 
+# *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh *.hxx 
+# *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.py *.f90
+
+FILE_PATTERNS          = *.h
+
+# The RECURSIVE tag can be used to turn specify whether or not subdirectories 
+# should be searched for input files as well. Possible values are YES and NO. 
+# If left blank NO is used.
+
+RECURSIVE              = YES
+
+# The EXCLUDE tag can be used to specify files and/or directories that should 
+# excluded from the INPUT source files. This way you can easily exclude a 
+# subdirectory from a directory tree whose root is specified with the INPUT tag.
+
+EXCLUDE                = 
+
+# The EXCLUDE_SYMLINKS tag can be used select whether or not files or 
+# directories that are symbolic links (a Unix filesystem feature) are excluded 
+# from the input.
+
+EXCLUDE_SYMLINKS       = NO
+
+# If the value of the INPUT tag contains directories, you can use the 
+# EXCLUDE_PATTERNS tag to specify one or more wildcard patterns to exclude 
+# certain files from those directories. Note that the wildcards are matched 
+# against the file with absolute path, so to exclude all test directories 
+# for example use the pattern */test/*
+
+EXCLUDE_PATTERNS       = 
+
+# The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names 
+# (namespaces, classes, functions, etc.) that should be excluded from the output. 
+# The symbol name can be a fully qualified name, a word, or if the wildcard * is used, 
+# a substring. Examples: ANamespace, AClass, AClass::ANamespace, ANamespace::*Test
+
+EXCLUDE_SYMBOLS        = 
+
+# The EXAMPLE_PATH tag can be used to specify one or more files or 
+# directories that contain example code fragments that are included (see 
+# the \include command).
+
+EXAMPLE_PATH           = 
+
+# If the value of the EXAMPLE_PATH tag contains directories, you can use the 
+# EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp 
+# and *.h) to filter out the source-files in the directories. If left 
+# blank all files are included.
+
+EXAMPLE_PATTERNS       = 
+
+# If the EXAMPLE_RECURSIVE tag is set to YES then subdirectories will be 
+# searched for input files to be used with the \include or \dontinclude 
+# commands irrespective of the value of the RECURSIVE tag. 
+# Possible values are YES and NO. If left blank NO is used.
+
+EXAMPLE_RECURSIVE      = NO
+
+# The IMAGE_PATH tag can be used to specify one or more files or 
+# directories that contain image that are included in the documentation (see 
+# the \image command).
+
+IMAGE_PATH             = 
+
+# The INPUT_FILTER tag can be used to specify a program that doxygen should 
+# invoke to filter for each input file. Doxygen will invoke the filter program 
+# by executing (via popen()) the command <filter> <input-file>, where <filter> 
+# is the value of the INPUT_FILTER tag, and <input-file> is the name of an 
+# input file. Doxygen will then use the output that the filter program writes 
+# to standard output.  If FILTER_PATTERNS is specified, this tag will be 
+# ignored.
+
+INPUT_FILTER           = 
+
+# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern 
+# basis.  Doxygen will compare the file name with each pattern and apply the 
+# filter if there is a match.  The filters are a list of the form: 
+# pattern=filter (like *.cpp=my_cpp_filter). See INPUT_FILTER for further 
+# info on how filters are used. If FILTER_PATTERNS is empty, INPUT_FILTER 
+# is applied to all files.
+
+FILTER_PATTERNS        = 
+
+# If the FILTER_SOURCE_FILES tag is set to YES, the input filter (if set using 
+# INPUT_FILTER) will be used to filter the input files when producing source 
+# files to browse (i.e. when SOURCE_BROWSER is set to YES).
+
+FILTER_SOURCE_FILES    = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to source browsing
+#---------------------------------------------------------------------------
+
+# If the SOURCE_BROWSER tag is set to YES then a list of source files will 
+# be generated. Documented entities will be cross-referenced with these sources. 
+# Note: To get rid of all source code in the generated output, make sure also 
+# VERBATIM_HEADERS is set to NO. If you have enabled CALL_GRAPH or CALLER_GRAPH 
+# then you must also enable this option. If you don't then doxygen will produce 
+# a warning and turn it on anyway
+
+SOURCE_BROWSER         = NO
+
+# Setting the INLINE_SOURCES tag to YES will include the body 
+# of functions and classes directly in the documentation.
+
+INLINE_SOURCES         = NO
+
+# Setting the STRIP_CODE_COMMENTS tag to YES (the default) will instruct 
+# doxygen to hide any special comment blocks from generated source code 
+# fragments. Normal C and C++ comments will always remain visible.
+
+STRIP_CODE_COMMENTS    = YES
+
+# If the REFERENCED_BY_RELATION tag is set to YES (the default) 
+# then for each documented function all documented 
+# functions referencing it will be listed.
+
+REFERENCED_BY_RELATION = YES
+
+# If the REFERENCES_RELATION tag is set to YES (the default) 
+# then for each documented function all documented entities 
+# called/used by that function will be listed.
+
+REFERENCES_RELATION    = YES
+
+# If the REFERENCES_LINK_SOURCE tag is set to YES (the default)
+# and SOURCE_BROWSER tag is set to YES, then the hyperlinks from
+# functions in REFERENCES_RELATION and REFERENCED_BY_RELATION lists will
+# link to the source code.  Otherwise they will link to the documentstion.
+
+REFERENCES_LINK_SOURCE = YES
+
+# If the USE_HTAGS tag is set to YES then the references to source code 
+# will point to the HTML generated by the htags(1) tool instead of doxygen 
+# built-in source browser. The htags tool is part of GNU's global source 
+# tagging system (see http://www.gnu.org/software/global/global.html). You 
+# will need version 4.8.6 or higher.
+
+USE_HTAGS              = NO
+
+# If the VERBATIM_HEADERS tag is set to YES (the default) then Doxygen 
+# will generate a verbatim copy of the header file for each class for 
+# which an include is specified. Set to NO to disable this.
+
+VERBATIM_HEADERS       = YES
+
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+
+# If the ALPHABETICAL_INDEX tag is set to YES, an alphabetical index 
+# of all compounds will be generated. Enable this if the project 
+# contains a lot of classes, structs, unions or interfaces.
+
+ALPHABETICAL_INDEX     = NO
+
+# If the alphabetical index is enabled (see ALPHABETICAL_INDEX) then 
+# the COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns 
+# in which this list will be split (can be a number in the range [1..20])
+
+COLS_IN_ALPHA_INDEX    = 5
+
+# In case all classes in a project start with a common prefix, all 
+# classes will be put under the same header in the alphabetical index. 
+# The IGNORE_PREFIX tag can be used to specify one or more prefixes that 
+# should be ignored while generating the index headers.
+
+IGNORE_PREFIX          = 
+
+#---------------------------------------------------------------------------
+# configuration options related to the HTML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_HTML tag is set to YES (the default) Doxygen will 
+# generate HTML output.
+
+GENERATE_HTML          = YES
+
+# The HTML_OUTPUT tag is used to specify where the HTML docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `html' will be used as the default path.
+
+HTML_OUTPUT            = html
+
+# The HTML_FILE_EXTENSION tag can be used to specify the file extension for 
+# each generated HTML page (for example: .htm,.php,.asp). If it is left blank 
+# doxygen will generate files with .html extension.
+
+HTML_FILE_EXTENSION    = .html
+
+# The HTML_HEADER tag can be used to specify a personal HTML header for 
+# each generated HTML page. If it is left blank doxygen will generate a 
+# standard header.
+
+HTML_HEADER            = 
+
+# The HTML_FOOTER tag can be used to specify a personal HTML footer for 
+# each generated HTML page. If it is left blank doxygen will generate a 
+# standard footer.
+
+HTML_FOOTER            = 
+
+# The HTML_STYLESHEET tag can be used to specify a user-defined cascading 
+# style sheet that is used by each HTML page. It can be used to 
+# fine-tune the look of the HTML output. If the tag is left blank doxygen 
+# will generate a default style sheet. Note that doxygen will try to copy 
+# the style sheet file to the HTML output directory, so don't put your own 
+# stylesheet in the HTML output directory as well, or it will be erased!
+
+HTML_STYLESHEET        = 
+
+# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes, 
+# files or namespaces will be aligned in HTML using tables. If set to 
+# NO a bullet list will be used.
+
+HTML_ALIGN_MEMBERS     = YES
+
+# If the GENERATE_HTMLHELP tag is set to YES, additional index files 
+# will be generated that can be used as input for tools like the 
+# Microsoft HTML help workshop to generate a compressed HTML help file (.chm) 
+# of the generated HTML documentation.
+
+GENERATE_HTMLHELP      = NO
+
+# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML 
+# documentation will contain sections that can be hidden and shown after the 
+# page has loaded. For this to work a browser that supports 
+# JavaScript and DHTML is required (for instance Mozilla 1.0+, Firefox 
+# Netscape 6.0+, Internet explorer 5.0+, Konqueror, or Safari).
+
+HTML_DYNAMIC_SECTIONS  = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the CHM_FILE tag can 
+# be used to specify the file name of the resulting .chm file. You 
+# can add a path in front of the file if the result should not be 
+# written to the html output directory.
+
+CHM_FILE               = 
+
+# If the GENERATE_HTMLHELP tag is set to YES, the HHC_LOCATION tag can 
+# be used to specify the location (absolute path including file name) of 
+# the HTML help compiler (hhc.exe). If non-empty doxygen will try to run 
+# the HTML help compiler on the generated index.hhp.
+
+HHC_LOCATION           = 
+
+# If the GENERATE_HTMLHELP tag is set to YES, the GENERATE_CHI flag 
+# controls if a separate .chi index file is generated (YES) or that 
+# it should be included in the master .chm file (NO).
+
+GENERATE_CHI           = NO
+
+# If the GENERATE_HTMLHELP tag is set to YES, the BINARY_TOC flag 
+# controls whether a binary table of contents is generated (YES) or a 
+# normal table of contents (NO) in the .chm file.
+
+BINARY_TOC             = NO
+
+# The TOC_EXPAND flag can be set to YES to add extra items for group members 
+# to the contents of the HTML help documentation and to the tree view.
+
+TOC_EXPAND             = NO
+
+# The DISABLE_INDEX tag can be used to turn on/off the condensed index at 
+# top of each HTML page. The value NO (the default) enables the index and 
+# the value YES disables it.
+
+DISABLE_INDEX          = NO
+
+# This tag can be used to set the number of enum values (range [1..20]) 
+# that doxygen will group on one line in the generated HTML documentation.
+
+ENUM_VALUES_PER_LINE   = 4
+
+# If the GENERATE_TREEVIEW tag is set to YES, a side panel will be
+# generated containing a tree-like index structure (just like the one that 
+# is generated for HTML Help). For this to work a browser that supports 
+# JavaScript, DHTML, CSS and frames is required (for instance Mozilla 1.0+, 
+# Netscape 6.0+, Internet explorer 5.0+, or Konqueror). Windows users are 
+# probably better off using the HTML help feature.
+
+GENERATE_TREEVIEW      = NO
+
+# If the treeview is enabled (see GENERATE_TREEVIEW) then this tag can be 
+# used to set the initial width (in pixels) of the frame in which the tree 
+# is shown.
+
+TREEVIEW_WIDTH         = 250
+
+#---------------------------------------------------------------------------
+# configuration options related to the LaTeX output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will 
+# generate Latex output.
+
+GENERATE_LATEX         = NO
+
+# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `latex' will be used as the default path.
+
+LATEX_OUTPUT           = latex
+
+# The LATEX_CMD_NAME tag can be used to specify the LaTeX command name to be 
+# invoked. If left blank `latex' will be used as the default command name.
+
+LATEX_CMD_NAME         = latex
+
+# The MAKEINDEX_CMD_NAME tag can be used to specify the command name to 
+# generate index for LaTeX. If left blank `makeindex' will be used as the 
+# default command name.
+
+MAKEINDEX_CMD_NAME     = makeindex
+
+# If the COMPACT_LATEX tag is set to YES Doxygen generates more compact 
+# LaTeX documents. This may be useful for small projects and may help to 
+# save some trees in general.
+
+COMPACT_LATEX          = NO
+
+# The PAPER_TYPE tag can be used to set the paper type that is used 
+# by the printer. Possible values are: a4, a4wide, letter, legal and 
+# executive. If left blank a4wide will be used.
+
+PAPER_TYPE             = a4wide
+
+# The EXTRA_PACKAGES tag can be to specify one or more names of LaTeX 
+# packages that should be included in the LaTeX output.
+
+EXTRA_PACKAGES         = 
+
+# The LATEX_HEADER tag can be used to specify a personal LaTeX header for 
+# the generated latex document. The header should contain everything until 
+# the first chapter. If it is left blank doxygen will generate a 
+# standard header. Notice: only use this tag if you know what you are doing!
+
+LATEX_HEADER           = 
+
+# If the PDF_HYPERLINKS tag is set to YES, the LaTeX that is generated 
+# is prepared for conversion to pdf (using ps2pdf). The pdf file will 
+# contain links (just like the HTML output) instead of page references 
+# This makes the output suitable for online browsing using a pdf viewer.
+
+PDF_HYPERLINKS         = NO
+
+# If the USE_PDFLATEX tag is set to YES, pdflatex will be used instead of 
+# plain latex in the generated Makefile. Set this option to YES to get a 
+# higher quality PDF documentation.
+
+USE_PDFLATEX           = NO
+
+# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \\batchmode. 
+# command to the generated LaTeX files. This will instruct LaTeX to keep 
+# running if errors occur, instead of asking the user for help. 
+# This option is also used when generating formulas in HTML.
+
+LATEX_BATCHMODE        = NO
+
+# If LATEX_HIDE_INDICES is set to YES then doxygen will not 
+# include the index chapters (such as File Index, Compound Index, etc.) 
+# in the output.
+
+LATEX_HIDE_INDICES     = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_RTF tag is set to YES Doxygen will generate RTF output 
+# The RTF output is optimized for Word 97 and may not look very pretty with 
+# other RTF readers or editors.
+
+GENERATE_RTF           = NO
+
+# The RTF_OUTPUT tag is used to specify where the RTF docs will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `rtf' will be used as the default path.
+
+RTF_OUTPUT             = rtf
+
+# If the COMPACT_RTF tag is set to YES Doxygen generates more compact 
+# RTF documents. This may be useful for small projects and may help to 
+# save some trees in general.
+
+COMPACT_RTF            = NO
+
+# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated 
+# will contain hyperlink fields. The RTF file will 
+# contain links (just like the HTML output) instead of page references. 
+# This makes the output suitable for online browsing using WORD or other 
+# programs which support those fields. 
+# Note: wordpad (write) and others do not support links.
+
+RTF_HYPERLINKS         = NO
+
+# Load stylesheet definitions from file. Syntax is similar to doxygen's 
+# config file, i.e. a series of assignments. You only have to provide 
+# replacements, missing definitions are set to their default value.
+
+RTF_STYLESHEET_FILE    = 
+
+# Set optional variables used in the generation of an rtf document. 
+# Syntax is similar to doxygen's config file.
+
+RTF_EXTENSIONS_FILE    = 
+
+#---------------------------------------------------------------------------
+# configuration options related to the man page output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_MAN tag is set to YES (the default) Doxygen will 
+# generate man pages
+
+GENERATE_MAN           = YES
+
+# The MAN_OUTPUT tag is used to specify where the man pages will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `man' will be used as the default path.
+
+MAN_OUTPUT             = man
+
+# The MAN_EXTENSION tag determines the extension that is added to 
+# the generated man pages (default is the subroutine's section .3)
+
+MAN_EXTENSION          = .3
+
+# If the MAN_LINKS tag is set to YES and Doxygen generates man output, 
+# then it will generate one additional man file for each entity 
+# documented in the real man page(s). These additional files 
+# only source the real man page, but without them the man command 
+# would be unable to find the correct page. The default is NO.
+
+MAN_LINKS              = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the XML output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_XML tag is set to YES Doxygen will 
+# generate an XML file that captures the structure of 
+# the code including all documentation.
+
+GENERATE_XML           = NO
+
+# The XML_OUTPUT tag is used to specify where the XML pages will be put. 
+# If a relative path is entered the value of OUTPUT_DIRECTORY will be 
+# put in front of it. If left blank `xml' will be used as the default path.
+
+XML_OUTPUT             = xml
+
+# The XML_SCHEMA tag can be used to specify an XML schema, 
+# which can be used by a validating XML parser to check the 
+# syntax of the XML files.
+
+XML_SCHEMA             = 
+
+# The XML_DTD tag can be used to specify an XML DTD, 
+# which can be used by a validating XML parser to check the 
+# syntax of the XML files.
+
+XML_DTD                = 
+
+# If the XML_PROGRAMLISTING tag is set to YES Doxygen will 
+# dump the program listings (including syntax highlighting 
+# and cross-referencing information) to the XML output. Note that 
+# enabling this will significantly increase the size of the XML output.
+
+XML_PROGRAMLISTING     = YES
+
+#---------------------------------------------------------------------------
+# configuration options for the AutoGen Definitions output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_AUTOGEN_DEF tag is set to YES Doxygen will 
+# generate an AutoGen Definitions (see autogen.sf.net) file 
+# that captures the structure of the code including all 
+# documentation. Note that this feature is still experimental 
+# and incomplete at the moment.
+
+GENERATE_AUTOGEN_DEF   = NO
+
+#---------------------------------------------------------------------------
+# configuration options related to the Perl module output
+#---------------------------------------------------------------------------
+
+# If the GENERATE_PERLMOD tag is set to YES Doxygen will 
+# generate a Perl module file that captures the structure of 
+# the code including all documentation. Note that this 
+# feature is still experimental and incomplete at the 
+# moment.
+
+GENERATE_PERLMOD       = NO
+
+# If the PERLMOD_LATEX tag is set to YES Doxygen will generate 
+# the necessary Makefile rules, Perl scripts and LaTeX code to be able 
+# to generate PDF and DVI output from the Perl module output.
+
+PERLMOD_LATEX          = NO
+
+# If the PERLMOD_PRETTY tag is set to YES the Perl module output will be 
+# nicely formatted so it can be parsed by a human reader.  This is useful 
+# if you want to understand what is going on.  On the other hand, if this 
+# tag is set to NO the size of the Perl module output will be much smaller 
+# and Perl will parse it just the same.
+
+PERLMOD_PRETTY         = YES
+
+# The names of the make variables in the generated doxyrules.make file 
+# are prefixed with the string contained in PERLMOD_MAKEVAR_PREFIX. 
+# This is useful so different doxyrules.make files included by the same 
+# Makefile don't overwrite each other's variables.
+
+PERLMOD_MAKEVAR_PREFIX = 
+
+#---------------------------------------------------------------------------
+# Configuration options related to the preprocessor   
+#---------------------------------------------------------------------------
+
+# If the ENABLE_PREPROCESSING tag is set to YES (the default) Doxygen will 
+# evaluate all C-preprocessor directives found in the sources and include 
+# files.
+
+ENABLE_PREPROCESSING   = YES
+
+# If the MACRO_EXPANSION tag is set to YES Doxygen will expand all macro 
+# names in the source code. If set to NO (the default) only conditional 
+# compilation will be performed. Macro expansion can be done in a controlled 
+# way by setting EXPAND_ONLY_PREDEF to YES.
+
+MACRO_EXPANSION        = NO
+
+# If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES 
+# then the macro expansion is limited to the macros specified with the 
+# PREDEFINED and EXPAND_AS_DEFINED tags.
+
+EXPAND_ONLY_PREDEF     = NO
+
+# If the SEARCH_INCLUDES tag is set to YES (the default) the includes files 
+# in the INCLUDE_PATH (see below) will be search if a #include is found.
+
+SEARCH_INCLUDES        = YES
+
+# The INCLUDE_PATH tag can be used to specify one or more directories that 
+# contain include files that are not input files but should be processed by 
+# the preprocessor.
+
+INCLUDE_PATH           = 
+
+# You can use the INCLUDE_FILE_PATTERNS tag to specify one or more wildcard 
+# patterns (like *.h and *.hpp) to filter out the header-files in the 
+# directories. If left blank, the patterns specified with FILE_PATTERNS will 
+# be used.
+
+INCLUDE_FILE_PATTERNS  = 
+
+# The PREDEFINED tag can be used to specify one or more macro names that 
+# are defined before the preprocessor is started (similar to the -D option of 
+# gcc). The argument of the tag is a list of macros of the form: name 
+# or name=definition (no spaces). If the definition and the = are 
+# omitted =1 is assumed. To prevent a macro definition from being 
+# undefined via #undef or recursively expanded use the := operator 
+# instead of the = operator.
+
+PREDEFINED             = 
+
+# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then 
+# this tag can be used to specify a list of macro names that should be expanded. 
+# The macro definition that is found in the sources will be used. 
+# Use the PREDEFINED tag if you want to use a different macro definition.
+
+EXPAND_AS_DEFINED      = 
+
+# If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then 
+# doxygen's preprocessor will remove all function-like macros that are alone 
+# on a line, have an all uppercase name, and do not end with a semicolon. Such 
+# function macros are typically used for boiler-plate code, and will confuse 
+# the parser if not removed.
+
+SKIP_FUNCTION_MACROS   = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to external references   
+#---------------------------------------------------------------------------
+
+# The TAGFILES option can be used to specify one or more tagfiles. 
+# Optionally an initial location of the external documentation 
+# can be added for each tagfile. The format of a tag file without 
+# this location is as follows: 
+#   TAGFILES = file1 file2 ... 
+# Adding location for the tag files is done as follows: 
+#   TAGFILES = file1=loc1 "file2 = loc2" ... 
+# where "loc1" and "loc2" can be relative or absolute paths or 
+# URLs. If a location is present for each tag, the installdox tool 
+# does not have to be run to correct the links.
+# Note that each tag file must have a unique name
+# (where the name does NOT include the path)
+# If a tag file is not located in the directory in which doxygen 
+# is run, you must also specify the path to the tagfile here.
+
+TAGFILES               = 
+
+# When a file name is specified after GENERATE_TAGFILE, doxygen will create 
+# a tag file that is based on the input files it reads.
+
+GENERATE_TAGFILE       = 
+
+# If the ALLEXTERNALS tag is set to YES all external classes will be listed 
+# in the class index. If set to NO only the inherited external classes 
+# will be listed.
+
+ALLEXTERNALS           = NO
+
+# If the EXTERNAL_GROUPS tag is set to YES all external groups will be listed 
+# in the modules index. If set to NO, only the current project's groups will 
+# be listed.
+
+EXTERNAL_GROUPS        = YES
+
+# The PERL_PATH should be the absolute path and name of the perl script 
+# interpreter (i.e. the result of `which perl').
+
+PERL_PATH              = /usr/bin/perl
+
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool   
+#---------------------------------------------------------------------------
+
+# If the CLASS_DIAGRAMS tag is set to YES (the default) Doxygen will 
+# generate a inheritance diagram (in HTML, RTF and LaTeX) for classes with base 
+# or super classes. Setting the tag to NO turns the diagrams off. Note that 
+# this option is superseded by the HAVE_DOT option below. This is only a 
+# fallback. It is recommended to install and use dot, since it yields more 
+# powerful graphs.
+
+CLASS_DIAGRAMS         = YES
+
+# You can define message sequence charts within doxygen comments using the \msc 
+# command. Doxygen will then run the mscgen tool (see http://www.mcternan.me.uk/mscgen/) to 
+# produce the chart and insert it in the documentation. The MSCGEN_PATH tag allows you to 
+# specify the directory where the mscgen tool resides. If left empty the tool is assumed to 
+# be found in the default search path.
+
+MSCGEN_PATH            = 
+
+# If set to YES, the inheritance and collaboration graphs will hide 
+# inheritance and usage relations if the target is undocumented 
+# or is not a class.
+
+HIDE_UNDOC_RELATIONS   = YES
+
+# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is 
+# available from the path. This tool is part of Graphviz, a graph visualization 
+# toolkit from AT&T and Lucent Bell Labs. The other options in this section 
+# have no effect if this option is set to NO (the default)
+
+HAVE_DOT               = NO
+
+# If the CLASS_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect inheritance relations. Setting this tag to YES will force the 
+# the CLASS_DIAGRAMS tag to NO.
+
+CLASS_GRAPH            = YES
+
+# If the COLLABORATION_GRAPH and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for each documented class showing the direct and 
+# indirect implementation dependencies (inheritance, containment, and 
+# class references variables) of the class with other documented classes.
+
+COLLABORATION_GRAPH    = YES
+
+# If the GROUP_GRAPHS and HAVE_DOT tags are set to YES then doxygen 
+# will generate a graph for groups, showing the direct groups dependencies
+
+GROUP_GRAPHS           = YES
+
+# If the UML_LOOK tag is set to YES doxygen will generate inheritance and 
+# collaboration diagrams in a style similar to the OMG's Unified Modeling 
+# Language.
+
+UML_LOOK               = NO
+
+# If set to YES, the inheritance and collaboration graphs will show the 
+# relations between templates and their instances.
+
+TEMPLATE_RELATIONS     = NO
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDE_GRAPH, and HAVE_DOT 
+# tags are set to YES then doxygen will generate a graph for each documented 
+# file showing the direct and indirect include dependencies of the file with 
+# other documented files.
+
+INCLUDE_GRAPH          = YES
+
+# If the ENABLE_PREPROCESSING, SEARCH_INCLUDES, INCLUDED_BY_GRAPH, and 
+# HAVE_DOT tags are set to YES then doxygen will generate a graph for each 
+# documented header file showing the documented files that directly or 
+# indirectly include this file.
+
+INCLUDED_BY_GRAPH      = YES
+
+# If the CALL_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will 
+# generate a call dependency graph for every global function or class method. 
+# Note that enabling this option will significantly increase the time of a run. 
+# So in most cases it will be better to enable call graphs for selected 
+# functions only using the \callgraph command.
+
+CALL_GRAPH             = NO
+
+# If the CALLER_GRAPH, SOURCE_BROWSER and HAVE_DOT tags are set to YES then doxygen will 
+# generate a caller dependency graph for every global function or class method. 
+# Note that enabling this option will significantly increase the time of a run. 
+# So in most cases it will be better to enable caller graphs for selected 
+# functions only using the \callergraph command.
+
+CALLER_GRAPH           = NO
+
+# If the GRAPHICAL_HIERARCHY and HAVE_DOT tags are set to YES then doxygen 
+# will graphical hierarchy of all classes instead of a textual one.
+
+GRAPHICAL_HIERARCHY    = YES
+
+# If the DIRECTORY_GRAPH, SHOW_DIRECTORIES and HAVE_DOT tags are set to YES 
+# then doxygen will show the dependencies a directory has on other directories 
+# in a graphical way. The dependency relations are determined by the #include
+# relations between the files in the directories.
+
+DIRECTORY_GRAPH        = YES
+
+# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images 
+# generated by dot. Possible values are png, jpg, or gif
+# If left blank png will be used.
+
+DOT_IMAGE_FORMAT       = png
+
+# The tag DOT_PATH can be used to specify the path where the dot tool can be 
+# found. If left blank, it is assumed the dot tool can be found in the path.
+
+DOT_PATH               = 
+
+# The DOTFILE_DIRS tag can be used to specify one or more directories that 
+# contain dot files that are included in the documentation (see the 
+# \dotfile command).
+
+DOTFILE_DIRS           = 
+
+# The MAX_DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of 
+# nodes that will be shown in the graph. If the number of nodes in a graph 
+# becomes larger than this value, doxygen will truncate the graph, which is 
+# visualized by representing a node as a red box. Note that doxygen if the number 
+# of direct children of the root node in a graph is already larger than 
+# MAX_DOT_GRAPH_NOTES then the graph will not be shown at all. Also note 
+# that the size of a graph can be further restricted by MAX_DOT_GRAPH_DEPTH.
+
+DOT_GRAPH_MAX_NODES    = 50
+
+# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the 
+# graphs generated by dot. A depth value of 3 means that only nodes reachable 
+# from the root by following a path via at most 3 edges will be shown. Nodes 
+# that lay further from the root node will be omitted. Note that setting this 
+# option to 1 or 2 may greatly reduce the computation time needed for large 
+# code bases. Also note that the size of a graph can be further restricted by 
+# DOT_GRAPH_MAX_NODES. Using a depth of 0 means no depth restriction.
+
+MAX_DOT_GRAPH_DEPTH    = 0
+
+# Set the DOT_TRANSPARENT tag to YES to generate images with a transparent 
+# background. This is disabled by default, which results in a white background. 
+# Warning: Depending on the platform used, enabling this option may lead to 
+# badly anti-aliased labels on the edges of a graph (i.e. they become hard to 
+# read).
+
+DOT_TRANSPARENT        = NO
+
+# Set the DOT_MULTI_TARGETS tag to YES allow dot to generate multiple output 
+# files in one run (i.e. multiple -o and -T options on the command line). This 
+# makes dot run faster, but since only newer versions of dot (>1.8.10) 
+# support this, this feature is disabled by default.
+
+DOT_MULTI_TARGETS      = NO
+
+# If the GENERATE_LEGEND tag is set to YES (the default) Doxygen will 
+# generate a legend page explaining the meaning of the various boxes and 
+# arrows in the dot generated graphs.
+
+GENERATE_LEGEND        = YES
+
+# If the DOT_CLEANUP tag is set to YES (the default) Doxygen will 
+# remove the intermediate dot files that are used to generate 
+# the various graphs.
+
+DOT_CLEANUP            = YES
+
+#---------------------------------------------------------------------------
+# Configuration::additions related to the search engine   
+#---------------------------------------------------------------------------
+
+# The SEARCHENGINE tag specifies whether or not a search engine should be 
+# used. If set to NO the values of all tags below this one will be ignored.
+
+SEARCHENGINE           = NO
Index: /tags/pap_merge_20090108/ppStack/Makefile.am
===================================================================
--- /tags/pap_merge_20090108/ppStack/Makefile.am	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/Makefile.am	(revision 21200)
@@ -0,0 +1,3 @@
+SUBDIRS = src
+
+CLEANFILES = *~ core core.*
Index: /tags/pap_merge_20090108/ppStack/README.txt
===================================================================
--- /tags/pap_merge_20090108/ppStack/README.txt	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/README.txt	(revision 21200)
@@ -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_merge_20090108/ppStack/autogen.sh
===================================================================
--- /tags/pap_merge_20090108/ppStack/autogen.sh	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/autogen.sh	(revision 21200)
@@ -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_merge_20090108/ppStack/configure.ac
===================================================================
--- /tags/pap_merge_20090108/ppStack/configure.ac	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/configure.ac	(revision 21200)
@@ -0,0 +1,35 @@
+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
+  Doxyfile
+])
+AC_OUTPUT
Index: /tags/pap_merge_20090108/ppStack/src/.cvsignore
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/.cvsignore	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/.cvsignore	(revision 21200)
@@ -0,0 +1,10 @@
+.deps
+Makefile
+Makefile.in
+.libs
+*.lo
+*.la
+ppStack
+config.h
+config.h.in
+stamp-h1
Index: /tags/pap_merge_20090108/ppStack/src/Makefile.am
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/Makefile.am	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/Makefile.am	(revision 21200)
@@ -0,0 +1,28 @@
+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		\
+	ppStackSources.c	\
+	ppStackThread.c
+
+noinst_HEADERS = 		\
+	ppStack.h
+
+CLEANFILES = *~
+
+clean-local:
+	-rm -f TAGS
+# Tags for emacs
+tags:
+	etags `find . -name \*.[ch] -print`
Index: /tags/pap_merge_20090108/ppStack/src/ppStack.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStack.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStack.c	(revision 21200)
@@ -0,0 +1,78 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.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;
+    }
+
+    (void) psTraceSetLevel("ppStack", 5);
+
+    if (!ppStackArgumentsSetup(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 (!psphotInit()) {
+        psErrorStackPrint(stderr, "Error initialising psphot.\n");
+        exitValue = PS_EXIT_PROG_ERROR;
+        goto die;
+    }
+
+    if (!ppStackCamera(config)) {
+        psErrorStackPrint(stderr, "Error setting up input files.\n");
+        exitValue = PS_EXIT_CONFIG_ERROR;
+        goto die;
+    }
+
+    if (!ppStackArgumentsParse(config)) {
+        psErrorStackPrint(stderr, "Error reading arguments.\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_merge_20090108/ppStack/src/ppStack.h
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStack.h	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStack.h	(revision 21200)
@@ -0,0 +1,162 @@
+#ifndef PP_STACK_H
+#define PP_STACK_H
+
+#define PPSTACK_RECIPE "PPSTACK"        // Name of the recipe
+#define PPSTACK_INSPECT_PIXELS "PPSTACK.PIXELS" // Name of rejected pixels metadata items
+
+#include <pslib.h>
+#include <psmodules.h>
+
+// Mask values for inputs
+typedef enum {
+    PPSTACK_MASK_MATCH  = 0x01,         // PSF-matching failed
+    PPSTACK_MASK_CHI2   = 0x02,         // Chi^2 too deviant
+    PPSTACK_MASK_REJECT = 0x04,         // Rejection failed
+    PPSTACK_MASK_BAD    = 0x08,         // Bad image (too many pixels rejected)
+    PPSTACK_MASK_ALL    = 0xff          // All errors
+} ppStackMask;
+
+// Thread for stacking chunks
+//
+// Each input file contributes a readout, into which is read a chunk from that file
+typedef struct {
+    psArray *readouts;                  // Input readouts to read and stack
+    bool read;                          // Has the scan been read?
+    bool busy;                          // Is the scan being processed?
+    int firstScan;                      // First row of the chunk to be read for this group
+    int lastScan;                       // Last row of the chunk to be read for this group
+} ppStackThread;
+
+// Allocator
+ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack
+    );
+
+// Data for threads
+typedef struct {
+    psArray *threads;                   // Threads doing stacking
+    int lastScan;                       // Last row that's been read
+    psArray *imageFits;                 // FITS file pointers for images
+    psArray *maskFits;                  // FITS file pointers for masks
+    psArray *weightFits;                // FITS file pointers for weights
+} ppStackThreadData;
+
+// Set up thread data
+ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells
+                                          const psArray *imageNames, // Names of images to read
+                                          const psArray *maskNames, // Names of masks to read
+                                          const psArray *weightNames, // Names of weight maps to read
+                                          const pmConfig *config // Configuration
+    );
+
+// Read chunk into the first available file thread
+ppStackThread *ppStackThreadRead(bool *status, // Status of read
+                                 ppStackThreadData *stack, // Stacks available for reading
+                                 pmConfig *config, // Configuration
+                                 int numChunk, // Chunk number (only for interest)
+                                 int overlap // Overlap between subsequent scans
+    );
+
+// Initialise the threads
+void ppStackThreadInit(void);
+
+// Setup command-line arguments
+bool ppStackArgumentsSetup(int argc, char *argv[], // Command-line arguments
+                           pmConfig *config  // Configuration
+    );
+
+// Parse command-line arguments
+bool ppStackArgumentsParse(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 psArray *psfs   // List of input PSFs
+    );
+
+// Re-do photometry on input sources
+//
+// Photometry for the sources is replaced by what is measured
+bool ppStackInputPhotometry(const pmReadout *ro, // Readout
+                            const psArray *sources, // Sources to photometer
+                            const pmConfig *config // Configuration
+    );
+
+// Perform stacking on a readout
+//
+// Returns an array of pixels to inspect for each input image
+psArray *ppStackReadoutInitial(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 match
+                               const psArray *kernels, // Array with array of kernels used in each PSF match
+                               const psVector *weightings, // Weighting factors for each image
+                               const psVector *addVariance // Additional variance for rejection
+    );
+
+// Thread entry point for ppStackReadoutInitial
+bool ppStackReadoutInitialThread(psThreadJob *job // Job to process
+    );
+
+// Concatenate inspection lists for each input image
+bool ppStackInspect(psThreadJob *job    // Job to process
+    );
+
+// Perform stacking on a readout
+bool ppStackReadoutFinal(const pmConfig *config,   // Configuration
+                         pmReadout *outRO,   // Output readout
+                         const psArray *readouts, // Input readouts
+                         const psArray *rejected, // Array with pixels rejected in each image
+                         const psVector *weightings // Weighting factors for each image
+    );
+
+// Thread entry point for ppStackReadoutFinal
+bool ppStackReadoutFinalThread(psThreadJob *job // Job to process
+    );
+
+// 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
+                  float *chi2,          // Chi^2 from the stamps
+                  float *weighting,     // Stack weighting (1/noise^2)
+                  psArray *sources,     // Array of sources
+                  const pmPSF *psf,     // Target PSF
+                  psRandom *rng,        // Random number generator
+                  const pmConfig *config // Configuration
+    );
+
+
+/// Calculate transparency differences between images
+///
+/// Corrects the source PSF photometry to a common system.  Return the sum of the exposure times.
+float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input
+                                 const pmFPAview *view, // View to readout
+                                 const pmConfig *config // Configuration
+    );
+
+#endif
Index: /tags/pap_merge_20090108/ppStack/src/ppStackArguments.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackArguments.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackArguments.c	(revision 21200)
@@ -0,0 +1,315 @@
+#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"
+            "\tSOURCES(STR):   Sources filename\n"
+            "\tWEIGHTING(F32): Relative weighting to be applied\n",
+            program);
+    fprintf(stderr, "\n");
+    psArgumentHelp(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, config->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(recipe, PS_LIST_TAIL, RECIPENAME, PS_META_REPLACE, 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, config->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(recipe, PS_LIST_TAIL, RECIPENAME, PS_META_REPLACE, 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, config->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; \
+        } \
+    } \
+    psImageMaskType value = pmConfigMaskGet(name, config); \
+    psMetadataAddImageMask(recipe, PS_LIST_TAIL, RECIPENAME, PS_META_REPLACE, NULL, value); \
+}
+
+// Get a string value from the command-line and add it to the target
+static bool valueArgStr(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, PS_META_REPLACE, NULL, value);
+    }
+    return false;
+}
+
+// Get a string value from the command-line or recipe and add it to the target
+static bool valueArgRecipeStr(psMetadata *arguments, // Command-line arguments
+                              psMetadata *recipe, // Recipe
+                              const char *argName, // Argument name in the command-line arguments
+                              const char *mdName, // Name for value in the metadata and recipe
+                              psMetadata *target // Target metadata to which to add value
+                              )
+{
+    psString value = psMetadataLookupStr(NULL, arguments, argName); // Value of interest
+    if (!value) {
+        value = psMetadataLookupStr(NULL, recipe, mdName);
+        if (!value) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s",
+                    mdName, PPSTACK_RECIPE);
+            return false;
+        }
+    }
+    return psMetadataAddStr(target, PS_LIST_TAIL, mdName, PS_META_REPLACE, NULL, value);
+}
+
+bool ppStackArgumentsSetup(int argc, char *argv[], pmConfig *config)
+{
+    assert(config);
+
+    // This capability makes things much faster when debugging
+    bool debugStack = false;            // Read old convolutions to debug the stacking?
+    int argNum = psArgumentGet(argc, argv, "-debug-stack"); // Argument number
+    if (argNum > 0) {
+        debugStack = true;
+        psArgumentRemove(argNum, &argc, argv);
+    }
+
+    pmConfigFileSetsMD(config->arguments, &argc, argv, "PPSTACK.SOURCES", "-sources", NULL);
+
+    if ((argNum = psArgumentGet(argc, argv, "-dumpconfig"))) {
+        psArgumentRemove(argNum, &argc, argv);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "DUMP_CONFIG", PS_META_REPLACE,
+                         "Filename for configuration dump", argv[argNum]);
+        psArgumentRemove(argNum, &argc, argv);
+    }
+
+
+    psMetadata *arguments = config->arguments; // 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, "-combine-sys", 0,
+                     "Relative systematic error in combination", NAN);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-val", 0, "Mask value of input bad pixels", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Mask value to give bad pixels", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-poor", 0, "Mask value to give poor pixels", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-poor-frac", 0, "Fraction of weight for poor pixels", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-deconv-limit", 0, "Maximum deconvolution fraction limit", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-image-rej", 0,
+                     "Pixel rejection fraction threshold for rejecting entire image", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-rows", 0, "Rows to read at once", 0);
+    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", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-psf-radius", 0, "Radius for PSF generation", NAN);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-psf-model", 0, "Model name for PSF generation", NULL);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-order", 0, "Spatial order for PSF generation", 0);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-variance", 0, "Use variance for rejection?", false);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-safe", 0,
+                      "Play safe with small numbers of pixels to combine?", false);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-radius", 0, "Radius (pixels) for matching sources", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-iter", 0, "Maximum iterations for zero point", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-tol", 0, "Tolerance for zero point iterations", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-zp-trans-iter", 0, "Iterations for transparency determination", 0);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-rej", 0, "Rejection threshold for transparency determination", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-trans-thresh", 0, "Threshold for transparency determination", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-rej", 0, "Rejection threshold for stars", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-limit", 0, "Limit on star rejection fraction for successful iteration", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys", 0, "Estimated systematic error", NAN);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-renorm", 0, "Renormalise variance maps?", false);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-mean", 0,
+                     "Statistic for mean in renormalisation", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-renorm-stdev", 0,
+                     "Statistic for stdev in renormalisation", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-renorm-width", 0, "Gaussian width for renormalisation", NAN);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-renorm-num", 0, "Number of samples for renormalisation", 0);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-weight", 0, "Suffix for temporary weight maps", NULL);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-temp-delete", 0,
+                      "Delete temporary files on completion?", false);
+    psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Number of threads to use", 0);
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "-psphot-visual", 0, "psphot visualisation", 0);
+
+    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(arguments, PS_LIST_TAIL, "STAMPS", 0, "Stamps file", stampsName);
+
+    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.");
+        return false;
+    }
+    psMetadataAddMetadata(arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs);
+    psFree(inputs);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "OUTPUT", 0, "Root name of the output image list", argv[2]);
+
+    valueArgStr(arguments, "-stats", "STATS", arguments);
+
+    int numThreads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads
+    if (numThreads > 0 && !psThreadPoolInit(numThreads)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", numThreads);
+        return false;
+    }
+
+    psMetadataAddBool(arguments, PS_LIST_TAIL, "PPSTACK.DEBUG.STACK", 0,
+                      "Read old convolved images to debug stack?", debugStack);
+
+    return true;
+}
+
+bool ppStackArgumentsParse(pmConfig *config)
+{
+    assert(config);
+
+    psMetadata *arguments = config->arguments; // Command-line arguments
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe
+    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("-combine-sys",    "COMBINE.SYS",    F32);
+    VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
+    VALUE_ARG_RECIPE_FLOAT("-image-rej",      "IMAGE.REJ",      F32);
+    VALUE_ARG_RECIPE_FLOAT("-deconv-limit",   "DECONV.LIMIT",   F32);
+    VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-poor-frac",      "POOR.FRACTION",  F32);
+
+    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);
+    valueArgRecipeStr(arguments, recipe, "-mask-bad",  "MASK.BAD",  recipe);
+    valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe);
+
+    VALUE_ARG_RECIPE_FLOAT("-zp-radius",       "ZP.RADIUS",       F32);
+    VALUE_ARG_RECIPE_INT("-zp-iter",           "ZP.ITER",         S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-zp-tol",          "ZP.TOL",          F32);
+    VALUE_ARG_RECIPE_INT("-zp-trans-iter",     "ZP.TRANS.ITER",   S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-zp-trans-rej",    "ZP.TRANS.REJ",    F32);
+    VALUE_ARG_RECIPE_FLOAT("-zp-trans-thresh", "ZP.TRANS.THRESH", F32);
+    VALUE_ARG_RECIPE_FLOAT("-zp-star-rej",     "ZP.STAR.REJ",     F32);
+    VALUE_ARG_RECIPE_FLOAT("-zp-star-limit",   "ZP.STAR.LIMIT",   F32);
+    VALUE_ARG_RECIPE_FLOAT("-zp-star-sys",     "ZP.STAR.SYS",     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);
+    valueArgRecipeStr(arguments, recipe, "-psf-model", "PSF.MODEL", recipe);
+
+    if (psMetadataLookupBool(NULL, arguments, "-photometry") ||
+        psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) {
+        psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE,
+                          "Do photometry on stacked image?", true);
+    }
+
+    if (psMetadataLookupBool(NULL, arguments, "-variance") ||
+        psMetadataLookupBool(NULL, recipe, "VARIANCE")) {
+        psMetadataAddBool(arguments, PS_LIST_TAIL, "VARIANCE", 0, "Use variance for rejection?", true);
+    }
+
+    if (psMetadataLookupBool(NULL, arguments, "-safe") ||
+        psMetadataLookupBool(NULL, recipe, "SAFE")) {
+        psMetadataAddBool(arguments, PS_LIST_TAIL, "SAFE", 0,
+                          "Play safe with small number of pixels to combine?", true);
+    }
+
+    if (psMetadataLookupBool(NULL, arguments, "-renorm") ||
+        psMetadataLookupBool(NULL, recipe, "RENORM")) {
+        psMetadataAddBool(arguments, PS_LIST_TAIL, "RENORM", 0, "Renormalise variance maps?", true);
+    }
+    VALUE_ARG_RECIPE_INT("-renorm-num", "RENORM.NUM", S32, 0);
+    VALUE_ARG_RECIPE_FLOAT("-renorm-width", "RENORM.WIDTH", F32);
+    valueArgRecipeStr(arguments, recipe, "-renorm-mean", "RENORM.MEAN", recipe);
+    valueArgRecipeStr(arguments, recipe, "-renorm-stdev", "RENORM.STDEV", recipe);
+
+    valueArgRecipeStr(arguments, recipe, "-temp-image",  "TEMP.IMAGE",  recipe);
+    valueArgRecipeStr(arguments, recipe, "-temp-mask",   "TEMP.MASK",   recipe);
+    valueArgRecipeStr(arguments, recipe, "-temp-weight", "TEMP.WEIGHT", recipe);
+
+    if (psMetadataLookupBool(NULL, arguments, "-temp-delete") ||
+        psMetadataLookupBool(NULL, recipe, "TEMP.DELETE")) {
+        psMetadataAddBool(arguments, PS_LIST_TAIL, "TEMP.DELETE", 0,
+                          "Delete temporary files on completion?", true);
+    }
+
+    psTrace("ppStack", 1, "Done parsing arguments\n");
+
+    // Dump configuration, now that's it's settled
+    bool status;
+    psString dump_file =  psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
+    if (dump_file) {
+        pmConfigCamerasCull(config, NULL);
+        pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS");
+
+        pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Input file
+        pmConfigDump(config, input->fpa, dump_file);
+    }
+
+    return true;
+
+ERROR:
+    return false;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackCamera.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackCamera.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackCamera.c	(revision 21200)
@@ -0,0 +1,388 @@
+#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
+        psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources
+
+        // 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;
+            }
+        }
+
+        // Add the psf file
+        if (!psf || strlen(psf) == 0) {
+            if (havePSFs) {
+                psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i);
+                return false;
+            }
+        } else {
+            if (!havePSFs && i != 0) {
+                psWarning("PSF not provided for all inputs --- ignoring.");
+            } else {
+                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, "PPSTACK.INPUT.PSF",
+                                                           "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, "PPSTACK.INPUT.PSF is not of type PSF");
+                    return false;
+                }
+                havePSFs = true;
+            }
+        }
+
+        // Add the sources file
+        {
+            psArray *sourcesFiles = psArrayAlloc(1); // Array of filenames for this FPA
+            sourcesFiles->data[0] = psMemIncrRefCounter(sources);
+            psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SOURCES.FILENAMES", PS_META_REPLACE,
+                               "Filenames of SOURCES files", sourcesFiles);
+            psFree(sourcesFiles);
+
+            bool status;
+            pmFPAfile *sourcesFile = pmFPAfileBindFromArgs(&status, imageFile, config,
+                                                           "PPSTACK.INPUT.SOURCES", "SOURCES.FILENAMES");
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)",
+                        i, sources);
+                return false;
+            }
+            if (sourcesFile->type != PM_FPA_FILE_CMF) {
+                psError(PS_ERR_IO, true, "PPSTACK.INPUT.SOURCES is not of type CMF");
+                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
+
+        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");
+    }
+    if (psMetadataLookup(config->arguments, "SOURCES.FILENAMES")) {
+        psMetadataRemoveKey(config->arguments, "SOURCES.FILENAMES");
+    }
+
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i);
+    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.PSF", 0, "Have PSFs available?", havePSFs);
+
+    // Output image
+    pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // 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;
+
+    if (!pmFPAAddSourceFromFormat(fpa, "Stack", output->format)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA.");
+        psFree(fpa);
+        return false;
+    }
+
+    // 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;
+    }
+
+    if (havePSFs) {
+        pmFPAfile *targetPSF = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.TARGET.PSF");
+        if (!targetPSF) {
+            psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.TARGET.PSF"));
+            return false;
+        }
+        if (targetPSF->type != PM_FPA_FILE_PSF) {
+            psError(PS_ERR_IO, true, "PPSTACK.TARGET.PSF is not of type PSF");
+            return false;
+        }
+        targetPSF->save = true;
+    }
+
+    // Output JPEGs
+    pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG1");
+    if (!jpeg1) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG1"));
+        return false;
+    }
+    if (jpeg1->type != PM_FPA_FILE_JPEG) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG1 is not of type JPEG");
+        return false;
+    }
+    jpeg1->save = true;
+    pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG2");
+    if (!jpeg2) {
+        psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG2"));
+        return false;
+    }
+    if (jpeg2->type != PM_FPA_FILE_JPEG) {
+        psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG2 is not of type JPEG");
+        return false;
+    }
+    jpeg2->save = true;
+
+
+    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);
+        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, recipe, "PHOTOMETRY") ||
+        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;
+        }
+    } else {
+        // Output PSF --- only required if photometry is not being performed
+        pmFPAfile *outPSF = pmFPAfileDefineOutputFromFile(config, output, "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;
+    }
+
+    return true;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackLoop.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackLoop.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackLoop.c	(revision 21200)
@@ -0,0 +1,1175 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <unistd.h>
+#include <string.h>
+#include <libgen.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <ppStats.h>
+
+#include "ppStack.h"
+
+//#define TESTING
+
+#define WCS_TOLERANCE 0.001             // Tolerance for WCS
+#define PIXELS_BUFFER 1024              // Initial size of pixel lists
+
+// Here follows lists of files for activation/deactivation at various stages.  Each must be NULL-terminated.
+
+// Files required in preparation for convolution
+static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL };
+
+// Files required for the convolution
+static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", NULL };
+
+// Output files for the combination
+static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
+                                "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL };
+
+// Files for photometry
+static char *photFiles[] = { "PSPHOT.INPUT", "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", NULL };
+
+static void memDump(const char *name)
+{
+    return;
+
+    static int num = 0;                 // Counter, to make files unique and give an idea of sequence
+
+    psString filename = NULL;           // Name of file
+    psStringAppend(&filename, "memdump_%s_%03d.txt", name, num);
+    FILE *memFile = fopen(filename, "w");
+    psFree(filename);
+
+    psMemBlock **leaks = NULL;
+    int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
+    fprintf(memFile, "# MemBlock Size Source\n");
+    unsigned long total = 0;            // Total memory used
+    for (int i = 0; i < numLeaks; i++) {
+        psMemBlock *mb = leaks[i];
+        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize,
+                mb->file, mb->lineno);
+        total += mb->userMemorySize;
+    }
+    fclose(memFile);
+    psFree(leaks);
+
+    fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0));
+    num++;
+}
+
+
+
+// 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
+static void 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);
+
+    for (int i = 0; files[i] != NULL; i++) {
+        pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file
+    }
+    return;
+}
+
+// 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;
+    }
+    psFree(view);
+    return true;
+}
+
+// Write an image to a FITS file
+static bool writeImage(const char *name, // Name of image
+                       psMetadata *header, // Header
+                       const psImage *image, // Image
+                       pmConfig *config // Configuration
+                       )
+{
+    assert(name);
+    assert(image);
+
+    psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name
+    psFits *fits = psFitsOpen(resolved, "w");
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved);
+        psFree(resolved);
+        return false;
+    }
+    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
+        psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved);
+        psFitsClose(fits);
+        psFree(resolved);
+        return false;
+    }
+    psFitsClose(fits);
+    psFree(resolved);
+    return true;
+}
+
+
+bool ppStackLoop(pmConfig *config)
+{
+    assert(config);
+
+    psTimerStart("PPSTACK_TOTAL");
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+    bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files?
+    const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images
+    if (!tempDir) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe");
+        return false;
+    }
+
+    psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, "OUTPUT")); // Name for temporary files
+    const char *tempName = basename(outputName);
+    if (!tempName) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files.");
+        return false;
+    }
+
+    const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images
+    const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks
+    const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for temp weight maps
+    if (!tempImage || !tempMask || !tempWeight) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.WEIGHT in recipe");
+        return false;
+    }
+
+    float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution
+    float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject
+                                                                     // before rejecting entire image
+    float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor"
+
+    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, 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
+
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int overlap = 2 * psMetadataLookupS32(NULL, ppsub,
+                                          "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans
+
+    if (!pmConfigMaskSetBits(NULL, NULL, config)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value.");
+        return false;
+    }
+
+    memDump("start");
+
+    // Preparation iteration: Load the sources, and get a target PSF model
+    psTrace("ppStack", 1, "Determining target PSF....\n");
+    psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching
+    pmPSF *targetPSF = NULL;            // Target PSF
+    float sumExposure = NAN;            // Sum of exposure times
+    if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) {
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, prepareFiles, true);
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            return false;
+        }
+
+        // Generate list of PSFs
+        psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD,
+                                                               "^PPSTACK.INPUT$"); // Iterator
+        psMetadataItem *fileItem; // Item from iteration
+        psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope
+        int numCols = 0, numRows = 0;   // Size of image
+        int index = 0;                  // Index for file
+        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(sourceLists);
+                psFree(fileIter);
+                psFree(psfs);
+                return false;
+            }
+            psfs->data[index] = psMemIncrRefCounter(psf);
+            psMetadataRemoveKey(chip->analysis, "PSPHOT.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(sourceLists);
+                psFree(fileIter);
+                psFree(psfs);
+                return false;
+            }
+            if (numCols == 0 && numRows == 0) {
+                numCols = naxis1;
+                numRows = naxis2;
+            }
+
+            pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources
+            psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources
+            if (!sources) {
+                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout.");
+                return NULL;
+            }
+            sourceLists->data[index] = psMemIncrRefCounter(sources);
+
+            // Re-do photometry if we don't trust the source lists
+            if (psMetadataLookupBool(NULL, recipe, "PHOT")) {
+                psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num);
+                pmFPAfileActivate(config->files, false, NULL);
+                fileActivationSingle(config, convolveFiles, true, index);
+                pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File
+                pmFPAview *view = filesIterateDown(config);
+                if (!view) {
+                    psFree(sourceLists);
+                    psFree(targetPSF);
+                    return false;
+                }
+
+                pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest
+
+                if (!ppStackInputPhotometry(ro, sources, config)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources");
+                    psFree(sourceLists);
+                    psFree(targetPSF);
+                    return false;
+                }
+
+                psFree(view);
+                if (!filesIterateUp(config)) {
+                    psFree(sourceLists);
+                    psFree(targetPSF);
+                    return false;
+                }
+                pmFPAfileActivate(config->files, false, NULL);
+                fileActivation(config, prepareFiles, true);
+            }
+
+            index++;
+        }
+        psFree(fileIter);
+
+        // Zero point calibration
+        sumExposure = ppStackSourcesTransparency(sourceLists, view, config);
+        if (!isfinite(sumExposure) || sumExposure <= 0) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences");
+            psFree(sourceLists);
+            psFree(targetPSF);
+            return false;
+        }
+
+        // Generate target PSF
+        targetPSF = ppStackPSF(config, numCols, numRows, psfs);
+        psFree(psfs);
+        if (!targetPSF) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+            psFree(sourceLists);
+            psFree(view);
+            return false;
+        }
+        psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN,
+                         "Target PSF for stack", targetPSF);
+
+        pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip
+        psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN,
+                         "Target PSF", targetPSF);
+        outChip->data_exists = true;
+
+        psFree(view);
+        if (!filesIterateUp(config)) {
+            return false;
+        }
+    }
+
+    memDump("psf");
+
+    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/%s.%d.%s", tempDir, tempName, i, tempImage);
+        psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask);
+        psStringAppend(&weightName, "%s/%s.%d.%s", tempDir, tempName, i, tempWeight);
+        psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, weightName);
+        imageNames->data[i] = imageName;
+        maskNames->data[i] = maskName;
+        weightNames->data[i] = weightName;
+    }
+    // Free the outputName that we grabbed above to set up tempName.
+    psFree(outputName);
+
+    // Generate convolutions and write them to disk
+    psTrace("ppStack", 1, "Convolving inputs to target PSF....\n");
+    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
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+    int numGood = 0;                    // Number of good frames
+    int numCols = 0, numRows = 0;       // Size of image
+    psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs
+    psVectorInit(inputMask, 0);
+    psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching
+    psVectorInit(matchChi2, NAN);
+    psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2)
+    psVectorInit(weightings, NAN);
+    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++) {
+        psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num);
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivationSingle(config, convolveFiles, true, i);
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
+        pmFPAview *view = filesIterateDown(config);
+        if (!view) {
+            psFree(sourceLists);
+            psFree(targetPSF);
+            psFree(rng);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(fpaList);
+            psFree(cellList);
+            return false;
+        }
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout
+        psFree(view);
+
+        if (numCols == 0 && numRows == 0) {
+            numCols = readout->image->numCols;
+            numRows = readout->image->numRows;
+        } else if (numCols != readout->image->numCols || numRows != readout->image->numRows) {
+            psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d",
+                    readout->image->numCols, readout->image->numRows, numCols, numRows);
+            psFree(sourceLists);
+            psFree(targetPSF);
+            psFree(rng);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(fpaList);
+            psFree(cellList);
+            return false;
+        }
+
+        // Background subtraction, scaling and normalisation is performed automatically by the image matching
+        psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction
+        psTimerStart("PPSTACK_MATCH");
+
+        if (!ppStackMatch(readout, &regions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i],
+                          sourceLists->data[i], targetPSF, rng, config)) {
+            psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i);
+            inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH;
+            psErrorClear();
+            continue;
+        }
+
+        if (stats) {
+            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
+                                                                PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel
+            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK,
+                             "Time to match PSF", psTimerMark("PPSTACK_MATCH"));
+            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK,
+                             "Mean deviation for stamps", kernels->mean);
+            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK,
+                             "RMS deviation for stamps", kernels->rms);
+            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
+                             "Number of stamps", kernels->numStamps);
+
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
+
+        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, config);
+        psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask
+        pmConfigMaskWriteHeader(config, maskHeader);
+        writeImage(maskNames->data[i],   maskHeader, readout->mask, config);
+        psFree(maskHeader);
+        writeImage(weightNames->data[i], hdu->header, readout->weight, config);
+
+        pmCell *inCell = readout->parent; // Input cell
+
+        psListAdd(cellList, PS_LIST_TAIL, inCell);
+        psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent);
+
+        cells->data[i] = psMemIncrRefCounter(inCell);
+        if (!filesIterateUp(config)) {
+            psFree(fpaList);
+            psFree(cellList);
+            return false;
+        }
+        numGood++;
+
+        memDump("match");
+    }
+    psFree(sourceLists);
+    psFree(targetPSF);
+    psFree(rng);
+
+    if (numGood == 0) {
+        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psFree(subKernels);
+        psFree(subRegions);
+        psFree(cells);
+        psFree(inputMask);
+        psFree(matchChi2);
+        psFree(fpaList);
+        psFree(cellList);
+        return false;
+    }
+
+
+    // Update concepts for output
+    {
+        pmFPAview view;                 // View for output
+        view.chip = view.cell = view.readout = 0;
+        pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell
+        pmFPA *outFPA = outCell->parent->parent; // Output FPA
+        pmConceptsAverageFPAs(outFPA, fpaList);
+        pmConceptsAverageCells(outCell, cellList, NULL, NULL, true);
+        psFree(fpaList);
+        psFree(cellList);
+
+        // Update the value of a concept
+#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
+            psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
+            psAssert(item, "Concept should be present"); \
+            psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
+            item->data.F32 = VALUE; \
+        }
+
+        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  sumExposure);
+        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", sumExposure);
+        UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
+    }
+
+
+    // Reject images out-of-hand on the basis of their match chi^2
+    {
+        psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                continue;
+            }
+            values->data.F32[values->n++] = matchChi2->data.F32[i];
+        }
+        assert(values->n == numGood);
+        if (!psVectorSortInPlace(values)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to sort vector.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(cells);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(values);
+            return false;
+        }
+        float median = numGood % 2 ? values->data.F32[numGood / 2] :
+            0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]);
+
+        float rms = 0.74 * (values->data.F32[numGood * 3 / 4] -
+                            values->data.F32[numGood / 4]); // Estimated RMS from interquartile range
+
+        psFree(values);
+
+        float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs)
+        if (isfinite(rej)) {
+            float thresh = median + rej * rms; // Threshold for rejection
+            psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f",
+                     median, rej, rms, thresh);
+
+            int numRej = 0;                 // Number rejected
+            numGood = 0;                    // Number of good images
+            for (int i = 0; i < num; i++) {
+              if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) {
+                    continue;
+                }
+                if (matchChi2->data.F32[i] > thresh) {
+                    numRej++;
+                    inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;
+                    psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f",
+                             i, matchChi2->data.F32[i]);
+                } else {
+                    numGood++;
+                }
+            }
+        }
+    }
+
+    if (numGood == 0) {
+        psError(PS_ERR_UNKNOWN, false, "No good images to combine.");
+        psFree(subKernels);
+        psFree(subRegions);
+        psFree(cells);
+        psFree(inputMask);
+        psFree(matchChi2);
+        return false;
+    }
+
+#ifdef TESTING
+    psTraceSetLevel("psModules.imcombine", 7);
+#endif
+
+
+    // Start threading
+    ppStackThreadInit();
+    ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, weightNames, config);
+
+    psTimerStart("PPSTACK_INITIAL");
+
+    memDump("preinitial");
+
+    // Stack the convolved files
+    psTrace("ppStack", 1, "Initial stack of convolved images....\n");
+    pmReadout *outRO = NULL;            // Output readout
+    pmFPAview *view = NULL;             // View to readout
+    psArray *inspect = NULL;            // Array of arrays of pixels to inspect
+    {
+        int row0, col0;                 // Offset for readout
+        int numCols, numRows;           // Size of readout
+        ppStackThread *thread = stack->threads->data[0]; // Representative thread
+        if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) {
+            psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(stack);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(cells);
+            return false;
+        }
+
+        pmFPAfileActivate(config->files, false, NULL);
+        fileActivation(config, combineFiles, true);
+        view = filesIterateDown(config);
+        if (!view) {
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(stack);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(cells);
+            return false;
+        }
+
+        pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell
+        outRO = pmReadoutAlloc(outCell); // Output readout
+
+        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
+        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+        if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to prepare output.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(stack);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(view);
+            psFree(outRO);
+            psFree(cells);
+            return false;
+        }
+
+        psFree(cells);
+
+        bool status;                    // Status of read
+        int numChunk;                   // Number of chunks
+        for (numChunk = 0; true; numChunk++) {
+            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap);
+            if (!status) {
+                // Something went wrong
+                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(stack);
+                psFree(inputMask);
+                psFree(matchChi2);
+                psFree(view);
+                psFree(outRO);
+                return false;
+            }
+            if (!thread) {
+                // Nothing more to read
+                break;
+            }
+
+            // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)
+            psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start
+            psArrayAdd(job->args, 1, thread);
+            psArrayAdd(job->args, 1, config);
+            psArrayAdd(job->args, 1, outRO);
+            psArrayAdd(job->args, 1, subRegions);
+            psArrayAdd(job->args, 1, subKernels);
+            psArrayAdd(job->args, 1, weightings);
+            psArrayAdd(job->args, 1, matchChi2);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(stack);
+                psFree(inputMask);
+                psFree(weightings);
+                psFree(matchChi2);
+                psFree(view);
+                psFree(outRO);
+                return false;
+            }
+            psFree(job);
+        }
+
+        if (!psThreadPoolWait(false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(stack);
+            psFree(inputMask);
+            psFree(matchChi2);
+            psFree(view);
+            psFree(outRO);
+            return false;
+        }
+        psFree(matchChi2);
+
+        // Harvest the jobs, gathering the inspection lists
+        inspect = psArrayAlloc(num);
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
+            inspect->data[i] = psArrayAllocEmpty(numChunk);
+        }
+        psThreadJob *job;               // Completed job
+        while ((job = psThreadJobGetDone())) {
+            psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0,
+                     "Job has incorrect type: %s", job->type);
+            psArray *results = job->results; // Results of job
+            for (int i = 0; i < num; i++) {
+                if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                    continue;
+                }
+                inspect->data[i] = psArrayAdd(inspect->data[i], 1, results->data[i]);
+            }
+            psFree(job);
+        }
+
+        memDump("initial");
+    }
+
+    if (stats) {
+        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_INITIAL", 0,
+                         "Time to make initial stack", psTimerMark("PPSTACK_INITIAL"));
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Time to make initial stack: %f sec", psTimerClear("PPSTACK_INITIAL"));
+
+    psTimerStart("PPSTACK_REJECT");
+
+    // Pixel rejection
+    psArray *rejected = psArrayAlloc(num);
+    {
+        int numRejected = 0;        // Number of inputs rejected completely
+
+        // Count images rejected out of hand
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                numRejected++;
+            }
+        }
+
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
+
+            psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start
+            psArrayAdd(job->args, 1, inspect);
+            PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(subKernels);
+                psFree(subRegions);
+                psFree(stack);
+                psFree(inputMask);
+                psFree(view);
+                psFree(outRO);
+                psFree(inspect);
+                psFree(rejected);
+                return false;
+            }
+            psFree(job);
+        }
+
+        if (!psThreadPoolWait(true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists.");
+            psFree(subKernels);
+            psFree(subRegions);
+            psFree(stack);
+            psFree(inputMask);
+            psFree(view);
+            psFree(outRO);
+            psFree(inspect);
+            psFree(rejected);
+            return false;
+        }
+
+        if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) {
+            pmStackRejectThreadsInit();
+        }
+
+        int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
+
+        // Reject bad pixels
+        for (int i = 0; i < num; i++) {
+            if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                continue;
+            }
+            psTimerStart("PPSTACK_REJECT");
+
+#ifdef TESTING
+            {
+                psImage *mask = psPixelsToMask(NULL, inspect->data[i],
+                                               psRegionSet(0, numCols - 1, 0, numRows - 1),
+                                               0xff); // Mask image
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "inspect_%03d.fits", i);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, mask, 0, NULL);
+                psFree(mask);
+                psFitsClose(fits);
+            }
+#endif
+
+            psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, stride,
+                                             subRegions->data[i], subKernels->data[i]); // Rejected pixels
+
+#ifdef TESTING
+            {
+                psImage *mask = psPixelsToMask(NULL, reject, psRegionSet(0, numCols - 1, 0, numRows - 1),
+                                               0xff); // Mask image
+                psString name = NULL;           // Name of image
+                psStringAppend(&name, "reject_%03d.fits", i);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, NULL, mask, 0, NULL);
+                psFree(mask);
+                psFitsClose(fits);
+            }
+#endif
+
+            psFree(inspect->data[i]);
+            inspect->data[i] = NULL;
+
+            if (!reject) {
+                psWarning("Rejection on image %d didn't work --- reject entire image.", i);
+                numRejected++;
+                inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT;
+            } else {
+                float frac = reject->n / (float)(numCols * numRows); // Pixel fraction
+                psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)",
+                         reject->n, i, frac * 100.0);
+                if (frac > imageRej) {
+                    psWarning("Image %d rejected completely because rejection fraction (%.3f) "
+                              "exceeds limit (%.3f)", i, frac, imageRej);
+                    psFree(reject);
+                    // reject == NULL means reject image completely
+                    reject = NULL;
+                    inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_BAD;
+                    numRejected++;
+                }
+            }
+
+            // Images without a list of rejected pixels (the list may be empty) are rejected completely
+            rejected->data[i] = reject;
+
+            if (stats) {
+                psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK,
+                                 "Time to perform rejection", psTimerMark("PPSTACK_REJECT"));
+                psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK,
+                                 "Number of pixels rejected", reject ? reject->n : 0);
+            }
+            psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i,
+                     psTimerClear("PPSTACK_REJECT"));
+        }
+
+        psFree(inspect);
+        psFree(subKernels);
+        psFree(subRegions);
+
+        if (numRejected >= num - 1) {
+            psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed.");
+            psFree(stack);
+            psFree(rejected);
+            psFree(inputMask);
+            psFree(view);
+            psFree(outRO);
+            return false;
+        }
+
+        if (stats) {
+            psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_IMAGES", 0,
+                             "Number of images rejected completely", numRejected);
+        }
+   }
+
+    psTimerStart("PPSTACK_FINAL");
+
+    memDump("reject");
+
+    // Final combination
+    psTrace("ppStack", 2, "Final stack of convolved images....\n");
+    {
+        stack->lastScan = 0;            // Reset read
+        bool status;                    // Status of read
+        for (int numChunk = 0; true; numChunk++) {
+            ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0);
+            if (!status) {
+                // Something went wrong
+                psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk);
+                psFree(stack);
+                psFree(rejected);
+                psFree(inputMask);
+                psFree(view);
+                psFree(outRO);
+                return false;
+            }
+            if (!thread) {
+                // Nothing more to read
+                break;
+            }
+
+            // call: ppStackReadoutFinal(config, outRO, readouts, rejected)
+            psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start
+            psArrayAdd(job->args, 1, thread);
+            psArrayAdd(job->args, 1, config);
+            psArrayAdd(job->args, 1, outRO);
+            psArrayAdd(job->args, 1, rejected);
+            psArrayAdd(job->args, 1, weightings);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(stack);
+                psFree(rejected);
+                psFree(inputMask);
+                psFree(view);
+                psFree(outRO);
+                return false;
+            }
+            psFree(job);
+        }
+
+        if (!psThreadPoolWait(true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to do final combination.");
+            psFree(stack);
+            psFree(inputMask);
+            psFree(rejected);
+            psFree(view);
+            psFree(outRO);
+            return false;
+        }
+    }
+
+    memDump("final");
+
+    if (stats) {
+        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_FINAL", 0,
+                         "Time to make final stack", psTimerMark("PPSTACK_FINAL"));
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Time to make final stack: %f sec", psTimerClear("PPSTACK_FINAL"));
+
+    psTrace("ppStack", 2, "Cleaning up after combination....\n");
+
+#if 0
+    // Ensure masked regions really look masked
+    {
+        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
+        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+        if (!pmReadoutMaskApply(outRO, maskBad)) {
+            psWarning("Unable to apply mask");
+        }
+    }
+#endif
+
+    // Close up
+    bool wcsDone = false;           // Have we done the WCS?
+    for (int i = 0; i < num; i++) {
+        if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+            continue;
+        }
+
+        ppStackThread *thread = stack->threads->data[0]; // Representative stack
+        pmReadout *inRO = thread->readouts->data[i]; // Template readout
+        if (inRO && !wcsDone) {
+            // Copy astrometry over
+            wcsDone = true;
+            pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU
+            pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU
+            pmChip *outChip = outRO->parent->parent; // Output chip
+            pmFPA *outFPA = outChip->parent; // Output FPA
+            if (!outHDU || !inHDU) {
+                psWarning("Unable to find HDU at FPA level to copy astrometry.");
+            } else {
+                if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) {
+                    psErrorClear();
+                    psWarning("Unable to read WCS astrometry from input FPA.");
+                    wcsDone = false;
+                } 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.");
+                        wcsDone = false;
+                    }
+                }
+            }
+        }
+
+        if (tempDelete) {
+            psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
+            psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
+            psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
+            if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 ||
+                unlink(weightResolved) == -1) {
+                psWarning("Unable to delete temporary files for image %d", i);
+            }
+            psFree(imageResolved);
+            psFree(maskResolved);
+            psFree(weightResolved);
+        }
+    }
+    psFree(imageNames);
+    psFree(maskNames);
+    psFree(weightNames);
+
+    psFree(inputMask);
+    psFree(stack);
+
+    memDump("cleanup");
+
+    // Generate binned JPEGs
+    {
+        int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level
+        int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level
+
+        // Target cells
+        pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1");
+        pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2");
+        psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+
+        pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts
+        if (!pmReadoutRebin(ro1, outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to bin output.");
+            psFree(ro1);
+            psFree(ro2);
+            psFree(outRO);
+            return false;
+        }
+        psFree(ro1);
+        psFree(ro2);
+    }
+
+    if (psMetadataLookupBool(&mdok, recipe, "RENORM")) {
+        // Statistics for renormalisation
+        psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                "RENORM.MEAN"));
+        psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                 "RENORM.STDEV"));
+        int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples
+        float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot
+        psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+        if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth,
+                                       renormMean, renormStdev, NULL)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
+            psFree(outRO);
+            return false;
+        }
+    }
+
+    if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
+        psTrace("ppStack", 1, "Photometering stacked image....\n");
+
+        psTimerStart("PPSTACK_PHOT");
+
+        fileActivation(config, combineFiles, false);
+        fileActivation(config, photFiles, true);
+        pmFPAview *photView = filesIterateDown(config);
+        if (!ppStackPhotometry(config, outRO, photView)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output.");
+            psFree(outRO);
+            psFree(photView);
+            return false;
+        }
+        psFree(photView);
+
+        fileActivation(config, combineFiles, true);
+
+        if (stats) {
+            pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File
+            pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources
+            psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
+            psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", sources->n);
+            psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0,
+                             "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
+        }
+        psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
+    }
+
+    psThreadPoolFinalize();
+
+    memDump("phot");
+
+    // Statistics on output
+    if (stats) {
+        psTrace("ppStack", 1, "Gathering statistics on stacked image....\n");
+        psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad
+        psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+        ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config);
+    }
+
+    memDump("stats");
+
+    // Put version information into the header
+    pmHDU *hdu = pmHDUFromCell(outRO->parent);
+    if (!hdu) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output.");
+        psFree(outRO);
+        psFree(view);
+        return false;
+    }
+    if (!hdu->header) {
+        hdu->header = psMetadataAlloc();
+    }
+    ppStackVersionMetadata(hdu->header);
+
+    psFree(outRO);
+    psFree(view);
+
+    // Write out summary statistics
+    if (stats) {
+        psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_STACK", 0,
+                         "Time to do photometry", psTimerClear("PPSTACK_TOTAL"));
+
+        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);
+    }
+
+    // Write out the output files
+    if (!filesIterateUp(config)) {
+        return false;
+    }
+
+    memDump("finish");
+
+    return true;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackMatch.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackMatch.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackMatch.c	(revision 21200)
@@ -0,0 +1,658 @@
+#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 MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
+#define FAKE_SIZE 1                     // Size of fake convolution kernel
+#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+                     PM_SOURCE_MODE_CR_LIMIT) // Mask to apply to input sources
+#define FAINT_SOURCE_FRAC 1.0e-4         // Set minimum flux to this fraction of faintest source flux
+
+//#define TESTING                         // Enable debugging output
+
+#ifdef TESTING
+// Read a FITS image
+static bool readImage(psImage **target, // Target for image
+                      const char *name, // Name of FITS file
+                      const pmConfig *config // Configuration
+    )
+{
+    psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
+    psFits *fits = psFitsOpen(resolved, "r");
+    psFree(resolved);
+    if (!fits) {
+        psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name);
+        return false;
+    }
+    psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
+    if (!image) {
+        psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name);
+        psFitsClose(fits);
+        return false;
+    }
+    psFitsClose(fits);
+
+    psFree(*target);
+    *target = image;
+
+    return true;
+}
+#endif
+
+// Get coordinates from a source
+static void coordsFromSource(float *x, float *y, const pmSource *source)
+{
+    assert(x && y);
+    assert(source);
+
+    if (source->modelPSF) {
+        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+        *x = source->peak->xf;
+        *y = source->peak->yf;
+    }
+    return;
+}
+
+static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
+                                   int exclusion // Exclusion zone, pixels
+    )
+{
+    psAssert(sources && sources->n > 0, "Require array of sources");
+    if (exclusion <= 0) {
+        return psMemIncrRefCounter(sources);
+    }
+
+    int num = sources->n;               // Number of sources
+    psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates
+    int numGood = 0;                    // Number of good sources
+    for (int i = 0; i < num; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
+        numGood++;
+    }
+    x->n = y->n = numGood;
+
+    psTree *tree = psTreePlant(2, 2, x, y); // kd tree
+
+    psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
+    psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source
+    int numFiltered = 0;                // Number of filtered sources
+    for (int i = 0; i < num; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source) {
+            continue;
+        }
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        coords->data.F64[0] = xSource;
+        coords->data.F64[1] = ySource;
+
+        long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
+        psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
+                coords->data.F64[0], coords->data.F64[1], numWithin);
+        if (numWithin == 1) {
+            // Only itself inside the exclusion zone
+            filtered = psArrayAdd(filtered, filtered->n, source);
+        } else {
+            numFiltered++;
+        }
+    }
+    psFree(coords);
+    psFree(tree);
+
+    psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
+
+    return filtered;
+}
+
+// Add background into the fake image
+// Based on ppSubBackground()
+static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
+                                     const pmConfig *config // Configuration
+    )
+{
+    psAssert(ro && ro->image, "Need readout image");
+    psAssert(config, "Need configuration");
+
+    psImage *image = ro->image;         // Image of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
+    psAssert(ppStackRecipe, "Need PPSTACK recipe");
+    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
+    psAssert(psphotRecipe, "Need PSPHOT recipe");
+
+    psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
+    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
+
+    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
+    psImageBinning *binning = psMetadataLookupPtr(NULL, psphotRecipe,
+                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
+    psAssert(binning, "Need binning parameters");
+    psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
+    if (!psImageUnbin(unbinned, binned, binning)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model");
+        psFree(unbinned);
+        return NULL;
+    }
+
+    // XXX should these really be here?? (probably not...)
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
+
+     return unbinned;
+}
+
+
+bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
+                  psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
+{
+    assert(readout);
+    assert(regions && !*regions);
+    assert(kernels && !*kernels);
+    assert(config);
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
+
+    // Look up appropriate values from the ppSub recipe
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
+
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
+    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
+    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    bool mdok;                          // Status of MD lookup
+    bool renorm = psMetadataLookupBool(&mdok, recipe, "RENORM"); // Renormalise variances?
+    float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
+    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+
+    if (!pmReadoutMaskNonfinite(readout, maskVal)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout.");
+        return false;
+    }
+
+    pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
+
+    static int numInput = -1;            // Index of input file
+    numInput++;
+#ifdef TESTING
+    // Read previously produced kernel
+    if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
+        const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+        assert(outName);
+        // Read convolution kernel
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+        psFree(filename);
+        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+        psFree(resolved);
+        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
+            psFitsClose(fits);
+            return false;
+        }
+        psFitsClose(fits);
+
+        // Add in variance factor
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+        if (!isfinite(vf)) {
+            vf = 1.0;
+        }
+        if (isfinite(vfItem->data.F32)) {
+            vfItem->data.F32 *= vf;
+        } else {
+            vfItem->data.F32 = vf;
+        }
+
+        // Read image, mask, weight
+        const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
+        const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
+        const char *tempWeight = psMetadataLookupStr(NULL, recipe, "TEMP.WEIGHT"); // Suffix for weight map
+        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
+        psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
+        psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
+        psStringAppend(&weightName, "%s.%d.%s", outName, numInput, tempWeight);
+
+        if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
+            !readImage(&readout->weight, weightName, config)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced image.");
+            psFree(imageName);
+            psFree(maskName);
+            psFree(weightName);
+            return false;
+        }
+        psFree(imageName);
+        psFree(maskName);
+        psFree(weightName);
+
+        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
+                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+        pmSubtractionAnalysis(readout->analysis, kernels, region,
+                              readout->image->numCols, readout->image->numRows);
+    } else {
+#endif
+
+        // Normal operations here
+        if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
+            assert(psf);
+            assert(sources);
+
+            int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
+            float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs
+            float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing
+            int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size
+            float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps
+            int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
+            int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
+            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
+            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
+            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
+                psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
+            psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
+            psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
+            int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius
+            int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order
+            int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel
+            float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction
+            bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters?
+            float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search
+            float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search
+            float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search
+            float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search
+            int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search
+            float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
+
+            // These values are specified specifically for stacking
+            const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
+
+            psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
+            if (optimum) {
+                optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
+            }
+
+#if 0
+            // Testing the normalisation of the fake image
+            {
+                pmReadout *test = pmReadoutAlloc(NULL); // Test readout
+                psArray *sources = psArrayAlloc(1);     // Array of sources
+                pmSource *source = pmSourceAlloc();     // Source
+                sources->data[0] = source;
+                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
+                source->psfMag = -13.0;
+                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
+                float sum = 0.0;
+                for (int y = 0; y < test->image->numRows; y++) {
+                    for (int x = 0; x < test->image->numCols; x++) {
+                        sum += test->image->data.F32[y][x];
+                    }
+                }
+                fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum));
+
+                psFits *fits = psFitsOpen("testphot.fits", "w");
+                psFitsWriteImage(fits, NULL, test->image, 0, NULL);
+                psFitsClose(fits);
+                exit(0);
+            }
+#endif
+
+            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
+
+            // For the sake of stamps, remove nearby sources
+            psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
+
+            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
+                                          stampSources, NULL, NULL, psf, NAN, footprint + size,
+                                          false, true)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
+                psFree(fake);
+                psFree(optWidths);
+                psFree(output);
+                return false;
+            }
+
+            // Add the background into the target image
+            psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
+            psBinaryOp(fake->image, fake->image, "+", bgImage);
+            psFree(bgImage);
+
+#ifdef TESTING
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "fake_%03d.fits", numInput);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
+                psFitsClose(fits);
+            }
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "real_%03d.fits", numInput);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
+                psFitsClose(fits);
+            }
+#endif
+
+            // Renormalise the variances if desired
+            if (renorm) {
+                // Statistics for renormalisation
+                psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                        "RENORM.MEAN"));
+                psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                         "RENORM.STDEV"));
+
+                if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
+                    psFree(output);
+                    return false;
+                }
+            }
+
+            if (threads > 0) {
+                pmSubtractionThreadsInit(readout, fake);
+            }
+
+            // Do the image matching
+            if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
+                                    threshold, stampSources, stampsName, type, size, order, widths, orders,
+                                    inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder,
+                                    optThresh, iter, rej, sysError, maskVal, maskBad, maskPoor, poorFrac,
+                                    badFrac, PM_SUBTRACTION_MODE_1)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
+                psFree(fake);
+                psFree(optWidths);
+                psFree(stampSources);
+                psFree(output);
+                return false;
+            }
+
+#ifdef TESTING
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "conv_%03d.fits", numInput);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
+                psFitsClose(fits);
+            }
+            {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
+                psString name = NULL;
+                psStringAppend(&name, "diff_%03d.fits", numInput);
+                psFits *fits = psFitsOpen(name, "w");
+                psFree(name);
+                psBinaryOp(fake->image, output->image, "-", fake->image);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
+                psFitsClose(fits);
+            }
+#endif
+
+            psFree(fake);
+            psFree(optWidths);
+            psFree(stampSources);
+
+            if (threads > 0) {
+                pmSubtractionThreadsFinalize(readout, fake);
+            }
+
+            // Set the variance factor
+            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+            float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
+            if (!isfinite(vf)) {
+                vf = 1.0;
+            }
+            if (isfinite(vfItem->data.F32)) {
+                vfItem->data.F32 *= vf;
+            } else {
+                vfItem->data.F32 = vf;
+            }
+
+            // 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);
+        } else {
+            // Fake the convolution
+            psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
+                             PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
+            psFree(region);
+            pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
+                                                                     PM_SUBTRACTION_MODE_1);
+            // Set solution to delta function
+            kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
+            psVectorInit(kernels->solution1, 0.0);
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = 1.0;
+            psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
+                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
+            psFree(kernels);
+        }
+
+#ifdef TESTING
+        // Write convolution kernel
+        {
+            const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
+            assert(outName);
+
+            psString filename = NULL;   // Output filename
+            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+            psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
+            psFree(filename);
+            psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
+            psFree(resolved);
+            pmReadoutWriteSubtractionKernels(output, fits);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
+    readout->analysis = psMetadataCopy(readout->analysis, output->analysis);
+
+// 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);
+
+    // Record chi^2
+    {
+        *chi2 = 0.0;
+        int num = 0;                    // Number of measurements of chi^2
+        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);
+            pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
+            *chi2 += kernels->mean;
+            num++;
+        }
+        psFree(iter);
+        *chi2 /= num;
+    }
+
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconv = psMetadataLookupF32(NULL, output->analysis,
+                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                  deconv, deconvLimit);
+        psFree(output);
+        return NULL;
+    }
+
+    // Renormalise the variances if desired
+    if (renorm) {
+        // Statistics for renormalisation
+        psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                "RENORM.MEAN"));
+        psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe,
+                                                                                 "RENORM.STDEV"));
+
+        if (!pmReadoutWeightRenormPixels(readout, maskBad, renormMean, renormStdev, rng)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances.");
+            psFree(output);
+            return false;
+        }
+    }
+
+    // Ensure the background value is zero
+    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
+    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
+        psWarning("Can't measure background for image.");
+        psErrorClear();
+    } else {
+        psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)",
+                 psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+        (void)psBinaryOp(readout->image, readout->image, "-",
+                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
+        *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
+        psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
+                         "Weighting by 1/noise^2 for stack",
+                         1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)));
+    }
+    psFree(bg);
+
+
+#if 0
+#define RADIUS 10                       // Radius of photometry
+#define MIN_ERR 0.05                    // Minimum photometric error, mag
+#define MAX_MAG -13                     // Maximum magnitude for source
+
+    // Ensure the normalisation is correct
+    // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
+    int numSources = sources->n;        // Number of sources
+    psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
+    psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
+    psVectorInit(ratioMask, 0xFF);
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Mask of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = sources->data[i]; // Source of interest
+        if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
+            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
+            continue;
+        }
+
+        float xSrc, ySrc;              // Source coordinates
+        coordsFromSource(&xSrc, &ySrc, source);
+        int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
+        int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
+        int numPix = 0;                 // Number of pixels
+        float sum = 0.0;                // Sum of pixels
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
+                    continue;
+                }
+                sum += image->data.F32[y][x];
+                numPix++;
+            }
+        }
+        if (sum >= 0 && numPix > 0) {
+            float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
+            ratio->data.F32[i] = mag - source->psfMag;
+            ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
+        }
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
+        psWarning("Unable to measure normalisation --- assuming correct.");
+    } else {
+        psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
+                 stats->robustMedian, stats->robustStdev);
+        float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
+        psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+    }
+    psFree(stats);
+    psFree(ratio);
+    psFree(ratioMask);
+#endif
+
+#ifdef TESTING
+    {
+        pmHDU *hdu = pmHDUFromCell(readout->parent);
+        psString name = NULL;
+        psStringAppend(&name, "convolved_%03d.fits", numInput);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    psFree(output);
+
+    return true;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackPSF.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackPSF.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackPSF.c	(revision 21200)
@@ -0,0 +1,31 @@
+#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 psArray *psfs)
+{
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF
+    float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF
+    const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF
+    int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF
+
+    // Solve for the target PSF
+    pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel,
+                               psfOrder, psfOrder);
+    if (!psf) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF.");
+        return NULL;
+    }
+
+    return psf;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackPhotometry.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackPhotometry.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackPhotometry.c	(revision 21200)
@@ -0,0 +1,157 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
+                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
+
+bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+
+    float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
+        return false;
+    }
+    float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
+        return false;
+    }
+    float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
+        return false;
+    }
+
+    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
+    if (!mdok || !maskValStr) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
+        return false;
+    }
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
+
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
+
+    // Measure background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    float bg = stats->robustMedian; // Background level
+    psFree(stats);
+    psFree(rng);
+
+    // Photometer sources
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int numSources = sources->n; // Number of sources
+    int numGood = 0;            // Number of good sources
+    float maxMag = -INFINITY;   // Maximum magnitude
+    for (int j = 0; j < numSources; j++) {
+        pmSource *source = sources->data[j]; // Source of interest
+        if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) {
+            source->psfMag = NAN;
+            continue;
+        }
+        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
+        int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
+        // Range of integration
+        int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius);
+        int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius);
+
+        // Integrate
+        double sumImage = 0.0, sumKernel = 0.0; // Sums from integration
+        int numBadPix = 0;         // Number of bad pixels
+        for (int v = yMin; v <= yMax; v++) {
+            float dy2 = PS_SQR(y - v); // Distance from centroid
+            for (int u = xMin; u <= xMax; u++) {
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) {
+                    numBadPix++;
+                    continue;
+                }
+                float dx2 = PS_SQR(x - u); // Distance from centroid
+                double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value
+                sumImage += (image->data.F32[v][u] - bg) * kernel;
+                sumKernel += kernel;
+            }
+        }
+
+        if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) {
+            source->psfMag = NAN;
+        } else {
+            source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius));
+            if (isfinite(source->psfMag)) {
+                numGood++;
+                maxMag = PS_MAX(maxMag, source->psfMag);
+            }
+        }
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
+
+    return true;
+}
+
+
+
+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 (psMetadataLookupBool(NULL, config->arguments, "-psphot-visual")) {
+        psphotSetVisual(true);
+    }
+
+    // Need to ensure aperture residual is not calculated
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
+    psMetadataItem *item = psMetadataLookup(recipe, "MEASURE.APTREND"); // Item determining aptrend
+    if (!item) {
+        psWarning("Unable to find MEASURE.APTREND in psphot recipe");
+        psErrorClear();
+    } else {
+        item->data.B = false;
+    }
+
+    // set maskValue and markValue in the psphot recipe
+    psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
+    psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
+    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
+    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue);
+
+    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();
+    }
+
+    if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
+        !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") ||
+        !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files.");
+        return false;
+    }
+
+    pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
+
+    return true;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackReadout.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackReadout.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackReadout.c	(revision 21200)
@@ -0,0 +1,279 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include <psmodules.h>
+#include <psphot.h>
+
+#include "ppStack.h"
+
+//#define TESTING                  // Write debugging output?
+
+bool ppStackReadoutInitialThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;          // Arguments
+    ppStackThread *thread = args->data[0]; // Thread
+    pmConfig *config = args->data[1];   // Configuration
+    pmReadout *outRO = args->data[2];   // Output readout
+    psArray *subRegions = args->data[3]; // Regions for PSF-matching
+    psArray *subKernels = args->data[4]; // Kernels for PSF-matching
+    psVector *weightings = args->data[5]; // Weightings (1/noise^2) for each image
+    psVector *addVariance = args->data[6]; // Additional variance when rejecting
+
+    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels,
+                                             weightings, addVariance);
+
+    job->results = inspect;
+    thread->busy = false;
+
+    return inspect ? true : false;
+}
+
+bool ppStackReadoutFinalThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;          // Arguments
+    ppStackThread *thread = args->data[0]; // Thread
+    pmConfig *config = args->data[1];   // Configuration
+    pmReadout *outRO = args->data[2];   // Output readout
+    psArray *rejected = args->data[3];  // Rejected pixels
+    psVector *weightings = args->data[4];  // Weighting (1/noise^2) for each image
+
+    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected,
+                                      weightings); // Status of operation
+
+    thread->busy = false;
+
+    return status;
+}
+
+
+bool ppStackInspect(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+
+    psArray *args = job->args;  // Input arguments
+    psArray *inspect = args->data[0]; // Array of pixel arrays
+    int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest
+
+    psArray *inputs = inspect->data[index]; // Array of interest
+    psPixels *output = NULL;    // Output pixel list
+    for (int i = 0; i < inputs->n; i++) {
+        psPixels *input = inputs->data[i]; // Input pixel list
+        if (!input || input->n == 0) {
+            continue;
+        }
+        output = psPixelsConcatenate(output, input);
+    }
+
+    if (!output) {
+        // If there are no pixels to inspect, then just fake it
+        output = psPixelsAllocEmpty(0);
+    }
+
+    psFree(inputs);
+    inspect->data[index] = output;
+
+    return true;
+}
+
+
+psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                               const psArray *regions, const psArray *kernels, const psVector *weightings,
+                               const psVector *addVariance)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(regions);
+    assert(kernels);
+    assert(readouts->n == regions->n);
+    assert(regions->n == kernels->n);
+    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
+    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
+    static int sectionNum = 0;          // Section number; for debugging outputs
+
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+    int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
+    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
+    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
+    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
+    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
+
+    psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
+    int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
+
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    for (int i = 0; i < num; i++) {
+        pmReadout *ro = readouts->data[i];
+        if (!ro) {
+            // Bad image
+            continue;
+        }
+
+#if 0
+        // This doesn't seem to work, so getting the weightings directly from a vector
+        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
+        if (!mdok || !isfinite(weighting)) {
+            psWarning("No weighting supplied for image %d --- set to unity.", i);
+            weighting = 1.0;
+        } else {
+            psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting);
+        }
+#endif
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
+    }
+
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys,
+                        true, useVariance, safe)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
+        psFree(stack);
+        return false;
+    }
+
+#ifdef TESTING
+    {
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_initial_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    // Save list of pixels to inspect
+    psArray *inspect = psArrayAlloc(num); // List of pixels to inspect
+    for (int i = 0; i < num; i++) {
+        pmStackData *data = stack->data[i]; // Data for this image
+        if (!data) {
+            continue;
+        }
+        pmReadout *readout = data->readout; // Readout of interest
+        if (!readout) {
+            continue;
+        }
+        inspect->data[i] = psMemIncrRefCounter(data->inspect);
+    }
+    psFree(stack);
+
+    sectionNum++;
+
+    return inspect;
+}
+
+
+
+bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
+                         const psArray *rejected, const psVector *weightings)
+{
+    assert(config);
+    assert(outRO);
+    assert(readouts);
+    assert(rejected);
+    assert(readouts->n == rejected->n);
+    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
+
+    static int sectionNum = 0;
+
+    // Get the recipe values
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
+
+    psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in
+    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
+    psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
+    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    int num = readouts->n;              // Number of inputs
+    psArray *stack = psArrayAlloc(num); // Array for stacking
+
+    int numGood = num;                  // Number of good inputs: images that haven't been completely rejected
+    for (int i = 0; i < num; i++) {
+        if (!rejected->data[i]) {
+            // Image completely rejected
+            numGood--;
+            continue;
+        }
+
+        pmReadout *ro = readouts->data[i];
+        assert(ro);
+
+#if 0
+        // This doesn't seem to work, so getting the weightings directly from a vector
+        bool mdok;                      // Status of MD lookup
+        float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight
+        if (!mdok || !isfinite(weighting)) {
+            psWarning("No WEIGHTING supplied for image %d --- set to unity.", i);
+            weighting = 1.0;
+        }
+#endif
+
+        // Ensure there is a mask, or pmStackCombine will complain
+        if (!ro->mask) {
+            ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_IMAGE_MASK);
+            psImageInit(ro->mask, 0);
+        }
+
+        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN);
+        data->reject = psMemIncrRefCounter(rejected->data[i]);
+        stack->data[i] = data;
+    }
+
+    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN,
+                        numGood != num, useVariance, false)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
+        psFree(stack);
+        return false;
+    }
+
+#ifdef TESTING
+    {
+        psString name = NULL;           // Name of image
+        psStringAppend(&name, "combined_final_%03d.fits", sectionNum);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
+    pmCell *outCell = outRO->parent;    // Output cell
+    pmChip *outChip = outCell->parent;  // Output chip
+
+    outRO->data_exists = true;
+    outCell->data_exists = true;
+    outChip->data_exists = true;
+
+    psFree(stack);
+
+    sectionNum++;
+
+    return true;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackSources.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackSources.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackSources.c	(revision 21200)
@@ -0,0 +1,193 @@
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+//#define TESTING                         // Enable debugging output
+
+// Size of fake image; set by hand because it's trouble to get it from other places
+#define FAKE_COLS 4861
+#define FAKE_ROWS 4913
+
+float ppStackSourcesTransparency(const psArray *sourceLists, const pmFPAview *view, const pmConfig *config)
+{
+    PS_ASSERT_ARRAY_NON_NULL(sourceLists, NAN);
+    PS_ASSERT_PTR_NON_NULL(view, NAN);
+    PS_ASSERT_PTR_NON_NULL(config, NAN);
+
+#ifdef TESTING
+    {
+        // Deliberately induce a major transparency difference
+        psArray *sources = sourceLists->data[1]; // Sources to correct
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i]; // Source of interest
+            if (!source) {
+                continue;
+            }
+            source->psfMag += 1.0;
+        }
+    }
+#endif
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    float radius = psMetadataLookupF32(NULL, recipe, "ZP.RADIUS"); // Radius (pixels) for matching sources
+    int iter = psMetadataLookupS32(NULL, recipe, "ZP.ITER"); // Maximum iterations
+    float tol = psMetadataLookupF32(NULL, recipe, "ZP.TOL"); // Tolerance for zero point iterations
+    int transIter = psMetadataLookupS32(NULL, recipe, "ZP.TRANS.ITER"); // Iterations for transparency
+    float transRej = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.REJ");// Rejection threshold for transparency
+    float transThresh = psMetadataLookupF32(NULL, recipe, "ZP.TRANS.THRESH"); // Threshold for transparency
+    float starRej = psMetadataLookupF32(NULL, recipe, "ZP.STAR.REJ"); // Rejection threshold for stars
+    float starLimit = psMetadataLookupF32(NULL, recipe, "ZP.STAR.LIMIT"); // Limit on star rejection fraction
+    float starSys = psMetadataLookupF32(NULL, recipe, "ZP.STAR.SYS"); // Estimated systematic error
+
+    psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms
+    if (!airmassZP) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe.");
+        return NAN;
+    }
+
+    int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs
+    psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n);
+
+    psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image
+    const char *filter = NULL;          // Filter name
+    float airmassTerm = NAN;            // Airmass term
+    float sumExpTime = 0.0;             // Sum of the exposure time
+    for (int i = 0; i < num; i++) {
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest
+        pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest
+
+#ifdef TESTING
+        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout
+        pmPSF *psf = psMetadataLookupPtr(NULL, config->arguments, "PSF.TARGET"); // PSF for fake image
+        pmReadoutFakeFromSources(fake, FAKE_COLS, FAKE_ROWS, sourceLists->data[i], NULL, NULL, psf, 5, 0, false, true);
+        psString name = NULL;
+        psStringAppend(&name, "start_%03d.fits", i);
+        psFits *fits = psFitsOpen(name, "w");
+        psFree(name);
+        psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+        psFitsClose(fits);
+        psFree(fake);
+#endif
+
+
+        float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time
+        float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass
+        const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name
+        if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 ||
+            !expFilter || strlen(expFilter) == 0) {
+            psError(PS_ERR_UNEXPECTED_NULL, false,
+                    "Unable to find exposure time (%f), airmass (%f) or filter (%s)",
+                    exptime, airmass, expFilter);
+            psFree(zp);
+            return NAN;
+        }
+
+        if (!filter) {
+            filter = expFilter;
+            airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter);
+            if (!isfinite(airmassTerm)) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, false,
+                        "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter);
+                psFree(zp);
+                return NAN;
+            }
+        } else if (strcmp(filter, expFilter) != 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter);
+            psFree(zp);
+            return NAN;
+        }
+
+        zp->data.F32[i] = airmassTerm * airmass - 2.5 * log10(exptime);
+        sumExpTime += exptime;
+    }
+
+    psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+    if (!matches) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
+        psFree(zp);
+        return NAN;
+    }
+    psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
+                                           transThresh, starRej, starSys); // Transparencies for each image
+
+#ifdef TESTING
+    {
+        // Dump the corrected magnitudes
+        FILE *outMatches = fopen("source_match.dat", "w"); // Output matches
+        psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes for each star
+        psVector *magErr = psVectorAlloc(num, PS_TYPE_F32); // Errors for each star
+        for (int i = 0; i < matches->n; i++) {
+            pmSourceMatch *match = matches->data[i]; // Match of interest
+            psVectorInit(mag, NAN);
+            psVectorInit(magErr, NAN);
+            for (int j = 0; j < match->num; j++) {
+                if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) {
+                    continue;
+                }
+                int index = match->image->data.U32[j]; // Image index
+                mag->data.F32[index] = match->mag->data.F32[j] - zp->data.F32[index] - trans->data.F32[index];
+                magErr->data.F32[index] = match->magErr->data.F32[j];
+            }
+            for (int j = 0; j < num; j++) {
+                fprintf(outMatches, "%f (%f) ", mag->data.F32[j], magErr->data.F32[j]);
+            }
+            fprintf(outMatches, "\n");
+        }
+        psFree(mag);
+        psFree(magErr);
+        fclose(outMatches);
+    }
+#endif
+
+    psFree(matches);
+    if (!trans) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies");
+        return NAN;
+    }
+
+    // M = m + c0 + c1 * airmass - 2.5log(t) + transparency
+    // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0
+    // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1
+    // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0)
+    // We don't need to know the magnitude zero point for the filter, since it cancels out
+    for (int i = 0; i < num; i++) {
+        psArray *sources = sourceLists->data[i]; // Sources of interest
+        float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i];
+        psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr);
+
+        for (int j = 0; j < sources->n; j++) {
+            pmSource *source = sources->data[j]; // Source of interest
+            if (!source) {
+                continue;
+            }
+            source->psfMag += magCorr;
+        }
+    }
+    psFree(trans);
+
+#ifdef TESTING
+    // Double check: all transparencies should be zero
+    {
+        psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches
+        if (!matches) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to match sources");
+            psFree(zp);
+            return NAN;
+        }
+        psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej,
+                                               transThresh, starRej, starSys);
+        for (int i = 0; i < num; i++) {
+            fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]);
+        }
+        psFree(trans);
+    }
+#endif
+
+    return sumExpTime;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackThread.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackThread.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackThread.c	(revision 21200)
@@ -0,0 +1,258 @@
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <unistd.h>
+#include <pslib.h>
+#include <psmodules.h>
+
+#include "ppStack.h"
+
+#define THREAD_WAIT 10000               // Microseconds to wait if thread is not available
+
+
+static void stackThreadFree(ppStackThread *thread)
+{
+    psFree(thread->readouts);
+    return;
+}
+
+ppStackThread *ppStackThreadAlloc(psArray *readouts)
+{
+    ppStackThread *thread = psAlloc(sizeof(ppStackThread));
+    psMemSetDeallocator(thread, (psFreeFunc)stackThreadFree);
+
+    thread->readouts = psMemIncrRefCounter(readouts);
+    thread->read = false;
+    thread->busy = false;
+    thread->firstScan = 0;
+    thread->lastScan = 0;
+    return thread;
+}
+
+static void stackThreadDataFree(ppStackThreadData *stack)
+{
+    psFree(stack->threads);
+    for (int i = 0; i < stack->imageFits->n; i++) {
+        psFitsClose(stack->imageFits->data[i]);
+        psFitsClose(stack->maskFits->data[i]);
+        psFitsClose(stack->weightFits->data[i]);
+        stack->imageFits->data[i] = stack->maskFits->data[i] = stack->weightFits->data[i] = NULL;
+    }
+    psFree(stack->imageFits);
+    psFree(stack->maskFits);
+    psFree(stack->weightFits);
+    return;
+}
+
+ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames,
+                                          const psArray *maskNames, const psArray *weightNames,
+                                          const pmConfig *config)
+{
+    PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, weightNames, NULL);
+
+    ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
+    psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree);
+
+    stack->lastScan = 0;
+
+    int numInputs = cells->n;           // Number of inputs
+    stack->imageFits  = psArrayAlloc(numInputs);
+    stack->maskFits   = psArrayAlloc(numInputs);
+    stack->weightFits = psArrayAlloc(numInputs);
+    for (int i = 0; i < numInputs; i++) {
+        if (!cells->data[i]) {
+            // Bad image
+            continue;
+        }
+        // Resolved names
+        psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false);
+        psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false);
+        psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false);
+        stack->imageFits->data[i] = psFitsOpen(imageResolved, "r");
+        stack->maskFits->data[i] = psFitsOpen(maskResolved, "r");
+        stack->weightFits->data[i] = psFitsOpen(weightResolved, "r");
+        psFree(imageResolved);
+        psFree(maskResolved);
+        psFree(weightResolved);
+        if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->weightFits->data[i]) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s",
+                    (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]);
+            return false;
+        }
+    }
+
+    int numThreads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
+
+    // Generate readouts for each input file in each file group
+    psArray *threads = psArrayAlloc(numThreads + 1); // A stack for each thread
+    for (int i = 0; i < threads->n; i++) {
+        psArray *readouts = psArrayAlloc(numInputs); // Input readouts
+        for (int j = 0; j < numInputs; j++) {
+            pmCell *cell = cells->data[j]; // Cell with data
+            if (!cell) {
+                continue;
+            }
+            readouts->data[j] = pmReadoutAlloc(cell);
+        }
+        threads->data[i] = ppStackThreadAlloc(readouts);
+        psFree(readouts);               // Drop reference
+    }
+    stack->threads = threads;
+
+    return stack;
+}
+
+
+ppStackThread *ppStackThreadRead(bool *status, ppStackThreadData *stack, pmConfig *config, int numChunk,
+                                 int overlap)
+{
+    assert(status);
+    PS_ASSERT_PTR_NON_NULL(stack, NULL);
+    PS_ASSERT_PTR_NON_NULL(config, NULL);
+
+    *status = true;
+
+    psArray *threads = stack->threads;  // Threads for reading
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe
+    if (!recipe) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE);
+        return NULL;
+    }
+    int rows = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of rows to read per chunk
+    if (rows <= 0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "ROWS is not set in the recipe.");
+        return NULL;
+    }
+
+    // Select an available group
+    while (true) {
+        // check for any groups which can read data
+        for (int j = 0; j < threads->n; j++) {
+            ppStackThread *thread = threads->data[j];
+            if (thread->read) {
+                continue;
+            }
+
+            thread->firstScan = stack->lastScan; // Overlap is taken care of in pmReadoutReadChunk
+            thread->lastScan = stack->lastScan + rows;
+            stack->lastScan = thread->lastScan;
+
+            psArray *readouts = thread->readouts;
+
+            psTimerStart ("ppStackReadChunk");
+
+            psTrace("ppStack", 2, "Reading data for chunk %d into group %d....\n", numChunk, j);
+            for (int i = 0; i < readouts->n; i++) {
+                pmReadout *ro = readouts->data[i]; // Input readout
+                if (!ro) {
+                    continue;
+                }
+
+                // override the recorded last scan
+                ro->thisImageScan  = thread->firstScan;
+                ro->thisWeightScan = thread->firstScan;
+                ro->thisMaskScan   = thread->firstScan;
+                ro->lastImageScan  = thread->lastScan;
+                ro->lastMaskScan   = thread->lastScan;
+                ro->lastWeightScan = thread->lastScan;
+                ro->forceScan      = true;
+
+                psFits *imageFits  = stack->imageFits->data[i]; // FITS file for image
+                psFits *maskFits   = stack->maskFits->data[i]; // FITS file for mask
+                psFits *weightFits = stack->weightFits->data[i]; // FITS file for weight
+
+
+                bool keepReading = false;
+                if (pmReadoutMore(ro, imageFits, 0, rows, config)) {
+                    keepReading = true;
+                    if (!pmReadoutReadChunk(ro, imageFits, 0, rows, overlap, config)) {
+                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT %d",
+                                numChunk, i);
+                        *status = false;
+                        return NULL;
+                    }
+                }
+
+                if (pmReadoutMoreMask(ro, maskFits, 0, rows, config)) {
+                    keepReading = true;
+                    if (!pmReadoutReadChunkMask(ro, maskFits, 0, rows, overlap, config)) {
+                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d",
+                                numChunk, i);
+                        *status = false;
+                        return NULL;
+                    }
+                }
+
+                if (pmReadoutMoreWeight(ro, weightFits, 0, rows, config)) {
+                    keepReading = true;
+                    if (!pmReadoutReadChunkWeight(ro, weightFits, 0, rows, overlap, config)) {
+                        psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.WEIGHT %d",
+                                numChunk, i);
+                        *status = false;
+                        return NULL;
+                    }
+                }
+                if (!keepReading) {
+                    return NULL;
+                }
+            }
+
+            thread->read = thread->busy = true;
+            return thread;
+        }
+
+        // Check for threads that are ready to read
+        bool wait = true;
+        for (int j = 0; j < threads->n; j++) {
+            ppStackThread *thread = threads->data[j];
+            if (thread->busy) {
+                continue;
+            }
+            thread->read = false;
+            wait = false;
+        }
+        if (wait) {
+            // No threads currently available
+            usleep(THREAD_WAIT);
+        }
+    }
+    return NULL;
+}
+
+
+void ppStackThreadInit(void)
+{
+    static bool threaded = false;       // Are we running threaded?
+    if (threaded) {
+        psAbort("Already running threaded.");
+    }
+    threaded = true;
+
+    {
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 7);
+        task->function = &ppStackReadoutInitialThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
+
+    {
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2);
+        task->function = &ppStackInspect;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
+
+    {
+        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 5);
+        task->function = &ppStackReadoutFinalThread;
+        psThreadTaskAdd(task);
+        psFree(task);
+    }
+
+    return;
+}
Index: /tags/pap_merge_20090108/ppStack/src/ppStackVersion.c
===================================================================
--- /tags/pap_merge_20090108/ppStack/src/ppStackVersion.c	(revision 21200)
+++ /tags/pap_merge_20090108/ppStack/src/ppStackVersion.c	(revision 21200)
@@ -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;
+}