Index: trunk/doc/pslib/psLibADD.tex =================================================================== --- trunk/doc/pslib/psLibADD.tex (revision 1554) +++ trunk/doc/pslib/psLibADD.tex (revision 1555) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.25 2004-08-17 01:56:41 eugene Exp $ +%%% $Id: psLibADD.tex,v 1.26 2004-08-17 03:09:16 price Exp $ \documentclass[panstarrs]{panstarrs} @@ -450,6 +450,4 @@ \paragraph{Non-linear Fitting: Levenberg-Marquardt Method} - -\TBD{describe LMM for psMinimize and psMinimizeChi2} For models in which the system of equations defined by the partial @@ -523,4 +521,56 @@ \end{center} + +\paragraph{Non-linear fitting: Powell's method} + +Powell's method is a type of ``Direction Set'' methods in +multi-dimensions for finding a local minimum. Given a starting point +(the ``best guess'' for the minimum) and a set of direction vectors, a +direction set method advances in the direction of the vectors, +determines a new direction vector by some method, and proceeds in this +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. + +We will use for our algorithm the modified version of Powell's +Quadratically Convergent Method, which is described below, adapted +from NR. + +\begin{enumerate} +\item Given a function in $N$ dimensions to minimize, $f$, and a best + guess for the minimum, point $P$ in $N$ dimensions, take an initial + set of $N$ vectors, $v_i$, to be the unit vectors. +\item Set point $Q = P$. +\item For each dimension in turn, move $Q$ \textit{only} in the + direction $v_i$ to minimize the function of interest. +\item Set vector $u = Q - P$. +\item Move $Q$ \textit{only} in the direction $u$ +\item Replace the vector along which the largest minimization was + made, $v_{i,\rm max}$, with $u$, except under either of the + following circumstances: + \begin{itemize} + \item If $f_QP \ge f_P$, then there is no point in keeping the new + vector, because there is no further minimization to be made in + that direction. + \item If $2 ( f_P - 2f_Q + f_{QP} ) \left[ ( f_P - f_Q ) - + \Delta_{\rm max} \right]^2 \ge ( f_P - f_{QP} )^2 \Delta_{\rm max}$, + then either the decrease in the function was not due to any single + direction, or we are close to the minimum. + \end{itemize} + where $f_P = f(P)$, $f_Q = f(Q)$, $f_{QP} = f(2Q - P)$, and + $\Delta_{\rm max} \ge 0$ is the magnitude of the minimization made + along $v_{i,\rm max}$. +\item Set $P$ to $Q$. +\item Return to step 3 until the change in this last move is less + than some specified tolerance, or a maximum number of iterations + has been reached. +\end{enumerate} + +In regards to minimizing the function only in a particular direction, +we shall adopt, as NR recommends, bracketing the minimum before +applying Brent's method, \tbd{which will be specified in detail +later}. + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -644,17 +694,12 @@ Correct time representation is critical in astronomical software. -PSLib uses the \code{psTime} structure to represent all time values. -This structure represents a time which is consists of seconds and -fractions of seconds in a time system defined by the \code{psTimeType} -element \code{type}. Two possible time systems are currently -available: TAI and UTC. Both are defined in terms of the reference -epoch 1970-01-01T00:00:00Z, but with minor modifications for -leap-seconds as needed. The first represenatation, TAI (International -Atomic Time), has seconds of uniform length and no leap seconds. The -exact zero reference is 1970/01/01,00:00:10 UTC. The second -representations is UTC, which has seconds of uniform length and -leap-seconds as needed to adjust it to remain within 0.9 seconds of -the Earth's rotation. It has a zero-point of exactly -1970/01/01,00:00:00 UTC. +PSLib uses the \code{psTime} structure to represent all time +values. This structure represents a time which is equivalent to TAI +(International Atomic Time) and has the following properties: +\begin{itemize} +\item it represents both seconds and microseconds +\item the seconds are continuous (no leap seconds) +\item the zero reference point is \tbd{1970/01/01,00:00:10} UTC. +\end{itemize} Julian Day (JD) and Modified Julian Day (MJD) are both continuous time @@ -713,19 +758,4 @@ \end{verbatim} -The conversion from a time and longitude to local mean sidereal time -is performed using the SLA Lib function \code{sla_GMST}. This -function requires the value $\Delta$ UT1 = UTC - UT1. The value of -$\Delta$ UT1 may be determined from the following site in real time: - -\code{ftp://maia.usno.navy.mil/ser7/ser7.dat} - -In addition, the long-term values may be determined from the table -found at: \code{ftp://maia.usno.navy.mil/ser7/finals.all}. See also -the web page \code{http://maia.usno.navy.mil/}. The most significant -accuracy requirements are for the current value when calculating the -LST. For this purpose, the table above (\code{ser7.dat}), which -provides predictions over a 2 month period, must be made available -locally to PSLib and updated regularly. - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -862,136 +892,47 @@ \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 +The following information is from \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$. - -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 - \phi_p) -\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 arg} (-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$: +Calabretta (1995, ADASS, 4, 233)}. + +Let the latitude be $\phi$ and the longitude $\theta$. The domains of +these are $-\pi < \phi \le \pi$ and $-\pi/2 \le \theta \le \pi/2$. + +For zenithal projections (e.g.\ Gnomonic and Orthographic) the +following hold: + +\begin{eqnarray} +x & = & R \sin (\phi) \\ +y & = & -R \cos (\phi) +\end{eqnarray} + +and + +\begin{eqnarray} +R & = & \sqrt{x^2 + y^2} \\ +\phi & = & {\rm arg} (-y,x) +\end{eqnarray} + +\paragraph{Gnomonic} + +The Gnomonic projection (``TAN'') is a zenithal projection. + +\begin{eqnarray} +R & = & \cot (\theta) 180^\circ/\pi \\ +\theta & = & \arctan (180^\circ/(\pi R)) +\end{eqnarray} + +\paragraph{Orthographic} + +The Orthographic projection (``SIN'') is a zenithal projection. + +\begin{eqnarray} +R & = & \cos (\theta) 180^\circ/\pi \\ +\theta & = & \arccos (\pi R / 180^\circ) +\end{eqnarray} + +\paragraph{Cartesian} + +The Cartesian projection (``CAR'') is a very simple cylindrical projection. \begin{eqnarray} @@ -1000,5 +941,5 @@ \end{eqnarray} -\subparagraph{Mercator} +\paragraph{Mercator} The Mercator projection (``MER'') is a cylindrical projection. @@ -1006,16 +947,16 @@ \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)} +y & = & \ln \left( \tan (45^\circ + \theta/2) \right) 180^\circ/\pi \\ +{\rm and}\hspace{1cm} \theta & = & 2 \arctan \left( e^{y\pi/180^\circ} \right) - 90^\circ +\end{eqnarray} + +\paragraph{Hammer-Aitoff} + +The Hammer-Aitoff projection is a general projection, and is defined: + +\begin{eqnarray} +x & = & 2 \alpha \cos (\theta) \sin (\phi/2) \\ +y & = & \alpha \sin \theta \\ +{\rm where}\hspace{1cm} \alpha^{-1} & \equiv & (180^\circ/\pi) \sqrt{\left(1 + \cos (\theta) \cos (\phi/2) \right) / 2} \end{eqnarray} @@ -1023,24 +964,7 @@ \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 \\ +\phi & = & 2 {\rm arg} (2z^2 - 1, xz \pi/360^\circ) \\ +\theta & = & \arcsin (yz\pi/180^\circ) \\ +{\rm where}\hspace{1cm} z & \equiv & \sqrt{1 - (x\pi/720^\circ)^2 - (y\pi/360^\circ)^2} \end{eqnarray}