Context Navigation

← Previous Changeset
Next Changeset →

Changeset 3564

Timestamp:

Mar 30, 2005, 11:14:48 AM (21 years ago)

Author:

eugene

Message:

serious reorg for release 13

Location:

trunk/doc/pslib

Files:

: 3 added
: 5 edited

ChangeLogADD.tex (modified) (5 diffs)
ChangeLogSDRS.tex (modified) (4 diffs)
psLibADD.tex (modified) (63 diffs)
psLibSDRS.tex (modified) (24 diffs)
psLibSDRS_Astrom.tex (modified) (4 diffs)
rotations.ps (added)
rotations.sxd (added)
transforms.ps (added)

Legend:

: Unmodified
: Added
: Removed

trunk/doc/pslib/ChangeLogADD.tex

-              r3438
+              r3564
 \subsection{Changes from version 06 to version 07}
+\subsection{Changes from version 06 (7 September 2004) to version 07 (24 November 2004)}
 \begin{itemize}
 …
 \end{itemize}
 \subsection{Changes from version 07 to version 08}
+\subsection{Changes from version 07 (24 November 2004) to version 08 (8 January 2005)}
 \begin{itemize}
 …
 \end{itemize}
 \subsection{Changes from version 08 to version 09}
+\subsection{Changes from version 08 (8 January 2005) to version 09 (14 February 2005)}
 \begin{itemize}
 …
   \code{PS_RESAMPLE_LAGRANGE}.
 \item Added section on FITS WCS.
+\end{itemize}
+\subsection{Changes from version 09 (14 February 2005) to version 10 (29 March 2005}
+\begin{itemize}
 \item Changes to the Time section:
 \begin{itemize}
 …
 \item Verbatim Fortran code for UT1 interpolation
 \end{itemize}
+\item section reorganization:
+\begin{itemize}
+\item renamed Astronomical Image Manipulations to Image Manipulations
+\item moved Image Manipulations to own subsection before Astronomy Utilities
+\item promoted all PSLib subsections to sections
 \end{itemize}
+\end{itemize}

trunk/doc/pslib/ChangeLogSDRS.tex

-              r3537
+              r3564
 %%% $Id: ChangeLogSDRS.tex,v 1.77 2005-03-29 03:42:16 price Exp $
+%%% $Id: ChangeLogSDRS.tex,v 1.78 2005-03-30 21:14:48 eugene Exp $
 \subsection{Changes from version 00 to version 01}
 …
 \end{itemize}
 \subsection{Changes from Revision 12 (9 February 2005) to Revision 13 (Present)}
+\subsection{Changes from Revision 12 (9 February 2005) to Revision 13 (30 March 2005)}
 \begin{itemize}
 …
   from a file.
 \item \code{psMetadataAddV} changed to use \code{va_list} parameter (bug 312).
+\item Added \code{psPixels} structure, \code{psPixelsToMask},
+  \code{psMaskToPixels}, \code{psPixelsConcatenate} and
+  \code{psPixelsTransform}.
 \item Modified \code{psImageTransform} in preparation for image combination.
+\item Added \code{psPixels} structure and related functions
 \item Added \code{psPlaneTransformDeriv}.
 \item Added \code{psImageGrowMask}.
+\item Added Earth Orientation Calculations Section
 \item Changes to the Time section:
 \begin{itemize}
 …
 \item Adding logical operations (and, or) to \code{psBinaryOp}.
 \end{itemize}
+\item Substantial reorganization:
+\begin{itemize}
+\item Moved Metadata, Database, and XML sections to new section
+\item Re-named Detector \& Sky Coordinates to Linear \& Spherical Coordinates
+\item Moved Exposure and Observatory information out of 'Astronomical Images'
+\item Moved Celestial Coordinate Systems out of 'Detector \& Sky Coordinates'
+\item Added Atmospheric Effects section, incorporating psGrommit and airmass functions from other sections)
+\item Moved Fixed Pattern out of Astronomical Images
+\end{itemize}
+\end{itemize}

trunk/doc/pslib/psLibADD.tex

-              r3448
+              r3564
 %%% $Id: psLibADD.tex,v 1.69 2005-03-18 20:40:14 jhoblitt Exp $
+%%% $Id: psLibADD.tex,v 1.70 2005-03-30 21:14:48 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 …
 \project{Pan-STARRS Image Processing Pipeline}
 \organization{Institute for Astronomy}
 \version{09}
+\version{10}
 \docnumber{PSDC-430-006}
+\setcounter{tocdepth}{5} % lowest level to be included in toc
 \newcommand\citealt{}
 …
 & 2004 May 21 & Added section on 2D Chebyshev fitting, then removed. \\ \hline
 & 2004 Jun 22 & modified stats specification \\ \hline
+--05 & ??? & ??? \\ \hline
+& 2004 Sep 7 & Frozen for PSLib-2 \\ \hline
+& 2004 Jul 13 & \\ \hline
+& 2004 Aug 16 & \\ \hline
+& 2004 Sep 01 & \\ \hline
+& 2004 Sep 07 & Frozen for PSLib-2 \\ \hline
 & 2004 Nov 24 & Frozen for Cycle 4 \\ \hline
 & 2005 Jan 21 & Draft for Cycle 5 \\ \hline
 & 2005 Feb 14 & Frozen for Cycle 5 \\ \hline
+& 2005 Mar 21 & Draft for Cycle 6 \\ \hline
 \RevisionsEnd
 …
 \pagenumbering{arabic}
 \section{Pan-STARRS Library PSLib}
 \subsection{Math Utilities}
 \subsubsection{Sorting}
+% \section{Pan-STARRS Library PSLib}
+\section{PSLib Math Utilities}
+\subsection{Sorting}
 A variety of sorting algorithms exist, with a wide range in speed for
 …
 \code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}.
 \subsubsection{Smoothing: Boxcar and Gaussian}
+\subsection{Smoothing: Boxcar and Gaussian}
 \label{smooth}
 …
 \end{equation}
 \subsubsection{Statistics}
+\subsection{Statistics}
 The general statistics function \code{psStats} performs a variety of
 …
 sample and robust estimators.
 \paragraph{Sample Statistics}
+\subsubsection{Sample Statistics}
 We define the following statistical terms, assuming there is a set of
 data elements $x_i$ with (standard) errors $\sigma_i$.
 \subparagraph{Mean}
+\paragraph{Mean}
 The simple mean is defined as:
 …
 \end{equation}
 \subparagraph{Weighted Mean}
+\paragraph{Weighted Mean}
 The weighted mean is defined as:
 …
 standard definition of the mean.
 \subparagraph{Median}
+\paragraph{Median}
 The median is defined as the value for which 50\% of the data values
 …
 calculating the sample median.
 \subparagraph{Upper and Lower Quartiles}
+\paragraph{Upper and Lower Quartiles}
 The upper and lower quartiles ($U_{\frac{1}{4}}$ and
 …
 the sample quartiles.
 \subparagraph{Standard Deviation}
+\paragraph{Standard Deviation}
 The standard deviation of the sample is given by:
 …
 \paragraph{Clipped Statistics}
+\subsubsection{Clipped Statistics}
 The clipped statistics are used to determine the mean and standard
 …
 \bar{x}| > k \sigma_i$.
 \paragraph{Robust Statistics}
+\subsubsection{Robust Statistics}
 The robust version of the statistics provides estimators of basic
 …
 quartiles are estimated in the same manner as above.
 \paragraph{Histograms}
+\subsubsection{Histograms}
 When calculating histograms in the presence of known errors in the
 …
 Note that the total adds to one --- the number of values added.
 \subsubsection{Matrix Operations}
+\subsection{Matrix Operations}
 In this section, we define the linear algebra operations performed on
 …
 \code{gsl_linalg_LU_decomp}.
 \paragraph{LU Decomposition}
+\subsubsection{LU Decomposition}
 \label{LUdecomp}
 …
 \end{equation}
 \paragraph{Calculate a matrix determinant}
+\subsubsection{Calculate a matrix determinant}
 The determinant $D$ of a matrix $a_{ij}$ is calculated from the
 …
 shall be used.
 \paragraph{Solving a Linear Equation}
+\subsubsection{Solving a Linear Equation}
 The LU decomposition of a matrix may be used to solve the
 …
 \end{eqnarray}
 \paragraph{Invert a matrix}
+\subsubsection{Invert a matrix}
 Inversion of a matrix using the LU decomposition is performed by
 …
 operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}.
 \paragraph{Perform matrix addition, subtraction and multiplication}
+\subsubsection{Perform matrix addition, subtraction and multiplication}
 Matrix binary arithmetic operations differ from image binary
 …
 \times$.
 \paragraph{Transpose a matrix}
+\subsubsection{Transpose a matrix}
 The transpose of a matrix is simply the reorganization of the matrix
 …
 where $M_{ij}$ is the matrix to be transposed.
 \paragraph{Convert a matrix to a vector}
+\subsubsection{Convert a matrix to a vector}
 Matrix-to-vector conversion is only defined for a matrix that has a
 …
 matrix is converted to a \code{PS_DIMEN_TRANV}-type vector.
 \subsubsection{Fitting}
 \paragraph{Chi-squared}
+\subsection{Fitting}
+\subsubsection{Chi-squared}
 \label{chisq}
 …
 \end{equation}
 \paragraph{General Polynomial Fitting}
+\subsubsection{General Polynomial Fitting}
 Given a set of data values $y_i$ with errors $\sigma_i$, related to
 …
 (section~\ref{LUdecomp}).
 \subsubsection{Non-linear Minimization}
 \paragraph{Levenberg-Marquardt Method}
+\subsection{Non-linear Minimization}
+\subsubsection{Levenberg-Marquardt Method}
 In the Levenberg-Marquardt Method (LMM; see NR \S 15.5), we make a
 …
 \paragraph{Powell's method}
+\subsubsection{Powell's method}
 Powell's method is a type of ``Direction Set'' methods in
 …
 manner until the advances along the vectors are smaller than some
 pre-defined tolerance.  Such direction set methods, including Powell's
 Quadratically Convergent method are discussed in NR\S10.5.
+Quadratically Convergent method are discussed in NR \S 10.5.
 We will use for our algorithm the modified version of Powell's
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsubsection{Polynomials}
+\subsection{Polynomials}
 \label{sec:polynomials}
 …
 $-1 < x < 1$.
 \paragraph{Multi-dimensional polynomials}
+\subsubsection{Multi-dimensional polynomials}
 Multi-dimensional polynomials shall be composed of multiplications of
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsubsection{(Fast) Fourier Transforms}
+\subsection{(Fast) Fourier Transforms}
 (Fast) Fourier Transforms (FFTs) shall be implemented using the
 …
 library}.
 \paragraph{FFTW Plans}
+\subsubsection{FFTW Plans}
 FFTW requires the user to create a ``plan'' for each transform size,
 …
 initialisation of the PSLib FFT functions and saved at the conclusion.
 \paragraph{Function mapping}
+\subsubsection{Function mapping}
 The forward and reverse transforms call the corresponding
 …
 place to avoid the need to pad the input array to hold the output.
 \paragraph{More Complicated Functions}
+\subsubsection{More Complicated Functions}
 \code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve
 …
 caller, and choose to normalise by $1/N^2$.
+\subsection{Image Manipulations}
+\subsubsection{Interpolation}
+Interpolation is needed in various image manipulation operations,
+including rotation and resampling.  We have specified a function to
+perform the interpolation using one of several possible interpolation
+methods, defined below.  It is important in the discussions that
+follow to remember that a pixel with column,row if $i,j$ has
+coordinate at the center of $i+0.5,j+0.5$ and corners with coordinates
+from $i,j$ to $i+1,j+1$.  Thus, the interpolation of a coordinate
+$x,y$ = 5.0,4.0 is a value midway between the four pixels with
+column,row of (5,4), (5,5), (6,4), (6,5).
+\paragraph{Nearest Pixel Interpolation ({\tt PS\_INTERPOLATE\_FLAT})}
+In this interpolation, the value of the closest pixel is returned.
+This is equivalent to pixel duplication or replication.
+\paragraph{Bilinear Interpolation ({\tt PS\_INTERPOLATE\_BILINEAR})}
+In this interpolation, the value at the coordinate is calculated using
+linear interpolation in two dimensions from the four nearest neighbor
+pixels.  The bilinear interpolation value at a coordinate $x,y$
+depends on the four nearest neighbor pixels and the fractional
+distance $fx,fy$ of the given coordinates from the centers of those
+four pixels.  Consider four neighboring pixels at column,row of $i,j$,
+$i+1,j$, $i,j+1$, and $i+1,j+1$ with pixel values $V_{0,0}$,
+$V_{1,0}$, $V_{0,1}$, $V_{1,1}$.  The value at $x,y$ is given by:
+\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \]
+This expression is more efficiently evaluated by factoring and
+calculating the expresion as:
+\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \]
+\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \]
+Note that the values of $f_x$ and $f_y$ require some care.  Given a
+coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 -
+int(f_x - 0.5)$.  For example, when interpolating the value at
+(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5)
+and the fractional coordinate values $f_x, f_y = 0.3, 0.7$.  The
+resulting coordinate would be contained within the pixel at column,row
+(5,5).
+\paragraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})}
+Because it would be slow to specify the size of the kernel
+dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8
+pixels in each dimension (a kernel of size 2 pixels in each dimension
+is handled by the bilinear interpolation).  These correspond to the
+options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and
+\code{PS_INTERPOLATE_LANCZOS4}, respectively.
+Given a position on the input image, $(x_0,y_0)$, a kernel is derived
+according to pixels local to the position:
+\begin{equation}
+  h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N)
+\end{equation}
+where
+\begin{eqnarray}
+  \delta x & = & x - x_0 \\
+  \delta y & = & y - y_0 \\
+  {\rm sinc}(z) & = & \sin(z)/z
+\end{eqnarray}
+and $N$ corresponds to the choice of kernel size.  For $N = 2$, the
+kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le
+$).  For $N = 3$, the kernel size is 6 pixels in each dimension
+(i.e., $-3 < \delta x \le 3$).  For $N = 4$, the kernel size is 8
+pixels in each dimension (i.e., $-4 < \delta x \le 4$).
+The interpolated value at the given position, $(x_0,y_0)$, is then
+simply the dot product of the kernel and the fluxes:
+\begin{equation}
+  f(x_0,y_0) = \sum_R f(x,y) h(x,y)
+\end{equation}
+where $R$ is the region defined by the kernel size, and $f(x,y)$ is
+the flux at the pixel position.
+For further information, see the
+\href{http://terapix.iap.fr/IMG/pdf/swarp.pdf}{SWarp manual}.
+\subsubsection{Image Cuts and Slices}
+Several functions specify operations which manipulate a collection of
+pixels to return a statistic on the pixel collection.  In the simplest
+case, these are trivial to define: if the boundaries of the region of
+interest are specified along integral pixel coordinates, then the
+pixels used to measure the statistic are always an exact integer.
+This is the case for the function \code{psImageSlice} which requires a
+starting coordinate which is an integer and a width in both dimensions
+which is an integer.  For the case of the functions \code{psImageCut}
+and \code{psImageRadialCut}, the situation is a bit more subtle.  In
+both of these cases, the region is unlikely to contain only whole
+pixels and some choices must be made.
+One posibility which we reject is to identify the fractional pixels
+which are overlapped by the region of interest and add that fraction
+of the pixel's flux when calculating the statistic of interest.  This
+is computationally intensive, and not necessarily well defined for all
+statistics.
+In PSLib, we instead identify the pixels overlapped by the region, use
+the complete set of pixel values, treating all pixels equally, and
+renormalize as needed.  To perform this, the region of interest is
+laid on top of the image pixels.  Any pixels which overlap the region
+are identified as part of the input sample.  The statistic (ie, sample
+mean, robust mode, etc), is then calculated on this collection of
+pixels.  If the output statistic is an average value, the measured
+value is reported.  If the output statistic is a sum value (sum of
+counts, sum of pixels), then the value is renormalized by the ratio of
+pixels used in the calculation to the pixel area of the region of
+interest.  For example, if the sum within a radial aperture is
+requested, the circle of the specified radius and center is placed on
+the pixel grid.  Any pixels which touch the circle are then placed in
+a list to be analysed.  The statistic of interest is the measured for
+this collection of pixels.  In the case of a circular aperture which
+is centered at the coordinate (2,2) and has a radius of 2, the number
+of pixels which are touched by the circle is 16, while the total pixel
+area of the circle is 12.57 square pixels.  In this case, the pixel
+sum is renormalized by the ratio (12.57/16.00).
+\paragraph{Radial Cuts}
+Consider an image with pixels $x_i,y_i$ and a reference coordinate
+$x_c, y_c$.  We want to construct a radial cut by measuring statistics
+for pixels in a sequence of radial annulii $r_s < r < r_e$.  For each
+annulus, we need to select the pixels which fall within this annulus.
+The coordinates of the center of pixel $i,j$ are $i+0.5,j+0.5$.  A
+given pixel has a distance from the reference coordinate of $dX = x_c
+- i - 0.5, dY = y_c - j - 0.5$.  The pixels to be used for a given
+radial annulus are all of those pixels for which $r_s < \sqrt{dX^2 +
+  dY^2} < r_e$.  This is more efficiently calculated by comparing the
+square of the radii and distances.  All pixels which satisfy the above
+condition are included in a specific annular radius.  All average
+quantities are calculated directly from the pixel ensemble
+statistics.
+\paragraph{Arbitrary Linear Cuts}
+Select the pixels which lie along a line following steps of 1 pixel
+length:
+\begin{verbatim}
+  dX = xe - xs;
+  dY = ye - ys;
+  L = hypot (dX, dY);
+  dX = dX / L;
+  dY = dY / L;
+  REALLOCATE (xvec[0].elements, float, MAX (L, 1));
+  REALLOCATE (yvec[0].elements, float, MAX (L, 1));
+  xvec[0].Nelements = L;
+  yvec[0].Nelements = L;
+  V = (float *)buf[0].matrix.buffer;
+  for (i = 0; i < L; i++) {
+    xi = xs + i*dX - 0.5;
+    yi = ys + i*dY - 0.5;
+    xvec[0].elements[i] = i;
+    yvec[0].elements[i] = V[xi + Nx*yi];
+  }
+\end{verbatim}
+\subsubsection{Image Rotation}
+Image rotation can be performed in two possible ways under different
+circumstances, identified in the following discussion.
+In the simplest case, the rotation angle is an integer multiple of 90
+degrees ($\pi/2$ rad).  In these cases, the input and output pixels
+have a one-to-one mapping.  If the input image has dimensions of $N_x,
+N_y$, then the output image will have dimensions of either $N_x, N_y$
+(for even multiples of 90 degrees) or $N_y, N_x$ (for odd multiples).
+If the angle of the rotation is not a multiple of 90, then the output
+pixels necessarily result from the interpolation of several input
+pixels.  In this case, for an input image of dimensions $N_x, N_y$ and
+rotation angle $\theta$, the output image has dimensions $Lx = |N_x
+\cos \theta| + |N_y \sin \theta|$ and $Ly = |N_x \sin \theta| + |N_y
+\cos \theta|$, each dimension rounded up to the nearest integer as
+needed.  Every pixel in the output image is in general derived from an
+interpolation over 4 neighboring pixels.  The coordinate of a pixel in
+the output image ($i,j$) corresponds to a fractional pixel coordinate
+($x,y$) in the input image according to:
+\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \]
+\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \]
+where the offset coordinate ($i_o,j_o$) depends on the sign of the
+sine of the angle $\theta$.  If the sign of that sine is positive, the
+offset coordinate is ($N_y\sin\theta$,0), otherwise it is
+(0,$-N_x\sin\theta$).
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Astronomy Utilities}
+\pagebreak
+\section{PSLib Astronomy Utilities}
 Most of the astronomy utilities will be implemented through wrapping
 …
 the next release}
 \subsubsection{Time}
+\subsection{Time}
 Correct time representation is \emph{critical} in astronomical software.  PSLib
 …
 ``1970-01-01T00:00:00Z'' UTC.
 \paragraph{Coordinated Universal Time (UTC)}
+\subsubsection{Coordinated Universal Time (UTC)}
 Coordinated Univeral Time (UTC) is defined by the International
 …
 timezone is forbidden.}
 \paragraph{International Atomic Time (TAI)}
+\subsubsection{International Atomic Time (TAI)}
 International Atomic Time or Temps Atomique International (TAI) is a system of
 …
 seconds since the UNIX epoch of ``1970-01-01T00:00:00Z''.
 \paragraph{Leap-seconds}
+\subsubsection{Leap-seconds}
 Leap seconds keep UTC within 0.9s of UT1.  The offset between TAI and
 …
 This data is available from: \code{ftp://maia.usno.navy.mil/ser7/tai-utc.dat}
 \paragraph{Gregorian dates to seconds}
+\subsubsection{Gregorian dates to seconds}
 The Perl code below, based on an algorithm described in the book ``Calendrical
 …
 \end{verbatim}
 Outputs year, month, day as \code{$y, $m, $d}.
+%$
 \emph{The above code was taken [and slightly altered] from
 …
 \paragraph{Universal Time (UT1)}
+\subsubsection{Universal Time (UT1)}
 \label{sec:ut1}
 …
 IERS publications references above, and should be interpolated in the same way.
 \paragraph{Julian Date and Modified Julian Date}
+\subsubsection{Julian Date and Modified Julian Date}
 The follow definitions of Julian Date (JD) and Modified Julian Date (MJD) was
 …
 http://www.iers.org/iers/earth/resolutions/UAI\_b1.html}.
 \subparagraph{Julian Date}
+\paragraph{Julian Date}
 \begin{verbatim}
 …
 \end{verbatim}
 \subparagraph{Modified Julian Date}
+\paragraph{Modified Julian Date}
 \begin{verbatim}
 …
 \end{verbatim}
 \subparagraph{JD and MJD conversion}
+\paragraph{JD and MJD conversion}
 Conversion between \code{psTime} values and MJD and JD are determined
 …
 \end{equation}
 \paragraph{Terrestrial Time (TT)}
+\subsubsection{Terrestrial Time (TT)}
 Terrestrial Time (TT) is defined as a fixed offset from TAI.
 …
 \end{equation}
 \paragraph{TT as Julian Centuries since J2000.0}
+\subsubsection{TT as Julian Centuries since J2000.0}
 The algorithm for calulating GMST requires TT formated in Julian centruies
 …
 \end{equation}
 \paragraph{UT1 as Julian Centuries since J2000.0}
+\subsubsection{UT1 as Julian Centuries since J2000.0}
 The algorithm for calulating GMST requires UT1 be formated in Julian centuries
 …
 \end{equation}
+\paragraph{Greenwich Mean Sidereal Time (GMST)}
+\subsubsection{Local Mean Sidereal Time (LMST)}
+Local Mean Sidereal Time (LMST) is Greenwich Mean Sideral Time (GMST) plus the
+observer's location in East longitude. Calculating LMST requires the input of
+Universal Time (UT1), Terrestrial Dynamical Time (TT) and a longitude (measured
+East of Greenwich).
+\begin{equation}
+LMST = GMST00(t_u, t) + longitude
+\end{equation}
+Gives $LMST$ in seconds.
+\subsubsection{Greenwich Mean Sidereal Time (GMST)}
 Greenwich Mean Sidereal Time (GMST) is caclulated from UT1 and TT.  This
 …
 Gives $GMST00$ in seconds.
+\paragraph{Longitude}
+\subsubsection{Longitude}
 Longitudes are often expressed in the form of decimal degrees while the
 …
 \end{equation}
+\paragraph{Local Mean Sidereal Time (LMST)}
+Local Mean Sidereal Time (LMST) is Greenwich Mean Sideral Time (GMST) plus the
+observer's location in East longitude. Calculating LMST requires the input of
+Universal Time (UT1), Terrestrial Dynamical Time (TT) and a longitude (measured
+East of Greenwich).
+\begin{equation}
+LMST = GMST00(t_u, t) + longitude
+\end{equation}
+Gives $LMST$ in seconds.
+\paragraph{Polar Motion}
+\subsubsection{Polar Motion}
+\tbd{move this to Earth Motion section}
 The polar coordinates, $x_p$ and $y_p$, required for the transformation from
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{Astronomical Image Manipulations}
+\paragraph{Interpolation}
+Interpolation is needed in various image manipulation operations,
+including rotation and resampling.  We have specified a function to
+perform the interpolation using one of several possible interpolation
+methods, defined below.  It is important in the discussions that
+follow to remember that a pixel with column,row if $i,j$ has
+coordinate at the center of $i+0.5,j+0.5$ and corners with coordinates
+from $i,j$ to $i+1,j+1$.  Thus, the interpolation of a coordinate
+$x,y$ = 5.0,4.0 is a value midway between the four pixels with
+column,row of (5,4), (5,5), (6,4), (6,5).
+\subparagraph{Nearest Pixel Interpolation ({\tt PS\_INTERPOLATE\_FLAT})}
+In this interpolation, the value of the closest pixel is returned.
+This is equivalent to pixel duplication or replication.
+\subparagraph{Bilinear Interpolation ({\tt PS\_INTERPOLATE\_BILINEAR})}
+In this interpolation, the value at the coordinate is calculated using
+linear interpolation in two dimensions from the four nearest neighbor
+pixels.  The bilinear interpolation value at a coordinate $x,y$
+depends on the four nearest neighbor pixels and the fractional
+distance $fx,fy$ of the given coordinates from the centers of those
+four pixels.  Consider four neighboring pixels at column,row of $i,j$,
+$i+1,j$, $i,j+1$, and $i+1,j+1$ with pixel values $V_{0,0}$,
+$V_{1,0}$, $V_{0,1}$, $V_{1,1}$.  The value at $x,y$ is given by:
+\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \]
+This expression is more efficiently evaluated by factoring and
+calculating the expresion as:
+\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \]
+\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \]
+Note that the values of $f_x$ and $f_y$ require some care.  Given a
+coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 -
+int(f_x - 0.5)$.  For example, when interpolating the value at
+(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5)
+and the fractional coordinate values $f_x, f_y = 0.3, 0.7$.  The
+resulting coordinate would be contained within the pixel at column,row
+(5,5).
+\subparagraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})}
+Because it would be slow to specify the size of the kernel
+dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8
+pixels in each dimension (a kernel of size 2 pixels in each dimension
+is handled by the bilinear interpolation).  These correspond to the
+options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and
+\code{PS_INTERPOLATE_LANCZOS4}, respectively.
+Given a position on the input image, $(x_0,y_0)$, a kernel is derived
+according to pixels local to the position:
+\begin{equation}
+  h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N)
+\end{equation}
+where
+\begin{eqnarray}
+  \delta x & = & x - x_0 \\
+  \delta y & = & y - y_0 \\
+  {\rm sinc}(z) & = & \sin(z)/z
+\end{eqnarray}
+and $N$ corresponds to the choice of kernel size.  For $N = 2$, the
+kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le
+$).  For $N = 3$, the kernel size is 6 pixels in each dimension
+(i.e., $-3 < \delta x \le 3$).  For $N = 4$, the kernel size is 8
+pixels in each dimension (i.e., $-4 < \delta x \le 4$).
+The interpolated value at the given position, $(x_0,y_0)$, is then
+simply the dot product of the kernel and the fluxes:
+\begin{equation}
+  f(x_0,y_0) = \sum_R f(x,y) h(x,y)
+\end{equation}
+where $R$ is the region defined by the kernel size, and $f(x,y)$ is
+the flux at the pixel position.
+For further information, see the
+\href{http://terapix.iap.fr/IMG/pdf/swarp.pdf}{SWarp manual}.
+\paragraph{Image Cuts and Slices}
+Several functions specify operations which manipulate a collection of
+pixels to return a statistic on the pixel collection.  In the simplest
+case, these are trivial to define: if the boundaries of the region of
+interest are specified along integral pixel coordinates, then the
+pixels used to measure the statistic are always an exact integer.
+This is the case for the function \code{psImageSlice} which requires a
+starting coordinate which is an integer and a width in both dimensions
+which is an integer.  For the case of the functions \code{psImageCut}
+and \code{psImageRadialCut}, the situation is a bit more subtle.  In
+both of these cases, the region is unlikely to contain only whole
+pixels and some choices must be made.
+One posibility which we reject is to identify the fractional pixels
+which are overlapped by the region of interest and add that fraction
+of the pixel's flux when calculating the statistic of interest.  This
+is computationally intensive, and not necessarily well defined for all
+statistics.
+In PSLib, we instead identify the pixels overlapped by the region, use
+the complete set of pixel values, treating all pixels equally, and
+renormalize as needed.  To perform this, the region of interest is
+laid on top of the image pixels.  Any pixels which overlap the region
+are identified as part of the input sample.  The statistic (ie, sample
+mean, robust mode, etc), is then calculated on this collection of
+pixels.  If the output statistic is an average value, the measured
+value is reported.  If the output statistic is a sum value (sum of
+counts, sum of pixels), then the value is renormalized by the ratio of
+pixels used in the calculation to the pixel area of the region of
+interest.  For example, if the sum within a radial aperture is
+requested, the circle of the specified radius and center is placed on
+the pixel grid.  Any pixels which touch the circle are then placed in
+a list to be analysed.  The statistic of interest is the measured for
+this collection of pixels.  In the case of a circular aperture which
+is centered at the coordinate (2,2) and has a radius of 2, the number
+of pixels which are touched by the circle is 16, while the total pixel
+area of the circle is 12.57 square pixels.  In this case, the pixel
+sum is renormalized by the ratio (12.57/16.00).
+\subparagraph{Radial Cuts}
+Consider an image with pixels $x_i,y_i$ and a reference coordinate
+$x_c, y_c$.  We want to construct a radial cut by measuring statistics
+for pixels in a sequence of radial annulii $r_s < r < r_e$.  For each
+annulus, we need to select the pixels which fall within this annulus.
+The coordinates of the center of pixel $i,j$ are $i+0.5,j+0.5$.  A
+given pixel has a distance from the reference coordinate of $dX = x_c
+- i - 0.5, dY = y_c - j - 0.5$.  The pixels to be used for a given
+radial annulus are all of those pixels for which $r_s < \sqrt{dX^2 +
+  dY^2} < r_e$.  This is more efficiently calculated by comparing the
+square of the radii and distances.  All pixels which satisfy the above
+condition are included in a specific annular radius.  All average
+quantities are calculated directly from the pixel ensemble
+statistics.
+\subparagraph{Arbitrary Linear Cuts}
+Select the pixels which lie along a line following steps of 1 pixel
+length:
+\subsection{2D transformations}
+In PSLib, we implement 2-dimensional transformations using
+\code{psPlaneTransform}, which contains a matrix of polynomial
+coefficients for each dimension.  Since we are using these to model
+the real world, where, for example, a particular point on the detector
+maps to a particular point on the sky, we consider only
+transformations that are ``one-to-one''.  This makes it possible to
+speak of inverse transformations, and of combining multiple
+transformations.
+Given a transformation, $f(x,y)$, the inverse transformation,
+$g(x,y)$, is that for which $g(f(x,y)) = (x,y)$ for $(x,y)$ over the
+range of interest (not necessarily the entire set of real numbers).
+Given two transformations, $f(x,y)$ and $g(x,y)$, the combined
+transformation is the transformation, $h(x,y) = g(f(x,y))$ for $(x,y)$
+over the range of interest (not necessarily the entire set of real
+numbers).
+Both of these operations are straightforward if the transformation is
+linear.  If the function $(u,v) = f(x,y)$ is:
+\begin{eqnarray}
+u & = & a + bx + cy \\
+v & = & d + ex + fy
+\end{eqnarray}
+then the inverse transformation $(x,y) = g(u,v)$ is:
+\begin{eqnarray}
+x & = & (-fa+cd)/\Delta + fu/\Delta - cv/\Delta \\
+y & = & (ae-bd)/\Delta - eu/\Delta + bv/\Delta
+\end{eqnarray}
+where $\Delta = bf - ce$ is the matrix determinant.  Given two
+functions $f_i(x,y)$ for $i=1,2$:
+\begin{eqnarray}
+u & = & a_i + b_i x + c_i y \\
+v & = & d_i + e_i x + f_i y
+\end{eqnarray}
+then the combined transformation, $(u,v) = f_2(f_1(x,y))$ is:
+\begin{eqnarray}
+u & = & (a_2 + b_2 a_1 + c_2 d_1) + (b_2 b_1 + c_2 e_1) x + (b_2 c_1 + c_2 f_1) y \\
+v & = & (d_2 + e_2 a_1 + f_2 d_1) + (e_2 b_1 + f_2 e_1) x + (e_2 c_1 + f_2 f_1) y
+\end{eqnarray}
+When the transformations are not linear, the inverse and combined
+transformations can be estimated by sampling a grid over the region of
+interest, calculating the transformation (or double transformation)
+for each sample, and using this information to derive the best fit
+transformation that produces the inverse or combined transformation.
+The inverse transformation should be of the same order as that of the
+forward transformation, while the combined transformation should be of
+the higher order of the two component transformations.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Spherical Rotations with Quaternions}
+\subsubsection{Quaternion Construction}
+The following describes the algorithms needed to implement 3-D
+rotations in terms of quaternions. A quaternion is an ordered set of
+four numbers, $\bar{q} = (q_0, q_1, q_2, q_3)$. A rotation of angle
+$\theta$ about the axis defined by the unit vector $(v_x, v_y, v_z)$
+has quaternion components:
+\begin{eqnarray}
+q_0 & = & v_x sin(\theta/2), \\
+q_1 & = & v_y sin(\theta/2), \\
+q_2 & = & v_z sin(\theta/2), and \\
+q_3 & = & cos(\theta/2). \\
+\end{eqnarray}
+Note that the sine and cosine are taken of the half angle of the
+rotation.  Note also that this implies that the quaternion components
+are normalized such that $|\bar{q}| \def q_0^2 + q_1^2 + q_2^2 + q_3^2
+= 1$.
+The 3-vector representation of the angle of the pole is determined
+from the coordinate of the pole ($\alpha_p, \delta_p$) by:
+\begin{eqnarray}
+v_x & = & \cos \delta_p \cos \alpha_p \\
+v_y & = & \cos \delta_p \sin \alpha_p \\
+v_x & = & \sin \delta_p \\
+\end{eqnarray}
+\subsubsection{Combining Two Rotations}
+Given two quaternions $\bar{p1}$ and $\bar{p2}$, there is a third
+quaternion, $\bar{p}$, which represents the result of first applying
+$\bar{p1}$, and then $\bar{p2}$. The components of $\bar{p}$ are given
+by:
 \begin{verbatim}
+  dX = xe - xs;
+  dY = ye - ys;
+  L = hypot (dX, dY);
+  dX = dX / L;
+  dY = dY / L;
+  REALLOCATE (xvec[0].elements, float, MAX (L, 1));
+  REALLOCATE (yvec[0].elements, float, MAX (L, 1));
+  xvec[0].Nelements = L;
+  yvec[0].Nelements = L;
+  V = (float *)buf[0].matrix.buffer;
+  for (i = 0; i < L; i++) {
+    xi = xs + i*dX - 0.5;
+    yi = ys + i*dY - 0.5;
+    xvec[0].elements[i] = i;
+    yvec[0].elements[i] = V[xi + Nx*yi];
+  }
+p_0 & = &  p2_3 p1_0 + p2_2 p1_1 - p2_1 p1_2 + p2_0 p1_3 \\
+p_1 & = & -p2_2 p1_0 + p2_3 p1_1 + p2_0 p1_2 + p2_1 p1_3 \\
+p_2 & = &  p2_1 p1_0 - p2_0 p1_1 + p2_3 p1_2 + p2_2 p1_3 \\
+p_3 & = & -p2_0 p1_0 - p2_1 p1_1 - p2_2 p1_2 + p2_3 p1_3 \\
 \end{verbatim}
+\paragraph{Image Rotation}
+Image rotation can be performed in two possible ways under different
+circumstances, identified in the following discussion.
+In the simplest case, the rotation angle is an integer multiple of 90
+degrees ($\pi/2$ rad).  In these cases, the input and output pixels
+have a one-to-one mapping.  If the input image has dimensions of $N_x,
+N_y$, then the output image will have dimensions of either $N_x, N_y$
+(for even multiples of 90 degrees) or $N_y, N_x$ (for odd multiples).
+If the angle of the rotation is not a multiple of 90, then the output
+pixels necessarily result from the interpolation of several input
+pixels.  In this case, for an input image of dimensions $N_x, N_y$ and
+rotation angle $\theta$, the output image has dimensions $Lx = |N_x
+\cos \theta| + |N_y \sin \theta|$ and $Ly = |N_x \sin \theta| + |N_y
+\cos \theta|$, each dimension rounded up to the nearest integer as
+needed.  Every pixel in the output image is in general derived from an
+interpolation over 4 neighboring pixels.  The coordinate of a pixel in
+the output image ($i,j$) corresponds to a fractional pixel coordinate
+($x,y$) in the input image according to:
+\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \]
+\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \]
+where the offset coordinate ($i_o,j_o$) depends on the sign of the
+sine of the angle $\theta$.  If the sign of that sine is positive, the
+offset coordinate is ($N_y\sin\theta$,0), otherwise it is
+(0,$-N_x\sin\theta$).
+\subsubsection{Celestial Coordinate Conversions}
+\subsubsection{Rotating a Vector}
+You may rotate a unit vector by first constructing a quaternion
+$\bar{p2}$, whose first three components are the components of the
+unit vector, and whose fourth component is zero. To rotate this vector
+by a quaternion $\bar{p1}$, you apply the formula above for combining
+two quaternions. The rotated vector is found in the first three
+components of the resulting quaternion, $\bar{p}$.
+\subsubsection{Rotation Matrix}
+The rotation matrix representation of a rotation may be derived
+directly from the quaternion representation.  The following formulae
+convert a quaternion to a rotation matrix:
+\begin{eqnarray}
+    rot_{x,x} & = &  q_0 q_0 - q_1 q_1 - q_2 q_2 + q_3 q_3 \\
+    rot_{y,y} & = & -q_0 q_0 + q_1 q_1 - q_2 q_2 + q_3 q_3 \\
+    rot_{z,z} & = & -q_0 q_0 - q_1 q_1 + q_2 q_2 + q_3 q_3 \\
+    rot_{x,y} & = & 2 (q_0 q_1 + q_2 q_3) \\
+    rot_{y,x} & = & 2 (q_0 q_1 - q_2 q_3) \\
+    rot_{x,z} & = & 2 (q_0 q_2 - q_1 q_3) \\
+    rot_{z,x} & = & 2 (q_0 q_2 + q_1 q_3) \\
+    rot_{y,z} & = & 2 (q_1 q_2 + q_0 q_3) \\
+    rot_{z,y} & = & 2 (q_1 q_2 - q_0 q_3) \\
+\end{eqnarray}
+\subsubsection{Conversion to Other Representations}
+You may convert a rotation matrix, m, to a quaternion, p, with the following
+code:
+\begin{verbatim}
+double diag_sum[3];
+int maxi;
+double recip;
+diag_sum[0]=1+m[0][0]-m[1][1]-m[2][2];
+diag_sum[1]=1-m[0][0]+m[1][1]-m[2][2];
+diag_sum[2]=1-m[0][0]-m[1][1]+m[2][2];
+diag_sum[3]=1+m[0][0]+m[1][1]+m[2][2];
+maxi=0;
+for(i=1;i<4;++i) {
+    if(diag_sum[i]>diag_sum[maxi]) maxi=i;
+}
+p[maxi]=0.5*sqrt(diag_sum[maxi]);
+recip=1./(4.*p[maxi]);
+if(maxi==0) {
+    p[1]=recip*(m[0][1]+m[1][0]);
+    p[2]=recip*(m[2][0]+m[0][2]);
+    p[3]=recip*(m[1][2]-m[2][1]);
+} else if(maxi==1) {
+    p[0]=recip*(m[0][1]+m[1][0]);
+    p[2]=recip*(m[1][2]+m[2][1]);
+    p[3]=recip*(m[2][0]-m[0][2]);
+} else if(maxi==2) {
+    p[0]=recip*(m[2][0]+m[0][2]);
+    p[1]=recip*(m[1][2]+m[2][1]);
+    p[3]=recip*(m[0][1]-m[1][0]);
+} else if(maxi==3) {
+    p[0]=recip*(m[1][2]-m[2][1]);
+    p[1]=recip*(m[2][0]-m[0][2]);
+    p[2]=recip*(m[0][1]-m[1][0]);
+}
+\end{verbatim}
+\subsection{Celestial Coordinate Conversions}
 Changes between spherical coordinate systems (ie, Ecliptic, Galactic,
 …
 the forward transformation.
 \paragraph{Galactic to ICRS}
+\subsubsection{Galactic to ICRS}
 The appropriate values, from the Hipparcos and Tycho Catalogues are:
 …
 \end{eqnarray}
 \paragraph{Ecliptic to ICRS}
+\subsubsection{Ecliptic to ICRS}
 The appropriate values, from Zombeck, are:
 …
 where $T$ is the time in Julian centuries since 1900.
 \paragraph{Precession}
+\subsubsection{Precession}
 The appropriate values, from Elixir, are:
 …
 \paragraph{Suggested test cases}
+\subsubsection{Suggested test cases}
 $(\alpha,\delta) = (0^\circ,0^\circ)$ transforms to Galactic
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{2D transformations}
+In PSLib, we implement 2-dimensional transformations using
+\code{psPlaneTransform}, which contains a matrix of polynomial
+coefficients for each dimension.  Since we are using these to model
+the real world, where, for example, a particular point on the detector
+maps to a particular point on the sky, we consider only
+transformations that are ``one-to-one''.  This makes it possible to
+speak of inverse transformations, and of combining multiple
+transformations.
+Given a transformation, $f(x,y)$, the inverse transformation,
+$g(x,y)$, is that for which $g(f(x,y)) = (x,y)$ for $(x,y)$ over the
+range of interest (not necessarily the entire set of real numbers).
+Given two transformations, $f(x,y)$ and $g(x,y)$, the combined
+transformation is the transformation, $h(x,y) = g(f(x,y))$ for $(x,y)$
+over the range of interest (not necessarily the entire set of real
+numbers).
+Both of these operations are straightforward if the transformation is
+linear.  If the function $(u,v) = f(x,y)$ is:
+\begin{eqnarray}
+u & = & a + bx + cy \\
+v & = & d + ex + fy
+\end{eqnarray}
+then the inverse transformation $(x,y) = g(u,v)$ is:
+\begin{eqnarray}
+x & = & (-fa+cd)/\Delta + fu/\Delta - cv/\Delta \\
+y & = & (ae-bd)/\Delta - eu/\Delta + bv/\Delta
+\end{eqnarray}
+where $\Delta = bf - ce$ is the matrix determinant.  Given two
+functions $f_i(x,y)$ for $i=1,2$:
+\begin{eqnarray}
+u & = & a_i + b_i x + c_i y \\
+v & = & d_i + e_i x + f_i y
+\end{eqnarray}
+then the combined transformation, $(u,v) = f_2(f_1(x,y))$ is:
+\begin{eqnarray}
+u & = & (a_2 + b_2 a_1 + c_2 d_1) + (b_2 b_1 + c_2 e_1) x + (b_2 c_1 + c_2 f_1) y \\
+v & = & (d_2 + e_2 a_1 + f_2 d_1) + (e_2 b_1 + f_2 e_1) x + (e_2 c_1 + f_2 f_1) y
+\end{eqnarray}
+When the transformations are not linear, the inverse and combined
+transformations can be estimated by sampling a grid over the region of
+interest, calculating the transformation (or double transformation)
+for each sample, and using this information to derive the best fit
+transformation that produces the inverse or combined transformation.
+The inverse transformation should be of the same order as that of the
+forward transformation, while the combined transformation should be of
+the higher order of the two component transformations.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{Projections}
+We implement three types of projections: {\em zenithal}, {\em
+cylindrical} and {\em pseudocylindrical}, each requiring slightly
+different handling.  Our representations are based on the treatment of
+projections presented by
+\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \&
+Calabretta (1995, ADASS, 4, 233)}.  In all of these projections, we
+are converting from a spherical coordinate $\alpha,\delta$ to a linear
+(2-D) coordinate $x_p,y_p$.  The projection is defined by the
+projection type, the projection center ($\alpha_p, \delta_p$) and the
+the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$).
+In the structure, \code{psProjection}, the projection type is defined
+by the element \code{type}, the projection center $\alpha_p,\delta_p$
+is defined by the elements \code{R,D}, and the plate scales,
+$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}.  The plate
+scales are applied independently to the $x$ and $y$ coordinates to
+convert them to the corresponding linear units (ie, pixels):
+%
+\begin{eqnarray}
+x_p & = & \rho_x x \\
+y_p & = & \rho_y y \\
+\end{eqnarray}
+%
+In the discussions below, we ignore this last step (or first step,
+depending on the direction of the conversion).
+\paragraph{Zenithal Projections}
+The {\em zenithal} projections are defined relative to a set of
+spherical coordinates with pole at the center of the projection
+($\alpha_p, \delta_p$), and which thus represents a coordinate system
+rotated relative to the coordinate system of $\alpha, \delta$.  In
+this spherical coordinate system, the coordinate of longitude is
+labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the
+latitude, measured from the pole, is labeled $\theta$ and has domain
+$0 \le \theta \le \pi$.  The coordinate frame of $\phi,\theta$ is
+defined so that $\phi_p$, the longitude of the target system pole, is
+.0.
+For an arbitrary projection center, it is necessary to convert the
+spherical coordinates to be projected ($\alpha,\delta$) to the
+projection spherical coordinate system coordinates ($\phi, \theta$).
+In practice, we construct the following useful trigonometric
+relationships between $\phi$ and $\theta$ which may be employed in the
+equations of $x,y$ below:
+%
+\begin{eqnarray}
+\sin \theta           & = & \sin \delta \sin \delta_p + \cos \delta \cos \delta_p \cos (\alpha - \alpha_p) \\
+\cos \theta \cos \phi & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \cos (\alpha - \alpha_p) \\
+\cos \theta \sin \phi & = & - \cos \delta \sin (\alpha - \alpha_p)
+\end{eqnarray}
+%
+For the inverse transformations, the equivalent relationships are:
+%
+\begin{eqnarray}
+\sin \delta                          & = & \sin \theta \sin \delta_p + \cos \theta \cos \delta_p \cos \phi \\
+\cos \delta \cos (\alpha - \alpha_p) & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \cos \phi \\
+\cos \delta \sin (\alpha - \alpha_p) & = & - \cos \theta \sin \phi
+\end{eqnarray}
+%
+For zenithal projections, the linear coordinates are related to
+$\phi,\theta$ by:
+%
+\begin{eqnarray}
+x & = & R_\theta \sin \phi \\
+y & = & -R_\theta \cos \phi
+\end{eqnarray}
+%
+and the inverse:
+%
+\begin{eqnarray}
+R_\theta & = & \sqrt{x^2 + y^2} \\
+\phi     & = & {\rm atan} (-y,x)
+\end{eqnarray}
+%
+The coordinates $x,y$ above are defined to be in angular units (ie,
+radians).
+From these relationships, we can calculate $\alpha, \delta$ as:
+%
+\begin{eqnarray}
+\alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\
+\delta            & = & \arcsin (\sin \delta) \\
+\end{eqnarray}
+%
+Note that if $(x,y) = (0,0)$, then $\alpha = \alpha_p, \delta = \delta_p$.
+\subparagraph{Gnomonic}
+The Gnomonic projection (``TAN'') is a zenithal projection with
+$R_\theta = \cot \theta$.  The resulting relationships for $(x,y)$ and
+for $\sin \theta, \cos \theta$ are:
+\begin{eqnarray}
+x           & = & \frac{\cos \theta \sin \phi}{\sin \theta} \\
+y           & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\
+\sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\
+\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\
+\end{eqnarray}
+where $\zeta = 1 / R_\theta$.
+\subparagraph{Orthographic}
+The Orthographic projection (``SIN'') is a zenithal projection with
+$R_\theta = \cos \theta$.  The resulting relationships for $(x,y)$ and
+for $\sin \theta, \cos \theta$ are:
+\begin{eqnarray}
+x           & = & \cos \theta \sin \phi \\
+y           & = & -\cos \theta \cos \phi \\
+\sin \theta & = & \sqrt{1 - R_\theta^2} \\
+\cos \theta & = & R_\theta \\
+\end{eqnarray}
+\paragraph{Cylindrical and Pseudocylindrical Projections}
+The {\em cylindrical} and {\em pseudocylindrical} projections are
+defined relative to a set of cylindrical coordinates whose pole is
+coincident with the pole of the spherical coordinates.  These
+projections are particularly used for full-sky representations, and
+are only defined for projection centers with $\delta_p = 0$.  In this
+spherical coordinate system, the coordinate of longitude is labeled
+$\phi$, and has domain of $-\pi < \phi \le \pi$, while the latitude,
+measured from the pole, is labeled $\theta$ and has domain $0 \le
+\theta \le \pi$.  The projection center longitude, $\alpha_p$
+corresponds to $\phi = 0$, thus the value of $\phi$ is determined as
+$\alpha - \alpha_p$ for all such projections.
+\subparagraph{Cartesian}
+The Cartesian projection (``CAR'') is a very simple cylindrical
+projection with the following relationships between $x,y$ and
+$\phi,\theta$:
+\begin{eqnarray}
+x & = & \phi \\
+y & = & \theta
+\end{eqnarray}
+\subparagraph{Mercator}
+The Mercator projection (``MER'') is a cylindrical projection.
+\begin{eqnarray}
+x & = & \phi \\
+y & = & \ln \left( \tan (\pi/4 + \theta/2) \right) \\
+{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^y \right) - \pi/2
+\end{eqnarray}
+\subparagraph{Hammer-Aitoff}
+The Hammer-Aitoff projection(``AIT'') is a pseudocylindrical projection, and is defined:
+\begin{eqnarray}
+x & = & 2 \zeta \cos \theta \sin \frac{\phi}{2} \\
+y & = & \zeta \sin \theta \\
+{\rm where}\hspace{1cm} \zeta^{-1} & \equiv & \sqrt{\frac{1}{2}\left(1 + \cos \theta \cos \frac{\phi}{2} \right)}
+\end{eqnarray}
+And in reverse:
+\begin{eqnarray}
+\phi & = & 2 {\rm \arctan} (2z^2 - 1, x z) \\
+\theta & = & \arcsin (yz) \\
+{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x/2)^2 - y^2}
+\end{eqnarray}
+\subparagraph{Parabolic}
+The Parabolic projection (``PAR'') is a pseudocylindrical projection, and is defined:
+\begin{eqnarray}
+x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\
+y & = & \pi \sin \frac{\theta}{3} \\
+\end{eqnarray}
+And in reverse:
+\begin{eqnarray}
+\theta & = & 3 \sin^{-1} \rho \\
+\phi   & = & \frac{x}{1 - 4\rho^2} \\
+{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\
+\end{eqnarray}
+\subsubsection{Offset}
+Coordinate offsets can be either spherical offsets or linear offsets.
+A spherical offset is performed by adding the components of the
+offset, after unit conversion, to the given position.  The resulting
+coordinates must be wrapped to within the allowed range ($-\pi$ to
+$\pi$, 0 to $2\pi$).
+A linear offset is defined to be a linear offset in a tangent
+projection centered on the starting coordinate with $y$ axis aligned
+with the local direction or increasing Declination.  This projection
+is undefined only for the coordinates exactly at the north and south
+poles, in which case the orientation is defined to have the $y$ axis
+parallel to the line of RA = 0.0.  The scale of the projection is 1.0
+(ie, 1 'pixel' is 1 radian) and the given offsets must the scaled
+based on the given offset units.
+Pseudo-code to implement the above for an offset:
+\begin{verbatim}
+psSphere *psSphereSetOffset (psSphere pos, psSphere offset) {
+  psPlane lin;
+  psSphere new;
+  psProjection proj;
+  proj.R = pos->r;
+  proj.D = pos->d;
+  proj.X = 0;
+  proj.Y = 0;
+  proj.type = PS_PROJ_TAN;
+  lin.x = offset.r;
+  lin.y = offset.d;
+  new = psDeproject (&lin, &proj);
+  return (new);
+}
+\end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{Tangent Plane to Sky}
+\subsection{Tangent Plane to Sky}
 \tbd{we will replace the SLALIB version of AOPPA with a new function}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{The One-to-Many Problem with Mosaic Cameras}
+\subsection{Sky to Tangent Plane (II)}
+This section describes the transformation between celestial coordinates
+(R.A., Dec.) and local terrestrial coordinates (Az, Alt). This transformation
+is broken down into a number of steps as described below.
+\paragraph{Reference Implementations}
+There are two reference implementatins for the code to account for the
+motion of the Earth in space. The first are the sample routines
+provided by the IERS to accompany chaper 5 of IERS Bulletin 32.  This
+document and the code can be downloaded from
+http://maia.usno.navy.mil/conv2003.html .  The second reference
+implementation is the SOFA software package managed by the IAU and
+available at http://www.iau-sofa.rl.ac.uk Only the 2003-04-29 version
+of SOFA should be used.  The IERS code requires a few of the rotation
+matrix utility routines from SOFA.
+Both implementations are in FORTRAN 77. The SOFA code has a more
+complex implementation of precession-nutation for backward
+compatibility with the pre 2003-01-01 conventions.  The IERS code
+includes some tricks to achieve greater precision in the fundamental
+arguments of nutation, which the SOFA code omits.  Therefore, the main
+reference for psLib should be the IERS code.  Note that the IERS code
+calculates the transform from terrestrial to celestial coordinates,
+while the SOFA code calculates its inverse.
+\subsubsection{Coordinate Systems}
+\begin{figure}
+\psfig{file=transforms.ps}
+\caption{Coordinates systems and the transformations between them}
+\end{figure}
+Figure X shows the transformation steps and intermediate coordinate systems
+between celestial and local terrestrial coordinate systems. The intermediate
+coordinate systems are defined below.
+\paragraph{ICRS}
+The official IAU-sanctioned celestial coordinate system is the
+International Celestial Reference System (ICRS). It is defined in terms of
+a number of radio sources whose positions have been measured using VLBI.
+It can be tied to the optical through the Hipparcos catalog. The ICRS has its
+origin at the solar system barycenter.
+\paragraph{GCRS}
+The Geocentric Celestial Reference System (GCRS) corresponds to the ICRS, but
+has its origin at the center of the Earth. The differences between the two
+systems are due to the velocity of the Earth (aberration), the position of
+the Earth (parallax), and general relativistic bending of light rays.
+There is no net rotation between the ICRS and the GCRS.
+\paragraph{ITRS}
+The International Terrestrial Reference System (ITRS) is a coordinate
+system which is fixed with respect to the Earth's crust.
+\paragraph{Intermediate Coordinate Systems - CIP, CEO, TEO}
+The transform between the GCRS and ITRS is conventionally
+decomposed into three parts in order to isolate the relatively rapid rotation
+of the Earth from the movement of the Earth's rotational axis in the GCRS
+and ITRS. All three sub-transforms are rigid rotations.
+This decomposition results in two intermediate coordinate systems. Both of
+these share the same pole, known as the Celestial Intermediate Pole (CIP).
+The CIP is defined by its motion in the GCRS to match the Tisserand
+mean axis of the Earth (Seidelmann 1982, Celesial Mechanics 27, 78-106),
+excluding motions with periods less than or equal
+to two days. The CIP approximates the angular momentum vector of the
+rotating Earth.
+The X axes of the intermediate coordinate systems are known as the
+Celestial and Terrestrial Ephemeris Origins. (CEO and TEO). Both are defined
+to be non-rotating origins. A non-rotating origin is a point on the equator
+whose instantaneous motion is always orthogonal to the equator
+(Kaplan 2003 IAU XXV Joint Discussion 16
+\footnote{http://aa.usno.navy.mil/kaplan/NROs\%5BJD16proc\%5D.pdf}).
+Thus the CEO is defined by its position in the GCRS at some epoch and by the
+motion of the CIP in the GCRS since that date. Similarly the TEO is
+defined by its position in the ITRS at some epoch and the motion of the
+CIP in the ITRS since that date.
+\subsubsection{ICRS - GCRS}
+The transformation between barycentric (ICRS) and geocentric (GCRS) coordinates
+involves two components. These are
+the general relativistic deflection of light rays by the Sun's gravity, and
+aberration, due to the orbital motion
+of the Earth.
+\paragraph{Gravitational Deflection}
+The Sun's gravity bends the path of light rays which pass near it.
+To first order, a light ray is deflected by an angle of $4GM/c^2r_0$ radians,
+where $G$ is the gravitational constant,
+$M$ is the mass of the Sun,
+$c$ is the speed of light, and
+$r_0$ is the point of closest approach to the light ray to the Sun.
+To the same order this is equal to the impact parameter - i.e. the point
+of closest approach if the light ray were not deflected. Note that
+$r_0/d = \tan(\theta)$, where $d$ is the distance from the Earth
+to the Sun, and $\theta$
+is the angular separation of the star from the center of the Sun.
+There is a maximum deflection of 1.75 arc seconds if we set
+$r_0$ to the radius of the sun.
+Since the Sun bends light rays toward it, a star appears shifted away from the sun in the sky.
+\paragraph{Aberration}
+Aberration is the apparent change in direction of a ray of light in the
+reference frame of a moving observer. Traditionally the aberration
+calculation has been done with a linear expansion of the full
+relativistic expression, often neglecting all but the linear term in
+$v/c$, since the relativistic terms are on the order of a miliarcsecond.
+However, the full relativistic expression poses no challenge for modern
+computers, so psLib will use the following procedure to calculate aberration.
+Suppose an observer has a velocity $\beta\hat{\beta}$, with respect to
+the Solar System barycenter, where $\beta$ is in units of the speed of
+light, and $\hat{\beta}$ is a unit vector. Suppose also that the unit vector
+$\hat{r}$ points toward a star in the barycenter frame of reference
+(i.e. the ``actual'' position).
+and $\hat{r}'$ gives the direction of the star in the observer's frame,
+(i.e. the apparent position).
+First, decompose $\hat{r}$ into components parallel and perpendicular to
+$\hat{\beta}$ by calculating
+$\mu = \hat{r}\cdot\hat{\beta}$ and
+$\vec{r}_\perp = \hat{r} - \mu \hat{\beta}$.
+Next, use the following expression for relativistic beaming, modified
+slightly from equation 4.8b of Rybicki and Lightman:
+\begin{equation}
+\mu' = \mu + \beta \frac{\mu^2 - 1}{1 - \beta\mu}
+\end{equation}
+where $\mu' = \hat{r}' \cdot \hat{\beta}$.
+Now, the component of $\hat{r}'$ perpendicular to $\hat{\beta}$
+(i.e. $\vec{r}_\perp'$) must point
+in the same direction as $\vec{r}_\perp$, but will have a different magnitude
+because $\hat{r}'$ is a unit vector. In other words,
+$\vec{r}_\perp' = a\vec{r}_\perp$, for some scalar $a$. So the next step is
+to calculate $a = \sqrt{(1-\mu'^2)/\vec{r}_\perp}$.
+Finally, reassemble the components of
+$\hat{r}' = \mu'\hat{\beta} + a \vec{r_\perp}$.
+\subsubsection{GCRS - ITRS}
+The transformation between geocentric celestial coordinates and terrestrial
+coordinates is a solid body rotation due to the motion of the Earth is space.
+This is conventionally broken down into three components to isolate the
+relatively rapid rotation of the Earth from the motion of its rotational axis.
+This section is largely a summary of
+Chapter 5 of IERS Technical Note 32 \footnote{http://maia.usno.navy.mil/conv2003.html}
+(hereafter IERS32),
+which is a description of the implementation of the Resoltions of the
+XXIVth General Assembly of the IAU, available from the same URL as above.
+These two documents describe a set of conventions which have been in effect
+since 2003-01-01. The conventions in effect before that date will not be
+implemented by psLib.
+\paragraph{Precession/Nutation}
+The transform between the GCRS and the CIP/CEO coordinate systems is described
+by the IAU 2000A precession-nutation model, which is accurate to the
+.2 mas level.
+For higher accuracy the user must apply corrections to the model, which are tabulated by the IERS.
+The IAU 2000A precession-nutation model may be calculated in the following
+way. First calculate the time $t$ as the number of Julian centuries since
+-01-01T12:00:00 TT.
+Next calculate the fundamental arguments of nutation using equations (40)
+and (41) of IERS32, reproduced below:
+\begin{eqnarray}
+F_1\equiv l\quad  =~&\ Mean\ Anomaly\ of\ the\ Moon \cr
+ =~& 134.96340251^\circ + 1717915923.2178'' t
+ + 31.8792'' t^2 + 0.051635'' t^3 - 0.00024470'' t^4,\cr
+F_2\equiv l'\quad =~&\ Mean\ Anomaly\ of\ the\ Sun\cr
+=~& 357.52910918^\circ + 129596581.0481'' t
+- 0.5532'' t^2 + 0.000136'' t^3 - 0.00001149'' t^4,\cr
+F_3\equiv F\quad  =~& L - \Omega\cr
+=~& 93.27209062^\circ + 1739527262.8478'' t - 12.7512'' t^2
+- 0.001037'' t^3 + 0.00000417'' t^4,\cr
+F_4\equiv D\quad  =~&\ Mean\ Elongation\ of\ the\ Moon\ from\ the\ Sun\cr
+=~& 297.85019547^\circ + 1602961601.2090'' t - 6.3706'' t^2
++ 0.006593'' t^3 - 0.00003169'' t^4,\cr
+F_5\equiv\Omega\quad  =~&\ Mean\ Longitude\ of\ the\ Ascending\ Node\ of\
+the\ Moon\cr
+=~& 125.04455501^\circ - 6962890.5431'' t + 7.4722'' t^2 + 0.007702'' t^3 - 0.00005939'' t^4 \cr
+F_6\ \equiv l_{Me}\quad    =~& 4.402 608 842 + 2608.7903 141 574\times t,\cr
+F_7\ \equiv l_{Ve}\quad    =~& 3.176 146 697 + 1021.3285 546 211 \times t,\cr
+F_8\ \equiv l_{E\ }\quad   =~& 1.753 470 314 + 628.3075 849 991 \times t,\cr
+F_9\equiv l_{Ma}\quad    =~& 6.203 480 913 + 334.0612 426 700 \times t,\cr
+F_{10}\equiv l_{Ju}\quad =~& 0.599 546 497 + 52.9690 962 641 \times t,\cr
+F_{11}\equiv l_{Sa}\quad =~& 0.874 016 757 + 21.3299 104 960 \times t,\cr
+F_{12}\equiv l_{Ur}\quad =~& 5.481 293 872 +  7.4781 598 567 \times t,\cr
+F_{13}\equiv l_{Ne}\quad =~& 5.311 886 287 +  3.8133 035 638 \times t,\cr
+F_{14}\equiv p_{a\ }\quad =~& 0.024 381 750 \times t + 0.000 005 386 91 \times t^2.
+\end{eqnarray}
+Next calculate the quantities $X$, $Y$, and $s$, using expressions of the form:
+\begin{equation}
+     \sum_{j} p_j t^j + \sum_{j}\sum_{i}[
+     (a_{{\rm s},j})_i t^j \sin ({\rm \scriptstyle {ARG_{i,j}}})
+   + (a_{{\rm c},j})_i t^j \cos ({\rm \scriptstyle {ARG_{i,j}}})]
+   ,
+\end{equation}
+where the $\rm \scriptstyle{ARG_{i,j}} = \sum_{k} w_{i,j,k} F_k$ represent linear
+combinations of the fundamental arguments of nutation.
+The constants $p_j$, $w_{i,j,k}$, $(a_{{\rm s},j})_i$, and $(a_{{\rm c},j})_i$
+are given in the ASCII files:
+tab5.2a.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2a.txt} (for $X$),
+tab5.2b.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2b.txt} (for $Y$), and
+tab5.2c.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2c.txt} (for $s+XY/2$).
+Note that the expansion is given for $s+XY/2$, since this series converges
+more rapidly than the one for $s$ alone.
+Each file contains a human-readable header, which includes the polynomial
+coeficients, $p_j$ under the heading ``Polynomial part''. The data part of the
+file lists the remaining constants, with rows cycling first through $i$, and
+then through $j$. There is a separate heading each time $j$ increments.
+Each row contains the following columns:
+\begin{itemize}
+\item col 1 - A running index of rows in the table.
+\item col 2 - The sine coeficients, $(a_{{\rm s},j})_i$
+\item col 3 - The cosine coeficients, $(a_{{\rm c},j})_i$
+\item cols 4 - 17 The weighting factors for the fundamental arguments of
+                  nutation, $w_{i,j,k}$.
+\end{itemize}
+A FORTRAN reference implementation for the precession/nutation model is available from the IERS
+\footnote{http://maia.usno.navy.mil/conv2000/chapter5/XYS2000A.f}.
+The psLib results should agree with the reference implementation to within
+the limits of numerical precision.
+Next, corrections to $X$, and $Y$ may be obtained from the IERS as part of
+Bulletin A, or B. It is recommended to use the values published daily in
+http://maia.usno.navy.mil/ser7/finals2000A.daily, which has the format
+described by http://maia.usno.navy.mil/ser7/readme.finals2000A. The
+quantities of interest are labeled dX and dY. Note that UT1$-$UTC and the
+polar motion values are obtained from this same table.
+By convention, nutation terms with periods of less
+than two days
+are accounted for by the corresponding polar motion. So it is sufficient to
+interpolate the corrections tabulated daily by the IERS, and take the result as
+instantaneous values.
+The final step is to use $X$, $Y$, and $s$ to calculate the rotation
+matrix from the CIP/CEO system to the GCRS using IERS32 equation (10),
+reproduced below:
+\begin{equation}
+\begin{pmatrix}1-aX^2& -aXY& X\cr -aXY& 1-aY^2& Y\cr -X& -Y&
+-a(X^2+Y^2)\cr
+\end{pmatrix} \cdot R_3(s),
+\end{equation}
+where $R_3$ denotes a rotation about the Z axis,
+$a = 1/(1+\sqrt{1 - X^2 + Y^2})$,
+and $X$ and $Y$ are expressed in radians.
+A FORTRAN reference implementation for this calculation is given
+by the IERS \footnote{http://maia.usno.navy.mil/conv2000/chapter5/BPN2000.f}.
+Note that above we gave the expression for the transform toward celestial
+coordinates (upward in figure X), in order to match the IERS reference code.
+The inverse transform may be found by inverting the resulting rotation.
+\paragraph{Rotation of the Earth}
+The transform from the CIP/CEO to CIP/TEO coordinate systems is a
+rotation about the CIP (i.e. the Z axis) by an angle known as the
+``Earth Rotation Angle''.
+By definition the Earth Rotation Angle is given by
+equation (13) of IERS32, reproduced below:
+\begin{equation}
+\theta(T_u)=2\pi(0.7790572732640 + 1.00273781191135448T_u),
+\end{equation}
+where $T_u$ is the Julian UT1 date minus 2451545.0 .
+\paragraph{Polar Motion}
+The motion of the CIP in the ITRS is known as ``polar motion''. Similarly to
+precession/nutation, the instantaneous position of the CIP in the
+ITRS is specified by the quantites $x_p$, and $y_p$, and a third quantity,
+$s'$, gives the position of the TEO with respect to the ITRS.
+The values of $x_p$ and $y_p$ are published daily by the IERS in
+http://maia.usno.navy.mil/ser7/finals2000A.daily, which has the format
+described by http://maia.usno.navy.mil/ser7/readme.finals2000A.
+The UT1$-$UTC, and the precession/nutation corrections (discussed elsewhere
+in this document) come from this same source.
+The polar motion coordinates should be interpolated using a third order
+polynomial, as described in
+IERS Gazette \#13 \footnote{http://maia.usno.navy.mil/iers-gaz13},
+which gives a
+FORTRAN reference implementation of the correct procedure.
+\tbd{reference to interpolation in this doc?}
+The values published by the IERS are smoothed to remove noise and
+variations on the timescale of a day or less. There are two sources of
+short timescale variations - tidal effects on the order of 0.1 milliarcseconds,
+and short period nutation terms on the order of 15 microarcseconds.
+Both of these effects may be modeled and added to the interpolated values
+for higher accuracy.
+The tidal effects should be included using the FORTRAN reference implementation
+of the Ray tidal model given in IERS Gazette \#13. This code should be
+mimiced to machine accuracy by psLib.
+By definition of the CIP, nutation terms with periods less than 2 days are
+not included in the IAU 2000A precession/nutation model.
+So these motions
+must be compensated for by their equivalent polar motions. These may
+be calculated using a form similar to that of the precession/nutation $X$,
+and $Y$. The constants to use are given in Table 5.1 of IERS32.
+Note that only the terms with periods less than 2 days should be used.
+The quantity $s'$ may be approximated with microarcsecond accuracy over this
+century by $s' = -4.7 \times 10^{-5} t$ in arcseconds. There is no need
+to apply short timescale corrections to $s'$.
+The transform from the ITRS to the CIP/TEO frame can be constructed by
+first rotating about the X axis by $y_p$, then rotating about the X axis by
+$x_p$, and finally rotating about the Z axis by $s'$.
+The IERS reference implementation for this is given in the subroutine
+POM2000 \footnote{http://maia.usno.navy.mil/conv2000/chapter5/POM2000.f}.
+Note that we describe the transform toward celestial coordinates (upward in
+figure X), in order to match the reference implementation.
+\subsubsection{ITRS - Alt/Az}
+\paragraph{Orientation of the Observer}
+An observer's astronomical longitude and latitude give the orientation of
+the local vertical with respect to the ITRS. Note that these coordinates
+can be approximated by the geographic longitude and latitude of the observatory,
+but their exact values must be calibrated from observation of stars
+with known coordinates in the ICRS.
+The transform from the ITRS to Az/Alt in the absence of atmospheric refraction
+is first a rotation about the Z axis by the observer's astronomical longitude,
+and then a rotation about the Y axis of 90 degrees minus the observer's
+astronomical latitude, followed by a rotation about the Z axis of 180 degrees
+so that North is zero azimuth.
+\paragraph{Atmospheric Refraction}
+\tbd{add in summary of Ken's paper}
+\subsection{Projections}
+We implement three types of projections: {\em zenithal}, {\em
+cylindrical} and {\em pseudocylindrical}, each requiring slightly
+different handling.  Our representations are based on the treatment of
+projections presented by
+\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \&
+Calabretta (1995, ADASS, 4, 233)}.  In all of these projections, we
+are converting from a spherical coordinate $\alpha,\delta$ to a linear
+(2-D) coordinate $x_p,y_p$.  The projection is defined by the
+projection type, the projection center ($\alpha_p, \delta_p$) and the
+the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$).
+In the structure, \code{psProjection}, the projection type is defined
+by the element \code{type}, the projection center $\alpha_p,\delta_p$
+is defined by the elements \code{R,D}, and the plate scales,
+$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}.  The plate
+scales are applied independently to the $x$ and $y$ coordinates to
+convert them to the corresponding linear units (ie, pixels):
+%
+\begin{eqnarray}
+x_p & = & \rho_x x \\
+y_p & = & \rho_y y \\
+\end{eqnarray}
+%
+In the discussions below, we ignore this last step (or first step,
+depending on the direction of the conversion).
+\subsubsection{Zenithal Projections}
+The {\em zenithal} projections are defined relative to a set of
+spherical coordinates with pole at the center of the projection
+($\alpha_p, \delta_p$), and which thus represents a coordinate system
+rotated relative to the coordinate system of $\alpha, \delta$.  In
+this spherical coordinate system, the coordinate of longitude is
+labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the
+latitude, measured from the pole, is labeled $\theta$ and has domain
+$0 \le \theta \le \pi$.  The coordinate frame of $\phi,\theta$ is
+defined so that $\phi_p$, the longitude of the target system pole, is
+.0.
+For an arbitrary projection center, it is necessary to convert the
+spherical coordinates to be projected ($\alpha,\delta$) to the
+projection spherical coordinate system coordinates ($\phi, \theta$).
+In practice, we construct the following useful trigonometric
+relationships between $\phi$ and $\theta$ which may be employed in the
+equations of $x,y$ below:
+%
+\begin{eqnarray}
+\sin \theta           & = & \sin \delta \sin \delta_p + \cos \delta \cos \delta_p \cos (\alpha - \alpha_p) \\
+\cos \theta \cos \phi & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \cos (\alpha - \alpha_p) \\
+\cos \theta \sin \phi & = & - \cos \delta \sin (\alpha - \alpha_p)
+\end{eqnarray}
+%
+For the inverse transformations, the equivalent relationships are:
+%
+\begin{eqnarray}
+\sin \delta                          & = & \sin \theta \sin \delta_p + \cos \theta \cos \delta_p \cos \phi \\
+\cos \delta \cos (\alpha - \alpha_p) & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \cos \phi \\
+\cos \delta \sin (\alpha - \alpha_p) & = & - \cos \theta \sin \phi
+\end{eqnarray}
+%
+For zenithal projections, the linear coordinates are related to
+$\phi,\theta$ by:
+%
+\begin{eqnarray}
+x & = & R_\theta \sin \phi \\
+y & = & -R_\theta \cos \phi
+\end{eqnarray}
+%
+and the inverse:
+%
+\begin{eqnarray}
+R_\theta & = & \sqrt{x^2 + y^2} \\
+\phi     & = & {\rm atan} (-y,x)
+\end{eqnarray}
+%
+The coordinates $x,y$ above are defined to be in angular units (ie,
+radians).
+From these relationships, we can calculate $\alpha, \delta$ as:
+%
+\begin{eqnarray}
+\alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\
+\delta            & = & \arcsin (\sin \delta) \\
+\end{eqnarray}
+%
+Note that if $(x,y) = (0,0)$, then $\alpha = \alpha_p, \delta = \delta_p$.
+\paragraph{Gnomonic}
+The Gnomonic projection (``TAN'') is a zenithal projection with
+$R_\theta = \cot \theta$.  The resulting relationships for $(x,y)$ and
+for $\sin \theta, \cos \theta$ are:
+\begin{eqnarray}
+x           & = & \frac{\cos \theta \sin \phi}{\sin \theta} \\
+y           & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\
+\sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\
+\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\
+\end{eqnarray}
+where $\zeta = 1 / R_\theta$.
+\paragraph{Orthographic}
+The Orthographic projection (``SIN'') is a zenithal projection with
+$R_\theta = \cos \theta$.  The resulting relationships for $(x,y)$ and
+for $\sin \theta, \cos \theta$ are:
+\begin{eqnarray}
+x           & = & \cos \theta \sin \phi \\
+y           & = & -\cos \theta \cos \phi \\
+\sin \theta & = & \sqrt{1 - R_\theta^2} \\
+\cos \theta & = & R_\theta \\
+\end{eqnarray}
+\subsubsection{Cylindrical and Pseudocylindrical Projections}
+The {\em cylindrical} and {\em pseudocylindrical} projections are
+defined relative to a set of cylindrical coordinates whose pole is
+coincident with the pole of the spherical coordinates.  These
+projections are particularly used for full-sky representations, and
+are only defined for projection centers with $\delta_p = 0$.  In this
+spherical coordinate system, the coordinate of longitude is labeled
+$\phi$, and has domain of $-\pi < \phi \le \pi$, while the latitude,
+measured from the pole, is labeled $\theta$ and has domain $0 \le
+\theta \le \pi$.  The projection center longitude, $\alpha_p$
+corresponds to $\phi = 0$, thus the value of $\phi$ is determined as
+$\alpha - \alpha_p$ for all such projections.
+\paragraph{Cartesian}
+The Cartesian projection (``CAR'') is a very simple cylindrical
+projection with the following relationships between $x,y$ and
+$\phi,\theta$:
+\begin{eqnarray}
+x & = & \phi \\
+y & = & \theta
+\end{eqnarray}
+\paragraph{Mercator}
+The Mercator projection (``MER'') is a cylindrical projection.
+\begin{eqnarray}
+x & = & \phi \\
+y & = & \ln \left( \tan (\pi/4 + \theta/2) \right) \\
+{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^y \right) - \pi/2
+\end{eqnarray}
+\paragraph{Hammer-Aitoff}
+The Hammer-Aitoff projection(``AIT'') is a pseudocylindrical projection, and is defined:
+\begin{eqnarray}
+x & = & 2 \zeta \cos \theta \sin \frac{\phi}{2} \\
+y & = & \zeta \sin \theta \\
+{\rm where}\hspace{1cm} \zeta^{-1} & \equiv & \sqrt{\frac{1}{2}\left(1 + \cos \theta \cos \frac{\phi}{2} \right)}
+\end{eqnarray}
+And in reverse:
+\begin{eqnarray}
+\phi & = & 2 {\rm \arctan} (2z^2 - 1, x z) \\
+\theta & = & \arcsin (yz) \\
+{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x/2)^2 - y^2}
+\end{eqnarray}
+\paragraph{Parabolic}
+The Parabolic projection (``PAR'') is a pseudocylindrical projection, and is defined:
+\begin{eqnarray}
+x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\
+y & = & \pi \sin \frac{\theta}{3} \\
+\end{eqnarray}
+And in reverse:
+\begin{eqnarray}
+\theta & = & 3 \sin^{-1} \rho \\
+\phi   & = & \frac{x}{1 - 4\rho^2} \\
+{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\
+\end{eqnarray}
+\subsection{Offset}
+Coordinate offsets can be either spherical offsets or linear offsets.
+A spherical offset is performed by adding the components of the
+offset, after unit conversion, to the given position.  The resulting
+coordinates must be wrapped to within the allowed range ($-\pi$ to
+$\pi$, 0 to $2\pi$).
+A linear offset is defined to be a linear offset in a tangent
+projection centered on the starting coordinate with $y$ axis aligned
+with the local direction or increasing Declination.  This projection
+is undefined only for the coordinates exactly at the north and south
+poles, in which case the orientation is defined to have the $y$ axis
+parallel to the line of RA = 0.0.  The scale of the projection is 1.0
+(ie, 1 'pixel' is 1 radian) and the given offsets must the scaled
+based on the given offset units.
+Pseudo-code to implement the above for an offset:
+\begin{verbatim}
+psSphere *psSphereSetOffset (psSphere pos, psSphere offset) {
+  psPlane lin;
+  psSphere new;
+  psProjection proj;
+  proj.R = pos->r;
+  proj.D = pos->d;
+  proj.X = 0;
+  proj.Y = 0;
+  proj.type = PS_PROJ_TAN;
+  lin.x = offset.r;
+  lin.y = offset.d;
+  new = psDeproject (&lin, &proj);
+  return (new);
+}
+\end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{The One-to-Many Problem with Mosaic Cameras}
 The \PS{} focal plane consists of several chips, so we will often want
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsubsection{General Astronomy Functions}
+\subsection{General Astronomy Functions}
 \tbd{we will provide a new airmass function}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsubsection{Positions of Major Solar System Objects}
+\subsection{Positions of Major Solar System Objects}
 \tbd{ephemerides code to replace this}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Missing and Todo}
+\tbd{define sunrise, sunset, sun position}
+\tbd{define moonrise, moonset, moon position, moon phase}
+\tbd{define planet functions}
+\tbd{clean up FITS I/O issues}
+\tbd{define Brent's method \& minimization bracketing}
+\pagebreak
 \section{Pan-STARRS Modules}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Missing and Todo}
+\tbd{define sunrise, sunset, sun position}
+\tbd{define moonrise, moonset, moon position, moon phase}
+\tbd{define planet functions}
+\tbd{clean up FITS I/O issues}
+\tbd{define Brent's method \& minimization bracketing}
 \appendix
 \section{Change Log}
 …
 \end{document}
-\section{Modules}
-\subsection{Image Processing Modules}
-\subsubsection{debias}
-\subsubsection{mask}
-\subsubsection{trim}
-\subsubsection{flatten}
-\subsubsection{sky/fringe subtract}
-\subsubsection{warp}
-\subsubsection{stack}
-\subsubsection{difference}
-\subsubsection{kernel convolution}
-\subsubsection{special stack}
-\subsection{Object Detection Modules}
-\subsubsection{find peaks}
-\subsubsection{background }
-\subsubsection{aperture photometry}
-\subsubsection{get shape}
-\subsection{Miscellaneous Modules}
-\section{Analysis Stages}
-\subsection{Phase 1}
-\subsection{Phase 2}
-\subsection{Phase 3}
-\subsection{Phase 4}
-\subsection{Cal 1}
-\subsection{Cal 2}
-\subsection{Cal 3}
-\subsection{Astrom Ref}
-\subsection{Photom Ref}
-\section{Architectual Components}

trunk/doc/pslib/psLibSDRS.tex

-              r3537
+              r3564
 %%% $Id: psLibSDRS.tex,v 1.191 2005-03-29 03:42:21 price Exp $
+%%% $Id: psLibSDRS.tex,v 1.192 2005-03-30 21:14:48 eugene Exp $
 \documentclass[panstarrs,spec]{panstarrs}
 …
 \project{Pan-STARRS Image Processing Pipeline}
 \organization{Institute for Astronomy}
 \version{12}
+\version{13}
 \docnumber{PSDC-430-007}
+% \setcounter{tocdepth}{5} % lowest level to be included in toc
 \setlength{\topsep}{-2pt}
 …
 \href{heasarc.gsfc.nasa.gov/docs/software/fitsio}{\tt heasarc.gsfc.nasa.gov/docs/software/fitsio}
+\item \tbd{SLALIB support is likely to be dropped} Many of the
+astronomy routines will wrap the StarLink Positional Astronomy
+libraries (SLALib):
+\item the StarLink Positional Astronomy libraries (SLALib) are a
+  useful reference:
 \href{star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{\tt
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\pagebreak
 \section{System Utilities}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\pagebreak
 \section{Basic Data Types and Collections}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\pagebreak
 \section{Data manipulation}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Image Pixel Lists}
+Usually an image mask is the best way to carry information about what
+pixels mean what.  However, in the case where the number of pixels in
+which we are interested is limited, it is more efficient to simply
+carry a list of pixels.  An example of this is in the image
+combination code, where we want to perform an operation on a
+relatively small fraction of pixels, and it is inefficient to go
+through an entire mask image checking each pixel.
+\begin{verbatim}
+typedef struct {
+    psVector *x;                        // x coordinate
+    psVector *y;                        // y coordinate
+} psPixels;
+\end{verbatim}
+Of course, the size of each of the vectors should match.  In the event
+that they do not match, any function which detects the problem shall
+generate a warning and use the size of the shorter of the vectors as
+the size.  The order in which the pixels are kept is not considered
+important.
+\begin{verbatim}
+psImage *psPixelsToMask(psImage *out, const psPixels *pixels, const psRegion *region, unsigned int maskVal);
+psPixels *psMaskToPixels(psPixels *out, const psImage *mask, unsigned int maskVal);
+\end{verbatim}
+\code{psPixelsToMask} shall return an image of type U8 with the
+\code{pixels} lying within the specified \code{region} set to the
+\code{maskVal}.  The \code{out} image shall be modified if supplied,
+or allocated and returned if \code{NULL}.  The size of the output
+image shall be \code{region->x1 - region->x0} by \code{region->y1 -
+region->y0}, with \code{out->x0 = region->x0} and \code{out->y0 =
+region->y0}.  In the event that either of \code{pixels} or
+\code{region} are \code{NULL}, the function shall generate an error
+and return \code{NULL}.
+\code{psMaskToPixels} shall return a \code{psPixels} consisting of the
+coordinates in the \code{mask} that match the \code{maskVal}.  The
+\code{out} pixel list shall be modified if supplied, or allocated and
+returned if \code{NULL}.  In hte event that \code{mask} is
+\code{NULL}, the function shall generate an error and return
+\code{NULL}.
+\begin{verbatim}
+psPixels *psPixelsConcatenate(psPixels *out, const psPixels *pixels);
+\end{verbatim}
+\code{psPixelsConcatenate} shall concatenate \code{pixels} onto
+\code{out}.  In the event that \code{out} is \code{NULL}, a new
+\code{psPixels} shall be allocated, and the contents of \code{pixels}
+simply copied in.  If \code{pixels} is \code{NULL}, the function shall
+generate an error and return \code{NULL}.  The function shall take
+care to ensure that there are no duplicate pixels in \code{out} (since
+the order in which the pixels are stored is not important, the values
+may be sorted, allowing the use of a faster algorithm than a linear
+scan).
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Image Regions}
+In many places, we need to refer to a rectangular area.  We define a
+structure to represent a rectangle:
+\begin{verbatim}
+typedef struct {
+  float x0;
+  float x1;
+  float y0;
+  float y1;
+} psRegion;
+psRegion *psRegionAlloc (float x0, float x1, float y0, float y1);
+\end{verbatim}
+\begin{verbatim}
+psRegion *psRegionFromString (char *region);
+\end{verbatim}
+This function converts the IRAF description of a region in the form
+\code{[x0:x1,y0:y1]}, used for header entries such as \code{BIASSEC},
+into the corresponding \code{psRegion} structure.
 \subsection{Vector and Image Arithmetic}
 \label{sec:arithmetic}
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\pagebreak
+\section{Rich Data Structures and I/O}
+\subsection{Metadata}
+\label{sec:metadata}
+\subsubsection{Conceptual Overview}
+Within PSLib, we provide a data structure to carry metadata and
+mechanisms to manipulate the metadata.  Metadata is a general concept
+that requires some discussion.  In any data analysis task, the
+ensemble of all possible data may be divided into two or three
+classes: there is the specific data of interest, there is data which
+is related or critical but not the primary data of interest, and there
+is all of the other data which may or may not be interesting.  For
+example, consider a simple 2D image obtained of a galaxy from a CCD
+camera on a telescope.  If you want to study the galaxy, the specific
+data of interest is the collection of pixels.  There are a variety of
+other pieces of data which are closely related and crucial to
+understanding the data in those pixels, such as the dimensions of the
+image, the coordinate system, the time of the image, the exposure
+time, and so forth.  Other data may be known which may be less
+critical to understanding the image, but which may be interesting or
+desired at a later date.  For example, the observer who took the
+image, the filter manufacturer, the humidity at the telescope, etc.
+Formally, all of the related data which describe the principal data of
+interest are metadata.  Note that which piece is the metadata and
+which is the data may depend on the context.  If you are examining the
+pixels in an image, the coordinate and flux of an object may be part
+of the metadata.  However, if you are analyzing a collection of
+objects extracted from an image, you may consider then pixel data
+simply part of the metadata associated with the list of objects.
+There are various ways to handle metadata vs data within a programming
+environment.  In C, it is convenient to use structures to group
+associated data together.  One possibility is to define the metadata
+as part of the associated data structure.  For example, the image data
+structure would have elements for all possible associated measurement.
+This approach is both cumbersome (because of the large number metadata
+types), impractical (because the full range of necessary metadata is
+difficult to know in advance) and inflexible (because any change in
+the collection of metadata requires addition of new structure elements
+and recompilation).
+An alternative is to place the metadata in a generic container and use
+lookup mechanisms to extract the appropriate metadata when needed.  An
+example of this is would be a text-based FITS header for an image read
+into a flat text buffer.  In this implementation, metadata lookup
+functions could return the current value of, for example, NAXIS1 (the
+number of columns of the image) by scanning through the header buffer.
+This method has the benefits of flexibility and simplicity of
+programming interface, but it has the disadvantage that all metadata
+is accessed though this lookup mechanism.  This may make the code less
+readable and it may slow down the access.
+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,
+\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
+new dimensions.  What should it do to the corresponding metadata?
+Clearly, it should change the corresponding value which defines
+\code{nX, nY}.  However, it is not quite so simple: there may be other
+metadata values which depend on those values.  These must also be
+changed appropriately.  What if the metadata element points to a
+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.  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
+structure have the advantage of being hardwired: every instance of
+that structure will have the same name for the same entry.  This is
+not necessarily the case with a more flexible metadata container.  The
+image exposure time is a notorious example in astronomy.  Different
+observatories use different header keywords (ie, metadata names) for
+the same concept of the exposure time (\code{EXPTIME},
+\code{EXPOSURE}, \code{OPENTIME}, \code{INTTIME}, etc).  Any system
+which operates on these metadata needs to address the issue of
+identifying these names.  This issue seems like an argument for
+hardwiring metadata in the structure, but in fact it does not present
+such a strong case.  If the metadata are hardwired, some function will
+still have to know how to interpret the various names to populate the
+structure.  The concept can still be localized with generic metadata
+containers by including abstract metadata names within the code which
+are tied to the various implementations-specific metadata names.
+\subsubsection{Metadata Representation}
+\begin{figure}
+\psfig{file=Metadata,width=6.5in}
+\caption{Metadata Structures\label{fig:metadata}}
+\end{figure}
+This section addresses the question of how \PS{} metadata should be
+represented in memory, not how it should be represented on disk.
+We define an item of metadata with the following structure:
+\filbreak
+\begin{verbatim}
+typedef struct {
+    int id;                             ///< unique ID for this item
+    char *name;                         ///< Name of item
+    psMetadataType type;                ///< type of this item
+    psElemType ptype;                   ///< primitive data type
+    const union {
+        psS32 S32;                      ///< integer data
+        psF32 F32;                      ///< floating-point data
+        psF64 F64;                      ///< double-precision data
+        void *V;                        ///< other type
+        psList *list;                   ///< psList entry
+        psMetadata *md;                 ///< psMetadata entry
+    } data;                             ///< value of metadata
+    char *comment;                      ///< optional comment ("", not NULL)
+} psMetadataItem;
+\end{verbatim}
+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.  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{V}.  A character string comment
+associated with this metadata item may be stored in the element
+\code{comment}. The \code{type} entry specifies how to interpret the
+type of the data being represented, given by the enumerated type
+\code{psMetadataType}:
+%
+\filbreak
+\begin{verbatim}
+typedef enum {                          ///< type of item.data is:
+    PS_META_PRIMITIVE,                  ///< primitive type: use item.ptype
+    PS_META_LIST,                       ///< psList; use item.data.list (used for non-unique data)
+    PS_META_META,                       ///< psMetadata: use item.data.list
+    PS_META_STR,                        ///< string (item.data.V)
+    PS_META_MATH,                       ///< psScalar, psVector, psImage (item.data.V)
+    PS_META_JPEG,                       ///< JPEG (item.data)
+    PS_META_PNG,                        ///< PNG (item.data)
+    PS_META_ASTROM,                     ///< astrometric coefficients (item.data)
+    PS_META_UNKNOWN,                    ///< other (item.data)
+    PS_META_NTYPE                       ///< Number of types; must be last
+} psMetadataType;
+\end{verbatim}
+If the data is a PSLib primitive data value, the primitive data type
+is given by the value of \code{ptype}.
+A collection of metadata is represented by the \code{psMetadata} structure:
+\begin{verbatim}
+typedef struct {
+    psList *list;                       ///< list of psMetadataItem
+    psHash *table;                      ///< hash table of the same metadata
+} 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
+(each \code{psMetadataItem} appears on both).  The first
+representation employs a doubly-linked list that allows the order of
+the metadata to be preserved (e.g., if FITS headers are read in a
+particular order, they should be written in the same order).  The
+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.  In such a case, the \code{psMetadata.list} structure contains
+the entries in their original sequence with duplicate keys.  The
+\code{psMetadata.hash} entries, which are required to have unique
+keys, would have a single entry with the keyword of the repeated key,
+with the value of \code{psMetadataType} set to \code{PS_META_LIST},
+and the \code{psMetadataItem.data} element pointing to a \code{psList}
+containing the actual entries.  If \code{psMetadataItemAlloc} is
+called with the type set to \code{PS_META_LIST}, such a repeated key
+is created.  If the data value passed to \code{psMetadataItemAlloc}
+(the quantity in ellipsis) is \code{NULL}, then an empty
+\code{psMetadataItem} with the given keyword is created to hold future
+entries of that keyword.
+The \code{psMetadataAdd} routine is required to check that all
+metadata names are unique unless the type is already qualified as
+\code{PS_META_LIST}; in this case the data are added to the
+corresponding \code{psMetadataItem.data} list.
+\subsubsection{Metadata APIs}
+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}-type formating codes.  The
+\code{comment} entry is a fixed string which is used for the comment
+associated with this metadata item.  The metadata data and the
+arguments to the \code{name} formatting codes are passed, in that
+order (metadata pointer first), 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, psMetadataType type, const char *comment, ...);
+psMetadataItem *psMetadataItemAllocV(const char *name, psMetadataType type, const char *comment, va_list list);
+\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);
+\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 return value defines the success
+(\code{true}) or failure of the operation.  The second function,
+\code{psMetadataAdd} takes a pointer or value which is interpreted by
+the function using variadic argument interpretation.  The third
+version is the \code{va_list} version of the second function.  All
+three functions take a parameter, \code{location}, which specifies
+where in the list to place the element, following the conventions for
+the \code{psList}.  The entry \code{mode} for \code{psMetadataAddItem}
+is a bit mask constructed by OR-ing the allowed option flags (eg,
+\code{PS_META_REPLACE}) which specifies minor variations on the
+behavior.  The \code{format} entry, which specifies both the metadata
+type and the optional flags, is constructed by bit-wise OR-ing the
+appropriate \code{psMetadataType} and allowed option option flags.
+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_LIST}).
+%
+\begin{verbatim}
+bool psMetadataAddItem(psMetadata *md, psMetadataItem *item, int location, int mode);
+bool psMetadataAdd(psMetadata *md, int location, const char *name, int format, const char *comment, ...);
+bool psMetadataAddV(psMetadata *md, int location, const char *name, int format, const char *comment,
+                    va_list list);
+\end{verbatim}
+The functions above take option flags which modify the behavior when
+metadata items are added to the metadata list.  These flags must be
+bit-exclusive of those used above for the \code{psMetadataTypes}.  The
+flags have the following meanings:
+\code{PS_META_DEFAULT}: This is the zero bit mask, to allow the
+default behavior for \code{psMetadataAddItem} above.  If this is OR-ed
+with a \code{psMetadataType}, the result is as if no OR-ing took
+place.
+\code{PS_META_REPLACE}: If the given metadata item exists in the
+metadata list, and is not of type \code{PS_META_LIST} or
+\code{PS_META_META} (ie, not a container type), then this entry is
+allowed to replace the existing entry.  If this mode bit is not set, a
+duplicate (non-container-type) entry is an error.
+\begin{verbatim}
+typedef enum {                          ///< option flags for psMetadata functions
+    PS_META_DEFAULT,                    ///< default behavior (0x0000) for use in mode above
+    PS_META_REPLACE,                    ///< allow entry to be replaced
+} psMetadataFlags;
+\end{verbatim}
+An example of code to use these metadata APIs to generate the
+structure seen in Figure~\ref{fig:metadata} is given below.
+\begin{verbatim}
+md = psMetadataAlloc();
+psMetadataAdd(md, PS_LIST_TAIL, "SIMPLE",   PS_META_BOOL, "basic fits",            TRUE);
+psMetadataAdd(md, PS_LIST_TAIL, "BLANK",    PS_META_S32,  "invalid pixel data",    -32768);
+psMetadataAdd(md, PS_LIST_TAIL, "DATE-OBS", PS_META_STR,  "observing date UT", "   2004-6-16");
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_LIST, "head of comment block", NULL);
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "DATA");
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "PARAMS");
+psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_META_F32,  "exposure time (sec)",   1.05);
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "FOO");
+cell = psMetadataAlloc();
+psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_META_STR,  "",                    "CCD00");
+psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR,  "",                    "BSEC-00");
+psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_META_STR,  "",                    "CHIP.00");
+psMetadataAdd(md,   PS_LIST_TAIL, "CELL.00",  PS_META_META, "",                    cell);
+cell = psMetadataAlloc();
+psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_META_STR,  "",                    "CCD01");
+psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR,  "",                    "BSEC-01");
+psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_META_STR,  "",                    "CHIP.01");
+psMetadataAdd(md,   PS_LIST_TAIL, "CELL.01",  PS_META_META, "",                    cell);
+\end{verbatim}
+The following code shows how to use the APIs to replace one of these values:
+\begin{verbatim}
+psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_META_F32 | PS_REPLACE,  "new exposure time (sec)",   2.05);
+\end{verbatim}
+Items may be removed from the metadata by specifying a key or a
+location in the list.  If the value of \code{name} is \code{NULL}, the
+value of \code{location} is used.  If the value of \code{name} is not
+\code{NULL}, then \code{location} must be set to
+\code{PS_LIST_UNKNOWN}.  If the key matches a metadata item, the item
+is removed from the metadata and \code{true} is returned; otherwise,
+\code{false} is returned.  If the key is not unique, then \emph{all}
+items corresponding to the key are removed, and \code{true} is
+returned.
+%
+\begin{verbatim}
+bool psMetadataRemove(psMetadata *md, int location, const char *key);
+\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 *md, const char *key);
+\end{verbatim}
+Several utility functions are provided for simple cases.  These
+functions perform the effort of casting the data to the appropriate
+type.  The numerical functions shall return 0.0 if their key is not
+found.  If the pointer value of \code{status} is not \code{NULL}, it
+is set to reflect the success or failure of the lookup.
+\begin{verbatim}
+void *psMetadataLookupPtr(bool *status, const psMetadata *md, const char *key);
+psS32 psMetadataLookupS32(bool *status, const psMetadata *md, const char *key);
+psF64 psMetadataLookupF64(bool *status, const psMetadata *md, const char *key);
+\end{verbatim}
+Items may be retrieved from the metadata by their entry position.  The
+value of which specifies the desired entry in the fashion of
+\code{psList}.
+\begin{verbatim}
+psMetadataItem *psMetadataGet(const psMetadata *md, int location);
+\end{verbatim}
+The metadata list component may be iterated over by using a
+\code{psListIterator} in a fashion equivalent to the usage for
+\code{psList}.  The iterator may be set to a location in the
+\code{psMetadata} list, and the user may get the previous or next item
+in the list relative to that location.  \code{psMetadataGetNext} has
+the ability to match the key using a POSIX regex, 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 in the metadata list, even those that are
+contained in a \code{PS_META_LIST}.  The value \code{iterator}
+specifies the iterator to be used.  In setting the iterator, the
+position of the iterator is defined by \code{location}, which follows
+the conventions of the \code{psList} iterators.
+\begin{verbatim}
+psListIterator *psMetadataIteratorAlloc(psMetadata *md, int location, bool mutable);
+bool psMetadataIteratorSet(psListIterator *iterator, int location);
+psMetadataItem *psMetadataGetAndIncrement(psListIterator *iterator, const char *regex);
+psMetadataItem *psMetadataGetAndDecrement(psListIterator *iterator, const char *regex);
+\end{verbatim}
+Metadata items may be printed to an open file descriptor based on a
+provided format.  The format string is an sprintf format statement
+with exactly one \% formatting command.  If the metadata item type is
+a numeric type, this formatting command must also be numeric, and type
+conversion performed to the value to match the format type.  If the
+metadata item type is a string, the formatting command must also be
+for a string (\%s type of command).  If the metadata type is any other
+data type, printing is not allowed.
+\begin{verbatim}
+bool psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *md);
+\end{verbatim}
+\subsubsection{Configuration files}
+\label{sec:configspec}
+It will be necessary for the \PS{} system, in order to load
+pre-defined settings, to parse a configuration file into a
+\code{psMetadata} structure.  This shall be performed by the
+function \code{psMetadataParseConfig}, as described below.
+\begin{verbatim}
+psMetadata *psMetadataParseConfig(psMetadata *md, int *nFail, const char *filename, bool overwrite);
+\end{verbatim}
+Given a metadata container, \code{md}, and the name of a configuration
+file, \code{filename}, \code{psMetadataParseConfig} shall parse the
+configuration file, placing the contained key/type/value/comment quads
+into the metadata, and returning a pointer to the metadata structure.
+The number of lines that failed to parse is returned in \code{nFail}.
+Multiple specifications of a key that haven't been declared (see
+below) are overwritten if and only if \code{overwrite} is \code{true}.
+If the metadata container is \code{NULL}, it shall be allocated.
+On error, the function shall return \code{NULL}.
+The configuration file shall consist of plain text with
+key/type/value/comment quads on separate lines.  Blank lines,
+including those consisting solely of whitespace (both spaces and
+tabs), shall be ignored, as shall lines that commence with the comment
+character (a hash mark, \code{#}), either immediately at the start of
+the line, or preceded by whitespace.  The key/type/value/comment quads
+shall all lie on a single line, separated by whitespace.
+The key shall be first, possibly preceded on the line by whitespace
+which should not form part of the key.
+Next, to assist the casting of the value, shall be a string
+identifying the type of the value, which shall correspond to one of
+the simple types supported in \code{psMetadata}:
+\code{STRING,BOOL,S32,F32,F64}; \code{STR} may be used to abbreviate
+\code{STRING}.
+\tbd{May, in the future, require more types, including U8,S16,C64,
+which will also necessitate updating the definition of psMetadata.}
+The value shall follow the type: strings may consist of multiple
+words, and shall have all leading and trailing whitespace removed;
+booleans shall simply be either \code{T} or \code{F}.
+Following the value may be an optional comment, preceded by a comment
+character (a hash mark, \code{#}), which in the case of a string
+value, serves to mark the end of the value, and for other types serves
+to identify the comment to the reader.  Only one comment character may
+be present on any single line (i.e., neither strings nor comments are
+permitted to contain the comment character).  The comment may consist
+of multiple words, and shall have leading and trailing whitespace
+removed.
+One wrinkle is the specification of vectors.  Keys for which the value
+is to be parsed as a vector shall be preceded immediately by a
+``vector symbol'', which we choose to be the ``at'' sign, \code{@}.
+In this case, the type shall be interpreted as the type for the
+vector, which may be any of the signed or unsigned integer or floating
+point types (\code{U8,U16,U32,U64,S8,S16,S32,S32,S64,F32,F64}) but not
+the complex floating point types; and the value shall consist of
+multiple numbers, separated either by a comma or whitespace.  These
+values shall populate a \code{psVector} of the appropriate type in the
+order in which they appear in the configuration file.
+\tbd{May add complex types, likely to be specified with values such as
+.23+4.56i in the future.}
+An additional hurdle is the specification of keys that may be
+non-unique (such as the \code{COMMENT} keyword in a FITS header).
+These keys shall be specified in the configuration file as non-unique
+by specifying the key at the start of the line (possibly preceded by
+whitespace) and specifying the type as a ``multiple symbol'', which we
+choose to be an asterisk, \code{*}.  No other data may be provided on
+this line, though a comment, preceeded by the comment marker, is
+valid.  A warning shall be produced when a key which has not been
+specified to be non-unique is repeated; in this case, the former value
+shall be overwritten if \code{overwrite} is \code{true}, otherwise the
+line shall be ignored and counted as one that could not be parsed.
+If a line does not conform to the rules laid out here, a warning shall
+be generated, it shall be ignored and counted as a line that could not
+be parsed.  The total number of lines that were not able to be parsed
+(including those that were ignored because \code{overwrite} is
+\code{false}, and any other parsing problems, but not including blank
+lines and comment lines) shall be returned by the function in the
+argument \code{nFail}.
+Here are some examples of lines of a valid configuration file:
+\filbreak
+\begin{verbatim}
+Double     F64     1.23456789      # This is a comment
+Float    F32 0.98765 # This is a comment too
+String  STR This is the string that forms the value #comment
+ # This is a comment line and is to be ignored
+boolean     BOOL    T # The value of `boolean' is `true'
+@primes U8  2,3 5 7,11,13 17 #   These are prime numbers
+comment MULTI # The rest of this line is ignored, but `comment' is set to be non-unique
+comment STR This
+comment STR     is
+comment STR       a
+comment STR        non-unique
+comment STR                  key
+Float F64 1.23456 # This generates a warning, and, if `overwrite' is `false', is ignored
+\end{verbatim}
+Of course, a real configuration file should look much nicer to humans
+than the above example, but PSLib must be able to parse such ugly
+files.
+We extend \code{psMetadataParseConfig} to allow a modest tree
+structure by defining a reserved keyword \code{TYPE}.  Any line in the
+config file which starts with the word \code{TYPE} shall be
+interpretted as defining a new valid type.  The defined type name
+follows the word \code{TYPE}, and is in turn followed by an arbitrary
+number of words.  These words are to be interpreted as the names of an
+embedded \code{psMetadata} entry, where the values are given on any
+line which (following the \code{TYPE} definition) employs the new type
+name.  For example, a new type may be defined as:
+\begin{verbatim}
+TYPE      CELL   EXTNAME   BIASSEC  CHIP
+CELL.00   CELL   CCD00     BSEC-00  CHIP.00
+CELL.01   CELL   CCD01     BSEC-01  CHIP.00
+\end{verbatim}
+When \code{psMetadataParseConfig} encounters the \code{TYPE} line, it
+should construct a \code{psMetadata} container and fill it with
+\code{psMetadataItems} having the names \code{EXTNAME, BIASSEC, CHIP},
+with type \code{PS_META_STR}, but data allocated.  When it next
+encounters an entry of type \code{CELL}, it should then use the given
+name (e.g., \code{CELL.00}) for the \code{psMetadataItem}, and copy
+the \code{psMetadata} data onto the \code{psMetadataItem.data.md}
+entry, filling in the values from the rest of the line (\code{CCD00,
+BSEC-00, CHIP.00}).  This hierarchical structure is illustrated in
+Figure~\ref{fig:metadata}.
+We further extend \code{psMetadataParseConfig} to allow the definition
+of a \code{psMetadata} entry using a sequence of successive lines to
+define the values of the \code{psMetadataItem} entries.  The initial
+line defines the new \code{psMetadata} entry and its name.  The
+following lines have the same format as the other metadata config file
+entries.  The sequence is terminated with a line with a single word
+\code{END}.  For example, a metadata entry may be defined as:
+\begin{verbatim}
+CELL      METADATA
+ EXTNAME   STR   CCD00
+ BIASSEC   STR   BSEC-00
+ CHIP      STR   CHIP.00
+ NCELL     S32   24
+END
+\end{verbatim}
+A series of test inputs is contained in
+\S\ref{sec:configtest}.
+\subsection{XML Functions}
+Within Pan-STARRS, we will use XML documents as a transport mechanism
+to carry data between programs and between IPP and other subsystems.
+Configuration information may be stored as well as XML documents, in
+addition to the text format discussed in the discussion on Metadata.
+XML is an extremely variable document format, and it is not currently
+the intention of PSLib to provide a complete PSLib version of XML
+operations.  Rather, a limited number of operations are defined to
+convert specific data structures to an appropriate XML document.  To
+maximize the simplicity of the XML APIs, we will use the convention
+that a single XML document to be parsed by PSLib shall contain only a
+single data structure.  Each of the XML APIs therefore take as input a
+reference to a complete XML document and return a PSLib data
+structure, or take a PSLib data structure and return a complete XML
+document.
+We start by defining a PSLib wrapper type which is a pointer to an XML
+document in memory.  We wrap the \code{libxml2} version of an XML
+document pointer for now:
+\begin{verbatim}
+typedef xmlDocPtr psXMLDoc;
+void psXMLDocFree(psXMLDoc *doc);
+\end{verbatim}
+The next pair of functions convert a \code{psMetadata} data structure
+to a complete \code{psXMLDoc} (in memory) and vice versa:
+\begin{verbatim}
+psXMLDoc *psMetadataToXMLDoc(const psMetadata *metadata);
+psMetadata *psXMLDocToMetadata(const psXMLDoc *doc);
+\end{verbatim}
+The next pair of functions loads the data in a named file into a
+complete \code{psXMLDoc} (in memory) and write out the \code{psXMLDoc}
+to a named file:
+\begin{verbatim}
+psXMLDoc *psXMLParseFile(const char *filename);
+int psXMLDocToFile(const psXMLDoc *doc, const char *filename);
+\end{verbatim}
+The next pair of functions accepts a block of memory and parses it
+into a complete \code{psXMLDoc} (also in memory), and vice versa:
+\begin{verbatim}
+psXMLDoc *psXMLParseMemory(const char *buffer, const int size);
+int psXMLDocToMemory(const psXMLDoc *doc, char *buffer);
+\end{verbatim}
+The next pair of functions read from and write to a file descriptor.
+The first converts the imcoming data to a complete \code{psXMLDoc}
+(also in memory), the second writes the \code{psXMLDoc} to the file
+descriptor:
+\begin{verbatim}
+psXMLDoc *psXMLParseFD(int fd);
+int psXMLDocToFD(const psXMLDoc *doc, int fd);
+\end{verbatim}
+\subsection{Database Functions}
+Many of the applications that PSLib will be used for will require
+access to a simple relational database.  PSLib includes generic
+database-independent interface mechanisms as part of its API set.  The
+most important aspect of PSLib's database support is to abstract as
+much database specific complexity as is feasible.  As almost all RDBMS
+provide at least a simple transactional model, commit and rollback
+support should be provided.
+Currently, only support for MySQL 4.1.x is required but other backends
+may be added as options in the future.  As a particular example which
+has implications for the database interaction model, support for
+SQLite may be required in the future.  Currently, the choice of
+backend database interface may be made as a compile option.  Details
+of the specified APIs in the discussion below refer to the relevant
+MySQL functions.
+Database errors must be trapped and placed onto the psError stack.
+The complete error message should be retrieved with the database's
+error function.
+\subsubsection{Managing the Database Connection}
+We specify a database handle which carries the information about the
+database connection:
+\begin{verbatim}
+    typedef struct {
+        MYSQL *mysql;
+    } psDB;
+\end{verbatim}
+The following collection of functions provides basic database functionality:
+\begin{verbatim}
+    // wraps mysql_init() & mysql_real_connect()
+    psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname);
+    // wraps mysql_close()
+    void psDBCleanup(psDB *dbh);
+    // wraps mysql_create_db()
+    bool psDBCreate(psDB *dbh, const char *dbname);
+    // wraps mysql_select_db()
+    bool psDBChange(psDB *dbh, const char *dbname);
+    // wraps mysql_drop_db()
+    bool psDBDrop(psDB *dbh, const char *dbname);
+\end{verbatim}
+For MySQL support, \code{psDBInit()} wraps \code{mysql_init()} and
+\code{mysql_real_connect()} in order to initialize a psDB structure and
+establish a database connection.  A null pointer should be returned on
+failure.
+When implementing support for SQLite, or other DB which is purely
+file-based, the \code{host}, \code{user}, and \code{passwd} arguments
+would be ignored while \code{dbname} would specify the path to the
+SQLite db file.
+\subsubsection{Interacting with Database Tables}
+The functions in this section perform high level interactions with the
+database tables.  All of them should behave ``atomically'' with
+respect to the state of the database.  Specifically, all interactions
+with the database should be done as a part of a transaction that is
+rolled-back on failure and committed only after all queries used by
+the API have been run.  In general, this API set attempts to treat a
+database table as a 2D matrix where columns can be represented by a
+\code{psVector} and rows as a \code{psMetadata} type.  A
+\code{psMetadata} collection is also used to define the columns of a
+table and as part of the query restrictions.
+\begin{verbatim}
+    bool psDBCreateTable(psDB *dbh, const char *tableName, psMetadata *md);
+\end{verbatim}
+This function generates and executes the SQL needed to create a table
+named \code{tableName}, with the column names and datatypes as
+described in \code{md}.  Each data item in the \code{psMetadata}
+collection represents a single table field.  The name of the field is
+given by the name of the \code{psMetadataItem} and the data type is
+give by the \code{psMetadataItem.type} and \code{psMetadataItem.ptype}
+entries.  A lookup table should be used to convert from PSLib types
+into MySQL compatible SQL data types.  For example, a
+\code{PS_META_STR} would map to an SQL99 varchar.  If the value of
+\code{type} is \code{PS_META_STR} then the \code{psMetadataItem.data}
+element is set to a string with the length for the field written as a
+text string.  The value of the \code{psMetadataItem.data} element is
+unused for the \code{PS_META_PRIMITIVE} types.  Other metadata types
+beyond \code{PS_META_STR} and \code{PS_META_PRIMITIVE} are not allowed
+in a table definition metadata collection.
+Database indexes can be specified setting the \code{comment} field to
+``\code{Primary Key}'' or ``\code{Key}''.  Comment are otherwise
+ignored.
+\begin{verbatim}
+    bool psDBDropTable(psDB *dbh, const char *tableName);
+\end{verbatim}
+This function deletes the specified table.
+\begin{verbatim}
+    psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, const psU64 limit);
+    psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col, psElemType pType, const psU64 limit);
+\end{verbatim}
+These functions generates and executes the SQL needed to select an entire
+column from a table or up to \code{limit} rows from it.  If \code{limit} is 0,
+the entire range is returned.  The database response is processed and a
+\code{psArray} of strings is returned.  The Num version of the function returns
+the data in a \code{psVector}, data cast to \code{pType}.  It returns an error
+(NULL) if the requested field is not a numerical type.
+\begin{verbatim}
+    psArray *psDBSelectRows(psDB *dbh, const char *tableName, psMetadata *where, const psU64 limit);
+\end{verbatim}
+This function returns rows from the specified table which match
+the restrictions given by \code{where}.  The restrictions are
+specified as field / value pairs.  The \code{psMetadata} collection
+where must consist of valid database fields, though the database query
+checking functions may be used to validate the fields as part of the
+query.  If \code{where} is \code{NULL}, then there are no restrictions
+on the rows selected.  The selected rows are returned as a
+\code{psArray} of \code{psMetadata} values, one per row.
+\begin{verbatim}
+    bool psDBInsertOneRow(psDB *dbh, const char *tableName, psMetadata *row);
+\end{verbatim}
+Insert the data from \code{row} into \code{tableName}.  It should be noted in
+the API reference that if fields are specified in \code{row} that do not exist
+in \code{tablename}, the insert will fail.
+\begin{verbatim}
+    bool psDBInsertRows(psDB *dbh, const char *tableName, psArray *rowSet);
+\end{verbatim}
+Similar to \code{psDBInsertOneRow()}, this function inserts many rows at once
+and is atomic for the complete set of rows.
+\begin{verbatim}
+    psArray *psDBDumpRows(psDB *dbh, const char *tableName);
+\end{verbatim}
+Fetch all rows as an psArray of psMetadata.
+\begin{verbatim}
+    psMetadata *psDBDumpCols(psDB *dbh, const char *tableName);
+\end{verbatim}
+Fetch all columns, as either a psVector or a psArray depnding on whether or not
+the column is numeric, and return them in a psMetadata structure where
+psMetadataItem.name contains the column's name.
+\begin{verbatim}
+    psS64 psDBUpdateRows(psDB *dbh, const char *tableName, psMetadata *where, psMetadata *values);
+\end{verbatim}
+Update the columns contained in \code{values} in the row(s) that have a field
+with the value indicated by \code{where} (note that this is only allows very
+limited use of SQL99's ``where'' semantics).  The number of rows modified is
+returned.  A negative value is return to indicate an error. If there are
+multiple psMetadataItems in \code{where} then each item should be considered as
+an additional constraint.  e.g.  ``where foo = x and where bar = y''
+\begin{verbatim}
+    psS64 psDBDeleteRows(psDB *dbh, const char *tableName, psMetadata *where);
+\end{verbatim}
+Delete the rows that are matched by \code{where} using the same semantics for
+\code{where} as in psDBUpdateRow().  A negative value is returned to indicate an
+error.
+\subsection{FITS I/O Functions}
+We need a variety of I/O functions between the disk and certain of our
+PSLib data structures.  We need the ability to access FITS headers,
+images and tables (both ASCII and Binary).  We define here the FITS
+I/O functions, all of which are currently specified as wrappers to
+functions within CFITSIO.  CFITSIO provides a wide range of utilities
+which we do not feel are particularly appropriate as part of a generic
+I/O library, such as assumptions about names which change the data
+interpretation, etc.  We are defining our calls to avoid the hidden
+'features'.  The CFITSIO functions which are wrapped should in general
+be the most basic versions.
+\begin{verbatim}
+typedef struct {
+    fitsfile fd;
+} psFits;
+\end{verbatim}
+We begin by defining a datatype to wrap the CFITSIO \code{fitsfile}
+structure.  This is necessary to allow repeated access to the data in
+a file without multiple open commands (which are expensive).
+\subsubsection{FITS File Manipulations}
+\begin{verbatim}
+psFits *psFitsAlloc(const char *filename);
+\end{verbatim}
+Opens a FITS file at positions the pointer to the PHU.
+\begin{verbatim}
+bool psFitsMoveExtName(psFits *fits, const char *extname);
+\end{verbatim}
+Positions the pointer to the beginning of the specified
+\code{extname}.  If the \code{extname} does not exist, the function
+shall fail.
+\begin{verbatim}
+bool psFitsMoveExtNum(psFits* fits, int extnum, bool relative);
+\end{verbatim}
+Moves the pointer to the beginning of the specified HDU number.  If
+\code{relative} is TRUE, \code{extnum} represents the number of HDUs
+relative to the current HDU.  The PHU is entry number 0, while the
+extended data segments start at number 1.
+\begin{verbatim}
+int psFitsGetExtNum(psFits* fits);
+\end{verbatim}
+Returns the current HDU number (i.e., file position).
+\begin{verbatim}
+int psFitsGetSize(psFits* fits);
+\end{verbatim}
+Returns the number of HDUs in the file.
+\begin{verbatim}
+psFitsType psFitsGetExtType(psFits* fits);
+\end{verbatim}
+Gets the current HDU's type (table or image).
+\subsubsection{FITS Header I/O Functions}
+\begin{verbatim}
+psMetadata *psFitsReadHeader(psMetadata *out, const psFits *fits);
+\end{verbatim}
+Read header data into a \code{psMetadata} structure.  The data is read
+from the current HDU pointed at by the \code{psFits *fits} entry.  If
+\code{out} is \code{NULL}, a new psMetadata is created.
+\begin{verbatim}
+psMetadata *psFitsReadHeaderSet (psFits *fits);
+\end{verbatim}
+Load a complete set of headers from the \code{psFits} file pointer.
+This function loads the headers from all extensions into a
+\code{psMetadata} collection, each entry of which is a pointer to a
+\code{psMetadata} structure containing the header data.  The metadata
+entry names are the \code{EXTNAME} values for each header (with the
+value of \code{PHU} for the primary header unit).  At the start of the
+operation, the file pointer is rewound to the beginning of the file.
+At the end, it is positioned where it started when the function was
+called.
+\begin{verbatim}
+bool psFitsWriteHeader(psMetadata *output, const psFits *fits);
+\end{verbatim}
+Write metadata into the header of a FITS image file.  The header is
+written at the current HDU.
+\subsubsection{FITS Image I/O Functions}
+\begin{verbatim}
+psImage *psFitsReadImage(psImage *output, psFits *fits, psRegion region, int z);
+\end{verbatim}
+Read an image or subimage from the \code{psFits} file pointer.  This
+function is a wrapper to the CFITSIO library function.  The input
+parameters allow a full image or a subimage to be read.  The region to
+be read is specified by \code{region}.  A negative value for either of
+\code{region.x1} or \code{region.y1} specifies the size of the region
+to be read counting down from the end of the array.
+If the native image is a cube, the value of z specifies the requested
+slice of the image.  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, or if the image on disk is zero- or
+one-dimensional.  This function need only read images of the native
+FITS image types (\code{psU8}, \code{psS16}, \code{psS32},
+\code{psF32}, \code{psF64}).  The user is expected to convert the data
+type as needed with \code{psImageCopy}.
+\begin{verbatim}
+bool psFitsUpdateImage(psFits *fits, const psImage *input, psRegion region, int z);
+\end{verbatim}
+\tbd{we have discussed this as the alternate name}
+Write an image section to the open \code{psFits} file pointer.  This
+operation may write a portion of an image over the existing bytes of
+an existing image.  Care must be taken to interpret \code{region},
+which specified the output pixels to be written / over-written.  If
+the combination of \code{region} and the size of \code{psImage *input}
+implies writing pixels outside the existing data area of the image,
+the function shall return an error (ie, if \code{region.x0 + image.nx
+>= NAXIS1}, \code{region.y0 + image.ny >= NAXIS2}, or \code{z >=
+NAXIS3}).  This function will only write images of the native FITS
+image types (\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32},
+\code{psF64}).  The user is expected to convert the data type as
+needed with \code{psImageCopy}.  The return value must be 0 for a
+successful operation and 1 for an error.
+\begin{verbatim}
+bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input, int depth);
+\end{verbatim}
+Create a new image based on the dimensions specified for the image and
+the requested depth.  The header and image data segments are written
+in the file at the current position of the \code{psFits} pointer.
+This function will only write images of the native FITS image types
+(\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, \code{psF64}).
+The user is expected to convert the data type as needed with
+\code{psImageCopy}.  The return value must be 0 for a successful
+operation and 1 for an error.
+\subsubsection{FITS Table I/O Functions}
+\begin{verbatim}
+psMetadata *psFitsReadTableRow (psFits *fits, int row);
+\end{verbatim}
+This function reads a single row of the table in the extension pointed
+at by the \code{psFits} file pointer.  The row number to be read is
+given by \code{row}.  The result is returned as a \code{psMetadata}
+collection with elements of the apporpriate types and keys
+corresponding to the table column names.  The function must apply the
+needed byte-swapping on the data in the row based on the description
+of the table data in the table header.  \tbr{we may need to be more
+flexible here: if we call this function repeatedly, it would be more
+efficient to pass the corresponding header or keep it somewhere (and
+the file pointer location, for that matter).}
+\begin{verbatim}
+void *psFitsReadTableRowRaw (int *nBytes, psFits *fits, int row);
+\end{verbatim}
+This function reads a single row of the table in the extension pointed
+at by the \code{psFits} file pointer.  The row number to be read is
+given by \code{row}.  The result is returned as collection of
+\code{nBytes} bytes allocated by the function.  The function must
+apply the needed byte-swapping on the data in the row based on the
+description of the table data in the table header.  \tbr{we may need
+to be more flexible here: if we call this function repeatedly, it
+would be more efficient to pass the corresponding header or keep it
+somewhere (and the file pointer location, for that matter).}
+\begin{verbatim}
+psArray *psFitsReadTableColumn (psFits *fits, char *colname);
+\end{verbatim}
+This function reads a single column of the table in the extension
+pointed at by the \code{psFits} file pointer.  The column is specified
+by the FITS table column key given by \code{row}.  The result is
+returned as a \code{psArray}, with the data from one row of the table
+column per array element.
+\begin{verbatim}
+psVector *psFitsReadTableColumnNum (psFits *fits, char *colname);
+\end{verbatim}
+This function reads a single column of the table in the extension
+pointed at by the \code{psFits} file pointer.  The column is specified
+by the FITS table column key given by \code{row} and must be of a
+numeric data type.  The result is returned as a \code{psVector} of the
+appropriate data type, with the data from one row of the table column
+per array element.
+\begin{verbatim}
+psArray *psFitsReadTableRaw (int *nBytes, psFits *fits);
+\end{verbatim}
+This function reads the entire data block from a table into the a
+\code{psArray}, with one element of the array per row.  The number of
+bytes per row is returned in \code{nBytes}.  The function must apply
+the needed byte-swapping on the data in each row based on the
+description of the table data in the table header.
+\begin{verbatim}
+psArray *psFitsReadTable (psFits *fits);
+\end{verbatim}
+This function reads the entire data block from a table into the a
+\code{psArray}, with one element of the array per row.  Each row is
+stored as a \code{psMetadata} collection as described above for
+\code{psFitsReadTableRow}.
+\begin{verbatim}
+bool psFitsWriteTable(psFits* fits, psMetadata *header, psArray* table);
+\end{verbatim}
+Accepts a \code{psArray} of \code{psMetadata} and writes it to the
+current HDU.  If the current HDU is not a table type, this will fail
+and return FALSE.
+\begin{verbatim}
+bool psFitsUpdateTable(psFits* fits, psMetadata *header, psMetadata* data, int row);
+\end{verbatim}
+Writes the \code{psMetadata} data to a FITS table at the specified row
+in the current HDU.  If the current HDU is not a table type, this will
+fail and return FALSE.
+\pagebreak
 \section{Astronomy-Specific Functions}
 …
 \begin{itemize}
 \item Dates and times
-\item Metadata
 \item Detector and sky positions
 \item Astronomy Image
+\item Astronomical Image Containers
 \item Astrometry
 \item Photometry
 …
 typedef enum {
     PS_IERS_A,                      ///< IERS Bulliten A
     PS_IERS_B,                      ///< IERS Bulliten B
 } psTimeBulliten;
+    PS_IERS_A,                      ///< IERS Bulletin A
+    PS_IERS_B,                      ///< IERS Bulletin B
+} psTimeBulletin;
 typedef struct {
     psS64            sec;           ///< seconds, negative values represent dates before 1970
     psU32            nsec;          ///< nanseconds
+    psU32            nsec;          ///< nanoseconds
     bool             leapsecond;    ///< if time falls on a UTC leapsecond
     psTimeType       type;          ///< type of time
 …
 This function may be used to convert between the various \code{psTimeType} time
 representations.  The \code{time} is modified and returned.  Conversion between
 all of the \emph{SI} length second systems should be implimented as first
+all of the \emph{SI} length second systems should be implimented by first
 converting to TAI and then to the destination system.  UT1 is a special case
 for conversion as it uses variable length secounds.  Conversation to UT1, via
 TAI, is allowed but conversion \emph{from} UT1 is currently forbidden.
+To convert to Local Mean Sidereal Time, it is necessary to provide the local
+longitude (specified in radians, positive East of Greenwich) as well:
+The following function converts to Local Mean Sidereal Time.  It is
+necessary to provide the local longitude (specified in radians,
+positive East of Greenwich) as well:
 \begin{verbatim}
 …
 \end{verbatim}
+The function may accept any of the \code{psTimeType} types with \emph{SI}
+length seconds.  The \code{time} is modified and returned.  Note that this
+function must supply the value $UT1-UTC$, which is available externally (see
+\code{psTimeGetUT1Delta()}) and should use the UT1 values interpolated from
+IERS bulliten B.  The following utility function encapsulates the PSLib
+mechanism to extract the value of $UT1-UTC$:
+\begin{verbatim}
+double psTimeGetUT1Delta(const psTime *time, psTimeBulliten bulliten);
+psSphere *psTimeGetPoleCoords(const psTime *time, psTimeBulliten bulliten);
+\end{verbatim}
+Leap seconds are added to UTC in order to keep it within $0.9s$ of UT1 (which
+is defined relative to the Earth's rotation, and hence is useful for
+astronomical purposes).
+The function may accept any of the \code{psTimeType} types with
+\emph{SI} length seconds.  The \code{time} is modified and returned.
+Note that this function requires the value $UT1-UTC$ (see
+\code{psTimeGetUT1Delta()}) and should use the UT1 values interpolated
+from IERS bulletin B.
+The following utility function encapsulates the PSLib mechanism to
+extract the value of $UT1-UTC$ from the IERS Time Tables:
+\begin{verbatim}
+double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin);
+\end{verbatim}
+The following function provides tidal corrections to UT1-UTC, using
+the Ray model of Simon et al (REF).
+\begin{verbatim}
+psTime *psTime_TideUT1Corr(const psTime *time);
+\end{verbatim}
+Leap seconds are added to UTC in order to keep it within $0.9s$ of UT1
+(which is defined relative to the Earth's rotation, and hence is
+useful for astronomical purposes).  The following function calculates
+the absolute number of leap seconds different between two times.
 \begin{verbatim}
 …
 \end{verbatim}
+This function calculates the absolute number of leap seconds different between
+two times.
+The following function accepts \code{PS_TIME_UTC} objects and
+determines if the time is potentaly a leapsecond, returning
+\code{TRUE} if so.
 \begin{verbatim}
 bool psTimeIsLeapSecond(const psTime *utc);
 \end{verbatim}
-This function only accepts \code{PS_TIME_UTC} objects and determines if the
-time is potentaly a leapsecond.
 \subsubsection{External Date and Time Formats}
 …
 \end{verbatim}
 The \code{psTimeToISO()} function will convertion \code{PS_TIME_UTC} objects
 with the \code{leapsecond} flag set to represent the number of seconds as
 ``:60''.
+The \code{psTimeToISO()} function converts \code{PS_TIME_UTC} objects
+with the \code{leapsecond} flag set to represent the number of seconds
+as ``:60''.
 \subsubsection{Date and Time Parsing}
 …
 \code{psS64} (the type of \code{psTime.sec}).
 Note that in both these functions, when handling UTC, that the difference
+between two times is not inclusive of leap seconds.  For example, if we add 30
+seconds to 1998-12-31T23:59:45Z, the result is 1999-01-01T00:00:14Z, since a
+leap second was introduced at 1999-01-01T00:00:00Z.
+Time math may only be done on the of \code{psTime} objects of the same type,
+and in the case of UT1, the functions shall internally convert the
 \code{psTime} inputs to TAI before performing the math; this is in order that
+leap seconds are accounted for.
+The type of the time returned by \code{psTimeMath} shall be the same as that of
+the input \code{time}.
+Note that in both these functions the difference between two times is
+the elapsed number of seconds, inclusing leap seconds.  For example,
+if we add 30 seconds to 1998-12-31T23:59:45Z, the result is
+-01-01T00:00:14Z, since a leap second was introduced at
+-01-01T00:00:00Z.
+Time math may only be done on \code{psTime} objects of the same type,
+and except in the case of UT1, the functions shall internally convert
+the \code{psTime} inputs to TAI before performing the math; this
+ensures that leap seconds are correctly counted.
+The type of the time returned by \code{psTimeMath} shall be the same
+as that of the input \code{time}.
 \subsubsection{Time Tables}
 …
 \begin{itemize}
 \item IERS Bulliten A \& B (1 year ago $\rightarrow$ now + $\sim$3 months)
+\item IERS Bulletin A \& B (1 year ago $\rightarrow$ now + $\sim$3 months)
 \begin{itemize}
 \item \code{ftp://maia.usno.navy.mil/ser7/finals.daily}
 \end{itemize}
 \item IERS Bulliten A \& B (1973-1-2 $\rightarrow$ now + $\sim$1 year)
+\item IERS Bulletin A \& B (1973-1-2 $\rightarrow$ now + $\sim$1 year)
 \begin{itemize}
 \item \code{ftp://maia.usno.navy.mil/ser7/finals.all}
 …
 format of these files shall be simple, for speed in reading.  Each line shall
 contain the date in MJD and the following values from both the A \& B IERS
 bullitens: $x_p$ (in arcseconds), $y_p$ (in arcseconds) and $\Delta$UT (in
+bulletins: $x_p$ (in arcseconds), $y_p$ (in arcseconds) and $\Delta$UT (in
 seconds).  This format must be readable by \code{psLookupTableRead}.  For
 example:
 …
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Regions}
+In many places, we need to refer to a rectangular area.  We define a
+structure to represent a rectangle:
+\begin{verbatim}
+typedef struct {
+  float x0;
+  float x1;
+  float y0;
+  float y1;
+} psRegion;
+psRegion *psRegionAlloc (float x0, float x1, float y0, float y1);
+\end{verbatim}
+\begin{verbatim}
+psRegion *psRegionFromString (char *region);
+\end{verbatim}
+This function converts the IRAF description of a region in the form
+\code{[x0:x1,y0:y1]}, used for header entries such as \code{BIASSEC},
+into the corresponding \code{psRegion} structure.
+\subsection{Metadata}
+\label{sec:metadata}
+\subsubsection{Conceptual Overview}
+Within PSLib, we provide a data structure to carry metadata and
+mechanisms to manipulate the metadata.  Metadata is a general concept
+that requires some discussion.  In any data analysis task, the
+ensemble of all possible data may be divided into two or three
+classes: there is the specific data of interest, there is data which
+is related or critical but not the primary data of interest, and there
+is all of the other data which may or may not be interesting.  For
+example, consider a simple 2D image obtained of a galaxy from a CCD
+camera on a telescope.  If you want to study the galaxy, the specific
+data of interest is the collection of pixels.  There are a variety of
+other pieces of data which are closely related and crucial to
+understanding the data in those pixels, such as the dimensions of the
+image, the coordinate system, the time of the image, the exposure
+time, and so forth.  Other data may be known which may be less
+critical to understanding the image, but which may be interesting or
+desired at a later date.  For example, the observer who took the
+image, the filter manufacturer, the humidity at the telescope, etc.
+Formally, all of the related data which describe the principal data of
+interest are metadata.  Note that which piece is the metadata and
+which is the data may depend on the context.  If you are examining the
+pixels in an image, the coordinate and flux of an object may be part
+of the metadata.  However, if you are analyzing a collection of
+objects extracted from an image, you may consider then pixel data
+simply part of the metadata associated with the list of objects.
+There are various ways to handle metadata vs data within a programming
+environment.  In C, it is convenient to use structures to group
+associated data together.  One possibility is to define the metadata
+as part of the associated data structure.  For example, the image data
+structure would have elements for all possible associated measurement.
+This approach is both cumbersome (because of the large number metadata
+types), impractical (because the full range of necessary metadata is
+difficult to know in advance) and inflexible (because any change in
+the collection of metadata requires addition of new structure elements
+and recompilation).
+An alternative is to place the metadata in a generic container and use
+lookup mechanisms to extract the appropriate metadata when needed.  An
+example of this is would be a text-based FITS header for an image read
+into a flat text buffer.  In this implementation, metadata lookup
+functions could return the current value of, for example, NAXIS1 (the
+number of columns of the image) by scanning through the header buffer.
+This method has the benefits of flexibility and simplicity of
+programming interface, but it has the disadvantage that all metadata
+is accessed though this lookup mechanism.  This may make the code less
+readable and it may slow down the access.
+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,
+\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
+new dimensions.  What should it do to the corresponding metadata?
+Clearly, it should change the corresponding value which defines
+\code{nX, nY}.  However, it is not quite so simple: there may be other
+metadata values which depend on those values.  These must also be
+changed appropriately.  What if the metadata element points to a
+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.  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
+structure have the advantage of being hardwired: every instance of
+that structure will have the same name for the same entry.  This is
+not necessarily the case with a more flexible metadata container.  The
+image exposure time is a notorious example in astronomy.  Different
+observatories use different header keywords (ie, metadata names) for
+the same concept of the exposure time (\code{EXPTIME},
+\code{EXPOSURE}, \code{OPENTIME}, \code{INTTIME}, etc).  Any system
+which operates on these metadata needs to address the issue of
+identifying these names.  This issue seems like an argument for
+hardwiring metadata in the structure, but in fact it does not present
+such a strong case.  If the metadata are hardwired, some function will
+still have to know how to interpret the various names to populate the
+structure.  The concept can still be localized with generic metadata
+containers by including abstract metadata names within the code which
+are tied to the various implementations-specific metadata names.
+\subsubsection{Metadata Representation}
+\begin{figure}
+\psfig{file=Metadata,width=6.5in}
+\caption{Metadata Structures\label{fig:metadata}}
+\end{figure}
+This section addresses the question of how \PS{} metadata should be
+represented in memory, not how it should be represented on disk.
+We define an item of metadata with the following structure:
+\filbreak
+\begin{verbatim}
+typedef struct {
+    int id;                             ///< unique ID for this item
+    char *name;                         ///< Name of item
+    psMetadataType type;                ///< type of this item
+    psElemType ptype;                   ///< primitive data type
+    const union {
+        psS32 S32;                      ///< integer data
+        psF32 F32;                      ///< floating-point data
+        psF64 F64;                      ///< double-precision data
+        void *V;                        ///< other type
+        psList *list;                   ///< psList entry
+        psMetadata *md;                 ///< psMetadata entry
+    } data;                             ///< value of metadata
+    char *comment;                      ///< optional comment ("", not NULL)
+} psMetadataItem;
+\end{verbatim}
+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.  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{V}.  A character string comment
+associated with this metadata item may be stored in the element
+\code{comment}. The \code{type} entry specifies how to interpret the
+type of the data being represented, given by the enumerated type
+\code{psMetadataType}:
+%
+\filbreak
+\begin{verbatim}
+typedef enum {                          ///< type of item.data is:
+    PS_META_PRIMITIVE,                  ///< primitive type: use item.ptype
+    PS_META_LIST,                       ///< psList; use item.data.list (used for non-unique data)
+    PS_META_META,                       ///< psMetadata: use item.data.list
+    PS_META_STR,                        ///< string (item.data.V)
+    PS_META_MATH,                       ///< psScalar, psVector, psImage (item.data.V)
+    PS_META_JPEG,                       ///< JPEG (item.data)
+    PS_META_PNG,                        ///< PNG (item.data)
+    PS_META_ASTROM,                     ///< astrometric coefficients (item.data)
+    PS_META_UNKNOWN,                    ///< other (item.data)
+    PS_META_NTYPE                       ///< Number of types; must be last
+} psMetadataType;
+\end{verbatim}
+If the data is a PSLib primitive data value, the primitive data type
+is given by the value of \code{ptype}.
+A collection of metadata is represented by the \code{psMetadata} structure:
+\begin{verbatim}
+typedef struct {
+    psList *list;                       ///< list of psMetadataItem
+    psHash *table;                      ///< hash table of the same metadata
+} 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
+(each \code{psMetadataItem} appears on both).  The first
+representation employs a doubly-linked list that allows the order of
+the metadata to be preserved (e.g., if FITS headers are read in a
+particular order, they should be written in the same order).  The
+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.  In such a case, the \code{psMetadata.list} structure contains
+the entries in their original sequence with duplicate keys.  The
+\code{psMetadata.hash} entries, which are required to have unique
+keys, would have a single entry with the keyword of the repeated key,
+with the value of \code{psMetadataType} set to \code{PS_META_LIST},
+and the \code{psMetadataItem.data} element pointing to a \code{psList}
+containing the actual entries.  If \code{psMetadataItemAlloc} is
+called with the type set to \code{PS_META_LIST}, such a repeated key
+is created.  If the data value passed to \code{psMetadataItemAlloc}
+(the quantity in ellipsis) is \code{NULL}, then an empty
+\code{psMetadataItem} with the given keyword is created to hold future
+entries of that keyword.
+The \code{psMetadataAdd} routine is required to check that all
+metadata names are unique unless the type is already qualified as
+\code{PS_META_LIST}; in this case the data are added to the
+corresponding \code{psMetadataItem.data} list.
+\subsubsection{Metadata APIs}
+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}-type formating codes.  The
+\code{comment} entry is a fixed string which is used for the comment
+associated with this metadata item.  The metadata data and the
+arguments to the \code{name} formatting codes are passed, in that
+order (metadata pointer first), 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, psMetadataType type, const char *comment, ...);
+psMetadataItem *psMetadataItemAllocV(const char *name, psMetadataType type, const char *comment, va_list list);
+\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);
+\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 return value defines the success
+(\code{true}) or failure of the operation.  The second function,
+\code{psMetadataAdd} takes a pointer or value which is interpreted by
+the function using variadic argument interpretation.  The third
+version is the \code{va_list} version of the second function.  All
+three functions take a parameter, \code{location}, which specifies
+where in the list to place the element, following the conventions for
+the \code{psList}.  The entry \code{mode} for \code{psMetadataAddItem}
+is a bit mask constructed by OR-ing the allowed option flags (eg,
+\code{PS_META_REPLACE}) which specifies minor variations on the
+behavior.  The \code{format} entry, which specifies both the metadata
+type and the optional flags, is constructed by bit-wise OR-ing the
+appropriate \code{psMetadataType} and allowed option option flags.
+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_LIST}).
+%
+\begin{verbatim}
+bool psMetadataAddItem(psMetadata *md, psMetadataItem *item, int location, int mode);
+bool psMetadataAdd(psMetadata *md, int location, const char *name, int format, const char *comment, ...);
+bool psMetadataAddV(psMetadata *md, int location, const char *name, int format, const char *comment,
+                    va_list list);
+\end{verbatim}
+The functions above take option flags which modify the behavior when
+metadata items are added to the metadata list.  These flags must be
+bit-exclusive of those used above for the \code{psMetadataTypes}.  The
+flags have the following meanings:
+\code{PS_META_DEFAULT}: This is the zero bit mask, to allow the
+default behavior for \code{psMetadataAddItem} above.  If this is OR-ed
+with a \code{psMetadataType}, the result is as if no OR-ing took
+place.
+\code{PS_META_REPLACE}: If the given metadata item exists in the
+metadata list, and is not of type \code{PS_META_LIST} or
+\code{PS_META_META} (ie, not a container type), then this entry is
+allowed to replace the existing entry.  If this mode bit is not set, a
+duplicate (non-container-type) entry is an error.
+\begin{verbatim}
+typedef enum {                          ///< option flags for psMetadata functions
+    PS_META_DEFAULT,                    ///< default behavior (0x0000) for use in mode above
+    PS_META_REPLACE,                    ///< allow entry to be replaced
+} psMetadataFlags;
+\end{verbatim}
+An example of code to use these metadata APIs to generate the
+structure seen in Figure~\ref{fig:metadata} is given below.
+\begin{verbatim}
+md = psMetadataAlloc();
+psMetadataAdd(md, PS_LIST_TAIL, "SIMPLE",   PS_META_BOOL, "basic fits",            TRUE);
+psMetadataAdd(md, PS_LIST_TAIL, "BLANK",    PS_META_S32,  "invalid pixel data",    -32768);
+psMetadataAdd(md, PS_LIST_TAIL, "DATE-OBS", PS_META_STR,  "observing date UT", "   2004-6-16");
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_LIST, "head of comment block", NULL);
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "DATA");
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "PARAMS");
+psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_META_F32,  "exposure time (sec)",   1.05);
+psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_META_STR,  "",                      "FOO");
+cell = psMetadataAlloc();
+psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_META_STR,  "",                    "CCD00");
+psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR,  "",                    "BSEC-00");
+psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_META_STR,  "",                    "CHIP.00");
+psMetadataAdd(md,   PS_LIST_TAIL, "CELL.00",  PS_META_META, "",                    cell);
+cell = psMetadataAlloc();
+psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_META_STR,  "",                    "CCD01");
+psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR,  "",                    "BSEC-01");
+psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_META_STR,  "",                    "CHIP.01");
+psMetadataAdd(md,   PS_LIST_TAIL, "CELL.01",  PS_META_META, "",                    cell);
+\end{verbatim}
+The following code shows how to use the APIs to replace one of these values:
+\begin{verbatim}
+psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_META_F32 | PS_REPLACE,  "new exposure time (sec)",   2.05);
+\end{verbatim}
+Items may be removed from the metadata by specifying a key or a
+location in the list.  If the value of \code{name} is \code{NULL}, the
+value of \code{location} is used.  If the value of \code{name} is not
+\code{NULL}, then \code{location} must be set to
+\code{PS_LIST_UNKNOWN}.  If the key matches a metadata item, the item
+is removed from the metadata and \code{true} is returned; otherwise,
+\code{false} is returned.  If the key is not unique, then \emph{all}
+items corresponding to the key are removed, and \code{true} is
+returned.
+%
+\begin{verbatim}
+bool psMetadataRemove(psMetadata *md, int location, const char *key);
+\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 *md, const char *key);
+\end{verbatim}
+Several utility functions are provided for simple cases.  These
+functions perform the effort of casting the data to the appropriate
+type.  The numerical functions shall return 0.0 if their key is not
+found.  If the pointer value of \code{status} is not \code{NULL}, it
+is set to reflect the success or failure of the lookup.
+\begin{verbatim}
+void *psMetadataLookupPtr(bool *status, const psMetadata *md, const char *key);
+psS32 psMetadataLookupS32(bool *status, const psMetadata *md, const char *key);
+psF64 psMetadataLookupF64(bool *status, const psMetadata *md, const char *key);
+\end{verbatim}
+Items may be retrieved from the metadata by their entry position.  The
+value of which specifies the desired entry in the fashion of
+\code{psList}.
+\begin{verbatim}
+psMetadataItem *psMetadataGet(const psMetadata *md, int location);
+\end{verbatim}
+The metadata list component may be iterated over by using a
+\code{psListIterator} in a fashion equivalent to the usage for
+\code{psList}.  The iterator may be set to a location in the
+\code{psMetadata} list, and the user may get the previous or next item
+in the list relative to that location.  \code{psMetadataGetNext} has
+the ability to match the key using a POSIX regex, 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 in the metadata list, even those that are
+contained in a \code{PS_META_LIST}.  The value \code{iterator}
+specifies the iterator to be used.  In setting the iterator, the
+position of the iterator is defined by \code{location}, which follows
+the conventions of the \code{psList} iterators.
+\begin{verbatim}
+psListIterator *psMetadataIteratorAlloc(psMetadata *md, int location, bool mutable);
+bool psMetadataIteratorSet(psListIterator *iterator, int location);
+psMetadataItem *psMetadataGetAndIncrement(psListIterator *iterator, const char *regex);
+psMetadataItem *psMetadataGetAndDecrement(psListIterator *iterator, const char *regex);
+\end{verbatim}
+Metadata items may be printed to an open file descriptor based on a
+provided format.  The format string is an sprintf format statement
+with exactly one \% formatting command.  If the metadata item type is
+a numeric type, this formatting command must also be numeric, and type
+conversion performed to the value to match the format type.  If the
+metadata item type is a string, the formatting command must also be
+for a string (\%s type of command).  If the metadata type is any other
+data type, printing is not allowed.
+\begin{verbatim}
+bool psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *md);
+\end{verbatim}
+\subsubsection{Configuration files}
+\label{sec:configspec}
+It will be necessary for the \PS{} system, in order to load
+pre-defined settings, to parse a configuration file into a
+\code{psMetadata} structure.  This shall be performed by the
+function \code{psMetadataParseConfig}, as described below.
+\begin{verbatim}
+psMetadata *psMetadataParseConfig(psMetadata *md, int *nFail, const char *filename, bool overwrite);
+\end{verbatim}
+Given a metadata container, \code{md}, and the name of a configuration
+file, \code{filename}, \code{psMetadataParseConfig} shall parse the
+configuration file, placing the contained key/type/value/comment quads
+into the metadata, and returning a pointer to the metadata structure.
+The number of lines that failed to parse is returned in \code{nFail}.
+Multiple specifications of a key that haven't been declared (see
+below) are overwritten if and only if \code{overwrite} is \code{true}.
+If the metadata container is \code{NULL}, it shall be allocated.
+On error, the function shall return \code{NULL}.
+The configuration file shall consist of plain text with
+key/type/value/comment quads on separate lines.  Blank lines,
+including those consisting solely of whitespace (both spaces and
+tabs), shall be ignored, as shall lines that commence with the comment
+character (a hash mark, \code{#}), either immediately at the start of
+the line, or preceded by whitespace.  The key/type/value/comment quads
+shall all lie on a single line, separated by whitespace.
+The key shall be first, possibly preceded on the line by whitespace
+which should not form part of the key.
+Next, to assist the casting of the value, shall be a string
+identifying the type of the value, which shall correspond to one of
+the simple types supported in \code{psMetadata}:
+\code{STRING,BOOL,S32,F32,F64}; \code{STR} may be used to abbreviate
+\code{STRING}.
+\tbd{May, in the future, require more types, including U8,S16,C64,
+which will also necessitate updating the definition of psMetadata.}
+The value shall follow the type: strings may consist of multiple
+words, and shall have all leading and trailing whitespace removed;
+booleans shall simply be either \code{T} or \code{F}.
+Following the value may be an optional comment, preceded by a comment
+character (a hash mark, \code{#}), which in the case of a string
+value, serves to mark the end of the value, and for other types serves
+to identify the comment to the reader.  Only one comment character may
+be present on any single line (i.e., neither strings nor comments are
+permitted to contain the comment character).  The comment may consist
+of multiple words, and shall have leading and trailing whitespace
+removed.
+One wrinkle is the specification of vectors.  Keys for which the value
+is to be parsed as a vector shall be preceded immediately by a
+``vector symbol'', which we choose to be the ``at'' sign, \code{@}.
+In this case, the type shall be interpreted as the type for the
+vector, which may be any of the signed or unsigned integer or floating
+point types (\code{U8,U16,U32,U64,S8,S16,S32,S32,S64,F32,F64}) but not
+the complex floating point types; and the value shall consist of
+multiple numbers, separated either by a comma or whitespace.  These
+values shall populate a \code{psVector} of the appropriate type in the
+order in which they appear in the configuration file.
+\tbd{May add complex types, likely to be specified with values such as
+.23+4.56i in the future.}
+An additional hurdle is the specification of keys that may be
+non-unique (such as the \code{COMMENT} keyword in a FITS header).
+These keys shall be specified in the configuration file as non-unique
+by specifying the key at the start of the line (possibly preceded by
+whitespace) and specifying the type as a ``multiple symbol'', which we
+choose to be an asterisk, \code{*}.  No other data may be provided on
+this line, though a comment, preceeded by the comment marker, is
+valid.  A warning shall be produced when a key which has not been
+specified to be non-unique is repeated; in this case, the former value
+shall be overwritten if \code{overwrite} is \code{true}, otherwise the
+line shall be ignored and counted as one that could not be parsed.
+If a line does not conform to the rules laid out here, a warning shall
+be generated, it shall be ignored and counted as a line that could not
+be parsed.  The total number of lines that were not able to be parsed
+(including those that were ignored because \code{overwrite} is
+\code{false}, and any other parsing problems, but not including blank
+lines and comment lines) shall be returned by the function in the
+argument \code{nFail}.
+Here are some examples of lines of a valid configuration file:
+\filbreak
+\begin{verbatim}
+Double     F64     1.23456789      # This is a comment
+Float    F32 0.98765 # This is a comment too
+String  STR This is the string that forms the value #comment
+ # This is a comment line and is to be ignored
+boolean     BOOL    T # The value of `boolean' is `true'
+@primes U8  2,3 5 7,11,13 17 #   These are prime numbers
+comment MULTI # The rest of this line is ignored, but `comment' is set to be non-unique
+comment STR This
+comment STR     is
+comment STR       a
+comment STR        non-unique
+comment STR                  key
+Float F64 1.23456 # This generates a warning, and, if `overwrite' is `false', is ignored
+\end{verbatim}
+Of course, a real configuration file should look much nicer to humans
+than the above example, but PSLib must be able to parse such ugly
+files.
+We extend \code{psMetadataParseConfig} to allow a modest tree
+structure by defining a reserved keyword \code{TYPE}.  Any line in the
+config file which starts with the word \code{TYPE} shall be
+interpretted as defining a new valid type.  The defined type name
+follows the word \code{TYPE}, and is in turn followed by an arbitrary
+number of words.  These words are to be interpreted as the names of an
+embedded \code{psMetadata} entry, where the values are given on any
+line which (following the \code{TYPE} definition) employs the new type
+name.  For example, a new type may be defined as:
+\begin{verbatim}
+TYPE      CELL   EXTNAME   BIASSEC  CHIP
+CELL.00   CELL   CCD00     BSEC-00  CHIP.00
+CELL.01   CELL   CCD01     BSEC-01  CHIP.00
+\end{verbatim}
+When \code{psMetadataParseConfig} encounters the \code{TYPE} line, it
+should construct a \code{psMetadata} container and fill it with
+\code{psMetadataItems} having the names \code{EXTNAME, BIASSEC, CHIP},
+with type \code{PS_META_STR}, but data allocated.  When it next
+encounters an entry of type \code{CELL}, it should then use the given
+name (e.g., \code{CELL.00}) for the \code{psMetadataItem}, and copy
+the \code{psMetadata} data onto the \code{psMetadataItem.data.md}
+entry, filling in the values from the rest of the line (\code{CCD00,
+BSEC-00, CHIP.00}).  This hierarchical structure is illustrated in
+Figure~\ref{fig:metadata}.
+We further extend \code{psMetadataParseConfig} to allow the definition
+of a \code{psMetadata} entry using a sequence of successive lines to
+define the values of the \code{psMetadataItem} entries.  The initial
+line defines the new \code{psMetadata} entry and its name.  The
+following lines have the same format as the other metadata config file
+entries.  The sequence is terminated with a line with a single word
+\code{END}.  For example, a metadata entry may be defined as:
+\begin{verbatim}
+CELL      METADATA
+ EXTNAME   STR   CCD00
+ BIASSEC   STR   BSEC-00
+ CHIP      STR   CHIP.00
+ NCELL     S32   24
+END
+\end{verbatim}
+A series of test inputs is contained in
+\S\ref{sec:configtest}.
+\subsection{XML Functions}
+Within Pan-STARRS, we will use XML documents as a transport mechanism
+to carry data between programs and between IPP and other subsystems.
+Configuration information may be stored as well as XML documents, in
+addition to the text format discussed in the discussion on Metadata.
+XML is an extremely variable document format, and it is not currently
+the intention of PSLib to provide a complete PSLib version of XML
+operations.  Rather, a limited number of operations are defined to
+convert specific data structures to an appropriate XML document.  To
+maximize the simplicity of the XML APIs, we will use the convention
+that a single XML document to be parsed by PSLib shall contain only a
+single data structure.  Each of the XML APIs therefore take as input a
+reference to a complete XML document and return a PSLib data
+structure, or take a PSLib data structure and return a complete XML
+document.
+We start by defining a PSLib wrapper type which is a pointer to an XML
+document in memory.  We wrap the \code{libxml2} version of an XML
+document pointer for now:
+\begin{verbatim}
+typedef xmlDocPtr psXMLDoc;
+void psXMLDocFree(psXMLDoc *doc);
+\end{verbatim}
+The next pair of functions convert a \code{psMetadata} data structure
+to a complete \code{psXMLDoc} (in memory) and vice versa:
+\begin{verbatim}
+psXMLDoc *psMetadataToXMLDoc(const psMetadata *metadata);
+psMetadata *psXMLDocToMetadata(const psXMLDoc *doc);
+\end{verbatim}
+The next pair of functions loads the data in a named file into a
+complete \code{psXMLDoc} (in memory) and write out the \code{psXMLDoc}
+to a named file:
+\begin{verbatim}
+psXMLDoc *psXMLParseFile(const char *filename);
+int psXMLDocToFile(const psXMLDoc *doc, const char *filename);
+\end{verbatim}
+The next pair of functions accepts a block of memory and parses it
+into a complete \code{psXMLDoc} (also in memory), and vice versa:
+\begin{verbatim}
+psXMLDoc *psXMLParseMemory(const char *buffer, const int size);
+int psXMLDocToMemory(const psXMLDoc *doc, char *buffer);
+\end{verbatim}
+The next pair of functions read from and write to a file descriptor.
+The first converts the imcoming data to a complete \code{psXMLDoc}
+(also in memory), the second writes the \code{psXMLDoc} to the file
+descriptor:
+\begin{verbatim}
+psXMLDoc *psXMLParseFD(int fd);
+int psXMLDocToFD(const psXMLDoc *doc, int fd);
+\end{verbatim}
+\subsection{Database Functions}
+Many of the applications that PSLib will be used for will require
+access to a simple relational database.  PSLib includes generic
+database-independent interface mechanisms as part of its API set.  The
+most important aspect of PSLib's database support is to abstract as
+much database specific complexity as is feasible.  As almost all RDBMS
+provide at least a simple transactional model, commit and rollback
+support should be provided.
+Currently, only support for MySQL 4.1.x is required but other backends
+may be added as options in the future.  As a particular example which
+has implications for the database interaction model, support for
+SQLite may be required in the future.  Currently, the choice of
+backend database interface may be made as a compile option.  Details
+of the specified APIs in the discussion below refer to the relevant
+MySQL functions.
+Database errors must be trapped and placed onto the psError stack.
+The complete error message should be retrieved with the database's
+error function.
+\subsubsection{Managing the Database Connection}
+We specify a database handle which carries the information about the
+database connection:
+\begin{verbatim}
+    typedef struct {
+        MYSQL *mysql;
+    } psDB;
+\end{verbatim}
+The following collection of functions provides basic database functionality:
+\begin{verbatim}
+    // wraps mysql_init() & mysql_real_connect()
+    psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname);
+    // wraps mysql_close()
+    void psDBCleanup(psDB *dbh);
+    // wraps mysql_create_db()
+    bool psDBCreate(psDB *dbh, const char *dbname);
+    // wraps mysql_select_db()
+    bool psDBChange(psDB *dbh, const char *dbname);
+    // wraps mysql_drop_db()
+    bool psDBDrop(psDB *dbh, const char *dbname);
+\end{verbatim}
+For MySQL support, \code{psDBInit()} wraps \code{mysql_init()} and
+\code{mysql_real_connect()} in order to initialize a psDB structure and
+establish a database connection.  A null pointer should be returned on
+failure.
+When implementing support for SQLite, or other DB which is purely
+file-based, the \code{host}, \code{user}, and \code{passwd} arguments
+would be ignored while \code{dbname} would specify the path to the
+SQLite db file.
+\subsubsection{Interacting with Database Tables}
+The functions in this section perform high level interactions with the
+database tables.  All of them should behave ``atomically'' with
+respect to the state of the database.  Specifically, all interactions
+with the database should be done as a part of a transaction that is
+rolled-back on failure and committed only after all queries used by
+the API have been run.  In general, this API set attempts to treat a
+database table as a 2D matrix where columns can be represented by a
+\code{psVector} and rows as a \code{psMetadata} type.  A
+\code{psMetadata} collection is also used to define the columns of a
+table and as part of the query restrictions.
+\begin{verbatim}
+    bool psDBCreateTable(psDB *dbh, const char *tableName, psMetadata *md);
+\end{verbatim}
+This function generates and executes the SQL needed to create a table
+named \code{tableName}, with the column names and datatypes as
+described in \code{md}.  Each data item in the \code{psMetadata}
+collection represents a single table field.  The name of the field is
+given by the name of the \code{psMetadataItem} and the data type is
+give by the \code{psMetadataItem.type} and \code{psMetadataItem.ptype}
+entries.  A lookup table should be used to convert from PSLib types
+into MySQL compatible SQL data types.  For example, a
+\code{PS_META_STR} would map to an SQL99 varchar.  If the value of
+\code{type} is \code{PS_META_STR} then the \code{psMetadataItem.data}
+element is set to a string with the length for the field written as a
+text string.  The value of the \code{psMetadataItem.data} element is
+unused for the \code{PS_META_PRIMITIVE} types.  Other metadata types
+beyond \code{PS_META_STR} and \code{PS_META_PRIMITIVE} are not allowed
+in a table definition metadata collection.
+Database indexes can be specified setting the \code{comment} field to
+``\code{Primary Key}'' or ``\code{Key}''.  Comment are otherwise
+ignored.
+\begin{verbatim}
+    bool psDBDropTable(psDB *dbh, const char *tableName);
+\end{verbatim}
+This function deletes the specified table.
+\begin{verbatim}
+    psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, const psU64 limit);
+    psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col, psElemType pType, const psU64 limit);
+\end{verbatim}
+These functions generates and executes the SQL needed to select an entire
+column from a table or up to \code{limit} rows from it.  If \code{limit} is 0,
+the entire range is returned.  The database response is processed and a
+\code{psArray} of strings is returned.  The Num version of the function returns
+the data in a \code{psVector}, data cast to \code{pType}.  It returns an error
+(NULL) if the requested field is not a numerical type.
+\begin{verbatim}
+    psArray *psDBSelectRows(psDB *dbh, const char *tableName, psMetadata *where, const psU64 limit);
+\end{verbatim}
+This function returns rows from the specified table which match
+the restrictions given by \code{where}.  The restrictions are
+specified as field / value pairs.  The \code{psMetadata} collection
+where must consist of valid database fields, though the database query
+checking functions may be used to validate the fields as part of the
+query.  If \code{where} is \code{NULL}, then there are no restrictions
+on the rows selected.  The selected rows are returned as a
+\code{psArray} of \code{psMetadata} values, one per row.
+\begin{verbatim}
+    bool psDBInsertOneRow(psDB *dbh, const char *tableName, psMetadata *row);
+\end{verbatim}
+Insert the data from \code{row} into \code{tableName}.  It should be noted in
+the API reference that if fields are specified in \code{row} that do not exist
+in \code{tablename}, the insert will fail.
+\begin{verbatim}
+    bool psDBInsertRows(psDB *dbh, const char *tableName, psArray *rowSet);
+\end{verbatim}
+Similar to \code{psDBInsertOneRow()}, this function inserts many rows at once
+and is atomic for the complete set of rows.
+\begin{verbatim}
+    psArray *psDBDumpRows(psDB *dbh, const char *tableName);
+\end{verbatim}
+Fetch all rows as an psArray of psMetadata.
+\begin{verbatim}
+    psMetadata *psDBDumpCols(psDB *dbh, const char *tableName);
+\end{verbatim}
+Fetch all columns, as either a psVector or a psArray depnding on whether or not
+the column is numeric, and return them in a psMetadata structure where
+psMetadataItem.name contains the column's name.
+\begin{verbatim}
+    psS64 psDBUpdateRows(psDB *dbh, const char *tableName, psMetadata *where, psMetadata *values);
+\end{verbatim}
+Update the columns contained in \code{values} in the row(s) that have a field
+with the value indicated by \code{where} (note that this is only allows very
+limited use of SQL99's ``where'' semantics).  The number of rows modified is
+returned.  A negative value is return to indicate an error. If there are
+multiple psMetadataItems in \code{where} then each item should be considered as
+an additional constraint.  e.g.  ``where foo = x and where bar = y''
+\begin{verbatim}
+    psS64 psDBDeleteRows(psDB *dbh, const char *tableName, psMetadata *where);
+\end{verbatim}
+Delete the rows that are matched by \code{where} using the same semantics for
+\code{where} as in psDBUpdateRow().  A negative value is returned to indicate an
+error.
+\subsection{FITS I/O Functions}
+We need a variety of I/O functions between the disk and certain of our
+PSLib data structures.  We need the ability to access FITS headers,
+images and tables (both ASCII and Binary).  We define here the FITS
+I/O functions, all of which are currently specified as wrappers to
+functions within CFITSIO.  CFITSIO provides a wide range of utilities
+which we do not feel are particularly appropriate as part of a generic
+I/O library, such as assumptions about names which change the data
+interpretation, etc.  We are defining our calls to avoid the hidden
+'features'.  The CFITSIO functions which are wrapped should in general
+be the most basic versions.
+\begin{verbatim}
+typedef struct {
+    fitsfile fd;
+} psFits;
+\end{verbatim}
+We begin by defining a datatype to wrap the CFITSIO \code{fitsfile}
+structure.  This is necessary to allow repeated access to the data in
+a file without multiple open commands (which are expensive).
+\subsubsection{FITS File Manipulations}
+\begin{verbatim}
+psFits *psFitsAlloc(const char *filename);
+\end{verbatim}
+Opens a FITS file at positions the pointer to the PHU.
+\begin{verbatim}
+bool psFitsMoveExtName(psFits *fits, const char *extname);
+\end{verbatim}
+Positions the pointer to the beginning of the specified
+\code{extname}.  If the \code{extname} does not exist, the function
+shall fail.
+\begin{verbatim}
+bool psFitsMoveExtNum(psFits* fits, int extnum, bool relative);
+\end{verbatim}
+Moves the pointer to the beginning of the specified HDU number.  If
+\code{relative} is TRUE, \code{extnum} represents the number of HDUs
+relative to the current HDU.  The PHU is entry number 0, while the
+extended data segments start at number 1.
+\begin{verbatim}
+int psFitsGetExtNum(psFits* fits);
+\end{verbatim}
+Returns the current HDU number (i.e., file position).
+\begin{verbatim}
+int psFitsGetSize(psFits* fits);
+\end{verbatim}
+Returns the number of HDUs in the file.
+\begin{verbatim}
+psFitsType psFitsGetExtType(psFits* fits);
+\end{verbatim}
+Gets the current HDU's type (table or image).
+\subsubsection{FITS Header I/O Functions}
+\begin{verbatim}
+psMetadata *psFitsReadHeader(psMetadata *out, const psFits *fits);
+\end{verbatim}
+Read header data into a \code{psMetadata} structure.  The data is read
+from the current HDU pointed at by the \code{psFits *fits} entry.  If
+\code{out} is \code{NULL}, a new psMetadata is created.
+\begin{verbatim}
+psMetadata *psFitsReadHeaderSet (psFits *fits);
+\end{verbatim}
+Load a complete set of headers from the \code{psFits} file pointer.
+This function loads the headers from all extensions into a
+\code{psMetadata} collection, each entry of which is a pointer to a
+\code{psMetadata} structure containing the header data.  The metadata
+entry names are the \code{EXTNAME} values for each header (with the
+value of \code{PHU} for the primary header unit).  At the start of the
+operation, the file pointer is rewound to the beginning of the file.
+At the end, it is positioned where it started when the function was
+called.
+\begin{verbatim}
+bool psFitsWriteHeader(psMetadata *output, const psFits *fits);
+\end{verbatim}
+Write metadata into the header of a FITS image file.  The header is
+written at the current HDU.
+\subsubsection{FITS Image I/O Functions}
+\begin{verbatim}
+psImage *psFitsReadImage(psImage *output, psFits *fits, psRegion region, int z);
+\end{verbatim}
+Read an image or subimage from the \code{psFits} file pointer.  This
+function is a wrapper to the CFITSIO library function.  The input
+parameters allow a full image or a subimage to be read.  The region to
+be read is specified by \code{region}.  A negative value for either of
+\code{region.x1} or \code{region.y1} specifies the size of the region
+to be read counting down from the end of the array.
+If the native image is a cube, the value of z specifies the requested
+slice of the image.  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, or if the image on disk is zero- or
+one-dimensional.  This function need only read images of the native
+FITS image types (\code{psU8}, \code{psS16}, \code{psS32},
+\code{psF32}, \code{psF64}).  The user is expected to convert the data
+type as needed with \code{psImageCopy}.
+\begin{verbatim}
+bool psFitsUpdateImage(psFits *fits, const psImage *input, psRegion region, int z);
+\end{verbatim}
+\tbd{we have discussed this as the alternate name}
+Write an image section to the open \code{psFits} file pointer.  This
+operation may write a portion of an image over the existing bytes of
+an existing image.  Care must be taken to interpret \code{region},
+which specified the output pixels to be written / over-written.  If
+the combination of \code{region} and the size of \code{psImage *input}
+implies writing pixels outside the existing data area of the image,
+the function shall return an error (ie, if \code{region.x0 + image.nx
+>= NAXIS1}, \code{region.y0 + image.ny >= NAXIS2}, or \code{z >=
+NAXIS3}).  This function will only write images of the native FITS
+image types (\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32},
+\code{psF64}).  The user is expected to convert the data type as
+needed with \code{psImageCopy}.  The return value must be 0 for a
+successful operation and 1 for an error.
+\begin{verbatim}
+bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input, int depth);
+\end{verbatim}
+Create a new image based on the dimensions specified for the image and
+the requested depth.  The header and image data segments are written
+in the file at the current position of the \code{psFits} pointer.
+This function will only write images of the native FITS image types
+(\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, \code{psF64}).
+The user is expected to convert the data type as needed with
+\code{psImageCopy}.  The return value must be 0 for a successful
+operation and 1 for an error.
+\subsubsection{FITS Table I/O Functions}
+\begin{verbatim}
+psMetadata *psFitsReadTableRow (psFits *fits, int row);
+\end{verbatim}
+This function reads a single row of the table in the extension pointed
+at by the \code{psFits} file pointer.  The row number to be read is
+given by \code{row}.  The result is returned as a \code{psMetadata}
+collection with elements of the apporpriate types and keys
+corresponding to the table column names.  The function must apply the
+needed byte-swapping on the data in the row based on the description
+of the table data in the table header.  \tbr{we may need to be more
+flexible here: if we call this function repeatedly, it would be more
+efficient to pass the corresponding header or keep it somewhere (and
+the file pointer location, for that matter).}
+\begin{verbatim}
+void *psFitsReadTableRowRaw (int *nBytes, psFits *fits, int row);
+\end{verbatim}
+This function reads a single row of the table in the extension pointed
+at by the \code{psFits} file pointer.  The row number to be read is
+given by \code{row}.  The result is returned as collection of
+\code{nBytes} bytes allocated by the function.  The function must
+apply the needed byte-swapping on the data in the row based on the
+description of the table data in the table header.  \tbr{we may need
+to be more flexible here: if we call this function repeatedly, it
+would be more efficient to pass the corresponding header or keep it
+somewhere (and the file pointer location, for that matter).}
+\begin{verbatim}
+psArray *psFitsReadTableColumn (psFits *fits, char *colname);
+\end{verbatim}
+This function reads a single column of the table in the extension
+pointed at by the \code{psFits} file pointer.  The column is specified
+by the FITS table column key given by \code{row}.  The result is
+returned as a \code{psArray}, with the data from one row of the table
+column per array element.
+\begin{verbatim}
+psVector *psFitsReadTableColumnNum (psFits *fits, char *colname);
+\end{verbatim}
+This function reads a single column of the table in the extension
+pointed at by the \code{psFits} file pointer.  The column is specified
+by the FITS table column key given by \code{row} and must be of a
+numeric data type.  The result is returned as a \code{psVector} of the
+appropriate data type, with the data from one row of the table column
+per array element.
+\begin{verbatim}
+psArray *psFitsReadTableRaw (int *nBytes, psFits *fits);
+\end{verbatim}
+This function reads the entire data block from a table into the a
+\code{psArray}, with one element of the array per row.  The number of
+bytes per row is returned in \code{nBytes}.  The function must apply
+the needed byte-swapping on the data in each row based on the
+description of the table data in the table header.
+\begin{verbatim}
+psArray *psFitsReadTable (psFits *fits);
+\end{verbatim}
+This function reads the entire data block from a table into the a
+\code{psArray}, with one element of the array per row.  Each row is
+stored as a \code{psMetadata} collection as described above for
+\code{psFitsReadTableRow}.
+\begin{verbatim}
+bool psFitsWriteTable(psFits* fits, psMetadata *header, psArray* table);
+\end{verbatim}
+Accepts a \code{psArray} of \code{psMetadata} and writes it to the
+current HDU.  If the current HDU is not a table type, this will fail
+and return FALSE.
+\begin{verbatim}
+bool psFitsUpdateTable(psFits* fits, psMetadata *header, psMetadata* data, int row);
+\end{verbatim}
+Writes the \code{psMetadata} data to a FITS table at the specified row
+in the current HDU.  If the current HDU is not a table type, this will
+fail and return FALSE.
+\subsection{Detector and Sky Coordinates}
+\subsection{Linear and Spherical Coordinates}
 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.
+while the second are angular coordinates which define a position on
+the sphere of the sky.  We put these into two structures,
+\code{psPlane} and \code{psSphere}, respectively.  Partitioning these
+two will enable error-checking.  An alternative representation for
+angular positions is the 3-D unit vector.  These are used in
+particular as part of spherical rotation calculations.  We define
+\code{psCube} to represent such an element.
+%
 \begin{verbatim}
 …
     double dErr;                        ///< Error in Dec
 } psSphere;
+typedef struct {
+    double x;                           ///< cos (DEC) cos (RA)
+    double y;                           ///< cos (DEC) sin (RA)
+    double z;                           ///< sin (DEC)
+    double xErr;                        ///< Error in x
+    double yErr;                        ///< Error in y
+    double zErr;                        ///< Error in z
+} psCube;
 \end{verbatim}
 …
 Constructors for these are straight-forward:
 \begin{verbatim}
 psPlane *psPlaneAlloc(void);
+psPlane  *psPlaneAlloc(void);
 psSphere *psSphereAlloc(void);
+psCube   *psCubeAlloc(void);
 \end{verbatim}
 Initialization of the structures is not necessary.
 …
 \code{NULL}, the function shall generate an error and return
 \code{NULL}.
-\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 cosDeltaP;                 ///< Cosine of target pole latitude in the source system
-    double sinDeltaP;                 ///< Sine of target pole latitude in the source system
-    double alphaP;                    ///< Longitude of the target system pole in the source system
-    double phiP;                      ///< Longitude of the ascending node in the target system
-} psSphereTransform;
-\end{verbatim}
-The constructor is defined as follows:
-\begin{verbatim}
-psSphereTransform *psSphereTransformAlloc(double alphaP, double deltaP, double phiP);
-\end{verbatim}
-where \code{alphaP} and \code{deltaP} define the coordinates in the
-input system of the north pole in the output system and \code{phiP}
-defines the longitude in the input system of the equatorial
-intersection between the two systems (e.g, the first point of Ares).
-The constructor must calculate the sines and cosines above.
-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., ICRS, Ecliptic and Galactic).  These are
-constructors as well as the above \code{psSphereTransformAlloc}.
+%
-\begin{verbatim}
-psSphereTransform *psSphereTransformICRSToEcliptic(psTime *time);
-psSphereTransform *psSphereTransformEclipticToICRS(psTime *time);
-psSphereTransform *psSphereTransformICRSToGalactic(void);
-psSphereTransform *psSphereTransformGalacticToICRS(void);
-\end{verbatim}
-We also require the ability to precess coordinates from one equinox to
-another.
-\begin{verbatim}
-psSphere *psSpherePrecess(psSphere *coords, const psTime *fromTime, const psTime *toTime);
-\end{verbatim}
-Given coordinates, \code{coords}, with equinox for \code{fromTime},
-the coordinates are precessed to equinox for \code{toTime}.  The
-\code{coords} are modified in-place.  Equinoxes shall be Julian
-equinoxes (as opposed to Bessellian).
-\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
-\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_NTYPE                       ///< Number of types; must be last
-} psProjectionType;
-\end{verbatim}
-The constructor is straight-forward:
-\begin{verbatim}
-psProjection *psProjectionAlloc(double R, double D, double Xs, double Ys, psProjectionType type);
-\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.  Both an offset
-mode and an offset unit may be defined.  The mode may be either
-\code{PS_SPHERICAL}, in which case the specified offset corresponds to
-an offset in angles, or it may be \code{PS_LINEAR}, in which case the
-offset corresponds to a linear offset in a local projection.  The
-offset unit may be in one of \code{PS_ARCSEC}, \code{PS_ARCMIN},
-\code{PS_DEGREE}, and \code{PS_RADIAN}, which specifies the units of
-the offset only.
-\begin{verbatim}
-psSphere *psSphereGetOffset(const psSphere *position1,
-                            const psSphere *position2,
-                            psSphereOffsetMode mode,
-                            psSphereOffsetUnit unit);
-psSphere *psSphereSetOffset(const psSphere *position,
-                            const psSphere *offset,
-                            psSphereOffsetMode mode,
-                            psSphereOffsetUnit unit);
-typedef enum {
-    PS_SPHERICAL;                       ///< Offset on a sphere
-    PS_LINEAR;                          ///< Linear offset
-} psSphereOffsetMode;
-typedef enum {
-    PS_ARCSEC;                          ///< Arcseconds
-    PS_ARCMIN;                          ///< Arcminutes
-    PS_DEGREE;                          ///< Degrees
-    PS_RADIAN;                          ///< Radians
-} psSphereOffsetUnit;
-\end{verbatim}
-Note that these should propagate the errors appropriately.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%%% Astronomical Images and Astrometry
-\include{psLibSDRS_Astrom}
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{Pixel Lists}
-Usually an image mask is the best way to carry information about what
-pixels mean what.  However, in the case where the number of pixels in
-which we are interested is limited, it is more efficient to simply
-carry a list of pixels.  An example of this is in the image
-combination code, where we want to perform an operation on a
-relatively small fraction of pixels, and it is inefficient to go
-through an entire mask image checking each pixel.
-\begin{verbatim}
-typedef struct {
-    psVector *x;                        // x coordinate
-    psVector *y;                        // y coordinate
-} psPixels;
-\end{verbatim}
-Of course, the size of each of the vectors should match.  In the event
-that they do not match, any function which detects the problem shall
-generate a warning and use the size of the shorter of the vectors as
-the size.  The order in which the pixels are kept is not considered
-important.
-\begin{verbatim}
-psImage *psPixelsToMask(psImage *out, const psPixels *pixels, const psRegion *region, unsigned int maskVal);
-psPixels *psMaskToPixels(psPixels *out, const psImage *mask, unsigned int maskVal);
-\end{verbatim}
-\code{psPixelsToMask} shall return an image of type U8 with the
-\code{pixels} lying within the specified \code{region} set to the
-\code{maskVal}.  The \code{out} image shall be modified if supplied,
-or allocated and returned if \code{NULL}.  The size of the output
-image shall be \code{region->x1 - region->x0} by \code{region->y1 -
-region->y0}, with \code{out->x0 = region->x0} and \code{out->y0 =
-region->y0}.  In the event that either of \code{pixels} or
-\code{region} are \code{NULL}, the function shall generate an error
-and return \code{NULL}.
-\code{psMaskToPixels} shall return a \code{psPixels} consisting of the
-coordinates in the \code{mask} that match the \code{maskVal}.  The
-\code{out} pixel list shall be modified if supplied, or allocated and
-returned if \code{NULL}.  In hte event that \code{mask} is
-\code{NULL}, the function shall generate an error and return
-\code{NULL}.
-\begin{verbatim}
-psPixels *psPixelsConcatenate(psPixels *out, const psPixels *pixels);
-\end{verbatim}
-\code{psPixelsConcatenate} shall concatenate \code{pixels} onto
-\code{out}.  In the event that \code{out} is \code{NULL}, a new
-\code{psPixels} shall be allocated, and the contents of \code{pixels}
-simply copied in.  If \code{pixels} is \code{NULL}, the function shall
-generate an error and return \code{NULL}.  The function shall take
-care to ensure that there are no duplicate pixels in \code{out} (since
-the order in which the pixels are stored is not important, the values
-may be sorted, allowing the use of a faster algorithm than a linear
-scan).
 \begin{verbatim}
 …
 and return \code{NULL}.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsubsection{Spherical Rotations}
+Spherical rotations represent coordinate transformation in 3-D, as
+well as the effects of precession and nutation.  We need spherical
+rotatations to convert between ICRS, Galactic and Ecliptic
+coordinates, and to determine Alt-Az coordinates for sources.  All of
+these basic spherical transformations represent rotations of the
+spherical coordinate reference.  We specify a general transformation
+function which takes a structure, \code{psSphereRot}, defining the
+transformation between two spherical coordinate systems.  The
+structure contains the elements of a quaternion to represent the
+spherical rotational.  We define two allocators for
+\code{psSphereRot}, one which defines the rotation in terms of the
+coordinate of the pole and the rotation about that pole.  The other
+defines the rotation from the elements of the quaternion.  We also
+specify functions to manipulate \code{psSphereRot} in several useful
+way.
+\begin{verbatim}
+typedef struct {
+    double q0;
+    double q1;
+    double q2;
+    double q3;
+} psSphereRot;
+\end{verbatim}
+The constructor is defined as follows:
+\begin{verbatim}
+psSphereRot *psSphereRotAlloc(double alphaP, double deltaP, double phiP);
+\end{verbatim}
+where \code{alphaP} and \code{deltaP} define the coordinates in the
+input system of the axis of rotation (the north pole of the output
+system), while \code{phiP} defines the rotation about that pole.  This
+last angle is also equal to 270\degree - $\phi_a$, where $\phi_a$ is
+the longitude in the output system of the ascending node (equatorial
+intersection between the two systems, e.g, the first point of Ares).
+The \code{psSphereRot} may also be constructed by supplying the
+elements of the quaternion to the following function:
+\begin{verbatim}
+psSphereRot *psSphereRotQuat(double q0, double q1, double q2, double q3);
+\end{verbatim}
+This function normalizes the quaternion, so the input elements need
+not be normalized.
+Spherical coordinates may be transformed by providing the
+transformation and the coordinate in the input system to
+\code{psSphereRot}.  The output pointer may be optionally supplied, or
+if \code{NULL}, is allocated by the function.
+\begin{verbatim}
+psSphere *psSphereRotApply(psSphere *out, const psSphereRot *transform, const psSphere *coord);
+\end{verbatim}
+The following function combines two rotations, to produce a single
+rotation which is the equivallent of applying the first rotation and
+then the second.  The output rotation may be supplied, or will be
+allocated if \code{NULL}.
+\begin{verbatim}
+psSphereRot *psSphereRotCombine(psSphereRot *out, psSphereRot *rot1, psSphereRot *rot2)
+\end{verbatim}
+The following function changes the given rotation to its inverse:
+\begin{verbatim}
+psSphereRot *psSphereRotInvert(psSphereRot *rot)
+\end{verbatim}
+The 3-vector representation of the angles (\code{psCube}) is needed to
+implement these functions, and is useful in other circumstances as
+well.  Two utility functions are provided to convert between the
+angular and 3-vector representations:
+\begin{verbatim}
+psSphere *psCubeToSphere(psCube *cube);
+psCube *psSphereToCube(psSphere *sphere);
+\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.  Both an offset
+mode and an offset unit may be defined.  The mode may be either
+\code{PS_SPHERICAL}, in which case the specified offset corresponds to
+an offset in angles, or it may be \code{PS_LINEAR}, in which case the
+offset corresponds to a linear offset in a local projection.  The
+offset unit may be in one of \code{PS_ARCSEC}, \code{PS_ARCMIN},
+\code{PS_DEGREE}, and \code{PS_RADIAN}, which specifies the units of
+the offset only.
+\begin{verbatim}
+psSphere *psSphereGetOffset(const psSphere *position1,
+                            const psSphere *position2,
+                            psSphereOffsetMode mode,
+                            psSphereOffsetUnit unit);
+psSphere *psSphereSetOffset(const psSphere *position,
+                            const psSphere *offset,
+                            psSphereOffsetMode mode,
+                            psSphereOffsetUnit unit);
+typedef enum {
+    PS_SPHERICAL;                       ///< Offset on a sphere
+    PS_LINEAR;                          ///< Linear offset
+} psSphereOffsetMode;
+typedef enum {
+    PS_ARCSEC;                          ///< Arcseconds
+    PS_ARCMIN;                          ///< Arcminutes
+    PS_DEGREE;                          ///< Degrees
+    PS_RADIAN;                          ///< Radians
+} psSphereOffsetUnit;
+\end{verbatim}
+Note that these should propagate the errors appropriately.
+\subsection{Celestial Coordinate Systems}
+The following functions simply return the appropriate
+\code{psSphereRot} to convert between predefined spherical
+coordinate systems (i.e., ICRS, Ecliptic and Galactic).  These are
+constructors as well as the above \code{psSphereRotAlloc}.
+%
+\begin{verbatim}
+psSphereRot *psSphereRotICRSToEcliptic(psTime *time);
+psSphereRot *psSphereRotEclipticToICRS(psTime *time);
+psSphereRot *psSphereRotICRSToGalactic(void);
+psSphereRot *psSphereRotGalacticToICRS(void);
+\end{verbatim}
+\subsection{Earth Orientation Calculations}
+One of the critical sets of calculations in astronomy is the sequence
+of steps needed to convert between the celestial coordinates of an
+object and the observed coordinates of the object.  This problem is
+best divided into two major components: transformation between the
+celestial sphere and coordinates relative to the surface of the solid
+earth, excluding the effects of the atmosphere, and compensation for
+the effects of the atmosphere.  In this section, we address the first
+of these two transformations: the Earth Orientation Calculations.
+The Earth Orientation Calculations are further subdivided into several
+steps, illustrated in Figure~\ref{CoordinateSystems} .  Celestial
+coordinates are defined in the International Celestial Reference
+System (ICRS), which has the solar barycenter as its reference
+position and velocity.  The next coordinate system is the Geocentric
+Celestial Reference System (GCRS), which uses the earth barycenter as
+a reference.  The transformation between these two includes the
+abberation due to the Earth's velocity, the parallax of the object,
+which depends on both the Earth's position and the distance to the
+object of interest, and the general relativistic correction for the
+bending of light as it approaches the Earth.
+The next set of transformations compenstate for the 3-D rotation of
+the Earth on various timescales, including the effects of precession,
+nutation, and simple solid-body rotation.  These calculations can be
+performed using different amounts of information for higher levels of
+precision.  Since the Earth's rotation is constantly affected by
+stochastic processes (weather, earthquakes, etc), these conversions
+are constantly modified by observations reported by authoratative
+sources such as the US Naval Observatory.  The target of this
+transformation is the International Terrestrial Reference System
+(ITRS), which is fixed with respect to the Earth's crust.  This
+transformation is subdivided into slow precesion and nutation
+(yielding the coordinate system CIP/CEO), followed by the Earth's
+rotation (yielding the coordinate system CIP/TEO), and finally
+corrections for the short-period motion of the Earth's pole.
+\subsubsection{Transformation from ICRS to GCRS}
+\tbd{we need a function to construct the direction and speed elements
+  given the time}.
+\tbd{supply the velocity as an un-normalized 3 vector?}
+\paragraph{Aberration}
+The following function calculates the \code{apparent} position of a
+star, given its \code{actual} position and the velocity vector of the
+observer, represented as a speed and a direction:
+\begin{verbatim}
+psAberration(psSphere *apparent, psSphere *actual, psSphere direction, double speed);
+\end{verbatim}
+The \code{actual} and \code{apparent} positions are represented as
+\code{psSphere} entries, as is the \code{direction} of motion.  The
+speed in that direction is given in units of the speed of light.
+\paragraph{Gravitational Deflection}
+\paragraph{Parallax}
+\begin{verbatim}
+double psEOC_ParallaxFactor(psSphere *coords, psTime *time);
+\end{verbatim}
+Calculate the parallax factor for the given position and time.
+\subsubsection{Transformation from GCRS to ITRS}
+\paragraph{Precession/Nutation}
+The following routine calculates the components of the rotation
+between the CEO and GCRS frames, $X$, $Y$, and $s$, using to the
+IAU2000A precession \& nutation model:
+%
+\begin{verbatim}
+psSphere *psEOC_PrecessionModel(double *s, const psTime *time)
+\end{verbatim}
+%
+The input to this function is the desired \code{time}, which may be
+represented in any format other than UT1.  This routine must give
+results identical to the IERS XYS2000A subroutine within the limits of
+machine accuracy.
+The following function provides interpolated corrections to $X$ and
+$Y$ from the tables provided by the IERS, just as it does for UT1 and
+polar motion.
+\begin{verbatim}
+psSphere *psEOC_GetPolarCorr(const psTime *time, psTimeBulletin bulletin);
+\end{verbatim}
+The polar correction is applied to the $X$ and $Y$ elements of the
+rotation to provide higher accuracy.  The spherical rotation term is
+generated by providing the three elements of the rotation to the
+following function:
+\begin{verbatim}
+psSphereRot *psSphereRot_CEOtoGCRS(double s, const psSphere *pole)
+\end{verbatim}
+The retulting \code{psSphereRot} may be used to determine the rotation
+from CIP/CEO to GCRS.  This function must give results identical to
+the IERS BPN2000, within the limits of machine accuracy.
+\paragraph{Earth Rotation}
+The following routine calculates the rotation of the Earth about the CIP:
+\begin{verbatim}
+psSphereRot *psSphereRot_TEOtoCEO(const psTime *time)
+\end{verbatim}
+The IERS code to create the comparable rotation is embedded in
+T2C2000, with the Earth Rotation Angle calculated by ERA2000.
+\paragraph{Polar Motion}
+The following function provides interpolated values of the polar
+motion components, $x_p$ and $y_p$, extracted from the IERS tables.
+\begin{verbatim}
+psSphere *psEOC_GetPoleCoords(const psTime *time, psTimeBulletin bulletin);
+\end{verbatim}
+The following function provides tidal corrections to the polar motion
+components, $x_p$ and $y_p$, using the Ray model of Simon et al (see
+ADD).
+\begin{verbatim}
+psSphere *psEOC_TidePolarCorr(const psTime *time);
+\end{verbatim}
+The following function provides the additional corrections due to nutation
+terms with periods less than or equal to two days:
+\begin{verbatim}
+psSphere *psEOC_NutationCorr(psTime *time);
+\end{verbatim}
+The following function should generate the \code{psSphereRot} transform from
+ITRS to CIP/TEO:
+\begin{verbatim}
+psSphereRot *psSphereRot_ITRStoTEO(psSphere pole, psTime *time);
+\end{verbatim}
+The time argument should be used to internally calculate $s'$.
+This function should give identical results to the IERS POM2000 subroutine.
+\subsubsection{Earth Orientation Wrappers}
+The following function generates the complete spherical rotation to
+account for precession between two times.  If \code{NULL} is provided
+for either time, it is assumed to have the reference equinox value of
+J2000.
+\begin{verbatim}
+psSphere *psSpherePrecess(const psTime *fromTime, const psTime *toTime, psPrecessMethod mode);
+\end{verbatim}
+The mode argument is used to specify the level of detail used in the
+calculation.
+\begin{verbatim}
+typedef enum {
+  PS_PRECESS_ROUGH,
+  PS_PRECESS_COMPLETE,
+  PS_PRECESS_IAU2000A,
+} psPrecessMethod;
+\end{verbatim}
+\subsection{Atmospheric Effects}
+\tbd{The ATM effects components should be deferred until we clean up
+  the refraction definitions}
+A-priori astrometric transformations between the telescope orientation
+(Alt/Az) and the predicted stellar coordinates above the atmosphere
+(DEC/HA) requires several pieces of information describing the current
+environmental conditions.  These quantities are consistent across an
+image, and may vary only slowly with time.  Pre-computing these
+quantities for exposures means that subsequent transformations are
+faster.  The structure below carries the environmental data of interest.
+For historical reasons, this structure is known colloquially as ``the
+Grommit''.
+\tbd{this structure needs to be modified to correspond to what we
+  actually need to carry around for the atmosphere functions}
+\tbd{provide a single Grommit to carry around all EOC + ATM
+  pre-calculated entries and a separate structure for ATM effect?}
+\begin{verbatim}
+typedef struct {
+    const double latitude;              ///< geodetic latitude (radians)
+    const double longitude;             ///< longitude + ... (radians)
+    const double height;                ///< height (HM)
+    const double abberationMag;         ///< magnitude of diurnal aberration vector
+    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 siderealTime;          ///< local apparent sidereal time (radians)
+} psGrommit;
+\end{verbatim}
+The \code{psGrommit} is calculated from telescope information for the
+particular exposure, \code{exp}:
+\begin{verbatim}
+psGrommit *psGrommitAlloc(const psExposure *exp);
+\end{verbatim}
+\tbd{these functions probably need to take the ATM structure}
+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, psTime *lst, float height);
+\end{verbatim}
+which returns the airmass for a given position and local sidereal time
+(\code{lst}).
+\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.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{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}
+The constructor for \code{psFixedPattern} shall be:
+\begin{verbatim}
+psFixedPattern *psFixedPatternAlloc(double x0,        double y0,
+                                    double xScale,    double yScale,
+                                    const psImage *x, const psImage *y);
+\end{verbatim}
+Here, \code{x0}, \code{y0}, \code{xScale} and \code{yScale} have the
+same meaning as in the \code{psFixedPattern} structure.  Note that the
+values of the fixed pattern offsets are specified as images, the
+values from which need to be copied into the \code{double **x} and
+\code{double **y} of \code{psFixedPattern}, and that the number of
+elements may be derived from the size of the images.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+\subsection{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
+\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_NTYPE                       ///< Number of types; must be last
+} psProjectionType;
+\end{verbatim}
+The constructor is straight-forward:
+\begin{verbatim}
+psProjection *psProjectionAlloc(double R, double D, double Xs, double Ys, psProjectionType type);
+\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}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 \subsection{Photometry}
 …
 M_{\rm pM} - pA, M_{\rm sP} - M_{\rm sM} - sA)$.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% Astronomical Images and Astrometry
+\include{psLibSDRS_Astrom}
 \subsection{Astronomical objects}
 …
 \appendix
+\pagebreak
 \section{Configuration File Test Inputs}
 \label{sec:configtest}

trunk/doc/pslib/psLibSDRS_Astrom.tex

-              r3181
+              r3564
 \code{TRUE} if they are all correctly assigned, otherwise \code{FALSE}.
 \subsubsection{Coordinate Transformations}
+\subsubsection{Detector Coordinate Transformations}
 \begin{figure}
 …
 \end{verbatim}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% add psExposure to metadata?
+% add grommit to metadata?
+% add photsystem data to metadata?
+% add statistics to metadata?
+\subsubsection{Observatory data}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\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.
+\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:  transformation to the sky using
+pre-computed coefficients 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{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 (const psPlane *coord, const psFPA *fpa);
+psChip *psChipInFPA (const psPlane *coord, const psFPA *fpa);
+psCell *psCellInChip(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);
+% astrometry comes from cell (no need for parent)
+\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);
+% astrometry comes from chip (no need for parent)
+\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, float color, float mag, const psFPA *fpa);
+% astrometry comes from FPA (no need for parent)
+\end{verbatim}
+which converts coordinates \code{in} on the specified focal plane
+\code{fpa} to tangent plane coordinates, applying the appropriate
+distortion terms.  The \code{color} and magnitude (\code{mag}) of the
+source is necessary in order to perform the distortion between the
+focal plane and the tangent plane.
+\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}.
+% astrometry comes from cell
+\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.
+% astrometry comes from cell,chip,fpa (PARENT IS NEEDED HERE)
+\begin{verbatim}
+psSphere *psCoordCellToSky(psSphere *out, const psPlane *in, float color, float mag, 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.  The \code{color} and magnitude (\code{mag}) of the source is
+necessary in order to perform the distortion between the focal plane
+and the tangent plane.
+% astrometry comes from cell (no need for parent)
+\begin{verbatim}
+psSphere *psCoordCellToSkyQuick(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.toTP}.  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, float color, float mag, const psFPA *fpa);
+\end{verbatim}
+which converts the tangent plane coordinates \code{in} to focal plane
+coordinates.  The \code{color} and magnitude (\code{mag}) of the
+source is necessary in order to perform the distortion between the
+focal plane and the tangent plane.
+\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, float color, float mag, 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.  The \code{color} and magnitude (\code{mag}) of the
+source is necessary in order to perform the distortion between the
+focal plane and the tangent plane.
+\begin{verbatim}
+psPlane *psCoordSkyToCellQuick(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.toTP}.
+Although the accuracy of this transformation is lower than the
+complete transformation above, the calculation is substantially faster
+as it only involves linear transformations.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\subsection{Astrometry and World Coordinate System}
+The FITS World Coordinate System (WCS) headers are commonly employed
+with astronomical images in order to relate pixels to celestial (or
+otherwise) coordinates.  Since it is a FITS standard, we must be able
+to read and write from WCS into our internal format.  For the time
+being, we will consider only celestial WCS (i.e., no spectral
+wavelength calibrations, etc).  Because WCS does not support the
+multiple layers that we have built for \PS{}, we will use a simple
+internal representation: a transformation, which handles any
+distortions (i.e., goes directly from the coordinate frame of the
+image to the tangent plane); and the projection.
+\begin{verbatim}
+bool psAstrometryReadWCS(psPlaneTransform **transform, // Output transformation
+                         psProjection **projection, // Output projection
+                         psMetadata *header // Input FITS header
+                         );
+bool psAstrometryWriteWCS(psMetadata *header, // Output FITS header
+                          psPlaneTransform *transform, // Input transformation
+                          psProjection *projection, // Input projection
+                          double color, // Mean color to use
+                          double magnitude, // Mean magnitude to use
+                          );
+bool psAstrometrySimplify(psPlaneTransform **transform, // Output transformation
+                          psProjection **projection, // Output projection
+                          psCell *cell // Cell for which to generate transform and projection
+                          );
+\end{verbatim}
+\code{pmReadAstrometry} shall parse the specified FITS \code{header},
+returning new instances of the \code{transform} and \code{projection}
+that represent the WCS.  The function shall return \code{true} if it
+was able to successfully generate the outputs; otherwise it shall
+return \code{false}.
+\code{pmWriteAstrometry} shall add WCS keywords to the supplied FITS
+\code{header} that implement the given \code{transform} and
+\code{projection}.  The function shall return \code{true} if it was
+able to successfully generate the output; otherwise it shall return
+\code{false}.
+\code{pmSimplifyAstrometry} shall take a \code{cell} and simplify the
+internal astrometric representation (\code{cell->toFPA} or equivalent,
+\code{cell->parent->parent->toTangentPlane} and
+\code{cell->parent->parent->grommit}) to a single \code{transform} and
+\code{projection}.  This allows the subsequent use of
+\code{pmWriteAstrometry} in the case that we have only the
+multi-layered \PS{} internal astrometric representation.  The function
+shall return \code{true} if it was able to successfully generate the
+output; otherwise it shall return \code{false}.
+\subsection{Observatory data}
 We need a container for the observatory data that doesn't change per
 …
 \end{verbatim}
 \subsubsection{Exposure information}
+\subsection{Exposure information}
 We need several quantities from the telescope in order to make a
 …
 \end{verbatim}
-\subsubsection{Environmental Information}
-A-priori astrometric transformations between the tangent plane and the
-sky require several pieces of information describing the current
-environmental conditions.  These quantities are consistent from image
-to image, and may vary only slowly with time.  Pre-computing these
-quantities for exposures means that subsequent transformations are
-faster.  The structure below carries the environment data of interest.
-For historical reasons, this structure is known colloquially as ``the
-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, \code{exp}:
-\begin{verbatim}
-psGrommit *psGrommitAlloc(const psExposure *exp);
-\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}
-The constructor for \code{psFixedPattern} shall be:
-\begin{verbatim}
-psFixedPattern *psFixedPatternAlloc(double x0,        double y0,
-                                    double xScale,    double yScale,
-                                    const psImage *x, const psImage *y);
-\end{verbatim}
-Here, \code{x0}, \code{y0}, \code{xScale} and \code{yScale} have the
-same meaning as in the \code{psFixedPattern} structure.  Note that the
-values of the fixed pattern offsets are specified as images, the
-values from which need to be copied into the \code{double **x} and
-\code{double **y} of \code{psFixedPattern}, and that the number of
-elements may be derived from the size of the images.
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\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.  \tbd{SLAlib functions will be replaced in the next
-release}.
-\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:  transformation to the sky using
-pre-computed coefficients 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{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 (const psPlane *coord, const psFPA *fpa);
-psChip *psChipInFPA (const psPlane *coord, const psFPA *fpa);
-psCell *psCellInChip(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);
-% astrometry comes from cell (no need for parent)
-\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);
-% astrometry comes from chip (no need for parent)
-\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, float color, float mag, const psFPA *fpa);
-% astrometry comes from FPA (no need for parent)
-\end{verbatim}
-which converts coordinates \code{in} on the specified focal plane
-\code{fpa} to tangent plane coordinates, applying the appropriate
-distortion terms.  The \code{color} and magnitude (\code{mag}) of the
-source is necessary in order to perform the distortion between the
-focal plane and the tangent plane.
-\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);
-% astrometry comes from 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, float color, float mag, const psCell *cell);
-% astrometry comes from cell,chip,fpa (PARENT IS NEEDED HERE)
-\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.  The \code{color} and magnitude (\code{mag}) of the source is
-necessary in order to perform the distortion between the focal plane
-and the tangent plane.
-\begin{verbatim}
-psSphere *psCoordCellToSkyQuick(psSphere *out, const psPlane *in, const psCell *cell);
-% astrometry comes from cell (no need for parent)
-\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.toTP}.  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, float color, float mag, const psFPA *fpa);
-\end{verbatim}
-which converts the tangent plane coordinates \code{in} to focal plane
-coordinates.  The \code{color} and magnitude (\code{mag}) of the
-source is necessary in order to perform the distortion between the
-focal plane and the tangent plane.
-\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, float color, float mag, 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.  The \code{color} and magnitude (\code{mag}) of the
-source is necessary in order to perform the distortion between the
-focal plane and the tangent plane.
-\begin{verbatim}
-psPlane *psCoordSkyToCellQuick(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.toTP}.
-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, psTime *lst, float height);
-\end{verbatim}
-which returns the airmass for a given position and local sidereal time
-(\code{lst}).
-\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{Astrometry and World Coordinate System}
-The FITS World Coordinate System (WCS) headers are commonly employed
-with astronomical images in order to relate pixels to celestial (or
-otherwise) coordinates.  Since it is a FITS standard, we must be able
-to read and write from WCS into our internal format.  For the time
-being, we will consider only celestial WCS (i.e., no spectral
-wavelength calibrations, etc).  Because WCS does not support the
-multiple layers that we have built for \PS{}, we will use a simple
-internal representation: a transformation, which handles any
-distortions (i.e., goes directly from the coordinate frame of the
-image to the tangent plane); and the projection.
-\begin{verbatim}
-bool psAstrometryReadWCS(psPlaneTransform **transform, // Output transformation
-                         psProjection **projection, // Output projection
-                         psMetadata *header // Input FITS header
-                         );
-bool psAstrometryWriteWCS(psMetadata *header, // Output FITS header
-                          psPlaneTransform *transform, // Input transformation
-                          psProjection *projection, // Input projection
-                          double color, // Mean color to use
-                          double magnitude, // Mean magnitude to use
-                          );
-bool psAstrometrySimplify(psPlaneTransform **transform, // Output transformation
-                          psProjection **projection, // Output projection
-                          psCell *cell // Cell for which to generate transform and projection
-                          );
-\end{verbatim}
-\code{pmReadAstrometry} shall parse the specified FITS \code{header},
-returning new instances of the \code{transform} and \code{projection}
-that represent the WCS.  The function shall return \code{true} if it
-was able to successfully generate the outputs; otherwise it shall
-return \code{false}.
-\code{pmWriteAstrometry} shall add WCS keywords to the supplied FITS
-\code{header} that implement the given \code{transform} and
-\code{projection}.  The function shall return \code{true} if it was
-able to successfully generate the output; otherwise it shall return
-\code{false}.
-\code{pmSimplifyAstrometry} shall take a \code{cell} and simplify the
-internal astrometric representation (\code{cell->toFPA} or equivalent,
-\code{cell->parent->parent->toTangentPlane} and
-\code{cell->parent->parent->grommit}) to a single \code{transform} and
-\code{projection}.  This allows the subsequent use of
-\code{pmWriteAstrometry} in the case that we have only the
-multi-layered \PS{} internal astrometric representation.  The function
-shall return \code{true} if it was able to successfully generate the
-output; otherwise it shall return \code{false}.

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