Index: trunk/doc/pslib/psLibADD.tex
===================================================================
--- trunk/doc/pslib/psLibADD.tex	(revision 6009)
+++ trunk/doc/pslib/psLibADD.tex	(revision 9785)
@@ -1,3 +1,3 @@
-%%% $Id: psLibADD.tex,v 1.93 2006-01-16 01:11:40 eugene Exp $
+%%% $Id: psLibADD.tex,v 1.94 2006-10-30 21:45:46 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 
@@ -366,13 +366,13 @@
 size is choosen to keep about 25 points per bin.  Thus, the bin size
 ($d\sigma$) is set to about:
-\[
+\begin{equation}
 d\sigma = 2.35 \sigma (N_{\rm bin}/N_{50}) = (25 \times 2.35) (\sigma/N_{50})
-\]
+\end{equation}
 With the limitation that $\sigma/d\sigma$ should be limited on one end
 to the value 1, and the other to the value 5.  The easiest way to set
 this limit is to define dN to be:
-\[
+\begin{equation}
 dN = (\sigma / d\sigma) = (N_{50} / N_{\rm bin}) / 2.35 = 0.017 * N_{50}
-\]
+\end{equation}
 \subsubsection{Histograms}
 
@@ -426,5 +426,5 @@
 T_2(x) & = & 2x^2 - 1 \\
 T_3(x) & = & 4x^3 - 3x \\
-T_4(x) & = & 8x^4 - 8x^2 + 1 \\
+T_4(x) & = & 8x^4 - 8x^2 + 1 
 \end{eqnarray}
 Chebyshev polynomials follow the recurrence relation:
@@ -437,5 +437,5 @@
 \begin{eqnarray}
 d_j  & = & 2xd_{j+1} - d_{j+2} + c_j \\
-f(x) & = & x*d_1 - d_2 + 1/2 c_0 \\
+f(x) & = & x*d_1 - d_2 + 1/2 c_0
 \end{eqnarray}
 
@@ -721,9 +721,9 @@
 $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 \]
+\begin{equation} 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 \end{equation}
 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 \]
+\begin{equation} r_x = V_{0,0} + (V_{1,0} - V_{0,0})f_x \end{equation}
+\begin{equation} V = r_x + (V_{0,1} + (V_{1,1} - V_{0,1})f_x - r_x)f_y \end{equation}
 
 Note that the values of $f_x$ and $f_y$ require some care.  Given a
@@ -876,6 +876,6 @@
 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 \]
+\begin{equation} x = (i - i_o)*\cos\theta + (j - j_o)*\sin\theta \end{equation}
+\begin{equation} y = (i_o - i)*\sin\theta + (j - j_o)*\cos\theta \end{equation}
 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
@@ -1092,5 +1092,5 @@
 P_0     & = & \left| C_0 \right|^2 / N^2 \\
 P_j     & = & \left( \left| C_j \right|^2 + \left| C_{N-j} \right|^2 \right)/ N^2 \\
-P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2 \\
+P_{N/2} & = & \left| C_{N/2} \right|^2 / N^2
 \end{eqnarray}
 where $j = 1, 2, \ldots, (N/2 - 1)$.
@@ -1098,4 +1098,156 @@
 Note that we leave the issue of ``windowing'' the data up to the
 caller, and choose to normalise by $1/N^2$.
+
+\subsection{Ellipse Representations}
+
+Images of astronomical objects may often be represented using a model
+consisting of a radial profile combined with an elliptical contour.
+Two common ways to measure such a shape are to fit a model to the
+light distribution or to measure second-order moments, perhaps with
+some weighting profile.  In the special case of a 2D Gaussian with an
+elliptical contour, these representations are equivalent.  The
+following discussion shows how to relate the fitted parameters and
+second-order moments of a elliptical Gaussian of an arbitrary
+orientation with the parameters of an unrotate elliptical Gaussian.
+
+Consider a 2D Gaussian with an elliptical contour.  If the ellipse is
+oriented with the major axis along the x-axis, then the formula for
+such a Gaussian may be written $f = exp(-z)$ where
+\begin{equation}
+\label{aligned-ellipse}
+z = \frac{x^2}{2\sigma_a^2} + \frac{y^2}{2\sigma_b^2}
+\end{equation}
+with $\sigma_a$ the semi-major axis and $\sigma_b$ the semi-minor axis
+of the 1$\sigma$ contour.  Given such a Gaussian, we may measure its
+second moments by integration, and find that the second moment tensor
+is
+\begin{equation} \left| \begin{array}{cc}
+\sigma_a^2 & 0 \\
+0 & \sigma_b^2 \\
+\end{array} \right| \end{equation}
+
+Now consider the same ellipse rotated to an arbitrary angle $\theta$.
+The formula for such a Gaussian may be written $f = exp(-z)$ where
+\begin{equation}
+\label{rotated-ellipse}
+z = \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + \sigma_{xy}xy
+\end{equation}
+Note that, in the above form of the equation, $\sigma_{xy}$ goes to 0
+as the ellipse is rotated to be aligned with the x (or y) axis.  Thus,
+in this representation, $\sigma_{xy}$ is well behaved, but does not
+have the same units of length that $\sigma_x$ or $\sigma_y$ have.  Our
+goal is to determine the relationships between the rotated and
+unrotated components of the Gaussian formula as well as the second
+moments.
+
+To determine the behavior of $\sigma_x$, etc, under rotation, we start
+with the aligned ellipse (\ref{aligned-ellipse}) and rotate the
+coordinate frame by an angle $-\theta$:
+\begin{equation} 
+\left( \begin{array}{c} x^\prime \\ y^\prime \end{array} \right) =
+\left| \begin{array}{cc} \cos \theta & \sin \theta \\ 
+                        -\sin \theta & \cos \theta 
+\end{array} \right|
+\left( \begin{array}{c} x \\ y \end{array} \right)
+\end{equation}
+where $x^\prime$ and $y^\prime$ are the coordinates for the unrotated
+(aligned) ellipse.  Applying this rotation to (\ref{aligned-ellipse}) yields:
+\begin{equation} 
+z = \frac{x^2 \cos^2 \theta + y^2 \sin^2 \theta + 2 x y \sin \theta \cos \theta}{2\sigma_a^2} +
+    \frac{x^2 \sin^2 \theta + y^2 \cos^2 \theta - 2 x y \sin \theta \cos \theta}{2\sigma_b^2} 
+\end{equation}
+Grouping these terms together, we find:
+\begin{equation} 
+z = \frac{x^2}{2}(\sigma_a^{-2} \cos^2 \theta + \sigma_b^{-2}\sin^2 \theta) + 
+    \frac{y^2}{2}(\sigma_b^{-2} \cos^2 \theta + \sigma_a^{-2}\sin^2 \theta) - 
+    \frac{xy}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
+\end{equation}
+We then associate the components of this equation with those of (\ref{rotated-ellipse}) and find:
+\begin{eqnarray}
+\sigma_x^{-2} & = & \sigma_a^{-2} \cos^2 \theta + \sigma_b^{-2}\sin^2 \theta \\
+\sigma_y^{-2} & = & \sigma_b^{-2} \cos^2 \theta + \sigma_a^{-2}\sin^2 \theta \\
+\sigma_{xy}   & = & \frac{1}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
+\end{eqnarray}
+Replacing $\cos^2$ and $\sin^2$ with the double-angle relationships, we find:
+\begin{eqnarray}
+\sigma_x^{-2} & = & \frac{1}{2}(\sigma_a^{-2} + \sigma_b^{-2}) - \frac{1}{2}(\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
+\sigma_y^{-2} & = & \frac{1}{2}(\sigma_a^{-2} + \sigma_b^{-2}) + \frac{1}{2}(\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
+\sigma_{xy}   & = & \frac{1}{2} \sin (2 \theta) (\sigma_b^{-2} - \sigma_a^{-2})
+\end{eqnarray}
+These formulae thus define the values of $\sigma_x$, $\sigma_y$, and
+$\sigma_{xy}$ given $\sigma_a$, $\sigma_b$, and $\theta$.  Note that
+in this equation and the one above, we represent the quantities in
+terms of $\sigma_b^{-2} - \sigma_a^{-2}$ which is always greater than
+0, thus attributing the sign of the equation to the $\sin$ or $\cos$
+term.  This is necessary to determine the angle in the proper quadrant
+using the arctangent below.
+
+With the above relationships, we may now form combinations that help
+us to solve for $\sigma_a$, $\sigma_b$, and $\theta$:
+\begin{eqnarray}
+f_1 = \sigma_y^{-2} + \sigma_x^{-2} & = & \sigma_b^{-2} + \sigma_a^{-2} \\
+f_2 = \sigma_y^{-2} - \sigma_x^{-2} & = & (\sigma_b^{-2} - \sigma_a^{-2}) \cos (2 \theta) \\
+f_3 = \sqrt{f_2^2 + 4\sigma_{xy}^2} & = & \sigma_b^{-2} - \sigma_a^{-2}
+\end{eqnarray}
+From these, we may derive the equations for $\sigma_a$, $\sigma_b$, and $\theta$:
+\begin{eqnarray}
+\theta & = & \frac{1}{2} \arg (-2 \sigma_{xy}, f_2) \\
+\sigma_a & = & \sqrt{\frac{2}{f_1 - f_3}} \\
+\sigma_b & = & \sqrt{\frac{2}{f_1 + f_3}}
+\end{eqnarray}
+
+The relationship between the rotated ($m_{i,j}$) and unrotated
+($M_{i,j}$) second moments, the latter being equal to $\sigma_a^2$ and
+$\sigma_b^2$, is derived in a similar fashion.  We start with the
+point that the second moment is rotated as a tensor:
+\begin{equation} 
+\left| 
+\begin{array}{cc}
+m_{x,x} & m_{x,y} \\
+m_{y,x} & m_{y,y} \\
+\end{array} \right| 
+=
+\left| 
+\begin{array}{cc}
++\cos \theta & +\sin \theta \\
+-\sin \theta & +\cos \theta \\
+\end{array} \right| 
+\left| 
+\begin{array}{cc}
+M_{x,x} & 0 \\
+0       & M_{y,y} \\
+\end{array} \right| 
+\left| 
+\begin{array}{cc}
++\cos \theta & -\sin \theta \\
++\sin \theta & +\cos \theta \\
+\end{array} \right| 
+\end{equation}
+Multiplying this out and substituting $\sigma_a^2$, $\sigma_b^2$ for $M_{x,x}$, $M_{y,y}$, we find:
+\begin{eqnarray}
+m_{x,x} & = & \sigma_a^{2} \cos^2 \theta + \sigma_b^{2}\sin^2 \theta \\
+m_{y,y} & = & \sigma_b^{2} \cos^2 \theta + \sigma_a^{2}\sin^2 \theta \\
+m_{x,y} & = & -\frac{1}{2} \sin (2 \theta) (\sigma_a^2 - \sigma_b^2)
+\end{eqnarray}
+Using the double-angle relationships, these become:
+\begin{eqnarray}
+m_{x,x} & = & \frac{1}{2}(\sigma_a^{2} + \sigma_b^{2}) + \frac{1}{2}(\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
+m_{y,y} & = & \frac{1}{2}(\sigma_a^{2} + \sigma_b^{2}) - \frac{1}{2}(\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
+m_{x,y} & = & -\frac{1}{2} \sin (2 \theta) (\sigma_a^{-2} - \sigma_b^{-2})
+\end{eqnarray}
+These three formulae define the second moments in terms of $\sigma_a$, $\sigma_b$, and $\theta$. 
+
+We define equivalent intermediate products to the above:
+\begin{eqnarray}
+g_1 = m_{x,x} + m_{y,y}   	 & = & \sigma_a^{2} + \sigma_b^{2} \\
+g_2 = m_{x,x} - m_{y,y}   	 & = & (\sigma_a^{2} - \sigma_b^{2}) \cos (2 \theta) \\
+g_3 = \sqrt{f_2^2 + 4 m_{x,y}^2} & = & \sigma_a^{2} - \sigma_b^{-2}
+\end{eqnarray}
+From these, we may derive the equations for $\sigma_a$, $\sigma_b$, and $\theta$:
+\begin{eqnarray}
+\theta   & = & \frac{1}{2} \arg (-2 m_{x,y}, g_2) \\
+\sigma_a & = & \sqrt{\frac{g_1 - g_3}{2}} \\
+\sigma_b & = & \sqrt{\frac{g_1 + g_3}{2}}
+\end{eqnarray}
 
 \section{PSLib Astronomy Utilities}
@@ -1793,5 +1945,5 @@
 \alpha_p & = & 180^\circ - 192.85948^\circ \\
 \delta_p & = & 90^\circ - 62.87175^\circ \\
-\phi_p & = & 90^\circ + 32.93192^\circ \\
+\phi_p & = & 90^\circ + 32.93192^\circ
 \end{eqnarray}
 
@@ -2768,5 +2920,5 @@
 \begin{eqnarray}
 x_p & = & \rho_x x \\
-y_p & = & \rho_y y \\
+y_p & = & \rho_y y
 \end{eqnarray}
 % 
@@ -2830,5 +2982,5 @@
 \begin{eqnarray}
 \alpha - \alpha_p & = & \arctan (\sin \alpha, \cos \alpha) \\
-\delta            & = & \arcsin (\sin \delta) \\
+\delta            & = & \arcsin (\sin \delta)
 \end{eqnarray}
 %
@@ -2845,5 +2997,5 @@
 y           & = & \frac{-\cos \theta \cos \phi}{\sin \theta} \\
 \sin \theta & = & \zeta / \sqrt{1 + \zeta^2} \\
-\cos \theta & = & 1 / \sqrt{1 + \zeta^2} \\
+\cos \theta & = & 1 / \sqrt{1 + \zeta^2}
 \end{eqnarray}
 
@@ -2860,5 +3012,5 @@
 y           & = & -\cos \theta \cos \phi \\
 \sin \theta & = & \sqrt{1 - R_\theta^2} \\
-\cos \theta & = & R_\theta \\
+\cos \theta & = & R_\theta
 \end{eqnarray}
 
@@ -2922,5 +3074,5 @@
 \begin{eqnarray}
 x & = & \phi \left( 2 \cos \frac{2 \theta}{3} - 1 \right) \\
-y & = & \pi \sin \frac{\theta}{3} \\
+y & = & \pi \sin \frac{\theta}{3}
 \end{eqnarray}
 
@@ -2930,5 +3082,5 @@
 \theta & = & 3 \sin^{-1} \rho \\
 \phi   & = & \frac{x}{1 - 4\rho^2} \\
-{\rm where}\hspace{1cm} \rho & \equiv & y/\pi \\
+{\rm where}\hspace{1cm} \rho & \equiv & y/\pi
 \end{eqnarray}
 
@@ -3031,11 +3183,11 @@
 This function is a two-dimensional Gaussian with an elliptical
 cross-section and a constant local background:
-\[
+\begin{equation}
 f(x,y) = Z_o e^{-z} + S_o
-\]
+\end{equation}
 where
-\[
+\begin{equation}
 z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
-\]
+\end{equation}
 
 Below is the relationship between the \code{psModel} parameters and
@@ -3080,11 +3232,11 @@
 This function is a polynomial approximation of a 2D Gaussian.  The
 function is very similar to the real Gaussian:
-\[
+\begin{equation}
 f(x,y) = Z_o (1 + z + z^2/2 + z^3/6)^{-1} + S_o
-\]
+\end{equation}
 where
-\[
+\begin{equation}
 z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
-\]
+\end{equation}
 
 Below is the relationship between the \code{psModel} parameters and
@@ -3133,11 +3285,11 @@
 than the Taylor series values of 1/2 and 1/6.  The
 function is very similar to the pseudo-Gaussian:
-\[
+\begin{equation}
 f(x,y) = Z_o (1 + z + B_2 (z^2/2 + B_3 z^3/6))^{-1} + S_o
-\]
+\end{equation}
 where
-\[
+\begin{equation}
 z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy}
-\]
+\end{equation}
 
 Below is the relationship between the \code{psModel} parameters and
@@ -3189,12 +3341,12 @@
 core, where the core has a different contour from the wings.  
 
-\[
+\begin{equation}
 f(x,y) = Z_{\rm pk} (1 + z_1 + z_2^M)^{-1} + Sky
-\]
+\end{equation}
 where
-\[
+\begin{equation}
 z_1 = \frac{x^2}{2\sigma_{x,in}^2} + \frac{y^2}{2\sigma_{y,in}^2} + x y \sigma_{xy,in}
 z_2 = \frac{x^2}{2\sigma_{x,out}^2} + \frac{y^2}{2\sigma_{y,out}^2} + x y \sigma_{xy,out}
-\]
+\end{equation}
 
 \begin{verbatim}