Index: trunk/doc/pslib/psLibADD.tex =================================================================== --- trunk/doc/pslib/psLibADD.tex (revision 318) +++ trunk/doc/pslib/psLibADD.tex (revision 335) @@ -43,16 +43,17 @@ ``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)$ in the best -case, 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., +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 orignially mades use of the quicksort -algorithm, but now uses \tbd{what algorithm?}. 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. - -For PSLib, the sorting functions shall be implemented via +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 @@ -73,5 +74,7 @@ neighbors in the form: -\[ g_i = \sum_{n=-N}^N c_n f_i \] +\begin{equation} + g_i = \sum_{n=-N}^N c_n f_i +\end{equation} % where the values of $c_n$ determine the filter type. For boxcar @@ -80,9 +83,11 @@ the ends of the data range to reduce the value of $c_n$ as fewer input data points may be used). For Gaussian smoothing, the crucial -parameter is $\sigma$, the standard deviation. The -value of $N$ should be chosen to be large enough, $N = 5\sigma$, and -the values of $c_n$ are then just the Gaussian curve: - -\[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \] +parameter is $\sigma$, the standard deviation. The value of $N$ +should be chosen to be large enough to sample the Gaussian, $N = +5\sigma$, and the values of $c_n$ are then just the Gaussian curve: + +\begin{equation} +c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} +\end{equation} \subsection{Statistics} @@ -97,5 +102,5 @@ 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. +sample and robust estimators. \subsubsection{Sample Statistics} @@ -106,5 +111,8 @@ \paragraph{Mean} -The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \] +The mean is defined as: +\begin{equation} +\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i +\end{equation} \paragraph{Median} @@ -115,5 +123,5 @@ 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 > 10000$ elements) as the basic +samples which are large (e.g., $N > 10^4$ elements) as the basic robust median is quicker and more accurate. @@ -129,5 +137,5 @@ 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 > 10000$ elements) as the robust quartiles +which are large (e.g., $N > 10^4$ elements) as the robust quartiles are quicker and more accurate. @@ -136,21 +144,26 @@ The standard deviation of the sample is given by: -\[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \] +\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: -\[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \] +\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} \subsubsection{Clipped Statistics} -The clipped statistics are used to determine the mean and $\sigma$ 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 $\sigma$ used to exclude outliers. Typical -values for both $N$ and $k$ are 3. The algorithm is as follows: +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} @@ -162,16 +175,17 @@ \item Repeat the following N times: \begin{enumerate} - \item Exclude all values $x_{i}$ for while $|x_{i} - \bar{x}| > k \sigma$. + \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 calulcated value of $\bar{x}$ is the clipped mean. -\item The last calulcated value of $\sigma$ is the clipped standard deviation. +\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} \subsubsection{Robust Statistics} -The robust version of the statistics provide estimators of basic +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 @@ -202,5 +216,5 @@ 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$. + to $U_{\frac{1}{4}}$; this is the robust mode, $\mbox{mode}_r$. \item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 8$. \item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to @@ -242,5 +256,6 @@ diagonal matrices that satisfy the relationship $A = L U$ where $L$ is a lower-diagonal matrix of the form: -\[ L = \left( +\begin{equation} +L = \left( \begin{matrix} \alpha_{11} & 0 & 0 & 0 \\ @@ -249,8 +264,9 @@ \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: -\[ U = \left( +(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} \\ @@ -259,15 +275,19 @@ 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: -\[ \beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} \] +\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: % -\[ \alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) \] +\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} @@ -275,5 +295,7 @@ The determinant $D$ of a matrix $a_{ij}$ is calculated from the product of the diagonal elements of the LU decomposition: -\[ D = P_{i=1}{N} U_{ii} \] +\begin{equation} +D = P_{i=1}{N} U_{ii} +\end{equation} This shall be calculated using the GSL function @@ -293,11 +315,15 @@ $\beta_{ij}$, the technique of back-substitution is used to solve the equation above. First solve the equation $L y = B$: -\[ 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) \] -% +\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$: -\[ 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) \] +\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} @@ -313,18 +339,24 @@ Matrix binary arithmetic operations differ from image binary -arithmetic operations in a fundamental way: For image operations, the +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 matrix. So, for +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 multiplication correspond to: +image and matrix multiplications correspond to: % -\[ \mbox{image}_{ij} = \alpha_{ij} \times \beta_{ij} \] -\[ \mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki} \] +\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 math function \code{psMatrixOp} shall implement the matrix -version of the binary arithmetical operations for the following -operators: $+, -, \times, /$ +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} @@ -336,17 +368,32 @@ given by % -\[ T_{ij} = M_{ji} \] +\begin{equation} +T_{ij} = M_{ji} +\end{equation} where $M_{ij}$ is the matrix to be transposed. \subsubsection{Convert a matrix to a vector} -\TBD{write me} + +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{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} - -\TBD{define for use with psGetArrayPolynomial} - -\TBD{change psGetArrayPolynomial to return a psPolynomial} Given a set of data values $y_i$ with errors $\sigma_i$, related to @@ -354,15 +401,19 @@ 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$ as defined above (\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 - -\[ \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} \] +$\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: -\[ \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} \] +\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 @@ -375,10 +426,10 @@ For models in which the system of equations defined by the partial derivatives cannot be solved with the linear technique, other options -are necessary. The Levenberg-Marquardt Method (LMM; see Numerical -Recipes for a discusion) may be used for these situations. In LMM, we -make a guess at the input parameters, measure the $\chi^2$, vary the -parameters by a particular choice based on the gradient, measure the -$\chi^2$ again, and adjust the parameters and the parameter varient -based on the results. +are necessary. The Levenberg-Marquardt Method (LMM; see NR \S 15.5) +may be used for these situations. In LMM, we make a guess at the +input parameters, measure the $\chi^2$, vary the parameters by a +particular choice based on the gradient, measure the $\chi^2$ again, +and adjust the parameters and the parameter varient based on the +results. Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent @@ -389,6 +440,8 @@ $\alpha_{j,k}$ at this parameter selection as follows: -\[ \beta_k = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \] -\[ \alpha_{j,k} = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2} \] +\begin{eqnarray} +\beta_k & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \\ +\alpha_{j,k} & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2} +\end{eqnarray} % We now define the new parameter guess for $a_k$ based on the gradient @@ -396,11 +449,15 @@ follows: -\[ A_{j,k} = \alpha_{j,k} ~ (j \ne k) \] -\[ A_{j,k} = (1 + \lambda) \alpha_{j,k} ~ (j = k) \] +\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: -\[ A_{j,k} \alpha^\prime_k = \beta_j \] +\begin{equation} +A_{j,k} a^\prime_k = \beta_j +\end{equation} % -where $\alpha^\prime_k$ represents our new attempt at a parameter +where $a^\prime_k$ represents our new attempt at a parameter guess. We use this parameter set to calculate $\chi^2$. If the new value of $\chi^2$ is lower than the previous guess, we accept this new @@ -437,4 +494,6 @@ \end{center} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \subsection{Polynomials} @@ -446,5 +505,5 @@ \begin{itemize} \item They are bounded on $-1 < x < 1$, with the maxima and minima -over this range being 1 and -1, respectively; +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. @@ -477,6 +536,9 @@ Multi-dimensional polynomials shall be composed of multiplications of -1D Chebyshev polynomials. - +1D Chebyshev polynomials, with the coefficients stored in tensors of +the appropriate rank. + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{(Fast) Fourier Transforms} @@ -500,5 +562,6 @@ 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. +initialisation of the PSLib FFT functions (\code{psFFTInit()}) and +saved at the conclusion (\code{psFFTDone()}). \subsubsection{Function mapping} @@ -515,5 +578,5 @@ These plans should be formulated using the \code{FFTW_ESTIMATE} flag, -which will allow FFTW to default to a short planning time if the +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. @@ -521,11 +584,11 @@ \subsubsection{More Complicated Functions} -The \code{psFFTFilter} and \code{psFFTFilterComplex} functions provide +The \code{psFFTFilter()} and \code{psFFTFilterComplex()} functions provide the means to apply a filter (purely real and imaginary multipliers, respectively) to the data in the Fourier plane. If the filter -function specified for \code{psFFTFilter} returns a real value, $r$, +function specified for \code{psFFTFilter()} returns a real value, $r$, then the corresponding value in the Fourier plane should be multiplied by $r$. If the real and imaginary filter functions specified for -\code{psFFTFilterComplex} return the values $r$ and $s$, respectively, +\code{psFFTFilterComplex()} return the values $r$ and $s$, respectively, then the corresponding value in the Fourier plane should be multiplied by the complex number $r + si$. @@ -535,10 +598,10 @@ Fourier transform is multiplied by the complex conjugate of the second Fourier transform to yield the Fourier transform of the -cross-correlation (NR 13.2). In the latter, the two Fourier +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 13.1). +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 13.4): +the elements of the power spectrum are (NR \S 13.4): \begin{eqnarray} P_0 & = & \left| C_0 \right|^2 / N^2 \\ @@ -550,4 +613,6 @@ 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} @@ -557,4 +622,6 @@ \href{http://star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{SLALIB Positional Astronomy Library}. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsubsection{Celestial Coordinate Conversions} @@ -571,4 +638,6 @@ \end{tabular} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \subsubsection{Projections} @@ -584,6 +653,6 @@ \begin{eqnarray} -x & = & R \sin \phi \\ -y & = & -R \cos \phi +x & = & R \sin (\phi) \\ +y & = & -R \cos (\phi) \end{eqnarray} @@ -600,5 +669,5 @@ \begin{eqnarray} -R & = & \cot \theta 180^\circ/\pi \\ +R & = & \cot (\theta) 180^\circ/\pi \\ \theta & = & \arctan (180^\circ/(\pi R)) \end{eqnarray} @@ -609,5 +678,5 @@ \begin{eqnarray} -R & = & \cos \theta 180^\circ/\pi \\ +R & = & \cos (\theta) 180^\circ/\pi \\ \theta & = & \arccos (\pi R / 180^\circ) \end{eqnarray} @@ -637,7 +706,7 @@ \begin{eqnarray} -x & = & 2 \alpha \cos \theta \sin (\phi/2) \\ +x & = & 2 \alpha \cos (\theta) \sin (\phi/2) \\ y & = & \alpha \sin \theta \\ -{\rm where}\hspace{1cm} \alpha^{-1} & \equiv & (180^\circ/\pi) \sqrt{\left(1 + \cos \theta \cos (\phi/2) \right) / 2} +{\rm where}\hspace{1cm} \alpha^{-1} & \equiv & (180^\circ/\pi) \sqrt{\left(1 + \cos (\theta) \cos (\phi/2) \right) / 2} \end{eqnarray} @@ -649,4 +718,109 @@ {\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x\pi/720^\circ)^2 - (y\pi/360^\circ)^2} \end{eqnarray} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Tangent Plane to Sky} + +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.) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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} + +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. + +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} + +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{Offsets} + + +