Context Navigation

← Previous Change
Next Change →

psLibADD.tex

Timestamp:

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

Author:

eugene

Message:

serious reorg for release 13

File:

: 1 edited

trunk/doc/pslib/psLibADD.tex (modified) (63 diffs)

Legend:

: Unmodified
: Added
: Removed

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}

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

Context Navigation

Changeset 3564 for trunk/doc/pslib/psLibADD.tex

Legend:

trunk/doc/pslib/psLibADD.tex

Download in other formats: