Index: /trunk/doc/pslib/ChangeLogADD.tex =================================================================== --- /trunk/doc/pslib/ChangeLogADD.tex (revision 3717) +++ /trunk/doc/pslib/ChangeLogADD.tex (revision 3718) @@ -29,5 +29,5 @@ \end{itemize} -\subsection{Changes from version 09 (14 February 2005) to version 10 (29 March 2005} +\subsection{Changes from version 09 (14 February 2005) to version 10 (29 March 2005)} \begin{itemize} @@ -48,4 +48,11 @@ \item promoted all PSLib subsections to sections \end{itemize} +\end{itemize} +\subsection{Changes from version 10 (29 March 2005) to version 11 (19 April 2005)} + +\begin{itemize} +\item defined earth orientation algorithms +\item provided quaternions for rotations +\item moved some sections to reflect order in SDRS (matrix, fftw) \end{itemize} Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 3717) +++ /trunk/doc/pslib/psLibADD.tex (revision 3718) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.70 2005-03-30 21:14:48 eugene Exp $ +%%% $Id: psLibADD.tex,v 1.71 2005-04-19 23:40:47 eugene Exp $ \documentclass[panstarrs]{panstarrs} @@ -14,8 +14,8 @@ \project{Pan-STARRS Image Processing Pipeline} \organization{Institute for Astronomy} -\version{10} +\version{11} \docnumber{PSDC-430-006} -\setcounter{tocdepth}{5} % lowest level to be included in toc +% \setcounter{tocdepth}{5} % lowest level to be included in toc \newcommand\citealt{} @@ -41,4 +41,5 @@ 09 & 2005 Feb 14 & Frozen for Cycle 5 \\ \hline 10 & 2005 Mar 21 & Draft for Cycle 6 \\ \hline +11 & 2005 Apr 19 & Frozen for Cycle 6 \\ \hline \RevisionsEnd @@ -347,154 +348,49 @@ Note that the total adds to one --- the number of values added. -\subsection{Matrix Operations} - -In this section, we define the linear algebra operations performed on -matrices. We have defined APIs to implement the following matrix -functions: - +\subsection{Polynomials} +\label{sec:polynomials} + +We will employ Chebyshev polynomials (NR \S 5.8) to approximate functions: +\begin{equation} +f(x) = \sum_{i=0}^{n} c_i T_i(x) +\end{equation} +These have some desirable features: \begin{itemize} -\item 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. +\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} -Many of these operations are implemented using LU decomposition. We -define LU decomposition in the following paragraph. Implementation of -LU decomposition shall make use of the GSL function -\code{gsl_linalg_LU_decomp}. - -\subsubsection{LU Decomposition} -\label{LUdecomp} - -We wish to decompose the matrix $A$ with elements $a_{ij}$ into -diagonal matrices that satisfy the relationship $A = L U$ where $L$ is -a lower-diagonal matrix of the form: -\begin{equation} -L = \left( -\begin{matrix} -\alpha_{11} & 0 & 0 & 0 \\ -\alpha_{21} & \alpha_{22} & 0 & 0 \\ -\alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\ -\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\ -\end{matrix} \right) -\end{equation} -% -(where all diagonal values $\alpha_{ii} = 1$) and $U$ is an upper-diagonal matrix of the form: -\begin{equation} -U = \left( -\begin{matrix} -\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\ -0 & \beta_{22} & \beta_{23} & \beta_{24} \\ -0 & 0 & \beta_{33} & \beta_{34} \\ -0 & 0 & 0 & \beta_{44} \\ -\end{matrix} \right) -\end{equation} - -We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following -this procedure. First, $\alpha_{ii} = 1$. For all values of $j = 1, -2, \dots, N$, follow the next steps: For each value of $i = 1, 2, -\dots, j$, solve for $\beta_{ij}$ using the relationship: -\begin{equation} -\beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} -\end{equation} -For the values of $i = j+1, j+2, \dots, N$, solve for the values of -$\alpha_{ij}$ with the following: -% -\begin{equation} -\alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) -\end{equation} - -\subsubsection{Calculate a matrix determinant} - -The determinant $D$ of a matrix $a_{ij}$ is calculated from the -product of the diagonal elements of the LU decomposition: -\begin{equation} -D = P_{i=1}{N} U_{ii} -\end{equation} - -This shall be calculated using the GSL function -\code{gsl_linalg_LU_det}. If the matrix is large or the magnitude of -the elements is large, the determinant may overflow the data type. In -these cases, the (natural) logarithm of the determinant may be -calculated instead (as the sum of the logarithms of the diagonal -elements). In this case, the GSL function \code{gsl_linalg_LU_lndet} -shall be used. - -\subsubsection{Solving a Linear Equation} - -The LU decomposition of a matrix may be used to solve the -matrix equation $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ - -Given the LU decomposition of a matrix into $\alpha_{ij}$ and -$\beta_{ij}$, the technique of back-substitution is used to solve the -equation above. First solve the equation $L y = B$: -\begin{eqnarray} -y_1 & = & \frac{b_1}{\alpha_{11}} \\ -y_i & = & \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right) -\end{eqnarray} - -The values of $y$ may be used to solve for $x_i$ using the -relationship $U x = b$: -\begin{eqnarray} -x_N & = & \frac{y_N}{\beta_{NN}} \\ -x_i & = & \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) -\end{eqnarray} - -\subsubsection{Invert a matrix} - -Inversion of a matrix using the LU decomposition is performed by -performing back-subsitution to solve a series of linear equations of -the form $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ where the values -of $B_j$ represent the columns of the identity matrix. The solution -vectors $x_j$ represent the columns of the inverse matrix. This -operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}. - -\subsubsection{Perform matrix addition, subtraction and multiplication} - -Matrix binary arithmetic operations differ from image binary -arithmetic operations in a fundamental way: For image operations, the -resulting pixel value is determined by performing the operation on the -two corresponding input operand pixels. For matrix operations, the -resulting element value results from the operation on the -corresponding pair of row and colum from the input matrices. So, for -example, given the input matrices $\alpha_{ij}$ and $\beta_{ij}$, the -image and matrix multiplications correspond to: -% -\begin{eqnarray} -\mbox{image}_{ij} & = & \alpha_{ij} \times \beta_{ij} \\ -\mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki} -\end{eqnarray} -% -The matrix and image operations for addition and subtraction are -identical: the operation is performed on the corresponding elements in -each matrix. Matrix division is not defined (to divide, the user -should first invert the matrix, and then multiply). The matrix math -function \code{psMatrixOp} shall implement the matrix version of the -binary arithmetical operations for the following operators: $+, -, -\times$. - -\subsubsection{Transpose a matrix} - -The transpose of a matrix is simply the reorganization of the matrix -elements by 'flipping' the matrix along a diagonal. For non-square -matrices with dimensions $N \times M$, the resulting matrix has -dimensions $M \times N$. The element of the output matrix $T_{ij}$ is -given by -% -\begin{equation} -T_{ij} = M_{ji} -\end{equation} -where $M_{ij}$ is the matrix to be transposed. - -\subsubsection{Convert a matrix to a vector} - -Matrix-to-vector conversion is only defined for a matrix that has a -size of one in at least one dimension. In that case, the elements -should be extracted, and placed in a vector, with care that the -\code{psType} is defined correctly --- a $1\times N$ matrix is -converted to a \code{PS_DIMEN_VECTOR}-type vector, while a $N\times 1$ -matrix is converted to a \code{PS_DIMEN_TRANV}-type vector. +The first few Chebyshev polynomials are: +\begin{eqnarray} +T_0(x) & = & 1 \\ +T_1(x) & = & x \\ +T_2(x) & = & 2x^2 - 1 \\ +T_3(x) & = & 4x^3 - 3x \\ +T_4(x) & = & 8x^4 - 8x^2 + 1 \\ +\end{eqnarray} +Chebyshev polynomials follow the recurrence relation: +\begin{equation} +T_{n+1} = 2xT_n - T_{n-1} +\end{equation} + +Practically, Chebyshev polynomials should be evaluated using Clenshaw's recurrence +formula (NR \S 5.5): +\begin{eqnarray} +d_j & = & 2xd_{j+1} - d_{j+2} + c_j \\ +f(x) & = & x*d_1 - d_2 + 1/2 c_0 \\ +\end{eqnarray} + +It shall be the responsibility of the user to convert the domain into the range +$-1 < x < 1$. + +\subsubsection{Multi-dimensional polynomials} + +Multi-dimensional polynomials shall be composed of multiplications of +1D Chebyshev polynomials, with the coefficients stored in tensors of +the appropriate rank. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Fitting} @@ -670,111 +566,4 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Polynomials} -\label{sec:polynomials} - -We will employ Chebyshev polynomials (NR \S 5.8) to approximate functions: -\begin{equation} -f(x) = \sum_{i=0}^{n} c_i T_i(x) -\end{equation} -These have some desirable features: -\begin{itemize} -\item They are bounded on $-1 < x < 1$, with the maxima and minima -over this range being 1 and $-1$, respectively; -\item Truncation of the higher-order terms leaves one with the most accurate -lower-order polynomial representation of the desired function. -\end{itemize} - -The first few Chebyshev polynomials are: -\begin{eqnarray} -T_0(x) & = & 1 \\ -T_1(x) & = & x \\ -T_2(x) & = & 2x^2 - 1 \\ -T_3(x) & = & 4x^3 - 3x \\ -T_4(x) & = & 8x^4 - 8x^2 + 1 \\ -\end{eqnarray} -Chebyshev polynomials follow the recurrence relation: -\begin{equation} -T_{n+1} = 2xT_n - T_{n-1} -\end{equation} - -Practically, Chebyshev polynomials should be evaluated using Clenshaw's recurrence -formula (NR \S 5.5): -\begin{eqnarray} -d_j & = & 2xd_{j+1} - d_{j+2} + c_j \\ -f(x) & = & x*d_1 - d_2 + 1/2 c_0 \\ -\end{eqnarray} - -It shall be the responsibility of the user to convert the domain into the range -$-1 < x < 1$. - -\subsubsection{Multi-dimensional polynomials} - -Multi-dimensional polynomials shall be composed of multiplications of -1D Chebyshev polynomials, with the coefficients stored in tensors of -the appropriate rank. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsection{(Fast) Fourier Transforms} - -(Fast) Fourier Transforms (FFTs) shall be implemented using the -\href{www.fftw.org}{{\em Fastest Fourier Transform in the West} (FFTW) -library}. - -\subsubsection{FFTW Plans} - -FFTW requires the user to create a ``plan'' for each transform size, -the time to create which is a function of the desired speed of the -transform --- faster transforms for a given size (and machine) may be -performed if more time is spent testing plans. - -In the \PS{} IPP, we will want to perform FFTs on images of common -sizes (e.g.\ $512 \times 512$) regularly. This means that we would -gain from determining an FFTW plan for each of these common sizes. -FFTW provides a binary, \code{fftw-wisdom} which may be used to -generate and save ``wisdom''. The location of the \code{wisdom} file -will be specified as a configuration variable for the IPP (defaulting -to \code{/etc/fftw/wisdom}). The \code{wisdom} should be read in upon -initialisation of the PSLib FFT functions and saved at the conclusion. - -\subsubsection{Function mapping} - -The forward and reverse transforms call the corresponding -FFTW function to plan the transform: - -\begin{tabular}{ll} - PSLib function & Major FFTW call \\ \hline -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - \code{psFFTForward()} & \code{fftw_plan_dft_r2c_2d()} \\ - \code{psFFTReverse()} & \code{fftw_plan_dft_c2r_2d()} \\ -\end{tabular} - -These plans should be formulated using the \code{FFTW_ESTIMATE} flag, -which will allow FFTW to default to a minimal planning time if the -wisdom has not been loaded. Transforms should be performed out of -place to avoid the need to pad the input array to hold the output. - -\subsubsection{More Complicated Functions} - -\code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve -multiplication of two Fourier transforms. In the former, the first -Fourier transform is multiplied by the complex conjugate of the second -Fourier transform to yield the Fourier transform of the -cross-correlation (NR \S 13.2). In the latter, the two Fourier -transforms are multiplied directly to yield the Fourier transform of -the convolution (NR \S 13.1). - -If the elements of the discrete Fourier transform are $C_k$, then the -the elements of the power spectrum are (NR \S 13.4): -\begin{eqnarray} -P_0 & = & \left| C_0 \right|^2 / N^2 \\ -P_j & = & \left( \left| C_j \right|^2 + \left| C_{N-j} \right|^2 \right)/ N^2 \\ -P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2 \\ -\end{eqnarray} -where $j = 1, 2, \ldots, (N/2 - 1)$. - -Note that we leave the issue of ``windowing'' the data up to the -caller, and choose to normalise by $1/N^2$. - \subsection{Image Manipulations} @@ -971,11 +760,218 @@ \pagebreak + +\subsection{Matrix Operations} + +In this section, we define the linear algebra operations performed on +matrices. We have defined APIs to implement the following matrix +functions: + +\begin{itemize} +\item Invert a matrix; +\item Calculate a matrix determinant; +\item Perform matrix addition, subtraction and multiplication; +\item Transpose a matrix; and +\item Convert a matrix to a vector. +\end{itemize} + +Many of these operations are implemented using LU decomposition. We +define LU decomposition in the following paragraph. Implementation of +LU decomposition shall make use of the GSL function +\code{gsl_linalg_LU_decomp}. + +\subsubsection{LU Decomposition} +\label{LUdecomp} + +We wish to decompose the matrix $A$ with elements $a_{ij}$ into +diagonal matrices that satisfy the relationship $A = L U$ where $L$ is +a lower-diagonal matrix of the form: +\begin{equation} +L = \left( +\begin{matrix} +\alpha_{11} & 0 & 0 & 0 \\ +\alpha_{21} & \alpha_{22} & 0 & 0 \\ +\alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\ +\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\ +\end{matrix} \right) +\end{equation} +% +(where all diagonal values $\alpha_{ii} = 1$) and $U$ is an upper-diagonal matrix of the form: +\begin{equation} +U = \left( +\begin{matrix} +\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\ +0 & \beta_{22} & \beta_{23} & \beta_{24} \\ +0 & 0 & \beta_{33} & \beta_{34} \\ +0 & 0 & 0 & \beta_{44} \\ +\end{matrix} \right) +\end{equation} + +We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following +this procedure. First, $\alpha_{ii} = 1$. For all values of $j = 1, +2, \dots, N$, follow the next steps: For each value of $i = 1, 2, +\dots, j$, solve for $\beta_{ij}$ using the relationship: +\begin{equation} +\beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} +\end{equation} +For the values of $i = j+1, j+2, \dots, N$, solve for the values of +$\alpha_{ij}$ with the following: +% +\begin{equation} +\alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) +\end{equation} + +\subsubsection{Calculate a matrix determinant} + +The determinant $D$ of a matrix $a_{ij}$ is calculated from the +product of the diagonal elements of the LU decomposition: +\begin{equation} +D = P_{i=1}{N} U_{ii} +\end{equation} + +This shall be calculated using the GSL function +\code{gsl_linalg_LU_det}. If the matrix is large or the magnitude of +the elements is large, the determinant may overflow the data type. In +these cases, the (natural) logarithm of the determinant may be +calculated instead (as the sum of the logarithms of the diagonal +elements). In this case, the GSL function \code{gsl_linalg_LU_lndet} +shall be used. + +\subsubsection{Solving a Linear Equation} + +The LU decomposition of a matrix may be used to solve the +matrix equation $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ + +Given the LU decomposition of a matrix into $\alpha_{ij}$ and +$\beta_{ij}$, the technique of back-substitution is used to solve the +equation above. First solve the equation $L y = B$: +\begin{eqnarray} +y_1 & = & \frac{b_1}{\alpha_{11}} \\ +y_i & = & \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right) +\end{eqnarray} + +The values of $y$ may be used to solve for $x_i$ using the +relationship $U x = b$: +\begin{eqnarray} +x_N & = & \frac{y_N}{\beta_{NN}} \\ +x_i & = & \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) +\end{eqnarray} + +\subsubsection{Invert a matrix} + +Inversion of a matrix using the LU decomposition is performed by +performing back-subsitution to solve a series of linear equations of +the form $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ where the values +of $B_j$ represent the columns of the identity matrix. The solution +vectors $x_j$ represent the columns of the inverse matrix. This +operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}. + +\subsubsection{Perform matrix addition, subtraction and multiplication} + +Matrix binary arithmetic operations differ from image binary +arithmetic operations in a fundamental way: For image operations, the +resulting pixel value is determined by performing the operation on the +two corresponding input operand pixels. For matrix operations, the +resulting element value results from the operation on the +corresponding pair of row and colum from the input matrices. So, for +example, given the input matrices $\alpha_{ij}$ and $\beta_{ij}$, the +image and matrix multiplications correspond to: +% +\begin{eqnarray} +\mbox{image}_{ij} & = & \alpha_{ij} \times \beta_{ij} \\ +\mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki} +\end{eqnarray} +% +The matrix and image operations for addition and subtraction are +identical: the operation is performed on the corresponding elements in +each matrix. Matrix division is not defined (to divide, the user +should first invert the matrix, and then multiply). The matrix math +function \code{psMatrixOp} shall implement the matrix version of the +binary arithmetical operations for the following operators: $+, -, +\times$. + +\subsubsection{Transpose a matrix} + +The transpose of a matrix is simply the reorganization of the matrix +elements by 'flipping' the matrix along a diagonal. For non-square +matrices with dimensions $N \times M$, the resulting matrix has +dimensions $M \times N$. The element of the output matrix $T_{ij}$ is +given by +% +\begin{equation} +T_{ij} = M_{ji} +\end{equation} +where $M_{ij}$ is the matrix to be transposed. + +\subsubsection{Convert a matrix to a vector} + +Matrix-to-vector conversion is only defined for a matrix that has a +size of one in at least one dimension. In that case, the elements +should be extracted, and placed in a vector, with care that the +\code{psType} is defined correctly --- a $1\times N$ matrix is +converted to a \code{PS_DIMEN_VECTOR}-type vector, while a $N\times 1$ +matrix is converted to a \code{PS_DIMEN_TRANV}-type vector. + +\subsection{(Fast) Fourier Transforms} + +(Fast) Fourier Transforms (FFTs) shall be implemented using the +\href{www.fftw.org}{{\em Fastest Fourier Transform in the West} (FFTW) +library}. + +\subsubsection{FFTW Plans} + +FFTW requires the user to create a ``plan'' for each transform size, +the time to create which is a function of the desired speed of the +transform --- faster transforms for a given size (and machine) may be +performed if more time is spent testing plans. + +In the \PS{} IPP, we will want to perform FFTs on images of common +sizes (e.g.\ $512 \times 512$) regularly. This means that we would +gain from determining an FFTW plan for each of these common sizes. +FFTW provides a binary, \code{fftw-wisdom} which may be used to +generate and save ``wisdom''. The location of the \code{wisdom} file +will be specified as a configuration variable for the IPP (defaulting +to \code{/etc/fftw/wisdom}). The \code{wisdom} should be read in upon +initialisation of the PSLib FFT functions and saved at the conclusion. + +\subsubsection{Function mapping} + +The forward and reverse transforms call the corresponding +FFTW function to plan the transform: + +\begin{tabular}{ll} + PSLib function & Major FFTW call \\ \hline +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \code{psFFTForward()} & \code{fftw_plan_dft_r2c_2d()} \\ + \code{psFFTReverse()} & \code{fftw_plan_dft_c2r_2d()} \\ +\end{tabular} + +These plans should be formulated using the \code{FFTW_ESTIMATE} flag, +which will allow FFTW to default to a minimal planning time if the +wisdom has not been loaded. Transforms should be performed out of +place to avoid the need to pad the input array to hold the output. + +\subsubsection{More Complicated Functions} + +\code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve +multiplication of two Fourier transforms. In the former, the first +Fourier transform is multiplied by the complex conjugate of the second +Fourier transform to yield the Fourier transform of the +cross-correlation (NR \S 13.2). In the latter, the two Fourier +transforms are multiplied directly to yield the Fourier transform of +the convolution (NR \S 13.1). + +If the elements of the discrete Fourier transform are $C_k$, then the +the elements of the power spectrum are (NR \S 13.4): +\begin{eqnarray} +P_0 & = & \left| C_0 \right|^2 / N^2 \\ +P_j & = & \left( \left| C_j \right|^2 + \left| C_{N-j} \right|^2 \right)/ N^2 \\ +P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2 \\ +\end{eqnarray} +where $j = 1, 2, \ldots, (N/2 - 1)$. + +Note that we leave the issue of ``windowing'' the data up to the +caller, and choose to normalise by $1/N^2$. + \section{PSLib Astronomy Utilities} - -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} \subsection{Time} @@ -1015,5 +1011,5 @@ Leap-seconds are declared by the International Earth Rotation and Reference Systems Service (IERS)\footnote{IERS website - http://www.iers.org/}. -leap-seconds only occur in the UTC time system and cannot be accurately +Leap-seconds only occur in the UTC time system and cannot be accurately predicted due to variations in the Earth's rotational period. To determine the number of leap-second in a given UTC date a table of leap-seconds as annouced by @@ -1021,8 +1017,14 @@ leap-second occurs. -For ease of conversion, UTC should be represented as the number of seconds -since the UNIX epoch of ``1970-01-01T00:00:00Z'', non-inclusive of leap-seconds. -\emph{Times will always be expressed in the 'UTC timezone'. Use of the local -timezone is forbidden.} +For ease of conversion, UTC should be represented as the number of +seconds since the UNIX epoch of ``1970-01-01T00:00:00Z'', +non-inclusive of leap-seconds. \tbd{what does this statement +actually mean? what is the source of time (gettimeofday? let's be +explicit here} + +\tbd{Times will always be expressed in the 'UTC timezone'. Use of +the local timezone is forbidden. -- this makes no sense given that we +define LST. In any case, this statement, or somethign equivalent, +belongs in the SDRS not the ADD} \subsubsection{International Atomic Time (TAI)} @@ -1040,5 +1042,6 @@ For ease of conversion, TAI should be represented as the number of -seconds since the UNIX epoch of ``1970-01-01T00:00:00Z''. +seconds since the UNIX epoch of ``1970-01-01T00:00:00Z''. \tbd{what +does this statement actually mean?} \subsubsection{Leap-seconds} @@ -1076,7 +1079,25 @@ 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} +file will be made configurable. This data is available from +USNO\footnote{ftp://maia.usno.navy.mil/ser7/tai-utc.dat}. + +\subsubsection{Universal Time (UT1)} +\label{sec:ut1} + +UT1 is directly tied to the rotation of the Earth. Historically, time +has been measured with respect to the rising and setting of the +Sun. However, in the modern era of atomic clocks, the rotation of the +Earth makes for a highly unstable time standard. Tidal effects, +changes in the angular momentum of the atmosphere, seasonal changes in +the polar ice caps, movement within the Earth's core, and other +effects all cause measurable changes in the Earth's rotation on a +daily basis. However, UT1 is still vitally important for determining +the orientation of the Earth with respect to the sky. UT1 is +calculated by applying the value UT1-UTC to the value of UTC. UT1 is +continuously measured by the International Earth Rotation Service, and +tabulated values of the offset of UT1 from UTC are published at +regular intervals, along with predicted future values. The process of +calculating the UT1-UTC offsets are discussed in the section on Earth +Orientation (Section~\ref{sec:ut1}). \subsubsection{Gregorian dates to seconds} @@ -1209,47 +1230,4 @@ -\subsubsection{Universal Time (UT1)} -\label{sec:ut1} - -UT1 is directly tied to the rotation of the Earth. Historically, time has been -measured with respect to the rising and setting of the Sun. However, in the -modern era of atomic clocks, the rotation of the Earth makes for a highly -unstable time standard. Tidal effects, changes in the angular momentum of the -atmosphere, seasonal changes in the polar ice caps, movement within the Earth's -core, and other effects all cause measurable changes in the Earth's rotation on -a daily basis. However, UT1 is still vitally important for determining the -orientation of the Earth with respect to the sky. - -Since 2003-01-01, UT1 has been defined as directly proportional to the Earth -Rotation Angle (see IERS Technical Note 32\footnote{IERS Technical Note 32 - -http://maia.usno.navy.mil/conv2003.html}). Previous to that date, a different -definition was in effect (see IERS Technical Note 21\footnote{IERS Technical -Note 21 - http://maia.usno.navy.mil/conventions.html}). We will always use the -post-2003 definition. - -UT1 is continuously measured by the International Earth Rotation Service, and -tabulated values of the offset of UT1 from UTC are published at regular -intervals, along with predicted future values. IERS Bulletin A gives ``rapid -response'' values necessary for real-time and near real-time data analysis -(such as Pan-STARRS Otis and IPP subsystems). Bulletin B gives the results of a -final, definitive data reduction. An amalgam of Bulletin A and B values is -published daily on the IERS website\footnote{IERS Bulletin A \& B - -http://maia.usno.navy.mil/ser7/finals2000A.daily} along with a desciption of -the format\footnote{IERS finals2000A.daily table format - -http://maia.usno.navy.mil/ser7/readme.finals2000A}. - -The UT1 offsets should be interpolated using the prescription of IERS Gazette -\#13\footnote{IERS Gazette \#13 - http://maia.usno.navy.mil/iers-gaz13}. This -entails using a third order polynomial to interpolate the table values, and -then applying a model for diurnal and semi-diurnal fluctuations due to tidal -effects. An example implimentation\footnote{interp.f - -ftp://maia.usno.navy.mil/dist/interp.f} written in Fortran is available. - -\verbatiminput{interp.f} - -The Polar motion X, and Y coordinates are also important for determining the -orientation of the Earth with respect to the sky. This is also given in the -IERS publications references above, and should be interpolated in the same way. - \subsubsection{Julian Date and Modified Julian Date} @@ -1369,22 +1347,4 @@ Gives $GMST00$ in seconds. - -\subsubsection{Longitude} - -Longitudes are often expressed in the form of decimal degrees while the -algorithm for calculating GMST outputs seconds. - -\begin{equation} -1\degree = 240s -\end{equation} - -\subsubsection{Polar Motion} -\tbd{move this to Earth Motion section} - -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. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{2D transformations} @@ -1442,5 +1402,21 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Spherical Rotations with Quaternions} +\subsection{Spherical Rotations} + +Spherical rotations may be implemented with a variety of mathematical +methods. A single rotation may be specified by three angles which +describe rotations about the principal axes. Figure~\ref{rotations} +shows the rotation angles for an arbitrary 3-D rotation. A single +rotation may be represented by defining this set of Euler angles. As +an alternative, the quaternion is a convenient construct with which to +represent a single rotation, and it has useful mathematical properties +when applying those rotations. + +\begin{figure} +\begin{center} +\psfig{file=rotations.ps} +\caption{Definition of the rotation angles\label{rotations}} +\end{center} +\end{figure} \subsubsection{Quaternion Construction} @@ -1459,5 +1435,5 @@ Note that the sine and cosine are taken of the half angle of the rotation. Note also that this implies that the quaternion components -are normalized such that $|\bar{q}| \def q_0^2 + q_1^2 + q_2^2 + q_3^2 +are normalized such that $|\bar{q}| \equiv q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$. @@ -1472,22 +1448,22 @@ \subsubsection{Combining Two Rotations} -Given two quaternions $\bar{p1}$ and $\bar{p2}$, there is a third +Given two quaternions $\bar{a}$ and $\bar{b}$, there is a third quaternion, $\bar{p}$, which represents the result of first applying -$\bar{p1}$, and then $\bar{p2}$. The components of $\bar{p}$ are given +$\bar{a}$, and then $\bar{b}$. The components of $\bar{p}$ are given by: -\begin{verbatim} -p_0 & = & p2_3 p1_0 + p2_2 p1_1 - p2_1 p1_2 + p2_0 p1_3 \\ -p_1 & = & -p2_2 p1_0 + p2_3 p1_1 + p2_0 p1_2 + p2_1 p1_3 \\ -p_2 & = & p2_1 p1_0 - p2_0 p1_1 + p2_3 p1_2 + p2_2 p1_3 \\ -p_3 & = & -p2_0 p1_0 - p2_1 p1_1 - p2_2 p1_2 + p2_3 p1_3 \\ -\end{verbatim} +\begin{eqnarray} +p_0 & = & b_3 a_0 + b_2 a_1 - b_1 a_2 + b_0 a_3 \\ +p_1 & = & -b_2 a_0 + b_3 a_1 + b_0 a_2 + b_1 a_3 \\ +p_2 & = & b_1 a_0 - b_0 a_1 + b_3 a_2 + b_2 a_3 \\ +p_3 & = & -b_0 a_0 - b_1 a_1 - b_2 a_2 + b_3 a_3 \\ +\end{eqnarray} \subsubsection{Rotating a Vector} You may rotate a unit vector by first constructing a quaternion -$\bar{p2}$, whose first three components are the components of the +$\bar{b}$, whose first three components are the components of the unit vector, and whose fourth component is zero. To rotate this vector -by a quaternion $\bar{p1}$, you apply the formula above for combining +by a quaternion $\bar{a}$, you apply the formula above for combining two quaternions. The rotated vector is found in the first three components of the resulting quaternion, $\bar{p}$. @@ -1503,11 +1479,8 @@ rot_{y,y} & = & -q_0 q_0 + q_1 q_1 - q_2 q_2 + q_3 q_3 \\ rot_{z,z} & = & -q_0 q_0 - q_1 q_1 + q_2 q_2 + q_3 q_3 \\ - rot_{x,y} & = & 2 (q_0 q_1 + q_2 q_3) \\ rot_{y,x} & = & 2 (q_0 q_1 - q_2 q_3) \\ - rot_{x,z} & = & 2 (q_0 q_2 - q_1 q_3) \\ rot_{z,x} & = & 2 (q_0 q_2 + q_1 q_3) \\ - rot_{y,z} & = & 2 (q_1 q_2 + q_0 q_3) \\ rot_{z,y} & = & 2 (q_1 q_2 - q_0 q_3) \\ @@ -1576,5 +1549,10 @@ \end{itemize} Note that $\theta_p$, the latitude of the source system pole in the -target system, is equal to $\delta_p$ by symmetry. +target system, is equal to $\delta_p$ by symmetry. Transformations +between arbitrary systems are specified by determining the appropriate +values of $\alpha_p$, $\delta_p$, and $\phi_p$, and then constructing +the quarternion for this transformation. + +\tbd{can we drop this, since we do this with the quaternion?} The relevant trigonometric relationships are: @@ -1627,5 +1605,6 @@ \subsubsection{Precession} -The appropriate values, from Elixir, are: +Approximate precession, good for modest sub-arcsecond precision, may +be rapidly calculated using the following rotation angles: \begin{eqnarray} \alpha_p & = & 90^\circ - (0^\circ.6406161\, T + 0^\circ.0000839\, T^2 + 0^\circ.0000050\, T^3) \\ @@ -1660,40 +1639,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{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.) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsection{Sky to Tangent Plane (II)} +\subsection{Sky to Tangent Plane} This section describes the transformation between celestial coordinates @@ -1725,11 +1669,13 @@ \begin{figure} -\psfig{file=transforms.ps} -\caption{Coordinates systems and the transformations between them} - +\begin{center} +\psfig{file=earthrot.ps} +\caption{Coordinates systems and the transformations between them\label{earthrot}} +\end{center} \end{figure} -Figure X shows the transformation steps and intermediate coordinate systems -between celestial and local terrestrial coordinate systems. The intermediate -coordinate systems are defined below. + +Figure~\ref{earthrot} shows the transformation steps and intermediate +coordinate systems between celestial and local terrestrial coordinate +systems. The intermediate coordinate systems are defined below. \paragraph{ICRS} @@ -1786,41 +1732,39 @@ \paragraph{Gravitational Deflection} -The Sun's gravity bends the path of light rays which pass near it. -To first order, a light ray is deflected by an angle of $4GM/c^2r_0$ radians, -where $G$ is the gravitational constant, -$M$ is the mass of the Sun, -$c$ is the speed of light, and -$r_0$ is the point of closest approach to the light ray to the Sun. -To the same order this is equal to the impact parameter - i.e. the point -of closest approach if the light ray were not deflected. Note that -$r_0/d = \tan(\theta)$, where $d$ is the distance from the Earth -to the Sun, and $\theta$ -is the angular separation of the star from the center of the Sun. - -There is a maximum deflection of 1.75 arc seconds if we set -$r_0$ to the radius of the sun. -Since the Sun bends light rays toward it, a star appears shifted away from the sun in the sky. +The Sun's gravity bends the path of light rays which pass near it. To +first order, a light ray is deflected by an angle of $4GM/c^2r_0$ +radians, where $G$ is the gravitational constant, $M$ is the mass of +the Sun, $c$ is the speed of light, and $r_0$ is the point of closest +approach to the light ray to the Sun. To the same order this is equal +to the impact parameter - i.e. the point of closest approach if the +light ray were not deflected. Note that $r_0/d = \tan(\theta)$, where +$d$ is the distance from the Earth to the Sun, and $\theta$ is the +angular separation of the star from the center of the Sun. + +There is a maximum deflection of 1.75 arc seconds if we set $r_0$ to +the radius of the sun. Since the Sun bends light rays toward it, a +star appears shifted away from the sun in the sky. \paragraph{Aberration} -Aberration is the apparent change in direction of a ray of light in the -reference frame of a moving observer. Traditionally the aberration +Aberration is the apparent change in direction of a ray of light in +the reference frame of a moving observer. Traditionally the aberration calculation has been done with a linear expansion of the full relativistic expression, often neglecting all but the linear term in -$v/c$, since the relativistic terms are on the order of a miliarcsecond. -However, the full relativistic expression poses no challenge for modern -computers, so psLib will use the following procedure to calculate aberration. +$v/c$, since the relativistic terms are on the order of a +miliarcsecond. However, the full relativistic expression poses no +challenge for modern computers, so psLib will use the following +procedure to calculate aberration. Suppose an observer has a velocity $\beta\hat{\beta}$, with respect to the Solar System barycenter, where $\beta$ is in units of the speed of -light, and $\hat{\beta}$ is a unit vector. Suppose also that the unit vector -$\hat{r}$ points toward a star in the barycenter frame of reference -(i.e. the ``actual'' position). -and $\hat{r}'$ gives the direction of the star in the observer's frame, -(i.e. the apparent position). - -First, decompose $\hat{r}$ into components parallel and perpendicular to -$\hat{\beta}$ by calculating -$\mu = \hat{r}\cdot\hat{\beta}$ and +light, and $\hat{\beta}$ is a unit vector. Suppose also that the unit +vector $\hat{r}$ points toward a star in the barycenter frame of +reference (i.e. the ``actual'' position). and $\hat{r}'$ gives the +direction of the star in the observer's frame, (i.e. the apparent +position). + +First, decompose $\hat{r}$ into components parallel and perpendicular +to $\hat{\beta}$ by calculating $\mu = \hat{r}\cdot\hat{\beta}$ and $\vec{r}_\perp = \hat{r} - \mu \hat{\beta}$. @@ -1849,26 +1793,23 @@ relatively rapid rotation of the Earth from the motion of its rotational axis. -This section is largely a summary of -Chapter 5 of IERS Technical Note 32 \footnote{http://maia.usno.navy.mil/conv2003.html} -(hereafter IERS32), -which is a description of the implementation of the Resoltions of the -XXIVth General Assembly of the IAU, available from the same URL as above. -These two documents describe a set of conventions which have been in effect -since 2003-01-01. The conventions in effect before that date will not be -implemented by psLib. - +This section is largely a summary of Chapter 5 of IERS Technical Note +32 \footnote{http://maia.usno.navy.mil/conv2003.html} (hereafter +IERS32), which is a description of the implementation of the +Resoltions of the XXIVth General Assembly of the IAU, available from +the same URL as above. These two documents describe a set of +conventions which have been in effect since 2003-01-01. The +conventions in effect before that date will not be implemented by +psLib. \paragraph{Precession/Nutation} -The transform between the GCRS and the CIP/CEO coordinate systems is described -by the IAU 2000A precession-nutation model, which is accurate to the -0.2 mas level. -For higher accuracy the user must apply corrections to the model, which are tabulated by the IERS. - - - -The IAU 2000A precession-nutation model may be calculated in the following -way. First calculate the time $t$ as the number of Julian centuries since -2000-01-01T12:00:00 TT. +The transform between the GCRS and the CIP/CEO coordinate systems is +described by the IAU 2000A precession-nutation model, which is +accurate to the 0.2 mas level. For higher accuracy the user must +apply corrections to the model, which are tabulated by the IERS. + +The IAU 2000A precession-nutation model may be calculated in the +following way. First calculate the time $t$ as the number of Julian +centuries since 2000-01-01T12:00:00 TT. Next calculate the fundamental arguments of nutation using equations (40) @@ -1935,22 +1876,23 @@ \end{itemize} - -A FORTRAN reference implementation for the precession/nutation model is available from the IERS +A FORTRAN reference implementation for the precession/nutation model +is available from the IERS \footnote{http://maia.usno.navy.mil/conv2000/chapter5/XYS2000A.f}. The psLib results should agree with the reference implementation to within the limits of numerical precision. -Next, corrections to $X$, and $Y$ may be obtained from the IERS as part of -Bulletin A, or B. It is recommended to use the values published daily in -http://maia.usno.navy.mil/ser7/finals2000A.daily, which has the format -described by http://maia.usno.navy.mil/ser7/readme.finals2000A. The -quantities of interest are labeled dX and dY. Note that UT1$-$UTC and the -polar motion values are obtained from this same table. - -By convention, nutation terms with periods of less -than two days -are accounted for by the corresponding polar motion. So it is sufficient to -interpolate the corrections tabulated daily by the IERS, and take the result as -instantaneous values. +Next, corrections to $X$, and $Y$ may be obtained from the IERS as +part of Bulletin A, or B. It is recommended to use the values +published daily by USNO in the table +\code{finals2000A.daily}\footnote{http://maia.usno.navy.mil/ser7/finals2000A.daily}, +which has the format described by +\code{readme.finals2000A}\footnote{http://maia.usno.navy.mil/ser7/readme.finals2000A}. The +quantities of interest are labeled dX and dY. Note that UT1$-$UTC and +the polar motion values are obtained from this same table. + +By convention, nutation terms with periods of less than two days are +accounted for by the corresponding polar motion. So it is sufficient +to interpolate the corrections tabulated daily by the IERS, and take +the result as instantaneous values. The final step is to use $X$, $Y$, and $s$ to calculate the rotation @@ -1987,52 +1929,109 @@ \paragraph{Polar Motion} -The motion of the CIP in the ITRS is known as ``polar motion''. Similarly to -precession/nutation, the instantaneous position of the CIP in the -ITRS is specified by the quantites $x_p$, and $y_p$, and a third quantity, -$s'$, gives the position of the TEO with respect to the ITRS. -The values of $x_p$ and $y_p$ are published daily by the IERS in -http://maia.usno.navy.mil/ser7/finals2000A.daily, which has the format -described by http://maia.usno.navy.mil/ser7/readme.finals2000A. -The UT1$-$UTC, and the precession/nutation corrections (discussed elsewhere -in this document) come from this same source. - -The polar motion coordinates should be interpolated using a third order -polynomial, as described in -IERS Gazette \#13 \footnote{http://maia.usno.navy.mil/iers-gaz13}, -which gives a -FORTRAN reference implementation of the correct procedure. - -\tbd{reference to interpolation in this doc?} +The motion of the CIP in the ITRS is known as ``polar +motion''. Similarly to precession/nutation, the instantaneous position +of the CIP in the ITRS is specified by the quantites $x_p$, and $y_p$, +and a third quantity, $s'$, gives the position of the TEO with respect +to the ITRS. The values of $x_p$ and $y_p$ are published daily by the +IERS\footnote{http://maia.usno.navy.mil/ser7/finals2000A.daily}, with +a format described by their +\code{readme.finals2000A}\footnote{http://maia.usno.navy.mil/ser7/readme.finals2000A}. +The UT1$-$UTC, and the precession/nutation corrections (discussed +elsewhere in this document) come from this same source. + +The polar motion coordinates should be interpolated using a third +order polynomial, as described in IERS Gazette \#13 +\footnote{http://maia.usno.navy.mil/iers-gaz13}, which gives a FORTRAN +reference implementation of the correct procedure. The values published by the IERS are smoothed to remove noise and variations on the timescale of a day or less. There are two sources of -short timescale variations - tidal effects on the order of 0.1 milliarcseconds, -and short period nutation terms on the order of 15 microarcseconds. -Both of these effects may be modeled and added to the interpolated values -for higher accuracy. - -The tidal effects should be included using the FORTRAN reference implementation -of the Ray tidal model given in IERS Gazette \#13. This code should be -mimiced to machine accuracy by psLib. - -By definition of the CIP, nutation terms with periods less than 2 days are -not included in the IAU 2000A precession/nutation model. -So these motions -must be compensated for by their equivalent polar motions. These may -be calculated using a form similar to that of the precession/nutation $X$, -and $Y$. The constants to use are given in Table 5.1 of IERS32. -Note that only the terms with periods less than 2 days should be used. - -The quantity $s'$ may be approximated with microarcsecond accuracy over this -century by $s' = -4.7 \times 10^{-5} t$ in arcseconds. There is no need -to apply short timescale corrections to $s'$. +short timescale variations - tidal effects on the order of 0.1 +milliarcseconds, and short period nutation terms on the order of 15 +microarcseconds. Both of these effects may be modeled and added to +the interpolated values for higher accuracy. + +The tidal effects should be included by using the Ray tidal model +given in IERS Gazette \#13. The definition of this correction is +provided below. + +By definition of the CIP, nutation terms with periods less than 2 days +are not included in the IAU 2000A precession/nutation model. So these +motions must be compensated for by their equivalent polar +motions. These may be calculated using a form similar to that of the +precession/nutation $X$, and $Y$. The constants to use are given in +Table 5.1 of IERS32. Note that only the terms with periods less than +2 days should be used. + +The quantity $s'$ may be approximated with microarcsecond accuracy +over this century by $s' = -4.7 \times 10^{-5} t$ in arcseconds. There +is no need to apply short timescale corrections to $s'$. The transform from the ITRS to the CIP/TEO frame can be constructed by -first rotating about the X axis by $y_p$, then rotating about the X axis by -$x_p$, and finally rotating about the Z axis by $s'$. -The IERS reference implementation for this is given in the subroutine -POM2000 \footnote{http://maia.usno.navy.mil/conv2000/chapter5/POM2000.f}. -Note that we describe the transform toward celestial coordinates (upward in -figure X), in order to match the reference implementation. +first rotating about the X axis by $y_p$, then rotating about the X +axis by $x_p$, and finally rotating about the Z axis by $s'$. The +IERS reference implementation for this is given in the subroutine +POM2000 +\footnote{http://maia.usno.navy.mil/conv2000/chapter5/POM2000.f}. +Note that we describe the transform toward celestial coordinates +(upward in Figure~\ref{earthrot}), in order to match the reference +implementation. + +\subsubsection{Universal Time (UT1)} +\label{sec:ut1} + +Since 2003-01-01, UT1 has been defined as directly proportional to the +Earth Rotation Angle (see IERS Technical Note 32\footnote{IERS +Technical Note 32 - http://maia.usno.navy.mil/conv2003.html}). +Previous to that date, a different definition was in effect (see IERS +Technical Note 21\footnote{IERS Technical Note 21 - +http://maia.usno.navy.mil/conventions.html}). We will always use the +post-2003 definition. + +UT1 is continuously measured by the International Earth Rotation Service, and +tabulated values of the offset of UT1 from UTC are published at regular +intervals, along with predicted future values. IERS Bulletin A gives ``rapid +response'' values necessary for real-time and near real-time data analysis +(such as Pan-STARRS Otis and IPP subsystems). Bulletin B gives the results of a +final, definitive data reduction. An amalgam of Bulletin A and B values is +published daily on the IERS website\footnote{IERS Bulletin A \& B - +http://maia.usno.navy.mil/ser7/finals2000A.daily} along with a desciption of +the format\footnote{IERS finals2000A.daily table format - +http://maia.usno.navy.mil/ser7/readme.finals2000A}. + +The UT1 offsets should be interpolated using the prescription of IERS +Gazette \#13\footnote{IERS Gazette \#13 - +http://maia.usno.navy.mil/iers-gaz13}. This entails using a third +order polynomial to interpolate the table values, and then applying a +model for diurnal and semi-diurnal fluctuations due to tidal effects. +An example implimentation\footnote{interp.f - +ftp://maia.usno.navy.mil/dist/interp.f} written in Fortran is +available. The interpolated value of $dT$ must then have the tidal +correction from the Ray Tidal Model applied. + +\subsubsection{Ray Tidal Model} + +The Ray Model tidal corrections to X, Y, and dT are given by the the +Fortran code listed below. The input information is the epoch of +interest in MJD, while the output results are the corrections $C_x$, +$C_y$, and $C_dt$, which are in turn added to the interpolated values +determined above. + +\verbatiminput{raymodel.f} + +The Polar motion X, and Y coordinates are also important for determining the +orientation of the Earth with respect to the sky. This is also given in the +IERS publications references above, and should be interpolated in the same way. + +\subsubsection{Longitude} + +Longitudes are often expressed in the form of decimal degrees while the +algorithm for calculating GMST outputs seconds. + +\begin{equation} +1\degree = 240s +\end{equation} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsubsection{ITRS - Alt/Az} @@ -2344,5 +2343,5 @@ \subsection{Positions of Major Solar System Objects} -\tbd{ephemerides code to replace this} +\tbd{ephemerides code to replace this?} The SLALIB function \code{SLA_RDPLAN} returns the apparent position of @@ -2358,4 +2357,5 @@ \pagebreak + \section{Pan-STARRS Modules}