Index: trunk/doc/pslib/psLibADD.tex =================================================================== --- trunk/doc/pslib/psLibADD.tex (revision 3448) +++ trunk/doc/pslib/psLibADD.tex (revision 3564) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.69 2005-03-18 20:40:14 jhoblitt Exp $ +%%% $Id: psLibADD.tex,v 1.70 2005-03-30 21:14:48 eugene Exp $ \documentclass[panstarrs]{panstarrs} @@ -14,6 +14,8 @@ \project{Pan-STARRS Image Processing Pipeline} \organization{Institute for Astronomy} -\version{09} +\version{10} \docnumber{PSDC-430-006} + +\setcounter{tocdepth}{5} % lowest level to be included in toc \newcommand\citealt{} @@ -31,9 +33,12 @@ 01 & 2004 May 21 & Added section on 2D Chebyshev fitting, then removed. \\ \hline 02 & 2004 Jun 22 & modified stats specification \\ \hline -03--05 & ??? & ??? \\ \hline -06 & 2004 Sep 7 & Frozen for PSLib-2 \\ \hline +03 & 2004 Jul 13 & \\ \hline +04 & 2004 Aug 16 & \\ \hline +05 & 2004 Sep 01 & \\ \hline +06 & 2004 Sep 07 & Frozen for PSLib-2 \\ \hline 07 & 2004 Nov 24 & Frozen for Cycle 4 \\ \hline 08 & 2005 Jan 21 & Draft for Cycle 5 \\ \hline 09 & 2005 Feb 14 & Frozen for Cycle 5 \\ \hline +10 & 2005 Mar 21 & Draft for Cycle 6 \\ \hline \RevisionsEnd @@ -69,9 +74,9 @@ \pagenumbering{arabic} -\section{Pan-STARRS Library PSLib} - -\subsection{Math Utilities} - -\subsubsection{Sorting} +% \section{Pan-STARRS Library PSLib} + +\section{PSLib Math Utilities} + +\subsection{Sorting} A variety of sorting algorithms exist, with a wide range in speed for @@ -103,5 +108,5 @@ \code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}. -\subsubsection{Smoothing: Boxcar and Gaussian} +\subsection{Smoothing: Boxcar and Gaussian} \label{smooth} @@ -135,5 +140,5 @@ \end{equation} -\subsubsection{Statistics} +\subsection{Statistics} The general statistics function \code{psStats} performs a variety of @@ -148,10 +153,10 @@ sample and robust estimators. -\paragraph{Sample Statistics} +\subsubsection{Sample Statistics} We define the following statistical terms, assuming there is a set of data elements $x_i$ with (standard) errors $\sigma_i$. -\subparagraph{Mean} +\paragraph{Mean} The simple mean is defined as: @@ -160,5 +165,5 @@ \end{equation} -\subparagraph{Weighted Mean} +\paragraph{Weighted Mean} The weighted mean is defined as: @@ -170,5 +175,5 @@ standard definition of the mean. -\subparagraph{Median} +\paragraph{Median} The median is defined as the value for which 50\% of the data values @@ -181,5 +186,5 @@ calculating the sample median. -\subparagraph{Upper and Lower Quartiles} +\paragraph{Upper and Lower Quartiles} The upper and lower quartiles ($U_{\frac{1}{4}}$ and @@ -196,5 +201,5 @@ the sample quartiles. -\subparagraph{Standard Deviation} +\paragraph{Standard Deviation} The standard deviation of the sample is given by: @@ -218,5 +223,5 @@ -\paragraph{Clipped Statistics} +\subsubsection{Clipped Statistics} The clipped statistics are used to determine the mean and standard @@ -255,5 +260,5 @@ \bar{x}| > k \sigma_i$. -\paragraph{Robust Statistics} +\subsubsection{Robust Statistics} The robust version of the statistics provides estimators of basic @@ -312,5 +317,5 @@ quartiles are estimated in the same manner as above. -\paragraph{Histograms} +\subsubsection{Histograms} When calculating histograms in the presence of known errors in the @@ -342,5 +347,5 @@ Note that the total adds to one --- the number of values added. -\subsubsection{Matrix Operations} +\subsection{Matrix Operations} In this section, we define the linear algebra operations performed on @@ -361,5 +366,5 @@ \code{gsl_linalg_LU_decomp}. -\paragraph{LU Decomposition} +\subsubsection{LU Decomposition} \label{LUdecomp} @@ -402,5 +407,5 @@ \end{equation} -\paragraph{Calculate a matrix determinant} +\subsubsection{Calculate a matrix determinant} The determinant $D$ of a matrix $a_{ij}$ is calculated from the @@ -418,5 +423,5 @@ shall be used. -\paragraph{Solving a Linear Equation} +\subsubsection{Solving a Linear Equation} The LU decomposition of a matrix may be used to solve the @@ -438,5 +443,5 @@ \end{eqnarray} -\paragraph{Invert a matrix} +\subsubsection{Invert a matrix} Inversion of a matrix using the LU decomposition is performed by @@ -447,5 +452,5 @@ operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}. -\paragraph{Perform matrix addition, subtraction and multiplication} +\subsubsection{Perform matrix addition, subtraction and multiplication} Matrix binary arithmetic operations differ from image binary @@ -471,5 +476,5 @@ \times$. -\paragraph{Transpose a matrix} +\subsubsection{Transpose a matrix} The transpose of a matrix is simply the reorganization of the matrix @@ -484,5 +489,5 @@ where $M_{ij}$ is the matrix to be transposed. -\paragraph{Convert a matrix to a vector} +\subsubsection{Convert a matrix to a vector} Matrix-to-vector conversion is only defined for a matrix that has a @@ -493,7 +498,7 @@ matrix is converted to a \code{PS_DIMEN_TRANV}-type vector. -\subsubsection{Fitting} - -\paragraph{Chi-squared} +\subsection{Fitting} + +\subsubsection{Chi-squared} \label{chisq} @@ -506,5 +511,5 @@ \end{equation} -\paragraph{General Polynomial Fitting} +\subsubsection{General Polynomial Fitting} Given a set of data values $y_i$ with errors $\sigma_i$, related to @@ -531,7 +536,7 @@ (section~\ref{LUdecomp}). -\subsubsection{Non-linear Minimization} - -\paragraph{Levenberg-Marquardt Method} +\subsection{Non-linear Minimization} + +\subsubsection{Levenberg-Marquardt Method} In the Levenberg-Marquardt Method (LMM; see NR \S 15.5), we make a @@ -612,5 +617,5 @@ -\paragraph{Powell's method} +\subsubsection{Powell's method} Powell's method is a type of ``Direction Set'' methods in @@ -621,5 +626,5 @@ manner until the advances along the vectors are smaller than some pre-defined tolerance. Such direction set methods, including Powell's -Quadratically Convergent method are discussed in NR\S10.5. +Quadratically Convergent method are discussed in NR \S 10.5. We will use for our algorithm the modified version of Powell's @@ -665,5 +670,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{Polynomials} +\subsection{Polynomials} \label{sec:polynomials} @@ -703,5 +708,5 @@ $-1 < x < 1$. -\paragraph{Multi-dimensional polynomials} +\subsubsection{Multi-dimensional polynomials} Multi-dimensional polynomials shall be composed of multiplications of @@ -711,5 +716,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{(Fast) Fourier Transforms} +\subsection{(Fast) Fourier Transforms} (Fast) Fourier Transforms (FFTs) shall be implemented using the @@ -717,5 +722,5 @@ library}. -\paragraph{FFTW Plans} +\subsubsection{FFTW Plans} FFTW requires the user to create a ``plan'' for each transform size, @@ -733,5 +738,5 @@ initialisation of the PSLib FFT functions and saved at the conclusion. -\paragraph{Function mapping} +\subsubsection{Function mapping} The forward and reverse transforms call the corresponding @@ -750,5 +755,5 @@ place to avoid the need to pad the input array to hold the output. -\paragraph{More Complicated Functions} +\subsubsection{More Complicated Functions} \code{psFFTCrossCorrelate()} and \code{psFFTConvolve()} both involve @@ -772,7 +777,199 @@ caller, and choose to normalise by $1/N^2$. +\subsection{Image Manipulations} + +\subsubsection{Interpolation} + +Interpolation is needed in various image manipulation operations, +including rotation and resampling. We have specified a function to +perform the interpolation using one of several possible interpolation +methods, defined below. It is important in the discussions that +follow to remember that a pixel with column,row if $i,j$ has +coordinate at the center of $i+0.5,j+0.5$ and corners with coordinates +from $i,j$ to $i+1,j+1$. Thus, the interpolation of a coordinate +$x,y$ = 5.0,4.0 is a value midway between the four pixels with +column,row of (5,4), (5,5), (6,4), (6,5). + +\paragraph{Nearest Pixel Interpolation ({\tt PS\_INTERPOLATE\_FLAT})} + +In this interpolation, the value of the closest pixel is returned. +This is equivalent to pixel duplication or replication. + +\paragraph{Bilinear Interpolation ({\tt PS\_INTERPOLATE\_BILINEAR})} + +In this interpolation, the value at the coordinate is calculated using +linear interpolation in two dimensions from the four nearest neighbor +pixels. The bilinear interpolation value at a coordinate $x,y$ +depends on the four nearest neighbor pixels and the fractional +distance $fx,fy$ of the given coordinates from the centers of those +four pixels. Consider four neighboring pixels at column,row of $i,j$, +$i+1,j$, $i,j+1$, and $i+1,j+1$ with pixel values $V_{0,0}$, +$V_{1,0}$, $V_{0,1}$, $V_{1,1}$. The value at $x,y$ is given by: +\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \] +This expression is more efficiently evaluated by factoring and +calculating the expresion as: +\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \] +\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \] + +Note that the values of $f_x$ and $f_y$ require some care. Given a +coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 - +int(f_x - 0.5)$. For example, when interpolating the value at +(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5) +and the fractional coordinate values $f_x, f_y = 0.3, 0.7$. The +resulting coordinate would be contained within the pixel at column,row +(5,5). + +\paragraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})} + +Because it would be slow to specify the size of the kernel +dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8 +pixels in each dimension (a kernel of size 2 pixels in each dimension +is handled by the bilinear interpolation). These correspond to the +options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and +\code{PS_INTERPOLATE_LANCZOS4}, respectively. + +Given a position on the input image, $(x_0,y_0)$, a kernel is derived +according to pixels local to the position: +\begin{equation} + h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N) +\end{equation} +where +\begin{eqnarray} + \delta x & = & x - x_0 \\ + \delta y & = & y - y_0 \\ + {\rm sinc}(z) & = & \sin(z)/z +\end{eqnarray} +and $N$ corresponds to the choice of kernel size. For $N = 2$, the +kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le +2$). For $N = 3$, the kernel size is 6 pixels in each dimension +(i.e., $-3 < \delta x \le 3$). For $N = 4$, the kernel size is 8 +pixels in each dimension (i.e., $-4 < \delta x \le 4$). + +The interpolated value at the given position, $(x_0,y_0)$, is then +simply the dot product of the kernel and the fluxes: +\begin{equation} + f(x_0,y_0) = \sum_R f(x,y) h(x,y) +\end{equation} +where $R$ is the region defined by the kernel size, and $f(x,y)$ is +the flux at the pixel position. + +For further information, see the +\href{http://terapix.iap.fr/IMG/pdf/swarp.pdf}{SWarp manual}. + +\subsubsection{Image Cuts and Slices} + +Several functions specify operations which manipulate a collection of +pixels to return a statistic on the pixel collection. In the simplest +case, these are trivial to define: if the boundaries of the region of +interest are specified along integral pixel coordinates, then the +pixels used to measure the statistic are always an exact integer. +This is the case for the function \code{psImageSlice} which requires a +starting coordinate which is an integer and a width in both dimensions +which is an integer. For the case of the functions \code{psImageCut} +and \code{psImageRadialCut}, the situation is a bit more subtle. In +both of these cases, the region is unlikely to contain only whole +pixels and some choices must be made. + +One posibility which we reject is to identify the fractional pixels +which are overlapped by the region of interest and add that fraction +of the pixel's flux when calculating the statistic of interest. This +is computationally intensive, and not necessarily well defined for all +statistics. + +In PSLib, we instead identify the pixels overlapped by the region, use +the complete set of pixel values, treating all pixels equally, and +renormalize as needed. To perform this, the region of interest is +laid on top of the image pixels. Any pixels which overlap the region +are identified as part of the input sample. The statistic (ie, sample +mean, robust mode, etc), is then calculated on this collection of +pixels. If the output statistic is an average value, the measured +value is reported. If the output statistic is a sum value (sum of +counts, sum of pixels), then the value is renormalized by the ratio of +pixels used in the calculation to the pixel area of the region of +interest. For example, if the sum within a radial aperture is +requested, the circle of the specified radius and center is placed on +the pixel grid. Any pixels which touch the circle are then placed in +a list to be analysed. The statistic of interest is the measured for +this collection of pixels. In the case of a circular aperture which +is centered at the coordinate (2,2) and has a radius of 2, the number +of pixels which are touched by the circle is 16, while the total pixel +area of the circle is 12.57 square pixels. In this case, the pixel +sum is renormalized by the ratio (12.57/16.00). + +\paragraph{Radial Cuts} + +Consider an image with pixels $x_i,y_i$ and a reference coordinate +$x_c, y_c$. We want to construct a radial cut by measuring statistics +for pixels in a sequence of radial annulii $r_s < r < r_e$. For each +annulus, we need to select the pixels which fall within this annulus. +The coordinates of the center of pixel $i,j$ are $i+0.5,j+0.5$. A +given pixel has a distance from the reference coordinate of $dX = x_c +- i - 0.5, dY = y_c - j - 0.5$. The pixels to be used for a given +radial annulus are all of those pixels for which $r_s < \sqrt{dX^2 + + dY^2} < r_e$. This is more efficiently calculated by comparing the +square of the radii and distances. All pixels which satisfy the above +condition are included in a specific annular radius. All average +quantities are calculated directly from the pixel ensemble +statistics. + +\paragraph{Arbitrary Linear Cuts} + +Select the pixels which lie along a line following steps of 1 pixel +length: + +\begin{verbatim} + + dX = xe - xs; + dY = ye - ys; + L = hypot (dX, dY); + dX = dX / L; + dY = dY / L; + + REALLOCATE (xvec[0].elements, float, MAX (L, 1)); + REALLOCATE (yvec[0].elements, float, MAX (L, 1)); + xvec[0].Nelements = L; + yvec[0].Nelements = L; + + V = (float *)buf[0].matrix.buffer; + for (i = 0; i < L; i++) { + xi = xs + i*dX - 0.5; + yi = ys + i*dY - 0.5; + xvec[0].elements[i] = i; + yvec[0].elements[i] = V[xi + Nx*yi]; + } +\end{verbatim} + +\subsubsection{Image Rotation} + +Image rotation can be performed in two possible ways under different +circumstances, identified in the following discussion. + +In the simplest case, the rotation angle is an integer multiple of 90 +degrees ($\pi/2$ rad). In these cases, the input and output pixels +have a one-to-one mapping. If the input image has dimensions of $N_x, +N_y$, then the output image will have dimensions of either $N_x, N_y$ +(for even multiples of 90 degrees) or $N_y, N_x$ (for odd multiples). + +If the angle of the rotation is not a multiple of 90, then the output +pixels necessarily result from the interpolation of several input +pixels. In this case, for an input image of dimensions $N_x, N_y$ and +rotation angle $\theta$, the output image has dimensions $Lx = |N_x +\cos \theta| + |N_y \sin \theta|$ and $Ly = |N_x \sin \theta| + |N_y +\cos \theta|$, each dimension rounded up to the nearest integer as +needed. Every pixel in the output image is in general derived from an +interpolation over 4 neighboring pixels. The coordinate of a pixel in +the output image ($i,j$) corresponds to a fractional pixel coordinate +($x,y$) in the input image according to: +\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \] +\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \] +where the offset coordinate ($i_o,j_o$) depends on the sign of the +sine of the angle $\theta$. If the sign of that sine is positive, the +offset coordinate is ($N_y\sin\theta$,0), otherwise it is +(0,$-N_x\sin\theta$). + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Astronomy Utilities} +\pagebreak +\section{PSLib Astronomy Utilities} Most of the astronomy utilities will be implemented through wrapping @@ -782,5 +979,5 @@ the next release} -\subsubsection{Time} +\subsection{Time} Correct time representation is \emph{critical} in astronomical software. PSLib @@ -799,5 +996,5 @@ ``1970-01-01T00:00:00Z'' UTC. -\paragraph{Coordinated Universal Time (UTC)} +\subsubsection{Coordinated Universal Time (UTC)} Coordinated Univeral Time (UTC) is defined by the International @@ -829,5 +1026,5 @@ timezone is forbidden.} -\paragraph{International Atomic Time (TAI)} +\subsubsection{International Atomic Time (TAI)} International Atomic Time or Temps Atomique International (TAI) is a system of @@ -845,5 +1042,5 @@ seconds since the UNIX epoch of ``1970-01-01T00:00:00Z''. -\paragraph{Leap-seconds} +\subsubsection{Leap-seconds} Leap seconds keep UTC within 0.9s of UT1. The offset between TAI and @@ -883,5 +1080,5 @@ This data is available from: \code{ftp://maia.usno.navy.mil/ser7/tai-utc.dat} -\paragraph{Gregorian dates to seconds} +\subsubsection{Gregorian dates to seconds} The Perl code below, based on an algorithm described in the book ``Calendrical @@ -1003,4 +1200,5 @@ \end{verbatim} Outputs year, month, day as \code{$y, $m, $d}. +%$ \emph{The above code was taken [and slightly altered] from @@ -1011,5 +1209,5 @@ -\paragraph{Universal Time (UT1)} +\subsubsection{Universal Time (UT1)} \label{sec:ut1} @@ -1054,5 +1252,5 @@ IERS publications references above, and should be interpolated in the same way. -\paragraph{Julian Date and Modified Julian Date} +\subsubsection{Julian Date and Modified Julian Date} The follow definitions of Julian Date (JD) and Modified Julian Date (MJD) was @@ -1062,5 +1260,5 @@ http://www.iers.org/iers/earth/resolutions/UAI\_b1.html}. -\subparagraph{Julian Date} +\paragraph{Julian Date} \begin{verbatim} @@ -1089,5 +1287,5 @@ \end{verbatim} -\subparagraph{Modified Julian Date} +\paragraph{Modified Julian Date} \begin{verbatim} @@ -1097,5 +1295,5 @@ \end{verbatim} -\subparagraph{JD and MJD conversion} +\paragraph{JD and MJD conversion} Conversion between \code{psTime} values and MJD and JD are determined @@ -1115,5 +1313,5 @@ \end{equation} -\paragraph{Terrestrial Time (TT)} +\subsubsection{Terrestrial Time (TT)} Terrestrial Time (TT) is defined as a fixed offset from TAI. @@ -1123,5 +1321,5 @@ \end{equation} -\paragraph{TT as Julian Centuries since J2000.0} +\subsubsection{TT as Julian Centuries since J2000.0} The algorithm for calulating GMST requires TT formated in Julian centruies @@ -1131,5 +1329,5 @@ \end{equation} -\paragraph{UT1 as Julian Centuries since J2000.0} +\subsubsection{UT1 as Julian Centuries since J2000.0} The algorithm for calulating GMST requires UT1 be formated in Julian centuries @@ -1140,5 +1338,18 @@ \end{equation} -\paragraph{Greenwich Mean Sidereal Time (GMST)} +\subsubsection{Local Mean Sidereal Time (LMST)} + +Local Mean Sidereal Time (LMST) is Greenwich Mean Sideral Time (GMST) plus the +observer's location in East longitude. Calculating LMST requires the input of +Universal Time (UT1), Terrestrial Dynamical Time (TT) and a longitude (measured +East of Greenwich). + +\begin{equation} +LMST = GMST00(t_u, t) + longitude +\end{equation} + +Gives $LMST$ in seconds. + +\subsubsection{Greenwich Mean Sidereal Time (GMST)} Greenwich Mean Sidereal Time (GMST) is caclulated from UT1 and TT. This @@ -1159,6 +1370,5 @@ Gives $GMST00$ in seconds. - -\paragraph{Longitude} +\subsubsection{Longitude} Longitudes are often expressed in the form of decimal degrees while the @@ -1169,18 +1379,6 @@ \end{equation} -\paragraph{Local Mean Sidereal Time (LMST)} - -Local Mean Sidereal Time (LMST) is Greenwich Mean Sideral Time (GMST) plus the -observer's location in East longitude. Calculating LMST requires the input of -Universal Time (UT1), Terrestrial Dynamical Time (TT) and a longitude (measured -East of Greenwich). - -\begin{equation} -LMST = GMST00(t_u, t) + longitude -\end{equation} - -Gives $LMST$ in seconds. - -\paragraph{Polar Motion} +\subsubsection{Polar Motion} +\tbd{move this to Earth Motion section} The polar coordinates, $x_p$ and $y_p$, required for the transformation from @@ -1190,196 +1388,178 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{Astronomical Image Manipulations} - -\paragraph{Interpolation} - -Interpolation is needed in various image manipulation operations, -including rotation and resampling. We have specified a function to -perform the interpolation using one of several possible interpolation -methods, defined below. It is important in the discussions that -follow to remember that a pixel with column,row if $i,j$ has -coordinate at the center of $i+0.5,j+0.5$ and corners with coordinates -from $i,j$ to $i+1,j+1$. Thus, the interpolation of a coordinate -$x,y$ = 5.0,4.0 is a value midway between the four pixels with -column,row of (5,4), (5,5), (6,4), (6,5). - -\subparagraph{Nearest Pixel Interpolation ({\tt PS\_INTERPOLATE\_FLAT})} - -In this interpolation, the value of the closest pixel is returned. -This is equivalent to pixel duplication or replication. - -\subparagraph{Bilinear Interpolation ({\tt PS\_INTERPOLATE\_BILINEAR})} - -In this interpolation, the value at the coordinate is calculated using -linear interpolation in two dimensions from the four nearest neighbor -pixels. The bilinear interpolation value at a coordinate $x,y$ -depends on the four nearest neighbor pixels and the fractional -distance $fx,fy$ of the given coordinates from the centers of those -four pixels. Consider four neighboring pixels at column,row of $i,j$, -$i+1,j$, $i,j+1$, and $i+1,j+1$ with pixel values $V_{0,0}$, -$V_{1,0}$, $V_{0,1}$, $V_{1,1}$. The value at $x,y$ is given by: -\[ V = (V_{0,0}(1 - f_x) + V_{1,0}f_x)(1 - f_y) + (V_{0,1}(1-f_x) + V_{1,1}f_x)f_y \] -This expression is more efficiently evaluated by factoring and -calculating the expresion as: -\[ r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \] -\[ V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \] - -Note that the values of $f_x$ and $f_y$ require some care. Given a -coordinate $x,y$, the value of $f_x$ is calculated as $f_x - 0.5 - -int(f_x - 0.5)$. For example, when interpolating the value at -(5.8.5.2), the relevant neighbor pixels are (5,4), (6,4), (5,5), (6,5) -and the fractional coordinate values $f_x, f_y = 0.3, 0.7$. The -resulting coordinate would be contained within the pixel at column,row -(5,5). - -\subparagraph{Sinc Interpolation ({\tt PS\_INTERPOLATE\_LANCZOS[234]})} - -Because it would be slow to specify the size of the kernel -dynamically, we specify three hard-coded kernel sizes: 4, 6 and 8 -pixels in each dimension (a kernel of size 2 pixels in each dimension -is handled by the bilinear interpolation). These correspond to the -options \code{PS_INTERPOLATE_LANCZOS2}, \code{PS_INTERPOLATE_LANCZOS3} and -\code{PS_INTERPOLATE_LANCZOS4}, respectively. - -Given a position on the input image, $(x_0,y_0)$, a kernel is derived -according to pixels local to the position: -\begin{equation} - h(x,y) = {\rm sinc}(\pi \delta x) {\rm sinc}(\pi \delta x / N) \rm{sinc}(\pi \delta y) \rm{sinc}(\pi \delta y / N) -\end{equation} -where -\begin{eqnarray} - \delta x & = & x - x_0 \\ - \delta y & = & y - y_0 \\ - {\rm sinc}(z) & = & \sin(z)/z -\end{eqnarray} -and $N$ corresponds to the choice of kernel size. For $N = 2$, the -kernel size is 4 pixels in each dimension (i.e., $-2 < \delta x \le -2$). For $N = 3$, the kernel size is 6 pixels in each dimension -(i.e., $-3 < \delta x \le 3$). For $N = 4$, the kernel size is 8 -pixels in each dimension (i.e., $-4 < \delta x \le 4$). - -The interpolated value at the given position, $(x_0,y_0)$, is then -simply the dot product of the kernel and the fluxes: -\begin{equation} - f(x_0,y_0) = \sum_R f(x,y) h(x,y) -\end{equation} -where $R$ is the region defined by the kernel size, and $f(x,y)$ is -the flux at the pixel position. - -For further information, see the -\href{http://terapix.iap.fr/IMG/pdf/swarp.pdf}{SWarp manual}. - -\paragraph{Image Cuts and Slices} - -Several functions specify operations which manipulate a collection of -pixels to return a statistic on the pixel collection. In the simplest -case, these are trivial to define: if the boundaries of the region of -interest are specified along integral pixel coordinates, then the -pixels used to measure the statistic are always an exact integer. -This is the case for the function \code{psImageSlice} which requires a -starting coordinate which is an integer and a width in both dimensions -which is an integer. For the case of the functions \code{psImageCut} -and \code{psImageRadialCut}, the situation is a bit more subtle. In -both of these cases, the region is unlikely to contain only whole -pixels and some choices must be made. - -One posibility which we reject is to identify the fractional pixels -which are overlapped by the region of interest and add that fraction -of the pixel's flux when calculating the statistic of interest. This -is computationally intensive, and not necessarily well defined for all -statistics. - -In PSLib, we instead identify the pixels overlapped by the region, use -the complete set of pixel values, treating all pixels equally, and -renormalize as needed. To perform this, the region of interest is -laid on top of the image pixels. Any pixels which overlap the region -are identified as part of the input sample. The statistic (ie, sample -mean, robust mode, etc), is then calculated on this collection of -pixels. If the output statistic is an average value, the measured -value is reported. If the output statistic is a sum value (sum of -counts, sum of pixels), then the value is renormalized by the ratio of -pixels used in the calculation to the pixel area of the region of -interest. For example, if the sum within a radial aperture is -requested, the circle of the specified radius and center is placed on -the pixel grid. Any pixels which touch the circle are then placed in -a list to be analysed. The statistic of interest is the measured for -this collection of pixels. In the case of a circular aperture which -is centered at the coordinate (2,2) and has a radius of 2, the number -of pixels which are touched by the circle is 16, while the total pixel -area of the circle is 12.57 square pixels. In this case, the pixel -sum is renormalized by the ratio (12.57/16.00). - -\subparagraph{Radial Cuts} - -Consider an image with pixels $x_i,y_i$ and a reference coordinate -$x_c, y_c$. We want to construct a radial cut by measuring statistics -for pixels in a sequence of radial annulii $r_s < r < r_e$. For each -annulus, we need to select the pixels which fall within this annulus. -The coordinates of the center of pixel $i,j$ are $i+0.5,j+0.5$. A -given pixel has a distance from the reference coordinate of $dX = x_c -- i - 0.5, dY = y_c - j - 0.5$. The pixels to be used for a given -radial annulus are all of those pixels for which $r_s < \sqrt{dX^2 + - dY^2} < r_e$. This is more efficiently calculated by comparing the -square of the radii and distances. All pixels which satisfy the above -condition are included in a specific annular radius. All average -quantities are calculated directly from the pixel ensemble -statistics. - -\subparagraph{Arbitrary Linear Cuts} - -Select the pixels which lie along a line following steps of 1 pixel -length: +\subsection{2D transformations} + +In PSLib, we implement 2-dimensional transformations using +\code{psPlaneTransform}, which contains a matrix of polynomial +coefficients for each dimension. Since we are using these to model +the real world, where, for example, a particular point on the detector +maps to a particular point on the sky, we consider only +transformations that are ``one-to-one''. This makes it possible to +speak of inverse transformations, and of combining multiple +transformations. + +Given a transformation, $f(x,y)$, the inverse transformation, +$g(x,y)$, is that for which $g(f(x,y)) = (x,y)$ for $(x,y)$ over the +range of interest (not necessarily the entire set of real numbers). + +Given two transformations, $f(x,y)$ and $g(x,y)$, the combined +transformation is the transformation, $h(x,y) = g(f(x,y))$ for $(x,y)$ +over the range of interest (not necessarily the entire set of real +numbers). + +Both of these operations are straightforward if the transformation is +linear. If the function $(u,v) = f(x,y)$ is: +\begin{eqnarray} +u & = & a + bx + cy \\ +v & = & d + ex + fy +\end{eqnarray} +then the inverse transformation $(x,y) = g(u,v)$ is: +\begin{eqnarray} +x & = & (-fa+cd)/\Delta + fu/\Delta - cv/\Delta \\ +y & = & (ae-bd)/\Delta - eu/\Delta + bv/\Delta +\end{eqnarray} +where $\Delta = bf - ce$ is the matrix determinant. Given two +functions $f_i(x,y)$ for $i=1,2$: +\begin{eqnarray} +u & = & a_i + b_i x + c_i y \\ +v & = & d_i + e_i x + f_i y +\end{eqnarray} +then the combined transformation, $(u,v) = f_2(f_1(x,y))$ is: +\begin{eqnarray} +u & = & (a_2 + b_2 a_1 + c_2 d_1) + (b_2 b_1 + c_2 e_1) x + (b_2 c_1 + c_2 f_1) y \\ +v & = & (d_2 + e_2 a_1 + f_2 d_1) + (e_2 b_1 + f_2 e_1) x + (e_2 c_1 + f_2 f_1) y +\end{eqnarray} + +When the transformations are not linear, the inverse and combined +transformations can be estimated by sampling a grid over the region of +interest, calculating the transformation (or double transformation) +for each sample, and using this information to derive the best fit +transformation that produces the inverse or combined transformation. +The inverse transformation should be of the same order as that of the +forward transformation, while the combined transformation should be of +the higher order of the two component transformations. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Spherical Rotations with Quaternions} + +\subsubsection{Quaternion Construction} + +The following describes the algorithms needed to implement 3-D +rotations in terms of quaternions. A quaternion is an ordered set of +four numbers, $\bar{q} = (q_0, q_1, q_2, q_3)$. A rotation of angle +$\theta$ about the axis defined by the unit vector $(v_x, v_y, v_z)$ +has quaternion components: +\begin{eqnarray} +q_0 & = & v_x sin(\theta/2), \\ +q_1 & = & v_y sin(\theta/2), \\ +q_2 & = & v_z sin(\theta/2), and \\ +q_3 & = & cos(\theta/2). \\ +\end{eqnarray} +Note that the sine and cosine are taken of the half angle of the +rotation. Note also that this implies that the quaternion components +are normalized such that $|\bar{q}| \def q_0^2 + q_1^2 + q_2^2 + q_3^2 += 1$. + +The 3-vector representation of the angle of the pole is determined +from the coordinate of the pole ($\alpha_p, \delta_p$) by: +\begin{eqnarray} +v_x & = & \cos \delta_p \cos \alpha_p \\ +v_y & = & \cos \delta_p \sin \alpha_p \\ +v_x & = & \sin \delta_p \\ +\end{eqnarray} + +\subsubsection{Combining Two Rotations} + +Given two quaternions $\bar{p1}$ and $\bar{p2}$, 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 +by: \begin{verbatim} - - dX = xe - xs; - dY = ye - ys; - L = hypot (dX, dY); - dX = dX / L; - dY = dY / L; - - REALLOCATE (xvec[0].elements, float, MAX (L, 1)); - REALLOCATE (yvec[0].elements, float, MAX (L, 1)); - xvec[0].Nelements = L; - yvec[0].Nelements = L; - - V = (float *)buf[0].matrix.buffer; - for (i = 0; i < L; i++) { - xi = xs + i*dX - 0.5; - yi = ys + i*dY - 0.5; - xvec[0].elements[i] = i; - yvec[0].elements[i] = V[xi + Nx*yi]; - } +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} -\paragraph{Image Rotation} - -Image rotation can be performed in two possible ways under different -circumstances, identified in the following discussion. - -In the simplest case, the rotation angle is an integer multiple of 90 -degrees ($\pi/2$ rad). In these cases, the input and output pixels -have a one-to-one mapping. If the input image has dimensions of $N_x, -N_y$, then the output image will have dimensions of either $N_x, N_y$ -(for even multiples of 90 degrees) or $N_y, N_x$ (for odd multiples). - -If the angle of the rotation is not a multiple of 90, then the output -pixels necessarily result from the interpolation of several input -pixels. In this case, for an input image of dimensions $N_x, N_y$ and -rotation angle $\theta$, the output image has dimensions $Lx = |N_x -\cos \theta| + |N_y \sin \theta|$ and $Ly = |N_x \sin \theta| + |N_y -\cos \theta|$, each dimension rounded up to the nearest integer as -needed. Every pixel in the output image is in general derived from an -interpolation over 4 neighboring pixels. The coordinate of a pixel in -the output image ($i,j$) corresponds to a fractional pixel coordinate -($x,y$) in the input image according to: -\[ x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \] -\[ y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \] -where the offset coordinate ($i_o,j_o$) depends on the sign of the -sine of the angle $\theta$. If the sign of that sine is positive, the -offset coordinate is ($N_y\sin\theta$,0), otherwise it is -(0,$-N_x\sin\theta$). - -\subsubsection{Celestial Coordinate Conversions} +\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 +unit vector, and whose fourth component is zero. To rotate this vector +by a quaternion $\bar{p1}$, you apply the formula above for combining +two quaternions. The rotated vector is found in the first three +components of the resulting quaternion, $\bar{p}$. + +\subsubsection{Rotation Matrix} + +The rotation matrix representation of a rotation may be derived +directly from the quaternion representation. The following formulae +convert a quaternion to a rotation matrix: + +\begin{eqnarray} + rot_{x,x} & = & q_0 q_0 - q_1 q_1 - q_2 q_2 + q_3 q_3 \\ + rot_{y,y} & = & -q_0 q_0 + q_1 q_1 - q_2 q_2 + q_3 q_3 \\ + rot_{z,z} & = & -q_0 q_0 - q_1 q_1 + q_2 q_2 + q_3 q_3 \\ + + rot_{x,y} & = & 2 (q_0 q_1 + q_2 q_3) \\ + rot_{y,x} & = & 2 (q_0 q_1 - q_2 q_3) \\ + + rot_{x,z} & = & 2 (q_0 q_2 - q_1 q_3) \\ + rot_{z,x} & = & 2 (q_0 q_2 + q_1 q_3) \\ + + rot_{y,z} & = & 2 (q_1 q_2 + q_0 q_3) \\ + rot_{z,y} & = & 2 (q_1 q_2 - q_0 q_3) \\ +\end{eqnarray} + +\subsubsection{Conversion to Other Representations} + +You may convert a rotation matrix, m, to a quaternion, p, with the following +code: + +\begin{verbatim} +double diag_sum[3]; +int maxi; +double recip; + +diag_sum[0]=1+m[0][0]-m[1][1]-m[2][2]; +diag_sum[1]=1-m[0][0]+m[1][1]-m[2][2]; +diag_sum[2]=1-m[0][0]-m[1][1]+m[2][2]; +diag_sum[3]=1+m[0][0]+m[1][1]+m[2][2]; + + +maxi=0; +for(i=1;i<4;++i) { + if(diag_sum[i]>diag_sum[maxi]) maxi=i; +} + + +p[maxi]=0.5*sqrt(diag_sum[maxi]); +recip=1./(4.*p[maxi]); + +if(maxi==0) { + p[1]=recip*(m[0][1]+m[1][0]); + p[2]=recip*(m[2][0]+m[0][2]); + p[3]=recip*(m[1][2]-m[2][1]); + +} else if(maxi==1) { + p[0]=recip*(m[0][1]+m[1][0]); + p[2]=recip*(m[1][2]+m[2][1]); + p[3]=recip*(m[2][0]-m[0][2]); + +} else if(maxi==2) { + p[0]=recip*(m[2][0]+m[0][2]); + p[1]=recip*(m[1][2]+m[2][1]); + p[3]=recip*(m[0][1]-m[1][0]); + +} else if(maxi==3) { + p[0]=recip*(m[1][2]-m[2][1]); + p[1]=recip*(m[2][0]-m[0][2]); + p[2]=recip*(m[0][1]-m[1][0]); +} +\end{verbatim} + +\subsection{Celestial Coordinate Conversions} Changes between spherical coordinate systems (ie, Ecliptic, Galactic, @@ -1426,5 +1606,5 @@ the forward transformation. -\paragraph{Galactic to ICRS} +\subsubsection{Galactic to ICRS} The appropriate values, from the Hipparcos and Tycho Catalogues are: @@ -1435,5 +1615,5 @@ \end{eqnarray} -\paragraph{Ecliptic to ICRS} +\subsubsection{Ecliptic to ICRS} The appropriate values, from Zombeck, are: @@ -1445,5 +1625,5 @@ where $T$ is the time in Julian centuries since 1900. -\paragraph{Precession} +\subsubsection{Precession} The appropriate values, from Elixir, are: @@ -1457,5 +1637,5 @@ -\paragraph{Suggested test cases} +\subsubsection{Suggested test cases} $(\alpha,\delta) = (0^\circ,0^\circ)$ transforms to Galactic @@ -1479,289 +1659,6 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{2D transformations} - -In PSLib, we implement 2-dimensional transformations using -\code{psPlaneTransform}, which contains a matrix of polynomial -coefficients for each dimension. Since we are using these to model -the real world, where, for example, a particular point on the detector -maps to a particular point on the sky, we consider only -transformations that are ``one-to-one''. This makes it possible to -speak of inverse transformations, and of combining multiple -transformations. - -Given a transformation, $f(x,y)$, the inverse transformation, -$g(x,y)$, is that for which $g(f(x,y)) = (x,y)$ for $(x,y)$ over the -range of interest (not necessarily the entire set of real numbers). - -Given two transformations, $f(x,y)$ and $g(x,y)$, the combined -transformation is the transformation, $h(x,y) = g(f(x,y))$ for $(x,y)$ -over the range of interest (not necessarily the entire set of real -numbers). - -Both of these operations are straightforward if the transformation is -linear. If the function $(u,v) = f(x,y)$ is: -\begin{eqnarray} -u & = & a + bx + cy \\ -v & = & d + ex + fy -\end{eqnarray} -then the inverse transformation $(x,y) = g(u,v)$ is: -\begin{eqnarray} -x & = & (-fa+cd)/\Delta + fu/\Delta - cv/\Delta \\ -y & = & (ae-bd)/\Delta - eu/\Delta + bv/\Delta -\end{eqnarray} -where $\Delta = bf - ce$ is the matrix determinant. Given two -functions $f_i(x,y)$ for $i=1,2$: -\begin{eqnarray} -u & = & a_i + b_i x + c_i y \\ -v & = & d_i + e_i x + f_i y -\end{eqnarray} -then the combined transformation, $(u,v) = f_2(f_1(x,y))$ is: -\begin{eqnarray} -u & = & (a_2 + b_2 a_1 + c_2 d_1) + (b_2 b_1 + c_2 e_1) x + (b_2 c_1 + c_2 f_1) y \\ -v & = & (d_2 + e_2 a_1 + f_2 d_1) + (e_2 b_1 + f_2 e_1) x + (e_2 c_1 + f_2 f_1) y -\end{eqnarray} - -When the transformations are not linear, the inverse and combined -transformations can be estimated by sampling a grid over the region of -interest, calculating the transformation (or double transformation) -for each sample, and using this information to derive the best fit -transformation that produces the inverse or combined transformation. -The inverse transformation should be of the same order as that of the -forward transformation, while the combined transformation should be of -the higher order of the two component transformations. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsubsection{Projections} - -We implement three types of projections: {\em zenithal}, {\em -cylindrical} and {\em pseudocylindrical}, each requiring slightly -different handling. Our representations are based on the treatment of -projections presented by -\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \& -Calabretta (1995, ADASS, 4, 233)}. In all of these projections, we -are converting from a spherical coordinate $\alpha,\delta$ to a linear -(2-D) coordinate $x_p,y_p$. The projection is defined by the -projection type, the projection center ($\alpha_p, \delta_p$) and the -the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$). - -In the structure, \code{psProjection}, the projection type is defined -by the element \code{type}, the projection center $\alpha_p,\delta_p$ -is defined by the elements \code{R,D}, and the plate scales, -$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}. The plate -scales are applied independently to the $x$ and $y$ coordinates to -convert them to the corresponding linear units (ie, pixels): -% -\begin{eqnarray} -x_p & = & \rho_x x \\ -y_p & = & \rho_y y \\ -\end{eqnarray} -% -In the discussions below, we ignore this last step (or first step, -depending on the direction of the conversion). - -\paragraph{Zenithal Projections} - -The {\em zenithal} projections are defined relative to a set of -spherical coordinates with pole at the center of the projection -($\alpha_p, \delta_p$), and which thus represents a coordinate system -rotated relative to the coordinate system of $\alpha, \delta$. In -this spherical coordinate system, the coordinate of longitude is -labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the -latitude, measured from the pole, is labeled $\theta$ and has domain -$0 \le \theta \le \pi$. The coordinate frame of $\phi,\theta$ is -defined so that $\phi_p$, the longitude of the target system pole, is -0.0. - -For an arbitrary projection center, it is necessary to convert the -spherical coordinates to be projected ($\alpha,\delta$) to the -projection spherical coordinate system coordinates ($\phi, \theta$). -In practice, we construct the following useful trigonometric -relationships between $\phi$ and $\theta$ which may be employed in the -equations of $x,y$ below: -% -\begin{eqnarray} -\sin \theta & = & \sin \delta \sin \delta_p + \cos \delta \cos \delta_p \cos (\alpha - \alpha_p) \\ -\cos \theta \cos \phi & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \cos (\alpha - \alpha_p) \\ -\cos \theta \sin \phi & = & - \cos \delta \sin (\alpha - \alpha_p) -\end{eqnarray} -% -For the inverse transformations, the equivalent relationships are: -% -\begin{eqnarray} -\sin \delta & = & \sin \theta \sin \delta_p + \cos \theta \cos \delta_p \cos \phi \\ -\cos \delta \cos (\alpha - \alpha_p) & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \cos \phi \\ -\cos \delta \sin (\alpha - \alpha_p) & = & - \cos \theta \sin \phi -\end{eqnarray} -% -For zenithal projections, the linear coordinates are related to -$\phi,\theta$ by: -% -\begin{eqnarray} -x & = & R_\theta \sin \phi \\ -y & = & -R_\theta \cos \phi -\end{eqnarray} -% -and the inverse: -% -\begin{eqnarray} -R_\theta & = & \sqrt{x^2 + y^2} \\ -\phi & = & {\rm atan} (-y,x) -\end{eqnarray} -% -The coordinates $x,y$ above are defined to be in angular units (ie, -radians). - -From these relationships, we can calculate $\alpha, \delta$ as: -% -\begin{eqnarray} -\alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\ -\delta & = & \arcsin (\sin \delta) \\ -\end{eqnarray} -% -Note that if $(x,y) = (0,0)$, then $\alpha = \alpha_p, \delta = \delta_p$. - -\subparagraph{Gnomonic} - -The Gnomonic projection (``TAN'') is a zenithal projection with -$R_\theta = \cot \theta$. The resulting relationships for $(x,y)$ and -for $\sin \theta, \cos \theta$ are: - -\begin{eqnarray} -x & = & \frac{\cos \theta \sin \phi}{\sin \theta} \\ -y & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\ -\sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\ -\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\ -\end{eqnarray} - -where $\zeta = 1 / R_\theta$. - -\subparagraph{Orthographic} - -The Orthographic projection (``SIN'') is a zenithal projection with -$R_\theta = \cos \theta$. The resulting relationships for $(x,y)$ and -for $\sin \theta, \cos \theta$ are: - -\begin{eqnarray} -x & = & \cos \theta \sin \phi \\ -y & = & -\cos \theta \cos \phi \\ -\sin \theta & = & \sqrt{1 - R_\theta^2} \\ -\cos \theta & = & R_\theta \\ -\end{eqnarray} - -\paragraph{Cylindrical and Pseudocylindrical Projections} - -The {\em cylindrical} and {\em pseudocylindrical} projections are -defined relative to a set of cylindrical coordinates whose pole is -coincident with the pole of the spherical coordinates. These -projections are particularly used for full-sky representations, and -are only defined for projection centers with $\delta_p = 0$. In this -spherical coordinate system, the coordinate of longitude is labeled -$\phi$, and has domain of $-\pi < \phi \le \pi$, while the latitude, -measured from the pole, is labeled $\theta$ and has domain $0 \le -\theta \le \pi$. The projection center longitude, $\alpha_p$ -corresponds to $\phi = 0$, thus the value of $\phi$ is determined as -$\alpha - \alpha_p$ for all such projections. - -\subparagraph{Cartesian} - -The Cartesian projection (``CAR'') is a very simple cylindrical -projection with the following relationships between $x,y$ and -$\phi,\theta$: - -\begin{eqnarray} -x & = & \phi \\ -y & = & \theta -\end{eqnarray} - -\subparagraph{Mercator} - -The Mercator projection (``MER'') is a cylindrical projection. - -\begin{eqnarray} -x & = & \phi \\ -y & = & \ln \left( \tan (\pi/4 + \theta/2) \right) \\ -{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^y \right) - \pi/2 -\end{eqnarray} - -\subparagraph{Hammer-Aitoff} - -The Hammer-Aitoff projection(``AIT'') is a pseudocylindrical projection, and is defined: - -\begin{eqnarray} -x & = & 2 \zeta \cos \theta \sin \frac{\phi}{2} \\ -y & = & \zeta \sin \theta \\ -{\rm where}\hspace{1cm} \zeta^{-1} & \equiv & \sqrt{\frac{1}{2}\left(1 + \cos \theta \cos \frac{\phi}{2} \right)} -\end{eqnarray} - -And in reverse: - -\begin{eqnarray} -\phi & = & 2 {\rm \arctan} (2z^2 - 1, x z) \\ -\theta & = & \arcsin (yz) \\ -{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x/2)^2 - y^2} -\end{eqnarray} - -\subparagraph{Parabolic} - -The Parabolic projection (``PAR'') is a pseudocylindrical projection, and is defined: - -\begin{eqnarray} -x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\ -y & = & \pi \sin \frac{\theta}{3} \\ -\end{eqnarray} - -And in reverse: - -\begin{eqnarray} -\theta & = & 3 \sin^{-1} \rho \\ -\phi & = & \frac{x}{1 - 4\rho^2} \\ -{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\ -\end{eqnarray} - -\subsubsection{Offset} - -Coordinate offsets can be either spherical offsets or linear offsets. - -A spherical offset is performed by adding the components of the -offset, after unit conversion, to the given position. The resulting -coordinates must be wrapped to within the allowed range ($-\pi$ to -$\pi$, 0 to $2\pi$). - -A linear offset is defined to be a linear offset in a tangent -projection centered on the starting coordinate with $y$ axis aligned -with the local direction or increasing Declination. This projection -is undefined only for the coordinates exactly at the north and south -poles, in which case the orientation is defined to have the $y$ axis -parallel to the line of RA = 0.0. The scale of the projection is 1.0 -(ie, 1 'pixel' is 1 radian) and the given offsets must the scaled -based on the given offset units. - -Pseudo-code to implement the above for an offset: - -\begin{verbatim} -psSphere *psSphereSetOffset (psSphere pos, psSphere offset) { - - psPlane lin; - psSphere new; - psProjection proj; - - proj.R = pos->r; - proj.D = pos->d; - proj.X = 0; - proj.Y = 0; - proj.type = PS_PROJ_TAN; - - lin.x = offset.r; - lin.y = offset.d; - - new = psDeproject (&lin, &proj); - return (new); -} -\end{verbatim} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsubsection{Tangent Plane to Sky} + +\subsection{Tangent Plane to Sky} \tbd{we will replace the SLALIB version of AOPPA with a new function} @@ -1798,5 +1695,596 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{The One-to-Many Problem with Mosaic Cameras} +\subsection{Sky to Tangent Plane (II)} + +This section describes the transformation between celestial coordinates +(R.A., Dec.) and local terrestrial coordinates (Az, Alt). This transformation +is broken down into a number of steps as described below. + +\paragraph{Reference Implementations} + +There are two reference implementatins for the code to account for the +motion of the Earth in space. The first are the sample routines +provided by the IERS to accompany chaper 5 of IERS Bulletin 32. This +document and the code can be downloaded from +http://maia.usno.navy.mil/conv2003.html . The second reference +implementation is the SOFA software package managed by the IAU and +available at http://www.iau-sofa.rl.ac.uk Only the 2003-04-29 version +of SOFA should be used. The IERS code requires a few of the rotation +matrix utility routines from SOFA. + +Both implementations are in FORTRAN 77. The SOFA code has a more +complex implementation of precession-nutation for backward +compatibility with the pre 2003-01-01 conventions. The IERS code +includes some tricks to achieve greater precision in the fundamental +arguments of nutation, which the SOFA code omits. Therefore, the main +reference for psLib should be the IERS code. Note that the IERS code +calculates the transform from terrestrial to celestial coordinates, +while the SOFA code calculates its inverse. + +\subsubsection{Coordinate Systems} + +\begin{figure} +\psfig{file=transforms.ps} +\caption{Coordinates systems and the transformations between them} + +\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. + +\paragraph{ICRS} +The official IAU-sanctioned celestial coordinate system is the +International Celestial Reference System (ICRS). It is defined in terms of +a number of radio sources whose positions have been measured using VLBI. +It can be tied to the optical through the Hipparcos catalog. The ICRS has its +origin at the solar system barycenter. + +\paragraph{GCRS} +The Geocentric Celestial Reference System (GCRS) corresponds to the ICRS, but +has its origin at the center of the Earth. The differences between the two +systems are due to the velocity of the Earth (aberration), the position of +the Earth (parallax), and general relativistic bending of light rays. +There is no net rotation between the ICRS and the GCRS. + +\paragraph{ITRS} +The International Terrestrial Reference System (ITRS) is a coordinate +system which is fixed with respect to the Earth's crust. + +\paragraph{Intermediate Coordinate Systems - CIP, CEO, TEO} +The transform between the GCRS and ITRS is conventionally +decomposed into three parts in order to isolate the relatively rapid rotation +of the Earth from the movement of the Earth's rotational axis in the GCRS +and ITRS. All three sub-transforms are rigid rotations. + +This decomposition results in two intermediate coordinate systems. Both of +these share the same pole, known as the Celestial Intermediate Pole (CIP). +The CIP is defined by its motion in the GCRS to match the Tisserand +mean axis of the Earth (Seidelmann 1982, Celesial Mechanics 27, 78-106), +excluding motions with periods less than or equal +to two days. The CIP approximates the angular momentum vector of the +rotating Earth. + +The X axes of the intermediate coordinate systems are known as the +Celestial and Terrestrial Ephemeris Origins. (CEO and TEO). Both are defined +to be non-rotating origins. A non-rotating origin is a point on the equator +whose instantaneous motion is always orthogonal to the equator +(Kaplan 2003 IAU XXV Joint Discussion 16 +\footnote{http://aa.usno.navy.mil/kaplan/NROs\%5BJD16proc\%5D.pdf}). +Thus the CEO is defined by its position in the GCRS at some epoch and by the +motion of the CIP in the GCRS since that date. Similarly the TEO is +defined by its position in the ITRS at some epoch and the motion of the +CIP in the ITRS since that date. + +\subsubsection{ICRS - GCRS} + +The transformation between barycentric (ICRS) and geocentric (GCRS) coordinates +involves two components. These are +the general relativistic deflection of light rays by the Sun's gravity, and +aberration, due to the orbital motion +of the Earth. + +\paragraph{Gravitational Deflection} + +The Sun's gravity bends the path of light rays which pass near it. +To first order, a light ray is deflected by an angle of $4GM/c^2r_0$ radians, +where $G$ is the gravitational constant, +$M$ is the mass of the Sun, +$c$ is the speed of light, and +$r_0$ is the point of closest approach to the light ray to the Sun. +To the same order this is equal to the impact parameter - i.e. the point +of closest approach if the light ray were not deflected. Note that +$r_0/d = \tan(\theta)$, where $d$ is the distance from the Earth +to the Sun, and $\theta$ +is the angular separation of the star from the center of the Sun. + +There is a maximum deflection of 1.75 arc seconds if we set +$r_0$ to the radius of the sun. +Since the Sun bends light rays toward it, a star appears shifted away from the sun in the sky. + +\paragraph{Aberration} + +Aberration is the apparent change in direction of a ray of light in the +reference frame of a moving observer. Traditionally the aberration +calculation has been done with a linear expansion of the full +relativistic expression, often neglecting all but the linear term in +$v/c$, since the relativistic terms are on the order of a miliarcsecond. +However, the full relativistic expression poses no challenge for modern +computers, so psLib will use the following procedure to calculate aberration. + +Suppose an observer has a velocity $\beta\hat{\beta}$, with respect to +the Solar System barycenter, where $\beta$ is in units of the speed of +light, and $\hat{\beta}$ is a unit vector. Suppose also that the unit vector +$\hat{r}$ points toward a star in the barycenter frame of reference +(i.e. the ``actual'' position). +and $\hat{r}'$ gives the direction of the star in the observer's frame, +(i.e. the apparent position). + +First, decompose $\hat{r}$ into components parallel and perpendicular to +$\hat{\beta}$ by calculating +$\mu = \hat{r}\cdot\hat{\beta}$ and +$\vec{r}_\perp = \hat{r} - \mu \hat{\beta}$. + +Next, use the following expression for relativistic beaming, modified +slightly from equation 4.8b of Rybicki and Lightman: +\begin{equation} +\mu' = \mu + \beta \frac{\mu^2 - 1}{1 - \beta\mu} +\end{equation} +where $\mu' = \hat{r}' \cdot \hat{\beta}$. + +Now, the component of $\hat{r}'$ perpendicular to $\hat{\beta}$ +(i.e. $\vec{r}_\perp'$) must point +in the same direction as $\vec{r}_\perp$, but will have a different magnitude +because $\hat{r}'$ is a unit vector. In other words, +$\vec{r}_\perp' = a\vec{r}_\perp$, for some scalar $a$. So the next step is +to calculate $a = \sqrt{(1-\mu'^2)/\vec{r}_\perp}$. + +Finally, reassemble the components of +$\hat{r}' = \mu'\hat{\beta} + a \vec{r_\perp}$. + + +\subsubsection{GCRS - ITRS} +The transformation between geocentric celestial coordinates and terrestrial +coordinates is a solid body rotation due to the motion of the Earth is space. +This is conventionally broken down into three components to isolate the +relatively rapid rotation of the Earth from the motion of its rotational axis. + +This section is largely a summary of +Chapter 5 of IERS Technical Note 32 \footnote{http://maia.usno.navy.mil/conv2003.html} +(hereafter IERS32), +which is a description of the implementation of the Resoltions of the +XXIVth General Assembly of the IAU, available from the same URL as above. +These two documents describe a set of conventions which have been in effect +since 2003-01-01. The conventions in effect before that date will not be +implemented by psLib. + + +\paragraph{Precession/Nutation} + +The transform between the GCRS and the CIP/CEO coordinate systems is described +by the IAU 2000A precession-nutation model, which is accurate to the +0.2 mas level. +For higher accuracy the user must apply corrections to the model, which are tabulated by the IERS. + + + +The IAU 2000A precession-nutation model may be calculated in the following +way. First calculate the time $t$ as the number of Julian centuries since +2000-01-01T12:00:00 TT. + +Next calculate the fundamental arguments of nutation using equations (40) +and (41) of IERS32, reproduced below: +\begin{eqnarray} +F_1\equiv l\quad =~&\ Mean\ Anomaly\ of\ the\ Moon \cr + =~& 134.96340251^\circ + 1717915923.2178'' t + + 31.8792'' t^2 + 0.051635'' t^3 - 0.00024470'' t^4,\cr +F_2\equiv l'\quad =~&\ Mean\ Anomaly\ of\ the\ Sun\cr +=~& 357.52910918^\circ + 129596581.0481'' t +- 0.5532'' t^2 + 0.000136'' t^3 - 0.00001149'' t^4,\cr +F_3\equiv F\quad =~& L - \Omega\cr +=~& 93.27209062^\circ + 1739527262.8478'' t - 12.7512'' t^2 +- 0.001037'' t^3 + 0.00000417'' t^4,\cr +F_4\equiv D\quad =~&\ Mean\ Elongation\ of\ the\ Moon\ from\ the\ Sun\cr +=~& 297.85019547^\circ + 1602961601.2090'' t - 6.3706'' t^2 ++ 0.006593'' t^3 - 0.00003169'' t^4,\cr +F_5\equiv\Omega\quad =~&\ Mean\ Longitude\ of\ the\ Ascending\ Node\ of\ +the\ Moon\cr +=~& 125.04455501^\circ - 6962890.5431'' t + 7.4722'' t^2 + 0.007702'' t^3 - 0.00005939'' t^4 \cr +F_6\ \equiv l_{Me}\quad =~& 4.402 608 842 + 2608.7903 141 574\times t,\cr +F_7\ \equiv l_{Ve}\quad =~& 3.176 146 697 + 1021.3285 546 211 \times t,\cr +F_8\ \equiv l_{E\ }\quad =~& 1.753 470 314 + 628.3075 849 991 \times t,\cr +F_9\equiv l_{Ma}\quad =~& 6.203 480 913 + 334.0612 426 700 \times t,\cr +F_{10}\equiv l_{Ju}\quad =~& 0.599 546 497 + 52.9690 962 641 \times t,\cr +F_{11}\equiv l_{Sa}\quad =~& 0.874 016 757 + 21.3299 104 960 \times t,\cr +F_{12}\equiv l_{Ur}\quad =~& 5.481 293 872 + 7.4781 598 567 \times t,\cr +F_{13}\equiv l_{Ne}\quad =~& 5.311 886 287 + 3.8133 035 638 \times t,\cr +F_{14}\equiv p_{a\ }\quad =~& 0.024 381 750 \times t + 0.000 005 386 91 \times t^2. +\end{eqnarray} + +Next calculate the quantities $X$, $Y$, and $s$, using expressions of the form: + +\begin{equation} + \sum_{j} p_j t^j + \sum_{j}\sum_{i}[ + (a_{{\rm s},j})_i t^j \sin ({\rm \scriptstyle {ARG_{i,j}}}) + + (a_{{\rm c},j})_i t^j \cos ({\rm \scriptstyle {ARG_{i,j}}})] + , +\end{equation} + +where the $\rm \scriptstyle{ARG_{i,j}} = \sum_{k} w_{i,j,k} F_k$ represent linear +combinations of the fundamental arguments of nutation. + +The constants $p_j$, $w_{i,j,k}$, $(a_{{\rm s},j})_i$, and $(a_{{\rm c},j})_i$ +are given in the ASCII files: +tab5.2a.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2a.txt} (for $X$), +tab5.2b.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2b.txt} (for $Y$), and +tab5.2c.txt \footnote{http://maia.usno.navy.mil/conv2000/chapter5/tab5.2c.txt} (for $s+XY/2$). +Note that the expansion is given for $s+XY/2$, since this series converges +more rapidly than the one for $s$ alone. + +Each file contains a human-readable header, which includes the polynomial +coeficients, $p_j$ under the heading ``Polynomial part''. The data part of the +file lists the remaining constants, with rows cycling first through $i$, and +then through $j$. There is a separate heading each time $j$ increments. +Each row contains the following columns: + +\begin{itemize} +\item col 1 - A running index of rows in the table. +\item col 2 - The sine coeficients, $(a_{{\rm s},j})_i$ +\item col 3 - The cosine coeficients, $(a_{{\rm c},j})_i$ +\item cols 4 - 17 The weighting factors for the fundamental arguments of + nutation, $w_{i,j,k}$. +\end{itemize} + + +A FORTRAN reference implementation for the precession/nutation model is available from the IERS +\footnote{http://maia.usno.navy.mil/conv2000/chapter5/XYS2000A.f}. +The psLib results should agree with the reference implementation to within +the limits of numerical precision. + +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. + +The final step is to use $X$, $Y$, and $s$ to calculate the rotation +matrix from the CIP/CEO system to the GCRS using IERS32 equation (10), +reproduced below: + +\begin{equation} +\begin{pmatrix}1-aX^2& -aXY& X\cr -aXY& 1-aY^2& Y\cr -X& -Y& +1-a(X^2+Y^2)\cr +\end{pmatrix} \cdot R_3(s), +\end{equation} +where $R_3$ denotes a rotation about the Z axis, +$a = 1/(1+\sqrt{1 - X^2 + Y^2})$, +and $X$ and $Y$ are expressed in radians. +A FORTRAN reference implementation for this calculation is given +by the IERS \footnote{http://maia.usno.navy.mil/conv2000/chapter5/BPN2000.f}. + +Note that above we gave the expression for the transform toward celestial +coordinates (upward in figure X), in order to match the IERS reference code. +The inverse transform may be found by inverting the resulting rotation. + +\paragraph{Rotation of the Earth} + +The transform from the CIP/CEO to CIP/TEO coordinate systems is a +rotation about the CIP (i.e. the Z axis) by an angle known as the +``Earth Rotation Angle''. +By definition the Earth Rotation Angle is given by +equation (13) of IERS32, reproduced below: +\begin{equation} +\theta(T_u)=2\pi(0.7790572732640 + 1.00273781191135448T_u), +\end{equation} +where $T_u$ is the Julian UT1 date minus 2451545.0 . + +\paragraph{Polar Motion} + +The motion of the CIP in the ITRS is known as ``polar motion''. Similarly to +precession/nutation, the instantaneous position of the CIP in the +ITRS is specified by the quantites $x_p$, and $y_p$, and a third quantity, +$s'$, 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 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'$. + +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. + +\subsubsection{ITRS - Alt/Az} + +\paragraph{Orientation of the Observer} + +An observer's astronomical longitude and latitude give the orientation of +the local vertical with respect to the ITRS. Note that these coordinates +can be approximated by the geographic longitude and latitude of the observatory, +but their exact values must be calibrated from observation of stars +with known coordinates in the ICRS. + +The transform from the ITRS to Az/Alt in the absence of atmospheric refraction +is first a rotation about the Z axis by the observer's astronomical longitude, +and then a rotation about the Y axis of 90 degrees minus the observer's +astronomical latitude, followed by a rotation about the Z axis of 180 degrees +so that North is zero azimuth. + +\paragraph{Atmospheric Refraction} + +\tbd{add in summary of Ken's paper} + +\subsection{Projections} + +We implement three types of projections: {\em zenithal}, {\em +cylindrical} and {\em pseudocylindrical}, each requiring slightly +different handling. Our representations are based on the treatment of +projections presented by +\href{http://www.cv.nrao.edu/fits/documents/wcs/wcs.all.ps}{Greisen \& +Calabretta (1995, ADASS, 4, 233)}. In all of these projections, we +are converting from a spherical coordinate $\alpha,\delta$ to a linear +(2-D) coordinate $x_p,y_p$. The projection is defined by the +projection type, the projection center ($\alpha_p, \delta_p$) and the +the plate scales in the $x_p$ and $y_p$ directions ($\rho_x,\rho_y$). + +In the structure, \code{psProjection}, the projection type is defined +by the element \code{type}, the projection center $\alpha_p,\delta_p$ +is defined by the elements \code{R,D}, and the plate scales, +$\rho_x,\rho_y$, are defined by the elements \code{Xs,Ys}. The plate +scales are applied independently to the $x$ and $y$ coordinates to +convert them to the corresponding linear units (ie, pixels): +% +\begin{eqnarray} +x_p & = & \rho_x x \\ +y_p & = & \rho_y y \\ +\end{eqnarray} +% +In the discussions below, we ignore this last step (or first step, +depending on the direction of the conversion). + +\subsubsection{Zenithal Projections} + +The {\em zenithal} projections are defined relative to a set of +spherical coordinates with pole at the center of the projection +($\alpha_p, \delta_p$), and which thus represents a coordinate system +rotated relative to the coordinate system of $\alpha, \delta$. In +this spherical coordinate system, the coordinate of longitude is +labeled $\phi$, and has domain of $-\pi < \phi \le \pi$, while the +latitude, measured from the pole, is labeled $\theta$ and has domain +$0 \le \theta \le \pi$. The coordinate frame of $\phi,\theta$ is +defined so that $\phi_p$, the longitude of the target system pole, is +0.0. + +For an arbitrary projection center, it is necessary to convert the +spherical coordinates to be projected ($\alpha,\delta$) to the +projection spherical coordinate system coordinates ($\phi, \theta$). +In practice, we construct the following useful trigonometric +relationships between $\phi$ and $\theta$ which may be employed in the +equations of $x,y$ below: +% +\begin{eqnarray} +\sin \theta & = & \sin \delta \sin \delta_p + \cos \delta \cos \delta_p \cos (\alpha - \alpha_p) \\ +\cos \theta \cos \phi & = & \sin \delta \cos \delta_p - \cos \delta \sin \delta_p \cos (\alpha - \alpha_p) \\ +\cos \theta \sin \phi & = & - \cos \delta \sin (\alpha - \alpha_p) +\end{eqnarray} +% +For the inverse transformations, the equivalent relationships are: +% +\begin{eqnarray} +\sin \delta & = & \sin \theta \sin \delta_p + \cos \theta \cos \delta_p \cos \phi \\ +\cos \delta \cos (\alpha - \alpha_p) & = & \sin \theta \cos \delta_p - \cos \theta \sin \delta_p \cos \phi \\ +\cos \delta \sin (\alpha - \alpha_p) & = & - \cos \theta \sin \phi +\end{eqnarray} +% +For zenithal projections, the linear coordinates are related to +$\phi,\theta$ by: +% +\begin{eqnarray} +x & = & R_\theta \sin \phi \\ +y & = & -R_\theta \cos \phi +\end{eqnarray} +% +and the inverse: +% +\begin{eqnarray} +R_\theta & = & \sqrt{x^2 + y^2} \\ +\phi & = & {\rm atan} (-y,x) +\end{eqnarray} +% +The coordinates $x,y$ above are defined to be in angular units (ie, +radians). + +From these relationships, we can calculate $\alpha, \delta$ as: +% +\begin{eqnarray} +\alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\ +\delta & = & \arcsin (\sin \delta) \\ +\end{eqnarray} +% +Note that if $(x,y) = (0,0)$, then $\alpha = \alpha_p, \delta = \delta_p$. + +\paragraph{Gnomonic} + +The Gnomonic projection (``TAN'') is a zenithal projection with +$R_\theta = \cot \theta$. The resulting relationships for $(x,y)$ and +for $\sin \theta, \cos \theta$ are: + +\begin{eqnarray} +x & = & \frac{\cos \theta \sin \phi}{\sin \theta} \\ +y & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\ +\sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\ +\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\ +\end{eqnarray} + +where $\zeta = 1 / R_\theta$. + +\paragraph{Orthographic} + +The Orthographic projection (``SIN'') is a zenithal projection with +$R_\theta = \cos \theta$. The resulting relationships for $(x,y)$ and +for $\sin \theta, \cos \theta$ are: + +\begin{eqnarray} +x & = & \cos \theta \sin \phi \\ +y & = & -\cos \theta \cos \phi \\ +\sin \theta & = & \sqrt{1 - R_\theta^2} \\ +\cos \theta & = & R_\theta \\ +\end{eqnarray} + +\subsubsection{Cylindrical and Pseudocylindrical Projections} + +The {\em cylindrical} and {\em pseudocylindrical} projections are +defined relative to a set of cylindrical coordinates whose pole is +coincident with the pole of the spherical coordinates. These +projections are particularly used for full-sky representations, and +are only defined for projection centers with $\delta_p = 0$. In this +spherical coordinate system, the coordinate of longitude is labeled +$\phi$, and has domain of $-\pi < \phi \le \pi$, while the latitude, +measured from the pole, is labeled $\theta$ and has domain $0 \le +\theta \le \pi$. The projection center longitude, $\alpha_p$ +corresponds to $\phi = 0$, thus the value of $\phi$ is determined as +$\alpha - \alpha_p$ for all such projections. + +\paragraph{Cartesian} + +The Cartesian projection (``CAR'') is a very simple cylindrical +projection with the following relationships between $x,y$ and +$\phi,\theta$: + +\begin{eqnarray} +x & = & \phi \\ +y & = & \theta +\end{eqnarray} + +\paragraph{Mercator} + +The Mercator projection (``MER'') is a cylindrical projection. + +\begin{eqnarray} +x & = & \phi \\ +y & = & \ln \left( \tan (\pi/4 + \theta/2) \right) \\ +{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^y \right) - \pi/2 +\end{eqnarray} + +\paragraph{Hammer-Aitoff} + +The Hammer-Aitoff projection(``AIT'') is a pseudocylindrical projection, and is defined: + +\begin{eqnarray} +x & = & 2 \zeta \cos \theta \sin \frac{\phi}{2} \\ +y & = & \zeta \sin \theta \\ +{\rm where}\hspace{1cm} \zeta^{-1} & \equiv & \sqrt{\frac{1}{2}\left(1 + \cos \theta \cos \frac{\phi}{2} \right)} +\end{eqnarray} + +And in reverse: + +\begin{eqnarray} +\phi & = & 2 {\rm \arctan} (2z^2 - 1, x z) \\ +\theta & = & \arcsin (yz) \\ +{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x/2)^2 - y^2} +\end{eqnarray} + +\paragraph{Parabolic} + +The Parabolic projection (``PAR'') is a pseudocylindrical projection, and is defined: + +\begin{eqnarray} +x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\ +y & = & \pi \sin \frac{\theta}{3} \\ +\end{eqnarray} + +And in reverse: + +\begin{eqnarray} +\theta & = & 3 \sin^{-1} \rho \\ +\phi & = & \frac{x}{1 - 4\rho^2} \\ +{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\ +\end{eqnarray} + +\subsection{Offset} + +Coordinate offsets can be either spherical offsets or linear offsets. + +A spherical offset is performed by adding the components of the +offset, after unit conversion, to the given position. The resulting +coordinates must be wrapped to within the allowed range ($-\pi$ to +$\pi$, 0 to $2\pi$). + +A linear offset is defined to be a linear offset in a tangent +projection centered on the starting coordinate with $y$ axis aligned +with the local direction or increasing Declination. This projection +is undefined only for the coordinates exactly at the north and south +poles, in which case the orientation is defined to have the $y$ axis +parallel to the line of RA = 0.0. The scale of the projection is 1.0 +(ie, 1 'pixel' is 1 radian) and the given offsets must the scaled +based on the given offset units. + +Pseudo-code to implement the above for an offset: + +\begin{verbatim} +psSphere *psSphereSetOffset (psSphere pos, psSphere offset) { + + psPlane lin; + psSphere new; + psProjection proj; + + proj.R = pos->r; + proj.D = pos->d; + proj.X = 0; + proj.Y = 0; + proj.type = PS_PROJ_TAN; + + lin.x = offset.r; + lin.y = offset.d; + + new = psDeproject (&lin, &proj); + return (new); +} +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{The One-to-Many Problem with Mosaic Cameras} The \PS{} focal plane consists of several chips, so we will often want @@ -1820,5 +2308,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{General Astronomy Functions} +\subsection{General Astronomy Functions} \tbd{we will provide a new airmass function} @@ -1854,5 +2342,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{Positions of Major Solar System Objects} +\subsection{Positions of Major Solar System Objects} \tbd{ephemerides code to replace this} @@ -1869,16 +2357,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Missing and Todo} - -\tbd{define sunrise, sunset, sun position} - -\tbd{define moonrise, moonset, moon position, moon phase} - -\tbd{define planet functions} - -\tbd{clean up FITS I/O issues} - -\tbd{define Brent's method \& minimization bracketing} - +\pagebreak \section{Pan-STARRS Modules} @@ -2225,4 +2702,17 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Missing and Todo} + +\tbd{define sunrise, sunset, sun position} + +\tbd{define moonrise, moonset, moon position, moon phase} + +\tbd{define planet functions} + +\tbd{clean up FITS I/O issues} + +\tbd{define Brent's method \& minimization bracketing} + \appendix \section{Change Log} @@ -2230,39 +2720,2 @@ \end{document} - - -\section{Modules} - -\subsection{Image Processing Modules} -\subsubsection{debias} -\subsubsection{mask} -\subsubsection{trim} -\subsubsection{flatten} -\subsubsection{sky/fringe subtract} -\subsubsection{warp} -\subsubsection{stack} -\subsubsection{difference} -\subsubsection{kernel convolution} -\subsubsection{special stack} - -\subsection{Object Detection Modules} -\subsubsection{find peaks} -\subsubsection{background } -\subsubsection{aperture photometry} -\subsubsection{get shape} - -\subsection{Miscellaneous Modules} - -\section{Analysis Stages} -\subsection{Phase 1} -\subsection{Phase 2} -\subsection{Phase 3} -\subsection{Phase 4} -\subsection{Cal 1} -\subsection{Cal 2} -\subsection{Cal 3} -\subsection{Astrom Ref} -\subsection{Photom Ref} - -\section{Architectual Components} -