Index: /trunk/doc/pslib/psLibSDRS.tex
===================================================================
--- /trunk/doc/pslib/psLibSDRS.tex	(revision 255)
+++ /trunk/doc/pslib/psLibSDRS.tex	(revision 255)
@@ -0,0 +1,917 @@
+\documentclass[panstarrs]{panstarrs}
+%\documentclass[panstarrs]{panstarrs}
+
+% basic document variables
+\title{Pan-STARRS IPP Library Design}
+\author{}
+\shorttitle{PSLib Design}
+\group{Pan-STARRS Algorithm Group}
+\project{Pan-STARRS Image Processing Pipeline}
+\organization{Institute for Astronomy}
+\version{DR}
+\docnumber{PSDC-xxx-xxx}
+% note the use of the docnumber & version number:
+% the complete PSDC document number is given by
+% \thedocnumber-\theversion
+
+\begin{document}
+\maketitle
+
+% -- Revision History --
+% provide explicit values for the old versions
+% use '\theversion' for the current version (set above)
+% use \hline between each table row
+\RevisionsStart
+% version     Date         Description
+\theversion & 2003 Mar 11 & Hacking \\
+\RevisionsEnd
+
+\pagebreak
+\tableofcontents
+
+\pagebreak 
+\pagenumbering{arabic}
+
+%Here's a reference to another document:
+%\href{file:utils#psDlist}{utils.pdf}; the \code{#psDlist} points to
+%\CODE|\hlabel{psDlist}| in the file \file{utils.tex} (this is like
+%a regular \code{\label}, but defines the hyperlink too).
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{\PS{} Library (PSLib)}
+
+PSLib is a library of basic functions required for IPP (and MOPS)
+operations.  We expect that it will be, to a large extent, the major
+workhorse of the IPP.
+
+PSLib will consist of a variety of data structures and associated APIs
+to perform common functions.  It will provide low-level capabilities
+in each of the following areas:
+\begin{itemize}
+\item System utilities;
+\item Data containers;
+\item Data manipulation; and
+\item Astronomy-specific tasks.
+\end{itemize}
+
+Below we sketch out approximately 120 APIs which specify the required
+capabilities.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{System Utilities}
+
+\note{INSERT HERE STUFF FROM RHL'S \file{utils.pdf}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Data Containers}
+
+We require general data containers, so that associated values (e.g.\
+the elements of an array) may be connected as a whole.  We require the
+following types of containers:
+\begin{itemize}
+\item Arrays;
+\item Doubly-linked lists; and
+\item Hashes.
+\end{itemize}
+
+\note{INSERT HERE STUFF FROM RHL'S \file{utils.pdf}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Data manipulation}
+
+We require general data manipulation functions, which will act upon
+data (in particular, arrays/vectors).  We require the following capabilities:
+\begin{itemize}
+\item Bit masks;
+\item Vector and image arithmetic;
+\item Sorting;
+\item Statistics;
+\item Matrix operations and linear algebra;
+\item (Fast) Fourier Transforms;
+\item General functions; and
+\item Minimisation and fitting routines.
+\end{itemize}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Bit masks}
+
+Bit masks are required in order to turn options on and off.  We require the
+capability to have a bit mask of arbitrary length (i.e., not limited by the
+length of a \code{long}, say).
+
+\begin{verbatim}
+/** A bitmask of arbitrary length. */
+typedef struct {
+    int n;                              //!< Number of chars that form the mask
+    char *bits;                         //!< The bits
+} psBitMask;
+\end{verbatim}
+
+We also require the corresponding constructor and destructor:
+
+\begin{verbatim}
+/** Constructor */
+psBitMask *psBitMaskAlloc(int n         //!< Number of bits required
+    );
+ 
+/** Destructor */
+void psBitMaskFree(psBitMask *restrict myMask //!< Bit mask to destroy
+    );
+\end{verbatim}
+
+Four basic operations on bit masks are required:
+\begin{itemize}
+\item Set a bit;
+\item Check if a bit is set;
+\item \code{OR} two bit masks; and
+\item \code{AND} two bit masks.
+\end{itemize}
+The corresponding APIs are defined below.
+
+\begin{verbatim}
+/** Set a bit mask */
+psBitMask *
+psBitMaskSet(psBitMask *outMask,	//!< Output bit mask or NULL
+	     const psBitMask *myMask,	//!< Input bit mask
+	     int bit			//!< Bit to set
+    );
+\end{verbatim}
+
+\begin{verbatim}
+/** Check a bit mask.  Returns true or false */
+int
+psBitMaskTest(const psBitMask *checkMask, //!< Bit mask to check
+	      int bit			//!< Bit to check
+    );
+\end{verbatim}
+
+\begin{verbatim}
+/** apply the given operator to two bit masks */
+psBitMask *
+psBitMaskOp(psBitMask *outMask,		//!< Output bit mask or NULL
+	    const psBitMask *restrict inMask1, //!< Input bit mask 1
+	    char *operator,		//!< bit mask operator (AND, OR, XOR)
+	    const psBitMask *restrict inMask2 //!< Input bit mask 2
+    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Vector and Image Arithmetic}
+
+We will need to be able to perform various operations on vectors and
+matrices, e.g.\ dividing one image by another, subtracting a vector
+from an image, etc.  Both binary operations and unary operations are
+required.
+
+\begin{verbatim}
+/** Perform a binary operation on two data items (psImage, psVector, psScalar).
+*/
+psType *
+psBinaryOp (void *out,                  ///< destination (may be NULL)
+            void *in1,                  ///< first input
+            char *operator,             ///< operator
+            void *in2                   ///< second input
+);
+\end{verbatim}
+
+\begin{verbatim}
+/** Perform a binary operation on two data items (psImage, psVector, psScalar).
+*/
+psType *
+psUnaryOp (void *out,                   ///< destination (may be NULL)
+           void *in,                    ///< input
+           char *operator,              ///< operator
+);
+\end{verbatim}
+
+Note that these functions should return the appropriate type (i.e.,
+the \code{psType} return type refers to \code{psVector} and
+\code{psImage} and \code{psScalar}).  It is expected that the
+implementation of these functions will employ pre-processor macros to
+perform the onerous task of creating the loops.
+
+It is desirable to use the same functions for both vectors and
+matrices, so inputs are \code{void*}; this necessitates that vectors
+and matrices each have a type element at a pre-determined and constant
+location in the \code{struct}.  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}
+typedef struct {
+    psType type;
+    union {
+	int i;
+	float f;
+	double d;
+	complex float c;
+    }
+} 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.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Sorting}
+
+We require the ability to sort an array of floating-point values:
+
+\begin{verbatim}
+/** Sort an array. Inputs not restrict-ed to allow sort in place */
+psFloatArray *
+psSort(psFloatArray *out,		//!< Sorted array to return. May be NULL
+       const psFloatArray *myArray	//!< Array to sort
+    );
+\end{verbatim}
+
+We also require the ability to sort one array based on another.  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:
+
+\begin{verbatim}
+/** Sort an array, along with some other stuff.  Returns an index array */
+psIntArray *
+psSortIndex(psIntArray *restrict out;	//!< Output index array (may be NULL)
+	    const psFloatArray *restrict myArray //!< Array to sort
+    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Statistics}
+
+We require a very general statistics function, which, given an array
+of floating-point values, will be able to calculate the following
+statistics:
+\begin{itemize}
+\item Sample mean;
+\item Sample median;
+\item Sample mode;
+\item Sample standard deviation;
+\item Sample upper and lower quartiles;
+\item Robust mean, associated error and number of values used to calculate;
+\item Robust median, associated error and number of values used to calculate;
+\item Robust mode, associated error and number of values used to calculate;
+\item Robust standard deviation;
+\item Robust upper and lower quartiles;
+\item Clipped mean, associated error and number of values used to calculate;
+\item Clipped standard deviation; and
+\item Minimum and maximum value in array.
+\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
+function:
+
+\begin{verbatim}
+/** Do Statistics on an array.  Returns a status value. */
+psStats *
+psArrayStats(const psFloatArray *restrict myArray, //!< Array to be analysed
+	     const psIntArray *restrict maskArray, //!< Ignore elements where (maskArray & maskVal) != 0
+						   //!< May be NULL
+	     unsigned int maskVal,	//!< Only mask elements with one of these bits set in maskArray
+	     psStats *stats		//!< stats structure defines stats to be calculated and how
+	     );
+\end{verbatim}
+
+
+We also require to be able to generate histograms, given a list of
+upper and lower bounds for each of the bins, and maximum and minimum
+values.
+
+\begin{verbatim}
+/** Histograms */
+typedef struct {
+    const psFloatArray *restrict lower;	//!< Lower bounds for the bins
+    const psFloatArray *restrict upper; //!< Upper bounds for the bins
+    psIntArray *nums;			//!< Number in each of the bins
+    const float minVal, maxVal;		//!< Minimum and maximum values
+    int minNum, maxNum;			//!< Number below the minimum and above the maximum
+} psHistogram;
+\end{verbatim}
+
+The constructors and destructor follow.  We specify two constructors,
+so that the bounds of the bins may either be specified explicitly, or
+implicitly through simply specifying an upper and lower limit along
+with the size of the bins.
+
+\begin{verbatim}
+/** Constructor */
+psHistogram *
+psHistogramAlloc(float lower,		//!< Lower limit for the bins
+	         float upper,		//!< Upper limit for the bins
+	         float size		//!< Size of the bins
+    );
+\end{verbatim}
+
+\begin{verbatim}
+/** Generic constructor */
+psHistogram *
+psHistogramAllocGeneric(const psFloatArray *restrict lower, //!< Lower bounds for the bins
+		        const psFloatArray *restrict upper, //!< Upper bounds for the bins
+		        float minVal,	//!< Minimum value
+		        float maxVal	//!< Maximum value
+    );
+\end{verbatim}
+
+\begin{verbatim}
+/** Destructor */
+void
+psHistogramFree(psHistogram *restrict myHist //!< Histogram to destroy
+    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Matrix operations and linear algebra}
+
+We require the ability to perform basic linear algebra on matrices, in
+order to solve equations:
+
+\begin{verbatim}
+/** A matrix */
+typedef struct {
+    int xSize, ySize;			//!< Dimensions in x and y
+    float *restrict *restrict value;	//!< Values in matrix
+} psMatrix;
+\end{verbatim}
+
+Constructor and destructor are:
+
+\begin{verbatim}
+/** Constructor */
+psMatrix *
+psMatrixAlloc(int Xdimen,			//!< x dimension of new matrix
+	      int Ydimen			//!< y dimension of new matrix
+	      );
+\end{verbatim}
+
+\begin{verbatim}
+/** Destructor */
+void
+psMatrixFree(psMatrix *restrict myMatrix	//!< Matrix to destroy
+	     );
+\end{verbatim}
+
+We require the following basic matrix operations:
+\begin{itemize}
+\item Invert a matrix;
+\item Calculate a matrix determinant;
+\item Perform matrix addition, subtraction and multiplication;
+\item Transpose a matrix; and
+\item Convert a matrix to a vector.
+\end{itemize}
+The corresponding APIs are:
+
+\begin{verbatim}
+/** Invert matrix.  Not using restrict, to allow inversion to be done in-place */
+psMatrix *
+psMatrixInvert(psMatrix *out,		//!< Matrix to return, or NULL
+	       const psMatrix *myMatrix, //!< Matrix to be inverted
+	       float *restrict determinant //!< Determinant to return, or NULL
+    );
+\end{verbatim}
+
+\begin{verbatim}
+/** Matrix determinant */
+float
+psMatrixDeterminant(const psMatrix *restrict myMatrix //!< Matrix to get determinant for
+		    );
+\end{verbatim}
+
+\begin{verbatim}
+/** Matrix operations */
+psMatrix *
+psMatrixOp(psMatrix *out,		//!< Matrix to return, or NULL
+	   const psMatrix *matrix1,	//!< Matrix 1
+	   const char *op,		//!< Operation to perform
+	   const psMatrix *matrix2	//!< Matrix 2
+	   );
+\end{verbatim}
+
+\begin{verbatim}
+/** Transpose Matrix */
+psMatrix *
+psMatrixTranspose(psMatrix *out,	//!< Matrix to return, or NULL
+		  const psMatrix *myMatrix //!< Matrix to transpose
+		  );
+\end{verbatim}
+
+\begin{verbatim}
+/** Convert matrix to vector.  Intended for a 1-d matrix. */
+psVector *
+psMatrixToVector(psMatrix *myMatrix	//!< Matrix to convert
+    );
+\end{verbatim}
+
+Note that \code{psMatrixOp} differs from \code{psBinaryOp} in that
+\code{psBinaryOp} acts on each element independently, while
+\code{psMatrixOp} performs classical matrix arithmetic, involving row
+and column operations.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{(Fast) Fourier Transforms}
+
+We require the ability to calculate the (fast) fourier transforms of
+both floating-point and complex floating-point type arrays, as well as
+the power spectrum.  The APIs are:
+
+\begin{verbatim}
+/** Return Fourier Transform of an array */
+psComplexArray *
+psRealFFT(psComplexArray *restrict out,	//!< Output array to be returned; may be NULL
+	  const psFloatArray *restrict myArray //!< Input array
+	  );
+\end{verbatim}
+
+\begin{verbatim}
+/** Return [inverse?] Fourier Transform of a complex array.  May NOT be done in-place (restrict-ed) */
+psComplexArray *
+psComplexFFT(psComplexArray *restrict out, //!< Output array to be returned; may be NULL
+	     const psComplexArray *restrict myArray, //!< Input array
+	     int sign			//!< +1 or -1 to indicate direction of FT
+	     );
+\end{verbatim}
+
+\begin{verbatim}
+/** Return Power spectrum of a array */
+psFloatArray *
+psPowerSpec(psFloatArray *restrict out,	//!< Output array to be returned
+	    const psFloatArray *restrict myArray //!< Input array
+	    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{General functions}
+
+We require two types of general functions which will be used in fitting:
+Gaussians and Polynomials.  The Gaussian is easily defined:
+
+\begin{verbatim}
+/** evaluate 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: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] */
+float
+psGaussian(float x,			//!< Value at which to evaluate
+	   float mean,			//!< Mean for the Gaussian
+	   float stddev			//!< Standard deviation for the Gaussian
+	   );
+\end{verbatim}
+
+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$,
+colour 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 //!< Number of terms in x and y
+    );
+void psDPolynomial2DFree(psDPolynomial2D *restrict myPoly //!< Polynomial to destroy
+    );
+\end{verbatim}
+
+And the polynomials may be evaluated:
+
+\begin{verbatim}
+/** Evaluate 3D polynomial (double precision) */
+double
+psEvalDPolynomial3D(double x,		//!< Value x at which to evaluate
+		    double y,		//!< Value y at which to evaluate
+		    double z,		//!< Value z at which to evaluate
+		    const psDPolynomial3D *restrict myPoly //!< Coefficients for the polynomial
+		    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Minimisation and fitting routines}
+
+We require a general minimisation routine, a routine that will
+specifically minimise $\chi^2$ given a list of data with associated
+errors, and a function that will analytically determine the polynomial
+that goes through all the specified points.  The APIs are:
+
+\begin{verbatim}
+/** Minimize a particular function */
+psFloatArray *
+psMinimize(float (*myFunction)(const psFloatArray *restrict),	//!< Function to minimize
+	   psFloatArray *restrict initialGuess //!< Initial guess
+	   );
+\end{verbatim}
+
+\begin{verbatim}
+/** Minimize chi^2 for input data */
+psFloatArray *
+psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict,
+				  const psFloatArray *restrict), //!< Model to fit; receives domain and
+								 //!< parameters
+	       const psFloatArray *restrict domain, //!< The domain values for the corresponding measurements
+	       const psFloatArray *restrict data, //!< Data to fit
+	       const psFloatArray *restrict errors, //!< Errors in the data
+	       psFloatArray *restrict initialGuess, //!< Initial guess
+	       const psIntArray *restrict guessMask //!< 1 = fit for parameter, 0 = hold parameter constant
+	       );
+\end{verbatim}
+
+\begin{verbatim}
+/** Derive a polynomial that goes through the points --- can be done analytically */
+psFloatArray *
+psGetArrayPolynomial(const psFloatArray *restrict ord, //!< Ordinates (or NULL to just use the indices)
+		     const psFloatArray *restrict coord //!< Coordinates
+    );
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Astronomy-Specific Functions}
+
+Some basic, relatively simple astronomy-specific functions are
+required which will serve as the foundation for building the Phase $N$
+modules.  These functions are not expected to cover every forseeable
+function, but will serve as the building blocks of more complicated
+processing functions.
+
+We require functions covering each of the following areas:
+\begin{itemize}
+\item Astrometry;
+\item Dates and times;
+\item Image handling;
+\item Metadata;
+\item Detector and sky positions;
+\item Astronomical objects; and
+\item Photometry.
+\end{itemize}
+These are each dealt with below.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Astrometry}
+
+Astrometry is a basic functionality required for the IPP that will be
+used repeatedly, both for low-precision (roughly where is my favourite
+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{Terminology}
+
+Some brief review of terminology would be useful so that previous
+definitions do not influence the understanding of this document.
+
+A ``frame'' is a read of the detector.
+
+A ``cell'' is defined as the smallest element of the detector readout;
+usually associated with an amplifier.  Correspondingly, each cell has
+its own overscan region.  There may be multiple frames in a cell if
+the cell was used to provide fast guiding.
+
+A ``chip'' is defined as a contiguous piece of silicon, and consists
+of a group of cells.
+
+A ``focal plane'' is defined as a mosaic of chips, and consists of a
+group of chips.
+
+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
+2k$.  There would be eight chips, each consisting of two cells, and
+the focal plane consists of these eight chips.
+
+\subsubsection{Coordinate frames}
+
+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 colour of the
+object, respectively; and
+\item TP $\longleftrightarrow$ Sky: SLALib transformation using a
+transform pre-computed for each pointing.
+
+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
+colour; the latter serves as a check for charge transfer
+inefficiencies, while the latter will correct chromatic refraction,
+both through the atmosphere and the corrector lenses.
+
+\textbf{[If the magnitude terms serve to check CTI, then shouldn't we
+put them in the cell <--> chip section?]}
+
+We require structures to contain each of the above transformations as
+well as the pixel data.
+
+\subsubsection{A Frame}
+
+A frame is the result of a single read of a cell (or a portion
+thereof).  It contains a pointer to the data pixels with a
+separate pointer to the overscan pixels, and additional pointers to
+the objects found in the frame, and the frame 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 Frame: a collection of pixels */
+typedef struct {
+    int x0, y0;				//!< Offset from the lower-left corner
+    psImage *image;                     ///< imaging area of cell 
+    psDlist *objects;			///< objects derived from cell
+    psImage *overscan;                  ///< bias region (subimage) of cell
+    psMetaDataSet *md;			//!< Frame-level metadata
+} psFrame;
+\end{verbatim}    
+
+
+\subsubsection{A Cell}
+
+A cell consists of one or more frames (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.
+
+\begin{verbatim}
+/** a Cell: a collection of frames.
+ */
+typedef struct {
+    int nFrames;			///< number of frames in this cell realization; each may have its own
+					///< image, objects and overscan.
+    struct psFrame *frames;		//!< Frames from the cell
+    psMetaDataSet *md;			///< Cell-level metadata
+
+    psCoordXform *cellToChip;		///< Transformations from cell coordinates to chip coordinates
+    psCoordXform *cellToFPA;		///< Transformations from cell coordinates to FPA coordinates
+
+    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 expect 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
+
+    psMetaDataSet *md;			///< Chip-level metadata
+    psCoordXform *chipToFPA;            ///< Transformations from chip coordinates to FPA coordinates
+
+    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
+colour of the object) in order to correct for optical distortions and
+the effects of the atmosphere; hence we expect that it is prudent to
+include a reverse transformation which will be derived from
+numerically inverting the forward transformation.  Since colours are
+involved in the transformation, it is necessary to specify the colour
+the transformation is defined for.  We also include some values to
+characterise 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
+
+    psMetaDataSet *md;			///< FPA-level metadata 
+    psDistortion *TPtoFP;		///< Transformation term from 
+    psDistortion *FPtoTP;		///< Transformation term from 
+    psFixedPattern *pattern;		//!< Fixed pattern residual offsets
+    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{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 {
+    double latitude;                    ///< geodetic latitude (radians)
+    double sinLat, cosLat;              ///< sine and cosine of geodetic latitude
+    double abberationMag;               ///< magnitude of diurnal aberration vector
+    double height;                      ///< height (HM)
+    double temperature;                 ///< ambient temperature (TDK)
+    double pressure;                    ///< pressure (PMB)
+    double humidity;                    ///< relative humidity (RH)
+    double wavelength;                  ///< wavelength (WL)
+    double lapseRate;                   ///< lapse rate (TLR)
+    double refractA, refractB;          ///< refraction constants A and B (radians)
+    double longitudeOffset;             ///< longitude + eqn of equinoxes + ``sidereal UT'' (radians)
+    double siderealTime;                ///< local apparent sidereal time (radians)
+} psGrommit;
+\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
+    double ra, dec;                     //!< Telescope boresight
+    double ha;                          //!< Hour angle
+    double zd;                          //!< Zenith distance
+    double az;                          //!< Azimuth
+    double lst;                         //!< Local Sidereal Time
+    float mjd;                          //!< MJD of observation
+    float rotAngle;                     //!< Rotator position angle
+    float temp;                         //!< Air temperature, for estimating refraction
+    float pressure;                     //!< Air pressure, for calculating refraction
+    float humidity;                     //!< Relative humidity, for calculating refraction
+    float exptime;                      //!< Exposure time
+    /* Derived quantities */
+    float posAngle;                     //!< Position angle
+    float parallactic;                  //!< Parallactic angle
+    float airmass;                      //!< Airmass, calculated from zenith distance
+    float pf;                           //!< Parallactic factor
+    char *cameraName;                   ///< name of camera which provided exposure
+    char *telescopeName;                ///< name of telescope which provided exposure
+    psGrommit *grommit;                 //!< Data needed to convert from the sky to the tangent plane
+} psExposure;
+\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
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Dates and times}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Image handling}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Metadata}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Detector and sky positions}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Astronomical objects}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Photometry}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\bibliographystyle{plain} \bibliography{panstarrs}
+
+\end{document}
+