Index: /tags/ADD-11/doc/pslib/psLibADD.tex =================================================================== --- /tags/ADD-11/doc/pslib/psLibADD.tex (revision 3773) +++ /tags/ADD-11/doc/pslib/psLibADD.tex (revision 3773) @@ -0,0 +1,2734 @@ +%%% $Id: psLibADD.tex,v 1.73 2005-04-27 19:59:04 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{11} +\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 +SLALIB Positional Astronomy Library & \code{http://star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html } \\ \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} + +The robust version of the statistics provides 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 +strongly bias the sample statistics. + +The robust statistics are calculated by constructing a histogram of +the values and performing the measurements on the histogram. The +choice of a bin size requires some care. If the data are integer +valued, the natural bin size is an integer. Otherwise, the bin should +be a fraction of an estimate of the standard deviation. Use the +$3\sigma$ clipped standard deviation as an estimator of the standard +deviation. The bin size shall be set at $\sigma_e / 10$. The +remaining steps of the algorithm are as follows: + +\begin{itemize} +\item Construct the histogram with the specified bin size. +\item Smooth the histogram by a Gaussian with $\sigma_s = \sigma_e / + 4$. +\item Find the bin with the peak value in the range $L_{\frac{1}{4}}$ + to $U_{\frac{1}{4}}$; this is the robust mode, $\mbox{mode}_r$. +\item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 4$. +\item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to + $\mbox{mode}_r + dL$. +\item The resulting fit parameters are the robust mean, +$\mbox{mean}_r$ and the robust standard deviation, $\sigma_r$. +\end{itemize} + +To determine the robust median, construct the cumulative histogram +from the histogram above. Select the bin which contains the 50th +percentile value and its two neighbors. Fit a quadratic to these +three points. The robust median value is the coordinate of the +quadratic which returns the 50\% value. For the upper and lower +quartile points, the same process should be used, choosing the three +bins in the vicinity of the upper and lower quartile points. + +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 mean, median, mode, standard deviation and +quartiles are estimated in the same manner as above. + +\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} + +\subsubsection{Levenberg-Marquardt Method} + +In the Levenberg-Marquardt Method (LMM; see NR \S 15.5), 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 only works if the second derivative of the function can be +considered negligible, as in the case of minimizing $\chi^2$. + +Given some ordinates, $x_i$, we would like to find the parameters, +$a_k$, of the function $f(x_i; a_k)$ which minimize $\chi^2$ for some +measurements, $y_i$ and associated errors, $\sigma_i$. We start with +a set of parameter guesses, $a_k$. We calculate the gradient +$\beta_k$ and the Hessian matrix $\alpha_{j,k}$ at this parameter +selection as follows: +\begin{eqnarray} +\beta_k & = & \frac{\partial \chi^2}{\partial a_k} \\ +\alpha_{j,k} & = & \sum_i \frac{1}{\sigma_i^2} \frac{\partial f(c_i)}{\partial a_k} \frac{\partial f(c_i)}{\partial a_j} +\end{eqnarray} + +We now define the new parameter guess for $a_k$ based on the gradient +and Hessian by defining $A_{j,k}$ as a variant on $\alpha_{j,k}$ as +follows: + +\begin{eqnarray} +A_{j,k} & = & \alpha_{j,k} ~ (j \ne k) \\ +A_{j,k} & = & (1 + \lambda) \alpha_{j,k} ~ (j = k) +\end{eqnarray} +% +and solve the system of equations represented by: +\begin{equation} +A_{j,k} a^\prime_k = \beta_j +\end{equation} +% +where $a^\prime_k$ represents our new attempt at a parameter guess. We +use this parameter set to evaluate the function. If the new value of +the function is lower than the previous guess, 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_k$ are the +best-fit parameters for the system. + +The covariance matrix, $C_{i,j}$, which is the inverse of the matrix +$\alpha_{j,k}$ allows simple calculation 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: +\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \] +This expression is more efficiently evaluated by factoring and +calculating the expresion as: +\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \] +\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \] + +Note that the values of $f_x$ and $f_y$ require some care. Given a +coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 - +int(f_x - 0.5)$. For example, when interpolating the value at +(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5) +and the fractional coordinate values $f_x, f_y = 0.3, 0.7$. The +resulting coordinate would be contained within the pixel at column,row +(5,5). + +\paragraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})} + +Because it would be slow to specify the size of the kernel +dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8 +pixels in each dimension (a kernel of size 2 pixels in each dimension +is handled by the bilinear interpolation). These correspond to the +options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and +\code{PS_INTERPOLATE_LANCZOS4}, respectively. + +Given a position on the input image, $(x_0,y_0)$, a kernel is derived +according to pixels local to the position: +\begin{equation} + h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N) +\end{equation} +where +\begin{eqnarray} + \delta x & = & x - x_0 \\ + \delta y & = & y - y_0 \\ + {\rm sinc}(z) & = & \sin(z)/z +\end{eqnarray} +and $N$ corresponds to the choice of kernel size. For $N = 2$, the +kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le +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: +\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \] +\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \] +where the offset coordinate ($i_o,j_o$) depends on the sign of the +sine of the angle $\theta$. If the sign of that sine is positive, the +offset coordinate is ($N_y\sin\theta$,0), otherwise it is +(0,$-N_x\sin\theta$). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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$. + +\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)} +\label{sec: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=rotations.ps} +\caption{Definition of the rotation angles\label{rotations}} +\end{center} +\end{figure} + +\subsubsection{Quaternion Construction} + +The following describes the algorithms needed to implement 3-D +rotations in terms of quaternions. A quaternion is an ordered set of +four numbers, $\bar{q} = (q_0, q_1, q_2, q_3)$. A rotation of angle +$\theta$ about the axis defined by the unit vector $(v_x, v_y, v_z)$ +has quaternion components: +\begin{eqnarray} +q_0 & = & v_x sin(\theta/2), \\ +q_1 & = & v_y sin(\theta/2), \\ +q_2 & = & v_z sin(\theta/2), and \\ +q_3 & = & cos(\theta/2). \\ +\end{eqnarray} +Note that the sine and cosine are taken of the half angle of the +rotation. Note also that this implies that the quaternion components +are normalized such that $|\bar{q}| \equiv q_0^2 + q_1^2 + q_2^2 + q_3^2 += 1$. + +The 3-vector representation of the angle of the pole is determined +from the coordinate of the pole ($\alpha_p, \delta_p$) by: +\begin{eqnarray} +v_x & = & \cos \delta_p \cos \alpha_p \\ +v_y & = & \cos \delta_p \sin \alpha_p \\ +v_x & = & \sin \delta_p \\ +\end{eqnarray} + +\subsubsection{Combining Two Rotations} + +Given two quaternions $\bar{a}$ and $\bar{b}$, there is a third +quaternion, $\bar{p}$, which represents the result of first applying +$\bar{a}$, and then $\bar{b}$. The components of $\bar{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} + +\subsubsection{Rotating a Vector} + +You may rotate a unit vector by first constructing a quaternion +$\bar{b}$, whose first three components are the components of the +unit vector, and whose fourth component is zero. To rotate this vector +by a quaternion $\bar{a}$, you apply the formula above for combining +two quaternions. The rotated vector is found in the first three +components of the resulting quaternion, $\bar{p}$. + +\subsubsection{Rotation Matrix} + +The rotation matrix representation of a rotation may be derived +directly from the quaternion representation. The following formulae +convert a quaternion to a rotation matrix: + +\begin{eqnarray} + rot_{x,x} & = & q_0 q_0 - q_1 q_1 - q_2 q_2 + q_3 q_3 \\ + rot_{y,y} & = & -q_0 q_0 + q_1 q_1 - q_2 q_2 + q_3 q_3 \\ + rot_{z,z} & = & -q_0 q_0 - q_1 q_1 + q_2 q_2 + q_3 q_3 \\ + rot_{x,y} & = & 2 (q_0 q_1 + q_2 q_3) \\ + rot_{y,x} & = & 2 (q_0 q_1 - q_2 q_3) \\ + rot_{x,z} & = & 2 (q_0 q_2 - q_1 q_3) \\ + rot_{z,x} & = & 2 (q_0 q_2 + q_1 q_3) \\ + rot_{y,z} & = & 2 (q_1 q_2 + q_0 q_3) \\ + rot_{z,y} & = & 2 (q_1 q_2 - q_0 q_3) \\ +\end{eqnarray} + +\subsubsection{Conversion to Other Representations} + +You may convert a rotation matrix, m, to a quaternion, p, with the following +code: + +\begin{verbatim} +double diag_sum[3]; +int maxi; +double recip; + +diag_sum[0]=1+m[0][0]-m[1][1]-m[2][2]; +diag_sum[1]=1-m[0][0]+m[1][1]-m[2][2]; +diag_sum[2]=1-m[0][0]-m[1][1]+m[2][2]; +diag_sum[3]=1+m[0][0]+m[1][1]+m[2][2]; + + +maxi=0; +for(i=1;i<4;++i) { + if(diag_sum[i]>diag_sum[maxi]) maxi=i; +} + + +p[maxi]=0.5*sqrt(diag_sum[maxi]); +recip=1./(4.*p[maxi]); + +if(maxi==0) { + p[1]=recip*(m[0][1]+m[1][0]); + p[2]=recip*(m[2][0]+m[0][2]); + p[3]=recip*(m[1][2]-m[2][1]); + +} else if(maxi==1) { + p[0]=recip*(m[0][1]+m[1][0]); + p[2]=recip*(m[1][2]+m[2][1]); + p[3]=recip*(m[2][0]-m[0][2]); + +} else if(maxi==2) { + p[0]=recip*(m[2][0]+m[0][2]); + p[1]=recip*(m[1][2]+m[2][1]); + p[3]=recip*(m[0][1]-m[1][0]); + +} else if(maxi==3) { + p[0]=recip*(m[1][2]-m[2][1]); + p[1]=recip*(m[2][0]-m[0][2]); + p[2]=recip*(m[0][1]-m[1][0]); +} +\end{verbatim} + +\subsection{Celestial Coordinate Conversions} + +Changes between spherical coordinate systems (ie, Ecliptic, Galactic, +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 & = & 282.85948^\circ \\ +\delta_p & = & 62.87175^\circ \\ +\phi_p & = & 32.93192^\circ \\ +\end{eqnarray} + +\subsubsection{Ecliptic to ICRS} + +The appropriate values, from Zombeck, are: +\begin{eqnarray} +\alpha_p & = & 0^\circ \\ +\delta_p & = & 23^\circ27'8''.26 - 46''.845\, T - 0''.0059\, T^2 + 0''.00181\, T^3 \\ +\phi_p & = & 0^\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 & = & 90^\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 & = & 90^\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. + +\paragraph{Reference Implementations} + +There are two reference implementatins for the code to account for the +motion of the Earth in space. The first are the sample routines +provided by the IERS to accompany chaper 5 of IERS Bulletin +32.\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=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{\mu^2 - 1}{1 - \beta\mu} +\end{equation} +where $\mu' = \hat{r}' \cdot \hat{\beta}$. + +Now, the component of $\hat{r}'$ perpendicular to $\hat{\beta}$ +(i.e. $\vec{r}_\perp'$) must point +in the same direction as $\vec{r}_\perp$, but will have a different magnitude +because $\hat{r}'$ is a unit vector. In other words, +$\vec{r}_\perp' = a\vec{r}_\perp$, for some scalar $a$. So the next step is +to calculate $a = \sqrt{(1-\mu'^2)/\vec{r}_\perp}$. + +Finally, reassemble the components of +$\hat{r}' = \mu'\hat{\beta} + a \vec{r_\perp}$. + + +\subsubsection{GCRS - ITRS} +The transformation between geocentric celestial coordinates and terrestrial +coordinates is a solid body rotation due to the motion of the Earth is space. +This is conventionally broken down into three components to isolate the +relatively rapid rotation of the Earth from the motion of its rotational axis. + +This section is largely a summary of Chapter 5 of IERS Technical Note +32\footnote{http://maia.usno.navy.mil/conv2003.html} (hereafter +IERS32), which is a description of the implementation of the +Resoltions of the XXIVth General Assembly of the IAU, available from +the same URL as above. These two documents describe a set of +conventions which have been in effect since 2003-01-01. The +conventions in effect before that date will not be implemented by +psLib. + +\paragraph{Precession/Nutation} + +The transform between the GCRS and the CIP/CEO coordinate systems is +described by the IAU 2000A precession-nutation model, which is +accurate to the 0.2 mas level. For higher accuracy the user must +apply corrections to the model, which are tabulated by the IERS. + +\subparagraph{IAU 200A 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 t^j \sin ({\rm \scriptstyle {ARG_{i,j}}}) + + (a_{{\rm c},j})_i t^j \cos ({\rm \scriptstyle {ARG_{i,j}}})] + , +\end{equation} + +where the $\rm \scriptstyle{ARG_{i,j}} = \sum_{k} w_{i,j,k} F_k$ represent linear +combinations of the fundamental arguments of nutation. + +The constants $p_j$, $w_{i,j,k}$, $(a_{{\rm s},j})_i$, and $(a_{{\rm c},j})_i$ +are given in the ASCII files: +tab5.2a.txt\footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2a.txt} (for $X$), +tab5.2b.txt\footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2b.txt} (for $Y$), and +tab5.2c.txt\footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2c.txt} (for $s+XY/2$). +Note that the expansion is given for $s+XY/2$, since this series converges +more rapidly than the one for $s$ alone. + +Each file contains a human-readable header, which includes the polynomial +coeficients, $p_j$ under the heading ``Polynomial part''. The data part of the +file lists the remaining constants, with rows cycling first through $i$, and +then through $j$. There is a separate heading each time $j$ increments. +Each row contains the following columns: + +\begin{itemize} +\item col 1 - A running index of rows in the table. +\item col 2 - The sine coeficients, $(a_{{\rm s},j})_i$ +\item col 3 - The cosine coeficients, $(a_{{\rm c},j})_i$ +\item cols 4 - 17 The weighting factors for the fundamental arguments of + nutation, $w_{i,j,k}$. +\end{itemize} + +A FORTRAN reference implementation for the precession/nutation model +is available from the +IERS.\footnote{http://maia.usno.navy.mil/conv2000/chapter5/XYS2000A.f} +The psLib results should agree with the reference implementation to +within the limits of numerical precision. + +\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.daily}\footnote{http://maia.usno.navy.mil/ser7/finals2000A.daily}, +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. Note that UT1$-$UTC and +the polar motion values are obtained from this same table. + +By convention, nutation terms with periods of less than two days are +accounted for by the corresponding polar motion. So it is sufficient +to interpolate the corrections tabulated daily by the IERS, and take +the result as instantaneous values. + +\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 X +axis by $x_p$, and finally rotating about the Z axis by $s'$. The +IERS reference implementation for this is given in the subroutine +POM2000 +\footnote{http://maia.usno.navy.mil/conv2000/chapter5/POM2000.f}. +Note that we describe the transform toward celestial coordinates +(upward in Figure~\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}} + +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} + +\tbd{add in summary of Ken's paper} + +\subsection{Projections} + +We implement three types of projections: {\em zenithal}, {\em +cylindrical} and {\em pseudocylindrical}, each requiring slightly +different handling. Our representations are based on the treatment of +projections presented by +\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \& +Calabretta (1995, ADASS, 4, 233)}. In all of these projections, we +are converting from a spherical coordinate $\alpha,\delta$ to a linear +(2-D) coordinate $x_p,y_p$. The projection is defined by the +projection type, the projection center ($\alpha_p, \delta_p$) and the +the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$). + +In the structure, \code{psProjection}, the projection type is defined +by the element \code{type}, the projection center $\alpha_p,\delta_p$ +is defined by the elements \code{R,D}, and the plate scales, +$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}. The plate +scales are applied independently to the $x$ and $y$ coordinates to +convert them to the corresponding linear units (ie, pixels): +% +\begin{eqnarray} +x_p & = & \rho_x x \\ +y_p & = & \rho_y y \\ +\end{eqnarray} +% +In the discussions below, we ignore this last step (or first step, +depending on the direction of the conversion). + +\subsubsection{Zenithal Projections} + +The {\em zenithal} projections are defined relative to a set of +spherical coordinates with pole at the center of the projection +($\alpha_p, \delta_p$), and which thus represents a coordinate system +rotated relative to the coordinate system of $\alpha, \delta$. In +this spherical coordinate system, the coordinate of longitude is +labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the +latitude, measured from the pole, is labeled $\theta$ and has domain +$0 \le \theta \le \pi$. The coordinate frame of $\phi,\theta$ is +defined so that $\phi_p$, the longitude of the target system pole, is +0.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} + +The airmass is calculated using the SLALIB function \code{sla_AIRMAS}. + +The parallactic angle is calculated using the SLALIB function \code{sla_PA}. + +%The parallax factors are calculated using the following formulae +%(Smart et al.\ 2003, A\&A, 404, 317): +%\begin{eqnarray} +%P_\xi & = & \cos \alpha \sin \lambda \cos \epsilon - \sin \alpha \cos \lambda \\ +%P_\eta & = & (\sin \epsilon \cos \delta - \cos \epsilon \sin \alpha \sin \delta) \sin \lambda - \cos \alpha \sin \delta \cos \lambda +%\end{eqnarray} +%where $\alpha$ is the Right Ascension, $\delta$ is the Declination, +%$\lambda$ is the solar longitude, and $\epsilon = 23^\circ 27'08''.26$ +%is the inclination of the ecliptic. The solar longitude is obtained +%from the ecliptic coordinates of the Sun. + +\tbd{we will provide a new mean-to-apparent conversion} + +To calculate the parallax factors, get the mean-to-apparent parameters +(\code{sla_MAPPA}) for a mean epoch of 2000.0, and, given the mean +position of interest, calculate the apparent position +(\code{sla_MAPQK}) for a parallax of 1.0 arcsec $(\alpha_1,\delta_1)$, +and a parallax of 0.0 arcsec $(\alpha_0,\delta_0)$. Then the parallax +factors in radians are: +\begin{eqnarray} +P_x & = & 3,600 (180^\circ/\pi) (\alpha_1 - \alpha_0) cos (\delta_0) \\ +P_y & = & 3,600 (180^\circ/\pi) (\delta_1 - \delta_0) +\end{eqnarray} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Positions of Major Solar System Objects} + +\tbd{ephemerides code to replace this?} + +The SLALIB function \code{SLA_RDPLAN} returns the apparent position of +a specified planet, or the Moon. + +To calculate the position of the Sun, use \code{sla_EVP} to get the +position of the earth relative to the Sun, and convert from the +cartesian coordinates to spherical using \code{sla_DCC2S}, and +calculate the position on the opposite side of the sphere ($\alpha +\rightarrow \alpha + 12 {\rm hrs}$ and $\delta \rightarrow -\delta$). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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: +\[ +f(x,y) = Z_o e^{-z} + S_o +\] +where +\[ +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} +\] + +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: +\[ +f(x,y) = Z_o (1 + z + z^2/2 + z^3/6)^{-1} + S_o +\] +where +\[ +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} +\] + +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: +\[ +f(x,y) = Z_o (1 + z + B_2 (z^2/2 + B_3 z^3/6))^{-1} + S_o +\] +where +\[ +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} +\] + +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. + +\[ +f(x,y) = Z_{\rm pk} (1 + z_1 + z_2^M)^{-1} + Sky +\] +where +\[ +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} +\] + +\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}