Context Navigation

← Previous Change
Next Change →

psLibSDRS.tex

Timestamp:

May 19, 2004, 5:41:50 PM (22 years ago)

Author:

eugene

Message:

extensive edits based on MHPCC comments, see ChangeLogSDRS.tex

File:

: 1 edited

trunk/doc/pslib/psLibSDRS.tex (modified) (86 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/doc/pslib/psLibSDRS.tex

-              r703
+              r747
 %%% $Id: psLibSDRS.tex,v 1.45 2004-05-16 04:50:56 eugene Exp $
+%%% $Id: psLibSDRS.tex,v 1.46 2004-05-20 03:41:49 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 …
 \setlength{\topsep}{-2pt}
 \begin{document}
 \maketitle
+\sloppy
 % -- Revision History --
 …
 DR & 2004 Mar 29 & Draft \\ \hline
 & 2004 Apr 1  & First version, sent to MHPCC \\ \hline
 & 2004 Apr 23 & Added section on error handling \\ \hline
+& 2004 May 19 & Extensive modifications, see appendix \\ \hline
 \RevisionsEnd
 …
 well-defined, concise operations which can be coded with only a modest
 number of lines.  PSLib is a library of basic functions required by
 the IPP, but should include many programming concepts which may be useful
+the IPP, and it includes many programming concepts which may be useful
 for other software projects, especially those which deal with
 astronomical data handling tasks.
 …
 on the functions and data types of the earlier entries.
 The installed base of code for PSLib consists of header files, the
+binary library code, \code{libpslib.a} and the shared-library
+equivalent, \code{libpslib.so} (or \code{libpslib.dylib} in the case
+of OS/X).  Assuming these components have been installed into the
+library and search path, PSLib may be used within a program by
+including the line \code{#include <pslib.h>} into the C code and
 linking with \code{-lpslib}.
+The installed code base for PSLib consists of header files, the binary
+library code, \code{libpslib.a} and the shared-library equivalent,
+\code{libpslib.so} (or \code{libpslib.dylib} in the case of OS/X).
+Assuming these components have been installed into the library and
+search path, PSLib may be used within a program by including the line
+\code{#include <pslib.h>} into the C code and linking with
+\code{-lpslib}.
 This document describes the data structures and details the functions
 calls. The specified data structures and functions follow the naming
 conventions detailed in the IPP Code Conventions (PSDC-430-004).  In
+particular, these coding conventions restrict the namespace used by
+the library functions by requiring that all globally visible symbols
+start with the two letters \code{ps}.  Further namespace organization
+is achieved by encouraging functions to be named in the form
+psNounVerbPhrase, where Noun is the data type of principle relevance
+and VerbPhase describes the operation applied to that data type.  For
+example, the function which copies an image (of type \code{psImage})
 is called \code{psImageCopy()}.
+conventions detailed in the IPP Software Requirements Specification
+(PSDC-430-005).  In particular, these coding conventions restrict the
+namespace used by the library functions by requiring that all globally
+visible symbols start with the two letters \code{ps}.  Further
+namespace organization is achieved by encouraging functions to be
+named in the form psNounVerbPhrase, where Noun is the data type of
+principle relevance and VerbPhase describes the operation applied to
+that data type.  For example, the function which copies an image (of
+type \code{psImage}) is called \code{psImageCopy()}.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 …
 document will be implemented through wrapping external libraries:
 \begin{itemize}
+\item Many of the matrix functions, some of the polynomial and
+minimisation functions functions will wrap the GNU Scientific Library
+(GSL; \href{www.gnu.org/software/gsl/}{\tt www.gnu.org/software/gsl});
+\item The sort functions will wrap the system \code{qsort} call;
+\item Some of the Fourier transform functions will wrap the Fastest Fourier
+Transform in the West library (FFTW; \href{www.fftw.org}{\tt www.fftw.org});
+\item The FITS functions will wrap the CFITSIO library
+(\href{heasarc.gsfc.nasa.gov/docs/software/fitsio/}{\tt
+heasarc.gsfc.nasa.gov/docs/software/fitsio/}); and
+\item Many of the matrix functions, some of the polynomial and some of
+the minimization functions functions should wrap the GNU Scientific
+Library (GSL):
+\href{www.gnu.org/software/gsl/}{\tt www.gnu.org/software/gsl});
+\item The sort functions should wrap the system \code{qsort} call
+\item Some of the Fourier transform functions should wrap the Fastest Fourier
+Transform in the West library (FFTW):
+\href{www.fftw.org}{\tt www.fftw.org}
+\item The FITS functions should wrap the CFITSIO library:
+\href{heasarc.gsfc.nasa.gov/docs/software/fitsio}{\tt heasarc.gsfc.nasa.gov/docs/software/fitsio}
 \item Many of the astronomy routines will wrap the StarLink Positional
+Astronomy libraries (SLALib;
+Astronomy libraries (SLALib):
 \href{star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{\tt
 star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}.
+\item \tbd{Some graphics library, possibly the SM library.}
+star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}
 \end{itemize}
 …
 \subsubsection{Introduction}
+The \PS{} software system will need a level of memory management
+placed between the operating system (\code{malloc}/\code{free}) and
+the high level routines (e.g.\ \code{psMetaDataAlloc}).  This layer is
+in addition to the possibility that specific heavily used data types
+may need their own special-purpose memory managers.  However, since we
+require that all user-level objects be allocated via associated
+\code{Alloc/Free} functions, we will easily be able to implement such
 functionality without impacting the facilities described here.
+PSLib needs a level of memory management placed between the operating
+system (\code{malloc}/\code{free}) and the high level routines (e.g.\
+\code{psMetadataAlloc}).  This layer is in addition to the possibility
+that specific heavily used data types may need their own
+special-purpose memory managers.  However, since we require that all
+user-level objects be allocated via associated \code{Alloc/Free}
+functions, we will easily be able to implement such functionality
+without impacting the facilities described here.
 \subsubsection{Rationale}
 …
   We wish to insulate ourselves from the details of the system-provided
   \code{malloc}.  There is no guarantee that the goals of the system
   architect align with those of the \PS{} processing
+  architect align with those of the PSLib or the IPP.
 \item
 …
   While it is possible to do this by linking with external libraries
   (e.g.\ \href{gnu.org}{Electric Fence}), it is convenient to do so
   within the \PS{} framework.
+  within the framework provided by PSLib.
 \item
 …
 \subsubsection{Memory Management}
 In the following sections, we specify the API set, and define the
 appropriate data structures, needed by the PSLib memory management
+In the following sections, we specify the API set and define the
+appropriate data structures needed by the PSLib memory management
 system in order to meet the requirements specified by the desiderata
 listed above.
+listed above.
 Within the PSLib memory management system, every allocated memory
 …
 memory segments.  The segment preceeding the user-memory contains data
 describing the allocated block, using the \code{psMemBlock} structure.
 The first and last elements of this structure \code{void} pointers
+The first and last elements of this structure are \code{void} pointers
 called \code{startblock} and \code{endblock}, which are assigned a
 special value, \code{PS_MEM_MAGIC}.  The segment following the
 user-memory block consists of a single \code{void} pointer, and is
 also assigned the special value of \code{PS_MEM_MAGIC}.  This address
 is pointed to by the structure elements \code{endpost}.
 In practice, these bounding memory blocks mean that when a user is
+is pointed to by the structure element \code{endpost}.
+In practice, these bounding memory blocks mean that when a user
 requests $N$ bytes of memory, the memory management system in fact
 allocates \code{N + sizeof(psMemBlock) + sizeof(void)} bytes, starting
 at a particular address, \code{ADDR}.  It then fills in the first
 \code{sizeof(psMemBlock)} bytes with the data of the \code{psMemBlock}
 structure, and the last \code{sizeof(void)} bytes with the
+structure, and the last \code{sizeof(void)} bytes with
 \code{PS_MEM_MAGIC}.  It returns to the user the pointer corresponding
 to the address \code{ADDR + sizeof(psMemBlock)}.  If the memory
 management system reallocates a block of memory, it must also allocate
 the additional space and fill in the boundary values.  If the memory
 management system is given a specific pointer for some operation, it is
+able to find the corresponding \code{psMemBlock} by simply subtracting
 \code{sizeof(psMemBlock)} from the pointer address.
+management system is given a specific pointer for some operation, it
+is able to find the corresponding \code{psMemBlock} by simply
+subtracting \code{sizeof(psMemBlock)} from the pointer address.
 The purpose of the three boundary markers is to catch corruption and
 …
 is one where the coder mis-counts the range and either fills the data
 just before the start of the valid memory or just after the end of the
 valid memory.  These actions will (hopefully) alter the boundary-post
+valid memory.  These actions are likely to alter the boundary-post
 values, which can be detected by the memory management system.  In the
 second case, hexadecimal dumps of large blocks of memory are easier to
 …
     const int lineno;                   ///< set from __LINE__ in e.g. p_psAlloc
     int refCounter;                     ///< how many times pointer is referenced
     const void **endpost;               ///< pointer to endpost, initialised to PS_MEM_MAGIC
     const void *endblock;              ///< initialised to PS_MEM_MAGIC
+    const void **endpost;               ///< initialised to PS_MEM_MAGIC
+    const void *endblock;               ///< initialised to PS_MEM_MAGIC
 } psMemBlock;
 \end{verbatim}
 …
 memory.  As discussed below, the basic free function, \code{psFree},
 is specified to free the memory block only if the reference counter is
 set to 1.  See the discussion of the \code{psDlist} and \code{psHash}
+data containers for an example of the usage.  Usage of this feature is
 strongly encouraged, but not enforced by the memory management system.
+set to 1.  See the discussion in section~\ref{sec:free} for an example
+of the usage.  Usage of this feature is strongly encouraged, but not
+enforced by the memory management system.
 In order to trace double frees and other memory errors, the memory
 …
 data are left behind.  If endpost points to the memory location
 immediately following the \code{psMemBlock} data, then the memory
 block has been freed.  This state shall be enforces by \code{psFree}.
+block has been freed.  This state must be enforced by \code{psFree}.
 \subsubsection{APIs for Allocating and Freeing}
+PSLib provides the following APIs to create and destroy memory blocks:
+PSLib must provide the following APIs to create and destroy memory
+blocks:
+%
 \begin{verbatim}
 …
 \code{psRealloc}, and \code{psFree} have identical semantics to the
 standard C library functions \code{malloc}, \code{realloc}, and
 \code{free}.  In fact, these functions shall be implemented as C
+\code{free}.  In fact, these functions should be implemented as C
 preprocessor macros which call the following private functions:
+%
 …
 In order to enforce the use of the PSLib versions, the header file
 shall take steps to ensure that code calling the functions
+must take steps to ensure that code calling the functions
 \code{malloc}, \code{calloc}, \code{realloc}, or \code{free} will not
 compile.  This may be achieved by defining preprocessor macros which
+mask these functions with invalid statements (\eg{} \code{#define malloc(S) for}).
+In exceptional cases, such as the memory management system itself,
+programmers may choose to override this prohibition by defining the
+symbol \code{PS_ALLOW_MALLOC}.  Application code will call
+\code{p_psAlloc,p_psRealloc,p_psFree} via the macros defined above.
+mask these functions with invalid statements (\eg{} \code{#define
+malloc(S) for}).  In exceptional cases, such as the memory management
+system itself, programmers may choose to override this prohibition by
+defining the symbol \code{PS_ALLOW_MALLOC}.  Application code will
+call \code{p_psAlloc}, \code{p_psRealloc}, or \code{p_psFree} via the
+macros defined above.
 The functions \code{psAlloc} and \code{psRealloc} shall never return a
+The functions \code{psAlloc} and \code{psRealloc} must never return a
 \code{NULL} pointer. If they are unable to provide the requested
 memory they should attempt to obtain the desired memory by calling the
+memory they must attempt to obtain the desired memory by calling the
 routine registered by \code{psMemExhaustedSetCallback} (see
 \S\ref{secMemAdvanced}), and if still unsuccessful, call
 …
 types, specified below.  The callbacks are set using functions with
 names of the form \code{psCallbackSet}.  In all cases, the
+`\code{Set}' routine takes a pointer to the desired callback
+function and returns a pointer to the one that was previously
+installed. If the function pointer is \code{NULL}, the default
+callback function is reinstalled.  We discuss the use of each of the
+four callbacks below.
+`\code{Set}' routine takes a pointer to the desired callback function
+and returns a pointer to the one that was previously installed. The
+defaults for each of these callbacks is \code{NULL}, in which case the
+corresponding callback is skipped.  If the function pointer passed to
+the functions above is \code{NULL}, the default behavior is set.  We
+discuss the use of each of the four callbacks below.
 \subsubsubsection{\tt psMemExhaustedCallback}
 …
 information.  No return value is accepted, and no specific operations
 are expected.  The callback is for informational purposes only.  Where
 practical and efficient, the memory manager shall call the routine
+practical and efficient, the memory manager must call the routine
 registered using \code{psMemProblemCallbackSet} whenever a corrupted block
 of memory is discovered.  For example, doubly-freed blocks can be
 …
 typedef long (*psMemFreeCallback)(const psMemBlock *ptr);
 psMemFreeCallback psMemFreeCallbackSet(psMemFreeCallback func);
+long psMemGetId(void);
 \end{verbatim}
+%
 …
 \code{psMemAllocateIDSet} accept the desired ID value and return the
 old value to the user.  The return values of the handlers installed by
 \code{psMemAllocateCallbackSet} and \code{psMemFreeCallbackSet} are used to
 increment the values of \code{p_psMemAllocateID} and
+\code{psMemAllocateCallbackSet} and \code{psMemFreeCallbackSet} are
+used to increment the values of \code{p_psMemAllocateID} and
 \code{p_psMemFreeID} respectively.  For example, a return value of
 \code{0} implies that the value is unchanged; if the value is \code{2}
 the callback will be called again when the memory ID counter has
 increased by two.  This functionality may be useful to check, for
 example, every 100th block allocated.  The function, \code{long psMemGetId(void);}
+example, every 100th block allocated.  The function \code{psMemGetId}
 returns the next identification number to be assigned to a memory
 block.  This function can be used to guide the choice of ID set with
 …
 If the argument \code{array} is non-\code{NULL}, then \code{**array}
 is set to an array of \code{psMemBlock *} pointers when the function
+is set to an array of pointers to \code{psMemBlock} when the function
 returns.  These pointers represent the blocks which have been
 allocated but not freed.  It is the caller's responsibility to free
 …
 The return value is the number of corrupted blocks detected. If the
 argument \code{abort_on_error} is true, \code{psMemCheckCorruption}
 shall call \code{psAbort} as soon as memory corruption is detected.
+must call \code{psAbort} as soon as memory corruption is detected.
 \subsubsection{Reference Counting}
 …
 and decrement it when those references are removed.  The memory
 management routines respect the use of the \code{refCounter} field:
 \code{psFree} will not free a block for which \code{refCounter != 1},
 and \code{psAlloc} will initialize the field to 1.  \code{psFree}
 must generate an error if \code{refCounter != 1}.  However, they do
 not (and cannot practically) enforce the use of the counters; this is
 a requirement of external APIs which intend to use the feature.
+\code{psFree} will not free a block for which \code{refCounter} is not
+, and \code{psAlloc} will initialize the field to 1.  \code{psFree}
+must generate an error if \code{refCounter} is not 1.  However, they
+do not (and cannot practically) enforce the use of the counters; this
+is a requirement of external APIs which intend to use the feature.
 Several APIs are provided to manage the reference counters.  These
 …
 The functions all take a pointer to the start of a user block of
 memory.  The first simply returns the value of the reference counter.
+The next two functions increment or decrement the reference counter,
+returning the pointer which was passed in. These functions must
+validate the memory pointer by determining the corresponding
+If \code{vptr} is \code{NULL}, this function must return a value of
+NULL.  The next two functions increment or decrement the reference
+counter, returning the pointer which was passed in. These functions
+must validate the memory pointer by determining the corresponding
 \code{psMemBlock.id} and checking for consistency in the internal
+memory block table (the table pointer for \code{psMemBlock.id} should
+be in the valid range and should correspond to the address of the
+\code{psMemBlock}).  For an example implementation of the
+\code{refCounter} facilities, see the discussion of \code{psDlist}
+(\S\ref{sec:psDlist}).
+memory block table (the table pointer for \code{psMemBlock.id} must be
+in the valid range and must correspond to the address of the
+\code{psMemBlock}).
 \subsubsection{Relation of Memory Management to Structures}
 \label{sec:free}
+In this document, we specify several C \code{struct}s.  It is expected
+that instances of, for example, \code{struct psMyType} will be
+constructed using \code{psMyTypeAlloc()} calls, and destroyed using
+Within PSLib and throughout the Pan-STARRS project, we specify a
+variety of rich data structures.  The IPP Software Requirements
+Specification states that structures should be defined with
+corresponding constructors and destructors.  Instances of, for
+example, \code{psMyType} should be constructed using
+\code{psMyTypeAlloc()} calls, and destroyed using
 \code{psMyTypeFree()} calls.  The allocator will allocate the required
 memory with \code{psAlloc} and increment the appropriate
 …
 \begin{verbatim}
+void psMyTypeFree(psMyType *myType ///< Object to destroy
+    )
+void psMyTypeFree(psMyType *myType)
+{
     /* No operation if object is NULL */
     if (myType == NULL) {
+    /* data is not defined */
+    if (psMemGetRefCounter(myType) < 1) {
         return;
+    }
     /* Only call psFree if reference counter is 1 */
     if (psMemGetRefCounter(myType) == 1) {
+        psSubFree (myType->sub);
         psFree(myType);
         return;
 …
 \end{verbatim}
+This allows, for example, the \code{psMyType} to be imported into the
+metadata (\S\ref{sec:metadata}) without the user worrying about the
+details of the memory allocation/deallocation:
+\begin{verbatim}
+void psFooMetadata(psMetadata *md)
+{
+    psFoo *foo = psFooAlloc();
+    (void) psMetaDataAppend(md, psMetaDataItemAlloc(PS_META_FOO,foo,"Comment","foo.bar"));
+    (void) psFooFree(foo);
+}
+\end{verbatim}
+In the above case, \code{foo} is created, stuffed into the metadata,
+and then the programmer follows the rule of ``for every \code{alloc},
+there is an equal and opposite free'' before the function returns.
+However, the metadata needs to carry around the \code{psFoo}, and so
+it is important that \code{psFooFree} does not free the memory for
+\code{foo}, but only decrements its \code{refCounter}.  Hence, at the
+conclusion of the function, the memory pointed to by \code{foo} in the
+course of the function remains allocated, and the corresponding
+\code{refCounter} is 1 (specifically, the reference in the metadata).
+Note that the element of \code{myType}, \code{myType.sub} is
+explicitly freed with its associated destructor.  If this element
+points to a data block with multiple references, this call would only
+decrement the counter.
 \subsection{Tracing and Logging}
 …
 top-level.  If the user needs to dig deeper into the code, the trace
 level should be set lower, and the more detailed messages could be
 examined.  In a case of a real, poorly-understood problem with the
+examined.  In a case of a serious, poorly-understood problem with the
 code, the trace threshold would be placed to the bottom and the
 lowest-level step-by-step messages would be printed.
 The PSLib tracing facility will provide the above functionality, along
+The PSLib tracing facility provides the above functionality, along
 with the ability to assign different trace levels to messages from
 different software components.  Each trace message when placed in the
 code is assigned to be part of a specific tracing 'facility', defined
+different software components.  Each trace message, when placed in the
+code, is assigned to be part of a specific tracing 'facility', defined
 in more detail below.  The trace level for that specific message is
 also set when the message is placed.  Each facility may have its trace
 …
 which returns the trace level of the named facility following the
 rules specified above.  A specified trace message (identified by
 \code{psTrace}) shall be printed if and only if
+\code{psTrace}) must be printed if and only if
 \code{psTraceGetLevel(facil)} returns a value greater than or equal to
 the value of \code{myLevel} for that message.  That is, a larger
 …
 hence is more verbose.
+PSLib will include a utility function for examining the current
+tracing levels of all facilities: \code{void psTracePrintLevels(void);}.
+This function will print the hierarchy of trace facilities along with
+the current trace level for each facility.  For example, a particular
+PSLib includes a utility function for examining the current tracing
+levels of all facilities:
+%
+\begin{verbatim}
+void psTracePrintLevels(void);
+\end{verbatim}.
+%
+This function prints the hierarchy of trace facilities along with the
+current trace level for each facility.  For example, a particular
 program may have a few facilities defined, along with their trace
 levels.  A call to \code{psTracePrintLevels} may produce a listing
 …
 flat-field image is foo.fits
   doing the divide
  got an invalid pixel value (NaN) at 500,20
+  got an invalid pixel value (NaN) at 500,20
   divide is done
 \end{verbatim}
 …
+%
 The availability of the tracing facility at run-time, shall be decided
+The availability of the tracing facility at run-time, must be decided
 at compilation: If the C pre-processor macro \code{PS_NO_TRACE} is
 defined, all trace code shall be replaced by empty space so that none
+defined, all trace code must be replaced by empty space so that none
 of the code is compiled.  This can be implemented via macro front-ends
 to private versions of the user APIs.  In addition, a function
 …
 with \code{psTraceSetDestination}:
 \begin{verbatim}
+void psTraceSetDestination(FILE *fp     // Open file pointer to write to
+    );
+void psTraceSetDestination(FILE *fp);
 \end{verbatim}
 If the \code{fp} is \code{NULL}, then the trace is sent to standard
 output.
+output, otherwise it is sent to the specified file pointer.
 \subsubsection{Message Logging}
 …
 \begin{verbatim}
 void psLogMsg(char *name, int myLevel, char *fmt, ...);
 void psLogVMsg(char *name, int myLevel, char *fmt, va_list ap);
+void psLogMsgV(char *name, int myLevel, char *fmt, va_list ap);
 \end{verbatim}
 where \code{name} is a word to describe the source of the message,
 …
 is a printf-style formatting statement defining the actual message,
 and is followed by the arguments to the formatting code.  The second
 form, \code{psLogVMsg} is an equivalent command which takes a
+form, \code{psLogMsgV} is an equivalent command which takes a
 \code{va_list} argument.
 …
 \end{verbatim}
+%
-The default level is set to \code{PS_LOG_INFO}.
 At any time, the program may set the current log level, the level
 …
 invoked with \code{psLogMsg} is only printed if its value of
 \code{myLevel} is less than the current value set by
+\code{psLogSetLevel}.
+\code{psLogSetLevel}.  The default log level is set to
+\code{PS_LOG_INFO}.
 Log messages are sent to the destination most recently set using:
+%
 \begin{verbatim}
+int psLogSetDestination(int dest);
+\end{verbatim}
+%
+The only values that are initially defined are \code{PS_LOG_TO_STDERR}
+and \code{PS_LOG_TO_STDOUT} to write to \code{stderr} and
+\code{stdout} respectively, and \code{PS_LOG_TO_NONE} to turn off
+logging.  \tbd{Log to a server via TCP/IP; will possibly set these
+values to negative integers and allow psLogSetDestination to take
+a positive-integer file descriptor.}.
+int psLogSetDestination(char *dest);
+\end{verbatim}
+%
+The destination string consists of a URL in the form
+\code{protocol:location}.  The \code{protocol} value may be
+\code{file}, to send the log to a local file named by the value of
+\code{location}.  Future expansion may allow the logger to send
+messages to an IP logger, with a protocol to be defined later.  Three
+other special values are allowed for \code{dest}: \code{stderr} and
+\code{stdout}, which write to \code{stderr} and \code{stdout}
+respectively, and \code{none} to turn off logging.
 The output format is controlled with the function:
 …
 \code{PS_LOG_ERROR}, \code{PS_LOG_WARN}, and \code{PS_LOG_INFO}
 respectively. Other levels are represented numerically (\code{5}
 etc.). The other two field, \code{name} and \code{msg} are the
+etc.). The other two fields, \code{name} and \code{msg}, are the
 arguments to \code{psLogMsg}; note that \code{name} has a fixed width
 of 15 characters. If \code{msg} doesn't end in a newline, a single
 …
 \hlabel{errorStack}
+\PS{} errors shall be raised using a function, \code{psError}, with
+the caller supplying a component name and error message.  It is
+desireable to be able to trace an error through the program so that
+the events that led to the error may be quickly and clearly
+identified.
+\begin{verbatim}
+/// Prints an error message and doesn't abort; returns code
+int psError(const char *name,           ///< Category of code that caused the error
+            psErrorCode code,           ///< class of error (equivalent to errno)
+            psErrorStatus status,       ///< is this a new error?
+            const char *fmt,            ///< Format
+            ...                         ///< Extra arguments to use format
+    );
+typedef enum {
+    PS_OLD_ERROR = 0,                   ///< This is an old error, and should append to the error stack
+    PS_NEW_ERROR = 1,                   ///< This is a new error and should clear the error stack
+} psErrorStatus;
+\PS{} errors must be raised using a function, \code{psError}, with the
+caller supplying a component name and error message.  It is desireable
+to be able to trace an error through the program so that the events
+that led to the error may be quickly and clearly identified.
+\code{psError} prints an error message and doesn't abort, but instead
+returns the error code.
+\begin{verbatim}
+int psError(const char *name, psErrorCode code, int status, const char *fmt, ...);
 \end{verbatim}
 The \code{name} is of the form \code{aaa.bbb.ccc} and identifies the
 component raising the error.  The \code{psErrorCode} is an enumerated
 type which lists the possible \textit{classes} of errors
+component raising the error.  The \code{code} is an enumerated type
+which lists the possible \textit{classes} of errors
 (e.g. \code{PS_ERR_IO}) that \PS{} code can generate (see section
+\ref{psErrorCodes}). \code{status} specifies whether this is a new
+error, or whether this call to \code{psError} is in response to an
+error that has already resulted in a call to \code{psError}.  The
+\ref{psErrorCodes}). The \code{new} argument takes a boolean which, if
+\code{TRUE} specifies that the error was set initially at this
+location, and if \code{FALSE} specifies that an error was passed to
+this location.  Raising new error should clear the error stack.  The
 final required argument, \code{fmt}, is a \code{printf}-style format
 that is passed to \code{psLogVMsg} with code \code{PS_LOG_ERROR}.
+The result of a call to \code{psError} shall be to push an error onto
+a stack; this stack is cleared if \code{psErrorStatus} is true, or by a call
+to \code{psErrorClear}.  The errors are defined as the following:
+that is passed to \code{psLogMsgV} with code \code{PS_LOG_ERROR}.  The
+result of a call to \code{psError} must be to push an error onto a
+stack; this stack is cleared if \code{psErrorStatus} is true, or by a
+call to \code{psErrorClear}.
+The errors on the error stack are defined as the following:
 \begin{verbatim}
 typedef struct {
 …
 \end{verbatim}
+The last error reported is available from \code{psLastError}; if no
+errors are current, a non-\code{NULL} \code{psErr} shall be returned
+with code \code{PS_ERR_NONE}.  Previous errors on the stack shall be
+returned by \code{psGetError} (a value of \code{0} passed to
+\code{psGetError} is equivalent to a call to \code{psLastError}).
+\begin{verbatim}
+const psErr *psErrorGet(int which);     ///< return specified error (or an "error" with code PS_ERR_NONE)
+const psErr *psErrorLast(void);         ///< return last error (or an "error" with code PS_ERR_NONE)
+\end{verbatim}
+The error stack may be completely cleared:
+\begin{verbatim}
+void psErrorClear(void);                ///< Clear the error stack
+The last error reported is available from \code{psErrorLast}; if no
+errors are current, a non-\code{NULL} \code{psErr} must be returned
+with code \code{PS_ERR_NONE}.  Previous errors on the stack must be
+returned by \code{psErrorGet} (a value of \code{0} passed to
+\code{psErrorGet} is equivalent to a call to \code{psErrorLast}).
+The error stack may be completely cleared with \code{psErrorClear}.
+%
+\begin{verbatim}
+const psErr *psErrorGet(int which);
+const psErr *psErrorLast(void);
+void psErrorClear(void);
 \end{verbatim}
 The entire error stack may be printed to an open file descriptor by
 calling \code{psErrorStackPrint} (or \code{psErrorVStackPrint}); if
+calling \code{psErrorStackPrint} (or \code{psErrorStackPrintV}); if
 and only if there are current errors, the printf-style string
 \code{fmt} is first printed to the file descriptor \code{fd}. In this
 printout, error codes shall be replaced by their string equivalents as
+printout, error codes must be replaced by their string equivalents as
 defined in the next section.  Note that these are also available in
 the \code{psErr} structure. The successive lines of the traceback
 should be indented by an additional space (see example).
+\code{psErrorVStackPrint} shall not invoke \code{va_end}.
+\begin{verbatim}
+void psErrorStackPrint(FILE *fd, const char *fmt, ...); ///< print the errorstack to this file descriptor
+void psErrorVStackPrint(FILE *fd, const char *fmt, va_list va); ///< print the errorstack to this file
+                                                                ///< descriptor
+\code{psErrorStackPrintV} must not invoke \code{va_end}.
+%
+\begin{verbatim}
+void psErrorStackPrint(FILE *fd, const char *fmt, ...);
+void psErrorStackPrintV(FILE *fd, const char *fmt, va_list va);
 \end{verbatim}
 Any \code{errorCode}s less then or equal to \code{PS_ERR_BASE} (see
 next section) shall be taken to be valid values of \code{errno}, and
 \code{psErrorStackPrint} shall print the value returned by
+next section) must be taken to be valid values of \code{errno}, and
+\code{psErrorStackPrint} must print the value returned by
 \code{strerror} if such error codes are encountered.
 The routine \code{psErrorCodeString} returns the string associated
 with an error code:
+\begin{verbatim}
+const char *psErrorCodeString(psErrorCode code);        ///< return the string associated with an error code.
+\end{verbatim}
+\begin{verbatim}
+const char *psErrorCodeString(psErrorCode code);
+\end{verbatim}
 \subsubsection{Error Codes}
 …
 Both error codes for PSLib and error codes for projects that use PSLib
 may be registered.  In the former case, the error codes shall be
 registered on initialisation, whereas in the latter case, it is
 required to explicitly register.
+may be registered.  In the former case, the error codes must be
+registered on initialisation, whereas in the latter case, they must be
+explicitly registered by the programmer.
 \paragraph{Registering error codes}
+An array of error codes may be registered with the PSLib error handler
+using a private function:
+\begin{verbatim}
+void p_psErrorRegister(const psErrorDescription *errors, ///< register a set of errors
+                       int nerror                        ///< number of errors
+    );
+\end{verbatim}
+Where the errors are represented internally as follows:
+PSLib and any project needed to use PSLib must define the necessary
+error codes and associated message strings.  An array of error codes
+may be registered with the PSLib error handler using a private
+function:
+\begin{verbatim}
+void psErrorRegisterSet(const psErrorDescription *errors, int nerror);
+\end{verbatim}
+where the errors are represented internally as follows:
 \begin{verbatim}
 typedef struct {
 …
 } psErrorDescription;
 \end{verbatim}
+Projects wishing to employ the PSLib error handler should define
+a function,
+\begin{verbatim}
+void prefixErrorRegister(void);
+\end{verbatim}
+where \code{prefix} is particular to the project.  For example, PSLib
+should have a function,
+\begin{verbatim}
+void psErrorRegister(void);
+\end{verbatim}
+In this case (i.e., error codes for PSLib), the function should be
+called upon initialisation.  As another example, MOPS should implement
+a function,
+\begin{verbatim}
+void mopsErrorRegister(void);
+\end{verbatim}
+which would be called explicitly when MOPS utilises PSLib.
+It is left to the external project to produce the appropriate
+\code{prefixErrorRegister()}, but they may find helpful the
+implementation discussed below for use in PSLib.
+There is a clear need to coordinate the choice of error numbers.  It
+is expected that error code ranges for different projects shall be
+managed by the Project Office.
+PSLib internal errors must be registered with the function
+psErrorRegister, which should be called as part of the program
+initialization to set up the error codes.  It is left to the external
+project to produce their own error registration functions which must
+also be called during initialization. There is a clear need to
+coordinate the choice of error numbers.  It is expected that error
+code ranges for different projects must be managed by the Project
+Office within Pan-STARRS.
 \paragraph{Error Codes for PSLib}
 For ease of maintenance, error codes for PSLib shall be defined by an
+For ease of maintenance, error codes for PSLib must be defined by an
 auxiliary file, conventionally named \file{psErrorCodes.dat}.  The
 format of this file shall consist of a number of lines, each of the
+format of this file must consist of a number of lines, each of the
 form:
 \begin{verbatim}
 …
 where \code{[ = value]} and the comma are optional, and no spaces are
 significant except in the STRING.  Comments extend from \code{#} to
 the end of the line (except that a \code{\#} shall be replaced by
+the end of the line (except that a \code{\#} must be replaced by
 \code{#} and not taken to start a comment). For example,
 \begin{verbatim}
 …
 The values \code{NONE = 0} and {UNKNOWN} must be present.
 A script, called from the Makefiles, shall generate two files,
+A script, called from the Makefiles, must generate two files,
 \file{psErrorCodes.h} and \file{psErrorCodes.c} from the input file
 \file{psErrorCodes.dat}.  \file{psErrorCodes.h} shall define an
+\file{psErrorCodes.dat}.  \file{psErrorCodes.h} must define an
 enumerated type \code{psErrorCode} with elements \code{PS_ERR_NAME}
 and values as specified in \file{psErrorCodes.dat}, e.g.
 …
 \end{verbatim}
+Any \code{errorCode}s less then or equal to \code{PS_ERR_BASE} shall be taken
+to be valid values of \code{errno}.
+\file{psErrorCodes.c} shall define the necessary function to register
+the error codes.
+\subsubsection{Example}
+The following example code:
+\begin{verbatim}
+#include "psLib.h"
+static int primary(int i)
+{
+    if (i != 0) {                       // let's pretend it's an I/O error
+        return psError("tst.error.primary", PS_ERR_IO, 1, "Primary error");
+    }
+    return 0;
+}
+static int middle(void)
+{
+    if (primary(1) != 0) {
+        return psError("tst.error.middle", PS_ERR_UNKNOWN, 0, "Secondary error");
+    }
+    return 0;
+}
+static int toplevel(void)
+{
+    if (middle() != 0) {
+        return psError("tst.error", PS_ERR_UNKNOWN, 0, "Toplevel error");
+    }
+    return 0;
+}
+int main(void)
+{
+    if (toplevel() != 0) {
+        psErrorStackPrint(stdout, "Traceback:\n");
+        if (psErrorLast()->code == PS_ERR_UNKNOWN) {
+            fprintf(stderr, "Last error is of unknown type\n");
+        }
+        if (psErrorGet(2)->code == PS_ERR_IO) {
+            fprintf(stderr, "Third oldest error is of type IO\n");
+        }
+    }
+    psErrorClear();
+    psErrorStackPrint(stdout, "Traceback:\n");
+    if (psErrorLast()->code == PS_ERR_NONE) {
+        fprintf(stderr, "No errors. Hurrah\n");
+    }
+    return 0;
+}
+\end{verbatim}
+Produces the following output:
+\begin{verbatim}
+Traceback:
+tst.error.primary              I/O error                      Primary error
+ tst.error.middle              unknown error                  Secondary error
+  tst.error                    unknown error                  Toplevel error
+Last error is of unknown type
+Third oldest error is of type IO
+No errors. Hurrah
+\end{verbatim}
+Any \code{errorCode}s less then or equal to \code{PS_ERR_BASE} must be
+taken to be valid values of \code{errno}.  \file{psErrorCodes.c} must
+define the necessary function to register the error codes.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 …
 \subsection{Abort}
 \code{psAbort}, shall call \code{psMsgLog} with a
+\code{psAbort}, must call \code{psMsgLog} with a
 level of \code{PS_LOG_ABORT}, and then call \code{abort}.
 …
 Throughout PSLib, we require a variety of structures which correspond
 to different mathematical data concepts.  For example, we have several
+data structures which correspond to one-dimensional arrays (vectors)
+of different data types (\code{int}, \code{float}, etc).  We also have
 different data structures which correspond to two-dimensional arrays
+(images or matrices), again with different data types for the
+individual elements.
+to different mathematical data concepts.  For example, we have a data
+structure which corresponds to one-dimensional arrays (vectors) of
+different data types (\code{int}, \code{float}, etc).  We also have a
+data structure which corresponds to two-dimensional arrays (images or
+matrices), again with different data types for the individual
+elements.
 A variety of functions perform operations which are equivalent for
 …
 operation if $x$ is a vector and $y$ is a matrix.  Nor does it matter
 mathematically that the element data types match; the sum of a float
 and an integer is a well-defined quantity!  One constraint should be
+and an integer is a well-defined quantity.  One constraint should be
 noted: the size of the elements in each dimension must match.  For
 example, if $x$ were a vector of 100 elements, but $y$ were a vector
 …
 Given that some functions should be able to operate equivalently (or
+identically) on a wide range of data types, it seems cumbersome to be
+forced into defining a large number of C functions to handle the
+different data types, just because we have different structures.
+Admittedly, some details of the function would have to vary for
+different data types, but since the basic function is the same, it
+would help both the user and programmer if the same function could be
+used for different data types.  We therefore define a mechanism which
+allows the C functions to accept a generic data type, and determine
+the type of the data on the basis of the data.  The mechanism uses the
+structure \code{psType}.
+Each of these equivalent data types is defined by a structure in which
+the first element is always of type \code{psType}.  This element
+defines both the dimensions of the array and the data type of each
+element.  The structure is as follows:
+identically) on a wide range of data types, we define a mechanism
+which allows the C functions to accept a generic data type, and
+determine the type of the data on the basis of the data.
+Supported data types must be defined by a structure in which
+the first element is always of type \code{psType}:
 \begin{verbatim}
 typedef struct {
 …
 where \code{psDimen dimen} defines the dimensionality of the data and
 \code{psElemType type} defines the data type of each element.  These
 two variable types are defined as structures:
+two variable types are defined as:
 \begin{verbatim}
 typedef enum {
 …
     PS_DIMEN_TRANSV,                    ///< A transposed vector
     PS_DIMEN_IMAGE,                     ///< An image (matrix)
     PS_DIMEN_OTHER                      ///< Something else that's not supported for arithmetic
+    PS_DIMEN_OTHER                      ///< Not supported for arithmetic
 } psDimen;
 \end{verbatim}
 …
     PS_TYPE_F32,                        ///< Floating point
     PS_TYPE_F64,                        ///< Double-precision floating point
     PS_TYPE_C32,                        ///< Complex numbers consisting of floating point
     PS_TYPE_OTHER,                      ///< Something else that's not supported for arithmetic
+    PS_TYPE_C32,                        ///< Complex numbers consisting of floats
+    PS_TYPE_OTHER,                      ///< Not supported for arithmetic
 } psElemType;
 \end{verbatim}
 …
     psType type;                        ///< vector data type and dimension
     const int n;                        ///< size of vector
     const int nalloc;                   ///< data region relative to parent
+    const int nalloc;                   ///< allocated data block
     union {
+        psF32 *arr;                     ///< Pointers to floating-point data (default)
+        psS8  *arr_S8;                  ///< Pointers to short-integer data
+        psS16 *arr_S16;                 ///< Pointers to short-integer data
+        psS32 *arr_S32;                 ///< Pointers to integer data
+        psS64 *arr_S64;                 ///< Pointers to long-integer data
+        psU8  *arr_U18;                 ///< Pointers to unsigned-short-integer data
+        psU16 *arr_U16;                 ///< Pointers to unsigned-short-integer data
+        psU32 *arr_U32;                 ///< Pointers to unsigned-integer data
+        psU64 *arr_U64;                 ///< Pointers to unsigned-long-integer data
+        psF32 *arr_F32;                 ///< Pointers to floating-point data
+        psF64 *arr_F64;                 ///< Pointers to double-precision data
+        psF32 *arr_C32;                 ///< Pointers to complex floating-point data
+        void **arr_v;
+    } arr;
+        psS8  *S8;                      ///< Pointers to short-integer data
+        psS16 *S16;                     ///< Pointers to short-integer data
+        psS32 *S32;                     ///< Pointers to integer data
+        psS64 *S64;                     ///< Pointers to long-integer data
+        psU8  *U18;                     ///< Pointers to unsigned-short-integer data
+        psU16 *U16;                     ///< Pointers to unsigned-short-integer data
+        psU32 *U32;                     ///< Pointers to unsigned-integer data
+        psU64 *U64;                     ///< Pointers to unsigned-long-integer data
+        psF32 *F32;                     ///< Pointers to floating-point data
+        psF64 *F64;                     ///< Pointers to double-precision data
+        psF32 *C32;                     ///< Pointers to complex floating-point data
+        void **void;
+    } data;
 } psVector;
 \end{verbatim}
 …
 (the number of elements); \code{nalloc} is the number of elements
 allocated ($nalloc \ge n$).  The allocated memory is available in the
 union \code{arr} which consists of pointers to each of the defined
+union \code{data} which consists of pointers to each of the defined
 primitive data types.  Note the parallelism in the names of the types,
 union elements, and the psElemType names.  This parallelism allows us
 to use automatic construction mechanisms effectively.  The data type
 is defined by the first element, \code{psType}.  The
 structure is associated with a constructor and a destructor:
+is defined by the first element, \code{psType}.  The structure is
+associated with a constructor and a destructor:
+%
 \begin{verbatim}
 psVector *psVectorAlloc(int nalloc, psElemType type);
 psVector *psVectorRealloc(const psVector *vector, int nalloc);
 void psVectorFree(psVector *restrict vector);
+void psVectorFree(psVector *restrict vector, void (*elemFree)(void *));
 \end{verbatim}
+%
 …
 \code{nalloc} is larger than the current value of \code{psVector.n},
 \code{psVector.n} is left intact.  If the value of \code{myArray} is
+\code{NULL}, then \code{psVectorRealloc} must return an error.
+\subsection{Arrays of Pointer Types}
+Arrays of pointer types need some additional specification.  We
+require an array of pointers of type \code{void}, with which we can
+carry around a collection of data of an arbitrary type which is more
+complicated than the simple numeric types above.  The structure is as
+follows:
+%
+\begin{verbatim}
+typedef struct {
+    psType type;                        ///< Type of data.  Must be first element
+    int nAlloc;                         ///< Total number of elements available
+    int n;                              ///< Number of elements in use
+    void **arr;                         ///< The array data
+} psVoidPtrArray;
+\end{verbatim}
+%
+There is also an equivalent set of constructors and destructor:
+%
+\begin{verbatim}
+psVoidPtrArray *psVoidPtrArrayAlloc(int nAlloc);
+psVoidPtrArray *psVoidPtrArrayRealloc(psVoidPtrArray *myArray, int nAlloc);
+void psVoidPtrArrayFree(psVoidPtrArray *restrict myArray, void (*elemFree)(void *));
+\end{verbatim}
+%
+The only difference with the numeric array types is the addition of a
+destructor function which is passed to \code{psVoidPrtArrayFree}.
+This function, which may be \code{NULL}, is called for each existing
+element of the array before the array itself is freed.  If the
+function is \code{NULL}, the elements are are not freed.
+The routine \code{psVoidPtrArrayFree} assumes that all pointers had
+their reference counters incremented when they were inserted onto the
+array.\footnote{\eg{} \code{va->arr[i] = psMemIncrRefCounter(ptr);}}
+If \code{psVoidPtrArrayFree}'s argument \code{elemFree} is NULL, the
+list should be deleted, but not the elements on it (an error should be
+raised if the \code{refCounter} is 1; otherwise their
+\code{refCounter}'s should be decremented) --- this is to account for
+an array of heterogeneous types.
+\code{NULL}, then \code{psVectorRealloc} must return an error.  In
+\code{psVectorFree}, the function \code{elemFree} is required for
+arrays of pointer types; it is the destructor appropriate to the data
+pointed to by the pointers.  This function, which may be \code{NULL},
+is called for each existing element of the array before the array
+itself is freed.  If the function is \code{NULL}, the elements are are
+not freed.  This function must not be defined for any data type except
+the \code{void} pointer array.
 \subsection{Simple Images}
 …
     const int x0, y0;                   ///< data region relative to parent
     union {
+        psF32 **rows;                   ///< Pointers to floating-point data (default)
+        psS8  **rows_S8;                ///< Pointers to char data
+        psS16 **rows_S16;               ///< Pointers to short-integer data
+        psS32 **rows_S32;               ///< Pointers to integer data
+        psS64 **rows_S64;               ///< Pointers to long-integer data
+        psU8  **rows_U8;                ///< Pointers to unsigned-char data
+        psU16 **rows_U16;               ///< Pointers to unsigned-short-integer data
+        psU32 **rows_U32;               ///< Pointers to unsigned-integer data
+        psU64 **rows_U64;               ///< Pointers to unsigned-long-integer data
+        psF32 **rows_F32;               ///< Pointers to floating-point data
+        psF64 **rows_F64;               ///< Pointers to double-precision data
+        psC32 **rows_C32;               ///< Pointers to complex floating-point data
+    } rows;
+        psS8  **S8;                     ///< Pointers to char data
+        psS16 **S16;                    ///< Pointers to short-integer data
+        psS32 **S32;                    ///< Pointers to integer data
+        psS64 **S64;                    ///< Pointers to long-integer data
+        psU8  **U8;                     ///< Pointers to unsigned-char data
+        psU16 **U16;                    ///< Pointers to unsigned-short-integer data
+        psU32 **U32;                    ///< Pointers to unsigned-integer data
+        psU64 **U64;                    ///< Pointers to unsigned-long-integer data
+        psF32 **F32;                    ///< Pointers to floating-point data
+        psF64 **F64;                    ///< Pointers to double-precision data
+        psC32 **C32;                    ///< Pointers to complex floating-point data
+    } data;
     const struct psImage *parent;       ///< parent, if a subimage
     int Nchildren;                      ///< number of subimages
     struct psImage **children;          ///< children of this region; array of Nchildren pointers
+    struct psImage **children;          ///< children of this region (Nchildren total)
 } psImage;
 \end{verbatim}
 …
 (\code{parent}).  Unless this is image is a child of another image
 (represents a subset of the pixels of another image), the image data
+is allocated in a contiguous block.
+Create an image of a specified width, height, and data type.  This
+function must allow any of the valid image data types and not restrict
+to the valid FITS BITPIX types.
+is allocated in a contiguous block.  We define the following
+supporting functions:
 \begin{verbatim}
 psImage *psImageAlloc (int width, int height, psElemType type);
 \end{verbatim}
+where \code{width} and \code{height} specify the size of the image and
+\code{type} specifies the data type and the image dimensionality
+(which must be 2).
+Create an image of a specified \code{width}, \code{height}, and data
+\code{type}.  This function must allow any of the valid image data
+types and not restrict to the valid FITS BITPIX types.  The image
+dimensionality must be 2.
+\begin{verbatim}
+void psImageFree(psImage *image);
+\end{verbatim}
 Free the memory associated with a specific image, including the pixel
 data. Free the children of the image if they exist.
-\begin{verbatim}
-void psImageFree(psImage *image);
-\end{verbatim}
-Free only the pixels for a specified image:
-\begin{verbatim}
-psImage *psImageFreePixels(psImage *image);
-\end{verbatim}
 \subsection{Doubly-linked lists}
 …
 The type \code{psDlist} represents the container of the list.  It has
 a pointer to the first element in the linked list (\code{head}), a
 pointer to the last element in the list (\code{tail}), \tbd{an entry
 for the current cursor location (\code{iter})}, and an entry to define
 the number of elements in the list (\code{n}).
+pointer to the last element in the list (\code{tail}), an entry for
+the current cursor location (\code{iter}), and an entry to define the
+number of elements in the list (\code{n}).
 The elements of the list are defined by the type \code{psDlistElem}:
 …
 (\code{next}), the previous element in the list (\code{prev}), and a
 \code{void} pointer to whatever data is represented by this list
+element.
+A list may be created with the function
+element.    The following supporting functions are required:
 \begin{verbatim}
 psDlist *psDlistAlloc(void *data);
 \end{verbatim}
+which may take a pointer to a data item, or it may take \code{NULL}.
+The allocator creates both the \code{psDlist} and the first element in
+the list, pointed to by both \code{psDlist.head} and
 \code{psDlist.tail}.  If the data entry is \code{NULL}, then an empty
+list, with both pointers are set to \code{NULL} should be created.
+An entry may be added to the list with the function:
+Create a list.  This function may take a pointer to a data item, or it
+may take \code{NULL}.  The allocator creates both the \code{psDlist}
+and the first element in the list, pointed to by both
+\code{psDlist.head} and \code{psDlist.tail}.  If the data entry is
+\code{NULL}, then an empty list, with both pointers are set to
+\code{NULL} should be created.
 \begin{verbatim}
 psDlist *psDlistAdd(psDlist *list, void *data, int where);
 \end{verbatim}
+which takes a pointer to the list and also returns a pointer to the
+list.  The returned pointer must be used as the value of
+\code{psDlist} may have changed.  The value of \code{where} specifies
+if the specified data item should be placed on the front of the list
+(\code{PS_DLIST_HEAD}), at the end of the list (\code{PS_DLIST_TAIL}),
 to add after (\code{PS_DLIST_NEXT}) or before (\code{PS_DLIST_PREV})
+the current element (specified by the iteration cursor), or an index
+that the new \code{data} should inhabit.
+A data item may be retrieved from the list with the function:
+Add an entry to the list with this function, which takes a pointer to
+the list and also returns a pointer to the list.  The returned pointer
+must be used as the value of \code{psDlist} may have changed.  The
+value of \code{where} specifies if the specified data item should be
+placed on the front of the list (\code{PS_DLIST_HEAD}), at the end of
+the list (\code{PS_DLIST_TAIL}), to add after (\code{PS_DLIST_NEXT})
+or before (\code{PS_DLIST_PREV}) the current element (specified by the
+iteration cursor), or an index that the new \code{data} should
+inhabit.
 \begin{verbatim}
 void *psDlistGet(psDlist *list, int which);
 \end{verbatim}
+The value of \code{which} may be the numerical index or it may be one
+of the special values: \code{PS_DLIST_HEAD}, \code{PS_DLIST_TAIL},
+A data item may be retrieved from the list with this function.  The
+value of \code{which} may be the numerical index or it may be one of
+the special values: \code{PS_DLIST_HEAD}, \code{PS_DLIST_TAIL},
 \code{PS_DLIST_PREV}, and \code{PS_DLIST_NEXT}, all of which are
 defined as negative integers, allowing \code{where} to also be the
 index of one of the data items.
-A data item may be removed from the list with the function:
 \begin{verbatim}
 void *psDlistRemove(psDlist *list, void *data, int which);
 \end{verbatim}
+The value of \code{which} may be the numerical index or it may be one
+of the special values: \code{PS_DLIST_HEAD}, \code{PS_DLIST_TAIL},
+A data item may be removed from the list with this function.  The
+value of \code{which} may be the numerical index or it may be one of
+the special values: \code{PS_DLIST_HEAD}, \code{PS_DLIST_TAIL},
 \code{PS_DLIST_PREV}, \code{PS_DLIST_UNKNOWN}, and
 \code{PS_DLIST_NEXT}, all of which are defined as negative integers.
 If the value of \code{which} is \code{PS_DLIST_UNKNOWN}, then the data
+item is identified by matching the pointer value with \code{void *data}.
+All data items placed onto lists (\code{psDlistAdd}) shall have their
+item is identified by matching the pointer value with \code{void
+*data}.
+All data items placed onto lists (\code{psDlistAdd}) must have their
 reference counters (section \ref{secMemRefcounter}) incremented.  When
 elements are removed from a list with \code{psDlistRemove}, they shall
+elements are removed from a list with \code{psDlistRemove}, they must
 have their reference counters decremented. The action of retrieving
 data from a list (with \code{psDlistGet}) shall not affect their
+data from a list (with \code{psDlistGet}) must not affect their
 reference counter.
-A complete list may be freed with the destructor:
 \begin{verbatim}
 void psDlistFree(psDlist *list, void (*elemFree)(void *));
 \end{verbatim}
+If the element destructor (\code{elemFree}) is \code{NULL}, the list
+should be deleted, but not the elements, although their
+\code{refcounter}s should be decremented.
+Two functions are available to convert between the \code{psDlist} and
+\code{psVoidPtrArray} containers:
+\begin{verbatim}
+psVoidPtrArray *psDlistToArray(psDlist *dlist);
+psDlist *psArrayToDlist(psVoidPtrArray *arr);
+\end{verbatim}
+These functions do not free the elements or destroy the input
+collection.  Rather, they increment the reference counter for each of
+the elements.
+Iteration over all elements of the list using the iteration cursor
+\code{iter} is provided by the functions:
+A complete list may be freed with this destructor.  If the element
+destructor (\code{elemFree}) is \code{NULL}, the list should be
+deleted, but not the elements, although their \code{refcounter}s
+should be decremented.
+\begin{verbatim}
+psVector *psDlistToVector(psDlist *dlist);
+psDlist *psVectorToDlist(psVector *vector);
+\end{verbatim}
+These two functions are available to convert between the
+\code{psDlist} and \code{psVector} containers.  These functions do not
+free the elements or destroy the input collection.  Rather, they
+increment the reference counter for each of the elements.
 \begin{verbatim}
 void psDlistSetIterator(psDlist *list, int where, int which);
 …
 void *psDlistGetPrev(psDlist *list, int which);
 \end{verbatim}
+The first of these functions uses the value of \code{where} to set the
+iteration cursor for the given list to the beginning
+\code{PS_DLIST_HEAD} or the end \code{PS_DLIST_TAIL} for the iterator
+specified by \code{which}.  The next two functions move the iteration
+cursor forward or backwards, returning the data item from the
+resulting list entry, or returning \code{NULL} at the end of the list.
+Explicit traversal of the list using the \code{psDlistElem}s
+\code{prev} and \code{next} pointers is also supported.
+Iteration over all elements of the list using the iteration cursor
+\code{iter} is provided by these functions.  The first of these
+functions uses the value of \code{where} to set the iteration cursor
+for the given list to the beginning \code{PS_DLIST_HEAD} or the end
+\code{PS_DLIST_TAIL} for the iterator specified by \code{which}.  The
+next two functions move the iteration cursor forward or backwards,
+returning the data item from the resulting list entry, or returning
+\code{NULL} at the end of the list.  Explicit traversal of the list
+using the \code{psDlistElem}s \code{prev} and \code{next} pointers is
+also supported.
+\begin{verbatim}
+psDlist *psDlistSort(psDlist *list, int (*compare)(const void *a, const void *b) );
+\end{verbatim}
 A list may be sorted using \code{psDlistSort}, which requires the
 specification of a comparison function to specify how the objects on
 the list should be sorted.  The motivation is primarily to be able to
 iterate on a sorted list of keys from a hash.
-\begin{verbatim}
-psDlist *psDlistSort(psDlist *list, int (*compare)(const void *a, const void *b) );
-\end{verbatim}
 \subsection{Hash Tables}
 …
 \section{Data manipulation}
+There are a number of data concepts which can be naturally represented
+in C as structures.  We require a variety of basic data manipulation
+functions which will act upon data (in particular, arrays/vectors).
+We require the following capabilities:
+We require a variety of basic data manipulation functions which will
+act upon data (in particular, arrays/vectors).  We require the
+following capabilities:
 \begin{itemize}
 \item Bit masks;
 …
 \subsection{Bitsets}
 Bitsets are required in order to turn options on and off.  We
+require the capability to have a bitset of arbitrary length (i.e.,
 not limited by the length of a \code{long}, say).  The
+\code{psBitset} structure is defined below.  Note that the entry
 \code{bits} is an array of type \code{char} storing the bits as bits
 of each byte in the array, with 8 bits available for each byte in the
+array.  Also note that the constructor is passed the number of
+required bits, which implies that \code{ceil(n / 8)} bytes must be
 allocated.  The bitset structure is define by:
+Bitsets are required in order to turn options on and off.  We require
+the capability to have a bitset of arbitrary length (i.e., not limited
+by the length of a \code{long}, say).  The \code{psBitset} structure
+is defined below.  Note that the entry \code{bits} is an array of type
+\code{char} storing the bits as bits of each byte in the array, with 8
+bits available for each byte in the array.  Also note that the
+constructor is passed the number of required bits, which implies that
+\code{ceil(n/8)} bytes must be allocated.  The bitset structure is
+define by:
 \begin{verbatim}
 typedef struct {
 …
 where \code{n} is the requested number of bits.
 Four basic operations on bitsets are required:
+Several basic operations on bitsets are required:
 \begin{itemize}
 \item Set a bit;
 \item Check if a bit is set; and
 \item \code{OR}, \code{AND} and \code{XOR} two bitsets.
+\item \code{NOT} a bitset.
 \end{itemize}
 The corresponding APIs are defined below:
 …
 \subsection{Sorting}
 We require the ability to sort an array of floating-point values.  The
+following function returns the array, sorted from the smallest (i.e.\
+most negative) value in the first element, and the largest (i.e.\ most
+positive) value in the last element.  The input array, \code{in}, may
+be sorted in-place if it is also specified as the \code{out}
 array. This function is specified for input types \code{psU8, psU16,
 psF32, psF64}.  The input and output vectors must have the same type.
+We require the ability to sort a vector.  The following function
+returns the vector, sorted from the smallest (i.e.\ most negative)
+value in the first element, and the largest (i.e.\ most positive)
+value in the last element.  The input vector, \code{in}, may be sorted
+in-place if it is also specified as the \code{out} vector. This
+function is specified for input types \code{psU8, psU16, psF32,
+psF64}.  The input and output vectors must have the same type.
 \begin{verbatim}
 …
 \end{verbatim}
 We also require the ability to sort one array based on another.  For
+We also require the ability to sort one vector based on another.  For
 example, we may want to sort both \code{x} and \code{y} by the value
 in \code{x}.  In order to facilitate this, we will have a sort
 function return an array containing the indices for the unsorted list
 in the order appropriate for the sorted array.  The output vector must
+function return a vector containing the indices for the unsorted list
+in the order appropriate for the sorted vector.  The output vector must
 be of type \code{psU32}.  This function is specified for input types
 \code{psU8, psU16, psF32, psF64}.
 …
 \end{verbatim}
 The sorted arrays may be accessed in the following manner:
 \begin{verbatim}
 indexArray = psSortIndex(NULL, x);
 for (int i = 0; i < indexArray.n; i++) {
     doMyFunc(x[indexArray.arr.arr_U32[i]], y[indexArray[i].arr.arr_U32]);
+The sorted vectors may be accessed in the following manner:
+\begin{verbatim}
+indexVector = psSortIndex(NULL, x);
+for (int i = 0; i < indexVector.n; i++) {
+    doMyFunc(x[indexVector.arr.arr_U32[i]], y[indexVector[i].arr.arr_U32]);
+}
 \end{verbatim}
 …
 \subsection{Statistics Functions}
 \subsubsection{Array Statistics}
 We require a very general statistics function, which, given an array
+\subsubsection{Vector Statistics}
+We require a very general statistics function, which, given a vector
 of floating-point values, will be able to calculate the following
 statistics:
 …
 \item Clipped mean and number of values used to calculate;
 \item Clipped standard deviation; and
 \item Minimum and maximum value in array.
+\item Minimum and maximum value in vector.
 \end{itemize}
 For definitions of each of these, see the accompanying Algorithms
 Definition Document (ADD), but in general, ``sample'' refers to the
 entire array, ``robust'' refers to fitting the distribution in the
 array with a model, and ``clipped'' refers to clipping the
 distribution.  Each of these shall be available from a single
+entire vector, ``robust'' refers to fitting the distribution in the
+vector with a model, and ``clipped'' refers to clipping the
+distribution.  Each of these must be available from a single
 function:
 \begin{verbatim}
 psStats *psArrayStats(const psVector *restrict in,
                       const psVector *restrict mask,
                       unsigned int maskVal,
                       psStats *stats);
+psStats *psVectorStats(const psVector *restrict in,
+                       const psVector *restrict mask,
+                       unsigned int maskVal,
+                       psStats *stats);
 \end{verbatim}
+%
 …
 in \code{mask}, so that the user may explicitly reject specific
 entries) and a \code{psStats} structure, which will be altered and
+returned.  The input vector may be of type \code{psU8}, \code{psU16},
+\code{psF32}, \code{psF64}; the mask must be of type \code{psU8}.
+returned.  The \code{psStats} structure includes several fields which
+are used for both input and output: \code{min} and \code{max} may be
+used to specify a value range for which the statistics are calculated.
+\code{binsize} specifies a choice for the robust statistics histogram
+bin size.  If these are to be used, the user should set the
+corresponding \code{options} bits \code{PS_STAT_USE_RANGE} or
+\code{PS_STAT_USE_BINSIZE}.  \code{clipSigma} specifies the number of
+standard deviations for which data should be clipped.  \code{clipIter}
+specifies the number of iterations which should be used for clipping.
+The defaults for these two numbers is both 3.  Since the sample
+statistics scale like $N\log N$, for large numbers of input data
+points, it is faster to use the robust statistics.  If the number of
+data points is large, \code{psStats} must revert to the robust
+calculation even if the user requested sample statistics.  The values
+should be returned in the \code{sample} fields, but the bit
+\code{PS_STAT_ROBUST_FOR_SAMPLE} in \code{options} must be set in this
+case.  The cutoff for this decision must be made on the basis of the
+value in \code{sampleLimit},which should have a default of \tbd{3e5}.
+Default input field values must be set by the \code{psStats}
+constructor.  The input vector may be of type \code{psU8},
+\code{psU16}, \code{psF32}, \code{psF64}; the mask must be of type
+\code{psU8}.
 The \code{psStats} structure is defined with entries for each of the
 …
     double sampleUQ;                    ///< upper quartile of sample
     double sampleLQ;                    ///< lower quartile of sample
+    double sampleLimit;                 ///< Number of datapoints to
     double robustMean;                  ///< robust mean of data
-    int    robustMeanNvalues;           ///< number of measurements used for robust mean
     double robustMedian;                ///< robust median of data
-    int    robustMedianNvalues;         ///< number of measurements used for robust median
     double robustMode;                  ///< Robust mode of data
-    int    robustModeNvalues;           ///< Number of measurements used for robust mode
     double robustStdev;                 ///< robust standard deviation of data
     double robustUQ;                    ///< robust upper quartile
     double robustLQ;                    ///< robust lower quartile
+    int    robustN50;                   ///< Number of points UQ-LQ
+    int    robustNfit;                  ///< Number of points in Gauss. fit
     double clippedMean;                 ///< Nsigma clipped mean
-    int    clippedMeanNvalues;          ///< number of data points used for clipped mean
     double clippedStdev;                ///< standard deviation after clipping
+    int    clippedNvalues;              ///< number of data points used for clipped mean
     double clipSigma;                   ///< Nsigma used for clipping; user input
     int    clipIter;                    ///< Number of clipping iterations; user input
     double min;                         ///< minimum data value in data
     double max;                         ///< maximum data value in data
     int    nValues;                     ///< number of data values in data
+    double min;                         ///< minimum data value in data; input range
+    double max;                         ///< maximum data value in data; input range
+    double binsize;                     ///< binsize for robust fit (input/output)
     psStatsOptions options;             ///< bitmask of calculated values
 } psStats;
 …
     PS_STAT_SAMPLE_MEDIAN         = 0x000002,
     PS_STAT_SAMPLE_STDEV          = 0x000004,
+    PS_STAT_SAMPLE_UQ             = 0x000008,
+    PS_STAT_SAMPLE_LQ             = 0x000010,
+    PS_STAT_ROBUST_MEAN           = 0x000020,
+    PS_STAT_ROBUST_MEAN_NVALUES   = 0x000040,
+    PS_STAT_ROBUST_MEDIAN         = 0x000080,
+    PS_STAT_ROBUST_MEDIAN_NVALUES = 0x000100,
+    PS_STAT_ROBUST_MODE           = 0x000200,
+    PS_STAT_ROBUST_MODE_NVALUES   = 0x000400,
+    PS_STAT_ROBUST_STDEV          = 0x000800,
+    PS_STAT_ROBUST_UQ             = 0x001000,
+    PS_STAT_ROBUST_LQ             = 0x002000,
+    PS_STAT_CLIPPED_MEAN          = 0x004000,
+    PS_STAT_CLIPPED_MEAN_NVALUES  = 0x008000,
+    PS_STAT_CLIPPED_MEAN_NSIGMA   = 0x010000,
+    PS_STAT_CLIPPED_STDEV         = 0x020000,
+    PS_STAT_MAX                   = 0x040000,
+    PS_STAT_MIN                   = 0x080000,
+    PS_STAT_NVALUES               = 0x100000
+    PS_STAT_SAMPLE_QUARTILE       = 0x000008,
+    PS_STAT_ROBUST_MEAN           = 0x000010,
+    PS_STAT_ROBUST_MEDIAN         = 0x000020,
+    PS_STAT_ROBUST_MODE           = 0x000040,
+    PS_STAT_ROBUST_STDEV          = 0x000080,
+    PS_STAT_ROBUST_QUARTILE       = 0x000100,
+    PS_STAT_CLIPPED_MEAN          = 0x000200,
+    PS_STAT_CLIPPED_STDEV         = 0x000400,
+    PS_STAT_MAX                   = 0x000800,
+    PS_STAT_MIN                   = 0x001000,
+    PS_STAT_USE_RANGE             = 0x002000,
+    PS_STAT_USE_BINSIZE           = 0x004000
+    PS_STAT_ROBUST_FOR_SAMPLE     = 0x008000
 } psStatsOptions;
 \end{verbatim}
 …
     const psVector *restrict bounds;    ///< Bounds for the bins
     psVector *nums;                     ///< Number in each of the bins
     int minNum, maxNum;                 ///< Number below the minimum and above the maximum
+    int minNum, maxNum;                 ///< Number below minimum / above maximum
     int uniform;                        ///< Is it a uniform distribution?
 } psHistogram;
 …
 \end{verbatim}
 The input vector may be of types \code{psU8, psU16, psF32, psF64}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Analytical functions}
+\subsubsection{General Polynomials}
+PSLib provides APIs to represent and interact with polynomials in up
+to four dimensions, with both floating-point and double-precision
+numbers.  In Pan-STARRS processing, the astrometry requirements push
+the need for at least four dimensions ($x$,$y$, color and magnitude)
+and double-precision (for milli-arcsec precision) versions.  We must
+also be able to calculate the errors in the fit coefficients, as well
+as be able to turn on and off each coefficient.  This leads us to
+define the following polynomial types:
+\begin{verbatim}
+/** One-dimensional polynomial */
+typedef struct {
+    int n;                              ///< Number of terms
+    float *restrict coeff;              ///< Coefficients
+    float *restrict coeffErr;           ///< Error in coefficients
+    char *restrict mask;                ///< Coefficient mask
+} psPolynomial1D;
+\end{verbatim}
+\begin{verbatim}
+/** Two-dimensional polynomial */
+typedef struct {
+    int nX, nY;                         ///< Number of terms in x and y
+    float *restrict *restrict coeff;    ///< Coefficients
+    float *restrict *restrict coeffErr; ///< Error in coefficients
+    char *restrict *restrict mask;      ///< Coefficients mask
+} psPolynomial2D;
+\end{verbatim}
+etc., up to four dimensions.  We also define double-precision versions:
+\begin{verbatim}
+/** Double-precision one-dimensional polynomial */
+typedef struct {
+    int n;                              ///< Number of terms
+    double *restrict coeff;             ///< Coefficients
+    double *restrict coeffErr;          ///< Error in coefficients
+    char *restrict mask;                ///< Coefficient mask
+} psDPolynomial1D;
+\end{verbatim}
+\begin{verbatim}
+/** Double-precision two-dimensional polynomial */
+typedef struct {
+    int nX, nY;                         ///< Number of terms in x and y
+    double *restrict *restrict coeff;   ///< Coefficients
+    double *restrict *restrict coeffErr; ///< Error in coefficients
+    char *restrict *restrict mask;      ///< Coefficients mask
+} psDPolynomial2D;
+\end{verbatim}
+etc.  In what follows, we only show the version for double-precision
+two-dimensionals; the others may be inferred following the standard
+naming convention exampled above.
+The constructor and destructor are:
+\begin{verbatim}
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY);
+void psDPolynomial2DFree(psDPolynomial2D *restrict myPoly);
+\end{verbatim}
+where \code{nX} and \code{nY} are the number of terms in x and y
+respectively.  The coefficients and errors are set initially to 0.0.
+To evaluate the polynomials at specific coordinates, we define:
+\begin{verbatim}
+double psDPolynomial2DEval(double x, double y, const psDPolynomial2D *restrict myPoly);
+\end{verbatim}
+\subsubsection{Gaussians}
+Gaussians are used extensively in any data-analysis system, in
+particular to represent a real population distribution.  We require
+a function to evaluate a Gaussian for a given coordinate and one which
+generates a Gaussian deviate; a collection of data points whose
+distribution obeys a specified Gaussian.
+The Gaussian evaluation is provide by:
+\begin{verbatim}
+float psGaussian(float x, float mean, float sigma, int normal);
+\end{verbatim}
+which evaluates a Gaussian with the given \code{mean} and \code{sigma}
+at the given coordinate \code{x}.  If \code{normal} is true, the
+Gaussian is normalized (total integral = 1), otherwise, the Gaussian
+is non-normalized (central peak value = 1).  The evaluated Gaussian
+is:
+\[ \frac{1}{\sqrt{2\pi\sigma^2}} exp^{-\frac{(x-mean)^2}{2\sigma^2}} \]
+In the case of the non-normalized Gaussian, the leading coefficient is
+dropped.
+A vector with a specified Gaussian deviate distribution is provide by:
+\begin{verbatim}
+psVector *psGaussianDev(float mean, float sigma, int Npts);
+\end{verbatim}
+which generates a vector (type \code{psF32}) of \code{Npts} elements
+whose distribution has the given \code{mean} and \code{sigma}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Minimization and fitting routines}
+We require a general minimization routine, a routine that will
+specifically minimize $\chi^2$ given a list of data with associated
+errors, and a function that will analytically determine the best
+polynomial fit by minimizing $\chi^2$.
+Consider a function \code{myFunction} which is a function of a
+collection of parameters \code{params} and coordiate vector
+\code{coord}, and returns a single floating point value.  We define
+\code{psMinimize}, which determines the parameters of
+\code{myFunction} which minimize the function for the coordinate
+\code{coord}.  The returned vector must be the same length as the
+vector argument to the input function.
+\begin{verbatim}
+float myFunction (psVector *params, psVector *coord);
+float myFuncDeriv (psVector *params, psVector *coord);
+psVector *psMinimize(psVector *restrict initialGuess,
+                     float (*myFunction)(psVector *, psVector *),
+                     float (*myFuncDeriv)(psVector *, psVector *),
+                     const psVector *restrict coord,
+                     const psVector *restrict paramMask);
+\end{verbatim}
+\code{psMinimize} determines and returns the vector that minimizes the
+specified function.  It takes as input a function, \code{myFunction},
+an initial guess for the vector that minimizes the function,
+\code{initialGuess}, the current coordiate \code{coord}, and an
+optional mask for the vector elements (function parameters) to
+minimize, \code{paramMask} (all parameters are fit if \code{NULL}).
+Note that \code{paramMask} must be of type \code{psU8}, while
+\code{params} must be of type \code{psF32}.  The optional function,
+\code{myFuncDeriv} returns the derivative of the function.
+\begin{verbatim}
+psVector *psMinimizeChi2(psVector *restrict initialGuess,
+                         float (*evalModel)(psVector *, psVector *),
+                         const psVector *restrict domain,
+                         const psVector *restrict data,
+                         const psVector *restrict errors,
+                         const psVector *restrict paramMask,
+                         float *ChiSq);
+\end{verbatim}
+\code{psMinimizeChi2} fits a model to observations by minimizing
+$\chi^2$, returing the best-fit parameters.  The input parameters are
+a function that evaluates the model for a specified domain, given the
+parameters, \code{evalModel}; a list of observations, (\code{domain},
+\code{data}, and \code{errors}); an initial guess at the best-fit
+parameters, \code{initialGuess} which is returned with the best-fit
+parameters, and an optional mask specifying which parameters are to be
+fit, \code{paramMask}, which must be of type \code{psU8}.  All
+parameters are fit if this vector is \code{NULL}.
+\begin{verbatim}
+psPolynomial1D *psVectorFitPolynomial1D(psPolynomial1D myPoly,
+                                     const psVector *restrict x,
+                                     const psVector *restrict y,
+                                     const psVector *restrict yErr);
+\end{verbatim}
+\code{psVectorFitPolynomial} returns the polynomial that best fits the
+observations.  The input parameters are a polynomial that specifies
+the fit order, \code{myPoly}, which will be altered and returned with
+the best-fit coefficients; and the observations, \code{x}, \code{y}
+and \code{yErr}.  The independent variable list, \code{x} may be
+\code{NULL}, in which case the vector index is used.  The dependent
+variable error, \code{yErr} may be null, in which case the solution is
+determined in the assumption that all data errors are equal.  This
+function must be valid for types \code{psU8}, \code{psU16},
+\code{psF32}, \code{psF64}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Image Operations}
+We require a variety of functions to manipulate these image
+structures, including creation, destruction, input, output, and
+various manipulations of the pixels.  The required functions are
+listed below, and fall into several categories.
+\subsubsection{Image Structure Manipulation}
+\begin{verbatim}
+psImage *psImageSubset(const psImage *image, int nx, int ny, int x0, int y0);
+\end{verbatim}
+Define a subimage of the specified area of the given image.  This
+function must return an error if the requested subset area lies
+outside of the parent image.  The argument \code{image} is the parent
+image, \code{nx,ny} specify the dimensions of the desired subraster,
+and \code{x0, y0} specify the starting pixel of the subraster.  The
+entire subraster must be contained within the raster of the parent
+image.  Note that the \code{refCounter} for the parent should be
+incremented.  This function must be defined for the following types:
+\code{psU8}, \code{psU16}, \code{psS8}, \code{psS16}, \code{psF32},
+\code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageCopy(psImage *output, const psImage *input, psElemType type);
+\end{verbatim}
+Create a copy of the specified image, converting the type in the
+process.  If the output target pointer is not NULL, place the result
+in the specified structure.  The output image data must be allocated
+as a single, contiguous block of memory.  The output image may not be
+the input image.  This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\subsubsection{Image Pixel Extractions}
+\begin{verbatim}
+psVector *psImageSlice(psVector *out, const psImage *input,
+                       int x, int y, int nx, int ny,
+                       int direction, const psStats *stats);
+\end{verbatim}
+Extract pixels from rectlinear region to a vector (array of floats).
+The output vector contains either \code{nx} or \code{ny} elements,
+based on the value of the direction: e.g., if \code{direction} is
+\tbd{+x}, there are \code{nx} elements.  The input region is collapsed
+in the perpendicular direction, and each element of the output vectors
+is derived from the statistics of the pixels at that direction
+coordinate.  The statistic used to derive the output vector value is
+specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psVector *psImageCut(psVector *out, const psImage *input,
+                     float xs, float ys, float xe, float ye,
+                     float dw, const psStats *stats);
+\end{verbatim}
+Extract pixels from an image along a line to a vector (array of
+floats).  The vector \code{(xs,ys)} - \code{(xe,ye)} forms the basis of
+the output vector.  Pixels are considered in a rectangular region of
+width \code{dw} about this vector.  The input region is collapsed in
+the perpendicular direction, and each element of the output vector
+represents pixel-sized boxes, where the value is derived from the
+statistics of the pixels interpolated along the perpendicular
+direction.  The statistic used to derive the output vector value is
+specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}
+\begin{verbatim}
+psVector *psImageRadialCut(psVector *out, const psImage *input, float x, float y,
+                           const psVector *radii, const psStats *stats);
+\end{verbatim}
+Extract radial annuli data to a vector.  A vector is constructed where
+each vector elements is derived from the statistics of the pixels
+which land in one of a sequence of annuli.  The annuli are centered on
+the image pixel coordinate \code{x,y}, and have width \code{dr}.  The
+number of annuli is $radius / dr$.  The statistic used to derive the
+output vector value is specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\subsubsection{Image Geometry Manipulation}
+\begin{verbatim}
+psImage *psImageRebin(psImage *out, const psImage *in,
+                      float scale, const psStats *stats);
+\end{verbatim}
+Rebin image to new scale.  A new image is constructed in which the
+dimensions are reduced by a factor of \code{scale} $\le 1$ (it is an
+error for \code{scale} $> 1$).  The \code{scale} is equal in each
+dimension.  The output image is generated from all input image pixels.
+Each pixel in the output image is derived from the statistics of the
+corresponding set of input image pixels based on the statistics
+specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageRotate(psImage *out, const psImage *input, float angle);
+\end{verbatim}
+Rotate the input image by given angle, specified in degrees.  The
+output image must contain all of the pixels from the input image in
+their new frame.  Pixels in the output image which do not map to input
+pixels should be set to \tbd{value}.  The center of rotation is always
+the center pixel of the image.  The rotation is specified in the sense
+that a positive angle is an anti-clockwise rotation.  This function
+must be defined for the following types: \code{psU8}, \code{psU16},
+\code{psS8}, \code{psS16}, \code{psF32}, \code{psF64}, \code{psC32},
+\code{psC64}.
+\begin{verbatim}
+psImage *psImageShift(psImage *out, const psImage *input,
+                      float dx, float dy, float exposed);
+\end{verbatim}
+Shift image by an arbitrary number of pixels (\code{dx,dy}) in either
+direction.  If the shift values are fractional, the output pixel
+values should interpolate between the input pixel values.  The output
+image has the same dimensions as the input image.  Pixels which fall
+off the edge of the output image are lost.  Newly exposed pixels are
+set to the value given by \code{exposed}.  This function must be
+defined for the following types: \code{psU8}, \code{psU16},
+\code{psS8}, \code{psS16}, \code{psF32}, \code{psF64}, \code{psC32},
+\code{psC64}.
+\begin{verbatim}
+psImage *psImageRoll(psImage *out, const psImage *input, int dx, int dy);
+\end{verbatim}
+Roll image by an integer number of pixels (\code{dx,dy}) in either
+direction.  The output image is the same dimensions as the input
+image.  Edge pixels wrap to the other side (no values are lost).  This
+function must be defined for the following types: \code{psU8},
+\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64},
+\code{psC32}, \code{psC64}.
+\subsubsection{Image Statistical Functions}
+\begin{verbatim}
+psStats *psImageGetStats(psStats *stats, const psImage *input);
+\end{verbatim}
+Determine statistics for image (or subimage).  The statistics to be
+determined are specified by \code{stats}.  This function must be
+defined for the following types: \code{psU8}, \code{psU16},
+\code{psF32}, \code{psF64}.
+\begin{verbatim}
+psHistogram *psImageHistogram(psHistogram *hist, const psImage *input);
+\end{verbatim}
+Construct a histogram from an image (or subimage).  The histogram to
+generate is specified by \code{psHistogram hist} (see
+section~\ref{sec:histograms}).  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psPolynomial2D *psImageFitPolynomial(psPolynomial2D *coeffs, const psImage *input);
+\end{verbatim}
+Fit a 2-D Chebychev polynomial surface to an image.  The input
+structure \code{coeffs} contains the desired order and terms of
+interest.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psImage *psImageEvalPolynomial(psImage *input, const psPolynomial2D *coeffs);
+\end{verbatim}
+Evaluate a 2-D polynomial surface for the image pixels.  Given the
+input polynomial coefficients, set the image pixel values on the basis
+of the polynomial function.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\subsubsection{Image I/O Functions}
+\begin{verbatim}
+psImage *psImageReadSection(psImage *output, int x0, int y0, int nx, int ny, int z,
+                            const char *extname, int extnum, const char *filename);
+\end{verbatim}
+Read an image or subimage from a named file.  This function is a
+wrapper to the FITS library function.  The input parameters allow a
+full image or a subimage to be read.  The starting pixel of the region
+is specified by \code{x,y}, while the dimensions of the requested
+region are specified by \code{nx,ny}.  A negative value of -1 for
+these two parameters specifies the full array of the requested image,
+less absolute value of the variable.  If the native image is a cube,
+the value of z specifies the requested slice of the image.  The data
+is read from the extension specified by extname (matching the EXTNAME
+keyword) or by the extnum value (with 0 representing the PHU, 1 the
+first extension, etc).  This function must call \code{psError} and
+return \code{NULL} if any of the specified parameters are out of range
+for the data in the image file, if the specified image file does not
+exist, or the image on disk is zero- or one-dimensional.
+\begin{verbatim}
+psImage *psImageFReadSection(psImage *output, int x, int y, int nx, int ny, int z,
+                             const char *extname, int extnum, FILE *f);
+\end{verbatim}
+Read an image or subimage from file descriptor.  The input parameters
+and their behavior for this function are identical with those in
+\code{psImageReadSection}.
+\begin{verbatim}
+psImage *psImageWriteSection(const psImage *input, int x, int y, int z,
+                             const char *extname, int extnum, const char *filename);
+\end{verbatim}
+Write an image section to the named file, which may exist.  This
+operation may write a portion of an image over the existing bytes of
+an existing image.  If the file does not exist, it should be created.
+If the specified extension does not exist, it should be created.  If
+an extension is specified and no PHU exists, a basic PHU should be
+created.
+\begin{verbatim}
+psImage *psImageFWriteSection(const psImage *input, int x, int y, int z,
+                              const char *extname, int extnum, FILE *f);
+\end{verbatim}
+Write an image section to file descriptor as above:
+\subsubsection{Image Pixel Manipulations}
+\begin{verbatim}
+int psImageClip(psImage *input, float min, float vmin, float max, float vmax);
+\end{verbatim}
+Clip image values outside of range to given values.  All pixels with
+values \code{< min} are set to the value \code{vmin}. All pixels with
+values \code{> max} are set to the value \code{vmax}. Returns the
+number of clipped pixels.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+int psImageClipNaN(psImage *input, float value);
+\end{verbatim}
+Clip \code{NaN} image pixels to given value.  Pixels with \code{NaN},
+\code{+Inf} or \code{-Inf} values are set to the specified value.
+Returns the number of clipped pixels.  This function must be defined
+for the following types: \code{psU8}, \code{psU16}, \code{psS8},
+\code{psS16}, \code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+int psImageOverlaySection(psImage *image, const psImage *overlay,
+                          int x0, int y0, const char *op);
+\end{verbatim}
+Overlay subregion of image with another image.  Replace the pixels in
+the \code{image} which correspond to the pixels in \code{overlay} with
+values derived from the values in \code{image} and \code{overlay}
+based on the given operator \code{op}.  Valid operators are \code{=}
+(set image value to overlay value), \code{+} (add overlay value to
+image value), \code{-} (subtract overlay from image), \code{*}
+(multiply overlay times image), \code{/} (divide image by overlay).
+This function must be defined for the following types: \code{psU8},
+\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64},
+\code{psC32}, \code{psC64}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Vector and Image Arithmetic}
+\label{sec:arithmetic}
+We will need to be able to perform various operations on vectors and
+images, e.g.\ dividing one image by another, subtracting a vector
+from an image, etc.  Both binary operations and unary operations are
+required.  To avoid the burden of memorizing a ton of APIs, we specify
+two generic APIs for the binary and unary operations.
+\begin{verbatim}
+psType *psBinaryOp (void *out, void *in1, char *op, void *in2);
+psType *psUnaryOp (void *out, void *in, char *op);
+\end{verbatim}
+These functions determine the type of the operands on the basis of
+their \code{psType} elements, which always are the first elements.
+Note that these functions return a pointer to the appropriate type for
+the operation.  Since the result may is cast to \code{psType}, the
+resulting type may be determined by examining the return value.  It is
+expected that the implementation of these functions will employ
+pre-processor macros to perform the onerous task of creating the
+loops.  Also note that \code{psVector} is equivalent to
+\code{psVector}.  An attempt to perform an arithmetic operation on
+an object of dimension \code{PS_DIMEN_OTHER} should produce an error.
+Binary operations between an image and a vector have a potential
+ambiguity --- do the vector elements correspond to the rows or the
+columns?  For this reason, we define two vector types: a ``vector''
+(\code{PS_DIMEN_VECTOR}), and a ``transposed vector''
+(\code{PS_DIMEN_TRANSV}).  We specify that a ``vector'', when involved
+in binary operations on an image, acts on the rows, while a
+``transposed vector'' in the same context acts on the columns.
+Vectors, when created, will be created as ``vectors'', but may be
+converted to ``transposed vectors'' using the following function:
+\begin{verbatim}
+psVector *psVectorTranspose(psVector *out, psVector *myVector);
+\end{verbatim}
+It is further desirable to allow scalar values to be used within these
+functions, which requires the following additions:
+\begin{verbatim}
+p_ps_Scalar *psScalar (double value);
+p_ps_Scalar *psScalarType (char *mode, ...);
+\end{verbatim}
+The first creates a psType-ed structure from a constant value, while
+the second creates a psType-ed structure for a specified type.  The
+structure which carries a scalar value is specified as the following
+private type, and is analogous to the \code{psVector} and
+\code{psImage} structures:
+\begin{verbatim}
+typedef struct {
+    psType type;                        ///< data type information
+    union {
+        psS32 S32;                      ///< integer value entry
+        psF32 F32;                      ///< float value entry
+        psF64 F64;                      ///< double value entry
+        psC32 C32;                      ///< complex value entry
+    } data;
+} p_psScalar;
+\end{verbatim}
+This allows one to write the following to take the sine of the square
+of all pixels in an image:
+\begin{verbatim}
+psImage A,B;
+B = psBinaryOp (NULL, A, "^", psScalar(2));
+(void) psUnaryOp(B, B, "sin");
+\end{verbatim}
+Note that the \code{psUnaryOp} is performed on \code{B} in-place.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 …
 \end{verbatim}
+Matrix arithmetic is supported through \code{psMatrixOp}.  Note that
+\code{psMatrixOp} differs from \code{psBinaryOp} in that
+multiplication with \code{psBinaryOp} acts on corresponding elements,
+while \code{psMatrixOp} performs classical matrix multiplication
+involving row and column operations.  Addition and subtraction by
+\code{psMatrixOp} are identical to that for \code{psBinaryOp} (since
+that is how matrix addition is defined).  Matrix division is not
+defined.  \tbd{why specify addition as an operator? why not only
+define operators which are different for matrices?}.  This function is
+Matrix multiplication is supported through \code{psMatrixMultiply}.  This function is
 specified for data types \code{psF32, psF64}.
+\begin{verbatim}
+psImage *psMatrixOp(psImage *out, const psImage *in1, const char *op, const psImage *in2);
+\begin{verbatim}
+psImage *psMatrixMultiply(psImage *out, const psImage *in1, const psImage *in2);
 \end{verbatim}
 …
 always \code{psC32}.  If the input vector is \code{psF32}, the
 direction must be forward.  Neither the forward or inverse transforms
 shall multiply by $1/N$ (or $1/N^{1/2}$), and so it falls to the
+must multiply by $1/N$ (or $1/N^{1/2}$), and so it falls to the
 responsibility of the user to multiply a vector that has been forward-
 and reverse-transformed by $1/N$.
 …
 \begin{verbatim}
 psImage *psImageFFT(const psImage *image, int direction);
-psImage *psImagePowerSpectrum(const psImage *image);
 psImage *psImageReal(psImage *out, const psImage *in);
 psImage *psImageReal(psImage *out, const psImage *in);
 psImage *psImageReal(psImage *out, const psImage *in);
+psImage *psImageImaginary(psImage *out, const psImage *in);
+psImage *psImageComplex(psImage *real, const psImage *imag);
 psImage *psImageConjugate(psImage *out, const psImage *in);
+\end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Analytical functions}
+We require two types of general functions which will be used in fitting:
+Gaussians and Polynomials.  The Gaussian is defined as:
+\begin{verbatim}
+float psGaussian(float x, float mean, float stdev);
+\end{verbatim}
+which evaluates a non-normalized Gaussian with the given mean and
+sigma at the given coordianate.  Note that this is not a Gaussian
+deviate.  The evaluated Gaussian is:
+\[ 1/(\sqrt(2\pi)\sigma) exp(-\frac{(x-mean)^2}{2\sigma^2}) \]
+For the polynomial, \PS{} astrometry requirements lead us to specify
+that we must be able to support at least four dimensions, in both
+floating-point (for speed) and double-precision (for milli-arcsec
+precision) versions.  This comes from the need to fit over $x$,$y$,
+color and magnitude simultaneously.  We must also be able to
+calculate the errors in the fit coefficients, as well as be able to
+turn on and off each coefficient.  This leads us to define polynomial
+types:
+\begin{verbatim}
+/** One-dimensional polynomial */
+typedef struct {
+    int n;                              ///< Number of terms
+    float *restrict coeff;              ///< Coefficients
+    float *restrict coeffErr;           ///< Error in coefficients
+    char *restrict mask;                ///< Coefficient mask
+} psPolynomial1D;
+\end{verbatim}
+\begin{verbatim}
+/** Two-dimensional polynomial */
+typedef struct {
+    int nX, nY;                         ///< Number of terms in x and y
+    float *restrict *restrict coeff;    ///< Coefficients
+    float *restrict *restrict coeffErr; ///< Error in coefficients
+    char *restrict *restrict mask;      ///< Coefficients mask
+} psPolynomial2D;
+\end{verbatim}
+etc., up to four dimensions.  We also define double-precision versions:
+\begin{verbatim}
+/** Double-precision one-dimensional polynomial */
+typedef struct {
+    int n;                              ///< Number of terms
+    double *restrict coeff;             ///< Coefficients
+    double *restrict coeffErr;          ///< Error in coefficients
+    char *restrict mask;                ///< Coefficient mask
+} psDPolynomial1D;
+\end{verbatim}
+\begin{verbatim}
+/** Double-precision two-dimensional polynomial */
+typedef struct {
+    int nX, nY;                         ///< Number of terms in x and y
+    double *restrict *restrict coeff;   ///< Coefficients
+    double *restrict *restrict coeffErr; ///< Error in coefficients
+    char *restrict *restrict mask;      ///< Coefficients mask
+} psDPolynomial2D;
+\end{verbatim}
+etc.  In what follows, we only show the version for double-precision
+two-dimensionals; the others may be inferred following the standard
+naming convention exampled above.
+The constructor and destructor are:
+\begin{verbatim}
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY);
+void psDPolynomial2DFree(psDPolynomial2D *restrict myPoly);
+\end{verbatim}
+where \code{nX} and \code{nY} are the number of terms in x and y
+respectively.
+To evaluate the polynomials at specific coordinates, we define:
+\begin{verbatim}
+double psEvalDPolynomial2D(double x, double y, const psDPolynomial2D *restrict myPoly);
+\end{verbatim}
+\tbd{Generate a vector of random Gaussian deviates}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Minimization and fitting routines}
+\tbd{I am concerned the minimization functions defined here are
+  missing details - gene}
+We require a general minimization routine, a routine that will
+specifically minimize $\chi^2$ given a list of data with associated
+errors, and a function that will analytically determine the best
+polynomial fit by minimizing $\chi^2$.
+First, we define \code{psMinimize}, which returns the vector which
+minimizes the specified function, which is itself a function of a
+single vector.  The returned vector must be the same length as the
+vector argument to the input function.
+\begin{verbatim}
+psVector *psMinimize(psVector *restrict initialGuess,
+                     float (*myFunction)(const psVector *restrict),
+                     float (*myFuncDeriv)(const psVector *restrict),
+                     const psVector *restrict paramMask);
+\end{verbatim}
+\code{psMinimize} determines and returns the vector that minimizes the
+specified function.  It takes as input a function, \code{myFunction},
+an initial guess for the vector that minimizes the function,
+\code{initialGuess}, and an optional mask for the vector elements
+(function parameters) to minimize, \code{paramMask} (all parameters
+are fit if \code{NULL}).  Note that \code{paramMask} must be of type
+\code{psU8}, while \tbd{what are valid types for the function
+parameter lists?}.  The optional function, \code{myFuncDeriv} returns
+the derivative of the function.  \tbd{is this sufficient? clear?}
+\begin{verbatim}
+psVector *psMinimizeChi2(psVector *restrict initialGuess,
+                         float (*evalModel)(const psVector *restrict, const psVector *restrict),
+                         const psVector *restrict domain,
+                         const psVector *restrict data,
+                         const psVector *restrict errors,
+                         const psVector *restrict paramMask);
+\end{verbatim}
+\code{psMinimizeChi2} fits a model to observations by minimizing
+$\chi^2$, returing the best-fit parameters.  The input parameters are
+a function that evaluates the model for a specified domain, given the
+parameters, \code{evalModel}; a list of observations, (\code{domain},
+\code{data}, and \code{errors}); an initial guess at the best-fit
+parameters, \code{initialGuess} which is returned with the best-fit
+parameters \tbd{how? unclear!}, and an optional mask specifying which
+parameters are to be fit, \code{paramMask}, which must be of type
+\code{psU8}.  All parameters are fit if this vector is \code{NULL}.
+\begin{verbatim}
+psPolynomial1D *psGetArrayPolynomial(psPolynomial1D myPoly,
+                                     const psVector *restrict x,
+                                     const psVector *restrict y,
+                                     const psVector *restrict yErr);
+\end{verbatim}
+\code{psGetArrayPolynomial} returns the polynomial that best fits the
+observations.  The input parameters are a polynomial that specifies
+the fit order, \code{myPoly}, which will be altered and returned with
+the best-fit coefficients; and the observations, \code{x}, \code{y}
+and \code{yErr}.  The independent variable list, \code{x} may be
+\code{NULL}, in which case the vector index is used.  The dependent
+variable error, \code{yErr} may be null, in which case the solution is
+determined in the assumption that all data errors are equal.
+\tbd{valid types?}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Image Operations}
+\subsubsection{Image Structure Manipulation}
+We require a variety of functions to manipulate these image
+structures, including creation, destruction, input, output, and
+various manipulations of the pixels.  The required functions are
+listed below, and fall into several categories.
+Define a subimage of the specified area of the given image.  This
+function must return an error if the requested subset area lies
+outside of the parent image.
+\begin{verbatim}
+psImage *psImageSubset(const psImage *image, int nx, int ny, int x0, int y0);
+\end{verbatim}
+where \code{image} is the parent image, \code{nx,ny} specify the
+dimensions of the desired subraster, and \code{x0, y0} specify the
+starting pixel of the subraster.  The entire subraster must be
+contained within the raster of the parent image.  Note that the
+\code{refCounter} for the parent should be incremented.  This function
+must be defined for the following types: \code{psU8}, \code{psU16},
+\code{psS8}, \code{psS16}, \code{psF32}, \code{psF64}, \code{psC32},
+\code{psC64}.
+Create a copy of the specified image, converting the type in the
+process.  If the output target pointer is not NULL, place the result
+in the specified structure.  The output image data must be allocated
+as a single, contiguous block of memory.  The output image may not be
+the input image.  This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageCopy(psImage *output, const psImage *input, psElemType type);
+\end{verbatim}
+\subsubsection{Image Pixel Extractions}
+Extract pixels from rectlinear region to a vector (array of floats).
+The output vector contains either \code{nx} or \code{ny} elements,
+based on the value of the direction: e.g., if \code{direction} is
+\tbd{+x}, there are \code{nx} elements.  The input region is collapsed
+in the perpendicular direction, and each element of the output vectors
+is derived from the statistics of the pixels at that direction
+coordinate.  The statistic used to derive the output vector value is
+specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psVector *psImageSlice(psVector *out, const psImage *input, int x, int y, int nx, int ny,
+                           int direction, const psStats *stats);
+\end{verbatim}
+Extract pixels from an image along a line to a vector (array of
+floats).  The vector \code{(xs,ys)} - \code{(xe,ye)} forms the basis of
+the output vector.  Pixels are considered in a rectangular region of
+width \code{dw} about this vector.  The input region is collapsed in
+the perpendicular direction, and each element of the output vector
+represents pixel-sized boxes, where the value is derived from the
+statistics of the pixels interpolated along the perpendicular
+direction.  The statistic used to derive the output vector value is
+specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}
+\begin{verbatim}
+psVector *psImageCut(psVector *out, const psImage *input, float xs, float ys, float xe, float ye,
+                         float dw, const psStats *stats);
+\end{verbatim}
+Extract radial annuli data to a vector.  A vector is constructed where
+each vector elements is derived from the statistics of the pixels
+which land in one of a sequence of annuli.  The annuli are centered on
+the image pixel coordinate \code{x,y}, and have width \code{dr}.  The
+number of annuli is $radius / dr$.  The statistic used to derive the
+output vector value is specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psVector *psImageRadialCut(psVector *out, const psImage *input, float x, float y,
+                               const psVector *radii, const psStats *stats);
+\end{verbatim}
+\subsubsection{Image Geometry Manipulation}
+Rebin image to new scale.  A new image is constructed in which the
+dimensions are reduced by a factor of \code{scale} $\le 1$ (it is an
+error for \code{scale} $> 1$).  The \code{scale} is equal in each
+dimension.  The output image is generated from all input image pixels.
+Each pixel in the output image is derived from the statistics of the
+corresponding set of input image pixels based on the statistics
+specified by \code{stats}.  This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageRebin(psImage *out, const psImage *input, float scale, const psStats *stats);
+\end{verbatim}
+Rotate the input image by given angle, specified in degrees.  The
+output image must contain all of the pixels from the input image in
+their new frame.  Pixels in the output image which do not map to input
+pixels should be set to \tbd{value}.  The center of rotation is always
+the center pixel of the image.  The rotation is specified in the sense
+that a positive angle is an anti-clockwise rotation.    This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageRotate(psImage *out, const psImage *input, float angle);
+\end{verbatim}
+Shift image by an arbitrary number of pixels (\code{dx,dy}) in either
+direction.  If the shift values are fractional, the output pixel
+values should interpolate between the input pixel values.  The output
+image has the same dimensions as the input image.  Pixels which fall
+off the edge of the output image are lost.  Newly exposed pixels are
+set to the value given by \code{exposed}.    This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageShift(psImage *out, const psImage *input, float dx, float dy, float exposed);
+\end{verbatim}
+Roll image by an integer number of pixels (\code{dx,dy}) in either
+direction.  The output image is the same dimensions as the input
+image.  Edge pixels wrap to the other side (no values are lost).  This function must be defined for the following
+types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+psImage *psImageRoll(psImage *out, const psImage *input, int dx, int dy);
+\end{verbatim}
+\subsubsection{Image Statistical Functions}
+Determine statistics for image (or subimage).  The statistics to be
+determined are specified by \code{stats}.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psStats *psImageGetStats(psStats *stats, const psImage *input);
+\end{verbatim}
+Construct a histogram from an image (or subimage).  The histogram to
+generate is specified by \code{psHistogram hist} (see
+section~\ref{sec:histograms}).  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psHistogram *psImageHistogram(psHistogram *hist, const psImage *input);
+\end{verbatim}
+Fit a 2-D Chebychev polynomial surface to an image.  The input
+structure \code{coeffs} contains the desired order and terms of
+interest.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psPolynomial2D *psImageFitPolynomial(psPolynomial2D *coeffs, const psImage *input);
+\end{verbatim}
+Evaluate a 2-D polynomial surface for the image pixels.  Given the
+input polynomial coefficients, set the image pixel values on the basis
+of the polynomial function.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psF32}, \code{psF64}.
+\begin{verbatim}
+psImage *psImageEvalPolynomial(psImage *input, const psPolynomial2D *coeffs);
+\end{verbatim}
+\subsubsection{Image I/O Functions}
+Read an image or subimage from a named file.  This function is a
+wrapper to the FITS library function.  The input parameters allow a
+full image or a subimage to be read.  The starting pixel of the region
+is specified by \code{x,y}, while the dimensions of the requested
+region are specified by \code{nx,ny}.  A negative value of -1 for
+these two parameters specifies the full array of the requested image,
+less absolute value of the variable.  If the native image is a cube,
+the value of z specifies the requested slice of the image.  The data
+is read from the extension specified by extname (matching the EXTNAME
+keyword) or by the extnum value (with 0 representing the PHU, 1 the
+first extension, etc).  This function must call \code{psError} and
+return \code{NULL} if any of the specified parameters are out of range
+for the data in the image file, if the specified image file does not
+exist, or the image on disk is zero- or one-dimensional.
+\begin{verbatim}
+psImage *psImageReadSection(psImage *output, int x0, int y0, int nx, int ny, int z,
+                            const char *extname, int extnum, const char *filename);
+\end{verbatim}
+Read an image or subimage from file descriptor.  The input parameters
+and their behavior for this function are identical with those in
+\code{psImageReadSection}.
+\begin{verbatim}
+psImage *psImageFReadSection(psImage *output, int x, int y, int nx, int ny, int z,
+                             const char *extname, int extnum, FILE *f);
+\end{verbatim}
+\tbd{The use of \code{FILE*} to carry around the file descriptor is to be reviewed.}
+Write an image section to the named file, which may exist.  This
+operation may write a portion of an image over the existing bytes of
+an existing image.  If the file does not exist, it should be created.
+If the specified extension does not exist, it should be created.  If
+an extension is specified and no PHU exists, a basic PHU should be
+created.
+\begin{verbatim}
+psImage *psImageWriteSection(const psImage *input, int x, int y, int z,
+                             const char *extname, int extnum, const char *filename);
+\end{verbatim}
+Write an image section to file descriptor as above:
+\begin{verbatim}
+psImage *psImageFWriteSection(const psImage *input, int x, int y, int z,
+                              const char *extname, int extnum, FILE *f);
+\end{verbatim}
+\tbd{The use of \code{FILE*} to carry around the file descriptor is to be reviewed.}
+Read header data from a FITS image file into a \code{psMetaData}
+structure (see section~\ref{sec:metadata}.  The \code{extname} and
+\code{extnum} parameters specify the extension of interest as above.
+If the named extension does not exist, the function should return an
+error.
+\begin{verbatim}
+psMetadata *psImageReadHeader(psMetadata *output, const char *extname, int extnum, const char *filename);
+\end{verbatim}
+Read header data from a FITS image file descriptor into a
+\code{psMetaData} structure.
+\begin{verbatim}
+psMetadata *psImageFReadHeader(psMetadata *output, const char *extname, int extnum, FILE *f);
+\end{verbatim}
+\subsubsection{Image Pixel Manipulations}
+Clip image values outside of range to given values.  All pixels with
+values \code{< min} are set to the value \code{vmin}. All pixels with
+values \code{> max} are set to the value \code{vmax}. Returns the
+number of clipped pixels.  This function must be defined for the
+following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
+\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+int psImageClip(psImage *input, float min, float vmin, float max, float vmax);
+\end{verbatim}
+Clip \code{NaN} image pixels to given value.  Pixels with \code{NaN},
+\code{+Inf} or \code{-Inf} values are set to the specified value.
+Returns the number of clipped pixels.  This function must be defined
+for the following types: \code{psU8}, \code{psU16}, \code{psS8},
+\code{psS16}, \code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.
+\begin{verbatim}
+int psImageClipNaN(psImage *input, float value);
+\end{verbatim}
+Overlay subregion of image with another image.  Replace the pixels in
+the \code{image} which correspond to the pixels in \code{overlay} with
+values derived from the values in \code{image} and \code{overlay}
+based on the given operator \code{op}.  Valid operators are \code{=}
+(set image value to overlay value), \code{+} (add overlay value to
+image value), \code{-} (subtract overlay from image), \code{*}
+(multiply overlay times image), \code{/} (divide image by overlay).
+This function must be defined for the following types: \code{psU8},
+\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64},
+\code{psC32}, \code{psC64}.
+\begin{verbatim}
+int psImageOverlaySection(psImage *image, const psImage *overlay, int x0, int y0, const char *op);
+\end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Vector and Image Arithmetic}
+\label{sec:arithmetic}
+We will need to be able to perform various operations on vectors and
+images, e.g.\ dividing one image by another, subtracting a vector
+from an image, etc.  Both binary operations and unary operations are
+required.  To avoid the burden of memorizing a ton of APIs, we specify
+two generic APIs for the binary and unary operations.
+\begin{verbatim}
+psType *psBinaryOp (void *out, void *in1, char *op, void *in2);
+psType *psUnaryOp (void *out, void *in, char *op);
+\end{verbatim}
+These functions determine the type of the operands on the basis of
+their \code{psType} elements, which always are the first elements.
+Note that these functions return a pointer to the appropriate type for
+the operation.  Since the result may is cast to \code{psType}, the
+resulting type may be determined by examining the return value.  It is
+expected that the implementation of these functions will employ
+pre-processor macros to perform the onerous task of creating the
+loops.  Also note that \code{psVector} is equivalent to
+\code{psVector}.  An attempt to perform an arithmetic operation on
+an object of dimension \code{PS_DIMEN_OTHER} should produce an error.
+Binary operations between an image and a vector have a potential
+ambiguity --- do the vector elements correspond to the rows or the
+columns?  For this reason, we define two vector types: a ``vector''
+(\code{PS_DIMEN_VECTOR}), and a ``transposed vector''
+(\code{PS_DIMEN_TRANSV}).  We specify that a ``vector'', when involved
+in binary operations on an image, acts on the rows, while a
+``transposed vector'' in the same context acts on the columns.
+Vectors, when created, will be created as ``vectors'', but may be
+converted to ``transposed vectors'' using the following function:
+\begin{verbatim}
+/** Transpose a vector.  Changes the type to a PS_DIMEN_TRANSV */
+psVector *psVectorTranspose(psVector *out, //!< Output vector, or NULL
+                            psVector *myVector //!< Vector to be transposed
+    );
+\end{verbatim}
+It is further desirable to allow scalar values to be used within these
+functions, which requires the following additions:
+\begin{verbatim}
+/** create a psType-ed structure from a constant value. */
+p_ps_Scalar *
+psScalar (double value);
+\end{verbatim}
+\begin{verbatim}
+/** create a psType-ed structure from a specified type  */
+p_ps_Scalar *
+psScalarType (char *mode,               ///< type description
+              ...                       ///< value (or values) of specified types
+);
+\end{verbatim}
+\begin{verbatim}
+/** private structure used to pass constant values into the math operators. */
+typedef struct {
+    psType type;                        ///< data type information
+    union {
+        int i;                          ///< integer value entry
+        float f;                        ///< float value entry
+        double d;                       ///< double value entry
+        complex float c;                ///< complex value entry
+    } val;
+} p_psScalar;
+\end{verbatim}
+This allows one to write the following to take the sine of the square
+of all pixels in an image:
+\begin{verbatim}
+psImage A,B;
+B = psBinaryOp (NULL, A, "^", psScalar(2));
+(void) psUnaryOp(B, B, "sin");
+\end{verbatim}
+Note that the \code{psUnaryOp} is performed on \code{B} in-place.
+psImage *psImagePowerSpectrum(const psImage *in);
+\end{verbatim}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 …
 \subsection{Dates and times}
+\textbf{[May be deferred.]}
+We require a collection of functions to manipulate time data.  These
+operations primarily consist of conversions between specific time
+formats.  PSLib wraps the libTAI structure \code{taia} for use as the
+native time format.
+\begin{verbatim}
+typedef struct {
+  struct taia t;
+} psTime;
+\end{verbatim}
+Two utility functions provide the current time in two formats, MJD
+(modified Julian day) and the Sidereal time.
+\begin{verbatim}
+psTime psGetMJD(void);
+psTime psGetSidereal(float mjd, float longitude);
+\end{verbatim}
+A collection of functions convert from the native \code{psTime} format
+to various external formats.  Note that ISO Time is represented by
+YYYY/MM/DD,HH:MM:SS.SSS.
+\begin{verbatim}
+char *psTimeToISOTime (psTime time);
+double psTimeToUTC (psTime time);
+double psTimeToMJD (psTime time);
+double psTimeToJD (psTime time);
+struct timeval *psTimeToTimeval (psTime time);
+struct tm *psTimeToTm (psTime time);
+\end{verbatim}
+A matching set of functions convert from the external formats to the
+native \code{psTime}:
+\begin{verbatim}
+psTime *psISOTimeToTime (char *input);
+psTime *psUTCToTime (double input);
+psTime *psMJDToTime (double input);
+psTime *psJDToTime (double input);
+psTime *psTimevalToTime (struct timeval *input);
+psTime *psTMtoTime (struct tm *input);
+\end{verbatim}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 …
 readable and it may slow down the access.
+We propose an intermediate solution to this problem.  We specify a
 flexible, generic metadata container and access methods.  Data types
 which require access to a general collection of metadata should
+PSLib implements an intermediate solution to this problem.  We specify
+a flexible, generic metadata container and access methods.  Data types
+which require association with a general collection of metadata should
 include an entry of this metadata type.  However, a subset of metadata
 concepts which are basic and frequently required may be placed in the
 coded structure elements.  This approach allows the code to refer to
 the basic metadata concepts as part of the data structure (ie,
+image.nx), but also allows us to provide access to any arbitrary
+metadata which may be generated.  As a practical matter, the choice of
+which entries are only in the metadata and which are part of the
+explicit structure elements is rather subjective.  Any data elements
+which are frequently used should be put in the structure; those which
+are only infrequently needed should be left in the generic metadata.
+\code{image.nx}), but also allows us to provide access to any
+arbitrary metadata which may be generated.  As a practical matter, the
+choice of which entries are only in the metadata and which are part of
+the explicit structure elements is rather subjective.  Any data
+elements which are frequently used should be put in the structure;
+those which are only infrequently needed should be left in the generic
+metadata.
 There are some points of caution which must be noted.  Any
 manipulation of the data should be reflected in the metadata where
 appropriate.  This is always an issue of concern.  For example,
 consider an image of dimensions \code{nX, nY}.  If a function extracts
 a subraster, it must change the values of \code{nX, nY} to match the
+consider an image of dimensions \code{nx, ny}.  If a function extracts
+a subraster, it must change the values of \code{nx, ny} to match the
 new dimensions.  What should it do to the corresponding metadata?
 Clearly, it should change the corresponding value which defines
 …
 copy of the metadata which may be shared by other representations of
 the image?  These must be treated differently because the change would
+invalidate those other references.
+invalidate those other references.  Care must be taken, therefore,
+when writing functions which operate on the data to consider all of
+the relevant metadata entries which must also be updated.
 A related issue is the definition of metadata names.  Entries in a
 …
 \subsubsection{Metadata Representation}
-\tbd{descriptions in this section need to be clarified}
 This section addresses the question of how \PS{} metadata should be
+represented in memory. We do not (yet) address how it should be
+represented on disk.
+We propose that an item of metadata be represented as a C structure with at least the following
+fields:
+\begin{verbatim}
+/** A struct to define a single item of metadata */
+represented in memory, not how it should be represented on disk.
+We define an item of metadata with the following structure:
+\begin{verbatim}
 typedef struct {
     const int id;                       ///< unique ID for this item
 …
     psMetadataFlags flags;              ///< flags associated with this item
     const union {
         double d;                       ///< double value
         float f;                        ///< floating value
         int i;                          ///< integer value
         void *v;                        ///< other type
     } val;                              ///< value of metadata
+        psS32 S32;                      ///< integer value
+        psF32 F32;                      ///< floating value
+        psF64 F64;                      ///< double value
+        void *void;                     ///< other type
+    } data;                             ///< value of metadata
     char *comment;                      ///< optional comment ("", not NULL)
     psDlist *restrict items;            ///< list of psMetadataItems with the same name
 …
 \end{verbatim}
+The \code{id} is a unique identifier for this item of metadata; experience
+shows that such tags are useful.
+The \code{psMetadataType type} entry specifies the type of the data
+being represented; the possibilities are listed in section \ref{metadataTypes}.
+\paragraph{Possible Types of Metadata}
+\label{metadataTypes}
+The possible types of metadata are identified by the enumerated type
+\code{psMetadataType} (see also table \ref{tabMetaDataTypes}); the initial defined values are:
+\begin{verbatim}
+/** Possible types of metadata. */
+The \code{id} is a unique identifier for this item of metadata;
+experience shows that such tags are useful.  The entry \code{name}
+specifies the name of the metadata item.  Metadata naming conventions
+for the Pan-STARRS IPP are specified in the IPP Software Requirements
+Specification (PSDC-430-005).  The value of the metadata is given by
+the union \code{data}, and may be of type \code{psS32}, \code{psF32},
+\code{psF64}, or an arbitrary rich structure pointed at by the
+\code{void} pointer \code{void}.  A character string comment
+associated with this metadata item may be stored in the element
+\code{comment}. The \code{type} entry specifies the type of the data
+being represented, given by the enumerated type \code{psMetadataType}:
+\begin{verbatim}
 typedef enum {                          ///< type of val is:
     PS_META_ITEM_SET = 0,               ///< NULL; metadata is in psMetadataType.items
+    PS_META_FLOAT,                      ///< float (.f)
+    PS_META_INT,                        ///< int (.i)
+    PS_META_STR,                        ///< string (.v)
+    PS_META_IMG,                        ///< image (.v)
+    PS_META_JPEG,                       ///< JPEG (.v)
+    PS_META_PNG,                        ///< PNG (.v)
+    PS_META_ASTROM,                     ///< astrometric coefficients (.v)
+    PS_META_UNKNOWN,                    ///< other (.v)
+    PS_META_S32,                        ///< int (.S32)
+    PS_META_F32,                        ///< float (.F32)
+    PS_META_F64,                        ///< double (.F64)
+    PS_META_STR,                        ///< string (.void)
+    PS_META_IMG,                        ///< image (.void)
+    PS_META_JPEG,                       ///< JPEG (.void)
+    PS_META_PNG,                        ///< PNG (.void)
+    PS_META_ASTROM,                     ///< astrometric coefficients (.void)
+    PS_META_UNKNOWN,                    ///< other (.void)
     PS_META_NTYPE                       ///< Number of types; must be last
 } psMetadataType;
 \end{verbatim}
+\begin{table}
+\begin{tabular}{llll}
+\textbf{Value} & \textbf{Type} & \textbf{member of union} & \textbf{Comments}\\
+\hline
+\code{PS_META_FLOAT}   & float    & f & value, not pointer, is stored \\
+\code{PS_META_INT}     & int      & i & value, not pointer, is stored \\
+\code{PS_META_STR}     & string   & v & value, not pointer to original, is stored \\
+\code{PS_META_IMG}     & psImage  & v & \\
+\code{PS_META_JPEG}    & JPEG     & v & \\
+\code{PS_META_PNG}     & PNG      & v & \\
+\code{PS_META_ASTROM}  & psAstrom & v & \\
+\code{PS_META_UNKNOWN} & other    & v & \\
+\code{PS_META_NTYPE}   & (none)   &   & The number of types defined
+\end{tabular}
+\begin{caption}{Supported Metadata Types}
+\label{tabMetaDataTypes}
+  Possible types of metadata
+\end{caption}
+\end{table}
+\subsubsection{Collections of Metadata}
+\begin{verbatim}
+/** A set of metadata */
+The \code{flags} entry specifies some optional characteristics of the
+metadata type, given by the enumerated type \code{psMetadataFlags}:
+\begin{verbatim}
+typedef enum {
+    PS_META_TYPE_MASK =     0xffff,     ///< the type enum must fit in this mask
+    PS_META_UNIQUE =       0x10000,     ///< the name must be unique (default)
+    PS_META_NON_UNIQUE =   0x20000,     ///< the name may be repeated
+} psMetadataFlags;
+\end{verbatim}
+A collection of metadata is represented by the \code{psMetadata} structure:
+\begin{verbatim}
 typedef struct {
     psDlist *restrict list;             ///< list of psMetadataItem
 …
 } psMetadata;
 \end{verbatim}
 The type \code{psMetadata} is a container class for metadata. Note
 that there are in fact \emph{two} representations of the metadata
 …
 second representation employs a hash table which allows fast look-up
 given a specific metadata keyword.
+Certain metadata names (such as the FITS keywords \code{COMMENT} and
+\code{HISTORY} in a FITS header) may be repeated with different
+values.  The \code{psMetadataAppend} routine is required to check that
+all metadata names are unique unless the type is qualified as
+\code{PS_META_NON_UNIQUE}; in this case a unique integer must be added
+to each name specified.
+\subsubsection{Metadata APIs}
+In this section, we explain the metadata APIs more fully.
+The allocator for \code{psMetadataItem} returns a full
+\code{psMetadataItem} ready for insertion into the \code{psMetadata}.
+The \code{name} entry specifies the name to use for this metadata
+item, and may include \code{sprintf}-stype formating codes.  The
+\code{format} entry, which specifies both the metadata type and the
+optional flags, is constructed by bit-wise or'ing the appropriate type
+and flag.  The \code{comment} entry is a fixed string which is used
+for the comment associated with this metadata item.  The arguments to
+the \code{name} formatting codes and the metadata data itself are
+passed to \code{psMetadataItemAlloc} as arguments following the
+comment string.  The data must be a pointer for any data types which
+are stored in the element \code{data.void}, while other data types are
+passed as numeric values.  The argument list must be interpreted
+appropriately by the \code{va_list} operators in the function.
+\begin{verbatim}
+psMetadataItem *psMetadataItemAlloc(const char *name, int format, const char *comment, ...);
+psMetadataItem *psMetadataItemAllocV(const char *name, int format, const char *comment, va_list list);
+\end{verbatim}
+The \code{psMetadataItem} destructor is specified below.  Note that
+the destructor for \code{psMetadataItem} must call the appropriate
+destructor for the \code{val} (recall that it is the duty of the
+\code{psMyTypeFree}s to decrement the \code{refCounter} and free the
+memory if and only if the \code{refCounter == 1} --- see
+\S\ref{sec:free}).
+\begin{verbatim}
+void psMetadataItemFree(psMetadataItem *item);
+\end{verbatim}
+The constructor for the collection of metadata, \code{psMetadata},
+simply returns an empty metadata container (employing the constructors
+for the doubly-linked list and hash table).  The destructor needs to
+free each of the \code{psMetadataItem}s using \code{psMetadataItemFree}.
+\begin{verbatim}
+psMetadata *psMetadataAlloc(void);
+void psMetadataFree(psMetadata *md);
+\end{verbatim}
+Items may be added to the metadata in one of two ways --- firstly, an
+item may be added by appending a \code{psMetadataItem} which has
+already been created; and secondly by directly providing the data to
+be appended.  In both cases, the \code{psMetadataItem} that is
+appended to the metadata is returned.  The second function,
+\code{psMetadataAppend} takes a pointer or value which is interpretted
+by the \code{va_list} operators in the function.
+%
+\begin{verbatim}
+psMetadataItem *psMetadataAppendItem(psMetadata *restrict md,
+                                     psMetadataItem *restrict item);
+psMetadataItem *psMetadataAppend(psMetadata *restrict md, const char *name,
+                                 int format, const char *comment, ...);
+\end{verbatim}
+Items may be removed from the metadata by specifying a key.  If the
+key matches a metadata item, the item is removed from the metadata and
+returned; otherwise, \code{NULL} is returned.  If the key is not
+unique, then \emph{all} items corresponding to the key are removed,
+and the first item is returned.  Care should be taken not to leak
+memory when appending an item for which the key already exists in the
+metadata (and is not \code{PS_META_NON_UNIQUE}).
+%
+\begin{verbatim}
+psMetadataItem *psMetadataRemove(psMetadata *restrict md, const char *restrict key);
+\end{verbatim}
+The metadata may be iterated over by (re-)setting the iterator for the
+appropriate \code{psMetadata}, and getting the next item.
+\code{psMetadataGetNext} has the ability to match the beginning of a
+key, e.g., if the user only wants to iterate through
+\code{IPP.machines.sky} and doesn't want to bother with
+\code{IPP.machines.detector}.  The iterator should iterate over every
+item of metadata --- even those that are non-unique.
+\begin{verbatim}
+void psMetadataSetIterator(psMetadata *md);
+psMetadataItem *psMetadataGetNext(psMetadata *restrict md, const char *restrict match);
+\end{verbatim}
+Items may be found within the metadata by providing a key.  In the
+event that the key is non-unique, the first item is returned.
+\begin{verbatim}
+psMetadataItem *psMetadataLookup(const psMetadata *restrict md,
+                                 const char *restrict key);
+\end{verbatim}
+Metadata items may be printed to an open file descriptor, optionally
+pre-pending a specified string.
+\begin{verbatim}
+void psMetadataItemPrint(FILE *fd, const psMetadataItem *restrict md,
+                         const char *prefix);
+\end{verbatim}
+\begin{verbatim}
+psMetadata *psMetadataReadHeader(psMetadata *output, const char *extname,
+                              int extnum, const char *filename);
+\end{verbatim}
+Read header data from a FITS image file into a \code{psMetadata}
+structure (see section~\ref{sec:metadata}.  The \code{extname} and
+\code{extnum} parameters specify the extension of interest as above.
+If the named extension does not exist, the function should return an
+error.
+\begin{verbatim}
+psMetadata *psMetadataFReadHeader(psMetadata *output, const char *extname,
+                               int extnum, FILE *f);
+\end{verbatim}
+Read header data from a FITS image file descriptor into a
+\code{psMetadata} structure.
+\subsection{Detector and sky positions}
+Both detector and sky positions will be used extensively in the IPP.
+The first are linear coordinates which conform to Euclidean geometry
+while the second are angular coordinates for which additional care
+must often be taken.  We put these into two structures, \code{psPlane}
+and \code{psSphere}, respectively.  Partitioning these two will enable
+error-checking.
+%
+\begin{verbatim}
+typedef struct {
+    double x;                           ///< x position
+    double y;                           ///< y position
+    double xErr;                        ///< Error in x position
+    double yErr;                        ///< Error in y position
+} psPlane;
+typedef struct {
+    double r;                           ///< RA
+    double d;                           ///< Dec
+    double rErr;                        ///< Error in RA
+    double dErr;                        ///< Error in Dec
+} psSphere;
+\end{verbatim}
+Three major classes of coordinate transformations are necessary.
+First, linear coordinates from one frame must be converted to linear
+coordinates in a different frame of references.  Simple transformations
+of this type are independent of other quantities of the positions --
+they are simply mapping between two linear spaces.  In practice, these
+transformations may often be a function of the magnitude or color of
+the imaged object.  The second type of conversion is the
+transformation of linear coordinates to angular coordinates and
+vice-versa.  This conversion depends on the desired projection, and
+may represent the real mapping performed by the telescope or may
+simply represent a convenient mechanism to display 3D coordinates in
+useful forms.  The third conversion of interest is the transformation
+of one set of spherical coordinates to another set.  Frequently in
+astronomy, these conversions consist only of rotations between the two
+spherical coordinates systems, where the coordinates of the pole and
+equatorial rotation between the two systems define the
+transformation.  Conversions between standard coordinate systems such
+as Galactic, Ecliptic, and various epochs of the Celestial coordinates
+are represented by these spherical transformations.
+\subsubsection{Linear Coordinate Transformations}
+We specify two types of transforms between coordinate systems.  The
+first consists simply of two 2D polynomials to transform both
+components -- the output coordinates depend only on the input
+coordinates and no other quantities of objects at those coordinates.
+The second consists of two 4D polynomials in which the output
+coordinates are also specified to be a function of the magnitude and
+color of the object with the given coordinates.  This type of
+coordinate transformation is necessary to represent the
+(color-dependent) optical distortions caused by the atmosphere and
+camera optics, and the possibly effects of charge transfer
+inefficiency.  We specify two structures to represent the coefficients
+of these transformations:
+\begin{verbatim}
+typedef struct {
+    psDPolynomial2D *x;
+    psDPolynomial2D *y;
+} psPlaneTransform;
+\end{verbatim}
+The \code{psDPolynomial2D} structures represent polynomials of
+arbitrary order as a function of two dimensions.  There is one of
+these structures for each of the two output dimensions.  As an
+example, consider the simple transformation from one linear coordinate
+frame ($x,y$), e.g., on a CCD, to a second frame ($p,q$), e.g., on a
+chip. If we have only first order terms in the transformation
+\code{psPlaneTransform T}, the new coordinates would be represented by
+the terms:
+%
+\begin{verbatim}
+p = T.x->coeff[0][0] + x*T.x->coeff[1][0] + y*T.x->coeff[0][1];
+q = T.y->coeff[0][0] + x*T.y->coeff[1][0] + y*T.y->coeff[0][1];
+\end{verbatim}
+%
+where we have excluded the basic cross-term ($x \times y$) by using
+the mask: \code{T.x->mask[1][1] = 0; T.y->mask[1][1] = 0;}
+The \code{psPlaneDistortion} represents an optical distortion.  The
+lowest two terms are the $x$ and $y$ axis of the target system.  The
+higher two terms may represent magnitude and color terms.
+\begin{verbatim}
+typedef struct {
+    psDPolynomial4D *x;
+    psDPolynomial4D *y;
+} psPlaneDistort;
+\end{verbatim}
+Like \code{psPlaneTransform}, \code{psPlaneDistort} contains two
+\code{psDPolynomial4D} structures representing polynomials of
+arbitrary order as a function of four, rather than two dimensions.
+There is one of these structures for each of the two output
+dimensions.  In this structure, the highest two dimensions could
+represent a magnitude and a color.  As an example, consider the simple
+transformation from one linear coordinate frame ($x,y$), e.g., on a
+CCD, of an object with magnitude and color ($m,c$) to a second frame
+($p,q$), e.g., the focal plane. If we have only first order terms in
+the transformation \code{psPlaneTransform T}, the new coordinates
+would be represented by the terms:
+%
+\begin{verbatim}
+p = T.x->coeff[0][0][0][0] + x*T.x->coeff[1][0][0][0] + y*T.x->coeff[0][1][0][0]
+  + m*T.x->coeff[0][0][1][0] + c*T.x->coeff[0][0][0][1]
+q = T.y->coeff[0][0][0][0] + x*T.y->coeff[1][0][0][0] + y*T.y->coeff[0][1][0][0]
+  + m*T.y->coeff[0][0][1][0] + c*T.y->coeff[0][0][0][1]
+\end{verbatim}
+%
+where we have again excluded the cross-term ($x \times y$) by using the
+mask.
+We require corresponding functions to apply the transformations to a
+specified coordinate \code{coords}:
+%
+\begin{verbatim}
+psPlane *psPlaneTransformApply (psPlane *out,
+                                const psPlaneTransform *transform,
+                                const psPlane *coords);
+psPlane *psPlaneDistortApply (psPlane *out,
+                              const psPlaneDistort *distort,
+                              const psPlane *coords,
+                              float term3, float term4);
+\end{verbatim}
+\subsubsection{Celestial Coordinate Conversions}
+We need to be able to convert between ICRS, Galactic and Ecliptic
+coordinates, and potentially between arbitrary spherical coordinate
+systems.  All of these basic spherical transformations represent
+rotations of the spherical coordinate reference.  We specify a general
+transformation function which takes a structure,
+\code{psSphereTransform}, defining the transformation between two
+spherical coordinate systems (the structure contains the sines and
+cosines of the angles involved so as to minimize computation time for
+repeated transformations).  We also define a function to generate
+\code{psSphereTransform}, based on the three angles
+describing the location of the pole and the relative equatorial
+rotations of the two systems.  We also specify special functions to
+return the \code{psSphereTransform} for transformations
+between standard coordinate systems.
+\begin{verbatim}
+typedef struct {
+    double sinNPlon;                     ///< sin of North Pole longitude
+    double cosNPlon;                     ///< cos of North Pole longitude
+    double sinNPlat;                     ///< sin of North Pole lattitude
+    double cosNPlat;                     ///< cos of North Pole lattitude
+    double sinZP;                        ///< sin of First PT of Ares lon
+    double cosZP;                        ///< cos of First PT of Ares lon
+} psSphereTransform;
+\end{verbatim}
+The constructor and destructor are defined as follows:
+\begin{verbatim}
+psSphereTransform *psSphereTransformAlloc(double NPlon, double NPlat, double ZP);
+void psSphereTransformFree(psSphereTransform *trans);
+\end{verbatim}
+where \code{NPlon} and \code{NPlat} define the coordinates in the
+input system of the north pole in the output system and \code{ZP}
+defines the longitude in the input system of the equatorial
+intersection between the two systems (e.g, the first point of Ares).
+Spherical coordinates may be transformed by providing the
+transformation and the coordinate in the input system to
+\code{psSphereTransform}:
+\begin{verbatim}
+psSphere *psSphereTransformApply(psSphere *out,
+                                 const psSphereTransform *transform,
+                                 const psSphere *coord);
+\end{verbatim}
+The following functions simply return the appropriate
+\code{psSphereTransform} to convert between predefined spherical
+coordinate systems (i.e., {\bf I}CRS, {\bf E}cliptic and {\bf
+G}alactic).
+%
+\begin{verbatim}
+psSphereTransform *psSphereTransformItoE(void);
+psSphereTransform *psSphereTransformEtoI(void);
+psSphereTransform *psSphereTransformItoG(void);
+psSphereTransform *psSphereTransformGtoI(void);
+\end{verbatim}
+\subsubsection{Projections}
+We require functions to convert between spherical and linear
+coordinate systems based on a variety of projections.  The required
+projections include:
+\begin{itemize}
+\item TAN
+\item SIN
+\item AIT
+\item PAR
+\item GLS
+\end{itemize}
+We specify the following structure \code{psProjection} to define the
+parameters of the projection:
+\begin{verbatim}
+typedef struct {
+    double R, D;                         ///< coordinates of projection center
+    double Xs, Ys;                       ///< plate-scale in X and Y directions
+    psProjectionType type;               ///< projection type
+} psProjection;
+\end{verbatim}
+The projection type is defined by the following enumerated type \code{psProjectionType}:
+\begin{verbatim}
+typedef enum {                          ///< type of val is:
+    PS_PROJ_TAN,                        ///< Tangent projection
+    PS_PROJ_SIN,                        ///< Sine projection
+    PS_PROJ_AIT,                        ///< Aitoff projection
+    PS_PROJ_PAR,                        ///< Par projection
+    PS_PROJ_GLS,                        ///< GLS projection
+    PS_PROJ_NTYPE                       ///< Number of types; must be last
+} psProjectionType;
+\end{verbatim}
+The following functions will project and deproject (respectively)
+spherical coordinates:
+\begin{verbatim}
+psPlane  *psProject(const psSphere *coord, const psProjection *projection);
+psSphere *psDeproject(const psPlane *coord, const psProjection *projection);
+\end{verbatim}
+\subsubsection{Offsets}
+We require a function to calculate the offset between two positions on
+the sky, as well as a function to apply an offset to a position.  The
+first determines the offset (RA,Dec) on the sky between two positions.
+The second applies the given offset to the coordinate. The offsets may
+be of type: \tbd{Linear, Spherical/Arcsec, Spherical/Degreees, etc.}
+\begin{verbatim}
+psSphere *psSphereGetOffset(const psSphere *restrict position1,
+                            const psSphere *restrict position2,
+                            const char *type);
+psSphere *psSphereSetOffset(const psSphere *restrict position,
+                            const psSphere *restrict offset,
+                            const char *type);
+\end{verbatim}
+Note that these should propagate the errors appropriately.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Astronomical Images}
+\subsubsection{Overview}
+Above, we have defined a basic container for a single 2D collection of
+pixels (\code{psImage}), along with basic operations to manipulate the
+image pixels.  For astronomical applications, this data structure is
+insufficient for two reasons.  First, it does provide sufficient
+additional metadata to describe the data in detail.  Second, astronomy
+applications frequent involve multiple, related images.  For
+PanSTARRS, and for general astronomical applications, we require a
+richer collection of data structures which describe a very general
+image concept.  We have defined several layers in the hierarchy which
+are necessary to describe the image data which will be produced by the
+PanSTARRS Gigapixel cameras as well as other standard astronomical
+images.
+A simple 2D image is a basic data unit for much of astronomical
+imaging.  If we consider various optical and IR array cameras, a
+single readout of the detector produces a collection of pixels
+measurements.  We define our lowest-level astronomical image
+structure, \code{psReadout}, to contain the pixels produced by a
+single readout of the detector, along with metadata needed to define
+that readout: the origin and binning of the image relative to the
+original detector pixels explicitly in the structure, and pointers to
+the general metadata and derived objects, if any.
+A single detector may produce more than one read which is associated.
+For example, infrared detectors frequently produce an image
+immediately after the detector is reset followed by an image after the
+basic exposure is complete.  Both readouts correspond to the same
+pixels, though the binning or rastering may be different between the
+two readouts.  Another example is the video sequence produced by the
+PanSTARRS Gigapix camera guide cells, each of which represents a
+series of many images from a subraster of pixels in the detector
+readout portion.  The second level of our image container hierarchy,
+\code{psCell}, consists of a collection of readouts from a single
+detector.
+In the PanSTARRS Gigapix camera, the basic readout region is a
+fraction of the full imaging area of a single CCD chip.  The chip is
+divided into 64 cells, any fraction of which may have been readout
+for a given exposure.  In other cameras, such as Megacam at CFHT, the
+individual CCDs have multiple amplifiers addressing contiguous
+portions of the detector.  In such cameras, each amplifier produces a
+separate collection of pixels.  In the third level of our image
+container hierarchy, the data structure \code{psChip} represents a
+collection of different cells.
+The top level of our image container hierarchy is a complete focal
+plane array (\code{psFPA}).  This structure represents the collection
+of chips in the camera, all of which are read out in a given
+exposure.
+For example, take a mosaic camera consisting of eight $2k\times 4k$
+CCDs, each of which is read out through two amplifiers.  Then there
+would be sixteen cells in total, each of which is presumably $2k\times
+k$.  There would be eight chips, each consisting of two cells, and
+the focal plane consists of these eight chips.
+As another example, consider an observation by PS1.  The focal plane
+would consist of 60 chips, each of which consist of 64 cells (or less;
+a few cells may be dead).  Some cells (those containing guide stars
+for the orthogonal transfer) will contain multiple readouts.
+These data structures represent containers with which to carry around
+the collection of related image data.  There is no requirement on the
+functions or the structures that each instance of one of these data
+structures represent the physical hardware.  For example, it is not
+necessary that an instance of \code{psFPA} always carry the data for
+all 60 (or 64) Gigapixel camera OTAs.  The usage of these structures
+is such that all astronomical operations which apply to a CCD image
+should be performed on an instance of \code{psFPA}.  If a particular
+circumstance only requires a single 2D image, then that is represented
+by an instance of \code{psFPA} with one \code{psChip}, which in turn
+has one \code{psCell}, which in turn has one \code{psReadout}.
+These container levels also include in their definition the information
+needed to transform the coordinates in one of the levels to the
+coordinate system relevant at the higher levels.
+\subsubsection{A Readout}
+A readout is the result of a single read of a cell (or a portion
+thereof).  It contains a pointer to the pixel data, and additional
+pointers to the objects found in the readout, and the readout
+metadata.  It also contains the offset from the lower-left corner of
+the chip, in the case that the CCD was windowed.
+\begin{verbatim}
+typedef struct {
+    const int x0, y0;                ///< Offset from the lower-left corner
+    const int nx, ny;                ///< Image binning
+    psImage *image;                  ///< imaging area of cell
+    psDlist *objects;                ///< objects derived from cell
+    psMetadata *md;                  ///< Readout-level metadata
+} psReadout;
+\end{verbatim}
+\subsubsection{A Cell}
+A cell consists of one or more readouts (usually only one except in the
+case that the cell has been used for fast guiding).  It also contains
+a pointer to the cell metadata, and a pointer to its parent chip.  On
+the astrometry side, it also contains coordinate transforms from the
+cell to the chip and, as a convenience, from the cell to the focal
+plane.  It is expected that these transforms will consist of two
+first-order 2D polynomials, simply specifying a translation, rotation
+and magnification; hence they are easily inverted, and there is no
+need to add reverse transformations.  We also add an additional
+transformation, which is intended to provide a ``quick and dirty''
+transform from the cell coordinates to the sky; this transformation
+not guaranteed to be as precise as the ``standard'' transformation of
+Cell $\rightarrow$ Chip $\rightarrow$ Focal Plane $\rightarrow$
+Tangent Plane $\rightarrow$ Sky, but will be faster.
+\begin{verbatim}
+typedef struct {
+    int nReadouts;                      ///< number of readouts in this cell
+    struct psReadout *readouts;         ///< Readouts from the cell
+    psMetadata *md;                     ///< Cell-level metadata
+    psPlaneTransform *cellToChip;       ///< Transformations from cell to chip coords
+    psPlaneTransform *cellToFPA;        ///< Transformations from cell to FPA coords
+    psPlaneTransform *cellToSky;        ///< Direct from cell to sky coords
+    struct psChip  *parentChip;         ///< chip which contains this cell
+} psCell;
+\end{verbatim}
+\subsubsection{A Chip}
+A chip consists of one or more cells (according to the number of
+amplifiers on the CCD).  It contains a pointer to the chip metadata,
+and a pointer to the parent focal plane.  For astrometry, it contains
+a coordinate transform from the chip to the focal plane.  It is
+expected that this transforms will consist of two second-order 2D
+polynomials; hence we think that it is prudent to include a reverse
+transformation which will be derived from numerically inverting the
+forward transformation.
+\begin{verbatim}
+typedef struct {
+    int nCells;                         ///< Number of Cells assigned
+    struct psCell *cells;               ///< Cells in the Chip
+    psMetadata *md;                     ///< Chip-level metadata
+    psPlaneTransform *chipToFPA;        ///< Transformations from chip to FPA coords
+    psPlaneTransform *FPAtoChip;        ///< Transformations from FPA to chip
+    struct psFPA *parentFPA;            ///< FPA which contains this chip
+} psChip;
+\end{verbatim}
+\subsubsection{A Focal Plane}
+A focal plane consists of one or more chips (according to the number
+of pieces of contiguous silicon).  It contains pointers to the focal
+plane metadata and the exposure information.  For astrometry, it
+contains a transformation from the focal plane to the tangent plane
+and the fixed pattern residuals.  It is expected that the
+transformation will consist of two 4D polynomials (i.e.\ a function of
+two coordinates in position, the magnitude of the object, and the
+color of the object) in order to correct for optical distortions and
+the effects of the atmosphere; hence we think that it is prudent to
+include a reverse transformation which will be derived from
+numerically inverting the forward transformation.  Since colors are
+involved in the transformation, it is necessary to specify the color
+the transformation is defined for.  We also include some values to
+characterize the quality of the transformation: the root mean square
+deviation for the x and y transformation fits, and the $\chi^2$ for
+the transformation fit.
+\begin{verbatim}
+typedef struct {
+    int nChips;                         ///< Number of Cells assigned
+    int nAlloc;                         ///< Number of Cells available
+    struct psChip *chips;               ///< Chips in the Focal Plane Array
+    psMetadata *md;                     ///< FPA-level metadata
+    psPlaneDistort *TPtoFP;             ///< Transformation term from
+    psPlaneDistort *FPtoTP;             ///< Transformation term from
+    psFixedPattern *pattern;            ///< Fixed pattern residual offsets
+    const psExposure *exp;              ///< information about this exposure
+    psPhotSystem colorPlus, colorMinus; ///< Colour reference
+    float rmsX, rmsY;                   ///< Dispersion in astrometric solution
+    float chi2;                         ///< chi^2 of astrometric solution
+} psFPA;
+\end{verbatim}
+\subsubsection{Exposure information}
+We need several quantities from the telescope in order to make a
+first guess at the astrometric solution.  From these quantities,
+further quantities can be derived and stored for later use.
+\begin{verbatim}
+typedef struct {
+    const double ra, dec;               ///< Telescope boresight
+    const double ha;                    ///< Hour angle
+    const double zd;                    ///< Zenith distance
+    const double az;                    ///< Azimuth
+    const double lst;                   ///< Local Sidereal Time
+    const float mjd;                    ///< MJD of observation
+    const float rotAngle;               ///< Rotator position angle
+    const float temp;                   ///< Air temperature, for estimating refraction
+    const float pressure;               ///< Air pressure, for calculating refraction
+    const float humidity;               ///< Relative humidity, for refraction
+    const float exptime;                ///< Exposure time
+    /* Derived quantities */
+    const float posAngle;               ///< Position angle
+    const float parallactic;            ///< Parallactic angle
+    const float airmass;                ///< Airmass, calculated from zenith distance
+    const float pf;                     ///< Parallactic factor
+    const char *cameraName;             ///< name of camera which provided exposure
+    const char *telescopeName;          ///< name of telescope which provided exposure
+} psExposure;
+\end{verbatim}
+\subsubsection{Constructors and Destructors}
+Each of the above structures needs an appropriate constructor and
+destructor.  Other than \code{psExposure}, which contains significant
+non-pointer types, the constructors should not take any arguments, and
+the destructors should only take the structure to be destroyed.
+The constructor for \code{psExposure} is specified below.
+\begin{verbatim}
+psExposure *
+psExposureAlloc(double ra, double dec,  ///< Telescope boresight
+                double ha,              ///< Hour angle
+                double zd,              ///< Zenith distance
+                double az,              ///< Azimuth
+                double lst,             ///< Local Sidereal Time
+                float mjd,              ///< MJD
+                float rotAngle,         ///< Rotator position angle
+                float temp,             ///< Temperature
+                float pressure,         ///< Pressure
+                float humidity,         ///< Relative humidity
+                float exptime);         ///< Exposure time
+\end{verbatim}
+\subsection{Astrometry}
+Astrometry is a basic functionality required for the IPP that will be
+used repeatedly, both for low-precision (roughly where is my favorite
+object?) and high-precision (what is the proper motion of this star?).
+As such, it must be flexible, yet robust.  Accordingly, we will wrap
+the StarLink Astronomy Libraries (SLALib), which has already been
+developed.
+\subsubsection{Coordinate frames}
+\label{sec:coordinateFrames}
+There are five coordinate frames that we need to worry about for the
+purposes of astrometry:
+\begin{itemize}
+\item Cell: $(x,y)$ in pixels --- raw coordinates;
+\item Chip: $(X,Y)$ in pixels --- the location on the silicon;
+\item Focal Plane: $(p,q)$ in microns --- the location on the focal plane;
+\item Tangent Plane: $(l,m)$ in arcsec from the telescope boresight; and
+\item Sky: (RA,Dec) --- ICRS.
+\end{itemize}
+The following steps are required to convert from the cell coordinates to
+the sky:
+\begin{itemize}
+\item Cell $\longleftrightarrow$ Chip: two 2D polynomials, $(X,Y) = f(x,y)$;
+\item Chip $\longleftrightarrow$ FP: two 2D polynomials, $(p,q) = g(X,Y)$;
+\item FP $\longleftrightarrow$ TP: two 4D polynomials, $(l,m) =
+h(p,q,m,c)$, where $m$ and $c$ are the magnitude and color of the
+object, respectively; and
+\item TP $\longleftrightarrow$ Sky: SLALib transformation using a
+transform pre-computed for each pointing.
+\end{itemize}
+Note that the transformation between the Focal Plane and the Tangent
+Plane is a four-dimensional polynomial, in order to account for any
+possible dependencies in the astrometry on the stellar magnitude and
+color; the former serves as a check for charge transfer
+inefficiencies, while the latter will correct chromatic refraction,
+both through the atmosphere and the corrector lenses.
+We require structures to contain each of the above transformations as
+well as the pixel data.
+\subsubsection{SLALib information}
+SLALib requires several elements to perform the transformations
+between the tangent plane and the sky.  Pre-computing these quantities
+for each exposure means that subsequent transformations are faster.
+For historical reasons, this structure is known colloquially as
+``Wallace's Grommit''.
+\begin{verbatim}
+typedef struct {
+    const double latitude;              ///< geodetic latitude (radians)
+    const double sinLat, cosLat;        ///< sine and cosine of geodetic latitude
+    const double abberationMag;         ///< magnitude of diurnal aberration vector
+    const double height;                ///< height (HM)
+    const double temperature;           ///< ambient temperature (TDK)
+    const double pressure;              ///< pressure (PMB)
+    const double humidity;              ///< relative humidity (RH)
+    const double wavelength;            ///< wavelength (WL)
+    const double lapseRate;             ///< lapse rate (TLR)
+    const double refractA, refractB;    ///< refraction constants A and B (radians)
+    const double longitudeOffset;       ///< longitude + ... (radians)
+    const double siderealTime;          ///< local apparent sidereal time (radians)
+} psGrommit;
+\end{verbatim}
+The \code{psGrommit} is calculated from telescope information for the
+particular exposure:
+\begin{verbatim}
+psGrommit *psGrommitAlloc(const psExposure *exp);
+void psGrommitFree(psGrommit *grommit);
+\end{verbatim}
+\subsubsection{Fixed Pattern}
+The fixed pattern is a correction to the general astrometric solution
+formed by summing the residuals from many observations.  The intent is
+to correct for higher-order distortions in the camera system on a
+coarse grid (larger than individual pixels, but smaller than a single
+cell).  Hence, in addition to the offsets, we need to specify the size
+and scale of the grid in $x$ and $y$, as well as the origin of the
+grid.
+\begin{verbatim}
+typedef struct {
+    int nX, nY;                         ///< Number of elements in x and y
+    double x0, y0;                      ///< Position of 0,0 corner on focal plane
+    double xScale, yScale;              ///< Scale of the grid
+    double **x, **y;                    ///< The grid of offsets in x and y
+} psFixedPattern;
+\end{verbatim}
+\subsubsection{Position Finding}
+We require functions to return the structure containing given
+coordinates.  For example, we want the chip that corresponds to the
+focal plane coordinates $(p,q) = (-1.234,+5.678)$.  These routines
+handle the one-to-many problem --- i.e., for one given focal plane
+coordinate, there are many chips that this coordinate may be
+correspond to; these functions will select the correct one.
+%
+\begin{verbatim}
+psCell *psCellInFPA (psCell *out, const psPlane *coord, const psFPA *fpa);
+psChip *psChipInFPA (psChip *out, const psPlane *coord, const psFPA *fpa);
+psCell *psCellInChip(psCell *out, const psPlane *coord, const psChip *chip);
+\end{verbatim}
+\subsubsection{Conversion Functions}
+We require functions to convert between the various coordinate frames
+(Section~\ref{sec:coordinateFrames}).  The hierarchy of the coordinate
+frames and the transformations between each are shown in
+Figure~\ref{fig:coco}.  The functions that employ the transformations
+are shown in Figure~\ref{fig:cocoFunc}.  In addition to
+transformations between each adjoining coordinate frame in the
+hierarchy, we also require higher-level functions to convert between
+the Cell and Sky coordinate frames; these will simply perform the
+intermediate steps.
+\begin{figure}
+\psfig{file=coordinateFrames,height=7in,angle=-90}
+\caption{The coordinate systems in the \PS{} IPP, and the relation
+between each by transformations contained in the appropriate
+structures.}
+\label{fig:coco}
+\end{figure}
+\begin{figure}
+\psfig{file=coordinateConv,height=7in,angle=-90}
+\caption{Conversion between coordinate systems by PSLib.}
+\label{fig:cocoFunc}
+\end{figure}
+We specify the following functions to convert between coordinates in
+one type of frame to another type of frame.  The first group consist
+of unambiguous transformations: from the coordinates in a low-level
+frame to the coordinates in the containing higher-level frame, of
+which only one exists.  In all of these functions, the output
+coordinate structure may be \code{NULL} or may be supplied by the
+calling function.  In the former case, the structure must be
+allocated; in the latter case, the supplied structure must be used.
+\begin{verbatim}
+psPlane *psCoordCelltoChip (psPlane *out, const psPlane *in, const psCell *cell);
+\end{verbatim}
+which converts coordindates \code{in} on the specified \code{cell} to
+the coordinates on the parent chip.
+\begin{verbatim}
+psPlane *psCoordChiptoFPA (psPlane *out, const psPlane *in, const psChip *chip);
+\end{verbatim}
+which converts the coordinates \code{in} on the specified \code{chip}
+to the coordinates on the parent FPA.
+\begin{verbatim}
+psPlane *psCoordFPAToTP(psPlane *out, const psPlane *in, const psFPA *fpa);
+\end{verbatim}
+which converts coordinates \code{in} on the specified focal plane
+\code{fpa} to tangent plane coordinates, applying the appropriate
+distortion terms.
+\begin{verbatim}
+psSphere *psCoordTPtoSky(psSphere *out, const psPlane *in, const psGrommit *grommit);
+\end{verbatim}
+which converts the tangent plane coordinates \code{in} to (RA,Dec) on
+the sky, based on the environmental information specified by
+\code{grommit}.
+\begin{verbatim}
+psPlane *psCoordCellToFPA(psPlane *out, const psPlane *in, const psCell *cell);
+\end{verbatim}
+which performs the single-step conversion between Cell coordinates
+\code{in} and FPA coordinates.
+\begin{verbatim}
+psSphere *psCoordCellToSky(psSphere *out, const psPlane *in, const psCell *cell);
+\end{verbatim}
+which converts coordinates on the specified cell to (RA,Dec).  This
+transformation must be performed using the intermediate stage
+transformations of Cell to Chip, Chip to FPA, FPA to Tangent Plane,
+Tangent Plane to Sky.  The information needed for each of these
+transformations is available in the \code{.parent} elements of
+\code{psCell} and \code{psChip}, and the \code{psFPA.exposure}
+element.
+\begin{verbatim}
+psSphere *psCoordCellToSkyQD(psSphere *out, const psPlane *in, const psCell *cell);
+\end{verbatim}
+which uses the 'quick-and-dirty' transformation to convert coordinates
+on the specified cell to (RA,Dec).  This transformation should use the
+locally linear transformation specified by the element
+\code{psCell.cellToSky}.  Although the accuracy of this transformation
+is lower than the complete transformation above, the calculation is
+substantially faster as it only involves linear transformations.
+The following functions convert from high-level frames to the
+coordinates of contained lower-level frames.
+\begin{verbatim}
+psPlane *psCoordSkyToTP(psPlane *out, const psSphere *in, const psGrommit *grommit);
+\end{verbatim}
+which converts (RA,Dec) coordinates \code{in} to tangent plane coords
+based on the enviromental information supplied by \code{grommit}.
+\begin{verbatim}
+psPlane *psCoordTPtoFPA(psPlane *out, const psPlane *in, const psFPA *fpa);
+\end{verbatim}
+which converts the tangent plane coordinates \code{in} to focal plane coordinates.
+\begin{verbatim}
+psPlane *psCoordFPAtoChip (psPlane *out, const psPlane *in, const psChip *chip);
+\end{verbatim}
+which converts the specified FPA coordinates \code{in} to the
+coordinates on the given Chip.  The specified chip need not contain
+the input coordinate.  To find the chip which contains a particular
+coordinate, the function \code{psChipInFPA}, defined above, should be
+used.
+\begin{verbatim}
+psPlane *psCoordChiptoCell (psPlane *out, const psPlane *in, const psCell *cell);
+\end{verbatim}
+which converts the specified Chip coordinate \code{in} to the
+coordinate on the given Cell.  The specified Cell need not contain the
+input coordinate.  To find the cell which contains a particular
+coordinate, the function \code{psCellInChip}, defined above, should be
+used.
+\begin{verbatim}
+psPlane *psCoordSkyToCell(psPlane *out, const psSphere *in, psCell *cell);
+\end{verbatim}
+which directly converts (RA,Dec) \code{in} to coordinates on the
+specified cell.  The specified cell need not contain the input
+coordinates.
+\begin{verbatim}
+psPlane *psCoordSkyToCellQD(psPlane *out, const psSphere *in, psCell *cell);
+\end{verbatim}
+which directly converts (RA,Dec) \code{in} to coordinates on the
+specified cell.  The specified cell need not contain the input
+coordinates.  This transformation should use the locally linear
+transformation specified by the element \code{psCell.cellToSky}.
+Although the accuracy of this transformation is lower than the
+complete transformation above, the calculation is substantially faster
+as it only involves linear transformations.
+\subsubsection{Additional functions}
+We require additional functions to perform general functions which
+will be useful for astrometry.  Given coordinates on the sky, we
+need to get the airmass, the parallactic angle, and an estimate of
+the atmospheric refraction.
+\begin{verbatim}
+float psGetAirmass(const psSphere *coord, double siderealTime, float height);
+\end{verbatim}
+which returns the airmass for a given position and sidereal time.
+\begin{verbatim}
+float psGetParallactic(const psSphere *coord, double siderealTime);
+\end{verbatim}
+which returns the parallactic angle for a given position and sidereal time.
+\begin{verbatim}
+float psGetRefraction(float colour,            ///< Colour of object
+                      psPhotSystem colorPlus,  ///< Colour reference
+                      psPhotSystem colorMinus, ///< Colour reference
+                      const psExposure *exp);  ///< Telescope pointing information
+\end{verbatim}
+which provides an estimate of the atmospheric refraction, along the parallactic angle.
+\begin{verbatim}
+double psGetParallaxFactor(const psExposure *exp)
+\end{verbatim}
+Calculate the parallax factor for the given exposure \tbd{why do we
+  need this?}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Photometry}
+Photometric observations are performed in an instrumental photometric
+system, and must be related to other photometric systems.  We
+require a data structure which defines a photometric system, as well
+as a structure to define the transformation between photometric
+systems.
+The photometric system is defined by the psPhotSystem structure.
+A photometric system is identified by a human-readable \code{name}
+(ie, SDSS.g, Landolt92.B, GPC1.OTA32.r).  Each photometric system is
+given a unique identifier \code{ID}.  Observations taken with a
+specific camera, detector, and filter represent their own photometric
+system, and it may be necessary to perform transformations between
+these systems.  Photometric systems associated with observations from
+a specific camera/detector/filter combination can be associated with
+those components.
+\begin{verbatim}
+typedef struct {
+    const int ID;                       ///< ID number for this photometric system
+    const char *name;                   ///< Name of photometric system
+    const char *camera;                 ///< Camera for photometric system
+    const char *filter;                 ///< Filter used for photometric system
+    const char *detector;               ///< Detector used for photometric system
+} psPhotSystem;
+\end{verbatim}
+The following structure defines the transformation between two
+photometric systems.
+\begin{verbatim}
+typedef struct {
+    psPhotSystem src;                   ///< Source photometric system
+    psPhotSystem dst;                   ///< Destination photometric system
+    psPhotSystem pP, pM;                ///< Primary color reference
+    psPhotSystem sP, sM;                ///< Secondary color reference
+    float pA, sA;                       ///< Color offset for references
+    psPolynomial3D transform;           ///< Transformation from source to destination
+} psPhotTransform;
+\end{verbatim}
+The transformation between two photometric systems may depend on the
+airmass of the observation and on the colors of the object of
+interest.  For a specific observation, such a transformations can be
+defined as a polynomial function of the color of the star and the
+airmass of the observations.  If sufficient data exists, the
+transformation between the photometric systems may include more than
+one color, constraining the curvature of the stellar spectral energy
+distributions.  This latter term may be significant for stars which
+are highly reddened, for example.  Derived photometric quantities may
+have been corrected for airmass variations, in which case only color
+terms may be measurable.  The structure defines the transformation
+between a source photometric system (\code{src}) and a target
+photometric system (\code{dst}).  The photometric system of a primary
+color is defined by \code{pP, pM} such that the color is constructed
+as $pP - pM$.  A secondary color is defined by \code{sP, sM}.  For
+both, a reference color is specified (\code{pA, sA}): the polynomial
+transformation terms refer to colors in the form $pP - pM - pA$.  The
+transformation is specified as a 3D polynomial.  For a star of
+magnitude $M_{\rm src}$ in the source photometric system, with
+additional magnitude information in the other systems $M_{\rm pP}$,
+$M_{\rm pM}$, $M_{\rm sP}$, $M_{\rm sM}$, observed at an airmass of
+$z$, the magnitude of the star in the target system $M_{\rm dst}$ is
+given by: $M_{\rm dst} = M_{\rm src} + transform(z, M_{\rm pP} -
+M_{\rm pM} - pA, M_{\rm sP} - M_{\rm sM} - sA)$.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Astronomical objects}
+\subsubsection{Positions of Major SS Objects}
+We require the ability to calculate the position of major Solar System
+objects, as well as Lunar phase.  These functions all take the
+specified modified julian date, \code{mjd}, and the latitude/longitude
+of the observer.
+\begin{verbatim}
+psSphere *psSunGetPos(psTime time);
+psTime *psSunGetRise (psTime *twi15, psTime *twi18, psTime time);
+psTime *psSunGetSet (psTime *twi15, psTime *twi18, psTime time);
+psSphere *psMoonGetPos(psTime time, psSphere location);
+psTime *psMoonGetRise (psTime *twi15, psTime *twi18, psTime time);
+psTime *psMoonGetSet (psTime *twi15, psTime *twi18, psTime time);
+float psGetMoonPhase(psTime time);
+psSphere *psPlanetGetPos(psTime time, psSphere location);
+psTime *psPlanetGetRise (psTime *twi15, psTime *twi18, psTime time);
+psTime *psPlanetGetSet (psTime *twi15, psTime *twi18, psTime time);
+\end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\appendix
+\pagebreak
+\section{API Summary: all functions}
+\subsection{System Utilities}
+\input{psSystemGroup.tex}
+\subsection{Data Containers}
+\input{psDataGroup.tex}
+\subsection{Math Utilities}
+\input{psMathGroup.tex}
+\subsection{Astronomy Functions}
+\input{psAstroGroup.tex}
+\pagebreak
+\section{API Summary: all structures}
+\input{psStructures.tex}
+\section{Revision Change Log}
+\input{ChangeLogSDRS.tex}
+\bibliographystyle{plain} \bibliography{panstarrs}
+\end{document}
 An example of the usage of the metadata APIs is as follows:
 …
 \null\qquad\qquad\code{IPP.phase1.ota12.biassec}.
 We shall set in place a system for assigning the top-level `domains'
+We must set in place a system for assigning the top-level `domains'
 to responsible individuals, and for gathering a complete list of
 all metadata names in use throughout the project.
 …
+    }
 \end{verbatim}
-It is an unfortunate fact that certain metadata keywords (such as
-\code{COMMENT} and \code{HISTORY} in a FITS header) may be repeated
-with different values.  The \code{psMetadataAppend} routine is
-required to check that all metadata names are unique unless the type
-is qualified as \code{PS_META_NON_UNIQUE}; in this case a unique
-integer will be added to each name that you specify. In this case, you
-may either delete individual element separately or as a complete set:
-\begin{verbatim}
-psMetadataItemFree(psMetadataRemove(ms, "lang.hello.0"));
-psMetadataItemFree(psMetadataRemove(ms, "lang.hello"));
-\end{verbatim}
-\subsubsection{Metadata APIs}
-In this section, we explain the metadata APIs more fully.
-The allocator for \code{psMetadataItem} returns a full
-\code{psMetadataItem} ready for insertion into the
-\code{psMetadata}.  It is important to note that, in order to
-support multiple types, the data being input must be a pointer, even
-if it is a \code{float} or \code{int}, for example.  The name of the
-\code{psMetadataItem} takes a \code{sprintf} format string, with the
-corresponding arguments following.
-Note that the destructor for \code{psMetadataItem} must call the
-appropriate destructor for the \code{val} (recall that it is the duty
-of the \code{psMyTypeFree}s to decrement the \code{refCounter} and
-free the memory if and only if the \code{refCounter == 1} --- see
-\S\ref{sec:free}).
-\begin{verbatim}
-/** Constructor */
-psMetadataItem *psMetadataItemAlloc(int typeFlags, ///< type of this piece of metadata + flags
-                                    const void *val, ///< value of new item N.b. a pointer even if the item
-                                                     ///< is of type e.g. int
-                                    const char *comment, ///< comment associated with item
-                                    const char *name, ///< name of new item of metadata (may be in sprintf
-                                                      ///< format)
-                                    ... ///< possible arguments for name format
-    );
-/** Destructor */
-void psMetadataItemFree(psMetadataItem *ms ///< piece of metadata to destroy
-    );
-\end{verbatim}
-The constructor for the collection of metadata, \code{psMetadata},
-simply returns an empty metadata container (employing the constructors
-for the doubly-linked list and hash table).  The destructor needs to
-free each of the \code{psMetadataItem}s using
-\code{psMetadataItemFree}.
-\begin{verbatim}
-/** Constructor */
-psMetadata *psMetadataAlloc(void);   ///< make a new set of metadata
-/** Destructor */
-void psMetadataFree(psMetadata *md ///< destroy a set of metadata
-    );
-\end{verbatim}
-Items may be added to the metadata in one of two ways --- firstly, an
-item may be added by appending a \code{psMetadataItem} which has
-already been created; and secondly by directly providing the data to
-be appended.  In both cases, the \code{psMetadataItem} that is
-appended to the metadata is returned.
-\begin{verbatim}
-/// Add item to the end of the metadata
-psMetadataItem *psMetadataAppendItem(psMetadata *restrict md, ///< metadata to add to
-                                     psMetadataItem *restrict item ///< Metatdata to add
-    );
-/// Add item to the end of the metadata.  Combines psMetadataItemAlloc and psMetadataAppendItem
-psMetadataItem *psMetadataAppend(psMetadata *restrict md, ///< Metadata to add to
-                                 int typeFlags ///< type of this piece of metadata + flags
-                                 const void *val, ///< value of new item N.b. a pointer even if the item
-                                                  ///< is of type e.g. int
-                                 const char *comment, ///< comment associated with item
-                                 const char *name, ///< name of new item of metadata (may be in sprintf
-                                                   ///< format)
-                                 ...    ///< possible arguments for name format
-    );
-\end{verbatim}
-Items may be removed from the metadata by specifying a key.  If the
-key matches a metadata item, the item is removed from the metadata and
-returned; otherwise, \code{NULL} is returned.  If the key is not
-unique, then \emph{all} items corresponding to the key are removed,
-and the \tbd{first} item is returned.
-Care should be taken not to leak memory when appending an item for
-which the key already exists in the metadata (and is not
-\code{PS_META_NON_UNIQUE}).
-\begin{verbatim}
-/// delete entry from the metadata
-psMetadataItem *psMetadataRemove(psMetadata *restrict md, ///< metadata to delete from
-                                 const char *restrict key ///< Key to delete
-    );
-\end{verbatim}
-The metadata may be iterated over by (re-)setting the iterator for the
-appropriate \code{psMetadata}, and getting the next item.
-\code{psMetadataGetNext} has the ability to match the beginning of a
-key, e.g., if the user only wants to iterate through
-\code{IPP.machines.sky} and doesn't want to bother with
-\code{IPP.machines.detector}.  The iterator should iterate over every
-item of metadata --- even those that are non-unique.
-\begin{verbatim}
-/// reset the iterator to the start of the list
-void psMetadataSetIterator(psMetadata *md ///< metadata to set iterator for
-    );
-/// get the next item in the sequence
-psMetadataItem *psMetadataGetNext(psMetadata *restrict md, ///< metadata to get from
-                                  const char *restrict match, ///< Match this
-                                  int which ///< Which iterator to use
-    );
-\end{verbatim}
-Items may be found within the metadata by providing a key.  In the
-event that the key is non-unique, the first item is returned.
-\begin{verbatim}
-/// find the metadata with the specified key
-psMetadataItem *psMetadataLookup(const psMetadata *restrict md, ///< metadata to look up
-                                 const char *restrict key ///< Key to find
-    );
-\end{verbatim}
-Metadata items may be printed to an open file descriptor, optionally
-pre-pending a specified string.
-\begin{verbatim}
-/// print metadata item to the specified stream
-void psMetadataItemPrint(FILE *fd,              ///< file descriptor to write to
-                         const psMetadataItem *restrict md, ///< item of metadata to print
-                         const char *prefix        ///< print this at the beginning of each line
-    );
-\end{verbatim}
-\subsection{Detector and sky positions}
-Both detector and sky positions will be used extensively in the IPP.
-The first are linear coordinates which conform to Euclidean geometry
-while the second are angular coordinates for which additional care
-must often be taken.  We put these into two structures,
-\code{psPlaneCoord} and \code{psSphereCoord}, respectively.
-Partitioning these two will enable error-checking.
-\begin{verbatim}
-/** A point in 2-D space, with errors. */
-typedef struct {
-    double x;                           ///< x position
-    double y;                           ///< y position
-    double xErr;                        ///< Error in x position
-    double yErr;                        ///< Error in y position
-} psPlaneCoord;
-/** A point on the surface of a sphere, with errors */
-typedef struct {
-    double r;                           ///< RA
-    double d;                           ///< Dec
-    double rErr;                        ///< Error in RA
-    double dErr;                        ///< Error in Dec
-} psSphereCoord;
-\end{verbatim}
-Three major classes of coordinate transformations are necessary.
-First, linear coordinates from one frame must be converted to linear
-coordinates in a different frame of references.  Simple transformations
-of this type are independent of other quantities of the positions --
-they are simply mapping between two linear spaces.  In practice, these
-transformations may often be a function of the magnitude or color of
-the imaged object.  The second type of conversion is the
-transformation of linear coordinates to angular coordinates and
-vice-versa.  This conversion depends on the desired projection, and
-may represent the real mapping performed by the telescope or may
-simply represent a convenient mechanism to display 3D coordinates in
-useful forms.  The third conversion of interest is the transformation
-of one set of spherical coordinates to another set.  Frequently in
-astronomy, these conversions consist only of rotations between the two
-spherical coordinates systems, where the coordinates of the pole and
-equatorial rotation between the two systems define the
-transformation.  Conversions between standard coordinate systems such
-as Galactic, Ecliptic, and various epochs of the Celestial coordinates
-are represented by these spherical transformations.
-\subsubsection{Linear Coordinate Transformations}
-We specify two types of transforms between coordinate systems.  The
-first consists simply of two 2D polynomials to transform both
-components -- the output coordinates depend only on the input
-coordinates and no other quantities of objects at those coordinates.
-The second consists of two 4D polynomials in which the output
-coordinates are also specified to be a function of the magnitude and
-color of the object with the given coordinates.  This type of
-coordinate transformation is necessary to represent the
-(color-dependent) optical distortions caused by the atmosphere and
-camera optics, and the possibly effects of charge transfer
-inefficiency.  We specify two structures to represent the coefficients
-of these transformations:
-\begin{verbatim}
-/** A polynomial transformation between coordinate frames.  This may be a linear relationship, or may
- *  represent a higher-order transformation.
- */
-typedef struct {
-    psDPolynomial2D *x;
-    psDPolynomial2D *y;
-} psPlaneCoordXform;
-\end{verbatim}
-The \code{psDPolynomial2D} structures represent polynomials of
-arbitrary order as a function of two dimensions.  There is one of
-these structures for each of the two output dimensions.  As an
-example, consider the simple transformation from one linear coordinate
-frame \code{x,y} (say a single CCD) to a second frame \code{p,q} (say,
-a second CCD image). If we have only first order terms in the
-transformation \code{psPlaneCoordXform T}, the new coordinates would be
-represented by the terms:
+%
-\begin{verbatim}
-p = T.x->coeff[0][0] + x*T.x->coeff[1][0] + y*T.x->coeff[0][1];
-q = T.y->coeff[0][0] + x*T.y->coeff[1][0] + y*T.y->coeff[0][1];
-\end{verbatim}
+%
-where we have excluded the basic cross-term (\code{x*y}) by using the
-mask: \code{T.x->mask[1][1] = 0; T.y->mask[1][1] = 0;}
-\begin{verbatim}
-/** The optical distortion terms.  The lowest two terms are the x and y axis of the target system.  The higher
- *  two terms represent magnitude and color terms.
- */
-typedef struct {
-    psDPolynomial4D *x;
-    psDPolynomial4D *y;
-} psPlaneDistortion;
-\end{verbatim}
-Like \code{psPlaneCoordXform}, \code{psPlaneDistortion} contains two
-\code{psDPolynomial4D} structures representing polynomials of
-arbitrary order as a function of four, rather than two dimensions.
-There is one of these structures for each of the two output
-dimensions.  In this structure, the highest two dimensions could
-represent a magnitude and a color.  As an example, consider the simple
-transformation from one linear coordinate frame \code{x,y} (say a
-single CCD) of an object with magnitude and color \code{m,c} to a
-second frame \code{p,q} (say, a second CCD image). If we have only
-first order terms in the transformation \code{psPlaneCoordXform T}, the new
-coordinates would be represented by the terms:
+%
-\begin{verbatim}
-p = T.x->coeff[0][0][0][0] + x*T.x->coeff[1][0][0][0] + y*T.x->coeff[0][1][0][0] + m*T.x->coeff[0][0][1][0]
-    + c*T.x->coeff[0][0][0][1]
-q = T.y->coeff[0][0][0][0] + x*T.y->coeff[1][0][0][0] + y*T.y->coeff[0][1][0][0] + m*T.y->coeff[0][0][1][0]
-    + c*T.y->coeff[0][0][0][1]
-\end{verbatim}
+%
-where we have again excluded the cross-terms (\code{x*y}) by using the
-mask.
-We require corresponding functions to apply the transformations:
+%
-\begin{verbatim}
-/** apply the coordinate transformation to the given coordinate */
-psPlaneCoord *psPlaneCoordXformApply (psPlaneCoord *out, ///< Output coordinates, or NULL
-                                      const psPlaneCoordXform *frame, ///< coordinate transformation
-                                      const psPlaneCoord *coords ///< input coordiate
-    );
-/** apply the optical distortion to the given coordinate, magnitude, color */
-psPlaneCoord *psPlaneDistortionApply (psPlaneCoord *out, ///< Output coordinates, or NULL
-                                      const psPlaneDistortion *pattern, ///< optical distortion pattern
-                                      const psPlaneCoord *coords, ///< input coordinate
-                                      float mag, ///< magnitude of object
-                                      float color ///< color of object
-    );
-\end{verbatim}
+%
-\subsubsection{Celestial Coordinate Conversions}
-We need to be able to convert between ICRS, Galactic and Ecliptic
-coordinates, and potentially between arbitrary spherical coordinate
-systems.  All of these basic spherical transformations represent
-rotations of the spherical coordinate reference.  We specify a general
-transformation function which takes a structure,
-\code{psSphereCoordTransformation}, defining the transformation
-between two spherical coordinate systems (the structure contains the
-sines and cosines of the angles involved so as to minimize computation
-time for repeated transformations).  We also define a function to
-generate \code{psSphereCoordTransformation}, based on the three angles
-describing the location of the pole and the relative equatorial
-rotations of the two systems.  We also specify special functions to
-return the \code{psSphereCoordTransformation} for transformations
-between standard coordinate systems.
-\begin{verbatim}
-/** General spherical transformation */
-typedef struct {
-    double sin1, sin2, sin3, cos1, cos2, cos3; ///< Sines and cosines for transformation
-} psSphereCoordTransformation;
-\end{verbatim}
-The constructor and destructor are defined as follows:
-\begin{verbatim}
-/** Constructor */
-psSphereCoordTransformation *
-psSphereCoordTransformationAlloc(double pole1, ///< First location of pole
-                                 double pole2, ///< Second location of pole
-                                 double rotation ///< Rotation between systems
-    );
-/** Destructor */
-void psSphereCoordTranformationFree(psSphereCoordTransformation *trans ///< Transformation to destroy
-    );
-\end{verbatim}
-Spherical coordinates may be transformed by providing the transformation to
-\code{psSphereCoordTransform}:
-\begin{verbatim}
-/** Apply general spherical transformation */
-psSphereCoord *
-psSphereCoordTransform(const psSphereCoord *coord, ///< Coordinates to convert
-                       psSphereCoordSystem *sys ///< System to use to convert
-    );
-\end{verbatim}
-The following functions simply return the appropriate
-\code{psSphereCoordTransformation} to convert between predefined
-spherical coordinate systems (i.e., ICRS, Ecliptic and Galactic).
-\begin{verbatim}
-/** Return transformation structure to convert ICRS to Ecliptic */
-psSphereCoordTransformation *psSphereCoordTransformationItoE(void);
-/** Return transformation structure to convert Ecliptic to ICRS */
-psSphereCoordTransformation *psSphereCoordTransformationEtoI(void);
-/** Return transformation structure to convert ICRS to Galactic */
-psSphereCoordTransformation *psSphereCoordTransformationItoG(void);
-/** Return transformation structure to convert Galactic to ICRS */
-psSphereCoordTransformation *psSphereCoordTransformationGtoI(void);
-\end{verbatim}
-\subsubsection{Projections}
-We require functions to convert between spherical and linear
-coordinate systems based on a variety of projections.  The required
-projections include:
-\begin{itemize}
-\item TAN
-\item SIN
-\item AIT
-\item PAR
-\item GLS
-\end{itemize}
-The following functions will project and deproject (respectively)
-spherical coordinates:
-\begin{verbatim}
-/** Project spherical system onto a plane */
-psPlaneCoord *
-psCoordProject(const psSphereCoord *coord, ///< Spherical coordinates to project
-               const char *projection   ///< Projection to use
-    );
-/** Deproject plane onto spherical system */
-psSphereCoord *
-psCoordDeproject(const psPlaneCoord *coord, ///< Plane coordinates to deproject
-                 const char *projection ///< Projection to use
-    );
-\end{verbatim}
-\subsubsection{Offsets}
-We require a function to calculate the offset between two positions on
-the sky, as well as a function to apply an offset to a position.
-\begin{verbatim}
-/** Get offset (RA,Dec) on the sky between two positions position1 and position2 may not be identical */
-psSphereCoord *
-psSphereCoordGetOffset(const psSphereCoord *restrict position1, ///< Position 1
-                       const psSphereCoord *restrict position2, ///< Position 2
-                       const char *type         ///< Type of offset: Linear, Spherical/Arcsec,
-                                                ///< Spherical/Degreees etc
-    );
-/** Apply an offset to a position */
-psSphereCoord *
-psSphereCoordApplyOffset(const psSphereCoord *restrict position, ///< Position
-                         const psSphereCoord *restrict offset, ///< Offset
-                         const char *type               ///< Type of offset: Linear, Spherical/Arcsec,
-                                                        ///< Spherical/Degreees etc
-    );
-\end{verbatim}
-Note that these should propagate the errors appropriately.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Astronomical Images}
-\subsubsection{Overview}
-Above, we have defined a basic container for a single 2D collection of
-pixels (\code{psImage}), along with basic operations to manipulate the
-image pixels.  For astronomical applications, this data structure is
-insufficient for two reasons.  First, it does provide sufficient
-additional metadata to describe the data in detail.  Second, astronomy
-applications frequent involve multiple, related images.  For
-PanSTARRS, and for general astronomical applications, we require a
-richer collection of data structures which describe a very general
-image concept.  We have defined several layers in the hierarchy which
-are necessary to describe the image data which will be produced by the
-PanSTARRS Gigapixel cameras as well as other standard astronomical
-images.
-A simple 2D image is a basic data unit for much of astronomical
-imaging.  If we consider various optical and IR array cameras, a
-single readout of the detector produces a collection of pixels
-measurements.  We define our lowest-level astronomical image
-structure, \code{psReadout}, to contain the pixels produced by a
-single readout of the detector, along with metadata needed to define
-that readout: the origin and binning of the image relative to the
-original detector pixels explicitly in the structure, and pointers to
-the general metadata and derived objects, if any.
-A single detector may produce more than one read which is associated.
-For example, infrared detectors frequently produce an image
-immediately after the detector is reset followed by an image after the
-basic exposure is complete.  Both readouts correspond to the same
-pixels, though the binning or rastering may be different between the
-two readouts.  Another example is the video sequence produced by the
-PanSTARRS Gigapix camera guide cells, each of which represents a
-series of many images from a subraster of pixels in the detector
-readout portion.  The second level of our image container hierarchy,
-\code{psCell}, consists of a collection of readouts from a single
-detector.
-In the PanSTARRS Gigapix camera, the basic readout region is a
-fraction of the full imaging area of a single CCD chip.  The chip is
-divided into 64 cells, any fraction of which may have been readout
-for a given exposure.  In other cameras, such as Megacam at CFHT, the
-individual CCDs have multiple amplifiers addressing contiguous
-portions of the detector.  In such cameras, each amplifier produces a
-separate collection of pixels.  In the third level of our image
-container hierarchy, the data structure \code{psChip} represents a
-collection of different cells.
-The top level of our image container hierarchy is a complete focal
-plane array (\code{psFPA}).  This structure represents the collection
-of chips in the camera, all of which are read out in a given
-exposure.
-For example, take a mosaic camera consisting of eight $2k\times 4k$
-CCDs, each of which is read out through two amplifiers.  Then there
-would be sixteen cells in total, each of which is presumably $2k\times
-k$.  There would be eight chips, each consisting of two cells, and
-the focal plane consists of these eight chips.
-As another example, consider an observation by PS1.  The focal plane
-would consist of 60 chips, each of which consist of 64 cells (or less;
-a few cells may be dead).  Some cells (those containing guide stars
-for the orthogonal transfer) will contain multiple readouts.
-These data structures represent containers with which to carry around
-the collection of related image data.  There is no requirement on the
-functions or the structures that each instance of one of these data
-structures represent the physical hardware.  For example, it is not
-necessary that an instance of \code{psFPA} always carry the data for
-all 60 (or 64) Gigapixel camera OTAs.  The usage of these structures
-is such that all astronomical operations which apply to a CCD image
-should be performed on an instance of \code{psFPA}.  If a particular
-circumstance only requires a single 2D image, then that is represented
-by an instance of \code{psFPA} with one \code{psChip}, which in turn
-has one \code{psCell}, which in turn has one \code{psReadout}.
-These container levels also include in their definition the information
-needed to transform the coordinates in one of the levels to the
-coordinate system relevant at the higher levels.
-\tbd{the following discussions include astrometry issues which should
-  be deferred to the astrometry section}
-\subsubsection{A Readout}
-A readout is the result of a single read of a cell (or a portion
-thereof).  It contains a pointer to the pixel data, and additional
-pointers to the objects found in the readout, and the readout
-metadata.  It also contains the offset from the lower-left corner of
-the chip, in the case that the CCD was windowed.
-\begin{verbatim}
-/** a Readout: a collection of pixels */
-typedef struct {
-    const int x0, y0;                ///< Offset from the lower-left corner
-    const int nx, ny;                ///< Image binning
-    psImage *image;                  ///< imaging area of cell
-    psDlist *objects;                ///< objects derived from cell
-    psMetadata *md;                  ///< Readout-level metadata
-} psReadout;
-\end{verbatim}
-\subsubsection{A Cell}
-A cell consists of one or more readouts (usually only one except in the
-case that the cell has been used for fast guiding).  It also contains
-a pointer to the cell metadata, and a pointer to its parent chip.  On
-the astrometry side, it also contains coordinate transforms from the
-cell to the chip and, as a convenience, from the cell to the focal
-plane.  It is expected that these transforms will consist of two
-first-order 2D polynomials, simply specifying a translation, rotation
-and magnification; hence they are easily inverted, and there is no
-need to add reverse transformations.  We also add an additional
-transformation, which is intended to provide a ``quick and dirty''
-transform from the cell coordinates to the sky; this transformation
-not guaranteed to be as precise as the ``standard'' transformation of
-Cell $\rightarrow$ Chip $\rightarrow$ Focal Plane $\rightarrow$
-Tangent Plane $\rightarrow$ Sky, but will be faster.
-\begin{verbatim}
-/** a Cell: a collection of readouts.
- */
-typedef struct {
-    int nReadouts;                      ///< number of readouts in this cell realization; each may have its
-                                        ///< own image, objects and overscan.
-    struct psReadout *readouts;         ///< Readouts from the cell
-    psMetadata *md;                     ///< Cell-level metadata
-    psPlaneCoordXform *cellToChip;      ///< Transformations from cell coordinates to chip coordinates
-    psPlaneCoordXform *cellToFPA;       ///< Transformations from cell coordinates to FPA coordinates
-    psPlaneCoordXform *cellToSky;       ///< Quick and Dirty transformations from cell coordinates to sky
-    struct psChip  *parentChip;         ///< chip which contains this cell
-} psCell;
-\end{verbatim}
-\subsubsection{A Chip}
-A chip consists of one or more cells (according to the number of
-amplifiers on the CCD).  It contains a pointer to the chip metadata,
-and a pointer to the parent focal plane.  For astrometry, it contains
-a coordinate transform from the chip to the focal plane.  It is
-expected that this transforms will consist of two second-order 2D
-polynomials; hence we think that it is prudent to include a reverse
-transformation which will be derived from numerically inverting the
-forward transformation.
-\begin{verbatim}
-/** a Chip: a collection of cells.  Not all valid cells in a chip need to be listed in an
- *  instance of psChip.
- */
-typedef struct {
-    int nCells;                         ///< Number of Cells assigned
-    struct psCell *cells;               ///< Cells in the Chip
-    psMetadata *md;                     ///< Chip-level metadata
-    psPlaneCoordXform *chipToFPA;       ///< Transformations from chip coordinates to FPA coordinates
-    psPlaneCoordXform *FPAtoChip;       ///< Transformations from FPA coordinates to chip
-    struct psFPA *parentFPA;            ///< FPA which contains this chip
-} psChip;
-\end{verbatim}
-\subsubsection{A Focal Plane}
-A focal plane consists of one or more chips (according to the number
-of pieces of contiguous silicon).  It contains pointers to the focal
-plane metadata and the exposure information.  For astrometry, it
-contains a transformation from the focal plane to the tangent plane
-and the fixed pattern residuals.  It is expected that the
-transformation will consist of two 4D polynomials (i.e.\ a function of
-two coordinates in position, the magnitude of the object, and the
-color of the object) in order to correct for optical distortions and
-the effects of the atmosphere; hence we think that it is prudent to
-include a reverse transformation which will be derived from
-numerically inverting the forward transformation.  Since colors are
-involved in the transformation, it is necessary to specify the color
-the transformation is defined for.  We also include some values to
-characterize the quality of the transformation: the root mean square
-deviation for the x and y transformation fits, and the $\chi^2$ for
-the transformation fit.
-\begin{verbatim}
-/** a Focal plane array: a collection of chips.  Not all chips in a camera need to be listed in an instance of
- *  psFPA.
- */
-typedef struct {
-    int nChips;                         ///< Number of Cells assigned
-    int nAlloc;                         ///< Number of Cells available
-    struct psChip *chips;               ///< Chips in the Focal Plane Array
-    psMetadata *md;                     ///< FPA-level metadata
-    psPlaneDistortion *TPtoFP;          ///< Transformation term from
-    psPlaneDistortion *FPtoTP;          ///< Transformation term from
-    psFixedPattern *pattern;            ///< Fixed pattern residual offsets
-    const psExposure *exp;              ///< information about this exposure
-    psPhotSystem colorPlus, colorMinus; ///< Colour reference
-    float rmsX, rmsY;                   ///< Dispersion in astrometric solution
-    float chi2;                         ///< chi^2 of astrometric solution
-} psFPA;
-\end{verbatim}
-\subsubsection{Exposure information}
-We need several quantities from the telescope in order to make a
-first guess at the astrometric solution.  From these quantities,
-further quantities can be derived and stored for later use.
-\begin{verbatim}
-/** Exposure information from the telescope */
-typedef struct {
-    // Telescope longitude, latitude and height are stored separately, since they don't change with pointing
-    const double ra, dec;               ///< Telescope boresight
-    const double ha;                    ///< Hour angle
-    const double zd;                    ///< Zenith distance
-    const double az;                    ///< Azimuth
-    const double lst;                   ///< Local Sidereal Time
-    const float mjd;                    ///< MJD of observation
-    const float rotAngle;               ///< Rotator position angle
-    const float temp;                   ///< Air temperature, for estimating refraction
-    const float pressure;               ///< Air pressure, for calculating refraction
-    const float humidity;               ///< Relative humidity, for calculating refraction
-    const float exptime;                ///< Exposure time
-    /* Derived quantities */
-    const float posAngle;               ///< Position angle
-    const float parallactic;            ///< Parallactic angle
-    const float airmass;                ///< Airmass, calculated from zenith distance
-    const float pf;                     ///< Parallactic factor
-    const char *cameraName;             ///< name of camera which provided exposure
-    const char *telescopeName;          ///< name of telescope which provided exposure
-} psExposure;
-\end{verbatim}
-\subsubsection{Constructors and Destructors}
-Each of the above structures needs an appropriate constructor and
-destructor.  Other than \code{psExposure}, which contains significant
-non-pointer types, the constructors should not take any arguments, and
-the destructors should only take the structure to be destroyed.
-The constructor for \code{psExposure} is specified below.
-\begin{verbatim}
-/** Constructor */
-psExposure *
-psExposureAlloc(double ra, double dec,  ///< Telescope boresight
-                double ha,              ///< Hour angle
-                double zd,              ///< Zenith distance
-                double az,              ///< Azimuth
-                double lst,             ///< Local Sidereal Time
-                float mjd,              ///< MJD
-                float rotAngle,         ///< Rotator position angle
-                float temp,             ///< Temperature
-                float pressure,         ///< Pressure
-                float humidity,         ///< Relative humidity
-                float exptime           ///< Exposure time
-                );
-\end{verbatim}
-\subsection{Astrometry}
-Astrometry is a basic functionality required for the IPP that will be
-used repeatedly, both for low-precision (roughly where is my favorite
-object?) and high-precision (what is the proper motion of this star?).
-As such, it must be flexible, yet robust.  Accordingly, we will wrap
-the StarLink Astronomy Libraries (SLALib), which has already been
-developed.
-\subsubsection{Coordinate frames}
-\label{sec:coordinateFrames}
-There are five coordinate frames that we need to worry about for the
-purposes of astrometry:
-\begin{itemize}
-\item Cell: $(x,y)$ in pixels --- raw coordinates;
-\item Chip: $(X,Y)$ in pixels --- the location on the silicon;
-\item Focal Plane: $(p,q)$ in microns --- the location on the focal plane;
-\item Tangent Plane: $(l,m)$ in arcsec from the telescope boresight; and
-\item Sky: (RA,Dec) --- ICRS.
-\end{itemize}
-The following steps are required to convert from the cell coordinates to
-the sky \tbd{convert Cell-$>$Chip and Chip-$>$FP to psDistortion}:
-\begin{itemize}
-\item Cell $\longleftrightarrow$ Chip: two 2D polynomials, $(X,Y) = f(x,y)$;
-\item Chip $\longleftrightarrow$ FP: two 2D polynomials, $(p,q) = g(X,Y)$;
-\item FP $\longleftrightarrow$ TP: two 4D polynomials, $(l,m) =
-h(p,q,m,c)$, where $m$ and $c$ are the magnitude and color of the
-object, respectively; and
-\item TP $\longleftrightarrow$ Sky: SLALib transformation using a
-transform pre-computed for each pointing.
-\end{itemize}
-Note that the transformation between the Focal Plane and the Tangent
-Plane is a four-dimensional polynomial, in order to account for any
-possible dependencies in the astrometry on the stellar magnitude and
-color; the former serves as a check for charge transfer
-inefficiencies, while the latter will correct chromatic refraction,
-both through the atmosphere and the corrector lenses.
-We require structures to contain each of the above transformations as
-well as the pixel data.
-\subsubsection{SLALib information}
-SLALib requires several elements to perform the transformations
-between the tangent plane and the sky.  Pre-computing these quantities
-for each exposure means that subsequent transformations are faster.
-For historical reasons, this structure is known colloquially as
-``Wallace's Grommit''.
-\begin{verbatim}
-/** Information needed (by SLALIB) to convert Apparent to Observed Position */
-typedef struct {
-    const double latitude;              ///< geodetic latitude (radians)
-    const double sinLat, cosLat;        ///< sine and cosine of geodetic latitude
-    const double abberationMag;         ///< magnitude of diurnal aberration vector
-    const double height;                ///< height (HM)
-    const double temperature;           ///< ambient temperature (TDK)
-    const double pressure;              ///< pressure (PMB)
-    const double humidity;              ///< relative humidity (RH)
-    const double wavelength;            ///< wavelength (WL)
-    const double lapseRate;             ///< lapse rate (TLR)
-    const double refractA, refractB;    ///< refraction constants A and B (radians)
-    const double longitudeOffset;       ///< longitude + eqn of equinoxes + ``sidereal UT'' (radians)
-    const double siderealTime;          ///< local apparent sidereal time (radians)
-} psGrommit;
-\end{verbatim}
-The \code{psGrommit} is calculated from telescope information for the
-particular exposure:
-\begin{verbatim}
-/** Constructor */
-psGrommit *
-psGrommitAlloc(const psExposure *exp    ///< Relevant exposure
-    );
-/** Destructor */
-void
-psGrommitFree(psGrommit *grommit        ///< Grommit to destroy
-    );
-\end{verbatim}
-\subsubsection{Fixed Pattern}
-The fixed pattern is a correction to the general astrometric solution
-formed by summing the residuals from many observations.  The intent is
-to correct for higher-order distortions in the camera system on a
-coarse grid (larger than individual pixels, but smaller than a single
-cell).  Hence, in addition to the offsets, we need to specify the size
-and scale of the grid in $x$ and $y$, as well as the origin of the
-grid.
-\begin{verbatim}
-/** The fixed pattern residual offsets.  These are specified via a coarse grid of x and y offsets. */
-typedef struct {
-    int nX, nY;                         ///< Number of elements in x and y
-    double x0, y0;                      ///< Position of the lower-left corner of the grid on the focal plane
-    double xScale, yScale;              ///< Scale of the grid
-    double **x, **y;                    ///< The grid of offsets in x and y
-} psFixedPattern;
-\end{verbatim}
-\subsubsection{Position Finding}
-We require functions to return the structure containing given
-coordinates.  For example, we want the chip that corresponds to the
-focal plane coordinates $(p,q) = (-1.234,+5.678)$.  These routines
-handle the one-to-many problem --- i.e., for one given focal plane
-coordinate, there are many chips that this coordinate may be
-correspond to; these functions will select the correct one.
-\begin{verbatim}
-/** returns Chip in FPA which contains the given FPA coordinate */
-psChip *
-psChipInFPA (psChip *out,               ///< Chip to return, or NULL
-             const psPlaneCoord *coord  ///< coordinate in FPA
-             const psFPA *fpa,          ///< FPA description
-             );
-/** returns Cell in Chip which contains the given chip coordinate */
-psCell *
-psCellInChip(psCell *out,               ///< Cell to return, or NULL
-             const psPlaneCoord *coord  ///< coordinate in chip
-             const psChip *chip,        ///< chip description
-             );
-/** Return the cell in FPA which contains the given FPA coordinates */
-psCell *
-psCellInFPA(psCell *out,                ///< Cell to return, or NULL
-            const psPlaneCoord *coord   ///< Coordinate in FPA
-            const psFPA *fpa,           ///< FPA description
-            );
-\end{verbatim}
-\subsubsection{Conversion Functions}
-We require functions to convert between the various coordinate frames
-(Section~\ref{sec:coordinateFrames}).  The hierarchy of the coordinate
-frames and the transformations between each are shown in
-Figure~\ref{fig:coco}.  The functions that employ the transformations
-are shown in Figure~\ref{fig:cocoFunc}.  In addition to
-transformations between each adjoining coordinate frame in the
-hierarchy, we also require higher-level functions to convert between
-the Cell and Sky coordinate frames; these will simply perform the
-intermediate steps.
-\begin{figure}
-\psfig{file=coordinateFrames,height=7in,angle=-90}
-\caption{The coordinate systems in the \PS{} IPP, and the relation
-between each by transformations contained in the appropriate
-structures.}
-\label{fig:coco}
-\end{figure}
-\begin{figure}
-\psfig{file=coordinateConv,height=7in,angle=-90}
-\caption{Conversion between coordinate systems by PSLib.}
-\label{fig:cocoFunc}
-\end{figure}
-The function prototypes are:
-\begin{verbatim}
-/** Convert (RA,Dec) to cell and cell coordinates */
-psPlaneCoord *
-psCoordSkyToCell(psPlaneCoord *out,     ///< Coordinates to return, or NULL
-                 psCell *cell,          ///< Cell to return
-                 const psSphereCoord *in, ///< Input coordinates
-                 const psFPA *fpa       ///< FPA description
-                 );
-/** Convert cell and cell coordinate to (RA,Dec) */
-psSphereCoord *
-psCoordCellToSky(psSphereCoord *out,    ///< Coordinates to return, or NULL
-                 const psPlaneCoord *coord ///< cell coordinates to transform
-                 const psCell *cell,    ///< Cell to get coordinates for
-                 );
-/** Quick and dirty cell to (RA,Dec) --- employs cellToSky transformation */
-psSphereCoord *
-psCoordCellToSkyQuick(psSphereCoord *out, ///< Coordinates to return, or NULL
-                      const psPlaneCoord *coord ///< cell coordinates to transform
-                      const psCell *cell, ///< Cell description
-                      );
-/** Convert (RA,Dec) to tangent plane coords */
-psPlaneCoord *
-psCoordSkyToTP(psPlaneCoord *out,       ///< Coordinates to return, or NULL
-               const psSphereCoord *coord ///< input Sky coordinate
-               const psGrommit *grommit, ///< Grommit for fast conversion
-               );
-/** Convert tangent plane coords to focal plane coordinates */
-psPlaneCoord *
-psCoordTPtoFPA(psPlaneCoord *out,       ///< Coordinates to return, or NULL
-               const psPlaneCoord *coord ///< input TP coordinate
-               const psFPA *fpa,        ///< FPA description
-               );
-/** converts the specified FPA coord to the coord on the given Chip */
-psPlaneCoord *
-psCoordFPAtoChip (psPlaneCoord *out,    ///< Coordinates to return, or NULL
-                  const psPlaneCoord *coord ///< input FPA coordinate
-                  const psChip *chip,   ///< Chip of interest
-                  );
-/** converts the specified Chip coord to the coord on the given Cell */
-psPlaneCoord *
-psCoordChiptoCell (psPlaneCoord *out,   ///< Coordinates to return, or NULL
-                   const psPlaneCoord *coord ///< input Chip coordinate
-                   const psCell *cell,  ///< Cell of interest
-                   );
-/** converts the specified Cell coord to the coord on the parent Chip */
-psPlaneCoord *
-psCoordCelltoChip (psPlaneCoord *out,   ///< Coordinates to return, or NULL
-                   const psPlaneCoord *coord ///< input Cell coordinate
-                   const psCell *cell,  ///< Cell description
-                   );
-/** converts the specified Chip coord to the coord on the parent FPA */
-psPlaneCoord *
-psCoordChiptoFPA (psPlaneCoord *out,            ///< Coordinates to return, or NULL
-                  const psPlaneCoord *coord     ///< input Chip coordinate
-                  const psChip *chip,   ///< Chip description
-                  );
-/** Convert focal plane coords to tangent plane coordinates */
-psPlaneCoord *
-psCoordFPAToTP(psPlaneCoord *out,               ///< Coordinates to return, or NULL
-               const psPlaneCoord *coord ///< input FPA coordinate
-               const psFPA *fpa,        ///< FPA description
-               );
-/** Convert tangent plane coords to (RA,Dec) */
-psSphereCoord *
-psCoordTPtoSky(psSphereCoord *out,      ///< Coordinates to return, or NULL
-               const psPlaneCoord *coord ///< input TP coordinate
-               const psGrommit *grommit, ///< Grommit for fast conversion
-               );
-/** Convert Cell coords to FPA coordinates */
-psPlaneCoord *
-psCoordCellToFPA(psPlaneCoord *out,     ///< Coordinates to return, or NULL
-                 const psPlaneCoord *coord ///< Input cell coordinates
-                 const psCell *cell,    ///< Cell description
-                 );
-\end{verbatim}
-\subsubsection{Additional functions}
-We require additional functions to perform general functions which
-will be useful for astrometry.  Given coordinates on the sky, we
-need to get the airmass, the parallactic angle, and an estimate of
-the atmospheric refraction.
-\begin{verbatim}
-/** Get the airmass for a given position and sidereal time */
-float
-psGetAirmass(const psSphereCoord *coord, ///< Position on the sky
-             double siderealTime,       ///< Sidereal time
-             float height               ///< Height above sea level
-             );
-\end{verbatim}
-\begin{verbatim}
-/** Get the parallactic angle for a given position and sidereal time */
-float
-psGetParallactic(const psSphereCoord *coord, ///< Position on the sky
-                 double siderealTime    ///< Sidereal time
-                 );
-\end{verbatim}
-\begin{verbatim}
-/** Estimate atmospheric refraction, along the parallactic */
-float
-psGetRefraction(float colour,           ///< Colour of object
-                psPhotSystem colorPlus, ///< Colour reference
-                psPhotSystem colorMinus, ///< Colour reference
-                const psExposure *exp   ///< Telescope pointing information, for airmass, temp and pressure
-                );
-\end{verbatim}
-\begin{verbatim}
-/** Calculate the parallax factor */
-double
-psGetParallaxFactor(const psExposure *exp ///< Exposure details
-    );
-\end{verbatim}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Photometry}
-Photometric observations are performed in an instrumental photometric
-system, and must be related to other photometric systems.  We
-require a data structure which defines a photometric system, as well
-as a structure to define the transformation between photometric
-systems.
-The photometric system is defined by the psPhotSystem structure.
-A photometric system is identified by a human-readable \code{name}
-(ie, SDSS.g, Landolt92.B, GPC1.OTA32.r).  Each photometric system is
-given a unique identifier \code{ID}.  Observations taken with a
-specific camera, detector, and filter represent their own photometric
-system, and it may be necessary to perform transformations between
-these systems.  Photometric systems associated with observations from
-a specific camera/detector/filter combination can be associated with
-those components.
-\begin{verbatim}
-typedef struct {
-    const int ID;                       ///< ID number for this photometric system
-    const char *name;                   ///< Name of photometric system
-    const char *camera;                 ///< Camera for photometric system
-    const char *filter;                 ///< Filter used for photometric system
-    const char *detector;               ///< Detector used for photometric system
-} psPhotSystem;
-\end{verbatim}
-The following structure defines the transformation between two
-photometric systems.
-\begin{verbatim}
-typedef struct {
-    psPhotSystem src;                   ///< Source photometric system
-    psPhotSystem dst;                   ///< Destination photometric system
-    psPhotSystem pP, pM;                ///< Primary colour reference
-    psPhotSystem sP, sM;                ///< Secondary colour reference
-    float pA, sA;                       ///< Colour offset for primary and secondary references
-    psPolynomial3D transform;           ///< Transformation from source to destination
-} psPhotTransform;
-\end{verbatim}
-The transformation between two photometric systems may depend on the
-airmass of the observation and on the colors of the object of
-interest.  For a specific observation, such a transformations can be
-defined as a polynomial function of the color of the star and the
-airmass of the observations.  If sufficient data exists, the
-transformation between the photometric systems may include more than
-one color, constraining the curvature of the stellar spectral energy
-distributions.  This latter term may be significant for stars which
-are highly reddened, for example.  Derived photometric quantities may
-have been corrected for airmass variations, in which case only color
-terms may be measurable.  The structure defines the transformation
-between a source photometric system (\code{src}) and a target
-photometric system (\code{dst}).  The photometric system of a primary
-color is defined by \code{pP, pM} such that the color is constructed
-as $pP - pM$.  A secondary color is defined by \code{sP, sM}.  For
-both, a reference color is specified (\code{pA, sA}): the polynomial
-transformation terms refer to colors in the form $pP - pM - pA$.  The
-transformation is specified as a 3D polynomial.  For a star of
-magnitude $M_{\rm src}$ in the source photometric system, with
-additional magnitude information in the other systems $M_{\rm pP}$,
-$M_{\rm pM}$, $M_{\rm sP}$, $M_{\rm sM}$, observed at an airmass of
-$z$, the magnitude of the star in the target system $M_{\rm dst}$ is
-given by: $M_{\rm dst} = M_{\rm src} + transform(z, M_{\rm pP} -
-M_{\rm pM} - pA, M_{\rm sP} - M_{\rm sM} - sA)$.
-\TBD{Really want a set of polynomials defined for specific colour
-ranges.}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Astronomical objects}
-\subsubsection{Positions of Major SS Objects}
-We require the ability to calculate the position of major Solar System
-objects, as well as Lunar phase.
-\begin{verbatim}
-/** Get Sun Position */
-psSphereCoord *
-psGetSunPos(float mjd                   ///< MJD to get position for
-    );
-/** Get Moon position */
-psSphereCoord *
-psGetMoonPos(float mjd,                 ///< MJD to get position for
-             double latitude,           ///< Latitude for apparent position
-             double longitude           ///< Longitude for apparent position
-    );
-/** Get Moon phase */
-float
-psGetMoonPhase(float mjd                ///< MJD to get phase for
-    );
-/** Get Planet positions */
-psSphereCoord *
-psGetSolarSystemPos(const char *solarSystemObject, ///< Named S.S. object
-                    float mjd           ///< MJD to get position for
-    );
-\end{verbatim}
-\tbd{The rest of this section is tentative}
-We may need a variety of other SLALib-type functions:
-\begin{itemize}
-\item Lunation calculations to/from psTime
-\item rise/set for specified object
-\item 12 \& 18 deg twilight for closest sunrise / sunset
-\item length of current night
-\end{itemize}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%\appendix
+%
-%\pagebreak
-%\section{API Summary: all functions}
+%
-%\subsection{System Utilities}
-%\input{psSystemGroup.tex}
+%
-%\subsection{Data Containers}
-%\input{psDataGroup.tex}
+%
-%\subsection{Math Utilities}
-%\input{psMathGroup.tex}
+%
-%\subsection{Astronomy Functions}
-%\input{psAstroGroup.tex}
+%
-%\pagebreak
-%\section{API Summary: all structures}
-%\input{psStructures.tex}
-\bibliographystyle{plain} \bibliography{panstarrs}
-\end{document}

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 747 for trunk/doc/pslib/psLibSDRS.tex

Legend:

trunk/doc/pslib/psLibSDRS.tex

Download in other formats: