Index: /tags/ADD-09/doc/pslib/ChangeLogADD.tex =================================================================== --- /tags/ADD-09/doc/pslib/ChangeLogADD.tex (revision 3214) +++ /tags/ADD-09/doc/pslib/ChangeLogADD.tex (revision 3214) @@ -0,0 +1,29 @@ + +\subsection{Changes from version 06 to version 07} + +\begin{itemize} +\item Reworked discussion about lookup tables for \code{psTime}. +\item Moved discussion of lookup tables for $\Delta$UT, $x_p$ and + $y_p$ into the SDRS, along with more thorough specification. +\item Added discussion about using errors in calculating statistics. +\item Fixed up LMM to be specific to $\chi^2$ fitting. +\item Generalised entry on Gaussian smoothing. +\item modifications to the document name and PSDC number (and refs). +\end{itemize} + +\subsection{Changes from version 07 to version 08} + +\begin{itemize} +\item Added short section on histograms in the presence of errors. +\item Added short note on inverse spherical transformations. +\item Added section on astronomical object models +\end{itemize} + +\subsection{Changes from version 08 to version 09} + +\begin{itemize} +\item Added section on inverse and combined transformations. +\item Added \code{PS_RESAMPLE_LANCZOS[234]}, dropped + \code{PS_RESAMPLE_LAGRANGE}. +\item Added section on FITS WCS. +\end{itemize} Index: /tags/ADD-09/doc/pslib/psLibADD.tex =================================================================== --- /tags/ADD-09/doc/pslib/psLibADD.tex (revision 3214) +++ /tags/ADD-09/doc/pslib/psLibADD.tex (revision 3214) @@ -0,0 +1,2142 @@ +%%% $Id: psLibADD.tex,v 1.63 2005-02-14 21:06:55 eugene Exp $ +\documentclass[panstarrs]{panstarrs} + +\usepackage{amsmath} + +% 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{09} +\docnumber{PSDC-430-006} + +\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--05 & ??? & ??? \\ \hline +06 & 2004 Sep 7 & 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 +\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} + +\subsection{Math Utilities} + +\subsubsection{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]]}. + +\subsubsection{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} + +\subsubsection{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. + +\paragraph{Sample Statistics} + +We define the following statistical terms, assuming there is a set of +data elements $x_i$ with (standard) errors $\sigma_i$. + +\subparagraph{Mean} + +The simple mean is defined as: +\begin{equation} +\bar{x} = \frac {1}{N} \sum^N_i x_i +\end{equation} + +\subparagraph{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. + +\subparagraph{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. + +\subparagraph{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. + +\subparagraph{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} + + +\paragraph{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$. + +\paragraph{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. + +\paragraph{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. + +\subsubsection{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}. + +\paragraph{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} + +\paragraph{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. + +\paragraph{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} + +\paragraph{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}. + +\paragraph{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$. + +\paragraph{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. + +\paragraph{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. + +\subsubsection{Fitting} + +\paragraph{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} + +\paragraph{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}). + +\subsubsection{Non-linear Minimization} + +\paragraph{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} + + +\paragraph{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\S10.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}. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{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$. + +\paragraph{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. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{(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}. + +\paragraph{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. + +\paragraph{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. + +\paragraph{More Complicated Functions} + +\code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve +multiplication of two Fourier transforms. In the former, the first +Fourier transform is multiplied by the complex conjugate of the second +Fourier transform to yield the Fourier transform of the +cross-correlation (NR \S 13.2). In the latter, the two Fourier +transforms are multiplied directly to yield the Fourier transform of +the convolution (NR \S 13.1). + +If the elements of the discrete Fourier transform are $C_k$, then the +the elements of the power spectrum are (NR \S 13.4): +\begin{eqnarray} +P_0 & = & \left| C_0 \right|^2 / N^2 \\ +P_j & = & \left( \left| C_j \right|^2 + \left| C_{N-j} \right|^2 \right)/ N^2 \\ +P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2 \\ +\end{eqnarray} +where $j = 1, 2, \ldots, (N/2 - 1)$. + +Note that we leave the issue of ``windowing'' the data up to the +caller, and choose to normalise by $1/N^2$. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Astronomy Utilities} + +Most of the astronomy utilities will be implemented through wrapping +the +\href{http://star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{SLALIB +Positional Astronomy Library}. \tbd{SLAlib support will be dropped in +the next release} + +\subsubsection{Time} + +Correct time representation is critical in astronomical software. +PSLib uses the \code{psTime} structure to represent time values. This +structure represents a time which is consists of seconds and fractions +of seconds in a time system defined by the \code{psTimeType} element +\code{type}. Two possible time systems are currently available: TAI +and UTC. Both are defined in terms of the reference epoch +1970-01-01T00:00:00Z, but with minor modifications for leap seconds as +needed. The first represenatation, TAI (International Atomic Time), +has seconds of uniform length and no leap seconds. The exact zero +reference is 1970/01/01,00:00:10 UTC. The second representation is +UTC, which has seconds of uniform length and leap seconds as needed to +adjust it to remain within 0.9 seconds of the Earth's rotation. It +has a zero-point of exactly 1970/01/01,00:00:00 UTC. + +\paragraph{Coordinated Universal Time (UTC)} + +Coordinated Univeral Time (UTC) is a system of time with SI length +seconds but attempts to stay within 1s of UT1. This is done by the +insertion of leap second whenever UTC-UT1 $\ge$ 0.9s. By definition +UTC-TAI is an integer number of seconds. UTC went into effect on +"1972-01-01T00:00:00" 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} + +Leapseconds are declared by the International Earth Rotation and +Reference Systems Service (IERS). Leapseconds 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 leapsecond +in a given UTC date a table of leapseconds as annouced by the IERS +must be consulted. This table will have to be updated each time a new +leapsecond occurs. + +For ease of conversion, UTC should be represented as the number of +seconds since the UNIX epoch of "1970-01-01T00:00:00". + +\paragraph{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) 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 leapsecons 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:00". + +\paragraph{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 +\code{ftp://maia.usno.navy.mil/ser7/tai-utc.dat} + +\paragraph{Gregorian dates to seconds} + +The below algorithm converts from Gregorian-formatted dates to +seconds since the UNIX epoch. + +\begin{verbatim} + Given year, month, day. + + ### Make month in range 3..14 (treat Jan & Feb as months 13..14 of prev year): + if ( month <= 2 ) + { + year -= ( temp = ( 14 - month ) / 12 ) + month += 12 * temp + } + else if ( month > 14 ) + { + year += ( temp = ( month - 3 ) / 12 ) + month -= 12 * temp + } + + ### make year positive + if ( year < 0 ) + { + day -= 146097 * ( temp = ( 399 - year ) / 400 ) + year += 400 * temp + } + + ### 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 + day += ( month * 367 - 1094 ) / 12 + year % 100 * 1461 / 4 + + ( year / 100 * 36524 + year / 400 ) - 306 + + unix = ( ( day - 1 ) * 86400 ) - 62135596800 + utc = unix - leapseconds(unix) +\end{verbatim} + +To go the other way: + +\begin{verbatim} + unix = utc + leapseconds(utc) + day = ( unix + 62135596800 ) / 86400 + temp = 0 + + ### add 306 days to make relative to Mar 1, 0; also adjust day to be + ### within a range (1..2**28-1) where our calculations will work + ### with 32bit ints + if ( day > 2**28 - 307 ) + { + ### avoid overflow if day close to maxint + temp = ( day - 146097 + 306 ) / 146097 + 1 + day -= temp * 146097 - 306 + } + else if ( ( day += 306 ) <= 0 ) + { + temp = -( -day / 146097 + 1 ) ### avoid ambiguity in C division of negatives + day -= temp * 146097 + } + + cent = ( day * 4 - 1 ) / 146097 ### calc number of centuries day is after 29 Feb of yr 0 + day -= cent * 146097 / 4 ### (4 centuries = 146097 days) + year = ( day * 4 - 1 ) / 1461 ### calc number of years into the century, + day -= year * 1461 / 4 ### again March-based (4 yrs =~ 146[01] days) + month = ( day * 12 + 1093 ) / 367 ### get the month (3..14 represent March through + day -= ( month * 367 - 1094 ) / 12 ### February of following year) + year += cent * 100 + temp * 400 ### get the real year, which is off by + if ( month > 12 ) ### one if month is January or February + { + year++ + month -= 12 + } + + + Output year, month, day. +\end{verbatim} + +(Above taken from \code{DateTime.pm} (C) 2003 Dave Rolsky, available +from \code{datetime.perl.org}.) + + +\paragraph{Universal Time (UT1)} +\label{sec:ut1} + +Univseral Time is a measure of the rotation angle of the Earth. When +corrected for polar motion it is referred to as UT1. This is distict +from UT0 which does not involve corrections for polar motion. UT1 may +be calculated from UTC through a table lookup of the appropriate value +of UTC - UT1. + +\paragraph{Julian Day and Modified Julian Day} + +Julian Day (JD) and Modified Julian Day (MJD) are both continuous time +representations, with one julian day interval having a length of 86400 +TAI seconds. MJD is equal to JD - 2400000.5 and has a zero point +equal to that of TAI, while JD has a zero point 0.5 off of TAI. +Conversion between \code{psTime} values and MJD and JD are determined +from: + +\begin{verbatim} +mjd = psTime.sec/86400.0 + psTime.usec/86400000000.0 + 40587.0; + jd = psTime.sec/86400.0 + psTime.usec/86400000000.0 + 2440587.5; +\end{verbatim} + +$2451545.0$ JD $= 51544.5$ MJD is equivalent to "2000-01-01T00:00:00". + +\begin{equation} +{\rm JD} = {\rm MJD} + 2400000.5 +\end{equation} + +\paragraph{Terrestrial Dynamical Time (TDT)} + +Terrestrial Dynamical Time (TDT) is defined as a fixed offset from +TAI. Its only purpose as far as we are concerned is for its utility +in obtaining the GMST. + +\begin{equation} +{\rm TDT} = {\rm TAI} + 32.184{\rm s} +\end{equation} + +\paragraph{TDT as Julian Centuries since J2000.0} + +The algorithm for calulating GMST requires TDT formated in Julian centruies +since the J2000.0 epoch. + +\begin{equation} +t_u = \frac{{\rm JD}_{\rm TDT} - 2451545.0}{36525} +\end{equation} + +\paragraph{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} + +\paragraph{Greenwich Mean Sidereal Time (GMST)} + +Greenwich Mean Sidereal Time (GMST) is caclulated from UT1 and TDT. +This algorithm requires that both time inputs are expressed as Julian +centuries since J2000.0. + +Here $t_u$ is UT1 expressed in Julian centuries since J2000.0, and $t$ +is TDT 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. + +\paragraph{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} + +\paragraph{Local Mean Sidereal Time (LMST)} + +Local Mean Sidereal Time (LMST) is Greenwich Mean Sideral Time (GMST) +plus the observer's location in East longitude. Calculating LMST +requires the input of Universal Time (UT1), Terrestrial Dynamical Time +(TDT) and a longitude (measured East of Greenwich). + +\begin{equation} +LMST = GMST00(t_u, t) + longitude +\end{equation} + +Gives $LMST$ in seconds. + +\paragraph{Polar Coordinates} + +The polar coordinates, $x_p$ and $y_p$, required for the +transformation from tangent plane to celestial coordiates (and hence +\code{psGrommit}), may be calculated through table lookups. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Astronomical Image Manipulations} + +\paragraph{Interpolation} + +Interpolation is needed in various image manipulation operations, +including rotation and resampling. We have specified a function to +perform the interpolation using one of several possible interpolation +methods, defined below. It is important in the discussions that +follow to remember that a pixel with column,row if $i,j$ has +coordinate at the center of $i+0.5,j+0.5$ and corners with coordinates +from $i,j$ to $i+1,j+1$. Thus, the interpolation of a coordinate +$x,y$ = 5.0,4.0 is a value midway between the four pixels with +column,row of (5,4), (5,5), (6,4), (6,5). + +\subparagraph{Nearest Pixel Interpolation ({\tt PS\_INTERPOLATE\_FLAT})} + +In this interpolation, the value of the closest pixel is returned. +This is equivalent to pixel duplication or replication. + +\subparagraph{Bilinear Interpolation ({\tt PS\_INTERPOLATE\_BILINEAR})} + +In this interpolation, the value at the coordinate is calculated using +linear interpolation in two dimensions from the four nearest neighbor +pixels. The bilinear interpolation value at a coordinate $x,y$ +depends on the four nearest neighbor pixels and the fractional +distance $fx,fy$ of the given coordinates from the centers of those +four pixels. Consider four neighboring pixels at column,row of $i,j$, +$i+1,j$, $i,j+1$, and $i+1,j+1$ with pixel values $V_{0,0}$, +$V_{1,0}$, $V_{0,1}$, $V_{1,1}$. The value at $x,y$ is given by: +\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \] +This expression is more efficiently evaluated by factoring and +calculating the expresion as: +\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \] +\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \] + +Note that the values of $f_x$ and $f_y$ require some care. Given a +coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 - +int(f_x - 0.5)$. For example, when interpolating the value at +(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5) +and the fractional coordinate values $f_x, f_y = 0.3, 0.7$. The +resulting coordinate would be contained within the pixel at column,row +(5,5). + +\subparagraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})} + +Because it would be slow to specify the size of the kernel +dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8 +pixels in each dimension (a kernel of size 2 pixels in each dimension +is handled by the bilinear interpolation). These correspond to the +options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and +\code{PS_INTERPOLATE_LANCZOS4}, respectively. + +Given a position on the input image, $(x_0,y_0)$, a kernel is derived +according to pixels local to the position: +\begin{equation} + h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N) +\end{equation} +where +\begin{eqnarray} + \delta x & = & x - x_0 \\ + \delta y & = & y - y_0 \\ + {\rm sinc}(z) & = & \sin(z)/z +\end{eqnarray} +and $N$ corresponds to the choice of kernel size. For $N = 2$, the +kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le +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}. + +\paragraph{Image Cuts and Slices} + +Several functions specify operations which manipulate a collection of +pixels to return a statistic on the pixel collection. In the simplest +case, these are trivial to define: if the boundaries of the region of +interest are specified along integral pixel coordinates, then the +pixels used to measure the statistic are always an exact integer. +This is the case for the function \code{psImageSlice} which requires a +starting coordinate which is an integer and a width in both dimensions +which is an integer. For the case of the functions \code{psImageCut} +and \code{psImageRadialCut}, the situation is a bit more subtle. In +both of these cases, the region is unlikely to contain only whole +pixels and some choices must be made. + +One posibility which we reject is to identify the fractional pixels +which are overlapped by the region of interest and add that fraction +of the pixel's flux when calculating the statistic of interest. This +is computationally intensive, and not necessarily well defined for all +statistics. + +In PSLib, we instead identify the pixels overlapped by the region, use +the complete set of pixel values, treating all pixels equally, and +renormalize as needed. To perform this, the region of interest is +laid on top of the image pixels. Any pixels which overlap the region +are identified as part of the input sample. The statistic (ie, sample +mean, robust mode, etc), is then calculated on this collection of +pixels. If the output statistic is an average value, the measured +value is reported. If the output statistic is a sum value (sum of +counts, sum of pixels), then the value is renormalized by the ratio of +pixels used in the calculation to the pixel area of the region of +interest. For example, if the sum within a radial aperture is +requested, the circle of the specified radius and center is placed on +the pixel grid. Any pixels which touch the circle are then placed in +a list to be analysed. The statistic of interest is the measured for +this collection of pixels. In the case of a circular aperture which +is centered at the coordinate (2,2) and has a radius of 2, the number +of pixels which are touched by the circle is 16, while the total pixel +area of the circle is 12.57 square pixels. In this case, the pixel +sum is renormalized by the ratio (12.57/16.00). + +\subparagraph{Radial Cuts} + +Consider an image with pixels $x_i,y_i$ and a reference coordinate +$x_c, y_c$. We want to construct a radial cut by measuring statistics +for pixels in a sequence of radial annulii $r_s < r < r_e$. For each +annulus, we need to select the pixels which fall within this annulus. +The coordinates of the center of pixel $i,j$ are $i+0.5,j+0.5$. A +given pixel has a distance from the reference coordinate of $dX = x_c +- i - 0.5, dY = y_c - j - 0.5$. The pixels to be used for a given +radial annulus are all of those pixels for which $r_s < \sqrt{dX^2 + + dY^2} < r_e$. This is more efficiently calculated by comparing the +square of the radii and distances. All pixels which satisfy the above +condition are included in a specific annular radius. All average +quantities are calculated directly from the pixel ensemble +statistics. + +\subparagraph{Arbitrary Linear Cuts} + +Select the pixels which lie along a line following steps of 1 pixel +length: + +\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} + +\paragraph{Image Rotation} + +Image rotation can be performed in two possible ways under different +circumstances, identified in the following discussion. + +In the simplest case, the rotation angle is an integer multiple of 90 +degrees ($\pi/2$ rad). In these cases, the input and output pixels +have a one-to-one mapping. If the input image has dimensions of $N_x, +N_y$, then the output image will have dimensions of either $N_x, N_y$ +(for even multiples of 90 degrees) or $N_y, N_x$ (for odd multiples). + +If the angle of the rotation is not a multiple of 90, then the output +pixels necessarily result from the interpolation of several input +pixels. In this case, for an input image of dimensions $N_x, N_y$ and +rotation angle $\theta$, the output image has dimensions $Lx = |N_x +\cos \theta| + |N_y \sin \theta|$ and $Ly = |N_x \sin \theta| + |N_y +\cos \theta|$, each dimension rounded up to the nearest integer as +needed. Every pixel in the output image is in general derived from an +interpolation over 4 neighboring pixels. The coordinate of a pixel in +the output image ($i,j$) corresponds to a fractional pixel coordinate +($x,y$) in the input image according to: +\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \] +\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \] +where the offset coordinate ($i_o,j_o$) depends on the sign of the +sine of the angle $\theta$. If the sign of that sine is positive, the +offset coordinate is ($N_y\sin\theta$,0), otherwise it is +(0,$-N_x\sin\theta$). + +\subsubsection{Celestial Coordinate Conversions} + +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. + +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. + +\paragraph{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} + +\paragraph{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. + +\paragraph{Precession} + +The appropriate values, from Elixir, are: +\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. + + +\paragraph{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)$. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsubsection{2D transformations} + +In PSLib, we implement 2-dimensional transformations using +\code{psPlaneTransform}, which contains a matrix of polynomial +coefficients for each dimension. Since we are using these to model +the real world, where, for example, a particular point on the detector +maps to a particular point on the sky, we consider only +transformations that are ``one-to-one''. This makes it possible to +speak of inverse transformations, and of combining multiple +transformations. + +Given a transformation, $f(x,y)$, the inverse transformation, +$g(x,y)$, is that for which $g(f(x,y)) = (x,y)$ for $(x,y)$ over the +range of interest (not necessarily the entire set of real numbers). + +Given two transformations, $f(x,y)$ and $g(x,y)$, the combined +transformation is the transformation, $h(x,y) = g(f(x,y))$ for $(x,y)$ +over the range of interest (not necessarily the entire set of real +numbers). + +Both of these operations are straightforward if the transformation is +linear. If the function $(u,v) = f(x,y)$ is: +\begin{eqnarray} +u & = & a + bx + cy \\ +v & = & d + ex + fy +\end{eqnarray} +then the inverse transformation $(x,y) = g(u,v)$ is: +\begin{eqnarray} +x & = & (-fa+cd)/\Delta + fu/\Delta - cv/\Delta \\ +y & = & (ae-bd)/\Delta - eu/\Delta + bv/\Delta +\end{eqnarray} +where $\Delta = bf - ce$ is the matrix determinant. Given two +functions $f_i(x,y)$ for $i=1,2$: +\begin{eqnarray} +u & = & a_i + b_i x + c_i y \\ +v & = & d_i + e_i x + f_i y +\end{eqnarray} +then the combined transformation, $(u,v) = f_2(f_1(x,y))$ is: +\begin{eqnarray} +u & = & (a_2 + b_2 a_1 + c_2 d_1) + (b_2 b_1 + c_2 e_1) x + (b_2 c_1 + c_2 f_1) y \\ +v & = & (d_2 + e_2 a_1 + f_2 d_1) + (e_2 b_1 + f_2 e_1) x + (e_2 c_1 + f_2 f_1) y +\end{eqnarray} + +When the transformations are not linear, the inverse and combined +transformations can be estimated by sampling a grid over the region of +interest, calculating the transformation (or double transformation) +for each sample, and using this information to derive the best fit +transformation that produces the inverse or combined transformation. +The inverse transformation should be of the same order as that of the +forward transformation, while the combined transformation should be of +the higher order of the two component transformations. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Projections} + +We implement three types of projections: {\em zenithal}, {\em +cylindrical} and {\em pseudocylindrical}, each requiring slightly +different handling. Our representations are based on the treatment of +projections presented by +\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \& +Calabretta (1995, ADASS, 4, 233)}. In all of these projections, we +are converting from a spherical coordinate $\alpha,\delta$ to a linear +(2-D) coordinate $x_p,y_p$. The projection is defined by the +projection type, the projection center ($\alpha_p, \delta_p$) and the +the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$). + +In the structure, \code{psProjection}, the projection type is defined +by the element \code{type}, the projection center $\alpha_p,\delta_p$ +is defined by the elements \code{R,D}, and the plate scales, +$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}. The plate +scales are applied independently to the $x$ and $y$ coordinates to +convert them to the corresponding linear units (ie, pixels): +% +\begin{eqnarray} +x_p & = & \rho_x x \\ +y_p & = & \rho_y y \\ +\end{eqnarray} +% +In the discussions below, we ignore this last step (or first step, +depending on the direction of the conversion). + +\paragraph{Zenithal Projections} + +The {\em zenithal} projections are defined relative to a set of +spherical coordinates with pole at the center of the projection +($\alpha_p, \delta_p$), and which thus represents a coordinate system +rotated relative to the coordinate system of $\alpha, \delta$. In +this spherical coordinate system, the coordinate of longitude is +labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the +latitude, measured from the pole, is labeled $\theta$ and has domain +$0 \le \theta \le \pi$. The coordinate frame of $\phi,\theta$ is +defined so that $\phi_p$, the longitude of the target system pole, is +0.0. + +For an arbitrary projection center, it is necessary to convert the +spherical coordinates to be projected ($\alpha,\delta$) to the +projection spherical coordinate system coordinates ($\phi, \theta$). +In practice, we construct the following useful trigonometric +relationships between $\phi$ and $\theta$ which may be employed in the +equations of $x,y$ below: +% +\begin{eqnarray} +\sin \theta & = & \sin \delta \sin \delta_p + \cos \delta \cos \delta_p \cos (\alpha - \alpha_p) \\ +\cos \theta \cos \phi & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \cos (\alpha - \alpha_p) \\ +\cos \theta \sin \phi & = & - \cos \delta \sin (\alpha - \alpha_p) +\end{eqnarray} +% +For the inverse transformations, the equivalent relationships are: +% +\begin{eqnarray} +\sin \delta & = & \sin \theta \sin \delta_p + \cos \theta \cos \delta_p \cos \phi \\ +\cos \delta \cos (\alpha - \alpha_p) & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \cos \phi \\ +\cos \delta \sin (\alpha - \alpha_p) & = & - \cos \theta \sin \phi +\end{eqnarray} +% +For zenithal projections, the linear coordinates are related to +$\phi,\theta$ by: +% +\begin{eqnarray} +x & = & R_\theta \sin \phi \\ +y & = & -R_\theta \cos \phi +\end{eqnarray} +% +and the inverse: +% +\begin{eqnarray} +R_\theta & = & \sqrt{x^2 + y^2} \\ +\phi & = & {\rm atan} (-y,x) +\end{eqnarray} +% +The coordinates $x,y$ above are defined to be in angular units (ie, +radians). + +From these relationships, we can calculate $\alpha, \delta$ as: +% +\begin{eqnarray} +\alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\ +\delta & = & \arcsin (\sin \delta) \\ +\end{eqnarray} +% +Note that if $(x,y) = (0,0)$, then $\alpha = \alpha_p, \delta = \delta_p$. + +\subparagraph{Gnomonic} + +The Gnomonic projection (``TAN'') is a zenithal projection with +$R_\theta = \cot \theta$. The resulting relationships for $(x,y)$ and +for $\sin \theta, \cos \theta$ are: + +\begin{eqnarray} +x & = & \frac{\cos \theta \sin \phi}{\sin \theta} \\ +y & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\ +\sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\ +\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\ +\end{eqnarray} + +where $\zeta = 1 / R_\theta$. + +\subparagraph{Orthographic} + +The Orthographic projection (``SIN'') is a zenithal projection with +$R_\theta = \cos \theta$. The resulting relationships for $(x,y)$ and +for $\sin \theta, \cos \theta$ are: + +\begin{eqnarray} +x & = & \cos \theta \sin \phi \\ +y & = & -\cos \theta \cos \phi \\ +\sin \theta & = & \sqrt{1 - R_\theta^2} \\ +\cos \theta & = & R_\theta \\ +\end{eqnarray} + +\paragraph{Cylindrical and Pseudocylindrical Projections} + +The {\em cylindrical} and {\em pseudocylindrical} projections are +defined relative to a set of cylindrical coordinates whose pole is +coincident with the pole of the spherical coordinates. These +projections are particularly used for full-sky representations, and +are only defined for projection centers with $\delta_p = 0$. In this +spherical coordinate system, the coordinate of longitude is labeled +$\phi$, and has domain of $-\pi < \phi \le \pi$, while the latitude, +measured from the pole, is labeled $\theta$ and has domain $0 \le +\theta \le \pi$. The projection center longitude, $\alpha_p$ +corresponds to $\phi = 0$, thus the value of $\phi$ is determined as +$\alpha - \alpha_p$ for all such projections. + +\subparagraph{Cartesian} + +The Cartesian projection (``CAR'') is a very simple cylindrical +projection with the following relationships between $x,y$ and +$\phi,\theta$: + +\begin{eqnarray} +x & = & \phi \\ +y & = & \theta +\end{eqnarray} + +\subparagraph{Mercator} + +The Mercator projection (``MER'') is a cylindrical projection. + +\begin{eqnarray} +x & = & \phi \\ +y & = & \ln \left( \tan (\pi/4 + \theta/2) \right) \\ +{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^y \right) - \pi/2 +\end{eqnarray} + +\subparagraph{Hammer-Aitoff} + +The Hammer-Aitoff projection(``AIT'') is a pseudocylindrical projection, and is defined: + +\begin{eqnarray} +x & = & 2 \zeta \cos \theta \sin \frac{\phi}{2} \\ +y & = & \zeta \sin \theta \\ +{\rm where}\hspace{1cm} \zeta^{-1} & \equiv & \sqrt{\frac{1}{2}\left(1 + \cos \theta \cos \frac{\phi}{2} \right)} +\end{eqnarray} + +And in reverse: + +\begin{eqnarray} +\phi & = & 2 {\rm \arctan} (2z^2 - 1, x z) \\ +\theta & = & \arcsin (yz) \\ +{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x/2)^2 - y^2} +\end{eqnarray} + +\subparagraph{Parabolic} + +The Parabolic projection (``PAR'') is a pseudocylindrical projection, and is defined: + +\begin{eqnarray} +x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\ +y & = & \pi \sin \frac{\theta}{3} \\ +\end{eqnarray} + +And in reverse: + +\begin{eqnarray} +\theta & = & 3 \sin^{-1} \rho \\ +\phi & = & \frac{x}{1 - 4\rho^2} \\ +{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\ +\end{eqnarray} + +\subsubsection{Offset} + +Coordinate offsets can be either spherical offsets or linear offsets. + +A spherical offset is performed by adding the components of the +offset, after unit conversion, to the given position. The resulting +coordinates must be wrapped to within the allowed range ($-\pi$ to +$\pi$, 0 to $2\pi$). + +A linear offset is defined to be a linear offset in a tangent +projection centered on the starting coordinate with $y$ axis aligned +with the local direction or increasing Declination. This projection +is undefined only for the coordinates exactly at the north and south +poles, in which case the orientation is defined to have the $y$ axis +parallel to the line of RA = 0.0. The scale of the projection is 1.0 +(ie, 1 'pixel' is 1 radian) and the given offsets must the scaled +based on the given offset units. + +Pseudo-code to implement the above for an offset: + +\begin{verbatim} +psSphere *psSphereSetOffset (psSphere pos, psSphere offset) { + + psPlane lin; + psSphere new; + psProjection proj; + + proj.R = pos->r; + proj.D = pos->d; + proj.X = 0; + proj.Y = 0; + proj.type = PS_PROJ_TAN; + + lin.x = offset.r; + lin.y = offset.d; + + new = psDeproject (&lin, &proj); + return (new); +} +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Tangent Plane to Sky} + +\tbd{we will replace the SLALIB version of AOPPA with a new function} + +Mappings between the tangent plane and the sky will be implemented +using SLALIB. + +To speed up the transformations, SLALIB allows the storage of +``apparent-to-observed parameters'', which will be contained in a +\code{psGrommit}. The \code{psGrommit} consists of the following: +\begin{enumerate} +\item geodetic latitude (radians) +\item sine and cosine of geodetic latitude +\item magnitude of diurnal aberration vector +\item height (metres) +\item ambient temperature (degrees K) +\item pressure (mB) +\item relative humidity (0--1) +\item wavelength ($\mu$m) +\item lapse rate (degrees K per metre) +\item refraction constants A and B (radians) +\item longitude + eqn of equinoxes + ``sidereal $\Delta$ UT'' (radians) +\item local apparent sidereal time (radians) +\end{enumerate} +These may be calculated using \code{sla_AOPPA}. Note that a +\code{psGrommit} is only appropriate for a single exposure. + +Once the \code{psGrommit} has been calculated, the functions +\code{sla_OAPQK} and \code{sla_AOPQK} convert from the tangent plane +to the sky, and the sky to the tangent plane, respectively. (Note +that ``observed'' in SLALIB refers to the tangent plane, and +``apparent'' refers to the apparent position on the sky.) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{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. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{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} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsubsection{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$). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Missing and Todo} + +\tbd{define sunrise, sunset, sun position} + +\tbd{define moonrise, moonset, moon position, moon phase} + +\tbd{define planet functions} + +\tbd{clean up FITS I/O issues} + +\tbd{define Brent's method \& minimization bracketing} + +\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. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\appendix +\section{Change Log} +\input{ChangeLogADD.tex} + +\end{document} + + +\section{Modules} + +\subsection{Image Processing Modules} +\subsubsection{debias} +\subsubsection{mask} +\subsubsection{trim} +\subsubsection{flatten} +\subsubsection{sky/fringe subtract} +\subsubsection{warp} +\subsubsection{stack} +\subsubsection{difference} +\subsubsection{kernel convolution} +\subsubsection{special stack} + +\subsection{Object Detection Modules} +\subsubsection{find peaks} +\subsubsection{background } +\subsubsection{aperture photometry} +\subsubsection{get shape} + +\subsection{Miscellaneous Modules} + +\section{Analysis Stages} +\subsection{Phase 1} +\subsection{Phase 2} +\subsection{Phase 3} +\subsection{Phase 4} +\subsection{Cal 1} +\subsection{Cal 2} +\subsection{Cal 3} +\subsection{Astrom Ref} +\subsection{Photom Ref} + +\section{Architectual Components} +