Changeset 3718

Timestamp:

Apr 19, 2005, 1:40:47 PM (21 years ago)

Author:

eugene

Message:

cleaned up ADD text

Location:

trunk/doc/pslib

Files:

• ChangeLogADD.tex (modified) (2 diffs)
• psLibADD.tex (modified) (26 diffs)
• rotations.ps (modified) ( previous)

trunk/doc/pslib/ChangeLogADD.tex

@@ -29 +29 @@
 \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 +48 @@
 \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}

trunk/doc/pslib/psLibADD.tex

@@ -1 @@
-%%% $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 +14 @@
 \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 +41 @@
 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 +348 @@
 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 +566 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\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 +760 @@
 
 \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 +1011 @@
 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 +1017 @@
 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 +1042 @@
 
 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 +1079 @@
 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 +1230 @@
 
 
-\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 +1347 @@
 
 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 +1402 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\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 +1435 @@
 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 +1448 @@
 \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 +1479 @@
     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 +1549 @@
 \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 +1605 @@
 \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 +1639 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\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 +1669 @@
 
 \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 +1732 @@
 \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 +1793 @@
 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 +1876 @@
 \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 +1929 @@
 \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 +2343 @@
 \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 +2357 @@
 
 \pagebreak 
+
 \section{Pan-STARRS Modules}