Index: /trunk/doc/pslib/ChangeLogADD.tex =================================================================== --- /trunk/doc/pslib/ChangeLogADD.tex (revision 3563) +++ /trunk/doc/pslib/ChangeLogADD.tex (revision 3564) @@ -1,4 +1,4 @@ -\subsection{Changes from version 06 to version 07} +\subsection{Changes from version 06 (7 September 2004) to version 07 (24 November 2004)} \begin{itemize} @@ -12,5 +12,5 @@ \end{itemize} -\subsection{Changes from version 07 to version 08} +\subsection{Changes from version 07 (24 November 2004) to version 08 (8 January 2005)} \begin{itemize} @@ -20,5 +20,5 @@ \end{itemize} -\subsection{Changes from version 08 to version 09} +\subsection{Changes from version 08 (8 January 2005) to version 09 (14 February 2005)} \begin{itemize} @@ -27,4 +27,9 @@ \code{PS_RESAMPLE_LAGRANGE}. \item Added section on FITS WCS. +\end{itemize} + +\subsection{Changes from version 09 (14 February 2005) to version 10 (29 March 2005} + +\begin{itemize} \item Changes to the Time section: \begin{itemize} @@ -36,3 +41,11 @@ \item Verbatim Fortran code for UT1 interpolation \end{itemize} + +\item section reorganization: +\begin{itemize} +\item renamed Astronomical Image Manipulations to Image Manipulations +\item moved Image Manipulations to own subsection before Astronomy Utilities +\item promoted all PSLib subsections to sections \end{itemize} + +\end{itemize} Index: /trunk/doc/pslib/ChangeLogSDRS.tex =================================================================== --- /trunk/doc/pslib/ChangeLogSDRS.tex (revision 3563) +++ /trunk/doc/pslib/ChangeLogSDRS.tex (revision 3564) @@ -1,3 +1,3 @@ -%%% $Id: ChangeLogSDRS.tex,v 1.77 2005-03-29 03:42:16 price Exp $ +%%% $Id: ChangeLogSDRS.tex,v 1.78 2005-03-30 21:14:48 eugene Exp $ \subsection{Changes from version 00 to version 01} @@ -464,5 +464,5 @@ \end{itemize} -\subsection{Changes from Revision 12 (9 February 2005) to Revision 13 (Present)} +\subsection{Changes from Revision 12 (9 February 2005) to Revision 13 (30 March 2005)} \begin{itemize} @@ -481,10 +481,12 @@ from a file. \item \code{psMetadataAddV} changed to use \code{va_list} parameter (bug 312). -\item Added \code{psPixels} structure, \code{psPixelsToMask}, - \code{psMaskToPixels}, \code{psPixelsConcatenate} and - \code{psPixelsTransform}. + \item Modified \code{psImageTransform} in preparation for image combination. + +\item Added \code{psPixels} structure and related functions \item Added \code{psPlaneTransformDeriv}. \item Added \code{psImageGrowMask}. +\item Added Earth Orientation Calculations Section + \item Changes to the Time section: \begin{itemize} @@ -508,2 +510,13 @@ \item Adding logical operations (and, or) to \code{psBinaryOp}. \end{itemize} + +\item Substantial reorganization: +\begin{itemize} +\item Moved Metadata, Database, and XML sections to new section +\item Re-named Detector \& Sky Coordinates to Linear \& Spherical Coordinates +\item Moved Exposure and Observatory information out of 'Astronomical Images' +\item Moved Celestial Coordinate Systems out of 'Detector \& Sky Coordinates' +\item Added Atmospheric Effects section, incorporating psGrommit and airmass functions from other sections) +\item Moved Fixed Pattern out of Astronomical Images +\end{itemize} +\end{itemize} Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 3563) +++ /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} - Index: /trunk/doc/pslib/psLibSDRS.tex =================================================================== --- /trunk/doc/pslib/psLibSDRS.tex (revision 3563) +++ /trunk/doc/pslib/psLibSDRS.tex (revision 3564) @@ -1,3 +1,3 @@ -%%% $Id: psLibSDRS.tex,v 1.191 2005-03-29 03:42:21 price Exp $ +%%% $Id: psLibSDRS.tex,v 1.192 2005-03-30 21:14:48 eugene Exp $ \documentclass[panstarrs,spec]{panstarrs} @@ -11,7 +11,8 @@ \project{Pan-STARRS Image Processing Pipeline} \organization{Institute for Astronomy} -\version{12} +\version{13} \docnumber{PSDC-430-007} +% \setcounter{tocdepth}{5} % lowest level to be included in toc \setlength{\topsep}{-2pt} @@ -137,7 +138,6 @@ \href{heasarc.gsfc.nasa.gov/docs/software/fitsio}{\tt heasarc.gsfc.nasa.gov/docs/software/fitsio} -\item \tbd{SLALIB support is likely to be dropped} Many of the -astronomy routines will wrap the StarLink Positional Astronomy -libraries (SLALib): +\item the StarLink Positional Astronomy libraries (SLALib) are a + useful reference: \href{star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{\tt @@ -183,4 +183,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\pagebreak \section{System Utilities} @@ -1193,4 +1194,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\pagebreak \section{Basic Data Types and Collections} @@ -1894,4 +1896,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\pagebreak \section{Data manipulation} @@ -3028,4 +3031,86 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsection{Image Pixel Lists} + +Usually an image mask is the best way to carry information about what +pixels mean what. However, in the case where the number of pixels in +which we are interested is limited, it is more efficient to simply +carry a list of pixels. An example of this is in the image +combination code, where we want to perform an operation on a +relatively small fraction of pixels, and it is inefficient to go +through an entire mask image checking each pixel. + +\begin{verbatim} +typedef struct { + psVector *x; // x coordinate + psVector *y; // y coordinate +} psPixels; +\end{verbatim} + +Of course, the size of each of the vectors should match. In the event +that they do not match, any function which detects the problem shall +generate a warning and use the size of the shorter of the vectors as +the size. The order in which the pixels are kept is not considered +important. + +\begin{verbatim} +psImage *psPixelsToMask(psImage *out, const psPixels *pixels, const psRegion *region, unsigned int maskVal); +psPixels *psMaskToPixels(psPixels *out, const psImage *mask, unsigned int maskVal); +\end{verbatim} + +\code{psPixelsToMask} shall return an image of type U8 with the +\code{pixels} lying within the specified \code{region} set to the +\code{maskVal}. The \code{out} image shall be modified if supplied, +or allocated and returned if \code{NULL}. The size of the output +image shall be \code{region->x1 - region->x0} by \code{region->y1 - +region->y0}, with \code{out->x0 = region->x0} and \code{out->y0 = +region->y0}. In the event that either of \code{pixels} or +\code{region} are \code{NULL}, the function shall generate an error +and return \code{NULL}. + +\code{psMaskToPixels} shall return a \code{psPixels} consisting of the +coordinates in the \code{mask} that match the \code{maskVal}. The +\code{out} pixel list shall be modified if supplied, or allocated and +returned if \code{NULL}. In hte event that \code{mask} is +\code{NULL}, the function shall generate an error and return +\code{NULL}. + +\begin{verbatim} +psPixels *psPixelsConcatenate(psPixels *out, const psPixels *pixels); +\end{verbatim} + +\code{psPixelsConcatenate} shall concatenate \code{pixels} onto +\code{out}. In the event that \code{out} is \code{NULL}, a new +\code{psPixels} shall be allocated, and the contents of \code{pixels} +simply copied in. If \code{pixels} is \code{NULL}, the function shall +generate an error and return \code{NULL}. The function shall take +care to ensure that there are no duplicate pixels in \code{out} (since +the order in which the pixels are stored is not important, the values +may be sorted, allowing the use of a faster algorithm than a linear +scan). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Image Regions} + +In many places, we need to refer to a rectangular area. We define a +structure to represent a rectangle: +\begin{verbatim} +typedef struct { + float x0; + float x1; + float y0; + float y1; +} psRegion; +psRegion *psRegionAlloc (float x0, float x1, float y0, float y1); +\end{verbatim} + +\begin{verbatim} +psRegion *psRegionFromString (char *region); +\end{verbatim} +This function converts the IRAF description of a region in the form +\code{[x0:x1,y0:y1]}, used for header entries such as \code{BIASSEC}, +into the corresponding \code{psRegion} structure. + \subsection{Vector and Image Arithmetic} \label{sec:arithmetic} @@ -3443,4 +3528,1002 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\pagebreak +\section{Rich Data Structures and I/O} + +\subsection{Metadata} +\label{sec:metadata} + +\subsubsection{Conceptual Overview} + +Within PSLib, we provide a data structure to carry metadata and +mechanisms to manipulate the metadata. Metadata is a general concept +that requires some discussion. In any data analysis task, the +ensemble of all possible data may be divided into two or three +classes: there is the specific data of interest, there is data which +is related or critical but not the primary data of interest, and there +is all of the other data which may or may not be interesting. For +example, consider a simple 2D image obtained of a galaxy from a CCD +camera on a telescope. If you want to study the galaxy, the specific +data of interest is the collection of pixels. There are a variety of +other pieces of data which are closely related and crucial to +understanding the data in those pixels, such as the dimensions of the +image, the coordinate system, the time of the image, the exposure +time, and so forth. Other data may be known which may be less +critical to understanding the image, but which may be interesting or +desired at a later date. For example, the observer who took the +image, the filter manufacturer, the humidity at the telescope, etc. + +Formally, all of the related data which describe the principal data of +interest are metadata. Note that which piece is the metadata and +which is the data may depend on the context. If you are examining the +pixels in an image, the coordinate and flux of an object may be part +of the metadata. However, if you are analyzing a collection of +objects extracted from an image, you may consider then pixel data +simply part of the metadata associated with the list of objects. + +There are various ways to handle metadata vs data within a programming +environment. In C, it is convenient to use structures to group +associated data together. One possibility is to define the metadata +as part of the associated data structure. For example, the image data +structure would have elements for all possible associated measurement. +This approach is both cumbersome (because of the large number metadata +types), impractical (because the full range of necessary metadata is +difficult to know in advance) and inflexible (because any change in +the collection of metadata requires addition of new structure elements +and recompilation). + +An alternative is to place the metadata in a generic container and use +lookup mechanisms to extract the appropriate metadata when needed. An +example of this is would be a text-based FITS header for an image read +into a flat text buffer. In this implementation, metadata lookup +functions could return the current value of, for example, NAXIS1 (the +number of columns of the image) by scanning through the header buffer. +This method has the benefits of flexibility and simplicity of +programming interface, but it has the disadvantage that all metadata +is accessed though this lookup mechanism. This may make the code less +readable and it may slow down the access. + +PSLib implements an intermediate solution to this problem. We specify +a flexible, generic metadata container and access methods. Data types +which require association with a general collection of metadata should +include an entry of this metadata type. However, a subset of metadata +concepts which are basic and frequently required may be placed in the +coded structure elements. This approach allows the code to refer to +the basic metadata concepts as part of the data structure (ie, +\code{image.nx}), but also allows us to provide access to any +arbitrary metadata which may be generated. As a practical matter, the +choice of which entries are only in the metadata and which are part of +the explicit structure elements is rather subjective. Any data +elements which are frequently used should be put in the structure; +those which are only infrequently needed should be left in the generic +metadata. + +There are some points of caution which must be noted. Any +manipulation of the data should be reflected in the metadata where +appropriate. This is always an issue of concern. For example, +consider an image of dimensions \code{nx, ny}. If a function extracts +a subraster, it must change the values of \code{nx, ny} to match the +new dimensions. What should it do to the corresponding metadata? +Clearly, it should change the corresponding value which defines +\code{nX, nY}. However, it is not quite so simple: there may be other +metadata values which depend on those values. These must also be +changed appropriately. What if the metadata element points to a +copy of the metadata which may be shared by other representations of +the image? These must be treated differently because the change would +invalidate those other references. Care must be taken, therefore, +when writing functions which operate on the data to consider all of +the relevant metadata entries which must also be updated. + +A related issue is the definition of metadata names. Entries in a +structure have the advantage of being hardwired: every instance of +that structure will have the same name for the same entry. This is +not necessarily the case with a more flexible metadata container. The +image exposure time is a notorious example in astronomy. Different +observatories use different header keywords (ie, metadata names) for +the same concept of the exposure time (\code{EXPTIME}, +\code{EXPOSURE}, \code{OPENTIME}, \code{INTTIME}, etc). Any system +which operates on these metadata needs to address the issue of +identifying these names. This issue seems like an argument for +hardwiring metadata in the structure, but in fact it does not present +such a strong case. If the metadata are hardwired, some function will +still have to know how to interpret the various names to populate the +structure. The concept can still be localized with generic metadata +containers by including abstract metadata names within the code which +are tied to the various implementations-specific metadata names. + +\subsubsection{Metadata Representation} + +\begin{figure} +\psfig{file=Metadata,width=6.5in} +\caption{Metadata Structures\label{fig:metadata}} +\end{figure} + +This section addresses the question of how \PS{} metadata should be +represented in memory, not how it should be represented on disk. + +We define an item of metadata with the following structure: +\filbreak +\begin{verbatim} +typedef struct { + int id; ///< unique ID for this item + char *name; ///< Name of item + psMetadataType type; ///< type of this item + psElemType ptype; ///< primitive data type + const union { + psS32 S32; ///< integer data + psF32 F32; ///< floating-point data + psF64 F64; ///< double-precision data + void *V; ///< other type + psList *list; ///< psList entry + psMetadata *md; ///< psMetadata entry + } data; ///< value of metadata + char *comment; ///< optional comment ("", not NULL) +} psMetadataItem; +\end{verbatim} + +The \code{id} is a unique identifier for this item of metadata; +experience shows that such tags are useful. The entry \code{name} +specifies the name of the metadata item. The value of the metadata is +given by the union \code{data}, and may be of type \code{psS32}, +\code{psF32}, \code{psF64}, or an arbitrary rich structure pointed at +by the \code{void} pointer \code{V}. A character string comment +associated with this metadata item may be stored in the element +\code{comment}. The \code{type} entry specifies how to interpret the +type of the data being represented, given by the enumerated type +\code{psMetadataType}: +% +\filbreak +\begin{verbatim} +typedef enum { ///< type of item.data is: + PS_META_PRIMITIVE, ///< primitive type: use item.ptype + PS_META_LIST, ///< psList; use item.data.list (used for non-unique data) + PS_META_META, ///< psMetadata: use item.data.list + PS_META_STR, ///< string (item.data.V) + PS_META_MATH, ///< psScalar, psVector, psImage (item.data.V) + PS_META_JPEG, ///< JPEG (item.data) + PS_META_PNG, ///< PNG (item.data) + PS_META_ASTROM, ///< astrometric coefficients (item.data) + PS_META_UNKNOWN, ///< other (item.data) + PS_META_NTYPE ///< Number of types; must be last +} psMetadataType; +\end{verbatim} +If the data is a PSLib primitive data value, the primitive data type +is given by the value of \code{ptype}. + +A collection of metadata is represented by the \code{psMetadata} structure: +\begin{verbatim} +typedef struct { + psList *list; ///< list of psMetadataItem + psHash *table; ///< hash table of the same metadata +} psMetadata; +\end{verbatim} +The type \code{psMetadata} is a container class for metadata. Note +that there are in fact \emph{two} representations of the metadata +(each \code{psMetadataItem} appears on both). The first +representation employs a doubly-linked list that allows the order of +the metadata to be preserved (e.g., if FITS headers are read in a +particular order, they should be written in the same order). The +second representation employs a hash table which allows fast look-up +given a specific metadata keyword. + +Certain metadata names (such as the FITS keywords \code{COMMENT} and +\code{HISTORY} in a FITS header) may be repeated with different +values. In such a case, the \code{psMetadata.list} structure contains +the entries in their original sequence with duplicate keys. The +\code{psMetadata.hash} entries, which are required to have unique +keys, would have a single entry with the keyword of the repeated key, +with the value of \code{psMetadataType} set to \code{PS_META_LIST}, +and the \code{psMetadataItem.data} element pointing to a \code{psList} +containing the actual entries. If \code{psMetadataItemAlloc} is +called with the type set to \code{PS_META_LIST}, such a repeated key +is created. If the data value passed to \code{psMetadataItemAlloc} +(the quantity in ellipsis) is \code{NULL}, then an empty +\code{psMetadataItem} with the given keyword is created to hold future +entries of that keyword. + +The \code{psMetadataAdd} routine is required to check that all +metadata names are unique unless the type is already qualified as +\code{PS_META_LIST}; in this case the data are added to the +corresponding \code{psMetadataItem.data} list. + +\subsubsection{Metadata APIs} + +The allocator for \code{psMetadataItem} returns a full +\code{psMetadataItem} ready for insertion into the \code{psMetadata}. +The \code{name} entry specifies the name to use for this metadata +item, and may include \code{sprintf}-type formating codes. The +\code{comment} entry is a fixed string which is used for the comment +associated with this metadata item. The metadata data and the +arguments to the \code{name} formatting codes are passed, in that +order (metadata pointer first), to \code{psMetadataItemAlloc} as +arguments following the comment string. The data must be a pointer +for any data types which are stored in the element \code{data.void}, +while other data types are passed as numeric values. The argument +list must be interpreted appropriately by the \code{va_list} operators +in the function. +\begin{verbatim} +psMetadataItem *psMetadataItemAlloc(const char *name, psMetadataType type, const char *comment, ...); +psMetadataItem *psMetadataItemAllocV(const char *name, psMetadataType type, const char *comment, va_list list); +\end{verbatim} + +The constructor for the collection of metadata, \code{psMetadata}, +simply returns an empty metadata container (employing the constructors +for the doubly-linked list and hash table). The destructor needs to +free each of the \code{psMetadataItem}s using \code{psMetadataItemFree}. +\begin{verbatim} +psMetadata *psMetadataAlloc(void); +\end{verbatim} + +Items may be added to the metadata in one of two ways --- firstly, an +item may be added by appending a \code{psMetadataItem} which has +already been created; and secondly by directly providing the data to +be appended. In both cases, the return value defines the success +(\code{true}) or failure of the operation. The second function, +\code{psMetadataAdd} takes a pointer or value which is interpreted by +the function using variadic argument interpretation. The third +version is the \code{va_list} version of the second function. All +three functions take a parameter, \code{location}, which specifies +where in the list to place the element, following the conventions for +the \code{psList}. The entry \code{mode} for \code{psMetadataAddItem} +is a bit mask constructed by OR-ing the allowed option flags (eg, +\code{PS_META_REPLACE}) which specifies minor variations on the +behavior. The \code{format} entry, which specifies both the metadata +type and the optional flags, is constructed by bit-wise OR-ing the +appropriate \code{psMetadataType} and allowed option option flags. +Care should be taken not to leak memory when appending an item for +which the key already exists in the metadata (and is not +\code{PS_META_LIST}). +% +\begin{verbatim} +bool psMetadataAddItem(psMetadata *md, psMetadataItem *item, int location, int mode); +bool psMetadataAdd(psMetadata *md, int location, const char *name, int format, const char *comment, ...); +bool psMetadataAddV(psMetadata *md, int location, const char *name, int format, const char *comment, + va_list list); +\end{verbatim} + +The functions above take option flags which modify the behavior when +metadata items are added to the metadata list. These flags must be +bit-exclusive of those used above for the \code{psMetadataTypes}. The +flags have the following meanings: + +\code{PS_META_DEFAULT}: This is the zero bit mask, to allow the +default behavior for \code{psMetadataAddItem} above. If this is OR-ed +with a \code{psMetadataType}, the result is as if no OR-ing took +place. + +\code{PS_META_REPLACE}: If the given metadata item exists in the +metadata list, and is not of type \code{PS_META_LIST} or +\code{PS_META_META} (ie, not a container type), then this entry is +allowed to replace the existing entry. If this mode bit is not set, a +duplicate (non-container-type) entry is an error. + +\begin{verbatim} +typedef enum { ///< option flags for psMetadata functions + PS_META_DEFAULT, ///< default behavior (0x0000) for use in mode above + PS_META_REPLACE, ///< allow entry to be replaced +} psMetadataFlags; +\end{verbatim} + +An example of code to use these metadata APIs to generate the +structure seen in Figure~\ref{fig:metadata} is given below. + +\begin{verbatim} +md = psMetadataAlloc(); + +psMetadataAdd(md, PS_LIST_TAIL, "SIMPLE", PS_META_BOOL, "basic fits", TRUE); +psMetadataAdd(md, PS_LIST_TAIL, "BLANK", PS_META_S32, "invalid pixel data", -32768); +psMetadataAdd(md, PS_LIST_TAIL, "DATE-OBS", PS_META_STR, "observing date UT", " 2004-6-16"); +psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_LIST, "head of comment block", NULL); +psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "DATA"); +psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "PARAMS"); +psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME", PS_META_F32, "exposure time (sec)", 1.05); +psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "FOO"); + +cell = psMetadataAlloc(); +psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME", PS_META_STR, "", "CCD00"); +psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR, "", "BSEC-00"); +psMetadataAdd(cell, PS_LIST_TAIL, "CHIP", PS_META_STR, "", "CHIP.00"); +psMetadataAdd(md, PS_LIST_TAIL, "CELL.00", PS_META_META, "", cell); + +cell = psMetadataAlloc(); +psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME", PS_META_STR, "", "CCD01"); +psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR, "", "BSEC-01"); +psMetadataAdd(cell, PS_LIST_TAIL, "CHIP", PS_META_STR, "", "CHIP.01"); +psMetadataAdd(md, PS_LIST_TAIL, "CELL.01", PS_META_META, "", cell); +\end{verbatim} + +The following code shows how to use the APIs to replace one of these values: +\begin{verbatim} +psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME", PS_META_F32 | PS_REPLACE, "new exposure time (sec)", 2.05); +\end{verbatim} + +Items may be removed from the metadata by specifying a key or a +location in the list. If the value of \code{name} is \code{NULL}, the +value of \code{location} is used. If the value of \code{name} is not +\code{NULL}, then \code{location} must be set to +\code{PS_LIST_UNKNOWN}. If the key matches a metadata item, the item +is removed from the metadata and \code{true} is returned; otherwise, +\code{false} is returned. If the key is not unique, then \emph{all} +items corresponding to the key are removed, and \code{true} is +returned. +% +\begin{verbatim} +bool psMetadataRemove(psMetadata *md, int location, const char *key); +\end{verbatim} + +Items may be found within the metadata by providing a key. In the +event that the key is non-unique, the first item is returned. +\begin{verbatim} +psMetadataItem *psMetadataLookup(const psMetadata *md, const char *key); +\end{verbatim} + +Several utility functions are provided for simple cases. These +functions perform the effort of casting the data to the appropriate +type. The numerical functions shall return 0.0 if their key is not +found. If the pointer value of \code{status} is not \code{NULL}, it +is set to reflect the success or failure of the lookup. +\begin{verbatim} +void *psMetadataLookupPtr(bool *status, const psMetadata *md, const char *key); +psS32 psMetadataLookupS32(bool *status, const psMetadata *md, const char *key); +psF64 psMetadataLookupF64(bool *status, const psMetadata *md, const char *key); +\end{verbatim} + +Items may be retrieved from the metadata by their entry position. The +value of which specifies the desired entry in the fashion of +\code{psList}. +\begin{verbatim} +psMetadataItem *psMetadataGet(const psMetadata *md, int location); +\end{verbatim} + +The metadata list component may be iterated over by using a +\code{psListIterator} in a fashion equivalent to the usage for +\code{psList}. The iterator may be set to a location in the +\code{psMetadata} list, and the user may get the previous or next item +in the list relative to that location. \code{psMetadataGetNext} has +the ability to match the key using a POSIX regex, e.g., if the user +only wants to iterate through \code{IPP.machines.sky} and doesn't want +to bother with \code{IPP.machines.detector}. The iterator should +iterate over every item in the metadata list, even those that are +contained in a \code{PS_META_LIST}. The value \code{iterator} +specifies the iterator to be used. In setting the iterator, the +position of the iterator is defined by \code{location}, which follows +the conventions of the \code{psList} iterators. +\begin{verbatim} +psListIterator *psMetadataIteratorAlloc(psMetadata *md, int location, bool mutable); +bool psMetadataIteratorSet(psListIterator *iterator, int location); +psMetadataItem *psMetadataGetAndIncrement(psListIterator *iterator, const char *regex); +psMetadataItem *psMetadataGetAndDecrement(psListIterator *iterator, const char *regex); +\end{verbatim} + +Metadata items may be printed to an open file descriptor based on a +provided format. The format string is an sprintf format statement +with exactly one \% formatting command. If the metadata item type is +a numeric type, this formatting command must also be numeric, and type +conversion performed to the value to match the format type. If the +metadata item type is a string, the formatting command must also be +for a string (\%s type of command). If the metadata type is any other +data type, printing is not allowed. +\begin{verbatim} +bool psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *md); +\end{verbatim} + +\subsubsection{Configuration files} +\label{sec:configspec} + +It will be necessary for the \PS{} system, in order to load +pre-defined settings, to parse a configuration file into a +\code{psMetadata} structure. This shall be performed by the +function \code{psMetadataParseConfig}, as described below. + +\begin{verbatim} +psMetadata *psMetadataParseConfig(psMetadata *md, int *nFail, const char *filename, bool overwrite); +\end{verbatim} + +Given a metadata container, \code{md}, and the name of a configuration +file, \code{filename}, \code{psMetadataParseConfig} shall parse the +configuration file, placing the contained key/type/value/comment quads +into the metadata, and returning a pointer to the metadata structure. +The number of lines that failed to parse is returned in \code{nFail}. +Multiple specifications of a key that haven't been declared (see +below) are overwritten if and only if \code{overwrite} is \code{true}. +If the metadata container is \code{NULL}, it shall be allocated. + +On error, the function shall return \code{NULL}. + +The configuration file shall consist of plain text with +key/type/value/comment quads on separate lines. Blank lines, +including those consisting solely of whitespace (both spaces and +tabs), shall be ignored, as shall lines that commence with the comment +character (a hash mark, \code{#}), either immediately at the start of +the line, or preceded by whitespace. The key/type/value/comment quads +shall all lie on a single line, separated by whitespace. + +The key shall be first, possibly preceded on the line by whitespace +which should not form part of the key. + +Next, to assist the casting of the value, shall be a string +identifying the type of the value, which shall correspond to one of +the simple types supported in \code{psMetadata}: +\code{STRING,BOOL,S32,F32,F64}; \code{STR} may be used to abbreviate +\code{STRING}. + +\tbd{May, in the future, require more types, including U8,S16,C64, +which will also necessitate updating the definition of psMetadata.} + +The value shall follow the type: strings may consist of multiple +words, and shall have all leading and trailing whitespace removed; +booleans shall simply be either \code{T} or \code{F}. + +Following the value may be an optional comment, preceded by a comment +character (a hash mark, \code{#}), which in the case of a string +value, serves to mark the end of the value, and for other types serves +to identify the comment to the reader. Only one comment character may +be present on any single line (i.e., neither strings nor comments are +permitted to contain the comment character). The comment may consist +of multiple words, and shall have leading and trailing whitespace +removed. + +One wrinkle is the specification of vectors. Keys for which the value +is to be parsed as a vector shall be preceded immediately by a +``vector symbol'', which we choose to be the ``at'' sign, \code{@}. +In this case, the type shall be interpreted as the type for the +vector, which may be any of the signed or unsigned integer or floating +point types (\code{U8,U16,U32,U64,S8,S16,S32,S32,S64,F32,F64}) but not +the complex floating point types; and the value shall consist of +multiple numbers, separated either by a comma or whitespace. These +values shall populate a \code{psVector} of the appropriate type in the +order in which they appear in the configuration file. + +\tbd{May add complex types, likely to be specified with values such as + 1.23+4.56i in the future.} + +An additional hurdle is the specification of keys that may be +non-unique (such as the \code{COMMENT} keyword in a FITS header). +These keys shall be specified in the configuration file as non-unique +by specifying the key at the start of the line (possibly preceded by +whitespace) and specifying the type as a ``multiple symbol'', which we +choose to be an asterisk, \code{*}. No other data may be provided on +this line, though a comment, preceeded by the comment marker, is +valid. A warning shall be produced when a key which has not been +specified to be non-unique is repeated; in this case, the former value +shall be overwritten if \code{overwrite} is \code{true}, otherwise the +line shall be ignored and counted as one that could not be parsed. + +If a line does not conform to the rules laid out here, a warning shall +be generated, it shall be ignored and counted as a line that could not +be parsed. The total number of lines that were not able to be parsed +(including those that were ignored because \code{overwrite} is +\code{false}, and any other parsing problems, but not including blank +lines and comment lines) shall be returned by the function in the +argument \code{nFail}. + +Here are some examples of lines of a valid configuration file: +\filbreak +\begin{verbatim} +Double F64 1.23456789 # This is a comment +Float F32 0.98765 # This is a comment too +String STR This is the string that forms the value #comment + + # This is a comment line and is to be ignored +boolean BOOL T # The value of `boolean' is `true' + +@primes U8 2,3 5 7,11,13 17 # These are prime numbers + +comment MULTI # The rest of this line is ignored, but `comment' is set to be non-unique +comment STR This +comment STR is +comment STR a +comment STR non-unique +comment STR key +Float F64 1.23456 # This generates a warning, and, if `overwrite' is `false', is ignored +\end{verbatim} + +Of course, a real configuration file should look much nicer to humans +than the above example, but PSLib must be able to parse such ugly +files. + +We extend \code{psMetadataParseConfig} to allow a modest tree +structure by defining a reserved keyword \code{TYPE}. Any line in the +config file which starts with the word \code{TYPE} shall be +interpretted as defining a new valid type. The defined type name +follows the word \code{TYPE}, and is in turn followed by an arbitrary +number of words. These words are to be interpreted as the names of an +embedded \code{psMetadata} entry, where the values are given on any +line which (following the \code{TYPE} definition) employs the new type +name. For example, a new type may be defined as: +\begin{verbatim} +TYPE CELL EXTNAME BIASSEC CHIP +CELL.00 CELL CCD00 BSEC-00 CHIP.00 +CELL.01 CELL CCD01 BSEC-01 CHIP.00 +\end{verbatim} + +When \code{psMetadataParseConfig} encounters the \code{TYPE} line, it +should construct a \code{psMetadata} container and fill it with +\code{psMetadataItems} having the names \code{EXTNAME, BIASSEC, CHIP}, +with type \code{PS_META_STR}, but data allocated. When it next +encounters an entry of type \code{CELL}, it should then use the given +name (e.g., \code{CELL.00}) for the \code{psMetadataItem}, and copy +the \code{psMetadata} data onto the \code{psMetadataItem.data.md} +entry, filling in the values from the rest of the line (\code{CCD00, +BSEC-00, CHIP.00}). This hierarchical structure is illustrated in +Figure~\ref{fig:metadata}. + +We further extend \code{psMetadataParseConfig} to allow the definition +of a \code{psMetadata} entry using a sequence of successive lines to +define the values of the \code{psMetadataItem} entries. The initial +line defines the new \code{psMetadata} entry and its name. The +following lines have the same format as the other metadata config file +entries. The sequence is terminated with a line with a single word +\code{END}. For example, a metadata entry may be defined as: +\begin{verbatim} +CELL METADATA + EXTNAME STR CCD00 + BIASSEC STR BSEC-00 + CHIP STR CHIP.00 + NCELL S32 24 +END +\end{verbatim} + +A series of test inputs is contained in +\S\ref{sec:configtest}. + +\subsection{XML Functions} + +Within Pan-STARRS, we will use XML documents as a transport mechanism +to carry data between programs and between IPP and other subsystems. +Configuration information may be stored as well as XML documents, in +addition to the text format discussed in the discussion on Metadata. +XML is an extremely variable document format, and it is not currently +the intention of PSLib to provide a complete PSLib version of XML +operations. Rather, a limited number of operations are defined to +convert specific data structures to an appropriate XML document. To +maximize the simplicity of the XML APIs, we will use the convention +that a single XML document to be parsed by PSLib shall contain only a +single data structure. Each of the XML APIs therefore take as input a +reference to a complete XML document and return a PSLib data +structure, or take a PSLib data structure and return a complete XML +document. + +We start by defining a PSLib wrapper type which is a pointer to an XML +document in memory. We wrap the \code{libxml2} version of an XML +document pointer for now: +\begin{verbatim} +typedef xmlDocPtr psXMLDoc; +void psXMLDocFree(psXMLDoc *doc); +\end{verbatim} + +The next pair of functions convert a \code{psMetadata} data structure +to a complete \code{psXMLDoc} (in memory) and vice versa: +\begin{verbatim} +psXMLDoc *psMetadataToXMLDoc(const psMetadata *metadata); +psMetadata *psXMLDocToMetadata(const psXMLDoc *doc); +\end{verbatim} + +The next pair of functions loads the data in a named file into a +complete \code{psXMLDoc} (in memory) and write out the \code{psXMLDoc} +to a named file: +\begin{verbatim} +psXMLDoc *psXMLParseFile(const char *filename); +int psXMLDocToFile(const psXMLDoc *doc, const char *filename); +\end{verbatim} + +The next pair of functions accepts a block of memory and parses it +into a complete \code{psXMLDoc} (also in memory), and vice versa: +\begin{verbatim} +psXMLDoc *psXMLParseMemory(const char *buffer, const int size); +int psXMLDocToMemory(const psXMLDoc *doc, char *buffer); +\end{verbatim} + +The next pair of functions read from and write to a file descriptor. +The first converts the imcoming data to a complete \code{psXMLDoc} +(also in memory), the second writes the \code{psXMLDoc} to the file +descriptor: +\begin{verbatim} +psXMLDoc *psXMLParseFD(int fd); +int psXMLDocToFD(const psXMLDoc *doc, int fd); +\end{verbatim} + +\subsection{Database Functions} + +Many of the applications that PSLib will be used for will require +access to a simple relational database. PSLib includes generic +database-independent interface mechanisms as part of its API set. The +most important aspect of PSLib's database support is to abstract as +much database specific complexity as is feasible. As almost all RDBMS +provide at least a simple transactional model, commit and rollback +support should be provided. + +Currently, only support for MySQL 4.1.x is required but other backends +may be added as options in the future. As a particular example which +has implications for the database interaction model, support for +SQLite may be required in the future. Currently, the choice of +backend database interface may be made as a compile option. Details +of the specified APIs in the discussion below refer to the relevant +MySQL functions. + +Database errors must be trapped and placed onto the psError stack. +The complete error message should be retrieved with the database's +error function. + +\subsubsection{Managing the Database Connection} + +We specify a database handle which carries the information about the +database connection: + +\begin{verbatim} + typedef struct { + MYSQL *mysql; + } psDB; +\end{verbatim} + +The following collection of functions provides basic database functionality: + +\begin{verbatim} + // wraps mysql_init() & mysql_real_connect() + psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname); + + // wraps mysql_close() + void psDBCleanup(psDB *dbh); + + // wraps mysql_create_db() + bool psDBCreate(psDB *dbh, const char *dbname); + + // wraps mysql_select_db() + bool psDBChange(psDB *dbh, const char *dbname); + + // wraps mysql_drop_db() + bool psDBDrop(psDB *dbh, const char *dbname); +\end{verbatim} + +For MySQL support, \code{psDBInit()} wraps \code{mysql_init()} and +\code{mysql_real_connect()} in order to initialize a psDB structure and +establish a database connection. A null pointer should be returned on +failure. + +When implementing support for SQLite, or other DB which is purely +file-based, the \code{host}, \code{user}, and \code{passwd} arguments +would be ignored while \code{dbname} would specify the path to the +SQLite db file. + +\subsubsection{Interacting with Database Tables} + +The functions in this section perform high level interactions with the +database tables. All of them should behave ``atomically'' with +respect to the state of the database. Specifically, all interactions +with the database should be done as a part of a transaction that is +rolled-back on failure and committed only after all queries used by +the API have been run. In general, this API set attempts to treat a +database table as a 2D matrix where columns can be represented by a +\code{psVector} and rows as a \code{psMetadata} type. A +\code{psMetadata} collection is also used to define the columns of a +table and as part of the query restrictions. + +\begin{verbatim} + bool psDBCreateTable(psDB *dbh, const char *tableName, psMetadata *md); +\end{verbatim} + +This function generates and executes the SQL needed to create a table +named \code{tableName}, with the column names and datatypes as +described in \code{md}. Each data item in the \code{psMetadata} +collection represents a single table field. The name of the field is +given by the name of the \code{psMetadataItem} and the data type is +give by the \code{psMetadataItem.type} and \code{psMetadataItem.ptype} +entries. A lookup table should be used to convert from PSLib types +into MySQL compatible SQL data types. For example, a +\code{PS_META_STR} would map to an SQL99 varchar. If the value of +\code{type} is \code{PS_META_STR} then the \code{psMetadataItem.data} +element is set to a string with the length for the field written as a +text string. The value of the \code{psMetadataItem.data} element is +unused for the \code{PS_META_PRIMITIVE} types. Other metadata types +beyond \code{PS_META_STR} and \code{PS_META_PRIMITIVE} are not allowed +in a table definition metadata collection. + +Database indexes can be specified setting the \code{comment} field to +``\code{Primary Key}'' or ``\code{Key}''. Comment are otherwise +ignored. + +\begin{verbatim} + bool psDBDropTable(psDB *dbh, const char *tableName); +\end{verbatim} + +This function deletes the specified table. + +\begin{verbatim} + psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, const psU64 limit); + psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col, psElemType pType, const psU64 limit); +\end{verbatim} + +These functions generates and executes the SQL needed to select an entire +column from a table or up to \code{limit} rows from it. If \code{limit} is 0, +the entire range is returned. The database response is processed and a +\code{psArray} of strings is returned. The Num version of the function returns +the data in a \code{psVector}, data cast to \code{pType}. It returns an error +(NULL) if the requested field is not a numerical type. + +\begin{verbatim} + psArray *psDBSelectRows(psDB *dbh, const char *tableName, psMetadata *where, const psU64 limit); +\end{verbatim} + +This function returns rows from the specified table which match +the restrictions given by \code{where}. The restrictions are +specified as field / value pairs. The \code{psMetadata} collection +where must consist of valid database fields, though the database query +checking functions may be used to validate the fields as part of the +query. If \code{where} is \code{NULL}, then there are no restrictions +on the rows selected. The selected rows are returned as a +\code{psArray} of \code{psMetadata} values, one per row. + +\begin{verbatim} + bool psDBInsertOneRow(psDB *dbh, const char *tableName, psMetadata *row); +\end{verbatim} + +Insert the data from \code{row} into \code{tableName}. It should be noted in +the API reference that if fields are specified in \code{row} that do not exist +in \code{tablename}, the insert will fail. + +\begin{verbatim} + bool psDBInsertRows(psDB *dbh, const char *tableName, psArray *rowSet); +\end{verbatim} + +Similar to \code{psDBInsertOneRow()}, this function inserts many rows at once +and is atomic for the complete set of rows. + +\begin{verbatim} + psArray *psDBDumpRows(psDB *dbh, const char *tableName); +\end{verbatim} + +Fetch all rows as an psArray of psMetadata. + +\begin{verbatim} + psMetadata *psDBDumpCols(psDB *dbh, const char *tableName); +\end{verbatim} + +Fetch all columns, as either a psVector or a psArray depnding on whether or not +the column is numeric, and return them in a psMetadata structure where +psMetadataItem.name contains the column's name. + +\begin{verbatim} + psS64 psDBUpdateRows(psDB *dbh, const char *tableName, psMetadata *where, psMetadata *values); +\end{verbatim} + +Update the columns contained in \code{values} in the row(s) that have a field +with the value indicated by \code{where} (note that this is only allows very +limited use of SQL99's ``where'' semantics). The number of rows modified is +returned. A negative value is return to indicate an error. If there are +multiple psMetadataItems in \code{where} then each item should be considered as +an additional constraint. e.g. ``where foo = x and where bar = y'' + +\begin{verbatim} + psS64 psDBDeleteRows(psDB *dbh, const char *tableName, psMetadata *where); +\end{verbatim} + +Delete the rows that are matched by \code{where} using the same semantics for +\code{where} as in psDBUpdateRow(). A negative value is returned to indicate an +error. + +\subsection{FITS I/O Functions} + +We need a variety of I/O functions between the disk and certain of our +PSLib data structures. We need the ability to access FITS headers, +images and tables (both ASCII and Binary). We define here the FITS +I/O functions, all of which are currently specified as wrappers to +functions within CFITSIO. CFITSIO provides a wide range of utilities +which we do not feel are particularly appropriate as part of a generic +I/O library, such as assumptions about names which change the data +interpretation, etc. We are defining our calls to avoid the hidden +'features'. The CFITSIO functions which are wrapped should in general +be the most basic versions. + +\begin{verbatim} +typedef struct { + fitsfile fd; +} psFits; +\end{verbatim} +We begin by defining a datatype to wrap the CFITSIO \code{fitsfile} +structure. This is necessary to allow repeated access to the data in +a file without multiple open commands (which are expensive). + +\subsubsection{FITS File Manipulations} + +\begin{verbatim} +psFits *psFitsAlloc(const char *filename); +\end{verbatim} + +Opens a FITS file at positions the pointer to the PHU. + +\begin{verbatim} +bool psFitsMoveExtName(psFits *fits, const char *extname); +\end{verbatim} + +Positions the pointer to the beginning of the specified +\code{extname}. If the \code{extname} does not exist, the function +shall fail. + +\begin{verbatim} +bool psFitsMoveExtNum(psFits* fits, int extnum, bool relative); +\end{verbatim} + +Moves the pointer to the beginning of the specified HDU number. If +\code{relative} is TRUE, \code{extnum} represents the number of HDUs +relative to the current HDU. The PHU is entry number 0, while the +extended data segments start at number 1. + +\begin{verbatim} +int psFitsGetExtNum(psFits* fits); +\end{verbatim} + +Returns the current HDU number (i.e., file position). + +\begin{verbatim} +int psFitsGetSize(psFits* fits); +\end{verbatim} + +Returns the number of HDUs in the file. + +\begin{verbatim} +psFitsType psFitsGetExtType(psFits* fits); +\end{verbatim} + +Gets the current HDU's type (table or image). + +\subsubsection{FITS Header I/O Functions} + +\begin{verbatim} +psMetadata *psFitsReadHeader(psMetadata *out, const psFits *fits); +\end{verbatim} +Read header data into a \code{psMetadata} structure. The data is read +from the current HDU pointed at by the \code{psFits *fits} entry. If +\code{out} is \code{NULL}, a new psMetadata is created. + +\begin{verbatim} +psMetadata *psFitsReadHeaderSet (psFits *fits); +\end{verbatim} +Load a complete set of headers from the \code{psFits} file pointer. +This function loads the headers from all extensions into a +\code{psMetadata} collection, each entry of which is a pointer to a +\code{psMetadata} structure containing the header data. The metadata +entry names are the \code{EXTNAME} values for each header (with the +value of \code{PHU} for the primary header unit). At the start of the +operation, the file pointer is rewound to the beginning of the file. +At the end, it is positioned where it started when the function was +called. + +\begin{verbatim} +bool psFitsWriteHeader(psMetadata *output, const psFits *fits); +\end{verbatim} +Write metadata into the header of a FITS image file. The header is +written at the current HDU. + +\subsubsection{FITS Image I/O Functions} + +\begin{verbatim} +psImage *psFitsReadImage(psImage *output, psFits *fits, psRegion region, int z); +\end{verbatim} +Read an image or subimage from the \code{psFits} file pointer. This +function is a wrapper to the CFITSIO library function. The input +parameters allow a full image or a subimage to be read. The region to +be read is specified by \code{region}. A negative value for either of +\code{region.x1} or \code{region.y1} specifies the size of the region +to be read counting down from the end of the array. + +If the native image is a cube, the value of z specifies the requested +slice of the image. This function must call \code{psError} and return +\code{NULL} if any of the specified parameters are out of range for +the data in the image file, or if the image on disk is zero- or +one-dimensional. This function need only read images of the native +FITS image types (\code{psU8}, \code{psS16}, \code{psS32}, +\code{psF32}, \code{psF64}). The user is expected to convert the data +type as needed with \code{psImageCopy}. + +\begin{verbatim} +bool psFitsUpdateImage(psFits *fits, const psImage *input, psRegion region, int z); +\end{verbatim} +\tbd{we have discussed this as the alternate name} +Write an image section to the open \code{psFits} file pointer. This +operation may write a portion of an image over the existing bytes of +an existing image. Care must be taken to interpret \code{region}, +which specified the output pixels to be written / over-written. If +the combination of \code{region} and the size of \code{psImage *input} +implies writing pixels outside the existing data area of the image, +the function shall return an error (ie, if \code{region.x0 + image.nx +>= NAXIS1}, \code{region.y0 + image.ny >= NAXIS2}, or \code{z >= +NAXIS3}). This function will only write images of the native FITS +image types (\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, +\code{psF64}). The user is expected to convert the data type as +needed with \code{psImageCopy}. The return value must be 0 for a +successful operation and 1 for an error. + +\begin{verbatim} +bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input, int depth); +\end{verbatim} +Create a new image based on the dimensions specified for the image and +the requested depth. The header and image data segments are written +in the file at the current position of the \code{psFits} pointer. +This function will only write images of the native FITS image types +(\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, \code{psF64}). +The user is expected to convert the data type as needed with +\code{psImageCopy}. The return value must be 0 for a successful +operation and 1 for an error. + +\subsubsection{FITS Table I/O Functions} + +\begin{verbatim} +psMetadata *psFitsReadTableRow (psFits *fits, int row); +\end{verbatim} +This function reads a single row of the table in the extension pointed +at by the \code{psFits} file pointer. The row number to be read is +given by \code{row}. The result is returned as a \code{psMetadata} +collection with elements of the apporpriate types and keys +corresponding to the table column names. The function must apply the +needed byte-swapping on the data in the row based on the description +of the table data in the table header. \tbr{we may need to be more +flexible here: if we call this function repeatedly, it would be more +efficient to pass the corresponding header or keep it somewhere (and +the file pointer location, for that matter).} + +\begin{verbatim} +void *psFitsReadTableRowRaw (int *nBytes, psFits *fits, int row); +\end{verbatim} +This function reads a single row of the table in the extension pointed +at by the \code{psFits} file pointer. The row number to be read is +given by \code{row}. The result is returned as collection of +\code{nBytes} bytes allocated by the function. The function must +apply the needed byte-swapping on the data in the row based on the +description of the table data in the table header. \tbr{we may need +to be more flexible here: if we call this function repeatedly, it +would be more efficient to pass the corresponding header or keep it +somewhere (and the file pointer location, for that matter).} + +\begin{verbatim} +psArray *psFitsReadTableColumn (psFits *fits, char *colname); +\end{verbatim} +This function reads a single column of the table in the extension +pointed at by the \code{psFits} file pointer. The column is specified +by the FITS table column key given by \code{row}. The result is +returned as a \code{psArray}, with the data from one row of the table +column per array element. + +\begin{verbatim} +psVector *psFitsReadTableColumnNum (psFits *fits, char *colname); +\end{verbatim} +This function reads a single column of the table in the extension +pointed at by the \code{psFits} file pointer. The column is specified +by the FITS table column key given by \code{row} and must be of a +numeric data type. The result is returned as a \code{psVector} of the +appropriate data type, with the data from one row of the table column +per array element. + +\begin{verbatim} +psArray *psFitsReadTableRaw (int *nBytes, psFits *fits); +\end{verbatim} +This function reads the entire data block from a table into the a +\code{psArray}, with one element of the array per row. The number of +bytes per row is returned in \code{nBytes}. The function must apply +the needed byte-swapping on the data in each row based on the +description of the table data in the table header. + +\begin{verbatim} +psArray *psFitsReadTable (psFits *fits); +\end{verbatim} +This function reads the entire data block from a table into the a +\code{psArray}, with one element of the array per row. Each row is +stored as a \code{psMetadata} collection as described above for +\code{psFitsReadTableRow}. + +\begin{verbatim} +bool psFitsWriteTable(psFits* fits, psMetadata *header, psArray* table); +\end{verbatim} +Accepts a \code{psArray} of \code{psMetadata} and writes it to the +current HDU. If the current HDU is not a table type, this will fail +and return FALSE. + +\begin{verbatim} +bool psFitsUpdateTable(psFits* fits, psMetadata *header, psMetadata* data, int row); +\end{verbatim} +Writes the \code{psMetadata} data to a FITS table at the specified row +in the current HDU. If the current HDU is not a table type, this will +fail and return FALSE. + +\pagebreak \section{Astronomy-Specific Functions} @@ -3457,7 +4540,6 @@ \begin{itemize} \item Dates and times -\item Metadata \item Detector and sky positions -\item Astronomy Image +\item Astronomical Image Containers \item Astrometry \item Photometry @@ -3487,11 +4569,11 @@ typedef enum { - PS_IERS_A, ///< IERS Bulliten A - PS_IERS_B, ///< IERS Bulliten B -} psTimeBulliten; + PS_IERS_A, ///< IERS Bulletin A + PS_IERS_B, ///< IERS Bulletin B +} psTimeBulletin; typedef struct { psS64 sec; ///< seconds, negative values represent dates before 1970 - psU32 nsec; ///< nanseconds + psU32 nsec; ///< nanoseconds bool leapsecond; ///< if time falls on a UTC leapsecond psTimeType type; ///< type of time @@ -3523,11 +4605,12 @@ This function may be used to convert between the various \code{psTimeType} time representations. The \code{time} is modified and returned. Conversion between -all of the \emph{SI} length second systems should be implimented as first +all of the \emph{SI} length second systems should be implimented by first converting to TAI and then to the destination system. UT1 is a special case for conversion as it uses variable length secounds. Conversation to UT1, via TAI, is allowed but conversion \emph{from} UT1 is currently forbidden. -To convert to Local Mean Sidereal Time, it is necessary to provide the local -longitude (specified in radians, positive East of Greenwich) as well: +The following function converts to Local Mean Sidereal Time. It is +necessary to provide the local longitude (specified in radians, +positive East of Greenwich) as well: \begin{verbatim} @@ -3535,19 +4618,27 @@ \end{verbatim} -The function may accept any of the \code{psTimeType} types with \emph{SI} -length seconds. The \code{time} is modified and returned. Note that this -function must supply the value $UT1-UTC$, which is available externally (see -\code{psTimeGetUT1Delta()}) and should use the UT1 values interpolated from -IERS bulliten B. The following utility function encapsulates the PSLib -mechanism to extract the value of $UT1-UTC$: - -\begin{verbatim} -double psTimeGetUT1Delta(const psTime *time, psTimeBulliten bulliten); -psSphere *psTimeGetPoleCoords(const psTime *time, psTimeBulliten bulliten); -\end{verbatim} - -Leap seconds are added to UTC in order to keep it within $0.9s$ of UT1 (which -is defined relative to the Earth's rotation, and hence is useful for -astronomical purposes). +The function may accept any of the \code{psTimeType} types with +\emph{SI} length seconds. The \code{time} is modified and returned. +Note that this function requires the value $UT1-UTC$ (see +\code{psTimeGetUT1Delta()}) and should use the UT1 values interpolated +from IERS bulletin B. + +The following utility function encapsulates the PSLib mechanism to +extract the value of $UT1-UTC$ from the IERS Time Tables: + +\begin{verbatim} +double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin); +\end{verbatim} + +The following function provides tidal corrections to UT1-UTC, using +the Ray model of Simon et al (REF). +\begin{verbatim} +psTime *psTime_TideUT1Corr(const psTime *time); +\end{verbatim} + +Leap seconds are added to UTC in order to keep it within $0.9s$ of UT1 +(which is defined relative to the Earth's rotation, and hence is +useful for astronomical purposes). The following function calculates +the absolute number of leap seconds different between two times. \begin{verbatim} @@ -3555,13 +4646,11 @@ \end{verbatim} -This function calculates the absolute number of leap seconds different between -two times. +The following function accepts \code{PS_TIME_UTC} objects and +determines if the time is potentaly a leapsecond, returning +\code{TRUE} if so. \begin{verbatim} bool psTimeIsLeapSecond(const psTime *utc); \end{verbatim} - -This function only accepts \code{PS_TIME_UTC} objects and determines if the -time is potentaly a leapsecond. \subsubsection{External Date and Time Formats} @@ -3581,7 +4670,7 @@ \end{verbatim} -The \code{psTimeToISO()} function will convertion \code{PS_TIME_UTC} objects -with the \code{leapsecond} flag set to represent the number of seconds as -``:60''. +The \code{psTimeToISO()} function converts \code{PS_TIME_UTC} objects +with the \code{leapsecond} flag set to represent the number of seconds +as ``:60''. \subsubsection{Date and Time Parsing} @@ -3622,17 +4711,17 @@ \code{psS64} (the type of \code{psTime.sec}). -Note that in both these functions, when handling UTC, that the difference -between two times is not inclusive of leap seconds. For example, if we add 30 -seconds to 1998-12-31T23:59:45Z, the result is 1999-01-01T00:00:14Z, since a -leap second was introduced at 1999-01-01T00:00:00Z. - -Time math may only be done on the of \code{psTime} objects of the same type, -and in the case of UT1, the functions shall internally convert the -\code{psTime} inputs to TAI before performing the math; this is in order that -leap seconds are accounted for. - -The type of the time returned by \code{psTimeMath} shall be the same as that of -the input \code{time}. - +Note that in both these functions the difference between two times is +the elapsed number of seconds, inclusing leap seconds. For example, +if we add 30 seconds to 1998-12-31T23:59:45Z, the result is +1999-01-01T00:00:14Z, since a leap second was introduced at +1999-01-01T00:00:00Z. + +Time math may only be done on \code{psTime} objects of the same type, +and except in the case of UT1, the functions shall internally convert +the \code{psTime} inputs to TAI before performing the math; this +ensures that leap seconds are correctly counted. + +The type of the time returned by \code{psTimeMath} shall be the same +as that of the input \code{time}. \subsubsection{Time Tables} @@ -3645,10 +4734,10 @@ \begin{itemize} -\item IERS Bulliten A \& B (1 year ago $\rightarrow$ now + $\sim$3 months) +\item IERS Bulletin A \& B (1 year ago $\rightarrow$ now + $\sim$3 months) \begin{itemize} \item \code{ftp://maia.usno.navy.mil/ser7/finals.daily} \end{itemize} -\item IERS Bulliten A \& B (1973-1-2 $\rightarrow$ now + $\sim$1 year) +\item IERS Bulletin A \& B (1973-1-2 $\rightarrow$ now + $\sim$1 year) \begin{itemize} \item \code{ftp://maia.usno.navy.mil/ser7/finals.all} @@ -3673,5 +4762,5 @@ format of these files shall be simple, for speed in reading. Each line shall contain the date in MJD and the following values from both the A \& B IERS -bullitens: $x_p$ (in arcseconds), $y_p$ (in arcseconds) and $\Delta$UT (in +bulletins: $x_p$ (in arcseconds), $y_p$ (in arcseconds) and $\Delta$UT (in seconds). This format must be readable by \code{psLookupTableRead}. For example: @@ -3756,1027 +4845,15 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Regions} - -In many places, we need to refer to a rectangular area. We define a -structure to represent a rectangle: -\begin{verbatim} -typedef struct { - float x0; - float x1; - float y0; - float y1; -} psRegion; -psRegion *psRegionAlloc (float x0, float x1, float y0, float y1); -\end{verbatim} - -\begin{verbatim} -psRegion *psRegionFromString (char *region); -\end{verbatim} -This function converts the IRAF description of a region in the form -\code{[x0:x1,y0:y1]}, used for header entries such as \code{BIASSEC}, -into the corresponding \code{psRegion} structure. - -\subsection{Metadata} -\label{sec:metadata} - -\subsubsection{Conceptual Overview} - -Within PSLib, we provide a data structure to carry metadata and -mechanisms to manipulate the metadata. Metadata is a general concept -that requires some discussion. In any data analysis task, the -ensemble of all possible data may be divided into two or three -classes: there is the specific data of interest, there is data which -is related or critical but not the primary data of interest, and there -is all of the other data which may or may not be interesting. For -example, consider a simple 2D image obtained of a galaxy from a CCD -camera on a telescope. If you want to study the galaxy, the specific -data of interest is the collection of pixels. There are a variety of -other pieces of data which are closely related and crucial to -understanding the data in those pixels, such as the dimensions of the -image, the coordinate system, the time of the image, the exposure -time, and so forth. Other data may be known which may be less -critical to understanding the image, but which may be interesting or -desired at a later date. For example, the observer who took the -image, the filter manufacturer, the humidity at the telescope, etc. - -Formally, all of the related data which describe the principal data of -interest are metadata. Note that which piece is the metadata and -which is the data may depend on the context. If you are examining the -pixels in an image, the coordinate and flux of an object may be part -of the metadata. However, if you are analyzing a collection of -objects extracted from an image, you may consider then pixel data -simply part of the metadata associated with the list of objects. - -There are various ways to handle metadata vs data within a programming -environment. In C, it is convenient to use structures to group -associated data together. One possibility is to define the metadata -as part of the associated data structure. For example, the image data -structure would have elements for all possible associated measurement. -This approach is both cumbersome (because of the large number metadata -types), impractical (because the full range of necessary metadata is -difficult to know in advance) and inflexible (because any change in -the collection of metadata requires addition of new structure elements -and recompilation). - -An alternative is to place the metadata in a generic container and use -lookup mechanisms to extract the appropriate metadata when needed. An -example of this is would be a text-based FITS header for an image read -into a flat text buffer. In this implementation, metadata lookup -functions could return the current value of, for example, NAXIS1 (the -number of columns of the image) by scanning through the header buffer. -This method has the benefits of flexibility and simplicity of -programming interface, but it has the disadvantage that all metadata -is accessed though this lookup mechanism. This may make the code less -readable and it may slow down the access. - -PSLib implements an intermediate solution to this problem. We specify -a flexible, generic metadata container and access methods. Data types -which require association with a general collection of metadata should -include an entry of this metadata type. However, a subset of metadata -concepts which are basic and frequently required may be placed in the -coded structure elements. This approach allows the code to refer to -the basic metadata concepts as part of the data structure (ie, -\code{image.nx}), but also allows us to provide access to any -arbitrary metadata which may be generated. As a practical matter, the -choice of which entries are only in the metadata and which are part of -the explicit structure elements is rather subjective. Any data -elements which are frequently used should be put in the structure; -those which are only infrequently needed should be left in the generic -metadata. - -There are some points of caution which must be noted. Any -manipulation of the data should be reflected in the metadata where -appropriate. This is always an issue of concern. For example, -consider an image of dimensions \code{nx, ny}. If a function extracts -a subraster, it must change the values of \code{nx, ny} to match the -new dimensions. What should it do to the corresponding metadata? -Clearly, it should change the corresponding value which defines -\code{nX, nY}. However, it is not quite so simple: there may be other -metadata values which depend on those values. These must also be -changed appropriately. What if the metadata element points to a -copy of the metadata which may be shared by other representations of -the image? These must be treated differently because the change would -invalidate those other references. Care must be taken, therefore, -when writing functions which operate on the data to consider all of -the relevant metadata entries which must also be updated. - -A related issue is the definition of metadata names. Entries in a -structure have the advantage of being hardwired: every instance of -that structure will have the same name for the same entry. This is -not necessarily the case with a more flexible metadata container. The -image exposure time is a notorious example in astronomy. Different -observatories use different header keywords (ie, metadata names) for -the same concept of the exposure time (\code{EXPTIME}, -\code{EXPOSURE}, \code{OPENTIME}, \code{INTTIME}, etc). Any system -which operates on these metadata needs to address the issue of -identifying these names. This issue seems like an argument for -hardwiring metadata in the structure, but in fact it does not present -such a strong case. If the metadata are hardwired, some function will -still have to know how to interpret the various names to populate the -structure. The concept can still be localized with generic metadata -containers by including abstract metadata names within the code which -are tied to the various implementations-specific metadata names. - -\subsubsection{Metadata Representation} - -\begin{figure} -\psfig{file=Metadata,width=6.5in} -\caption{Metadata Structures\label{fig:metadata}} -\end{figure} - -This section addresses the question of how \PS{} metadata should be -represented in memory, not how it should be represented on disk. - -We define an item of metadata with the following structure: -\filbreak -\begin{verbatim} -typedef struct { - int id; ///< unique ID for this item - char *name; ///< Name of item - psMetadataType type; ///< type of this item - psElemType ptype; ///< primitive data type - const union { - psS32 S32; ///< integer data - psF32 F32; ///< floating-point data - psF64 F64; ///< double-precision data - void *V; ///< other type - psList *list; ///< psList entry - psMetadata *md; ///< psMetadata entry - } data; ///< value of metadata - char *comment; ///< optional comment ("", not NULL) -} psMetadataItem; -\end{verbatim} - -The \code{id} is a unique identifier for this item of metadata; -experience shows that such tags are useful. The entry \code{name} -specifies the name of the metadata item. The value of the metadata is -given by the union \code{data}, and may be of type \code{psS32}, -\code{psF32}, \code{psF64}, or an arbitrary rich structure pointed at -by the \code{void} pointer \code{V}. A character string comment -associated with this metadata item may be stored in the element -\code{comment}. The \code{type} entry specifies how to interpret the -type of the data being represented, given by the enumerated type -\code{psMetadataType}: -% -\filbreak -\begin{verbatim} -typedef enum { ///< type of item.data is: - PS_META_PRIMITIVE, ///< primitive type: use item.ptype - PS_META_LIST, ///< psList; use item.data.list (used for non-unique data) - PS_META_META, ///< psMetadata: use item.data.list - PS_META_STR, ///< string (item.data.V) - PS_META_MATH, ///< psScalar, psVector, psImage (item.data.V) - PS_META_JPEG, ///< JPEG (item.data) - PS_META_PNG, ///< PNG (item.data) - PS_META_ASTROM, ///< astrometric coefficients (item.data) - PS_META_UNKNOWN, ///< other (item.data) - PS_META_NTYPE ///< Number of types; must be last -} psMetadataType; -\end{verbatim} -If the data is a PSLib primitive data value, the primitive data type -is given by the value of \code{ptype}. - -A collection of metadata is represented by the \code{psMetadata} structure: -\begin{verbatim} -typedef struct { - psList *list; ///< list of psMetadataItem - psHash *table; ///< hash table of the same metadata -} psMetadata; -\end{verbatim} -The type \code{psMetadata} is a container class for metadata. Note -that there are in fact \emph{two} representations of the metadata -(each \code{psMetadataItem} appears on both). The first -representation employs a doubly-linked list that allows the order of -the metadata to be preserved (e.g., if FITS headers are read in a -particular order, they should be written in the same order). The -second representation employs a hash table which allows fast look-up -given a specific metadata keyword. - -Certain metadata names (such as the FITS keywords \code{COMMENT} and -\code{HISTORY} in a FITS header) may be repeated with different -values. In such a case, the \code{psMetadata.list} structure contains -the entries in their original sequence with duplicate keys. The -\code{psMetadata.hash} entries, which are required to have unique -keys, would have a single entry with the keyword of the repeated key, -with the value of \code{psMetadataType} set to \code{PS_META_LIST}, -and the \code{psMetadataItem.data} element pointing to a \code{psList} -containing the actual entries. If \code{psMetadataItemAlloc} is -called with the type set to \code{PS_META_LIST}, such a repeated key -is created. If the data value passed to \code{psMetadataItemAlloc} -(the quantity in ellipsis) is \code{NULL}, then an empty -\code{psMetadataItem} with the given keyword is created to hold future -entries of that keyword. - -The \code{psMetadataAdd} routine is required to check that all -metadata names are unique unless the type is already qualified as -\code{PS_META_LIST}; in this case the data are added to the -corresponding \code{psMetadataItem.data} list. - -\subsubsection{Metadata APIs} - -The allocator for \code{psMetadataItem} returns a full -\code{psMetadataItem} ready for insertion into the \code{psMetadata}. -The \code{name} entry specifies the name to use for this metadata -item, and may include \code{sprintf}-type formating codes. The -\code{comment} entry is a fixed string which is used for the comment -associated with this metadata item. The metadata data and the -arguments to the \code{name} formatting codes are passed, in that -order (metadata pointer first), to \code{psMetadataItemAlloc} as -arguments following the comment string. The data must be a pointer -for any data types which are stored in the element \code{data.void}, -while other data types are passed as numeric values. The argument -list must be interpreted appropriately by the \code{va_list} operators -in the function. -\begin{verbatim} -psMetadataItem *psMetadataItemAlloc(const char *name, psMetadataType type, const char *comment, ...); -psMetadataItem *psMetadataItemAllocV(const char *name, psMetadataType type, const char *comment, va_list list); -\end{verbatim} - -The constructor for the collection of metadata, \code{psMetadata}, -simply returns an empty metadata container (employing the constructors -for the doubly-linked list and hash table). The destructor needs to -free each of the \code{psMetadataItem}s using \code{psMetadataItemFree}. -\begin{verbatim} -psMetadata *psMetadataAlloc(void); -\end{verbatim} - -Items may be added to the metadata in one of two ways --- firstly, an -item may be added by appending a \code{psMetadataItem} which has -already been created; and secondly by directly providing the data to -be appended. In both cases, the return value defines the success -(\code{true}) or failure of the operation. The second function, -\code{psMetadataAdd} takes a pointer or value which is interpreted by -the function using variadic argument interpretation. The third -version is the \code{va_list} version of the second function. All -three functions take a parameter, \code{location}, which specifies -where in the list to place the element, following the conventions for -the \code{psList}. The entry \code{mode} for \code{psMetadataAddItem} -is a bit mask constructed by OR-ing the allowed option flags (eg, -\code{PS_META_REPLACE}) which specifies minor variations on the -behavior. The \code{format} entry, which specifies both the metadata -type and the optional flags, is constructed by bit-wise OR-ing the -appropriate \code{psMetadataType} and allowed option option flags. -Care should be taken not to leak memory when appending an item for -which the key already exists in the metadata (and is not -\code{PS_META_LIST}). -% -\begin{verbatim} -bool psMetadataAddItem(psMetadata *md, psMetadataItem *item, int location, int mode); -bool psMetadataAdd(psMetadata *md, int location, const char *name, int format, const char *comment, ...); -bool psMetadataAddV(psMetadata *md, int location, const char *name, int format, const char *comment, - va_list list); -\end{verbatim} - -The functions above take option flags which modify the behavior when -metadata items are added to the metadata list. These flags must be -bit-exclusive of those used above for the \code{psMetadataTypes}. The -flags have the following meanings: - -\code{PS_META_DEFAULT}: This is the zero bit mask, to allow the -default behavior for \code{psMetadataAddItem} above. If this is OR-ed -with a \code{psMetadataType}, the result is as if no OR-ing took -place. - -\code{PS_META_REPLACE}: If the given metadata item exists in the -metadata list, and is not of type \code{PS_META_LIST} or -\code{PS_META_META} (ie, not a container type), then this entry is -allowed to replace the existing entry. If this mode bit is not set, a -duplicate (non-container-type) entry is an error. - -\begin{verbatim} -typedef enum { ///< option flags for psMetadata functions - PS_META_DEFAULT, ///< default behavior (0x0000) for use in mode above - PS_META_REPLACE, ///< allow entry to be replaced -} psMetadataFlags; -\end{verbatim} - -An example of code to use these metadata APIs to generate the -structure seen in Figure~\ref{fig:metadata} is given below. - -\begin{verbatim} -md = psMetadataAlloc(); - -psMetadataAdd(md, PS_LIST_TAIL, "SIMPLE", PS_META_BOOL, "basic fits", TRUE); -psMetadataAdd(md, PS_LIST_TAIL, "BLANK", PS_META_S32, "invalid pixel data", -32768); -psMetadataAdd(md, PS_LIST_TAIL, "DATE-OBS", PS_META_STR, "observing date UT", " 2004-6-16"); -psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_LIST, "head of comment block", NULL); -psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "DATA"); -psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "PARAMS"); -psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME", PS_META_F32, "exposure time (sec)", 1.05); -psMetadataAdd(md, PS_LIST_TAIL, "COMMENT", PS_META_STR, "", "FOO"); - -cell = psMetadataAlloc(); -psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME", PS_META_STR, "", "CCD00"); -psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR, "", "BSEC-00"); -psMetadataAdd(cell, PS_LIST_TAIL, "CHIP", PS_META_STR, "", "CHIP.00"); -psMetadataAdd(md, PS_LIST_TAIL, "CELL.00", PS_META_META, "", cell); - -cell = psMetadataAlloc(); -psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME", PS_META_STR, "", "CCD01"); -psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_META_STR, "", "BSEC-01"); -psMetadataAdd(cell, PS_LIST_TAIL, "CHIP", PS_META_STR, "", "CHIP.01"); -psMetadataAdd(md, PS_LIST_TAIL, "CELL.01", PS_META_META, "", cell); -\end{verbatim} - -The following code shows how to use the APIs to replace one of these values: -\begin{verbatim} -psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME", PS_META_F32 | PS_REPLACE, "new exposure time (sec)", 2.05); -\end{verbatim} - -Items may be removed from the metadata by specifying a key or a -location in the list. If the value of \code{name} is \code{NULL}, the -value of \code{location} is used. If the value of \code{name} is not -\code{NULL}, then \code{location} must be set to -\code{PS_LIST_UNKNOWN}. If the key matches a metadata item, the item -is removed from the metadata and \code{true} is returned; otherwise, -\code{false} is returned. If the key is not unique, then \emph{all} -items corresponding to the key are removed, and \code{true} is -returned. -% -\begin{verbatim} -bool psMetadataRemove(psMetadata *md, int location, const char *key); -\end{verbatim} - -Items may be found within the metadata by providing a key. In the -event that the key is non-unique, the first item is returned. -\begin{verbatim} -psMetadataItem *psMetadataLookup(const psMetadata *md, const char *key); -\end{verbatim} - -Several utility functions are provided for simple cases. These -functions perform the effort of casting the data to the appropriate -type. The numerical functions shall return 0.0 if their key is not -found. If the pointer value of \code{status} is not \code{NULL}, it -is set to reflect the success or failure of the lookup. -\begin{verbatim} -void *psMetadataLookupPtr(bool *status, const psMetadata *md, const char *key); -psS32 psMetadataLookupS32(bool *status, const psMetadata *md, const char *key); -psF64 psMetadataLookupF64(bool *status, const psMetadata *md, const char *key); -\end{verbatim} - -Items may be retrieved from the metadata by their entry position. The -value of which specifies the desired entry in the fashion of -\code{psList}. -\begin{verbatim} -psMetadataItem *psMetadataGet(const psMetadata *md, int location); -\end{verbatim} - -The metadata list component may be iterated over by using a -\code{psListIterator} in a fashion equivalent to the usage for -\code{psList}. The iterator may be set to a location in the -\code{psMetadata} list, and the user may get the previous or next item -in the list relative to that location. \code{psMetadataGetNext} has -the ability to match the key using a POSIX regex, e.g., if the user -only wants to iterate through \code{IPP.machines.sky} and doesn't want -to bother with \code{IPP.machines.detector}. The iterator should -iterate over every item in the metadata list, even those that are -contained in a \code{PS_META_LIST}. The value \code{iterator} -specifies the iterator to be used. In setting the iterator, the -position of the iterator is defined by \code{location}, which follows -the conventions of the \code{psList} iterators. -\begin{verbatim} -psListIterator *psMetadataIteratorAlloc(psMetadata *md, int location, bool mutable); -bool psMetadataIteratorSet(psListIterator *iterator, int location); -psMetadataItem *psMetadataGetAndIncrement(psListIterator *iterator, const char *regex); -psMetadataItem *psMetadataGetAndDecrement(psListIterator *iterator, const char *regex); -\end{verbatim} - -Metadata items may be printed to an open file descriptor based on a -provided format. The format string is an sprintf format statement -with exactly one \% formatting command. If the metadata item type is -a numeric type, this formatting command must also be numeric, and type -conversion performed to the value to match the format type. If the -metadata item type is a string, the formatting command must also be -for a string (\%s type of command). If the metadata type is any other -data type, printing is not allowed. -\begin{verbatim} -bool psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *md); -\end{verbatim} - -\subsubsection{Configuration files} -\label{sec:configspec} - -It will be necessary for the \PS{} system, in order to load -pre-defined settings, to parse a configuration file into a -\code{psMetadata} structure. This shall be performed by the -function \code{psMetadataParseConfig}, as described below. - -\begin{verbatim} -psMetadata *psMetadataParseConfig(psMetadata *md, int *nFail, const char *filename, bool overwrite); -\end{verbatim} - -Given a metadata container, \code{md}, and the name of a configuration -file, \code{filename}, \code{psMetadataParseConfig} shall parse the -configuration file, placing the contained key/type/value/comment quads -into the metadata, and returning a pointer to the metadata structure. -The number of lines that failed to parse is returned in \code{nFail}. -Multiple specifications of a key that haven't been declared (see -below) are overwritten if and only if \code{overwrite} is \code{true}. -If the metadata container is \code{NULL}, it shall be allocated. - -On error, the function shall return \code{NULL}. - -The configuration file shall consist of plain text with -key/type/value/comment quads on separate lines. Blank lines, -including those consisting solely of whitespace (both spaces and -tabs), shall be ignored, as shall lines that commence with the comment -character (a hash mark, \code{#}), either immediately at the start of -the line, or preceded by whitespace. The key/type/value/comment quads -shall all lie on a single line, separated by whitespace. - -The key shall be first, possibly preceded on the line by whitespace -which should not form part of the key. - -Next, to assist the casting of the value, shall be a string -identifying the type of the value, which shall correspond to one of -the simple types supported in \code{psMetadata}: -\code{STRING,BOOL,S32,F32,F64}; \code{STR} may be used to abbreviate -\code{STRING}. - -\tbd{May, in the future, require more types, including U8,S16,C64, -which will also necessitate updating the definition of psMetadata.} - -The value shall follow the type: strings may consist of multiple -words, and shall have all leading and trailing whitespace removed; -booleans shall simply be either \code{T} or \code{F}. - -Following the value may be an optional comment, preceded by a comment -character (a hash mark, \code{#}), which in the case of a string -value, serves to mark the end of the value, and for other types serves -to identify the comment to the reader. Only one comment character may -be present on any single line (i.e., neither strings nor comments are -permitted to contain the comment character). The comment may consist -of multiple words, and shall have leading and trailing whitespace -removed. - -One wrinkle is the specification of vectors. Keys for which the value -is to be parsed as a vector shall be preceded immediately by a -``vector symbol'', which we choose to be the ``at'' sign, \code{@}. -In this case, the type shall be interpreted as the type for the -vector, which may be any of the signed or unsigned integer or floating -point types (\code{U8,U16,U32,U64,S8,S16,S32,S32,S64,F32,F64}) but not -the complex floating point types; and the value shall consist of -multiple numbers, separated either by a comma or whitespace. These -values shall populate a \code{psVector} of the appropriate type in the -order in which they appear in the configuration file. - -\tbd{May add complex types, likely to be specified with values such as - 1.23+4.56i in the future.} - -An additional hurdle is the specification of keys that may be -non-unique (such as the \code{COMMENT} keyword in a FITS header). -These keys shall be specified in the configuration file as non-unique -by specifying the key at the start of the line (possibly preceded by -whitespace) and specifying the type as a ``multiple symbol'', which we -choose to be an asterisk, \code{*}. No other data may be provided on -this line, though a comment, preceeded by the comment marker, is -valid. A warning shall be produced when a key which has not been -specified to be non-unique is repeated; in this case, the former value -shall be overwritten if \code{overwrite} is \code{true}, otherwise the -line shall be ignored and counted as one that could not be parsed. - -If a line does not conform to the rules laid out here, a warning shall -be generated, it shall be ignored and counted as a line that could not -be parsed. The total number of lines that were not able to be parsed -(including those that were ignored because \code{overwrite} is -\code{false}, and any other parsing problems, but not including blank -lines and comment lines) shall be returned by the function in the -argument \code{nFail}. - -Here are some examples of lines of a valid configuration file: -\filbreak -\begin{verbatim} -Double F64 1.23456789 # This is a comment -Float F32 0.98765 # This is a comment too -String STR This is the string that forms the value #comment - - # This is a comment line and is to be ignored -boolean BOOL T # The value of `boolean' is `true' - -@primes U8 2,3 5 7,11,13 17 # These are prime numbers - -comment MULTI # The rest of this line is ignored, but `comment' is set to be non-unique -comment STR This -comment STR is -comment STR a -comment STR non-unique -comment STR key -Float F64 1.23456 # This generates a warning, and, if `overwrite' is `false', is ignored -\end{verbatim} - -Of course, a real configuration file should look much nicer to humans -than the above example, but PSLib must be able to parse such ugly -files. - -We extend \code{psMetadataParseConfig} to allow a modest tree -structure by defining a reserved keyword \code{TYPE}. Any line in the -config file which starts with the word \code{TYPE} shall be -interpretted as defining a new valid type. The defined type name -follows the word \code{TYPE}, and is in turn followed by an arbitrary -number of words. These words are to be interpreted as the names of an -embedded \code{psMetadata} entry, where the values are given on any -line which (following the \code{TYPE} definition) employs the new type -name. For example, a new type may be defined as: -\begin{verbatim} -TYPE CELL EXTNAME BIASSEC CHIP -CELL.00 CELL CCD00 BSEC-00 CHIP.00 -CELL.01 CELL CCD01 BSEC-01 CHIP.00 -\end{verbatim} - -When \code{psMetadataParseConfig} encounters the \code{TYPE} line, it -should construct a \code{psMetadata} container and fill it with -\code{psMetadataItems} having the names \code{EXTNAME, BIASSEC, CHIP}, -with type \code{PS_META_STR}, but data allocated. When it next -encounters an entry of type \code{CELL}, it should then use the given -name (e.g., \code{CELL.00}) for the \code{psMetadataItem}, and copy -the \code{psMetadata} data onto the \code{psMetadataItem.data.md} -entry, filling in the values from the rest of the line (\code{CCD00, -BSEC-00, CHIP.00}). This hierarchical structure is illustrated in -Figure~\ref{fig:metadata}. - -We further extend \code{psMetadataParseConfig} to allow the definition -of a \code{psMetadata} entry using a sequence of successive lines to -define the values of the \code{psMetadataItem} entries. The initial -line defines the new \code{psMetadata} entry and its name. The -following lines have the same format as the other metadata config file -entries. The sequence is terminated with a line with a single word -\code{END}. For example, a metadata entry may be defined as: -\begin{verbatim} -CELL METADATA - EXTNAME STR CCD00 - BIASSEC STR BSEC-00 - CHIP STR CHIP.00 - NCELL S32 24 -END -\end{verbatim} - -A series of test inputs is contained in -\S\ref{sec:configtest}. - -\subsection{XML Functions} - -Within Pan-STARRS, we will use XML documents as a transport mechanism -to carry data between programs and between IPP and other subsystems. -Configuration information may be stored as well as XML documents, in -addition to the text format discussed in the discussion on Metadata. -XML is an extremely variable document format, and it is not currently -the intention of PSLib to provide a complete PSLib version of XML -operations. Rather, a limited number of operations are defined to -convert specific data structures to an appropriate XML document. To -maximize the simplicity of the XML APIs, we will use the convention -that a single XML document to be parsed by PSLib shall contain only a -single data structure. Each of the XML APIs therefore take as input a -reference to a complete XML document and return a PSLib data -structure, or take a PSLib data structure and return a complete XML -document. - -We start by defining a PSLib wrapper type which is a pointer to an XML -document in memory. We wrap the \code{libxml2} version of an XML -document pointer for now: -\begin{verbatim} -typedef xmlDocPtr psXMLDoc; -void psXMLDocFree(psXMLDoc *doc); -\end{verbatim} - -The next pair of functions convert a \code{psMetadata} data structure -to a complete \code{psXMLDoc} (in memory) and vice versa: -\begin{verbatim} -psXMLDoc *psMetadataToXMLDoc(const psMetadata *metadata); -psMetadata *psXMLDocToMetadata(const psXMLDoc *doc); -\end{verbatim} - -The next pair of functions loads the data in a named file into a -complete \code{psXMLDoc} (in memory) and write out the \code{psXMLDoc} -to a named file: -\begin{verbatim} -psXMLDoc *psXMLParseFile(const char *filename); -int psXMLDocToFile(const psXMLDoc *doc, const char *filename); -\end{verbatim} - -The next pair of functions accepts a block of memory and parses it -into a complete \code{psXMLDoc} (also in memory), and vice versa: -\begin{verbatim} -psXMLDoc *psXMLParseMemory(const char *buffer, const int size); -int psXMLDocToMemory(const psXMLDoc *doc, char *buffer); -\end{verbatim} - -The next pair of functions read from and write to a file descriptor. -The first converts the imcoming data to a complete \code{psXMLDoc} -(also in memory), the second writes the \code{psXMLDoc} to the file -descriptor: -\begin{verbatim} -psXMLDoc *psXMLParseFD(int fd); -int psXMLDocToFD(const psXMLDoc *doc, int fd); -\end{verbatim} - -\subsection{Database Functions} - -Many of the applications that PSLib will be used for will require -access to a simple relational database. PSLib includes generic -database-independent interface mechanisms as part of its API set. The -most important aspect of PSLib's database support is to abstract as -much database specific complexity as is feasible. As almost all RDBMS -provide at least a simple transactional model, commit and rollback -support should be provided. - -Currently, only support for MySQL 4.1.x is required but other backends -may be added as options in the future. As a particular example which -has implications for the database interaction model, support for -SQLite may be required in the future. Currently, the choice of -backend database interface may be made as a compile option. Details -of the specified APIs in the discussion below refer to the relevant -MySQL functions. - -Database errors must be trapped and placed onto the psError stack. -The complete error message should be retrieved with the database's -error function. - -\subsubsection{Managing the Database Connection} - -We specify a database handle which carries the information about the -database connection: - -\begin{verbatim} - typedef struct { - MYSQL *mysql; - } psDB; -\end{verbatim} - -The following collection of functions provides basic database functionality: - -\begin{verbatim} - // wraps mysql_init() & mysql_real_connect() - psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname); - - // wraps mysql_close() - void psDBCleanup(psDB *dbh); - - // wraps mysql_create_db() - bool psDBCreate(psDB *dbh, const char *dbname); - - // wraps mysql_select_db() - bool psDBChange(psDB *dbh, const char *dbname); - - // wraps mysql_drop_db() - bool psDBDrop(psDB *dbh, const char *dbname); -\end{verbatim} - -For MySQL support, \code{psDBInit()} wraps \code{mysql_init()} and -\code{mysql_real_connect()} in order to initialize a psDB structure and -establish a database connection. A null pointer should be returned on -failure. - -When implementing support for SQLite, or other DB which is purely -file-based, the \code{host}, \code{user}, and \code{passwd} arguments -would be ignored while \code{dbname} would specify the path to the -SQLite db file. - -\subsubsection{Interacting with Database Tables} - -The functions in this section perform high level interactions with the -database tables. All of them should behave ``atomically'' with -respect to the state of the database. Specifically, all interactions -with the database should be done as a part of a transaction that is -rolled-back on failure and committed only after all queries used by -the API have been run. In general, this API set attempts to treat a -database table as a 2D matrix where columns can be represented by a -\code{psVector} and rows as a \code{psMetadata} type. A -\code{psMetadata} collection is also used to define the columns of a -table and as part of the query restrictions. - -\begin{verbatim} - bool psDBCreateTable(psDB *dbh, const char *tableName, psMetadata *md); -\end{verbatim} - -This function generates and executes the SQL needed to create a table -named \code{tableName}, with the column names and datatypes as -described in \code{md}. Each data item in the \code{psMetadata} -collection represents a single table field. The name of the field is -given by the name of the \code{psMetadataItem} and the data type is -give by the \code{psMetadataItem.type} and \code{psMetadataItem.ptype} -entries. A lookup table should be used to convert from PSLib types -into MySQL compatible SQL data types. For example, a -\code{PS_META_STR} would map to an SQL99 varchar. If the value of -\code{type} is \code{PS_META_STR} then the \code{psMetadataItem.data} -element is set to a string with the length for the field written as a -text string. The value of the \code{psMetadataItem.data} element is -unused for the \code{PS_META_PRIMITIVE} types. Other metadata types -beyond \code{PS_META_STR} and \code{PS_META_PRIMITIVE} are not allowed -in a table definition metadata collection. - -Database indexes can be specified setting the \code{comment} field to -``\code{Primary Key}'' or ``\code{Key}''. Comment are otherwise -ignored. - -\begin{verbatim} - bool psDBDropTable(psDB *dbh, const char *tableName); -\end{verbatim} - -This function deletes the specified table. - -\begin{verbatim} - psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, const psU64 limit); - psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col, psElemType pType, const psU64 limit); -\end{verbatim} - -These functions generates and executes the SQL needed to select an entire -column from a table or up to \code{limit} rows from it. If \code{limit} is 0, -the entire range is returned. The database response is processed and a -\code{psArray} of strings is returned. The Num version of the function returns -the data in a \code{psVector}, data cast to \code{pType}. It returns an error -(NULL) if the requested field is not a numerical type. - -\begin{verbatim} - psArray *psDBSelectRows(psDB *dbh, const char *tableName, psMetadata *where, const psU64 limit); -\end{verbatim} - -This function returns rows from the specified table which match -the restrictions given by \code{where}. The restrictions are -specified as field / value pairs. The \code{psMetadata} collection -where must consist of valid database fields, though the database query -checking functions may be used to validate the fields as part of the -query. If \code{where} is \code{NULL}, then there are no restrictions -on the rows selected. The selected rows are returned as a -\code{psArray} of \code{psMetadata} values, one per row. - -\begin{verbatim} - bool psDBInsertOneRow(psDB *dbh, const char *tableName, psMetadata *row); -\end{verbatim} - -Insert the data from \code{row} into \code{tableName}. It should be noted in -the API reference that if fields are specified in \code{row} that do not exist -in \code{tablename}, the insert will fail. - -\begin{verbatim} - bool psDBInsertRows(psDB *dbh, const char *tableName, psArray *rowSet); -\end{verbatim} - -Similar to \code{psDBInsertOneRow()}, this function inserts many rows at once -and is atomic for the complete set of rows. - -\begin{verbatim} - psArray *psDBDumpRows(psDB *dbh, const char *tableName); -\end{verbatim} - -Fetch all rows as an psArray of psMetadata. - -\begin{verbatim} - psMetadata *psDBDumpCols(psDB *dbh, const char *tableName); -\end{verbatim} - -Fetch all columns, as either a psVector or a psArray depnding on whether or not -the column is numeric, and return them in a psMetadata structure where -psMetadataItem.name contains the column's name. - -\begin{verbatim} - psS64 psDBUpdateRows(psDB *dbh, const char *tableName, psMetadata *where, psMetadata *values); -\end{verbatim} - -Update the columns contained in \code{values} in the row(s) that have a field -with the value indicated by \code{where} (note that this is only allows very -limited use of SQL99's ``where'' semantics). The number of rows modified is -returned. A negative value is return to indicate an error. If there are -multiple psMetadataItems in \code{where} then each item should be considered as -an additional constraint. e.g. ``where foo = x and where bar = y'' - -\begin{verbatim} - psS64 psDBDeleteRows(psDB *dbh, const char *tableName, psMetadata *where); -\end{verbatim} - -Delete the rows that are matched by \code{where} using the same semantics for -\code{where} as in psDBUpdateRow(). A negative value is returned to indicate an -error. - -\subsection{FITS I/O Functions} - -We need a variety of I/O functions between the disk and certain of our -PSLib data structures. We need the ability to access FITS headers, -images and tables (both ASCII and Binary). We define here the FITS -I/O functions, all of which are currently specified as wrappers to -functions within CFITSIO. CFITSIO provides a wide range of utilities -which we do not feel are particularly appropriate as part of a generic -I/O library, such as assumptions about names which change the data -interpretation, etc. We are defining our calls to avoid the hidden -'features'. The CFITSIO functions which are wrapped should in general -be the most basic versions. - -\begin{verbatim} -typedef struct { - fitsfile fd; -} psFits; -\end{verbatim} -We begin by defining a datatype to wrap the CFITSIO \code{fitsfile} -structure. This is necessary to allow repeated access to the data in -a file without multiple open commands (which are expensive). - -\subsubsection{FITS File Manipulations} - -\begin{verbatim} -psFits *psFitsAlloc(const char *filename); -\end{verbatim} - -Opens a FITS file at positions the pointer to the PHU. - -\begin{verbatim} -bool psFitsMoveExtName(psFits *fits, const char *extname); -\end{verbatim} - -Positions the pointer to the beginning of the specified -\code{extname}. If the \code{extname} does not exist, the function -shall fail. - -\begin{verbatim} -bool psFitsMoveExtNum(psFits* fits, int extnum, bool relative); -\end{verbatim} - -Moves the pointer to the beginning of the specified HDU number. If -\code{relative} is TRUE, \code{extnum} represents the number of HDUs -relative to the current HDU. The PHU is entry number 0, while the -extended data segments start at number 1. - -\begin{verbatim} -int psFitsGetExtNum(psFits* fits); -\end{verbatim} - -Returns the current HDU number (i.e., file position). - -\begin{verbatim} -int psFitsGetSize(psFits* fits); -\end{verbatim} - -Returns the number of HDUs in the file. - -\begin{verbatim} -psFitsType psFitsGetExtType(psFits* fits); -\end{verbatim} - -Gets the current HDU's type (table or image). - -\subsubsection{FITS Header I/O Functions} - -\begin{verbatim} -psMetadata *psFitsReadHeader(psMetadata *out, const psFits *fits); -\end{verbatim} -Read header data into a \code{psMetadata} structure. The data is read -from the current HDU pointed at by the \code{psFits *fits} entry. If -\code{out} is \code{NULL}, a new psMetadata is created. - -\begin{verbatim} -psMetadata *psFitsReadHeaderSet (psFits *fits); -\end{verbatim} -Load a complete set of headers from the \code{psFits} file pointer. -This function loads the headers from all extensions into a -\code{psMetadata} collection, each entry of which is a pointer to a -\code{psMetadata} structure containing the header data. The metadata -entry names are the \code{EXTNAME} values for each header (with the -value of \code{PHU} for the primary header unit). At the start of the -operation, the file pointer is rewound to the beginning of the file. -At the end, it is positioned where it started when the function was -called. - -\begin{verbatim} -bool psFitsWriteHeader(psMetadata *output, const psFits *fits); -\end{verbatim} -Write metadata into the header of a FITS image file. The header is -written at the current HDU. - -\subsubsection{FITS Image I/O Functions} - -\begin{verbatim} -psImage *psFitsReadImage(psImage *output, psFits *fits, psRegion region, int z); -\end{verbatim} -Read an image or subimage from the \code{psFits} file pointer. This -function is a wrapper to the CFITSIO library function. The input -parameters allow a full image or a subimage to be read. The region to -be read is specified by \code{region}. A negative value for either of -\code{region.x1} or \code{region.y1} specifies the size of the region -to be read counting down from the end of the array. - -If the native image is a cube, the value of z specifies the requested -slice of the image. This function must call \code{psError} and return -\code{NULL} if any of the specified parameters are out of range for -the data in the image file, or if the image on disk is zero- or -one-dimensional. This function need only read images of the native -FITS image types (\code{psU8}, \code{psS16}, \code{psS32}, -\code{psF32}, \code{psF64}). The user is expected to convert the data -type as needed with \code{psImageCopy}. - -\begin{verbatim} -bool psFitsUpdateImage(psFits *fits, const psImage *input, psRegion region, int z); -\end{verbatim} -\tbd{we have discussed this as the alternate name} -Write an image section to the open \code{psFits} file pointer. This -operation may write a portion of an image over the existing bytes of -an existing image. Care must be taken to interpret \code{region}, -which specified the output pixels to be written / over-written. If -the combination of \code{region} and the size of \code{psImage *input} -implies writing pixels outside the existing data area of the image, -the function shall return an error (ie, if \code{region.x0 + image.nx ->= NAXIS1}, \code{region.y0 + image.ny >= NAXIS2}, or \code{z >= -NAXIS3}). This function will only write images of the native FITS -image types (\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, -\code{psF64}). The user is expected to convert the data type as -needed with \code{psImageCopy}. The return value must be 0 for a -successful operation and 1 for an error. - -\begin{verbatim} -bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input, int depth); -\end{verbatim} -Create a new image based on the dimensions specified for the image and -the requested depth. The header and image data segments are written -in the file at the current position of the \code{psFits} pointer. -This function will only write images of the native FITS image types -(\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, \code{psF64}). -The user is expected to convert the data type as needed with -\code{psImageCopy}. The return value must be 0 for a successful -operation and 1 for an error. - -\subsubsection{FITS Table I/O Functions} - -\begin{verbatim} -psMetadata *psFitsReadTableRow (psFits *fits, int row); -\end{verbatim} -This function reads a single row of the table in the extension pointed -at by the \code{psFits} file pointer. The row number to be read is -given by \code{row}. The result is returned as a \code{psMetadata} -collection with elements of the apporpriate types and keys -corresponding to the table column names. The function must apply the -needed byte-swapping on the data in the row based on the description -of the table data in the table header. \tbr{we may need to be more -flexible here: if we call this function repeatedly, it would be more -efficient to pass the corresponding header or keep it somewhere (and -the file pointer location, for that matter).} - -\begin{verbatim} -void *psFitsReadTableRowRaw (int *nBytes, psFits *fits, int row); -\end{verbatim} -This function reads a single row of the table in the extension pointed -at by the \code{psFits} file pointer. The row number to be read is -given by \code{row}. The result is returned as collection of -\code{nBytes} bytes allocated by the function. The function must -apply the needed byte-swapping on the data in the row based on the -description of the table data in the table header. \tbr{we may need -to be more flexible here: if we call this function repeatedly, it -would be more efficient to pass the corresponding header or keep it -somewhere (and the file pointer location, for that matter).} - -\begin{verbatim} -psArray *psFitsReadTableColumn (psFits *fits, char *colname); -\end{verbatim} -This function reads a single column of the table in the extension -pointed at by the \code{psFits} file pointer. The column is specified -by the FITS table column key given by \code{row}. The result is -returned as a \code{psArray}, with the data from one row of the table -column per array element. - -\begin{verbatim} -psVector *psFitsReadTableColumnNum (psFits *fits, char *colname); -\end{verbatim} -This function reads a single column of the table in the extension -pointed at by the \code{psFits} file pointer. The column is specified -by the FITS table column key given by \code{row} and must be of a -numeric data type. The result is returned as a \code{psVector} of the -appropriate data type, with the data from one row of the table column -per array element. - -\begin{verbatim} -psArray *psFitsReadTableRaw (int *nBytes, psFits *fits); -\end{verbatim} -This function reads the entire data block from a table into the a -\code{psArray}, with one element of the array per row. The number of -bytes per row is returned in \code{nBytes}. The function must apply -the needed byte-swapping on the data in each row based on the -description of the table data in the table header. - -\begin{verbatim} -psArray *psFitsReadTable (psFits *fits); -\end{verbatim} -This function reads the entire data block from a table into the a -\code{psArray}, with one element of the array per row. Each row is -stored as a \code{psMetadata} collection as described above for -\code{psFitsReadTableRow}. - -\begin{verbatim} -bool psFitsWriteTable(psFits* fits, psMetadata *header, psArray* table); -\end{verbatim} -Accepts a \code{psArray} of \code{psMetadata} and writes it to the -current HDU. If the current HDU is not a table type, this will fail -and return FALSE. - -\begin{verbatim} -bool psFitsUpdateTable(psFits* fits, psMetadata *header, psMetadata* data, int row); -\end{verbatim} -Writes the \code{psMetadata} data to a FITS table at the specified row -in the current HDU. If the current HDU is not a table type, this will -fail and return FALSE. - -\subsection{Detector and Sky Coordinates} +\subsection{Linear and Spherical Coordinates} Both detector and sky positions will be used extensively in the IPP. The first are linear coordinates which conform to Euclidean geometry -while the second are angular coordinates for which additional care -must often be taken. We put these into two structures, \code{psPlane} -and \code{psSphere}, respectively. Partitioning these two will enable -error-checking. +while the second are angular coordinates which define a position on +the sphere of the sky. We put these into two structures, +\code{psPlane} and \code{psSphere}, respectively. Partitioning these +two will enable error-checking. An alternative representation for +angular positions is the 3-D unit vector. These are used in +particular as part of spherical rotation calculations. We define +\code{psCube} to represent such an element. % \begin{verbatim} @@ -4794,4 +4871,13 @@ double dErr; ///< Error in Dec } psSphere; + +typedef struct { + double x; ///< cos (DEC) cos (RA) + double y; ///< cos (DEC) sin (RA) + double z; ///< sin (DEC) + double xErr; ///< Error in x + double yErr; ///< Error in y + double zErr; ///< Error in z +} psCube; \end{verbatim} @@ -4818,6 +4904,7 @@ Constructors for these are straight-forward: \begin{verbatim} -psPlane *psPlaneAlloc(void); +psPlane *psPlaneAlloc(void); psSphere *psSphereAlloc(void); +psCube *psCubeAlloc(void); \end{verbatim} Initialization of the structures is not necessary. @@ -4985,224 +5072,4 @@ \code{NULL}, the function shall generate an error and return \code{NULL}. - - -\subsubsection{Celestial Coordinate Conversions} - -We need to be able to convert between ICRS, Galactic and Ecliptic -coordinates, and potentially between arbitrary spherical coordinate -systems. All of these basic spherical transformations represent -rotations of the spherical coordinate reference. We specify a general -transformation function which takes a structure, -\code{psSphereTransform}, defining the transformation between two -spherical coordinate systems (the structure contains the sines and -cosines of the angles involved so as to minimize computation time for -repeated transformations). We also define a function to generate -\code{psSphereTransform}, based on the three angles -describing the location of the pole and the relative equatorial -rotations of the two systems. We also specify special functions to -return the \code{psSphereTransform} for transformations -between standard coordinate systems. - -\begin{verbatim} -typedef struct { - double cosDeltaP; ///< Cosine of target pole latitude in the source system - double sinDeltaP; ///< Sine of target pole latitude in the source system - double alphaP; ///< Longitude of the target system pole in the source system - double phiP; ///< Longitude of the ascending node in the target system -} psSphereTransform; -\end{verbatim} - -The constructor is defined as follows: -\begin{verbatim} -psSphereTransform *psSphereTransformAlloc(double alphaP, double deltaP, double phiP); -\end{verbatim} -where \code{alphaP} and \code{deltaP} define the coordinates in the -input system of the north pole in the output system and \code{phiP} -defines the longitude in the input system of the equatorial -intersection between the two systems (e.g, the first point of Ares). -The constructor must calculate the sines and cosines above. - -Spherical coordinates may be transformed by providing the -transformation and the coordinate in the input system to -\code{psSphereTransform}: -\begin{verbatim} -psSphere *psSphereTransformApply(psSphere *out, - const psSphereTransform *transform, - const psSphere *coord); -\end{verbatim} - -The following functions simply return the appropriate -\code{psSphereTransform} to convert between predefined spherical -coordinate systems (i.e., ICRS, Ecliptic and Galactic). These are -constructors as well as the above \code{psSphereTransformAlloc}. -% -\begin{verbatim} -psSphereTransform *psSphereTransformICRSToEcliptic(psTime *time); -psSphereTransform *psSphereTransformEclipticToICRS(psTime *time); -psSphereTransform *psSphereTransformICRSToGalactic(void); -psSphereTransform *psSphereTransformGalacticToICRS(void); -\end{verbatim} - -We also require the ability to precess coordinates from one equinox to -another. - -\begin{verbatim} -psSphere *psSpherePrecess(psSphere *coords, const psTime *fromTime, const psTime *toTime); -\end{verbatim} - -Given coordinates, \code{coords}, with equinox for \code{fromTime}, -the coordinates are precessed to equinox for \code{toTime}. The -\code{coords} are modified in-place. Equinoxes shall be Julian -equinoxes (as opposed to Bessellian). - -\subsubsection{Projections} - -We require functions to convert between spherical and linear -coordinate systems based on a variety of projections. The required -projections include: -\begin{itemize} -\item TAN -\item SIN -\item AIT -\item PAR -\end{itemize} - -We specify the following structure \code{psProjection} to define the -parameters of the projection: -\begin{verbatim} -typedef struct { - double R, D; ///< coordinates of projection center - double Xs, Ys; ///< plate-scale in X and Y directions - psProjectionType type; ///< projection type -} psProjection; -\end{verbatim} - -The projection type is defined by the following enumerated type \code{psProjectionType}: -\begin{verbatim} -typedef enum { ///< type of val is: - PS_PROJ_TAN, ///< Tangent projection - PS_PROJ_SIN, ///< Sine projection - PS_PROJ_AIT, ///< Aitoff projection - PS_PROJ_PAR, ///< Par projection - PS_PROJ_NTYPE ///< Number of types; must be last -} psProjectionType; -\end{verbatim} - -The constructor is straight-forward: -\begin{verbatim} -psProjection *psProjectionAlloc(double R, double D, double Xs, double Ys, psProjectionType type); -\end{verbatim} - -The following functions will project and deproject (respectively) -spherical coordinates: - -\begin{verbatim} -psPlane *psProject(const psSphere *coord, const psProjection *projection); -psSphere *psDeproject(const psPlane *coord, const psProjection *projection); -\end{verbatim} - -\subsubsection{Offsets} -We require a function to calculate the offset between two positions on -the sky, as well as a function to apply an offset to a position. The -first determines the offset (RA,Dec) on the sky between two positions. -The second applies the given offset to the coordinate. Both an offset -mode and an offset unit may be defined. The mode may be either -\code{PS_SPHERICAL}, in which case the specified offset corresponds to -an offset in angles, or it may be \code{PS_LINEAR}, in which case the -offset corresponds to a linear offset in a local projection. The -offset unit may be in one of \code{PS_ARCSEC}, \code{PS_ARCMIN}, -\code{PS_DEGREE}, and \code{PS_RADIAN}, which specifies the units of -the offset only. - -\begin{verbatim} -psSphere *psSphereGetOffset(const psSphere *position1, - const psSphere *position2, - psSphereOffsetMode mode, - psSphereOffsetUnit unit); - -psSphere *psSphereSetOffset(const psSphere *position, - const psSphere *offset, - psSphereOffsetMode mode, - psSphereOffsetUnit unit); - -typedef enum { - PS_SPHERICAL; ///< Offset on a sphere - PS_LINEAR; ///< Linear offset -} psSphereOffsetMode; - -typedef enum { - PS_ARCSEC; ///< Arcseconds - PS_ARCMIN; ///< Arcminutes - PS_DEGREE; ///< Degrees - PS_RADIAN; ///< Radians -} psSphereOffsetUnit; -\end{verbatim} -Note that these should propagate the errors appropriately. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%% Astronomical Images and Astrometry -\include{psLibSDRS_Astrom} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsection{Pixel Lists} - -Usually an image mask is the best way to carry information about what -pixels mean what. However, in the case where the number of pixels in -which we are interested is limited, it is more efficient to simply -carry a list of pixels. An example of this is in the image -combination code, where we want to perform an operation on a -relatively small fraction of pixels, and it is inefficient to go -through an entire mask image checking each pixel. - -\begin{verbatim} -typedef struct { - psVector *x; // x coordinate - psVector *y; // y coordinate -} psPixels; -\end{verbatim} - -Of course, the size of each of the vectors should match. In the event -that they do not match, any function which detects the problem shall -generate a warning and use the size of the shorter of the vectors as -the size. The order in which the pixels are kept is not considered -important. - -\begin{verbatim} -psImage *psPixelsToMask(psImage *out, const psPixels *pixels, const psRegion *region, unsigned int maskVal); -psPixels *psMaskToPixels(psPixels *out, const psImage *mask, unsigned int maskVal); -\end{verbatim} - -\code{psPixelsToMask} shall return an image of type U8 with the -\code{pixels} lying within the specified \code{region} set to the -\code{maskVal}. The \code{out} image shall be modified if supplied, -or allocated and returned if \code{NULL}. The size of the output -image shall be \code{region->x1 - region->x0} by \code{region->y1 - -region->y0}, with \code{out->x0 = region->x0} and \code{out->y0 = -region->y0}. In the event that either of \code{pixels} or -\code{region} are \code{NULL}, the function shall generate an error -and return \code{NULL}. - -\code{psMaskToPixels} shall return a \code{psPixels} consisting of the -coordinates in the \code{mask} that match the \code{maskVal}. The -\code{out} pixel list shall be modified if supplied, or allocated and -returned if \code{NULL}. In hte event that \code{mask} is -\code{NULL}, the function shall generate an error and return -\code{NULL}. - -\begin{verbatim} -psPixels *psPixelsConcatenate(psPixels *out, const psPixels *pixels); -\end{verbatim} - -\code{psPixelsConcatenate} shall concatenate \code{pixels} onto -\code{out}. In the event that \code{out} is \code{NULL}, a new -\code{psPixels} shall be allocated, and the contents of \code{pixels} -simply copied in. If \code{pixels} is \code{NULL}, the function shall -generate an error and return \code{NULL}. The function shall take -care to ensure that there are no duplicate pixels in \code{out} (since -the order in which the pixels are stored is not important, the values -may be sorted, allowing the use of a faster algorithm than a linear -scan). \begin{verbatim} @@ -5228,5 +5095,446 @@ and return \code{NULL}. -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\subsubsection{Spherical Rotations} + +Spherical rotations represent coordinate transformation in 3-D, as +well as the effects of precession and nutation. We need spherical +rotatations to convert between ICRS, Galactic and Ecliptic +coordinates, and to determine Alt-Az coordinates for sources. All of +these basic spherical transformations represent rotations of the +spherical coordinate reference. We specify a general transformation +function which takes a structure, \code{psSphereRot}, defining the +transformation between two spherical coordinate systems. The +structure contains the elements of a quaternion to represent the +spherical rotational. We define two allocators for +\code{psSphereRot}, one which defines the rotation in terms of the +coordinate of the pole and the rotation about that pole. The other +defines the rotation from the elements of the quaternion. We also +specify functions to manipulate \code{psSphereRot} in several useful +way. + +\begin{verbatim} +typedef struct { + double q0; + double q1; + double q2; + double q3; +} psSphereRot; +\end{verbatim} + +The constructor is defined as follows: +\begin{verbatim} +psSphereRot *psSphereRotAlloc(double alphaP, double deltaP, double phiP); +\end{verbatim} +where \code{alphaP} and \code{deltaP} define the coordinates in the +input system of the axis of rotation (the north pole of the output +system), while \code{phiP} defines the rotation about that pole. This +last angle is also equal to 270\degree - $\phi_a$, where $\phi_a$ is +the longitude in the output system of the ascending node (equatorial +intersection between the two systems, e.g, the first point of Ares). + +The \code{psSphereRot} may also be constructed by supplying the +elements of the quaternion to the following function: +\begin{verbatim} +psSphereRot *psSphereRotQuat(double q0, double q1, double q2, double q3); +\end{verbatim} +This function normalizes the quaternion, so the input elements need +not be normalized. + +Spherical coordinates may be transformed by providing the +transformation and the coordinate in the input system to +\code{psSphereRot}. The output pointer may be optionally supplied, or +if \code{NULL}, is allocated by the function. + +\begin{verbatim} +psSphere *psSphereRotApply(psSphere *out, const psSphereRot *transform, const psSphere *coord); +\end{verbatim} + +The following function combines two rotations, to produce a single +rotation which is the equivallent of applying the first rotation and +then the second. The output rotation may be supplied, or will be +allocated if \code{NULL}. + +\begin{verbatim} +psSphereRot *psSphereRotCombine(psSphereRot *out, psSphereRot *rot1, psSphereRot *rot2) +\end{verbatim} + +The following function changes the given rotation to its inverse: + +\begin{verbatim} +psSphereRot *psSphereRotInvert(psSphereRot *rot) +\end{verbatim} + +The 3-vector representation of the angles (\code{psCube}) is needed to +implement these functions, and is useful in other circumstances as +well. Two utility functions are provided to convert between the +angular and 3-vector representations: +\begin{verbatim} +psSphere *psCubeToSphere(psCube *cube); +psCube *psSphereToCube(psSphere *sphere); +\end{verbatim} + +\subsubsection{Offsets} +We require a function to calculate the offset between two positions on +the sky, as well as a function to apply an offset to a position. The +first determines the offset (RA,Dec) on the sky between two positions. +The second applies the given offset to the coordinate. Both an offset +mode and an offset unit may be defined. The mode may be either +\code{PS_SPHERICAL}, in which case the specified offset corresponds to +an offset in angles, or it may be \code{PS_LINEAR}, in which case the +offset corresponds to a linear offset in a local projection. The +offset unit may be in one of \code{PS_ARCSEC}, \code{PS_ARCMIN}, +\code{PS_DEGREE}, and \code{PS_RADIAN}, which specifies the units of +the offset only. + +\begin{verbatim} +psSphere *psSphereGetOffset(const psSphere *position1, + const psSphere *position2, + psSphereOffsetMode mode, + psSphereOffsetUnit unit); + +psSphere *psSphereSetOffset(const psSphere *position, + const psSphere *offset, + psSphereOffsetMode mode, + psSphereOffsetUnit unit); + +typedef enum { + PS_SPHERICAL; ///< Offset on a sphere + PS_LINEAR; ///< Linear offset +} psSphereOffsetMode; + +typedef enum { + PS_ARCSEC; ///< Arcseconds + PS_ARCMIN; ///< Arcminutes + PS_DEGREE; ///< Degrees + PS_RADIAN; ///< Radians +} psSphereOffsetUnit; +\end{verbatim} +Note that these should propagate the errors appropriately. + +\subsection{Celestial Coordinate Systems} + +The following functions simply return the appropriate +\code{psSphereRot} to convert between predefined spherical +coordinate systems (i.e., ICRS, Ecliptic and Galactic). These are +constructors as well as the above \code{psSphereRotAlloc}. +% +\begin{verbatim} +psSphereRot *psSphereRotICRSToEcliptic(psTime *time); +psSphereRot *psSphereRotEclipticToICRS(psTime *time); +psSphereRot *psSphereRotICRSToGalactic(void); +psSphereRot *psSphereRotGalacticToICRS(void); +\end{verbatim} + +\subsection{Earth Orientation Calculations} + +One of the critical sets of calculations in astronomy is the sequence +of steps needed to convert between the celestial coordinates of an +object and the observed coordinates of the object. This problem is +best divided into two major components: transformation between the +celestial sphere and coordinates relative to the surface of the solid +earth, excluding the effects of the atmosphere, and compensation for +the effects of the atmosphere. In this section, we address the first +of these two transformations: the Earth Orientation Calculations. + +The Earth Orientation Calculations are further subdivided into several +steps, illustrated in Figure~\ref{CoordinateSystems} . Celestial +coordinates are defined in the International Celestial Reference +System (ICRS), which has the solar barycenter as its reference +position and velocity. The next coordinate system is the Geocentric +Celestial Reference System (GCRS), which uses the earth barycenter as +a reference. The transformation between these two includes the +abberation due to the Earth's velocity, the parallax of the object, +which depends on both the Earth's position and the distance to the +object of interest, and the general relativistic correction for the +bending of light as it approaches the Earth. + +The next set of transformations compenstate for the 3-D rotation of +the Earth on various timescales, including the effects of precession, +nutation, and simple solid-body rotation. These calculations can be +performed using different amounts of information for higher levels of +precision. Since the Earth's rotation is constantly affected by +stochastic processes (weather, earthquakes, etc), these conversions +are constantly modified by observations reported by authoratative +sources such as the US Naval Observatory. The target of this +transformation is the International Terrestrial Reference System +(ITRS), which is fixed with respect to the Earth's crust. This +transformation is subdivided into slow precesion and nutation +(yielding the coordinate system CIP/CEO), followed by the Earth's +rotation (yielding the coordinate system CIP/TEO), and finally +corrections for the short-period motion of the Earth's pole. + +\subsubsection{Transformation from ICRS to GCRS} + +\tbd{we need a function to construct the direction and speed elements + given the time}. + +\tbd{supply the velocity as an un-normalized 3 vector?} + +\paragraph{Aberration} +The following function calculates the \code{apparent} position of a +star, given its \code{actual} position and the velocity vector of the +observer, represented as a speed and a direction: +\begin{verbatim} +psAberration(psSphere *apparent, psSphere *actual, psSphere direction, double speed); +\end{verbatim} +The \code{actual} and \code{apparent} positions are represented as +\code{psSphere} entries, as is the \code{direction} of motion. The +speed in that direction is given in units of the speed of light. + +\paragraph{Gravitational Deflection} + +\paragraph{Parallax} + +\begin{verbatim} +double psEOC_ParallaxFactor(psSphere *coords, psTime *time); +\end{verbatim} +Calculate the parallax factor for the given position and time. + +\subsubsection{Transformation from GCRS to ITRS} + +\paragraph{Precession/Nutation} + +The following routine calculates the components of the rotation +between the CEO and GCRS frames, $X$, $Y$, and $s$, using to the +IAU2000A precession \& nutation model: +% +\begin{verbatim} +psSphere *psEOC_PrecessionModel(double *s, const psTime *time) +\end{verbatim} +% +The input to this function is the desired \code{time}, which may be +represented in any format other than UT1. This routine must give +results identical to the IERS XYS2000A subroutine within the limits of +machine accuracy. + +The following function provides interpolated corrections to $X$ and +$Y$ from the tables provided by the IERS, just as it does for UT1 and +polar motion. + +\begin{verbatim} +psSphere *psEOC_GetPolarCorr(const psTime *time, psTimeBulletin bulletin); +\end{verbatim} + +The polar correction is applied to the $X$ and $Y$ elements of the +rotation to provide higher accuracy. The spherical rotation term is +generated by providing the three elements of the rotation to the +following function: +\begin{verbatim} +psSphereRot *psSphereRot_CEOtoGCRS(double s, const psSphere *pole) +\end{verbatim} +The retulting \code{psSphereRot} may be used to determine the rotation +from CIP/CEO to GCRS. This function must give results identical to +the IERS BPN2000, within the limits of machine accuracy. + +\paragraph{Earth Rotation} + +The following routine calculates the rotation of the Earth about the CIP: +\begin{verbatim} +psSphereRot *psSphereRot_TEOtoCEO(const psTime *time) +\end{verbatim} +The IERS code to create the comparable rotation is embedded in +T2C2000, with the Earth Rotation Angle calculated by ERA2000. + +\paragraph{Polar Motion} + +The following function provides interpolated values of the polar +motion components, $x_p$ and $y_p$, extracted from the IERS tables. +\begin{verbatim} +psSphere *psEOC_GetPoleCoords(const psTime *time, psTimeBulletin bulletin); +\end{verbatim} + +The following function provides tidal corrections to the polar motion +components, $x_p$ and $y_p$, using the Ray model of Simon et al (see +ADD). +\begin{verbatim} +psSphere *psEOC_TidePolarCorr(const psTime *time); +\end{verbatim} + +The following function provides the additional corrections due to nutation +terms with periods less than or equal to two days: +\begin{verbatim} +psSphere *psEOC_NutationCorr(psTime *time); +\end{verbatim} + +The following function should generate the \code{psSphereRot} transform from +ITRS to CIP/TEO: +\begin{verbatim} +psSphereRot *psSphereRot_ITRStoTEO(psSphere pole, psTime *time); +\end{verbatim} +The time argument should be used to internally calculate $s'$. +This function should give identical results to the IERS POM2000 subroutine. + +\subsubsection{Earth Orientation Wrappers} + +The following function generates the complete spherical rotation to +account for precession between two times. If \code{NULL} is provided +for either time, it is assumed to have the reference equinox value of +J2000. +\begin{verbatim} +psSphere *psSpherePrecess(const psTime *fromTime, const psTime *toTime, psPrecessMethod mode); +\end{verbatim} +The mode argument is used to specify the level of detail used in the +calculation. + +\begin{verbatim} +typedef enum { + PS_PRECESS_ROUGH, + PS_PRECESS_COMPLETE, + PS_PRECESS_IAU2000A, +} psPrecessMethod; +\end{verbatim} + +\subsection{Atmospheric Effects} + +\tbd{The ATM effects components should be deferred until we clean up + the refraction definitions} + +A-priori astrometric transformations between the telescope orientation +(Alt/Az) and the predicted stellar coordinates above the atmosphere +(DEC/HA) requires several pieces of information describing the current +environmental conditions. These quantities are consistent across an +image, and may vary only slowly with time. Pre-computing these +quantities for exposures means that subsequent transformations are +faster. The structure below carries the environmental data of interest. +For historical reasons, this structure is known colloquially as ``the +Grommit''. + +\tbd{this structure needs to be modified to correspond to what we + actually need to carry around for the atmosphere functions} + +\tbd{provide a single Grommit to carry around all EOC + ATM + pre-calculated entries and a separate structure for ATM effect?} + +\begin{verbatim} +typedef struct { + const double latitude; ///< geodetic latitude (radians) + const double longitude; ///< longitude + ... (radians) + const double height; ///< height (HM) + const double abberationMag; ///< magnitude of diurnal aberration vector + const double temperature; ///< ambient temperature (TDK) + const double pressure; ///< pressure (PMB) + const double humidity; ///< relative humidity (RH) + const double wavelength; ///< wavelength (WL) + const double lapseRate; ///< lapse rate (TLR) + const double refractA, refractB; ///< refraction constants A and B (radians) + const double siderealTime; ///< local apparent sidereal time (radians) +} psGrommit; +\end{verbatim} + +The \code{psGrommit} is calculated from telescope information for the +particular exposure, \code{exp}: +\begin{verbatim} +psGrommit *psGrommitAlloc(const psExposure *exp); +\end{verbatim} + +\tbd{these functions probably need to take the ATM structure} + +We require additional functions to perform general functions which +will be useful for astrometry. Given coordinates on the sky, we +need to get the airmass, the parallactic angle, and an estimate of +the atmospheric refraction. + +\begin{verbatim} +float psGetAirmass(const psSphere *coord, psTime *lst, float height); +\end{verbatim} +which returns the airmass for a given position and local sidereal time +(\code{lst}). + +\begin{verbatim} +float psGetParallactic(const psSphere *coord, double siderealTime); +\end{verbatim} +which returns the parallactic angle for a given position and sidereal time. + +\begin{verbatim} +float psGetRefraction(float colour, ///< Colour of object + psPhotSystem colorPlus, ///< Colour reference + psPhotSystem colorMinus, ///< Colour reference + const psExposure *exp); ///< Telescope pointing information +\end{verbatim} +which provides an estimate of the atmospheric refraction, along the parallactic angle. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Fixed Pattern} + +The fixed pattern is a correction to the general astrometric solution +formed by summing the residuals from many observations. The intent is +to correct for higher-order distortions in the camera system on a +coarse grid (larger than individual pixels, but smaller than a single +cell). Hence, in addition to the offsets, we need to specify the size +and scale of the grid in $x$ and $y$, as well as the origin of the +grid. + +\begin{verbatim} +typedef struct { + int nX, nY; ///< Number of elements in x and y + double x0, y0; ///< Position of 0,0 corner on focal plane + double xScale, yScale; ///< Scale of the grid + double **x, **y; ///< The grid of offsets in x and y +} psFixedPattern; +\end{verbatim} + +The constructor for \code{psFixedPattern} shall be: +\begin{verbatim} +psFixedPattern *psFixedPatternAlloc(double x0, double y0, + double xScale, double yScale, + const psImage *x, const psImage *y); +\end{verbatim} +Here, \code{x0}, \code{y0}, \code{xScale} and \code{yScale} have the +same meaning as in the \code{psFixedPattern} structure. Note that the +values of the fixed pattern offsets are specified as images, the +values from which need to be copied into the \code{double **x} and +\code{double **y} of \code{psFixedPattern}, and that the number of +elements may be derived from the size of the images. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 + +\subsection{Projections} + +We require functions to convert between spherical and linear +coordinate systems based on a variety of projections. The required +projections include: +\begin{itemize} +\item TAN +\item SIN +\item AIT +\item PAR +\end{itemize} + +We specify the following structure \code{psProjection} to define the +parameters of the projection: +\begin{verbatim} +typedef struct { + double R, D; ///< coordinates of projection center + double Xs, Ys; ///< plate-scale in X and Y directions + psProjectionType type; ///< projection type +} psProjection; +\end{verbatim} + +The projection type is defined by the following enumerated type \code{psProjectionType}: +\begin{verbatim} +typedef enum { ///< type of val is: + PS_PROJ_TAN, ///< Tangent projection + PS_PROJ_SIN, ///< Sine projection + PS_PROJ_AIT, ///< Aitoff projection + PS_PROJ_PAR, ///< Par projection + PS_PROJ_NTYPE ///< Number of types; must be last +} psProjectionType; +\end{verbatim} + +The constructor is straight-forward: +\begin{verbatim} +psProjection *psProjectionAlloc(double R, double D, double Xs, double Ys, psProjectionType type); +\end{verbatim} + +The following functions will project and deproject (respectively) +spherical coordinates: + +\begin{verbatim} +psPlane *psProject(const psSphere *coord, const psProjection *projection); +psSphere *psDeproject(const psPlane *coord, const psProjection *projection); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 \subsection{Photometry} @@ -5298,5 +5606,8 @@ M_{\rm pM} - pA, M_{\rm sP} - M_{\rm sM} - sA)$. -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Astronomical Images and Astrometry +\include{psLibSDRS_Astrom} \subsection{Astronomical objects} @@ -5329,4 +5640,5 @@ \appendix +\pagebreak \section{Configuration File Test Inputs} \label{sec:configtest} Index: /trunk/doc/pslib/psLibSDRS_Astrom.tex =================================================================== --- /trunk/doc/pslib/psLibSDRS_Astrom.tex (revision 3563) +++ /trunk/doc/pslib/psLibSDRS_Astrom.tex (revision 3564) @@ -270,5 +270,5 @@ \code{TRUE} if they are all correctly assigned, otherwise \code{FALSE}. -\subsubsection{Coordinate Transformations} +\subsubsection{Detector Coordinate Transformations} \begin{figure} @@ -322,12 +322,271 @@ \end{verbatim} -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% add psExposure to metadata? -% add grommit to metadata? -% add photsystem data to metadata? -% add statistics to metadata? - -\subsubsection{Observatory data} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Astrometry} + +Astrometry is a basic functionality required for the IPP that will be +used repeatedly, both for low-precision (roughly where is my favorite +object?) and high-precision (what is the proper motion of this star?). +As such, it must be flexible, yet robust. + +\subsubsection{Coordinate frames} +\label{sec:coordinateFrames} + +There are five coordinate frames that we need to worry about for the +purposes of astrometry: +\begin{itemize} +\item Cell: $(x,y)$ in pixels --- raw coordinates; +\item Chip: $(X,Y)$ in pixels --- the location on the silicon; +\item Focal Plane: $(p,q)$ in microns --- the location on the focal plane; +\item Tangent Plane: $(l,m)$ in arcsec from the telescope boresight; and +\item Sky: (RA,Dec) --- ICRS. +\end{itemize} + +The following steps are required to convert from the cell coordinates to +the sky: +\begin{itemize} +\item Cell $\longleftrightarrow$ Chip: two 2D polynomials, $(X,Y) = f(x,y)$; +\item Chip $\longleftrightarrow$ FP: two 2D polynomials, $(p,q) = g(X,Y)$; +\item FP $\longleftrightarrow$ TP: two 4D polynomials, $(l,m) = +h(p,q,m,c)$, where $m$ and $c$ are the magnitude and color of the +object, respectively; and +\item TP $\longleftrightarrow$ Sky: transformation to the sky using +pre-computed coefficients for each pointing. +\end{itemize} + +Note that the transformation between the Focal Plane and the Tangent +Plane is a four-dimensional polynomial, in order to account for any +possible dependencies in the astrometry on the stellar magnitude and +color; the former serves as a check for charge transfer +inefficiencies, while the latter will correct chromatic refraction, +both through the atmosphere and the corrector lenses. + +We require structures to contain each of the above transformations as +well as the pixel data. + +\subsubsection{Position Finding} + +We require functions to return the structure containing given +coordinates. For example, we want the chip that corresponds to the +focal plane coordinates $(p,q) = (-1.234,+5.678)$. These routines +handle the one-to-many problem --- i.e., for one given focal plane +coordinate, there are many chips that this coordinate may be +correspond to; these functions will select the correct one. +% +\begin{verbatim} +psCell *psCellInFPA (const psPlane *coord, const psFPA *fpa); +psChip *psChipInFPA (const psPlane *coord, const psFPA *fpa); +psCell *psCellInChip(const psPlane *coord, const psChip *chip); +\end{verbatim} + +\subsubsection{Conversion Functions} + +We require functions to convert between the various coordinate frames +(Section~\ref{sec:coordinateFrames}). The hierarchy of the coordinate +frames and the transformations between each are shown in +Figure~\ref{fig:coco}. The functions that employ the transformations +are shown in Figure~\ref{fig:cocoFunc}. In addition to +transformations between each adjoining coordinate frame in the +hierarchy, we also require higher-level functions to convert between +the Cell and Sky coordinate frames; these will simply perform the +intermediate steps. + +\begin{figure} +\psfig{file=coordinateFrames,height=7in,angle=-90} +\caption{The coordinate systems in the \PS{} IPP, and the relation +between each by transformations contained in the appropriate +structures.} +\label{fig:coco} +\end{figure} + +\begin{figure} +\psfig{file=coordinateConv,height=7in,angle=-90} +\caption{Conversion between coordinate systems by PSLib.} +\label{fig:cocoFunc} +\end{figure} + +We specify the following functions to convert between coordinates in +one type of frame to another type of frame. The first group consist +of unambiguous transformations: from the coordinates in a low-level +frame to the coordinates in the containing higher-level frame, of +which only one exists. In all of these functions, the output +coordinate structure may be \code{NULL} or may be supplied by the +calling function. In the former case, the structure must be +allocated; in the latter case, the supplied structure must be used. + +\begin{verbatim} +psPlane *psCoordCellToChip (psPlane *out, const psPlane *in, const psCell *cell); +% astrometry comes from cell (no need for parent) +\end{verbatim} +which converts coordindates \code{in} on the specified \code{cell} to +the coordinates on the parent chip. + +\begin{verbatim} +psPlane *psCoordChipToFPA (psPlane *out, const psPlane *in, const psChip *chip); +% astrometry comes from chip (no need for parent) +\end{verbatim} +which converts the coordinates \code{in} on the specified \code{chip} +to the coordinates on the parent FPA. + +\begin{verbatim} +psPlane *psCoordFPAToTP(psPlane *out, const psPlane *in, float color, float mag, const psFPA *fpa); +% astrometry comes from FPA (no need for parent) +\end{verbatim} +which converts coordinates \code{in} on the specified focal plane +\code{fpa} to tangent plane coordinates, applying the appropriate +distortion terms. The \code{color} and magnitude (\code{mag}) of the +source is necessary in order to perform the distortion between the +focal plane and the tangent plane. + +\begin{verbatim} +psSphere *psCoordTPToSky(psSphere *out, const psPlane *in, const psGrommit *grommit); +\end{verbatim} +which converts the tangent plane coordinates \code{in} to (RA,Dec) on +the sky, based on the environmental information specified by +\code{grommit}. + +% astrometry comes from cell +\begin{verbatim} +psPlane *psCoordCellToFPA(psPlane *out, const psPlane *in, const psCell *cell); +\end{verbatim} +which performs the single-step conversion between Cell coordinates +\code{in} and FPA coordinates. + +% astrometry comes from cell,chip,fpa (PARENT IS NEEDED HERE) +\begin{verbatim} +psSphere *psCoordCellToSky(psSphere *out, const psPlane *in, float color, float mag, const psCell *cell); +\end{verbatim} +which converts coordinates on the specified cell to (RA,Dec). This +transformation must be performed using the intermediate stage +transformations of Cell to Chip, Chip to FPA, FPA to Tangent Plane, +Tangent Plane to Sky. The information needed for each of these +transformations is available in the \code{.parent} elements of +\code{psCell} and \code{psChip}, and the \code{psFPA.exposure} +element. The \code{color} and magnitude (\code{mag}) of the source is +necessary in order to perform the distortion between the focal plane +and the tangent plane. + +% astrometry comes from cell (no need for parent) +\begin{verbatim} +psSphere *psCoordCellToSkyQuick(psSphere *out, const psPlane *in, const psCell *cell); +\end{verbatim} +which uses the 'quick-and-dirty' transformation to convert coordinates +on the specified cell to (RA,Dec). This transformation should use the +locally linear transformation specified by the element +\code{psCell.toTP}. Although the accuracy of this transformation +is lower than the complete transformation above, the calculation is +substantially faster as it only involves linear transformations. + +The following functions convert from high-level frames to the +coordinates of contained lower-level frames. + +\begin{verbatim} +psPlane *psCoordSkyToTP(psPlane *out, const psSphere *in, const psGrommit *grommit); +\end{verbatim} +which converts (RA,Dec) coordinates \code{in} to tangent plane coords +based on the enviromental information supplied by \code{grommit}. + +\begin{verbatim} +psPlane *psCoordTPToFPA(psPlane *out, const psPlane *in, float color, float mag, const psFPA *fpa); +\end{verbatim} +which converts the tangent plane coordinates \code{in} to focal plane +coordinates. The \code{color} and magnitude (\code{mag}) of the +source is necessary in order to perform the distortion between the +focal plane and the tangent plane. + +\begin{verbatim} +psPlane *psCoordFPAToChip (psPlane *out, const psPlane *in, const psChip *chip); +\end{verbatim} +which converts the specified FPA coordinates \code{in} to the +coordinates on the given Chip. The specified chip need not contain +the input coordinate. To find the chip which contains a particular +coordinate, the function \code{psChipInFPA}, defined above, should be +used. + +\begin{verbatim} +psPlane *psCoordChipToCell (psPlane *out, const psPlane *in, const psCell *cell); +\end{verbatim} +which converts the specified Chip coordinate \code{in} to the +coordinate on the given Cell. The specified Cell need not contain the +input coordinate. To find the cell which contains a particular +coordinate, the function \code{psCellInChip}, defined above, should be +used. + +\begin{verbatim} +psPlane *psCoordSkyToCell(psPlane *out, const psSphere *in, float color, float mag, psCell *cell); +\end{verbatim} +which directly converts (RA,Dec) \code{in} to coordinates on the +specified cell. The specified cell need not contain the input +coordinates. The \code{color} and magnitude (\code{mag}) of the +source is necessary in order to perform the distortion between the +focal plane and the tangent plane. + +\begin{verbatim} +psPlane *psCoordSkyToCellQuick(psPlane *out, const psSphere *in, psCell *cell); +\end{verbatim} +which directly converts (RA,Dec) \code{in} to coordinates on the +specified cell. The specified cell need not contain the input +coordinates. This transformation should use the locally linear +transformation specified by the element \code{psCell.toTP}. +Although the accuracy of this transformation is lower than the +complete transformation above, the calculation is substantially faster +as it only involves linear transformations. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Astrometry and World Coordinate System} + +The FITS World Coordinate System (WCS) headers are commonly employed +with astronomical images in order to relate pixels to celestial (or +otherwise) coordinates. Since it is a FITS standard, we must be able +to read and write from WCS into our internal format. For the time +being, we will consider only celestial WCS (i.e., no spectral +wavelength calibrations, etc). Because WCS does not support the +multiple layers that we have built for \PS{}, we will use a simple +internal representation: a transformation, which handles any +distortions (i.e., goes directly from the coordinate frame of the +image to the tangent plane); and the projection. + +\begin{verbatim} +bool psAstrometryReadWCS(psPlaneTransform **transform, // Output transformation + psProjection **projection, // Output projection + psMetadata *header // Input FITS header + ); +bool psAstrometryWriteWCS(psMetadata *header, // Output FITS header + psPlaneTransform *transform, // Input transformation + psProjection *projection, // Input projection + double color, // Mean color to use + double magnitude, // Mean magnitude to use + ); +bool psAstrometrySimplify(psPlaneTransform **transform, // Output transformation + psProjection **projection, // Output projection + psCell *cell // Cell for which to generate transform and projection + ); +\end{verbatim} + +\code{pmReadAstrometry} shall parse the specified FITS \code{header}, +returning new instances of the \code{transform} and \code{projection} +that represent the WCS. The function shall return \code{true} if it +was able to successfully generate the outputs; otherwise it shall +return \code{false}. + +\code{pmWriteAstrometry} shall add WCS keywords to the supplied FITS +\code{header} that implement the given \code{transform} and +\code{projection}. The function shall return \code{true} if it was +able to successfully generate the output; otherwise it shall return +\code{false}. + +\code{pmSimplifyAstrometry} shall take a \code{cell} and simplify the +internal astrometric representation (\code{cell->toFPA} or equivalent, +\code{cell->parent->parent->toTangentPlane} and +\code{cell->parent->parent->grommit}) to a single \code{transform} and +\code{projection}. This allows the subsequent use of +\code{pmWriteAstrometry} in the case that we have only the +multi-layered \PS{} internal astrometric representation. The function +shall return \code{true} if it was able to successfully generate the +output; otherwise it shall return \code{false}. + +\subsection{Observatory data} We need a container for the observatory data that doesn't change per @@ -350,5 +609,5 @@ \end{verbatim} -\subsubsection{Exposure information} +\subsection{Exposure information} We need several quantities from the telescope in order to make a @@ -388,369 +647,2 @@ \end{verbatim} -\subsubsection{Environmental Information} - -A-priori astrometric transformations between the tangent plane and the -sky require several pieces of information describing the current -environmental conditions. These quantities are consistent from image -to image, and may vary only slowly with time. Pre-computing these -quantities for exposures means that subsequent transformations are -faster. The structure below carries the environment data of interest. -For historical reasons, this structure is known colloquially as ``the -Grommit''. - -\begin{verbatim} -typedef struct { - const double latitude; ///< geodetic latitude (radians) - const double sinLat, cosLat; ///< sine and cosine of geodetic latitude - const double abberationMag; ///< magnitude of diurnal aberration vector - const double height; ///< height (HM) - const double temperature; ///< ambient temperature (TDK) - const double pressure; ///< pressure (PMB) - const double humidity; ///< relative humidity (RH) - const double wavelength; ///< wavelength (WL) - const double lapseRate; ///< lapse rate (TLR) - const double refractA, refractB; ///< refraction constants A and B (radians) - const double longitudeOffset; ///< longitude + ... (radians) - const double siderealTime; ///< local apparent sidereal time (radians) -} psGrommit; -\end{verbatim} - -The \code{psGrommit} is calculated from telescope information for the -particular exposure, \code{exp}: -\begin{verbatim} -psGrommit *psGrommitAlloc(const psExposure *exp); -\end{verbatim} - -\subsubsection{Fixed Pattern} - -The fixed pattern is a correction to the general astrometric solution -formed by summing the residuals from many observations. The intent is -to correct for higher-order distortions in the camera system on a -coarse grid (larger than individual pixels, but smaller than a single -cell). Hence, in addition to the offsets, we need to specify the size -and scale of the grid in $x$ and $y$, as well as the origin of the -grid. - -\begin{verbatim} -typedef struct { - int nX, nY; ///< Number of elements in x and y - double x0, y0; ///< Position of 0,0 corner on focal plane - double xScale, yScale; ///< Scale of the grid - double **x, **y; ///< The grid of offsets in x and y -} psFixedPattern; -\end{verbatim} - -The constructor for \code{psFixedPattern} shall be: -\begin{verbatim} -psFixedPattern *psFixedPatternAlloc(double x0, double y0, - double xScale, double yScale, - const psImage *x, const psImage *y); -\end{verbatim} -Here, \code{x0}, \code{y0}, \code{xScale} and \code{yScale} have the -same meaning as in the \code{psFixedPattern} structure. Note that the -values of the fixed pattern offsets are specified as images, the -values from which need to be copied into the \code{double **x} and -\code{double **y} of \code{psFixedPattern}, and that the number of -elements may be derived from the size of the images. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsection{Astrometry} - -Astrometry is a basic functionality required for the IPP that will be -used repeatedly, both for low-precision (roughly where is my favorite -object?) and high-precision (what is the proper motion of this star?). -As such, it must be flexible, yet robust. Accordingly, we will wrap -the StarLink Astronomy Libraries (SLALib), which has already been -developed. \tbd{SLAlib functions will be replaced in the next -release}. - -\subsubsection{Coordinate frames} -\label{sec:coordinateFrames} - -There are five coordinate frames that we need to worry about for the -purposes of astrometry: -\begin{itemize} -\item Cell: $(x,y)$ in pixels --- raw coordinates; -\item Chip: $(X,Y)$ in pixels --- the location on the silicon; -\item Focal Plane: $(p,q)$ in microns --- the location on the focal plane; -\item Tangent Plane: $(l,m)$ in arcsec from the telescope boresight; and -\item Sky: (RA,Dec) --- ICRS. -\end{itemize} - -The following steps are required to convert from the cell coordinates to -the sky: -\begin{itemize} -\item Cell $\longleftrightarrow$ Chip: two 2D polynomials, $(X,Y) = f(x,y)$; -\item Chip $\longleftrightarrow$ FP: two 2D polynomials, $(p,q) = g(X,Y)$; -\item FP $\longleftrightarrow$ TP: two 4D polynomials, $(l,m) = -h(p,q,m,c)$, where $m$ and $c$ are the magnitude and color of the -object, respectively; and -\item TP $\longleftrightarrow$ Sky: transformation to the sky using -pre-computed coefficients for each pointing. -\end{itemize} - -Note that the transformation between the Focal Plane and the Tangent -Plane is a four-dimensional polynomial, in order to account for any -possible dependencies in the astrometry on the stellar magnitude and -color; the former serves as a check for charge transfer -inefficiencies, while the latter will correct chromatic refraction, -both through the atmosphere and the corrector lenses. - -We require structures to contain each of the above transformations as -well as the pixel data. - -\subsubsection{Position Finding} - -We require functions to return the structure containing given -coordinates. For example, we want the chip that corresponds to the -focal plane coordinates $(p,q) = (-1.234,+5.678)$. These routines -handle the one-to-many problem --- i.e., for one given focal plane -coordinate, there are many chips that this coordinate may be -correspond to; these functions will select the correct one. -% -\begin{verbatim} -psCell *psCellInFPA (const psPlane *coord, const psFPA *fpa); -psChip *psChipInFPA (const psPlane *coord, const psFPA *fpa); -psCell *psCellInChip(const psPlane *coord, const psChip *chip); -\end{verbatim} - -\subsubsection{Conversion Functions} - -We require functions to convert between the various coordinate frames -(Section~\ref{sec:coordinateFrames}). The hierarchy of the coordinate -frames and the transformations between each are shown in -Figure~\ref{fig:coco}. The functions that employ the transformations -are shown in Figure~\ref{fig:cocoFunc}. In addition to -transformations between each adjoining coordinate frame in the -hierarchy, we also require higher-level functions to convert between -the Cell and Sky coordinate frames; these will simply perform the -intermediate steps. - -\begin{figure} -\psfig{file=coordinateFrames,height=7in,angle=-90} -\caption{The coordinate systems in the \PS{} IPP, and the relation -between each by transformations contained in the appropriate -structures.} -\label{fig:coco} -\end{figure} - -\begin{figure} -\psfig{file=coordinateConv,height=7in,angle=-90} -\caption{Conversion between coordinate systems by PSLib.} -\label{fig:cocoFunc} -\end{figure} - -We specify the following functions to convert between coordinates in -one type of frame to another type of frame. The first group consist -of unambiguous transformations: from the coordinates in a low-level -frame to the coordinates in the containing higher-level frame, of -which only one exists. In all of these functions, the output -coordinate structure may be \code{NULL} or may be supplied by the -calling function. In the former case, the structure must be -allocated; in the latter case, the supplied structure must be used. - -\begin{verbatim} -psPlane *psCoordCellToChip (psPlane *out, const psPlane *in, const psCell *cell); -% astrometry comes from cell (no need for parent) -\end{verbatim} -which converts coordindates \code{in} on the specified \code{cell} to -the coordinates on the parent chip. - -\begin{verbatim} -psPlane *psCoordChipToFPA (psPlane *out, const psPlane *in, const psChip *chip); -% astrometry comes from chip (no need for parent) -\end{verbatim} -which converts the coordinates \code{in} on the specified \code{chip} -to the coordinates on the parent FPA. - -\begin{verbatim} -psPlane *psCoordFPAToTP(psPlane *out, const psPlane *in, float color, float mag, const psFPA *fpa); -% astrometry comes from FPA (no need for parent) -\end{verbatim} -which converts coordinates \code{in} on the specified focal plane -\code{fpa} to tangent plane coordinates, applying the appropriate -distortion terms. The \code{color} and magnitude (\code{mag}) of the -source is necessary in order to perform the distortion between the -focal plane and the tangent plane. - -\begin{verbatim} -psSphere *psCoordTPToSky(psSphere *out, const psPlane *in, const psGrommit *grommit); -\end{verbatim} -which converts the tangent plane coordinates \code{in} to (RA,Dec) on -the sky, based on the environmental information specified by -\code{grommit}. - -\begin{verbatim} -psPlane *psCoordCellToFPA(psPlane *out, const psPlane *in, const psCell *cell); -% astrometry comes from cell -\end{verbatim} -which performs the single-step conversion between Cell coordinates -\code{in} and FPA coordinates. - -\begin{verbatim} -psSphere *psCoordCellToSky(psSphere *out, const psPlane *in, float color, float mag, const psCell *cell); -% astrometry comes from cell,chip,fpa (PARENT IS NEEDED HERE) -\end{verbatim} -which converts coordinates on the specified cell to (RA,Dec). This -transformation must be performed using the intermediate stage -transformations of Cell to Chip, Chip to FPA, FPA to Tangent Plane, -Tangent Plane to Sky. The information needed for each of these -transformations is available in the \code{.parent} elements of -\code{psCell} and \code{psChip}, and the \code{psFPA.exposure} -element. The \code{color} and magnitude (\code{mag}) of the source is -necessary in order to perform the distortion between the focal plane -and the tangent plane. - -\begin{verbatim} -psSphere *psCoordCellToSkyQuick(psSphere *out, const psPlane *in, const psCell *cell); -% astrometry comes from cell (no need for parent) -\end{verbatim} -which uses the 'quick-and-dirty' transformation to convert coordinates -on the specified cell to (RA,Dec). This transformation should use the -locally linear transformation specified by the element -\code{psCell.toTP}. Although the accuracy of this transformation -is lower than the complete transformation above, the calculation is -substantially faster as it only involves linear transformations. - -The following functions convert from high-level frames to the -coordinates of contained lower-level frames. - -\begin{verbatim} -psPlane *psCoordSkyToTP(psPlane *out, const psSphere *in, const psGrommit *grommit); -\end{verbatim} -which converts (RA,Dec) coordinates \code{in} to tangent plane coords -based on the enviromental information supplied by \code{grommit}. - -\begin{verbatim} -psPlane *psCoordTPToFPA(psPlane *out, const psPlane *in, float color, float mag, const psFPA *fpa); -\end{verbatim} -which converts the tangent plane coordinates \code{in} to focal plane -coordinates. The \code{color} and magnitude (\code{mag}) of the -source is necessary in order to perform the distortion between the -focal plane and the tangent plane. - -\begin{verbatim} -psPlane *psCoordFPAToChip (psPlane *out, const psPlane *in, const psChip *chip); -\end{verbatim} -which converts the specified FPA coordinates \code{in} to the -coordinates on the given Chip. The specified chip need not contain -the input coordinate. To find the chip which contains a particular -coordinate, the function \code{psChipInFPA}, defined above, should be -used. - -\begin{verbatim} -psPlane *psCoordChipToCell (psPlane *out, const psPlane *in, const psCell *cell); -\end{verbatim} -which converts the specified Chip coordinate \code{in} to the -coordinate on the given Cell. The specified Cell need not contain the -input coordinate. To find the cell which contains a particular -coordinate, the function \code{psCellInChip}, defined above, should be -used. - -\begin{verbatim} -psPlane *psCoordSkyToCell(psPlane *out, const psSphere *in, float color, float mag, psCell *cell); -\end{verbatim} -which directly converts (RA,Dec) \code{in} to coordinates on the -specified cell. The specified cell need not contain the input -coordinates. The \code{color} and magnitude (\code{mag}) of the -source is necessary in order to perform the distortion between the -focal plane and the tangent plane. - -\begin{verbatim} -psPlane *psCoordSkyToCellQuick(psPlane *out, const psSphere *in, psCell *cell); -\end{verbatim} -which directly converts (RA,Dec) \code{in} to coordinates on the -specified cell. The specified cell need not contain the input -coordinates. This transformation should use the locally linear -transformation specified by the element \code{psCell.toTP}. -Although the accuracy of this transformation is lower than the -complete transformation above, the calculation is substantially faster -as it only involves linear transformations. - -\subsubsection{Additional functions} - -We require additional functions to perform general functions which -will be useful for astrometry. Given coordinates on the sky, we -need to get the airmass, the parallactic angle, and an estimate of -the atmospheric refraction. - -\begin{verbatim} -float psGetAirmass(const psSphere *coord, psTime *lst, float height); -\end{verbatim} -which returns the airmass for a given position and local sidereal time -(\code{lst}). - -\begin{verbatim} -float psGetParallactic(const psSphere *coord, double siderealTime); -\end{verbatim} -which returns the parallactic angle for a given position and sidereal time. - -\begin{verbatim} -float psGetRefraction(float colour, ///< Colour of object - psPhotSystem colorPlus, ///< Colour reference - psPhotSystem colorMinus, ///< Colour reference - const psExposure *exp); ///< Telescope pointing information -\end{verbatim} -which provides an estimate of the atmospheric refraction, along the parallactic angle. - -\begin{verbatim} -double psGetParallaxFactor(const psExposure *exp) -\end{verbatim} -Calculate the parallax factor for the given exposure. \tbd{Why do we need this?}. - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\subsection{Astrometry and World Coordinate System} - -The FITS World Coordinate System (WCS) headers are commonly employed -with astronomical images in order to relate pixels to celestial (or -otherwise) coordinates. Since it is a FITS standard, we must be able -to read and write from WCS into our internal format. For the time -being, we will consider only celestial WCS (i.e., no spectral -wavelength calibrations, etc). Because WCS does not support the -multiple layers that we have built for \PS{}, we will use a simple -internal representation: a transformation, which handles any -distortions (i.e., goes directly from the coordinate frame of the -image to the tangent plane); and the projection. - -\begin{verbatim} -bool psAstrometryReadWCS(psPlaneTransform **transform, // Output transformation - psProjection **projection, // Output projection - psMetadata *header // Input FITS header - ); -bool psAstrometryWriteWCS(psMetadata *header, // Output FITS header - psPlaneTransform *transform, // Input transformation - psProjection *projection, // Input projection - double color, // Mean color to use - double magnitude, // Mean magnitude to use - ); -bool psAstrometrySimplify(psPlaneTransform **transform, // Output transformation - psProjection **projection, // Output projection - psCell *cell // Cell for which to generate transform and projection - ); -\end{verbatim} - -\code{pmReadAstrometry} shall parse the specified FITS \code{header}, -returning new instances of the \code{transform} and \code{projection} -that represent the WCS. The function shall return \code{true} if it -was able to successfully generate the outputs; otherwise it shall -return \code{false}. - -\code{pmWriteAstrometry} shall add WCS keywords to the supplied FITS -\code{header} that implement the given \code{transform} and -\code{projection}. The function shall return \code{true} if it was -able to successfully generate the output; otherwise it shall return -\code{false}. - -\code{pmSimplifyAstrometry} shall take a \code{cell} and simplify the -internal astrometric representation (\code{cell->toFPA} or equivalent, -\code{cell->parent->parent->toTangentPlane} and -\code{cell->parent->parent->grommit}) to a single \code{transform} and -\code{projection}. This allows the subsequent use of -\code{pmWriteAstrometry} in the case that we have only the -multi-layered \PS{} internal astrometric representation. The function -shall return \code{true} if it was able to successfully generate the -output; otherwise it shall return \code{false}. -