IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 335


Ignore:
Timestamp:
Mar 31, 2004, 3:54:50 PM (22 years ago)
Author:
Paul Price
Message:

Minor edits.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/doc/pslib/psLibADD.tex

    r318 r335  
    4343``heapsort'', ``quicksort'', and ``mergesort'' (see, e.g.,
    4444\citealt{knuth}, \citealt{sedgewick}, \citealt{press}).  These three
    45 sorting algorithms all run in a time of $O(n \log n)$ in the best
    46 case, but have different worse-case run times.  Implementations of all
    47 three sorting algorithms are described in references above and in
    48 various places online (e.g.,
     45sorting algorithms all run in a time of $O(n \log n)$ on the average,
     46but have different worse-case run times.  Implementations of all three
     47sorting algorithms are described in references above and in various
     48places online (e.g.,
    4949http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm).
    50 The linux \code{qsort} function orignially mades use of the quicksort
    51 algorithm, but now uses \tbd{what algorithm?}.  The GSL function
    52 \code{gsl_heapsort} uses the heapsort algorithm.  Both of these
    53 implemenations require a pointer to the comparison function, which may
    54 result in a slower code. 
    55 
    56 For PSLib, the sorting functions shall be implemented via
     50The linux \code{qsort} function uses a heapsort if it can allocate a
     51sufficiently large temporary buffer, and otherwise resorts to a
     52quicksort in-place.  The GSL function \code{gsl_heapsort} uses the
     53heapsort algorithm.  Both of these implemenations require a pointer to
     54the comparison function, which may result in a slower code than if the
     55comparison were done within the sort function.
     56
     57For PSLib, the sorting functions shall be implemented via the system
    5758\code{qsort}.  The function \code{psSort} shall return the sorted
    5859result of the input array without over-writing the input array.  The
    … …  
    7374neighbors in the form:
    7475
    75 \[ g_i = \sum_{n=-N}^N c_n f_i \]
     76\begin{equation}
     77 g_i = \sum_{n=-N}^N c_n f_i
     78\end{equation}
    7679%
    7780where the values of $c_n$ determine the filter type.  For boxcar
    … …  
    8083the ends of the data range to reduce the value of $c_n$ as fewer input
    8184data points may be used).  For Gaussian smoothing, the crucial
    82 parameter is $\sigma$, the standard deviation.  The
    83 value of $N$ should be chosen to be large enough, $N = 5\sigma$, and
    84 the values of $c_n$ are then just the Gaussian curve:
    85 
    86 \[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \]
     85parameter is $\sigma$, the standard deviation.  The value of $N$
     86should be chosen to be large enough to sample the Gaussian, $N =
     875\sigma$, and the values of $c_n$ are then just the Gaussian curve:
     88
     89\begin{equation}
     90c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}}
     91\end{equation}
    8792
    8893\subsection{Statistics}
    … …  
    97102use may vary depending on the context.  In addition, certain
    98103sigma-clipped values are defined as an intermediate choice between the
    99 sample and robust estimators. 
     104sample and robust estimators.
    100105
    101106\subsubsection{Sample Statistics}
    … …  
    106111\paragraph{Mean}
    107112
    108 The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \]
     113The mean is defined as:
     114\begin{equation}
     115\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i
     116\end{equation}
    109117
    110118\paragraph{Median}
    … …  
    115123number of values is odd, and by the average of the two middle values
    116124if the number of values is even.  This median should be avoided for
    117 samples which are large (e.g., $N > 10000$ elements) as the basic
     125samples which are large (e.g., $N > 10^4$ elements) as the basic
    118126robust median is quicker and more accurate. 
    119127
    … …  
    129137lower quartile, it is sufficient to provide the closest integer entry
    130138to these values.  The sample quartiles should be avoided for samples
    131 which are large (e.g., $N > 10000$ elements) as the robust quartiles
     139which are large (e.g., $N > 10^4$ elements) as the robust quartiles
    132140are quicker and more accurate.
    133141
    … …  
    136144The standard deviation of the sample is given by:
    137145
    138 \[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \]
     146\begin{equation}
     147\sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}}
     148\end{equation}
    139149
    140150To minimize the numerical rounding error, this should be calculated
    141151numerically as:
    142152
    143 \[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \]
     153\begin{equation}
     154\sigma = \sqrt{\frac{1}{N - 1} \left[ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} \left(\sum_{i = 1}^{N} (x_i - \bar{x})\right)^2 \ \right]}
     155\end{equation}
    144156
    145157\subsubsection{Clipped Statistics}
    146158
    147 The clipped statistics are used to determine the mean and $\sigma$ of
    148 a distribution in the presence of outliers. The clipped statistics are
    149 quicker than the complete robust statistical estimators and more
    150 reliable than the sample statistics.  The clipped statistics algorithm
    151 produces the clipped mean and clipped standard deviation.  The
    152 algorithm has 2 parameters: $N$, the number of iterations and $k$, the
    153 multiplying factor of $\sigma$ used to exclude outliers.  Typical
    154 values for both $N$ and $k$ are 3.  The algorithm is as follows:
     159The clipped statistics are used to determine the mean and standard
     160deviation of a distribution in the presence of outliers. The clipped
     161statistics are quicker than the complete robust statistical estimators
     162and more reliable than the sample statistics.  The clipped statistics
     163algorithm produces the clipped mean and clipped standard deviation.
     164The algorithm has 2 parameters: $N$, the number of iterations and $k$,
     165the multiplying factor of the standard deviation used to exclude
     166outliers.  Typical values for both $N$ and $k$ are 3.  The algorithm
     167is as follows:
    155168
    156169\begin{enumerate}
    … …  
    162175\item Repeat the following N times:
    163176  \begin{enumerate}
    164   \item Exclude all values $x_{i}$ for while $|x_{i} - \bar{x}| > k \sigma$.
     177  \item Exclude all values $x_{i}$ for which $|x_{i} - \bar{x}| > k \sigma$.
    165178  \item Compute the mean and standard deviation of the sub-sample.
    166179  \item Use the new mean for $\bar{x}$.
    167180  \item Use the new standard deviations for $\sigma$.
    168181\end{enumerate}
    169 \item The last calulcated value of $\bar{x}$ is the clipped mean.
    170 \item The last calulcated value of $\sigma$ is the clipped standard deviation.
     182\item The last calculated value of $\bar{x}$ is the clipped mean.
     183\item The last calculated value of $\sigma$ is the clipped standard
     184  deviation.
    171185\end{enumerate}
    172186
    173187\subsubsection{Robust Statistics}
    174188
    175 The robust version of the statistics provide estimators of basic
     189The robust version of the statistics provides estimators of basic
    176190statistical concepts which are reliable even for data samples with
    177191significant contamination.  A typical case is the situation in which
    … …  
    202216  4$. 
    203217\item Find the bin with the peak value in the range $L_{\frac{1}{4}}$
    204   to $U_{\frac{1}{4}}$; this is the robust mode $\mbox{mode}_r$. 
     218  to $U_{\frac{1}{4}}$; this is the robust mode, $\mbox{mode}_r$. 
    205219\item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 8$.
    206220\item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to
    … …  
    242256diagonal matrices that satisfy the relationship $A = L U$ where $L$ is
    243257a lower-diagonal matrix of the form:
    244 \[ L = \left(
     258\begin{equation}
     259L = \left(
    245260\begin{matrix}
    246261\alpha_{11} & 0           & 0           & 0 \\
    … …  
    249264\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\
    250265\end{matrix} \right)
    251 \]
     266\end{equation}
    252267%
    253 (where all diagonal values $\alpha){ii} = 1$) and $U$ is an upper-diagonal matrix of the form:
    254 \[ U = \left(
     268(where all diagonal values $\alpha_{ii} = 1$) and $U$ is an upper-diagonal matrix of the form:
     269\begin{equation}
     270U = \left(
    255271\begin{matrix}
    256272\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\
    … …  
    2592750          & 0          & 0          & \beta_{44} \\
    260276\end{matrix} \right)
    261 \]
    262 %
     277\end{equation}
     278
    263279We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following
    264280this procedure.  First, $\alpha_{ii} = 1$.  For all values of $j = 1,
    2652812, \dots, N$, follow the next steps: For each value of $i = 1, 2,
    266282\dots, j$, solve for $\beta_{ij}$ using the relationship:
    267 \[ \beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} \]
     283\begin{equation}
     284\beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj}
     285\end{equation}
    268286For the values of $i = j+1, j+2, \dots, N$, solve for the values of
    269287$\alpha_{ij}$ with the following:
    270288%
    271 \[ \alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) \]
     289\begin{equation}
     290\alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right)
     291\end{equation}
    272292
    273293\subsubsection{Calculate a matrix determinant}
    … …  
    275295The determinant $D$ of a matrix $a_{ij}$ is calculated from the
    276296product of the diagonal elements of the LU decomposition:
    277 \[ D = P_{i=1}{N} U_{ii} \]
     297\begin{equation}
     298D = P_{i=1}{N} U_{ii}
     299\end{equation}
    278300
    279301This shall be calculated using the GSL function
    … …  
    293315$\beta_{ij}$, the technique of back-substitution is used to solve the
    294316equation above. First solve the equation $L y = B$:
    295 \[ y_1 = \frac{b_1}{\alpha_{11}}\]
    296 \[ y_i = \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right) \]
    297 %
     317\begin{eqnarray}
     318y_1 & = & \frac{b_1}{\alpha_{11}} \\
     319y_i & = & \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right)
     320\end{eqnarray}
     321
    298322The values of $y$ may be used to solve for $x_i$ using the
    299323relationship $U x = b$:
    300 \[ x_N = \frac{y_N}{\beta_{NN}}\]
    301 \[ x_i = \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) \]
     324\begin{eqnarray}
     325x_N & = & \frac{y_N}{\beta_{NN}} \\
     326x_i & = & \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right)
     327\end{eqnarray}
    302328
    303329\subsubsection{Invert a matrix}
    … …  
    313339
    314340Matrix binary arithmetic operations differ from image binary
    315 arithmetic operations in a fundamental way:  For image operations, the
     341arithmetic operations in a fundamental way: For image operations, the
    316342resulting pixel value is determined by performing the operation on the
    317343two corresponding input operand pixels.  For matrix operations, the
    318344resulting element value results from the operation on the
    319 corresponding pair of row and colum from the input matrix.  So, for
     345corresponding pair of row and colum from the input matrices.  So, for
    320346example, given the input matrices $\alpha_{ij}$ and $\beta_{ij}$, the
    321 image and matrix multiplication correspond to:
     347image and matrix multiplications correspond to:
    322348%
    323 \[ \mbox{image}_{ij} = \alpha_{ij} \times \beta_{ij} \]
    324 \[ \mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki} \]
     349\begin{eqnarray}
     350\mbox{image}_{ij} & = & \alpha_{ij} \times \beta_{ij} \\
     351\mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki}
     352\end{eqnarray}
    325353%
    326 The matrix math function \code{psMatrixOp} shall implement the matrix
    327 version of the binary arithmetical operations for the following
    328 operators: $+, -, \times, /$
     354The matrix and image operations for addition and subtraction are
     355identical: the operation is performed on the corresponding elements in
     356each matrix.  Matrix division is not defined (to divide, the user
     357should first invert the matrix, and then multiply).  The matrix math
     358function \code{psMatrixOp} shall implement the matrix version of the
     359binary arithmetical operations for the following operators: $+, -,
     360\times$.
    329361
    330362\subsubsection{Transpose a matrix}
    … …  
    336368given by
    337369%
    338 \[ T_{ij} = M_{ji} \]
     370\begin{equation}
     371T_{ij} = M_{ji}
     372\end{equation}
    339373where $M_{ij}$ is the matrix to be transposed.
    340374
    341375\subsubsection{Convert a matrix to a vector}
    342 \TBD{write me}
     376
     377Matrix-to-vector conversion is only defined for a matrix that has a
     378size of one in at least one dimension.  In that case, the elements
     379should be extracted, and placed in a vector, with care that the
     380\code{psType} is defined correctly --- a $1\times N$ matrix is
     381converted to a \code{PS_DIMEN_VECTOR}-type vector, while a $N\times 1$
     382matrix is converted to a \code{PS_DIMEN_TRANV}-type vector.
    343383
    344384\subsection{Fitting}
    345385
     386\subsubsection{Chi-squared}
     387\label{chisq}
     388
     389Given a set of N ordinates, $x_i$, measured coordinates, $y_i$, with
     390errors, $\sigma_i$, and a model function, $f(x_i)$, then $\chi^2$
     391(``chi-squared'') is defined:
     392
     393\begin{equation}
     394\chi^2 = \sum_{i=0}^{N-1} \left( \frac{y_i - f(x_i)}{\sigma_i} \right)^2
     395\end{equation}
     396
    346397\subsubsection{General Polynomial Fitting}
    347 
    348 \TBD{define for use with psGetArrayPolynomial}
    349 
    350 \TBD{change psGetArrayPolynomial to return a psPolynomial}
    351398
    352399Given a set of data values $y_i$ with errors $\sigma_i$, related to
    … …  
    354401of $N_{par}$ free parameters of the form $y = \sum_{j=0}^{N_{par}}
    355402\alpha_j x^j$.  The model is determined by minimizing the value of
    356 $\chi^2$ as defined above (\ref{chisq}), and the minimization is
    357 determined by solving for the set of parameters $\alpha_j$ for which
    358 the partial derivatives of the $\chi^2$ with respect to those
    359 parameters are zero.  The partial derivatives are
    360 
    361 \[ \frac{\partial \chi^2}{\partial \alpha_k} = -2 \sum_i (y_i - \sum_j x_i^j \alpha_j) \frac{x_i^k}{\sigma_i^2} \]
     403$\chi^2$ (\ref{chisq}), and the minimization is determined by solving
     404for the set of parameters $\alpha_j$ for which the partial derivatives
     405of the $\chi^2$ with respect to those parameters are zero.  The
     406partial derivatives are
     407
     408\begin{equation}
     409\frac{\partial \chi^2}{\partial \alpha_k} = -2 \sum_i (y_i - \sum_j x_i^j \alpha_j) \frac{x_i^k}{\sigma_i^2}
     410\end{equation}
    362411
    363412Setting the derivatives to zero, this may be reduced to the following
    364413matrix equation:
    365414
    366 \[ \sum_j \alpha_j \sum_i \frac{x_i^k x_i^j}{\sigma_i^2} = \sum_i \frac{x_i^k y_i}{\sigma_i^2} \]
     415\begin{equation}
     416\sum_j \alpha_j \sum_i \frac{x_i^k x_i^j}{\sigma_i^2} = \sum_i \frac{x_i^k y_i}{\sigma_i^2}
     417\end{equation}
    367418
    368419This matrix equation may be solved with LU Decomposition
    … …  
    375426For models in which the system of equations defined by the partial
    376427derivatives cannot be solved with the linear technique, other options
    377 are necessary.  The Levenberg-Marquardt Method (LMM; see Numerical
    378 Recipes for a discusion) may be used for these situations.  In LMM, we
    379 make a guess at the input parameters, measure the $\chi^2$, vary the
    380 parameters by a particular choice based on the gradient, measure the
    381 $\chi^2$ again, and adjust the parameters and the parameter varient
    382 based on the results.
     428are necessary.  The Levenberg-Marquardt Method (LMM; see NR \S 15.5)
     429may be used for these situations.  In LMM, we make a guess at the
     430input parameters, measure the $\chi^2$, vary the parameters by a
     431particular choice based on the gradient, measure the $\chi^2$ again,
     432and adjust the parameters and the parameter varient based on the
     433results.
    383434
    384435Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent
    … …  
    389440$\alpha_{j,k}$ at this parameter selection as follows:
    390441
    391 \[ \beta_k = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \]
    392 \[ \alpha_{j,k} = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2} \]
     442\begin{eqnarray}
     443\beta_k & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \\
     444\alpha_{j,k} & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2}
     445\end{eqnarray}
    393446%
    394447We now define the new parameter guess for $a_k$ based on the gradient
    … …  
    396449follows:
    397450
    398 \[ A_{j,k} = \alpha_{j,k} ~ (j \ne k) \]
    399 \[ A_{j,k} = (1 + \lambda) \alpha_{j,k} ~ (j = k) \]
     451\begin{eqnarray}
     452A_{j,k} & = & \alpha_{j,k} ~ (j \ne k) \\
     453A_{j,k} & = & (1 + \lambda) \alpha_{j,k} ~ (j = k)
     454\end{eqnarray}
    400455%
    401456and solve the system of equations represented by:
    402 \[ A_{j,k} \alpha^\prime_k = \beta_j \]
     457\begin{equation}
     458A_{j,k} a^\prime_k = \beta_j
     459\end{equation}
    403460%
    404 where $\alpha^\prime_k$ represents our new attempt at a parameter
     461where $a^\prime_k$ represents our new attempt at a parameter
    405462guess. We use this parameter set to calculate $\chi^2$.  If the new
    406463value of $\chi^2$ is lower than the previous guess, we accept this new
    … …  
    437494\end{center}
    438495
     496%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     497
    439498\subsection{Polynomials}
    440499
    … …  
    446505\begin{itemize}
    447506\item They are bounded on $-1 < x < 1$, with the maxima and minima
    448 over this range being 1 and -1, respectively;
     507over this range being 1 and $-1$, respectively;
    449508\item Truncation of the higher-order terms leaves one with the most accurate
    450509lower-order polynomial representation of the desired function.
    … …  
    477536
    478537Multi-dimensional polynomials shall be composed of multiplications of
    479 1D Chebyshev polynomials.
    480 
     5381D Chebyshev polynomials, with the coefficients stored in tensors of
     539the appropriate rank.
     540
     541
     542%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    481543
    482544\subsection{(Fast) Fourier Transforms}
    … …  
    500562will be specified as a configuration variable for the IPP (defaulting
    501563to \code{/etc/fftw/wisdom}).  The \code{wisdom} should be read in upon
    502 initialisation of the PSLib FFT functions and saved at the conclusion.
     564initialisation of the PSLib FFT functions (\code{psFFTInit()}) and
     565saved at the conclusion (\code{psFFTDone()}).
    503566
    504567\subsubsection{Function mapping}
    … …  
    515578
    516579These plans should be formulated using the \code{FFTW_ESTIMATE} flag,
    517 which will allow FFTW to default to a short planning time if the
     580which will allow FFTW to default to a minimal planning time if the
    518581wisdom has not been loaded.  Transforms should be performed out of
    519582place to avoid the need to pad the input array to hold the output.
    … …  
    521584\subsubsection{More Complicated Functions}
    522585
    523 The \code{psFFTFilter} and \code{psFFTFilterComplex} functions provide
     586The \code{psFFTFilter()} and \code{psFFTFilterComplex()} functions provide
    524587the means to apply a filter (purely real and imaginary multipliers,
    525588respectively) to the data in the Fourier plane.  If the filter
    526 function specified for \code{psFFTFilter} returns a real value, $r$,
     589function specified for \code{psFFTFilter()} returns a real value, $r$,
    527590then the corresponding value in the Fourier plane should be multiplied
    528591by $r$.  If the real and imaginary filter functions specified for
    529 \code{psFFTFilterComplex} return the values $r$ and $s$, respectively,
     592\code{psFFTFilterComplex()} return the values $r$ and $s$, respectively,
    530593then the corresponding value in the Fourier plane should be multiplied
    531594by the complex number $r + si$.
    … …  
    535598Fourier transform is multiplied by the complex conjugate of the second
    536599Fourier transform to yield the Fourier transform of the
    537 cross-correlation (NR 13.2).  In the latter, the two Fourier
     600cross-correlation (NR \S 13.2).  In the latter, the two Fourier
    538601transforms are multiplied directly to yield the Fourier transform of
    539 the convolution (NR 13.1).
     602the convolution (NR \S 13.1).
    540603
    541604If the elements of the discrete Fourier transform are $C_k$, then the
    542 the elements of the power spectrum are (NR 13.4):
     605the elements of the power spectrum are (NR \S 13.4):
    543606\begin{eqnarray}
    544607P_0     & = & \left| C_0 \right|^2 / N^2 \\
    … …  
    550613Note that we leave the issue of ``windowing'' the data up to the
    551614caller, and choose to normalise by $1/N^2$.
     615
     616%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    552617
    553618\subsection{Astronomy Utilities}
    … …  
    557622\href{http://star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{SLALIB
    558623Positional Astronomy Library}.
     624
     625%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    559626
    560627\subsubsection{Celestial Coordinate Conversions}
    … …  
    571638\end{tabular}
    572639
     640%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     641
    573642\subsubsection{Projections}
    574643
    … …  
    584653
    585654\begin{eqnarray}
    586 x & = & R \sin \phi \\
    587 y & = & -R \cos \phi
     655x & = & R \sin (\phi) \\
     656y & = & -R \cos (\phi)
    588657\end{eqnarray}
    589658
    … …  
    600669
    601670\begin{eqnarray}
    602 R & = & \cot \theta 180^\circ/\pi \\
     671R & = & \cot (\theta) 180^\circ/\pi \\
    603672\theta & = & \arctan (180^\circ/(\pi R))
    604673\end{eqnarray}
    … …  
    609678
    610679\begin{eqnarray}
    611 R & = & \cos \theta 180^\circ/\pi \\
     680R & = & \cos (\theta) 180^\circ/\pi \\
    612681\theta & = & \arccos (\pi R / 180^\circ)
    613682\end{eqnarray}
    … …  
    637706
    638707\begin{eqnarray}
    639 x & = & 2 \alpha \cos \theta \sin (\phi/2) \\
     708x & = & 2 \alpha \cos (\theta) \sin (\phi/2) \\
    640709y & = & \alpha \sin \theta \\
    641 {\rm where}\hspace{1cm} \alpha^{-1} & \equiv & (180^\circ/\pi) \sqrt{\left(1 + \cos \theta \cos (\phi/2) \right) / 2}
     710{\rm where}\hspace{1cm} \alpha^{-1} & \equiv & (180^\circ/\pi) \sqrt{\left(1 + \cos (\theta) \cos (\phi/2) \right) / 2}
    642711\end{eqnarray}
    643712
    … …  
    649718{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x\pi/720^\circ)^2 - (y\pi/360^\circ)^2}
    650719\end{eqnarray}
     720
     721%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     722
     723\subsection{Tangent Plane to Sky}
     724
     725Mappings between the tangent plane and the sky will be implemented
     726using SLALIB.
     727
     728To speed up the transformations, SLALIB allows the storage of
     729``apparent-to-observed parameters'', which will be contained in a
     730\code{psGrommit}.  The \code{psGrommit} consists of the following:
     731\begin{enumerate}
     732\item geodetic latitude (radians)
     733\item sine and cosine of geodetic latitude
     734\item magnitude of diurnal aberration vector
     735\item height (metres)
     736\item ambient temperature (degrees K)
     737\item pressure (mB)
     738\item relative humidity (0--1)
     739\item wavelength ($\mu$m)
     740\item lapse rate (degrees K per metre)
     741\item refraction constants A and B (radians)
     742\item longitude + eqn of equinoxes + ``sidereal $\Delta$ UT'' (radians)
     743\item local apparent sidereal time (radians)
     744\end{enumerate}
     745These may be calculated using \code{sla_AOPPA}.  Note that a
     746\code{psGrommit} is only appropriate for a single exposure.
     747
     748Once the \code{psGrommit} has been calculated, the functions
     749\code{sla_OAPQK} and \code{sla_AOPQK} convert from the tangent plane
     750to the sky, and the sky to the tangent plane, respectively.  (Note
     751that ``observed'' in SLALIB refers to the tangent plane, and
     752``apparent'' refers to the apparent position on the sky.)
     753
     754%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     755
     756\subsection{The One-to-Many Problem with Mosaic Cameras}
     757
     758The \PS{} focal plane consists of several chips, so we will often want
     759to identify which chip a source lies on, when we know the coordinates
     760in the focal plane.  This is an example of the one-to-many problem
     761(one coordinate, many chips that it may lie on).
     762
     763If this needs to be repeated for only one (or a small number of) focal
     764plane coordinates, then the fastest method is to simply convert the
     765focal plane coordinates to chip coordinates for each of the chips, and
     766determine for which of the chips the chip coordinates are valid (i.e.\
     767on the chip).
     768
     769On the other hand, if this needs to be repeated for many source focal
     770plane coordinates, then it is most efficient to convert the centers of
     771each of the chips to coordinates on the focal plane and to use the
     772distance of the source to each of the centers to optimise which chips
     773are tested first.  This saves testing many chips for every source.
     774
     775
     776%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     777
     778\subsection{General Astronomy Functions}
     779
     780The airmass is calculated using the SLALIB function \code{sla_AIRMAS}.
     781
     782The parallactic angle is calculated using the SLALIB function \code{sla_PA}.
     783
     784%The parallax factors are calculated using the following formulae
     785%(Smart et al.\ 2003, A\&A, 404, 317):
     786%\begin{eqnarray}
     787%P_\xi & = & \cos \alpha \sin \lambda \cos \epsilon - \sin \alpha \cos \lambda \\
     788%P_\eta & = & (\sin \epsilon \cos \delta - \cos \epsilon \sin \alpha \sin \delta) \sin \lambda - \cos \alpha \sin \delta \cos \lambda
     789%\end{eqnarray}
     790%where $\alpha$ is the Right Ascension, $\delta$ is the Declination,
     791%$\lambda$ is the solar longitude, and $\epsilon = 23^\circ 27'08''.26$
     792%is the inclination of the ecliptic.  The solar longitude is obtained
     793%from the ecliptic coordinates of the Sun.
     794
     795To calculate the parallax factors, get the mean-to-apparent parameters
     796(\code{sla_MAPPA}) for a mean epoch of 2000.0, and, given the mean
     797position of interest, calculate the apparent position
     798(\code{sla_MAPQK}) for a parallax of 1.0 arcsec $(\alpha_1,\delta_1)$,
     799and a parallax of 0.0 arcsec $(\alpha_0,\delta_0)$.  Then the parallax
     800factors in radians are:
     801\begin{eqnarray}
     802P_x & = & 3,600 (180^\circ/\pi) (\alpha_1 - \alpha_0) cos (\delta_0) \\
     803P_y & = & 3,600 (180^\circ/\pi) (\delta_1 - \delta_0)
     804\end{eqnarray}
     805
     806%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     807
     808\subsection{Positions of Major Solar System Objects}
     809
     810The SLALIB function \code{SLA_RDPLAN} returns the apparent position of
     811a specified planet, or the Moon.
     812
     813To calculate the position of the Sun, use \code{sla_EVP} to get the
     814position of the earth relative to the Sun, and convert from the
     815cartesian coordinates to spherical using \code{sla_DCC2S}, and
     816calculate the position on the opposite side of the sphere ($\alpha
     817\rightarrow \alpha + 12 {\rm hrs}$ and $\delta \rightarrow -\delta$).
     818
     819%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     820
     821\subsection{Offsets}
     822
     823
     824
    651825
    652826
Note: See TracChangeset for help on using the changeset viewer.