%%% $Id: psLibADD.tex,v 1.95 2006-12-09 00:30:57 eugene Exp $
\documentclass[panstarrs]{panstarrs}

\usepackage{amsmath}
\usepackage{verbatim}

% basic document variables
\title{Pan-STARRS PS-1 Image Processing Pipeline}
\subtitle{Algorithm Design Description}
\shorttitle{IPP ADD}
\author{Eugene Magnier, Paul Price, Joshua Hoblitt, Robert Lupton}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{12}
\docnumber{PSDC-430-006}

% \setcounter{tocdepth}{5} % lowest level to be included in toc

\newcommand\citealt{}

\begin{document}
\maketitle

% -- Revision History --
% provide explicit values for the old versions
% use '\theversion' for the current version (set above)
% use \hline between each table row
\RevisionsStart
% version     Date         Description
00 & 2004 Mar 11 & Hacking \\ \hline
01 & 2004 May 21 & Added section on 2D Chebyshev fitting, then removed. \\ \hline
02 & 2004 Jun 22 & modified stats specification \\ \hline
03 & 2004 Jul 13 & \\ \hline
04 & 2004 Aug 16 & \\ \hline
05 & 2004 Sep 01 & \\ \hline
06 & 2004 Sep 07 & Frozen for PSLib-2 \\ \hline
07 & 2004 Nov 24 & Frozen for Cycle 4 \\ \hline
08 & 2005 Jan 21 & Draft for Cycle 5 \\ \hline
09 & 2005 Feb 14 & Frozen for Cycle 5 \\ \hline
10 & 2005 Apr 19 & Frozen for Cycle 6 \\ \hline
11 & 2005 Apr 27 & Update for Cycle 6 \\ \hline
\RevisionsEnd

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\DocumentsInternalSection
PSDC-230-001  &   PS-1 Design Reference Mission \\ \hline
PSDC-430-004  &   Pan-STARRS PS-1 IPP C Code Conventions \\ \hline
PSDC-430-005  &   Pan-STARRS PS-1 IPP Software Requirements Specification \\ \hline
PSDC-430-006  &   Pan-STARRS PS-1 IPP Algorithm Design Document \\ \hline
PSDC-430-011  &   Pan-STARRS PS-1 IPP System/Subsystem Design Description \\ \hline
\DocumentsExternalSection
Posix Standard                      & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\ \hline
Numerical Recipes (NR)              & Press, Teukolsky, Vetterline, Flannery \\ \hline
Knuth, D.E.                         & Sorting and Searching; The Art of Computer Programming \\ \hline
Sedgewick, R.                       & Algorithms, Ch. 8 \\ \hline
Sorting Summary                     & \code{http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm } \\ \hline
GSL                                 & \\ \hline
FFTW                                & (Fastest Fourier Transform in the West) \\ 
                                    & \code{http://www.fftw.org} \\ \hline
FITS Projection Article             & {Greisen \& Calabretta (1995, ADASS, 4, 233)} \\
                                    & \code{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps} \\ \hline
Hipparcos and Tycho Catalogues      & \code{http://astro.estec.esa.nl/Hipparcos/CATALOGUE_VOL1/catalog_vol1.html} \\ \hline
Zombeck                             & ``Handbook of Space Astronomy and Astrophysics'', second edition, \\ 
                                    & \code{http://ads.harvard.edu/books/hsaa/toc.html} \\ \hline
Reingold \& Dershowitz              & ``Calendrical Calculations: The Millenium Edition'', Cambridge University Press, 2002. \\
\hline
\DocumentsEnd

\tableofcontents
\pagebreak 
\pagenumbering{arabic}

% \section{Pan-STARRS Library PSLib}

\section{PSLib Math Utilities}

\subsection{Sorting}

A variety of sorting algorithms exist, with a wide range in speed for
different set sizes.  Three of the best sorting algorithms are
``heapsort'', ``quicksort'', and ``mergesort'' (see, e.g.,
\citealt{knuth}, \citealt{sedgewick}, \citealt{press}).  These three
sorting algorithms all run in a time of $O(n \log n)$ on the average,
but have different worse-case run times.  Implementations of all three
sorting algorithms are described in references above and in various
places online (e.g.,
http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm).
The linux \code{qsort} function uses a heapsort if it can allocate a
sufficiently large temporary buffer, and otherwise resorts to a
quicksort in-place.  The GSL function \code{gsl_heapsort} uses the
heapsort algorithm.  Both of these implemenations require a pointer to
the comparison function, which may result in a slower code than if the
comparison were done within the sort function.

For PSLib, the sorting functions shall be implemented via the system
\code{qsort}.  The function \code{psSort} shall return the sorted
result of the input array without over-writing the input array.  The
function \code{psSortIndex} shall return an integer index to the
sequence of the input array without overwriting the input array.
Given the following line of code:
\begin{verbatim}
out = psSortIndex (NULL,&in);}
\end{verbatim}
the elements of the array \code{out} are in the sequence
\code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}.

\subsection{Smoothing: Boxcar and Gaussian}
\label{smooth}

Smoothing may occasionally be perfomed on data.  We present the
algorithms for two typical versions: boxcar and Gaussian smoothing.
In both smoothing techniques, given a series of data values $f_i(x_i)$
where $x_i$ are the values of the corresponding to the center of the
bin, the smoothed values $g_i(x_i)$ are determined by calculating a
linear combination based on the input data point and its nearest $2N$
neighbors in the form:

\begin{equation}
g_i = \sum_{j=i_{\rm min}}^{i_{\rm max}} c_{ij} f_j
\end{equation}
%
where the values of $c_{ij}$ determine the filter type.  For boxcar
smoothing, the values $c_{ij}$ are constant and scaled to maintain the
zeroth moment of the data (care must be taken at the ends of the data
range to reduce the value of $c_{ij}$ as fewer input data points may
be used).  For Gaussian smoothing, the crucial parameter is
$\sigma_x$, the standard deviation.  The values of $i_{\rm min}$ and
$i_{\rm max}$ are functions of the standard deviation: $i_{\rm min}$
corresponds to the bin in which $x_i - N\sigma_x$ is found; similarly
$i_{\rm max}$ is the bin corresponding to $x_i + N\sigma_x$.  The
value of $N$ should be chosen to be large enough to sample the
Gaussian, $N = 5$.  The values of $c_{ij}$ are then just the Gaussian
curve:

\begin{equation}
c_{ij} = \frac{e^{\frac{-(x_j - x_i)^2}{2\sigma_x^2}}}{\sqrt{2\pi\sigma_x^2}}
\end{equation}

\subsection{Statistics}

The general statistics function \code{psStats} performs a variety of
statistical calculations on a collection of floating point numbers.
These statistics include both sample values, in which the formal
statistic is calculated for the complete sample, and robust values, in
which the statistic is estimated on the basis of an assumption of an
underlying population.  While more reliable in general, the robust
estimators may be somewhat slower than the sample statisic, so their
use may vary depending on the context.  In addition, certain
sigma-clipped values are defined as an intermediate choice between the
sample and robust estimators.

\subsubsection{Sample Statistics}

We define the following statistical terms, assuming there is a set of
data elements $x_i$ with (standard) errors $\sigma_i$.

\paragraph{Mean}

The simple mean is defined as:
\begin{equation}
\bar{x} = \frac {1}{N} \sum^N_i x_i
\end{equation}

\paragraph{Weighted Mean}

The weighted mean is defined as:
\begin{equation}
\bar{x} = \sum_i \frac{x_i}{\sigma_i^2} \ / \ \sum_i \frac{1}{\sigma_i^2}
\end{equation}

In the event that all the errors are identical, this reduces to the
standard definition of the mean.

\paragraph{Median}

The median is defined as the value for which 50\% of the data values
are larger and 50\% are smaller.  For a sample, the values are sorted
and the median is given by the value of the middle element, if the
number of values is odd, and by the average of the two middle values
if the number of values is even.  This median should be avoided for
samples which are large (e.g., $N > 10^4$ elements) as the basic
robust median is quicker and more accurate.  Errors are ignored when
calculating the sample median.

\paragraph{Upper and Lower Quartiles}

The upper and lower quartiles ($U_{\frac{1}{4}}$ and
$L_{\frac{1}{4}}$) are first and last quarter equivalents to the
median.  The upper quartile is the value for which 25\% of the data
are larger and 75\% are smaller.  The lower quartile is the value for
which 75\% of the data are larger and 25\% are smaller.  For a sample,
the values are sorted and the quartiles are given by the values of
elements $N/4$ and $3N/4$.  For the purpose of a sample upper and
lower quartile, it is sufficient to provide the closest integer entry
to these values.  The sample quartiles should be avoided for samples
which are large (e.g., $N > 10^4$ elements) as the robust quartiles
are quicker and more accurate.  Errors are ignored when calculating
the sample quartiles.

\paragraph{Standard Deviation}

The standard deviation of the sample is given by:

\begin{equation}
\sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}}
\end{equation}

To minimize the numerical rounding error, this should be calculated
numerically as:

\begin{equation}
\sigma = \sqrt{\frac{1}{N - 1} \left[ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} \left(\sum_{i = 1}^{N} (x_i - \bar{x})\right)^2 \ \right]}
\end{equation}

If the errors are known, then the sample standard deviation is:

\begin{equation}
\sigma = \left( \sum_i \frac{1}{\sigma_i^2} \right) ^{-1/2}
\end{equation}


\subsubsection{Clipped Statistics}

The clipped statistics are used to determine the mean and standard
deviation of a distribution in the presence of outliers. The clipped
statistics are quicker than the complete robust statistical estimators
and more reliable than the sample statistics.  The clipped statistics
algorithm produces the clipped mean and clipped standard deviation.
The algorithm has 2 parameters: $N$, the number of iterations and $k$,
the multiplying factor of the standard deviation used to exclude
outliers.  Typical values for both $N$ and $k$ are 3.  The algorithm
is as follows:

\begin{enumerate}
\item Compute the sample median. The number of data points must be
      limited to 10000; the input dataset must be randomly subsampled
      if more data points are used.
\item Compute the sample standard deviation.
\item Use the sample median as the first estimator of the mean, $\bar{x}$.
\item Use the sample standard deviation as the first estimator of the
  true standard deviation, $\sigma$.
\item Repeat the following N times:
  \begin{enumerate}
  \item Exclude all values $x_{i}$ for which $|x_{i} - \bar{x}| > k \sigma$.
  \item Compute the mean and standard deviation of the sub-sample.
  \item Use the new mean for $\bar{x}$.
  \item Use the new standard deviations for $\sigma$.
\end{enumerate}
\item The last calculated value of $\bar{x}$ is the clipped mean.
\item The last calculated value of $\sigma$ is the clipped standard
  deviation.
\end{enumerate}

If the errors in the input values are known, then the clips are made
on the basis of the errors in the input values instead of the standard
deviation of the sample: values are excluded for which $|x_i -
\bar{x}| > k \sigma_i$.

\subsubsection{Robust Statistics}

Robust statistics algorithms provide estimators of basic statistical
concepts which are reliable even for data samples with significant
contamination.  A typical case is the situation in which the data of
interest represent a primary population of interest with a
single-valued mean and standard deviation and a secondary population
of data with a substantially different distribution.  For example, an
image of an uncrowded night-time field may consist of a sparse
collection of stars and an overall background level.  The majority of
pixels have the background value, with some variance due to noise
sources, but many are significantly higher (contributed by stars) or
significantly lower (dead pixels).  If we want to measure the mean of
the background, a robust mean is necessary as the outliers will bias
the sample statistics.  We define two levels of robust statistical
estimators: robust statistics using the cumulative histogram, and
statistics measured by fitting the differential histogram.

The robust statistics from the cumulative are calculated by
constructing a cumulative histogram of the values and performing
simple measurements of the histogram data distribution.  The use of
the cumulative histogram by this algorithm reduces the sensitivity to
the exact bin size.  The initial bin size is set to 1/1000 of the
total data range. The data values are found using quadratic
interpolation between the bin of interest and its two neighbors.

\begin{itemize}
\item Construct the histogram with the specified bin size.
\item Construct the cumulative histogram from the specific histogram
\item Find the bin which contains the 50\% data point.  
\item Interpolate to the exact 50\% position: this is the robust
  histogram median.
\item Find the bins which contains the 15.8655\% and 84.1345\% data
  points.  
\item Interpolate to find these two positions exactly: these are the
  $\pm 1\sigma$ positions.
\item Determine $\sigma$ as 1/2 of the distance between these
  positions.
\item If the measured $\sigma$ is less than 2 times the bin size,
  exclude points which are more than 25 bins from the median,
  recalculate the bin size, and perform the algorithm again.
\item Find the bins which contains the 25\% and 75\% data points.
\item Interpolate to find these two positions exactly: these are the
  upper and lower quartile positions.
\end{itemize}

If the errors in the input values are known, then the same approach is
used, except that the histograms become probability density functions
(PDFs).  In this case, the input values are spread out, so that they
do not simply contribute a single unit to the histogram, but rather
contribute a fraction of a value, equivalent to the weight.  In the
interests of speed, a boxcar PDF may be used to represent each input
value (as opposed to a Gaussian), where the boxcar width is equal to
$2 \sqrt{2 \ln 2}$ times the error and each input value contributes
constant area.  Then the robust median and standard deviation are
estimated in the same manner as above.

\subsubsection{Fitted Statistics}

The fitted statistics algorithm starts with the histogram used for the
robust statistics and determines the population statistics by fitting
a Gaussian model to the histogram.  The algorithm is as follows:

\begin{itemize}
\item Perform the Robust Histogram Statistics algorithm above,
  yielding an estimated standard deviation, $\sigma$.

\item Generate a new histogram for the data sample setting a bin size,
  $d\sigma$, based on the estimated standard deviation and the number
  of data points in the inner 50 percentile ($N_{\rm 50}$) as follows:

\begin{itemize}
\item let $dN = (\sigma / d\sigma) = 0.017 N_{50}$
\item limit $dN$ to the range 1 to 5.
\item set the bin size $d\sigma = \sigma / dN$
\end{itemize}

\item Smooth the resulting histogram with a Gaussian with $\sigma_x$ =
  1 bin in this new histogram.
\item Find the bin with the peak value in the range $\pm 2 \sigma$ of
  the robust histogram median.

\item Fit a Gaussian to the bins in the range $\pm 20 \sigma$ of the
  robust histogram median. Limit the fit range to the data range, if
  the latter is less then $\pm 20 \sigma$.  If the data range is small
  compared to the estimated $\sigma$, fit at least 4 bins of the
  hisgram centered on the robust histogram median.

\item The robust mean $\mbox{mean}_r$ is derived directly from the
  fitted Gaussian mean.  
\item The robust standard deviation, $\sigma_r$, is determined by
  subtracting the smoothing scale in quadrature: $\sigma_r^2 =
  \sigma_{\rm fit}^2 - \sigma_s^2$
\end{itemize}

To explain the choice of the histogram bin size: a histogram of a
Gaussian distribution with bin size $d\sigma$ will have approximately
$(2.35 \sigma/d\sigma)$ bins covering the range LQ to UQ.  Thus, the
average number of points per bin ($N_{\rm bin}$) in that interval will
be $N_{50} / (2.35 \sigma/d\sigma)$.  The value of $d\sigma$ should be
no larger than $\sigma$, regardless of the number of points, to avoid
too much undersampling.  The value of $d\sigma$ should also be no
smaller than $5\sigma$, again regardless of the number of points, to
avoid excessive oversample.  Intermediate to those two values, the bin
size is choosen to keep about 25 points per bin.  Thus, the bin size
($d\sigma$) is set to about:
\begin{equation}
d\sigma = 2.35 \sigma (N_{\rm bin}/N_{50}) = (25 \times 2.35) (\sigma/N_{50})
\end{equation}
With the limitation that $\sigma/d\sigma$ should be limited on one end
to the value 1, and the other to the value 5.  The easiest way to set
this limit is to define dN to be:
\begin{equation}
dN = (\sigma / d\sigma) = (N_{50} / N_{\rm bin}) / 2.35 = 0.017 * N_{50}
\end{equation}
\subsubsection{Histograms}

When calculating histograms in the presence of known errors in the
input values, the approach described above for the robust statistics
is used (i.e., the histograms become probability density functions).

An example may help here.  Say we have our histogram bounds being 0,
1, 2, 3, 4, 5; and our value is $2.5 \pm 0.5$.  Then, the width of the
contribution is $0.5 \times 2.35... \approx 1.175$.  Half the width is
0.5875, so we will treat this value as a boxcar from $2.5 - 0.5875$ to
$2.5 + 0.5875$.

Consequently, the bins 0 to 1 and 4 to 5 get no value, because none of
the boxcar overlaps.  The bin 1 to 2 gets 0.0875, because that's the
fraction of the boxcar that overlaps with it; same thing with the bin
3 to 4.  The bin 2 to 3 gets 0.825 because that's the fraction of the
boxcar that overlaps with it.  So the single value $2.5 \pm 0.5$ makes
the following histogram:

\begin{tabular}{lr}
Bin & Value \\ \hline
0--1 & 0 \\
1--2 & 0.0875 \\
2--3 & 0.8250 \\
3--4 & 0.0875 \\
4--5 & 0 \\
\end{tabular}

Note that the total adds to one --- the number of values added.

\subsection{Polynomials}
\label{sec:polynomials}

We will employ Chebyshev polynomials (NR \S 5.8) to approximate functions:
\begin{equation}
f(x) = \sum_{i=0}^{n} c_i T_i(x)
\end{equation}
These have some desirable features:
\begin{itemize}
\item They are bounded on $-1 < x < 1$, with the maxima and minima
over this range being 1 and $-1$, respectively;
\item Truncation of the higher-order terms leaves one with the most accurate
lower-order polynomial representation of the desired function.
\end{itemize}

The first few Chebyshev polynomials are:
\begin{eqnarray}
T_0(x) & = & 1 \\
T_1(x) & = & x \\
T_2(x) & = & 2x^2 - 1 \\
T_3(x) & = & 4x^3 - 3x \\
T_4(x) & = & 8x^4 - 8x^2 + 1 
\end{eqnarray}
Chebyshev polynomials follow the recurrence relation:
\begin{equation}
T_{n+1} = 2xT_n - T_{n-1}
\end{equation}

Practically, Chebyshev polynomials should be evaluated using Clenshaw's recurrence
formula (NR \S 5.5):
\begin{eqnarray}
d_j  & = & 2xd_{j+1} - d_{j+2} + c_j \\
f(x) & = & x*d_1 - d_2 + 1/2 c_0
\end{eqnarray}

It shall be the responsibility of the user to convert the domain into the range
$-1 < x < 1$.

\subsubsection{Multi-dimensional polynomials}

Multi-dimensional polynomials shall be composed of multiplications of
1D Chebyshev polynomials, with the coefficients stored in tensors of
the appropriate rank.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Fitting}

\subsubsection{Chi-squared}
\label{chisq}

Given a set of N ordinates, $x_i$, measured coordinates, $y_i$, with
errors, $\sigma_i$, and a model function, $f(x_i)$, then $\chi^2$
(``chi-squared'') is defined:

\begin{equation}
\chi^2 = \sum_{i=0}^{N-1} \left( \frac{y_i - f(x_i)}{\sigma_i} \right)^2
\end{equation}

\subsubsection{General Polynomial Fitting}

Given a set of data values $y_i$ with errors $\sigma_i$, related to
independent data values $x_i$, we would like to fit a polynomial model
of $N_{par}$ free parameters of the form $y = \sum_{j=0}^{N_{par}}
\alpha_j x^j$.  The model is determined by minimizing the value of
$\chi^2$ (\ref{chisq}), and the minimization is determined by solving
for the set of parameters $\alpha_j$ for which the partial derivatives
of the $\chi^2$ with respect to those parameters are zero.  The
partial derivatives are

\begin{equation}
\frac{\partial \chi^2}{\partial \alpha_k} = -2 \sum_i (y_i - \sum_j x_i^j \alpha_j) \frac{x_i^k}{\sigma_i^2}
\end{equation}

Setting the derivatives to zero, this may be reduced to the following
matrix equation:

\begin{equation}
\sum_j \alpha_j \sum_i \frac{x_i^k x_i^j}{\sigma_i^2} = \sum_i \frac{x_i^k y_i}{\sigma_i^2}
\end{equation}

This matrix equation may be solved with LU Decomposition
(section~\ref{LUdecomp}).

\subsection{Non-linear Minimization}

Non-linear minimization techniques use an iterative approach to find a
minimization when an analytical inversion is impractical or not
possible.  These techniques use a starting guess for the parameters of
interest, and make a sequence of new guess parameters based on the
properties of the function at the previous position.  If the new
parameters yield a reduced function value, the new position is used as
the starting position for the next iteration.  Otherwise, the guess
must be modified and another attempt is made.  Convergence may be
determined based on the absolute amount of change in the function
value, or by comparison with the expectation for a linear system.

The two common techniques used to construct a guess parameter set are
the 'steepest descent method' and the 'Gauss-Newton method'.  In the
first case, the guess is selected some distance along the local
gradient.  In the second case, a local Taylor expansion of the
function is used to construct a linear model for the function, and the
new guess is chosen to minimize that linear model.  The methods
discussed below make use of combinations of these two methods.  Aside
from the differences in their guess steps, the two methods differ in
using first derivatives of the function supplied by the user, or by
locally calculating the first derivatives.

Mathematically, we would like to choose the parameter set $a_m$ to
minimize a function of those parameters $F(a_m)$.  We iterate by
choosing a new parameter set $a^\prime_m = a_m + \delta_m$ based on
the behavior of the function at $a_m$.  

The steepest descent method chooses a step direction of $\bar{\delta}
= -\nabla F$, or defining $g_m$ as a component of the gradient,
$\delta_m = -g_m$.  The Gauss-Newton method uses a Taylor expansion of
the function to solve for the step: $(\nabla^2 F) \bar{\delta} =
-\nabla F$.  Defining an element of the Hessian matrix $H_{m,n}$ as a
component of the second derivatives, we can write the Gauss-Newton
step as $\delta_m = - H_{m,n}^{-1} g_n$.

\subsubsection{Levenberg-Marquardt Method}

In the Levenberg-Marquardt Method (LMM; see NR \S 15.5, Madsen et al),
we make a guess at the input parameters, evaluate the function of
interest, vary the parameters by a particular choice based on the
gradient, evaluate the function again, and adjust the parameters and
the parameter varient based on the results.  The LMM requires the
second derivative of the function to be negligible, as in the case of
minimizing $\chi^2$.

Consider the chi-square function. Given some ordinates, $x_i$, we
would like to find the parameters, $a_m$, of the function which
minimize $\chi^2$ for some measurements, $y_i$ and associated errors,
$\sigma_i$:
\begin{eqnarray}
\chi^2 (\bar{a})  & = & \sum_i \frac{1}{\sigma_i^2} \left( y_i - p(x_i;a_m) \right)^2 
\end{eqnarray}

We simplify this as:
\begin{eqnarray}
p_i (a_m)         & = & p(x_i;a_m) \\
f_i (a_m)         & = & \frac{1}{\sigma_i} (y_i - p_i) \\
\chi^2 (\bar{a})  & = & \sum_i f_i^2 
\end{eqnarray}

We write the minimization function $F(a_m) = \frac{1}{2} \chi^2$ to
avoid the various extra factors of 2.  We can now write out the needed
derivatives in terms of $f_i$:
\begin{eqnarray}
F(a_m)            & = & \frac{1}{2} \sum_i f_i^2 \\
\nabla F(a_m)     & = & \sum_i f_i \frac{\partial f_i}{\partial a_m} \\
\nabla^2 F(a_m)   & = & \sum_i \frac{\partial f_i}{\partial a_m} \frac{\partial f_i}{\partial a_n}
\end{eqnarray}

where we have dropped the second-derivatives of the function in the
representation of $\nabla^2 F(a_m)$.  Since $\frac{\partial
f_i}{\partial a_m} = -\frac{1}{\sigma_i}\frac{\partial p_i}{\partial
a_m}$, we can write these in terms of the derivatives of $p_i$ only:
\begin{eqnarray}
\nabla F(a_m)     & = & -\sum_i \frac{f_i}{\sigma_i} \frac{\partial p_i}{\partial a_m} \\
\nabla^2 F(a_m)   & = &  \sum_i \frac{1}{\sigma_i^2} \frac{\partial p_i}{\partial a_m} \frac{\partial p_i}{\partial a_n}
\end{eqnarray}

Writing these in matrix representation, and replacing $f_i$, we have:
\begin{eqnarray}
-g_m              & = & \sum_i \frac{(y_i - p_i)}{\sigma_i^2} \frac{\partial p_i}{\partial a_m} \\
H_{m,n}           & = & \sum_i \frac{1}{\sigma_i^2} \frac{\partial p_i}{\partial a_m} \frac{\partial p_i}{\partial a_n}
\end{eqnarray}

In the Levenberg-Marquart Method, we define a new guess using a
combination of the Steepest Descent and Gauss-Newton methods discussed
above.  We replace the Hessian matrix above with $A_{m,n}$ as a
variant on $H_{m,n}$ as follows:

\begin{eqnarray}
A_{m,n} & = & H_{m,n} \mbox{if} (j \ne k) \\
A_{m,n} & = & H_{m,n}(1 + \lambda) \mbox{if} (j = k)
\end{eqnarray}
%
and solve the system of equations represented by:
\begin{equation}
A_{m,n} \delta_n = -g_m
\end{equation}
%
The new parameter guess is then found from this value with
$\alpha^\prime_n = \alpha_n + \delta_n$. We use this parameter set to
evaluate the function.

To assess the quality of the new parameter set, we compare the change
in $\chi^2$ with the change expected from the linear model (the Taylor
expansion).  If the linear model were correct, we would have expected
a change (a reduction) in $\chi^2$ of $\Delta =
\frac{\lambda}{2}\sum\delta_m^2 + \frac{1}{2}\sum\delta_m g_m$.  We
use the 'gain ratio' $\rho = \frac{\chi^2_{\rm old} - \chi^2_{\rm
new}}{\Delta}$ to judge the new step. If $\rho > 0$, we accept this
new set of parameters and decrease $\lambda$ by a factor of 10,
otherwise we keep the old set, and increase the value of $\lambda$ by
a factor of 10.  We repeat this process until the value of the
function changes by much less than the tolerance.  The resulting
values of $a_m$ are the best-fit parameters for the system.

The covariance matrix, $C_{i,j}$, which is the inverse of the matrix
$H_{m,n}$ provides an estimate of the confidence limits of the
parameters.

%If the errors are normally distributed, the formal errors on the
%parameters are then calculated by setting $\lambda = 0$ and
%calculating the covarience matrix $C_{i,j}$, the inverse of the matrix
%$\alpha_{j,k}$.
%The independent 68.3\% confidence limit on parameter $a_k$ is then
%$\sqrt{C_{k,k}}$.  Confidence contours for sets of parameters may be
%defined as well by the function $\Delta = \delta\bar{a} P_{j,k}^{-1}
%\delta\bar{a}$ where $P_{j,k}$ is the projected matrix of $C_{j,k}$,
%ie those rows and columns of $C_{j,k}$ associated with the parameters
%of interest, the vector $\delta\bar{a}$.  The value of $\Delta$ is
%given by the table below for specific confidence limits and numbers of
%parameters.  Note that it is necessary to be able to calculate both
%the function as well as its derivative for any combination of
%parameters and dependent variables.
%
%\begin{center}
%\begin{tabular}{|l|r|r|r|}
%\hline
%{\bf P} & \multicolumn{3}{c|}{\bf $N_{par}$} \\
%        & 1    & 2    & 3    \\
%\hline
%68.3\%  & 1.00 & 2.30 & 3.53 \\
%95.4\%  & 4.00 & 6.17 & 8.02 \\
%99.73\% & 9.00 & 11.8 & 14.2 \\
%\hline
%\end{tabular}
%\end{center}

\subsubsection{Powell's method}

Powell's method is a type of ``Direction Set'' methods in
multi-dimensions for finding a local minimum.  Given a starting point
(the ``best guess'' for the minimum) and a set of direction vectors, a
direction set method advances in the direction of the vectors,
determines a new direction vector by some method, and proceeds in this
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 \S 10.5.

We will use for our algorithm the modified version of Powell's
Quadratically Convergent Method, which is described below, adapted
from NR.

\begin{enumerate}
\item Given a function in $N$ dimensions to minimize, $f$, and a best
  guess for the minimum, point $P$ in $N$ dimensions, take an initial
  set of $N$ vectors, $v_i$, to be the unit vectors.
\item Set point $Q = P$.
\item For each dimension in turn, move $Q$ \textit{only} in the
  direction $v_i$ to minimize the function of interest.
\item Set vector $u = Q - P$.
\item Move $Q$ \textit{only} in the direction $u$
\item Replace the vector along which the largest minimization was
  made, $v_{i,\rm max}$, with $u$, except under either of the
  following circumstances:
  \begin{itemize}
  \item If $f_QP \ge f_P$, then there is no point in keeping the new
    vector, because there is no further minimization to be made in
    that direction.
  \item If $2 ( f_P - 2f_Q + f_{QP} ) \left[ ( f_P - f_Q ) -
  \Delta_{\rm max} \right]^2 \ge ( f_P - f_{QP} )^2 \Delta_{\rm max}$,
  then either the decrease in the function was not due to any single
  direction, or we are close to the minimum.
  \end{itemize}
  where $f_P = f(P)$, $f_Q = f(Q)$, $f_{QP} = f(2Q - P)$, and
  $\Delta_{\rm max} \ge 0$ is the magnitude of the minimization made
  along $v_{i,\rm max}$.
\item Set $P$ to $Q$.
\item Return to step 3 until the change in this last move is less
  than some specified tolerance, or a maximum number of iterations
  has been reached.
\end{enumerate}

In regards to minimizing the function only in a particular direction,
we shall adopt, as NR recommends, bracketing the minimum before
applying Brent's method, \tbd{which will be specified in detail
later}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\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:
\begin{equation} 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 \end{equation}
This expression is more efficiently evaluated by factoring and
calculating the expresion as:
\begin{equation} r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \end{equation}
\begin{equation} V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \end{equation}

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
2$).  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:
\begin{equation} x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \end{equation}
\begin{equation} y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \end{equation}
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$).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\pagebreak 

\subsection{Matrix Operations}

In this section, we define the linear algebra operations performed on
matrices.  We have defined APIs to implement the following matrix
functions:

\begin{itemize}
\item Invert a matrix;
\item Calculate a matrix determinant;
\item Perform matrix addition, subtraction and multiplication;
\item Transpose a matrix; and
\item Convert a matrix to a vector.
\end{itemize}

Many of these operations are implemented using LU decomposition.  We
define LU decomposition in the following paragraph.  Implementation of
LU decomposition shall make use of the GSL function
\code{gsl_linalg_LU_decomp}.

\subsubsection{LU Decomposition}
\label{LUdecomp}

We wish to decompose the matrix $A$ with elements $a_{ij}$ into
diagonal matrices that satisfy the relationship $A = L U$ where $L$ is
a lower-diagonal matrix of the form:
\begin{equation}
L = \left(
\begin{matrix}
\alpha_{11} & 0           & 0           & 0 \\
\alpha_{21} & \alpha_{22} & 0           & 0 \\
\alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\
\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\
\end{matrix} \right)
\end{equation}
%
(where all diagonal values $\alpha_{ii} = 1$) and $U$ is an upper-diagonal matrix of the form:
\begin{equation}
U = \left(
\begin{matrix}
\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\
0          & \beta_{22} & \beta_{23} & \beta_{24} \\
0          & 0          & \beta_{33} & \beta_{34} \\
0          & 0          & 0          & \beta_{44} \\
\end{matrix} \right)
\end{equation}

We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following
this procedure.  First, $\alpha_{ii} = 1$.  For all values of $j = 1,
2, \dots, N$, follow the next steps: For each value of $i = 1, 2,
\dots, j$, solve for $\beta_{ij}$ using the relationship:
\begin{equation}
\beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj}
\end{equation}
For the values of $i = j+1, j+2, \dots, N$, solve for the values of
$\alpha_{ij}$ with the following:
%
\begin{equation}
\alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right)
\end{equation}

\subsubsection{Calculate a matrix determinant}

The determinant $D$ of a matrix $a_{ij}$ is calculated from the
product of the diagonal elements of the LU decomposition: 
\begin{equation}
D = P_{i=1}{N} U_{ii}
\end{equation}

This shall be calculated using the GSL function
\code{gsl_linalg_LU_det}.  If the matrix is large or the magnitude of
the elements is large, the determinant may overflow the data type.  In
these cases, the (natural) logarithm of the determinant may be
calculated instead (as the sum of the logarithms of the diagonal
elements).  In this case, the GSL function \code{gsl_linalg_LU_lndet}
shall be used.

\subsubsection{Solving a Linear Equation}

The LU decomposition of a matrix may be used to solve the
matrix equation $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $

Given the LU decomposition of a matrix into $\alpha_{ij}$ and
$\beta_{ij}$, the technique of back-substitution is used to solve the
equation above. First solve the equation $L y = B$:
\begin{eqnarray}
y_1 & = & \frac{b_1}{\alpha_{11}} \\
y_i & = & \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right)
\end{eqnarray}

The values of $y$ may be used to solve for $x_i$ using the
relationship $U x = b$:
\begin{eqnarray}
x_N & = & \frac{y_N}{\beta_{NN}} \\
x_i & = & \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right)
\end{eqnarray}

\subsubsection{Invert a matrix}

Inversion of a matrix using the LU decomposition is performed by
performing back-subsitution to solve a series of linear equations of
the form $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ where the values
of $B_j$ represent the columns of the identity matrix.  The solution
vectors $x_j$ represent the columns of the inverse matrix.  This
operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}.

\subsubsection{Perform matrix addition, subtraction and multiplication}

Matrix binary arithmetic operations differ from image binary
arithmetic operations in a fundamental way: For image operations, the
resulting pixel value is determined by performing the operation on the
two corresponding input operand pixels.  For matrix operations, the
resulting element value results from the operation on the
corresponding pair of row and colum from the input matrices.  So, for
example, given the input matrices $\alpha_{ij}$ and $\beta_{ij}$, the
image and matrix multiplications correspond to:
%
\begin{eqnarray}
\mbox{image}_{ij} & = & \alpha_{ij} \times \beta_{ij} \\
\mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki}
\end{eqnarray}
%
The matrix and image operations for addition and subtraction are
identical: the operation is performed on the corresponding elements in
each matrix.  Matrix division is not defined (to divide, the user
should first invert the matrix, and then multiply).  The matrix math
function \code{psMatrixOp} shall implement the matrix version of the
binary arithmetical operations for the following operators: $+, -,
\times$.

\subsubsection{Transpose a matrix}

The transpose of a matrix is simply the reorganization of the matrix
elements by 'flipping' the matrix along a diagonal.  For non-square
matrices with dimensions $N \times M$, the resulting matrix has
dimensions $M \times N$.  The element of the output matrix $T_{ij}$ is
given by
%
\begin{equation}
T_{ij} = M_{ji}
\end{equation}
where $M_{ij}$ is the matrix to be transposed.

\subsubsection{Convert a matrix to a vector}

Matrix-to-vector conversion is only defined for a matrix that has a
size of one in at least one dimension.  In that case, the elements
should be extracted, and placed in a vector, with care that the
\code{psType} is defined correctly --- a $1\times N$ matrix is
converted to a \code{PS_DIMEN_VECTOR}-type vector, while a $N\times 1$
matrix is converted to a \code{PS_DIMEN_TRANV}-type vector.

\subsection{(Fast) Fourier Transforms}

(Fast) Fourier Transforms (FFTs) shall be implemented using the
\href{www.fftw.org}{{\em Fastest Fourier Transform in the West} (FFTW)
library}.

\subsubsection{FFTW Plans}

FFTW requires the user to create a ``plan'' for each transform size,
the time to create which is a function of the desired speed of the
transform --- faster transforms for a given size (and machine) may be
performed if more time is spent testing plans.

In the \PS{} IPP, we will want to perform FFTs on images of common
sizes (e.g.\ $512 \times 512$) regularly.  This means that we would
gain from determining an FFTW plan for each of these common sizes.
FFTW provides a binary, \code{fftw-wisdom} which may be used to
generate and save ``wisdom''.  The location of the \code{wisdom} file
will be specified as a configuration variable for the IPP (defaulting
to \code{/etc/fftw/wisdom}).  The \code{wisdom} should be read in upon
initialisation of the PSLib FFT functions and saved at the conclusion.

\subsubsection{Function mapping}

The forward and reverse transforms call the corresponding
FFTW function to plan the transform:

\begin{tabular}{ll}
  PSLib function        & Major FFTW call \\ \hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  \code{psFFTForward()} & \code{fftw_plan_dft_r2c_2d()} \\
  \code{psFFTReverse()} & \code{fftw_plan_dft_c2r_2d()} \\
\end{tabular}

These plans should be formulated using the \code{FFTW_ESTIMATE} flag,
which will allow FFTW to default to a minimal planning time if the
wisdom has not been loaded.  Transforms should be performed out of
place to avoid the need to pad the input array to hold the output.

\subsubsection{More Complicated Functions}

\code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve
multiplication of two Fourier transforms.  In the former, the first
Fourier transform is multiplied by the complex conjugate of the second
Fourier transform to yield the Fourier transform of the
cross-correlation (NR \S 13.2).  In the latter, the two Fourier
transforms are multiplied directly to yield the Fourier transform of
the convolution (NR \S 13.1).

If the elements of the discrete Fourier transform are $C_k$, then the
the elements of the power spectrum are (NR \S 13.4):
\begin{eqnarray}
P_0     & = & \left| C_0 \right|^2 / N^2 \\
P_j     & = & \left( \left| C_j \right|^2 + \left| C_{N-j} \right|^2 \right)/ N^2 \\
P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2
\end{eqnarray}
where $j = 1, 2, \ldots, (N/2 - 1)$.

Note that we leave the issue of ``windowing'' the data up to the
caller, and choose to normalise by $1/N^2$.

\subsection{Ellipse Representations}

Images of astronomical objects may often be represented using a model
consisting of a radial profile combined with an elliptical contour.
Two common ways to measure such a shape are to fit a model to the
light distribution or to measure second-order moments, perhaps with
some weighting profile.  In the special case of a 2D Gaussian with an
elliptical contour, these representations are equivalent.  The
following discussion shows how to relate the fitted parameters and
second-order moments of a elliptical Gaussian of an arbitrary
orientation with the parameters of an unrotate elliptical Gaussian.

Consider a 2D Gaussian with an elliptical contour.  If the ellipse is
oriented with the major axis along the x-axis, then the formula for
such a Gaussian may be written $f = exp(-z)$ where
\begin{equation}
\label{aligned-ellipse}
z = \frac{x^2}{2\sigma_a^2} + \frac{y^2}{2\sigma_b^2}
\end{equation}
with $\sigma_a$ the semi-major axis and $\sigma_b$ the semi-minor axis
of the 1$\sigma$ contour.  Given such a Gaussian, we may measure its
second moments by integration, and find that the second moment tensor
is
\begin{equation} \left| \begin{array}{cc}
\sigma_a^2 & 0 \\
0 & \sigma_b^2 \\
\end{array} \right| \end{equation}

Now consider the same ellipse rotated to an arbitrary angle $\theta$.
The formula for such a Gaussian may be written $f = exp(-z)$ where
\begin{equation}
\label{rotated-ellipse}
z = \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + \sigma_{xy}xy
\end{equation}
Note that, in the above form of the equation, $\sigma_{xy}$ goes to 0
as the ellipse is rotated to be aligned with the x (or y) axis.  Thus,
in this representation, $\sigma_{xy}$ is well behaved, but does not
have the same units of length that $\sigma_x$ or $\sigma_y$ have.  Our
goal is to determine the relationships between the rotated and
unrotated components of the Gaussian formula as well as the second
moments.

To determine the behavior of $\sigma_x$, etc, under rotation, we start
with the aligned ellipse (\ref{aligned-ellipse}) and rotate the
coordinate frame by an angle $-\theta$:
\begin{equation} 
\left( \begin{array}{c} x^\prime \\ y^\prime \end{array} \right) =
\left| \begin{array}{cc} \cos \theta & -\sin \theta \\ 
                         \sin \theta & \cos \theta 
\end{array} \right|
\left( \begin{array}{c} x \\ y \end{array} \right)
\end{equation}
where $x^\prime$ and $y^\prime$ are the coordinates for the unrotated
(aligned) ellipse.  Applying this rotation to (\ref{aligned-ellipse}) yields:
\begin{equation} 
z = \frac{x^2 \cos^2 \theta + y^2 \sin^2 \theta - 2 x y \sin \theta \cos \theta}{2\sigma_a^2} +
    \frac{x^2 \sin^2 \theta + y^2 \cos^2 \theta + 2 x y \sin \theta \cos \theta}{2\sigma_b^2} 
\end{equation}
Grouping these terms together, we find:
\begin{equation} 
z = \frac{x^2}{2}(\sigma_a^{-2} \cos^2 \theta + \sigma_b^{-2}\sin^2 \theta) + 
    \frac{y^2}{2}(\sigma_b^{-2} \cos^2 \theta + \sigma_a^{-2}\sin^2 \theta) + 
    \frac{xy}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
\end{equation}
We then associate the components of this equation with those of (\ref{rotated-ellipse}) and find:
\begin{eqnarray}
\sigma_x^{-2} & = & \sigma_a^{-2} \cos^2 \theta + \sigma_b^{-2}\sin^2 \theta \\
\sigma_y^{-2} & = & \sigma_b^{-2} \cos^2 \theta + \sigma_a^{-2}\sin^2 \theta \\
\sigma_{xy}   & = & \frac{1}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
\end{eqnarray}
Replacing $\cos^2$ and $\sin^2$ with the double-angle relationships, we find:
\begin{eqnarray}
\sigma_x^{-2} & = & \frac{1}{2}(\sigma_a^{-2} + \sigma_b^{-2}) - \frac{1}{2}(\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
\sigma_y^{-2} & = & \frac{1}{2}(\sigma_a^{-2} + \sigma_b^{-2}) + \frac{1}{2}(\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
\sigma_{xy}   & = & \frac{1}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
\end{eqnarray}
These formulae thus define the values of $\sigma_x$, $\sigma_y$, and
$\sigma_{xy}$ given $\sigma_a$, $\sigma_b$, and $\theta$.  Note that
in this equation and the one above, we represent the quantities in
terms of $\sigma_b^{-2} - \sigma_a^{-2}$ which is always greater than
0, thus attributing the sign of the equation to the $\sin$ or $\cos$
term.  This is necessary to determine the angle in the proper quadrant
using the arctangent below.

With the above relationships, we may now form combinations that help
us to solve for $\sigma_a$, $\sigma_b$, and $\theta$:
\begin{eqnarray}
f_1 = \sigma_y^{-2} + \sigma_x^{-2} & = & \sigma_b^{-2} + \sigma_a^{-2} \\
f_2 = \sigma_y^{-2} - \sigma_x^{-2} & = & (\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
f_3 = \sqrt{f_2^2 + 4\sigma_{xy}^2} & = & \sigma_b^{-2} - \sigma_a^{-2}
\end{eqnarray}
From these, we may derive the equations for $\sigma_a$, $\sigma_b$, and $\theta$:
\begin{eqnarray}
\theta & = & \frac{1}{2} \arg (2 \sigma_{xy}, f_2) \\
\sigma_a & = & \sqrt{\frac{2}{f_1 - f_3}} \\
\sigma_b & = & \sqrt{\frac{2}{f_1 + f_3}}
\end{eqnarray}

The relationship between the rotated ($m_{i,j}$) and unrotated
($M_{i,j}$) second moments, the latter being equal to $\sigma_a^2$ and
$\sigma_b^2$, is derived in a similar fashion.  We start with the
point that the second moment is rotated as a tensor:
\begin{equation} 
\left| 
\begin{array}{cc}
m_{x,x} & m_{x,y} \\
m_{y,x} & m_{y,y} \\
\end{array} \right| 
=
\left| 
\begin{array}{cc}
+\cos \theta & -\sin \theta \\
+\sin \theta & +\cos \theta \\
\end{array} \right| 
\left| 
\begin{array}{cc}
M_{x,x} & 0 \\
0       & M_{y,y} \\
\end{array} \right| 
\left| 
\begin{array}{cc}
+\cos \theta & +\sin \theta \\
-\sin \theta & +\cos \theta \\
\end{array} \right| 
\end{equation}
Multiplying this out and substituting $\sigma_a^2$, $\sigma_b^2$ for $M_{x,x}$, $M_{y,y}$, we find:
\begin{eqnarray}
m_{x,x} & = & \sigma_a^{2} \cos^2 \theta + \sigma_b^{2}\sin^2 \theta \\
m_{y,y} & = & \sigma_b^{2} \cos^2 \theta + \sigma_a^{2}\sin^2 \theta \\
m_{x,y} & = & \frac{1}{2} \sin (2 \theta) (\sigma_a^2 - \sigma_b^2)
\end{eqnarray}
Using the double-angle relationships, these become:
\begin{eqnarray}
m_{x,x} & = & \frac{1}{2}(\sigma_a^{2} + \sigma_b^{2}) + \frac{1}{2}(\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
m_{y,y} & = & \frac{1}{2}(\sigma_a^{2} + \sigma_b^{2}) - \frac{1}{2}(\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
m_{x,y} & = & \frac{1}{2} \sin (2 \theta) (\sigma_a^{2} - \sigma_b^{2})
\end{eqnarray}
These three formulae define the second moments in terms of $\sigma_a$, $\sigma_b$, and $\theta$. 

We define equivalent intermediate products to the above:
\begin{eqnarray}
g_1 = m_{x,x} + m_{y,y}   	 & = & \sigma_a^{2} + \sigma_b^{2} \\
g_2 = m_{x,x} - m_{y,y}   	 & = & (\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
g_3 = \sqrt{g_2^2 + 4 m_{x,y}^2} & = & \sigma_a^{2} - \sigma_b^{2}
\end{eqnarray}
From these, we may derive the equations for $\sigma_a$, $\sigma_b$, and $\theta$:
\begin{eqnarray}
\theta   & = & \frac{1}{2} \arg (2 m_{x,y}, g_2) \\
\sigma_a & = & \sqrt{\frac{g_1 + g_3}{2}} \\
\sigma_b & = & \sqrt{\frac{g_1 - g_3}{2}}
\end{eqnarray}

\section{PSLib Astronomy Utilities}

\subsection{Time}

Correct time representation is \emph{critical} in astronomical software.  PSLib
uses the \code{psTime} structure to represent time values.  This structure
represents a time which consists of seconds and nanoseconds in a time
system defined by the \code{psTimeType} element \code{type}.  All available
time-systems are defined in terms of the reference epoch
``1970-01-01T00:00:00Z'' (Gregorian\footnote{Gregorian Calendar -
http://en.wikipedia.org/wiki/Gregorian\_calendar}), but with minor
modifications, as needed, for features such as leap-seconds.  The first
represenatation, TAI (International Atomic Time), has seconds of uniform length
(SI seconds) and no leap-seconds.  The exact zero reference is
``1970-01-01T00:00:10Z'' UTC.  The second representation is UTC, which has
seconds of uniform length and leap-seconds as needed to adjust it to remain
within $0.9s$ of the Earth's rotation.  It has a zero-point of exactly
``1970-01-01T00:00:00Z'' UTC.

\subsubsection{Coordinated Universal Time (UTC)}

Coordinated Univeral Time (UTC) is defined by the International
Telecommunication Union (ITU)\footnote{ITU website -
http://www.itu.int/home/index.html}.  It is a system of time with SI length
seconds but attempts to stay within $1s$ of UT1.  This is done by the insertion
of a ``leap-second'' whenever $\lvert UTC-UT1 \rvert \ge 0.9s$.  By
definition\footnote{UTC definition -
http://www.cl.cam.ac.uk/~mgk25/volatile/ITU-R-TF.460-4.pdf}, $UTC-TAI$ is an
integer number of seconds.  UTC went into effect on ``1972-01-01T00:00:00Z''
and is defined as being $TAI-UTC = 10s$ on that date.  For dates prior to
``1972-01-01'' a fixed offset of 10s relative to TAI will be assumed.

\begin{equation}
{\rm UTC} = {\rm TAI} - 10{\rm s} - {\rm leapseconds}
\end{equation}

Leap-seconds are declared by the International Earth Rotation and Reference
Systems Service (IERS)\footnote{IERS website - http://www.iers.org/}.
Leap-seconds only occur in the UTC time system and cannot be accurately
predicted due to variations in the Earth's rotational period.  To determine the
number of leap-second in a given UTC date a table of leap-seconds as annouced by
the IERS must be consulted.  This table will have to be updated each time a new
leap-second occurs.

For ease of conversion, UTC should be represented as the number of
seconds since the UNIX epoch of ``1970-01-01T00:00:00Z'',
non-inclusive of leap-seconds.  \tbd{what does this statement actually
mean? what is the source of time (gettimeofday?  let's be explicit
here}

\tbd{Times will always be expressed in the 'UTC timezone'.  Use of the
local timezone is forbidden.  -- this makes no sense given that we
define LST.  In any case, this statement, or somethign equivalent,
belongs in the SDRS not the ADD}

\subsubsection{International Atomic Time (TAI)}

International Atomic Time or Temps Atomique International (TAI) is a system of
time defined by the Bureau International des Poids et Mesures
(BIPM)\footnote{BIPM website - http://www.bipm.fr/} with SI length seconds as
measured at sea level.  To convert from UTC to TAI add the base delta of $10s$
and all of the accumulated leap-seconds since ``1972-01-01'' up until the UTC
date being converted.

\begin{equation}
{\rm TAI} = {\rm UTC} + 10{\rm s} + {\rm leapseconds}
\end{equation}

For ease of conversion, TAI should be represented as the number of
seconds since the UNIX epoch of ``1970-01-01T00:00:00Z''. \tbd{what
does this statement actually mean?}

\subsubsection{Leap-seconds}

Leap seconds keep UTC within 0.9s of UT1.  The offset between TAI and
UTC must be looked up from tables.  Jumps in the offset correspond to
leap seconds.

\begin{verbatim}
 1972 JUL  1 =JD 2441499.5  TAI-UTC=  11.0       S + (MJD - 41317.) X 0.0 S
 1973 JAN  1 =JD 2441683.5  TAI-UTC=  12.0       S + (MJD - 41317.) X 0.0 S
 1974 JAN  1 =JD 2442048.5  TAI-UTC=  13.0       S + (MJD - 41317.) X 0.0 S
 1975 JAN  1 =JD 2442413.5  TAI-UTC=  14.0       S + (MJD - 41317.) X 0.0 S
 1976 JAN  1 =JD 2442778.5  TAI-UTC=  15.0       S + (MJD - 41317.) X 0.0 S
 1977 JAN  1 =JD 2443144.5  TAI-UTC=  16.0       S + (MJD - 41317.) X 0.0 S
 1978 JAN  1 =JD 2443509.5  TAI-UTC=  17.0       S + (MJD - 41317.) X 0.0 S
 1979 JAN  1 =JD 2443874.5  TAI-UTC=  18.0       S + (MJD - 41317.) X 0.0 S
 1980 JAN  1 =JD 2444239.5  TAI-UTC=  19.0       S + (MJD - 41317.) X 0.0 S
 1981 JUL  1 =JD 2444786.5  TAI-UTC=  20.0       S + (MJD - 41317.) X 0.0 S
 1982 JUL  1 =JD 2445151.5  TAI-UTC=  21.0       S + (MJD - 41317.) X 0.0 S
 1983 JUL  1 =JD 2445516.5  TAI-UTC=  22.0       S + (MJD - 41317.) X 0.0 S
 1985 JUL  1 =JD 2446247.5  TAI-UTC=  23.0       S + (MJD - 41317.) X 0.0 S
 1988 JAN  1 =JD 2447161.5  TAI-UTC=  24.0       S + (MJD - 41317.) X 0.0 S
 1990 JAN  1 =JD 2447892.5  TAI-UTC=  25.0       S + (MJD - 41317.) X 0.0 S
 1991 JAN  1 =JD 2448257.5  TAI-UTC=  26.0       S + (MJD - 41317.) X 0.0 S
 1992 JUL  1 =JD 2448804.5  TAI-UTC=  27.0       S + (MJD - 41317.) X 0.0 S
 1993 JUL  1 =JD 2449169.5  TAI-UTC=  28.0       S + (MJD - 41317.) X 0.0 S
 1994 JUL  1 =JD 2449534.5  TAI-UTC=  29.0       S + (MJD - 41317.) X 0.0 S
 1996 JAN  1 =JD 2450083.5  TAI-UTC=  30.0       S + (MJD - 41317.) X 0.0 S
 1997 JUL  1 =JD 2450630.5  TAI-UTC=  31.0       S + (MJD - 41317.) X 0.0 S
 1999 JAN  1 =JD 2451179.5  TAI-UTC=  32.0       S + (MJD - 41317.) X 0.0 S
\end{verbatim}

For the present time, it should be assumed that this table resides on
local disk in a known location (i.e., there is no need that it is
downloaded from the internet by PSLib).  Later, the location of this
file will be made configurable. This data is available from
USNO\footnote{ftp://maia.usno.navy.mil/ser7/tai-utc.dat}.

\subsubsection{Universal Time (UT1)}

UT1 is directly tied to the rotation of the Earth.  Historically, time
has been measured with respect to the rising and setting of the
Sun. However, in the modern era of atomic clocks, the rotation of the
Earth makes for a highly unstable time standard. Tidal effects,
changes in the angular momentum of the atmosphere, seasonal changes in
the polar ice caps, movement within the Earth's core, and other
effects all cause measurable changes in the Earth's rotation on a
daily basis.  However, UT1 is still vitally important for determining
the orientation of the Earth with respect to the sky.  UT1 is
calculated by applying the value UT1-UTC to the value of UTC.  UT1 is
continuously measured by the International Earth Rotation Service, and
tabulated values of the offset of UT1 from UTC are published at
regular intervals, along with predicted future values.  The process of
calculating the UT1-UTC offsets are discussed in the section on Earth
Orientation (Section~\ref{sec:ut1}).

\subsubsection{Gregorian dates to seconds}

The Perl code below, based on an algorithm described in the book ``Calendrical
Calculations''\footnote{Calendrical Calculations -
http://emr.cs.iit.edu/home/reingold/calendar-book/second-edition/} and modified
to return seconds, converts from Gregorian-formatted dates to seconds since the
UNIX epoch.

Given year, month, day as \code{$y, $m, $d}.
\begin{verbatim}
    use integer;

    my $adj;

    # make month in range 3..14 (treat Jan & Feb as months 13..14 of
    # prev year)
    if ( $m <= 2 )
    {
        $y -= ( $adj = ( 14 - $m ) / 12 );
        $m += 12 * $adj;
    }
    elsif ( $m > 14 )
    {
        $y += ( $adj = ( $m - 3 ) / 12 );
        $m -= 12 * $adj;
    }

    # make year positive (oh, for a use integer 'sane_div'!)
    if ( $y < 0 )
    {
        $d -= 146097 * ( $adj = ( 399 - $y ) / 400 );
        $y += 400 * $adj;
    }

    # add: day of month, days of previous 0-11 month period that began
    # w/March, days of previous 0-399 year period that began w/March
    # of a 400-multiple year), days of any 400-year periods before
    # that, and 306 days to adjust from Mar 1, year 0-relative to Jan
    # 1, year 1-relative (whew)

    $d += ( $m * 367 - 1094 ) / 12 + $y % 100 * 1461 / 4 +
          ( $y / 100 * 36524 + $y / 400 ) - 306;

    # convert from count of days to seconds since the UNIX epoch
    $unix = ( ( $d - 1 ) * 86400 ) - 62135596800;
    $utc = $unix - leapseconds($unix);
\end{verbatim}
Outputs seconds as \code{$utc}.

To go the other way:

Given the number of seconds since the UNIX epoch as \code{$utc}.
\begin{verbatim}
    use integer;

    my $unix = $utc + leapseconds( $utc )
    $d = ( unix + 62135596800 ) / 86400
 
    my $rd = $d;

    my $yadj = 0;
    my ( $c, $y, $m );

    # add 306 days to make relative to Mar 1, 0; also adjust $d to be
    # within a range (1..2**28-1) where our calculations will work
    # with 32bit ints
    if ( $d > 2**28 - 307 )
    {
        # avoid overflow if $d close to maxint
        $yadj = ( $d - 146097 + 306 ) / 146097 + 1;
        $d -= $yadj * 146097 - 306;
    }
    elsif ( ( $d += 306 ) <= 0 )
    {
        $yadj =
          -( -$d / 146097 + 1 );    # avoid ambiguity in C division of negatives        $d -= $yadj * 146097;
    }

    $c = ( $d * 4 - 1 ) / 146097;   # calc # of centuries $d is after 29 Feb of yr 0
    $d -= $c * 146097 / 4;          # (4 centuries = 146097 days)
    $y = ( $d * 4 - 1 ) / 1461;     # calc number of years into the century,
    $d -= $y * 1461 / 4;            # again March-based (4 yrs =~ 146[01] days)
    $m = ( $d * 12 + 1093 ) / 367;  # get the month (3..14 represent March through
    $d -= ( $m * 367 - 1094 ) / 12; # February of following year)
    $y += $c * 100 + $yadj * 400;   # get the real year, which is off by
    ++$y, $m -= 12 if $m > 12;      # one if month is January or February

    if ( $_[0] )
    {
        my $dow;

        if ( $rd < -6 )
        {
            $dow = ( $rd + 6 ) % 7;
            $dow += $dow ? 8 : 1;
        }
        else
        {
            $dow = ( ( $rd + 6 ) % 7 ) + 1;
        }

        my $doy =
            $class->_end_of_last_month_day_of_year( $y, $m );

        $doy += $d;

        my $quarter;
        {
            no integer;
            $quarter = int( ( 1 / 3.1 ) * $m ) + 1;
        }

        my $qm = ( 3 * $quarter ) - 2;

        my $doq =
            ( $doy -
              $class->_end_of_last_month_day_of_year( $y, $qm )
            );
\end{verbatim}
Outputs year, month, day as \code{$y, $m, $d}.
%$

\emph{The above code was taken [and slightly altered] from
\code{DateTime.pm}\footnote{DateTime.pm -
http://search.cpan.org/~drolsky/DateTime/} (C)  2003 Dave Rolsky.
Please see the DateTime project website\footnote{DateTime project -
http://datetime.perl.org} for further details.}


\subsubsection{Julian Date and Modified Julian Date}

The follow definitions of Julian Date (JD) and Modified Julian Date (MJD) was
taken from, ``RESOLUTION B1: ON THE USE OF JULIAN DATES'' of ``The XXIIIrd
International Astronomical Union General Assembly''\footnote{RESOLUTION B1: ON
THE USE OF JULIAN DATES -
http://www.iers.org/iers/earth/resolutions/UAI\_b1.html}.

\paragraph{Julian Date}

\begin{verbatim}
1. Julian day number (JDN)

The Julian day number associated with the solar day is the number assigned to a
day in a continuous count of days beginning with the Julian day number 0
assigned to the day starting at Greenwich mean noon on 1 January 4713 BC,
Julian proleptic calendar -4712.

2. Julian Date (JD)

The Julian Date (JD) of any instant is the Julian day number for the preceding
noon plus the fraction of the day since that instant. A Julian Date begins at
12h 0m 0s and is composed of 86400 seconds. To determine time intervals in a
uniform time system it is necessary to express the JD in a uniform time scale.
For that purpose it is recommended that JD be specified as SI seconds in
Terrestrial Time (TT) where the length of day is 86,400 SI seconds.

In some cases it may be necessary to specify Julian Date using a different time
scale. (See Seidelmann, 1992, for an explanation of the various time scales in
use). The time scale used should be indicated when required such as JD(UT1). It
should be noted that time intervals calculated from differences of Julian Dates
specified in non-uniform time scales, such as UTC, may need to be corrected for
changes in time scales (e.g. leap seconds).
\end{verbatim}

\paragraph{Modified Julian Date}

\begin{verbatim}
"that for those cases where it is convenient to employ a day beginning at
midnight, the Modified Julian Date (equivalent to the Julian Date minus 2 400
000.5) be used"
\end{verbatim}

\paragraph{JD and MJD conversion}

Conversion between \code{psTime} values and MJD and JD are determined
from:

Where \code{psTime} is a \code{PS_TIME_TAI}.
\begin{verbatim}
mjd = psTime.sec/86400.0 + psTime.nsec/86400000000000.0 + 40587.0;
 jd = psTime.sec/86400.0 + psTime.nsec/86400000000000.0 + 2440587.5;
\end{verbatim}

For reference $2451545.0$ JD $= 51544.5$ MJD is equivalent to
``2000-01-01T00:00:00Z''.

\begin{equation}
{\rm JD} = {\rm MJD} + 2400000.5
\end{equation}

\subsubsection{Terrestrial Time (TT)}

Terrestrial Time (TT) is defined as a fixed offset from TAI.

\begin{equation}
{\rm TT} = {\rm TAI} + 32.184{\rm s}
\end{equation}

\subsubsection{TT as Julian Centuries since J2000.0}

The algorithm for calulating GMST requires TT formated in Julian centruies
since the J2000.0 epoch.

\begin{equation} t_u = \frac{{\rm JD}_{\rm TT} - 2451545.0}{36525}
\end{equation}

\subsubsection{UT1 as Julian Centuries since J2000.0}

The algorithm for calulating GMST requires UT1 be formated in Julian centuries
since the J2000.0 epoch.

\begin{equation}
t = \frac{{\rm JD}_{\rm UT1} - 2451545.0}{36525}
\end{equation}

\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
algorithm requires that both time inputs are expressed as Julian centuries
since J2000.0 \footnote{Expressions to implement the IAU 2000 definition of UT1
- http://www.edpsciences.org/articles//aa/abs/2003/30/aa3487/aa3487.html}.

Here $t_u$ is UT1 expressed in Julian centuries since J2000.0, and $t$
is TT expressed in Julian centuries since J2000.0.

\begin{eqnarray}
{\rm GMST00}(t_u, t) & = & UT1 + 24110.5493771\\
& & + 8639877.3173760\, t_u + 307.4771600\, t\\
& & + 0.0931118\, t^2 - 0.0000062\, t^3\\
& & + 0.0000013\, t^4
\end{eqnarray}

Gives $GMST00$ in seconds.

\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}

Spherical rotations may be implemented with a variety of mathematical
methods.  A single rotation may be specified by three angles which
describe rotations about the principal axes.  Figure~\ref{rotations}
shows the rotation angles for an arbitrary 3-D rotation.  A single
rotation may be represented by defining this set of Euler angles.  As
an alternative, the quaternion is a convenient construct with which to
represent a single rotation, and it has useful mathematical properties
when applying those rotations.

\begin{figure}
\begin{center}
\psfig{file=pics/rotations}
\caption{Definition of the rotation angles\label{rotations}.  Three
rotations are performed in series: first, a rotation of $\alpha_p$ is
made about the $z$ axis; second, a rotation of $\delta_p$ (not $90 -
\delta_p$ as shown) is made about the modified $y$ axis, $y'$;
finally, a rotation of $\phi_p$ is made about the modified $z$ axis,
$z''$.  Note that it is the coordinate system that rotates, not the
position of interest.}
\end{center}
\end{figure}


\subsubsection{Quaternions}

A quaternion is an ordered set of four numbers, $q = (q_0, q_1, q_2,
q_3)$, which is useful for specifying rotations.  A quaternion is made
up of a three-vector which specifies an axis about which to rotate,
and a scalar which specifies the amount of rotation.  In the
following, we call the final value, $q_3$, the scalar value; note that
other sources (e.g., MathWorld) may choose to call the first value the
scalar value.

The conjugate of a quaterion, $q = (q_0, q_1, q_2, q_3)$, is $\bar{q}
= (-q_0, -q_1, -q_2, q_3)$.  Note that the vector components are
negated, but not the scalar component.

\subsubsection{Quaternion for a position}

Given an angular position on the sky, $(\alpha, \delta)$, we can
construct a quaternion by treating it as a unit vector in cartesian
space:
\begin{eqnarray}
p_0 & = & \cos \delta \cos \alpha \\
p_1 & = & \cos \delta \sin \alpha \\
p_2 & = & \sin \delta
\end{eqnarray}
and we set the scalar value to zero, $p_3 = 0$.

Given a quaternion, $p$, we can calculate the position using the
inverse of the above equations:
\begin{eqnarray}
\phi & = & \arctan(p_1, p_0) \\
\theta & = & \arcsin(p_2)
\end{eqnarray}
where $\phi$ is the longitude and $\theta$ is the latitude.  Note that
in this case, we neglect the scalar component of the quaternion --- it
should be zero --- and that we implicitly assume that the quaternion
is normalised to unity.

\subsubsection{Quaternion for a rotation}

A rotation of angle $\theta$ about the axis defined by the unit vector
$(v_x, v_y, v_z)$ is specified by a quaternion with components:
\begin{eqnarray}
r_0 & = & v_x \sin(\theta/2) \\
r_1 & = & v_y \sin(\theta/2) \\
r_2 & = & v_z \sin(\theta/2) \\
r_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 $|r| \equiv r_0^2 + r_1^2 + r_2^2 + r_3^2
= 1$.

\subsubsection{Multiplication of quaternions}

Given two quaternions $a$ and $b$, there is a third quaternion, $p =
ab$.  The components of $p$ are given by:

\begin{eqnarray}
p_0 & = &  b_3 a_0 + b_2 a_1 - b_1 a_2 + b_0 a_3 \\
p_1 & = & -b_2 a_0 + b_3 a_1 + b_0 a_2 + b_1 a_3 \\
p_2 & = &  b_1 a_0 - b_0 a_1 + b_3 a_2 + b_2 a_3 \\
p_3 & = & -b_0 a_0 - b_1 a_1 - b_2 a_2 + b_3 a_3
\end{eqnarray}

Note that quaternion multiplication is associative (whether you do the
left pair or the right pair first doesn't matter):
\begin{equation}
(ab)c = a(bc)
\end{equation}
but not commutative (you can't switch the order of the operands):
\begin{equation}
abc \ne acb
\end{equation}

\subsubsection{Rotating a Vector}

Rotation of a position is performed by constructing the quaternion for
the position, $p$, and the rotation, $r$, according to the above
equations, and calculating the product:
\begin{equation}
q = r p \bar{r}
\end{equation}
$q$ is the quaternion of the result.  Note the use of the conjugate of
the rotation quaternion.

A general rotation may be specified by three individual rotations
about a predefined set of axes.  We choose to specify rotations around
the $z$, $y$ and $z$ axes, in that order.  The amount of rotation
around each of these axes are known as Euler angles.  Given the Euler
angles of a rotation, the rotation may be performed by rotating in
turn about the designated axes.  Euler angles are specified below for
the various rotations required.  To use them, the following rotation
quaternions are used:

\begin{itemize}
\item First, about the Z axis:
\begin{eqnarray}
r_0 & = & 0 \\
r_1 & = & 0 \\
r_2 & = & \sin(\alpha_p/2) \\
r_3 & = & \cos(\alpha_p/2)
\end{eqnarray}
\item Second, about the Y axis:
\begin{eqnarray}
s_0 & = & 0 \\
s_1 & = & \sin(\delta_p/2) \\
s_2 & = & 0 \\
s_3 & = & \cos(\delta_p/2)
\end{eqnarray}
\item Finally, about the Z axis again:
\begin{eqnarray}
t_0 & = & 0 \\
t_1 & = & 0 \\
t_2 & = & \sin(\phi_p/2) \\
t_3 & = & \cos(\phi_p/2)
\end{eqnarray}
\end{itemize}

These three quaternions may be multiplied together to yield the
quaternion of the combined rotation: $tsr$ (note the order --- $r$ is
done first, so it is nearest the position quaternion, etc.).

\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,
and ICRS, or Celestial, coordinates) is equivalent to a rotation in
3D.  Given two coordinate system, $\alpha,\delta$ and $\phi,\theta$
which differ by only a rotation, the transformation between these two
systems are defined by the following three parameters:
\begin{itemize}
\item $\alpha_p$ : the longitude of the target system pole in the
  source system
\item $\delta_p$ : the latitutde of the target system pole in the
  source system.  
\item $\phi_p$ : the longitude of the ascending node in the target system
\end{itemize}
Note that $\theta_p$, the latitude of the source system pole in the
target system, is equal to $\delta_p$ by symmetry.  Transformations
between arbitrary systems are specified by determining the appropriate
values of $\alpha_p$, $\delta_p$, and $\phi_p$, and then constructing
the quarternion for this transformation.

The relevant trigonometric relationships are:
%
\begin{eqnarray}
\sin \theta                        & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \sin (\alpha - \alpha_p) \\
\cos \theta \sin (\phi - \phi_p)   & = & \cos \delta \cos \delta_p \sin (\alpha - \alpha_p) + \sin \delta \sin \delta_p \\
\cos \theta \cos (\phi - \phi_p)   & = & \cos \delta \cos (\alpha - \alpha_p)
\end{eqnarray}
%
and for the inverse transformations, the equivalent relationships are:
%
\begin{eqnarray}
\sin \delta                          & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \sin (\phi - \phi_p) \\
\cos \delta \sin (\alpha - \alpha_p) & = & \cos \theta \cos \delta_p \sin (\phi - \phi_p) + \sin \theta \sin \delta_p \\
\cos \delta \cos (\alpha - \alpha_p) & = & \cos \theta \cos (\phi - \phi_p)
\end{eqnarray}
Since $\theta$ and $\delta$ have domains of $-\pi/2, \pi/2$, the value
of these angles are found by applying the arcsin to the sine of these
angles ($\theta = \arcsin \sin \theta$) which is always single-valued
and defined.  The value of $\alpha-\alpha_p$ may be found from
\code{atan2(y,x)}, where $y = \cos \delta \sin (\alpha - \alpha_p)$
and $x = \cos \delta \cos (\alpha - \alpha_p)$; and similarly for
$\phi-\phi_p$.

Note that the symmetry between the two sets of above equations means
that inverse transformations can be made simply by switching
$\alpha_p$ and $\phi_p$, changing the sign of $\delta_p$, and using
the forward transformation.

\subsubsection{Galactic to ICRS}

The appropriate values, from the Hipparcos and Tycho Catalogues are:
\begin{eqnarray}
\alpha_p & = & 180^\circ - 192.85948^\circ \\
\delta_p & = & 90^\circ - 62.87175^\circ \\
\phi_p & = & 90^\circ + 32.93192^\circ
\end{eqnarray}

\subsubsection{Ecliptic to ICRS}

The appropriate values, from Zombeck, are:
\begin{eqnarray}
\alpha_p & = & 270^\circ \\
\delta_p & = & 23^\circ27'8''.26 - 46''.845\, T - 0''.0059\, T^2 + 0''.00181\, T^3 \\
\phi_p & = & 90^\circ
\end{eqnarray}
where $T$ is the time in Julian centuries since 1900.

\subsubsection{Precession}

Approximate precession, good for modest sub-arcsecond precision, may
be rapidly calculated using the following rotation angles:
\begin{eqnarray}
\alpha_p & = & 180^\circ + (0^\circ.6406161\, T + 0^\circ.0000839\, T^2 + 0^\circ.0000050\, T^3) \\
\delta_p & = & 0^\circ .5567530\, T - 0^\circ.0001185\, T^2 - 0^\circ.0000116\, T^3 \\
\phi_p & = & 180^\circ + 0^\circ.6406161\, T + 0^\circ.0003041\, T^2 + 0^\circ.0000051\, T^3
\end{eqnarray}
where $T$ is $($MJD$_{\rm out} -$ MJD$_{\rm in})/36525$ is the
difference between the two epochs, in Julian centuries.  This
precession form shall be used to implement \code{PS_PRECESS_ROUGH}.

\subsubsection{Suggested test cases}

$(\alpha,\delta) = (0^\circ,0^\circ)$ transforms to Galactic
coordinates $(l,b) = (96.337272^\circ,-60.188553^\circ)$, and Ecliptic
coordinates $(\lambda,\beta) = (0^\circ,0^\circ)$.

$(\alpha,\delta) = (0^\circ,90^\circ)$ transforms to Galactic coordinates
$(l,b) = (122.93192^\circ,27.12825^\circ)$, and Ecliptic coordinates
at J2000.0 (i.e., $T=1$), $(\lambda,\beta) =
(90^\circ,66.560719^\circ)$.

$(\alpha,\delta) = (180^\circ,30^\circ)$ transforms to Galactic
coordinates $(l,b) = (195.639488^\circ,78.353806^\circ)$, and Ecliptic
coordinates at J2100.0 (i.e., $T=2$), $(\lambda,\beta) =
(167.072470^\circ,27.308813^\circ)$.

For each of the above input coordinates, precessing from J2100 to
J1900 gives the following output coordinates:
$(357.437^\circ,-1.113^\circ)$, $(358.719^\circ,88.886^\circ)$ and
$(177.423^\circ,31.113^\circ)$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Sky to Tangent Plane}

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.

\subsubsection{Reference Implementations}

There are two reference implementations 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.\footnote{http://maia.usno.navy.mil/conv2003.html} The second
reference implementation is the SOFA software package managed by the
IAU.\footnote{http://www.iau-sofa.rl.ac.uk} Only the 2003-04-29
version of SOFA should be considered.  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.  This code may be using as
a comparison for testing purposes.

\subsubsection{Coordinate Systems}

\begin{figure}
\begin{center}
\psfig{file=pics/earthrot.ps}
\caption{Coordinates systems and the transformations between them\label{earthrot}}
\end{center}
\end{figure}

Figure~\ref{earthrot} 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{1 - \mu^2}{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^2}$.

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 0.2 mas level.  For higher accuracy the user must
apply corrections to the model, which are tabulated by the IERS.

\subparagraph{IAU 2000A Precession/Nutation Model : {\tt psEOC\_PrecessionModel}}

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 2000-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 \sin ({\rm \scriptstyle {ARG_{i,j}}}) t^j 
   + (a_{{\rm c},j})_i \cos ({\rm \scriptstyle {ARG_{i,j}}})] t^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.

\subparagraph{Corrections to the Model : {\tt psEOC\_PrecessionCorr}}

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
by USNO in the table
\code{finals2000A.data}\footnote{http://maia.usno.navy.mil/ser7/finals2000A.data},
which has the format described by
\code{readme.finals2000A}\footnote{http://maia.usno.navy.mil/ser7/readme.finals2000A}. The
quantities of interest are labeled dX and dY, which should be
interpolated using Lagrangian interpolation.  Note that UT1$-$UTC and
the polar motion values are obtained from this same table.

There is no correction to apply to $S$.

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.

\subparagraph{Spherical Rotation from Polar Coordinates : {\tt psSphereRot\_CEOtoGCRS}}

In order to relate the values $X$, $Y$, and $s$ to the rotation
components, the rotation matrix below must be used.  The definitions
of $X$, $Y$, and $s$ transform from the CIP/CEO system to the GCRS
using IERS32 equation (10), reproduced below:

\begin{equation}
\label{CEOtoGCRS}
\begin{pmatrix}1-aX^2& -aXY& X\cr -aXY& 1-aY^2& Y\cr -X& -Y&
1-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~\ref{earthrot}), in order to
match the IERS reference code.  The inverse transform may be found by
inverting the resulting rotation.

\paragraph{Earth Rotation}

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'$, which give the position of the TEO with
respect to the ITRS.  The values of $x_p$ and $y_p$ are published
daily by the
IERS,\footnote{http://maia.usno.navy.mil/ser7/finals2000A.daily} with
a format described by their
\code{readme.finals2000A}\footnote{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.

\subparagraph{Polar Motion from Bulletin : {\tt psEOC\_GetPolarMotion}}

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.

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 by using the Ray tidal model
given in IERS Gazette \#13. The definition of this correction is
provided below (Section~\ref{Raymodel}).

\subparagraph{Polar Motion Nutation Correction : {\tt psEOC\_NutationCorr}}

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'$.

\subparagraph{Spherical Rotation from Polar Motion : {\tt psSphereRot\_ITRStoTEO}}

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 Y
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~\ref{earthrot}), in order to match the reference
implementation.

\subsubsection{Universal Time (UT1)}
\label{sec:ut1}

Since 2003-01-01, UT1 has been defined as directly proportional to the
Earth Rotation Angle (see IERS Technical Note 32\footnote{IERS
Technical Note 32 - http://maia.usno.navy.mil/conv2003.html}).
Previous to that date, a different definition was in effect (see IERS
Technical Note 21\footnote{IERS Technical Note 21 -
http://maia.usno.navy.mil/conventions.html}).  We will always use the
post-2003 definition.

UT1 is continuously measured by the International Earth Rotation Service, and
tabulated values of the offset of UT1 from UTC are published at regular
intervals, along with predicted future values.  IERS Bulletin A gives ``rapid
response'' values necessary for real-time and near real-time data analysis
(such as Pan-STARRS Otis and IPP subsystems). Bulletin B gives the results of a
final, definitive data reduction.  An amalgam of Bulletin A and B values is
published daily on the IERS website\footnote{IERS Bulletin A \& B -
http://maia.usno.navy.mil/ser7/finals2000A.daily} along with a desciption of
the format\footnote{IERS finals2000A.daily table format -
http://maia.usno.navy.mil/ser7/readme.finals2000A}.

The UT1 offsets should be interpolated using the prescription of IERS
Gazette \#13\footnote{IERS Gazette \#13 -
http://maia.usno.navy.mil/iers-gaz13}.  This entails using a third
order polynomial to interpolate the table values, and then applying a
model for diurnal and semi-diurnal fluctuations due to tidal effects.
An example implimentation\footnote{interp.f -
ftp://maia.usno.navy.mil/dist/interp.f} written in Fortran is
available.  The interpolated value of $dT$ must then have the tidal
correction from the Ray Tidal Model applied.

\subsubsection{Ray Tidal Model : {\tt psEOC\_PolarTideCorr}}
\label{Raymodel}

The Ray Model tidal corrections to X, Y, and dT are given by the the
Fortran code listed below.  The input information is the epoch of
interest in MJD, while the output results are the corrections $C_x$,
$C_y$, and $C_dt$, which are in turn added to the interpolated values
determined above.

\verbatiminput{raymodel.f}

The Polar motion X, and Y coordinates are also important for determining the
orientation of the Earth with respect to the sky. This is also given in the
IERS publications references above, and should be interpolated in the same way.

\subsubsection{Longitude}

Longitudes are often expressed in the form of decimal degrees while the
algorithm for calculating GMST outputs seconds.

\begin{equation}
1\degree = 240s
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\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}
\newcommand\citep[1]{\textit{#1}}
\newcommand\citet[1]{\textit{#1}}

{\em The following discussion is adapted from an article by Ken Chambers}

The hypsometric structure and index of refraction of the Earth's
atmosphere produces:
\begin{itemize}
\item atmospheric refraction, resulting in an apparent positional
displacement of astronomical objects towards smaller topocentric
zenith distances,
\item chromatic dispersion, along great circles intersecting the
topocentric zenith with shorter wavelengths having smaller zenith
distances; and
\item extinction, from scattering and absorption of light by
atmospheric gases (including water vapor) and aerosols.
\end{itemize}

Atmospheric refraction $R(\lambda) = z_{vac} - z_1$ is the difference
between the topocentric zenith angle {\it in vacuo}, where the index
of refraction is $n \equiv 1$, and the observed refracted zenith angle
at the observatory $z_1$.  (All subscripts ``1" in this discussion
indicate the local values of quantities at a given observatory site,
and all units are SI.)  There are various ways to express the equation
for zenith angle refraction; common ones include the ``general
equation", e.g. \citep{s1936,s1962}, the Auer \& Standish
transformation \citep{as1979,as2000}, and the Saastamoinen
approximation \citep{saas73a,saas73b,saas73c}.  All are derived from
the ``refractive invariant'', the fact that along the refracted light
path described by the locus of points $r(z,n)$, where $r$ is the
radial distance from the center of the Earth and $z$ and $n$ are the
local values of the refracted zenith angle and the index of refraction
of air respectively, the product $ n r {\rm sin} z = \rm {constant} $
remains invariant.

We have revisted the Saastamoinen approximation in search of an
improved analytical expression accurate to $\sim 1$ mas up to zenith
angles of 75 degrees.  Although different approximations are made at
different stages, both the general equation and the Saastamoinen
approximation lead to an expression for the refraction $R(\lambda)$ in
the form of a series in odd powers of tan$z_1$:
\begin{equation}
 R(\lambda) = R_1 {\rm tan} z_1 + R_3 {\rm tan}^3 z_1 + 
              R_5 {\rm tan}^5 z_1 + R_7 {\rm tan}^7 z_1 + 
              R_9 {\rm tan}^9 z_1 + \dots                    
\end{equation}
where we have chosen to label the coefficients with a subscript reflecting
the exponent of the tan$z_1$ terms.
We find (c.f. Stone 1996)
%
\footnote{ The $R_1$ coefficient and the first two terms in the $R_2$
          coefficient in equation 2 above give results similar to that
          of Stone's equation (4), if it is modified such that the
          $\kappa$ factor multiplies only his $\beta$, rather than the
          whole coefficient as is done for both coefficients.  We
          suspect this is a typographical error; his equation equation
          as written gives much greater residuals when compared to the
          Pulkovo Refraction Tables than he reports in his Table 2,
          whereas the agreement is comparable, if his equation is
          modified as described, and computed for the Pulkovo baseline
          model of $\lambda = 590$nm, $15\deg C$, 101325 Pa, zero
          water vapor, and mean sea level at latitude 45 degrees.  An
          exact comparison is not possible because the residuals in
          his Table 2 are averages computed from a wide but
          unspecified range of meteorological conditions.  }
%
\begin{equation}
 R(\lambda) = \gamma_1(1 - h_1/r_1){\rm tan} z_1 + 
                [\gamma_1(\gamma_1 / 2  - h_1/r_1) + \delta_1]
                {\rm tan}^3 z_1                                  
\end{equation}
in radians, where we define for convenience the variable 
$$ \gamma_1 = n_1 -1 $$ where $n_1(\lambda)$ is the index of
refraction of air at the observatory, which is dependent on the
observatory's altitude and the meteorological conditions at the time
of the observation (see Section 3 below).  Each of the other variables
Eq.(2) take detailed discussions which, for the sake of clarity, are
divided into separate subsections.

\subparagraph{Observatory height}

The height of the observatory from the geometric center 
of the Earth is 
\begin{equation}
\begin{array}{ll}
 r_1 = r_e + r_h, \qquad  r_h \approx r_n + r_o,                    
\end{array}
\end{equation}
where 
$r_e$ is the local radius of the reference ellipsoid,
$r_h$ is the local height above the reference ellipsoid,
$r_n$ is the local geoid height, normal to the reference ellipsoid,  
$r_o$ is the orthometric height, which is the height above the geoid
      (or Mean Sea Level) in the direction of normal gravity
      (i.e. a local plumb line). 

\subparagraph{The magnitude of normal gravity at the observatory }
                                                                                
The local magnitude of normal gravity
%
\footnote{ Stone uses a common
  \citep[e.g.][]{allen1973,seid1992,allen2001} expression for the
  normal gravity as a function of lattitude and altitude, which
  apparently first appeared in \citet{lamb1949}.  However its
  derivation is unclear -- cited as {\it ``Some notes on the
  calculation of the geopotential", unpublished manuscript},
  \citet{lamb1949}.  }
%
$g_1$
is the acceleration due to the combination
of the gravity of the ellipsoid, the local mass distribution,
and the centripetal acceleration from the Earth's rotation,
the vector sum being directed opposite to the local zenith by definition,
being close but not identical to the normal to the ellipsoidal, but
not intersecting the the geometric center of the Earth.
i.e. the atmospheric topocentric zenith differs from the geocentric
astronomical zenith, the impact of this difference on calculating
the atmospheric refraction is minimal, see \citep{seid1992}. 
The acceleration at the observatory is 
\begin{equation}
 g_1(r_h,\phi) = g(\phi)\left[ 1 - 2( 1 + f + m_r - 2 f \sin^2\phi)
                 \left(\frac{r_h}{a} \right) +
               3 \left(\frac{r_h}{a} \right)^2 \right]
\end{equation}
where
\begin{equation}
 g(\phi) = g_e \left( \frac{1 + k_s \sin^2 \phi}{\sqrt{1 - \epsilon^2 \sin^2 \phi}} \right)
\end{equation}
and $\phi$ is the latitude of the observatory and the other constants
are given in Table 2.
\begin{table}[!ht]
\caption{WGS-84 World Geodetic System Reference Ellipsoid for GPS}
\smallskip
\begin{center}
{\small
\begin{tabular}{lr}
%\tableline
\noalign{\smallskip}
Semi-major axis $a$                              & 6356752.3142 $m$\\
Flattening $f =(a-b)/a$                          & 0.003352811 \\
Eccentricity $\epsilon = \sqrt{a^2 +b^2}/a$      & 0.081819\\
Polar gravity $g_p$                              & 9.8321849378 $m s^{-2}$ \\
Equatorial gravity $g_e$                         & 9.7803253359 $m s^{-2}$ \\
Somigliana's Constant $k_s = ((b/a)(g_p/g_e)-1)$ & 0.001931853 \\
Angular velocity of the Earth $\omega$           & 0.00007292115 $rad/s$\\
$GM_{earth}$ including atmosphere         & 3986004.418$\times 10^8 m^3
s^{-2}$\\
Gravity ratio $m_r = \omega^2 a^2 b /(GM)$       & 0.003449787\\
\noalign{\smallskip}
%\tableline
\end{tabular}
}
\end{center}
\end{table}


\subparagraph{The scale height above the observatory} 

The scale height of the atmosphere above the observatory is
\begin{equation}
 h_1 = \mathcal{Z}_1\mathcal{R}T_1/g_1 M_a[1 - x_w(1-M_w/M_a)],  
\end{equation}
where
$T_1$ is the local air temperature at the observatory in degrees Kelvin,
$\mathcal{R} = 8.314472 {\rm J \ mol^{-1} K^{-1}}$ 
is the gas constant,
\begin{equation}
\begin{array}{ll} \mathcal{Z}_1 = 1 - (P_1/T_1)\left[a_0 + a_1t_1 + a_2t_1^2 + 
                  (b_0 + b_1)t_1 x_w + (c_0+c_1t_1)x_w^2 \right] \\
                  \qquad \qquad +(P_1/T_1)^2(d+ex_w^2)  
                \end{array}
\end{equation}
is the compressibility of moist air at local air temperature 
$t_1 = T_1 - 273.15$  
in degrees Celsius. 
The values of the other constants in $\mathcal{Z}_1$ are given in
Appendix A. 
The quantity  
$M_a = 0.0289635 + 1.2001 \times 10^{-8}(x_{c1} - 400)$ 
kg/mole
is the molar mass of dry air with a $CO_2$ concentration $x_{c1}$ 
in $\mu$mol/mol,
$M_w = 0.018015 $ kg/mole is the molar mass of water vapor, 
and
$x_w$ is the molar fraction of water vapor in moist air, which 
depends on the local humidity.   

To calculate both the scale height $h_1$ and 
the index of refraction of moist air at the observatory $n_1$
we need to calculate $x_{w1}$, the molar fraction of water vapor
from local measurements of the relative humidity or, much better, the
dew point. Simple formulas such as \citet{davis1992}
are only suitable above $0\deg C$, whereas at Mauna Kea observatories the
temperature is often below $0\deg C$, and occasionally the surrounding 
ground is covered in snow or ice, which also alters saturation vapor pressure.
Thus, we adopt the best available equation for the 
saturation water vapor pressure $p_{sv}$, 
the IAPWS equations \citep{huang1998}. 
Haung's equations are:

\begin{equation} 
\begin{array}{ll}
   \Omega = T + K_9/(T - K_{10}) \qquad A = \Omega^2 +K_1\Omega + K_2 \\ 
   B = K_3\Omega^2 + K_4\Omega + K_5 \qquad C = K_6\Omega^2 + K_7\Omega + K_8 \\
   X = - B + \sqrt{B^2 - 4AC} \qquad p_{sv}(t) = 10^6(2C/X)^4  \\
\end{array}
\end{equation} 
For saturation vapor pressure over ice or snow, use 
\begin{equation} 
\begin{array}{ll}
   \Theta = T/273.16  \qquad Y=A_1(1-\Theta^{-1.5})+A_2(1-\Theta^{-1.25}) \\ 
   p_{sv}(t) = 611.657 e^Y  \qquad  \\
\end{array}
\end{equation} 
The constants are given in Table 1.  
\begin{table}[!ht]
\caption{Constants for Compressibility and Humidity Equations} 
\smallskip
\begin{center}
{\small
\begin{tabular}{lcl}
%\tableline
\noalign{\smallskip}
 $a_0=1.58123\times 10^{-6}$KPa$^{-1}$     &  &  $K_1= 1.16705214528E+03$ \\  
 $a_1=-2.9331\times 10^{-8}$Pa$^{-1}$      &  &  $K_2= -7.24213167032E+05$ \\  
 $a_2=1.1043\times 10^{-8}$K$^-1$Pa$^{-1}$ &  &  $K_3= -1.70738469401e+01$ \\  
 $b_0=5.707\times 10^{-6}$KPa$^{-1}$       &  &  $K_4= 1.20208247025E+04$ \\  
 $b_1=-2.051\times 10^{-8}$Pa$^{-1}$       &  &  $K_5= -3.23255503223E+06$ \\  
 $c_0=1.9898\times 10^{-4}$KPa$^{-1}$      &  &  $K_6= 1.49151086135E+01$ \\  
 $c_1=-2.376\times 10^{-6}$Pa$^{-1}$       &  &  $K_7= -4.82326573616E+03$ \\  
 $d  =1.83\times 10^{-11}$K$^2$Pa$^{-2}$   &  &  $K_8= 4.05113405421E+05 $ \\  
 $e  =-0.765\times 10^{-8}$$^2$KPa$^{-2}$  &  &  $K_9= -2.38555575678E-01$ \\  
                                           &  & $K_{10}=6.50175348448E+02 $ \\  
                                           &  & $A_1=-13.928169$           \\  
                                           &  & $A_2=34.7078238$           \\  
\noalign{\smallskip}
%\tableline
\end{tabular}
}
\end{center}
\end{table}

\noindent Now, to calculate the mole fraction of water vapor $x_w$, we need the
so called ``enhancement factor"
\begin{equation}
   f(p,t)= a' + b' p + c' t^2 
\end{equation}
where $a' = 1.00062$,$b'= 3.14 \times 10^{-8}$, and 
$c' = 5.60\times 10^{-7}$, and where $p$ and $t$ are the air pressure in
Pascals and air temperature. 

If you have the more precise measurement of the dew point $t_d$ (or frost point)then
\begin{equation} 
x_w = f(p,t) \times p_{sv}(t_d)/p.
\end{equation}
On the other hand if you only have the relative humidity $RH$, a less accurate
expression is 
\begin{equation} 
x_w = (RH/100) \times f(p,t) \times p_{sv}(t_d)/p.
\end{equation}

\subparagraph{The index of refraction of moist air at the observatory}

The Ciddor equation for the index of refraction of moist air 
\citep{ciddor1996} has been adopted by the International Association
of Geodesy (IAG) as the standard as it is believed to provide the most
accurate results under the largest
range of wavelength, temperature, and humidity conditions (300 to 1690 nm;
-40 to 100 C; 0-80\% RH). Note for astronomy, the air temperature at the
tropopause in the 1976 standard atmosphere is 216.6 K, below the stated range ofvalidity. Astronomical observers often observe up to 90\% RH, where
water droplets can form and change the effective index of refraction.
Note that developments in the equation for the index of refraction of moist
air have been poorly tracked in the astronomical literature and it is critical to examine in
every application what equation is actually being used. \footnote{
Edlen's (1953) original fit to the available data covered the wavelength
range 2752 to 6440 /AA with reasonable residuals. His 1953 constants still
survive in the astronomical literature in the equations of \citet{allen1973};
%Allen (1973);
\citet{stone1996}; \citet{allen2001}; 
%Stone(1996); Allen(2001), 
Roe (2002 - who extrapolates the 1953 equation to
K band to correct for dichoric adaptive optics) and in the computer codes
SLALIB and ZEEMAX. None discuss the range of validity. However, Elden himself
revised them \citep{elden1966},
%(Elden 1966), 
and these were further discussed and updated by
\citet{pr1972, bd1993, bd1994, ciddor1996, bp1998}. 
\citet{rueg1998} and \citet{sz2004} make convincing arguments 
for the Ciddor equation, and the latter's approach is followed here.}  

The index of refraction of air at standard temperature and pressure
is \citep{ciddor1996} 
\begin{equation}
\gamma_{as} = 10^{-8} \left( { \left[ { k_1 \over k_0 - \sigma^2} \right] 
                    + \left[ { k_3\over k_2 - \sigma^2} \right] } \right),
\end{equation} 
where $\sigma = 1/\lambda$ is the wavenumber of wavelength of light in microns.
Adjusting for the (annually varying and secularly increasing) value of 
atmospheric $CO_2$ concentration $x_{CO2}$ in units of $\mu$mole/mole, 
the expression for dry air becomes  
\begin{equation}
\gamma_{axs} = \gamma_{as} 
               \left[ 1 + 5.34\times10^{-7} 
               (x_{CO2} - 450\ \mu{\rm mole/mole} )\right]. 
\end{equation} 
For water vapor under standard conditions, Ciddor finds 
\begin{equation}
\gamma_{ws} = 1.022 \times 10^{-8}
              \left[\omega_0 + \omega_1 \sigma^2 + 
              \omega_2 \sigma^4 + \omega_3 \sigma^6 \right]. 
\end{equation} 
Following \citep{owens1967}, 
the indicies can be combined in proportion to their densities,
thus the index of refraction of moist air at the observatory 
$n_1 = \gamma_1 +1 $ is given by     
\begin{equation}
\gamma_1 = (\rho_a/ \rho_{axs})\gamma_{axs} + 
           (\rho_w/ \rho_{ws})\gamma_{ws},  
\end{equation} 
where
\begin{equation} 
\begin{array}{ll}
   \rho_{a}   = (1 - x_w) P_1 M_a/({\cal Z}_m {\cal R} T_1) \\ 
   \rho_{w}   = x_w P_1 M_w/( {\cal Z}_m {\cal R} T_1) \\ 
   \rho_{axs} = P_{STP} M_a/( {\cal Z}_a {\cal R} T_{STP}) \\ 
\end{array}
\end{equation} 
and $\rho_{ws}$ is given in Table 3.  

\begin{table}[!ht]
\caption{Constants in the Ciddor Eq. for index of refraction of moist air}  
\smallskip
\begin{center}
{\small
\begin{tabular}{lll}
%\tableline
\noalign{\smallskip}
$k_0 = 23.0185\ \mu {\rm m}^{-2}$       &
$\omega_0 = 295.235\  \mu {\rm m}^{-2}$ & 
$P_{STP} = 101325$ Pa                   \\ 

$k_1 = 5792105\ \mu{\rm m}^{-2} $       &
$\omega_1 = 2.6422\   \mu {\rm m}^{-2}$ & 
$T_{STP} = 288.15 $ K                   \\

$k_2 = 57.362\ \mu{\rm m}^{-2}  $       & 
$\omega_2 = -0.03238\ \mu {\rm m}^{-4}$ & 
${\cal Z}_{a} = 0.9995922115$           \\

$k_3 = 167917\ \mu {\rm m}^{-2} $       & 
$\omega_3 = 0.004028\ \mu {\rm m}^{-6}$ &
$\rho_{ws} = 0.00985938\ {\rm kg\ m}^3$ \\
\noalign{\smallskip}
%\tableline
\end{tabular}
}
\end{center}
\end{table}

\subparagraph{The tropopause term in the equation of refraction}

The final term in the Refraction Equation (2) 
is $\delta_1$, which comes from our re-derivation of the 
Saastamoinen approximation: 
\begin{equation}
 \delta_1 = 5 \left( {h_1 \over {r_1 T_1}}\right)^2
              \left[ {{{\gamma_{ft}T_{ft}^2 - \gamma_t T_t^2} \over
                                {1 - (h_1\beta/T_1)}}} 
                                 + \gamma_tT_t^2 \right]   
\end{equation}
where $\beta$ is the lapse rate of the troposphere, 
and the index of refraction of the free troposphere is given by 
\begin{equation}
 \gamma_{ft} = \gamma_t(T_{ft}/T_t) ^{-(T_1/h_1 \beta)-1}       
\end{equation}
and the index of refraction of the tropopause is given by  
\begin{equation}
 \gamma_t = \gamma_1 
           {\rm exp} \left[ {T_1 (r_t - r_1) \over T_t h_1 }\right]    
\end{equation}
where $\beta$ is the lapse rate in K/m,  
the temperature of the tropopause $T_t$, and height of the tropopause $r_t$ 
are all determined from contemporaneous meterological data 
(radiosonde or modern forecast models). Then the temperature of the 
free troposphere is given by 
\begin{equation}
 T_{ft} = T_t - \beta(r_t - r_1).    
\end{equation}

\subparagraph{Calculating the atmospheric refraction from both 
  the observed and true zenith angle}

The monochromatic refraction can now be calculated for any given wavelength
$\lambda_{air}$
(formally only within the range of validity - 300 to 1670 nm) given
the altitude of the observatory
$h_1$; 
contemporaneous meteorological measurements at the observatory of   
air temperaure $T_1$ (K); 
atmospheric pressure $P_1$ (Pa);
percent relative humidity $RH$, or preferably dew point temperature $t_d$ ($\deg$C);
as well as a small set of additional meteorological data:   
the lapse rate of the troposphere $\beta$ (K/m);
the radius of the tropopause $r_t = r_e+h_t$(m), where
$h_t$ is the height above MSL of the tropopause; 
the temperature of the tropopause $T_t$(K);
and the atmospheric concentration of CO$_2$ 
$x_{CO2}$ ($\mu$mole/mole). 
The characterization of the troposphere and tropopause can 
be determined either by interpolation between radiosonde measurements 
or from hourly updated meteorological models available at many observatory 
sites, or, worst case, simply adopting the 1976 Standard Atmosphere 
values:\footnote{
\citet{seid1992}
%Seidelmann 1992 
contains a typographical error quoting this value in units of K/km} 
$\beta = -0.0065$ K/m, 
$h_t = 11000$m,
$T_t = 216.6$ K and assuming 
$x_{C02} = 375 \mu$mole/mole.\footnote{
Closed rooms have higher CO$_2$ concentration, thus the STP laboratory measurements 
have concentrations near 450 $\mu$mole/mole (the prefered unit to parts per million per volume). The secular
increase of atmospheric CO$_2$ in the industrial age is well documented, with 
an annual cycle superimposed due to terrestrial biomass and ocean exchange. 
}


If the number of electrons being created from the illumination from the source in
the interval of wavelength $d\lambda$ is $N_{\lambda}d\lambda$, 
then the mean refraction is  
\begin{equation}
\bar R = {{\int R(\lambda) N_{\lambda} d \lambda} \over 
         {\int N_{\lambda} d\lambda}}  
\end{equation}

\subsubsection{Air Mass and Extinction}    

By Laplace's theorem, the monochromatic airmass (mass per unit area
along the refracted path) is 
\begin{equation}
M(z_1) = (P_1/g_1) R(\lambda) / {\rm sin} z_1
\end{equation}
in kg/m$^2$. Thus
\begin{equation}
 M(z_1) = (P_1/g_1) \left( \gamma_1(1 - h_1/r_1){\rm sec} z_1 + 
                [\gamma_1(\gamma_1 / 2  - h_1/r_1) + \delta_1]
                {\rm sec}^3 z_1 \right).                                  
\end{equation}
This is generally normalized to $P_1/g_1$, i.e. in spite of the
daily changes in barometric pressure, and thus daily changes in the 
true mass of air over the observatory, the resulting change in extinction is
generally treated as a drift in photometric zeropoint.  
For a survey program like Pan-STARRS, one could instead normalize to a 
standard barometric pressure (i.e. the altitute pressure), and thus
during a low in atmospheric pressure, the airmass would be less than
1 at zenith, and during periods of high pressure the airmass would
be greater than 1 at zenith. From the variation with zeropoint and
temperature and barometric pressure, this would remove most of the
observed variation in zeropoint in the CFHT legacy program. 
(Magnier, private comunication). 

The mean airmass is then
\begin{equation}
\bar M(z_1) = \left( P_1 \over g_1 \right) \int 
               \left( \gamma_1(1 - h_1/r_1){\rm sec} z_1 + 
                [\gamma_1(\gamma_1 / 2  - h_1/r_1) + \delta_1]
                {\rm sec}^3 z_1 \right) N_{\lambda} d \lambda.                                   
\end{equation}
and depends weakly on the filter bandpass. Use of this more
accurate expression for airmass should lead to improved extinction
corrections at high airmass.  

\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.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
to identify which chip a source lies on, when we know the coordinates
in the focal plane.  This is an example of the one-to-many problem
(one coordinate, many chips that it may lie on).

If this needs to be repeated for only one (or a small number of) focal
plane coordinates, then the fastest method is to simply convert the
focal plane coordinates to chip coordinates for each of the chips, and
determine for which of the chips the chip coordinates are valid (i.e.\
on the chip).

On the other hand, if this needs to be repeated for many source focal
plane coordinates, then it is most efficient to convert the centers of
each of the chips to coordinates on the focal plane and to use the
distance of the source to each of the centers to optimise which chips
are tested first.  This saves testing many chips for every source.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{General Astronomy Functions}

\tbd{we will provide a new airmass function, and a new mean-to-apparent conversion}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Positions of Major Solar System Objects}

\tbd{ephemerides code?}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\pagebreak 

\section{Pan-STARRS Modules}

\subsection{Object Models}

In the discussions below, the following relationships between the
given C variables and the \code{pmSource} entries should be used:
\begin{itemize}
\item \code{Xo} is \code{pmMomemt.x}
\item \code{Yo} is \code{pmMomemt.y}
\item \code{SigmaX} is \code{pmMomemt.Sx}
\item \code{SigmaY} is \code{pmMomemt.Sy}
\item \code{Zo} is \code{pmPeak.counts - pmMoment.Sky}
\item \code{So} is \code{pmMomemt.Sky}
\end{itemize}

\subsubsection{Real 2D Gaussian}

This function is a two-dimensional Gaussian with an elliptical
cross-section and a constant local background:
\begin{equation}
f(x,y) = Z_o e^{-z} + S_o
\end{equation}
where
\begin{equation}
z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
\end{equation}

Below is the relationship between the \code{psModel} parameters and
the function parameters, sample C-code implementing the function
efficiently, and the value of the derivatives. 

\begin{verbatim}
  param[0] = So;
  param[1] = Zo;
  param[2] = Xo;
  param[3] = Yo;
  param[4] = sqrt(2) / SigmaX;
  param[5] = sqrt(2) / SigmaY;
  param[6] = Sxy;

  X = x[0] - param[2];
  Y = x[1] - param[3];
  
  px = param[4]*X;
  py = param[5]*Y;

  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
  r = exp(-z);
  f = param[1]*r + param[0];
  /* f is the function value */

  q = param[1]*r;
  deriv[0] = +1;
  deriv[1] = +r;
  deriv[2] = q*(2*px*param[4] + param[6]*Y);
  deriv[3] = q*(2*py*param[5] + param[6]*X);
  deriv[4] = -2*q*px*X;
  deriv[5] = -2*q*py*Y;
  deriv[6] = -q*X*Y;
\end{verbatim}

The intial guess for the Gaussian parameters may be taken from the
moments, peak value, and local sky.

\subsubsection{Pseudo-Gaussian}

This function is a polynomial approximation of a 2D Gaussian.  The
function is very similar to the real Gaussian:
\begin{equation}
f(x,y) = Z_o (1 + z + z^2/2 + z^3/6)^{-1} + S_o
\end{equation}
where
\begin{equation}
z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
\end{equation}

Below is the relationship between the \code{psModel} parameters and
the function parameters, sample C-code implementing the function
efficiently, and the value of the derivatives:

\begin{verbatim}
  param[0] = So;
  param[1] = Zo;
  param[2] = Xo;
  param[3] = Yo;
  param[4] = sqrt(2) / SigmaX;
  param[5] = sqrt(2) / SigmaY;
  param[6] = Sxy;

  X = x[0] - param[2];
  Y = x[1] - param[3];
  
  px = param[4]*X;
  py = param[5]*Y;

  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
  t = 1 + z + 0.5*z*z;
  r = 1.0 / (t*(1 + z/3)); /* ~ exp (-Z) */
  f = param[1]*r + param[0];
  /* f is the function value */

  /* note difference from a pure gaussian: q = param[1]*r */
  q = param[1]*r*r*t;
  deriv[0] = +1;
  deriv[1] = +r;
  deriv[2] = q*(2*px*param[4] + param[6]*Y);
  deriv[3] = q*(2*py*param[5] + param[6]*X);
  deriv[4] = -2*q*px*X;
  deriv[5] = -2*q*py*Y;
  deriv[6] = -q*X*Y;
\end{verbatim}

The intial guess for the Gaussian parameters may be taken from the
moments, peak value, and local sky.

\subsubsection{Waussian}

The Waussian is a modified polynomial approximation of a 2D Gaussian,
with non-linear polynomial terms having variable coefficients, rather
than the Taylor series values of 1/2 and 1/6.  The
function is very similar to the pseudo-Gaussian:
\begin{equation}
f(x,y) = Z_o (1 + z + B_2 (z^2/2 + B_3 z^3/6))^{-1} + S_o
\end{equation}
where
\begin{equation}
z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
\end{equation}

Below is the relationship between the \code{psModel} parameters and
the function parameters, sample C-code implementing the function
efficiently, and the value of the derivatives.  Note the fudge factors
of 100 in the derivatives of $B_2$ and $B_3$: these are included to
slow the variation of these parameters, which are otherwise very
sensitive to small errors.

\begin{verbatim}
  param[0] = So;
  param[1] = Zo;
  param[2] = Xo;
  param[3] = Yo;
  param[4] = Sx;
  param[5] = Sy;
  param[6] = Sxy;
  param[7] = B2;
  param[8] = B3;

  X = x - param[2];
  Y = y - param[2];
  
  px = param[4]*X;
  py = param[5]*Y;

  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
  t = 0.5*z*z*(1 + param[8]*z/3);
  r = 1.0 / (1 + z + param[7]*t); /* ~ exp (-Z) */
  f = param[1]*r + param[0];

  /* note difference from gaussian: q = param[1]*r */
  q = param[1]*r*r*(1 + param[7]*z*(1 + param[8]*z/2));
  deriv[0] = +1;
  deriv[1] = +r;
  deriv[2] = q*(2*px*param[4] + param[6]*Y);
  deriv[3] = q*(2*py*param[5] + param[6]*X);
  deriv[4] = -2*q*px*X;
  deriv[5] = -2*q*py*Y;
  deriv[6] = -q*X*Y;
  deriv[7] = -100*param[1]*r*r*t;
  deriv[8] = -100*param[1]*r*r*param[7]*(z*z*z)/6;
  /* the values of 100 dampen the swing of param[7,8] */
\end{verbatim}

\subsubsection{Twisted Gaussian}

This function describes an object with power-law wings and a flattened
core, where the core has a different contour from the wings.  

\begin{equation}
f(x,y) = Z_{\rm pk} (1 + z_1 + z_2^M)^{-1} + Sky
\end{equation}
where
\begin{equation}
z_1 = \frac{x^2}{2\sigma_{x,in}^2} + \frac{y^2}{2\sigma_{y,in}^2} + x y \sigma_{xy,in}
z_2 = \frac{x^2}{2\sigma_{x,out}^2} + \frac{y^2}{2\sigma_{y,out}^2} + x y \sigma_{xy,out}
\end{equation}

\begin{verbatim}
  param[0]  = So;
  param[1]  = Zo;
  param[2]  = Xo;
  param[3]  = Yo;
  param[4]  = SxInner;
  param[5]  = SyInner;
  param[6]  = SxyInner;
  param[7]  = SxOuter;
  param[8]  = SyOuter;
  param[9]  = SxyOuter;
  param[10] = N;

  X = x - param[2];
  Y = y - param[3];
  
  px1 = param[4]*X;
  py1 = param[5]*Y;
  px2 = param[7]*X;
  py2 = param[8]*Y;

  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + param[4]*X*Y;
  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + param[9]*X*Y;

  r  = 1.0 / (1 + z1 + pow(z2,param[10]));
  f  = param[5]*r + param[6];

  q1 = param[5]*SQ(r);
  q2 = param[5]*SQ(r)*param[10]*pow(z2,(param[10]-1));

  deriv[0] = +1;
  deriv[1] = +r;
  deriv[2] = q1*(2*px1*param[4] + param[6]*Y) + q2*(2*px2*param[7] + param[9]*Y);
  deriv[3] = q1*(2*py1*param[5] + param[6]*X) + q2*(2*py2*param[8] + param[9]*X);

  /* these fudge factors impede the growth of param[4] beyond param[7] */
  f1 = fabs(param[7]) / fabs(param[4]);
  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
  deriv[4] = -2*q1*px1*X*f2;

  /* these fudge factors impede the growth of param[5] beyond param[8] */
  f1 = fabs(param[8]) / fabs(param[5]);
  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
  deriv[5] = -2*q1*py1*Y*f2;

  deriv[6] = -q1*X*Y;

  deriv[7] = -2*q2*px2*X;
  deriv[8] = -2*q2*py2*Y;
  deriv[9] = -q2*X*Y;
  deriv[10] = -q1*ln(z2);
\end{verbatim}

The intial guess for the Gaussian parameters may be taken from the
moments, peak value, and local sky.

\tbd{future galaxy models to be implemented}

\begin{verbatim}
float Sersic()
  param[0] = So;
  param[1] = Zo;
  param[2] = Xo;
  param[3] = Yo;
  param[4] = Sx;
  param[5] = Sy;
  param[6] = Sxy;
  param[7] = Nexp;

float SersicBulge()
  param[0]  So;
  param[1]  Zo;
  param[2]  Xo;
  param[3]  Yo;
  param[4]  SxInner;
  param[5]  SyInner;
  param[6]  SxyInner;
  param[7]  Zd;
  param[8]  SxOuter;
  param[9]  SyOuter;
  param[10] = SxyOuter;
  param[11] = Nexp;
\end{verbatim}

\subsection{WCS Translation}

The FITS World Coordinate System (WCS) standard is specified in two
papers: Greisen \& Calabretta, 2002, A\&A, 375, 1061 (Paper I); and
Calabretta \& Greisen, 2002, A\&A, 375, 1077 (Paper II).  Two further
papers (Papers III and IV) are available as drafts, and have not yet
been accepted.  All these papers, in their most up-to-date form, are
available from
\href{http://www.atnf.csiro.au/people/mcalabre/WCS/index.html}{Mark
Calabretta's page}.

Papers I and II together lay out a system for converting pixel
coordinates to celestial coordinates (RA, Dec).  However, these assume
that linear transformations, followed by projection or deprojection
are sufficient, whereas we do not expect that this is adequate to
describe the \PS{} focal plane.  Paper III has to do with spectral
coordinates, and does not concern us.  Paper IV has a proposed
mechanism for dealing with distortions which appears to be adequate to
our needs.  While the formalism has not been officially approved, and
the syntax may change, the current version of the paper provides a
means for translating the multilayered \PS{} system to FITS.
Consequently, we will use the current version (the version we consider
here is dated 22 April 2004) and update any syntax changes later as
required.

Paper IV proposes two separate distortion corrections --- before and
after a linear transformation.  Given pixel coordinates, a distortion
correction is made yielding ``corrected pixel coordinates''.  This
first distortion allows correction of the individual detector (for
example, if the detector is not flat, as may be the case for MegaCam).
A linear transformation is then performed to ``intermediate pixel
coordinates'', accounting for translation, rotation and scale.  The
second distortion correction to ``corrected intermediate pixel
coordinates'' is performed, which corrects for optical distortion.
The resultant is then scaled and deprojected onto the sky, yielding
the celestial coordinates.

\subsubsection{Implementation}

The first distortion will correct for tilts or bends of the detectors
so that pixels are in the same plane as the focal plane
(\code{psCell.toChip}).  The linear transformation will correct for
the position on the focal plane (\code{psChip.toFPA}).  The second
distortion will correct for distortions in the optics
(\code{psFPA.toTangentPlane}).  The projection will be a simple
gnomonic (``TAN'') projection.  Thus, the Paper IV formalism satisfies
our needs.

Paper IV also goes far beyond our needs, by providing several types of
distortion functions.  We are interested (at least, for now) solely in
simple polynomial distortion functions, as this is what is currently
implemented for \PS{} (i.e., we are not interested in cubic spline,
B-spline or lookup tables; nor are we interested in the use of
auxiliary variables, though this may change in the future).  In
particular, the proposed WCS system cannot handle the color and
magnitude dependence currently built into psLib's \code{psDistortion},
so we will need to specify a mean color and magnitude at which to
evaluate the spatial polynomials.

In the event that the multiple layers (\code{psCell.toChip}
$\rightarrow$ \code{psChip.toFPA} $\rightarrow$
\code{psFPA.toTangentPlane}) are not available, the short-cut
transformation (\code{psCell.toFPA}) can be used as the first WCS
distortion.  If the only available information is the ``quick and
dirty'' transformation (\code{psCell.toSky}), this may be implemented
using the linear transformation without any of the WCS distortions,
followed by a linear `projection' (``LIN'').

Reading the WCS into a psFPA can be done in the reverse order of
writing the WCS.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\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}
\input{ChangeLogADD.tex}

\end{document}
