
% declare the overall style of the paper.  The \documentstyle lines
% that are relevant for the AASTeX macros are shown; one is uncommented out
% so that the file can be processed.

%\documentstyle[12pt,aasms4]{article}
\documentstyle[11pt,aaspp4]{article}
%\documentstyle[aas2pp4]{article}

% The eqsecnum style changes the way equations are numbered.  Normally,
% equations are just numbered sequentially through the entire paper.
% If eqsecnum appears in the \documentstyle command, equation numbers will
% be sequential through each section, and will be formatted "(sec-eqn)",
% where sec is the current section number and eqn is the number of the
% equation within that section.  The eqsecnum option can be used with
% any substyle.

%\documentstyle[11pt,eqsecnum,aaspp4]{article}

% Authors are permitted to use the fonts provided by the American Mathematical
% Society, if they are available to them on their local system.  These fonts
% are not part of the AASTeX macro package or the regular TeX distribution.

%\documentstyle[12pt,amssym,aasms4]{article}

% Here's some slug-line data.  The receipt and acceptance dates will be 
% filled in by the editorial staff with the appropriate dates.  Rules will 
% appear on the title page of the manuscript until these are uncommented 
% out by the editorial staff.

%\received{4 August 1988}
%\accepted{23 September 1988}
%\journalid{337}{15 January 1989}
%\articleid{11}{14}


% Authors may supply running head information, if they wish to do so, although
% this may be modified by the editorial offices.  The left head contains a
% list of authors, usually three allowed---otherwise use et al.  The right
% head is a modified title of up to roughly 44 characters.  Running heads
% are not printed.

\lefthead{Newberg}
\righthead{SDSS QSO Target Selection}

% This is the end of the "preamble".  Now we wish to start with the
% real material for the paper, which we indicate with \begin{document}.
% Following the \begin{document} command is the front matter for the
% paper, viz., the title, author and address data, the abstract, and
% any keywords or subject headings that are relevant.

\begin{document}

%\title{QSO Target Selection in the SDSS}
\title{The Algorithm to Select Outliers From the Stellar Locus in Three Colors}

\author{Heidi Jo Newberg}
\affil{Fermi National Accelerator Laboratory, Box 500,
    Batavia, IL 60510 \\
    Electronic mail: heidi@fnal.gov}

% Notice that each of these authors has alternate affiliations, which
% are identified by the \altaffilmark after each name.  The actual alternate
% affiliation information is typeset in footnotes at the bottom of the
% first page, and the text itself is specified in \altaffiltext commands.
% There is a separate \altaffiltext for each alternate affiliation
% indicated above.

% The abstract environment prints out the receipt and acceptance dates
% if they are relevant for the journal style.  For the aasms style, they
% will print out as horizontal rules for the editorial staff to type
% on, so long as the author does not include \received and \accepted
% commands.  This should not be done, since \received and \accepted dates
% are not known to the author.

%\begin{abstract}
%Put abstract here.
%\end{abstract}
% The different journals have different requirements for keywords.  The
% keywords.apj file, found on aas.org in the pubs/aastex-misc directory, 
% contains a list of keywords used with the ApJ and Letters.  These are 
% usually assigned by the editor, but authors may include them in their 
% manuscripts if they wish. 
%\keywords{quasars, stars:general, stars:statistics}

% That's it for the front matter.  On to the main body of the paper.
% We'll only put in tutorial remarks at the beginning of each section
% so you can see entire sections together.

% In the first two sections, you should notice the use of the LaTeX \cite
% command to identify citations.  The citations are tied to the
% reference list via symbolic KEYs.  We have chosen the first three
% characters of the first author's name plus the last two numeral of the
% year of publication.  The corresponding reference has a \bibitem
% command in the reference list below.
%
% Please see the AASTeX manual for a more complete discussion on how to make
% \cite-\bibitem work for you.   

%\section{Introduction}
%
%SDSS, other QSO searches, etc.
%
%\section {Color selection of QSOs}
%
%\subsection{Parameterization of the stellar locus}
%
%Review of Newberg and Yanny (1997).
%
\section {Overview}

Given a parameterization of the locus of stars in SDSS colors, we would like
to determine whether an individual catalog entry is consistent with being a
star.  Since there are so many more stars than QSOs in the region containing
the stellar locus, and since the locus of QSOs is diffuse, one cannot recover
much of the stellar locus region of multicolor space by asking whether the
object is consistent with being a QSO.  If it is consistent with being a
star, it is probably a star.

For the individual catalog entry, we require the $u'-g', g'-r',$ and $r'-i'$
colors.  Additionally, we must know the variances and the covariances in each
color.  The target selection algorithm currently assumes that the errors in
each filter are free from systematics, so that $cov(a-b, b-c) = -var(b)$.
The variances for a given color is given by $var(a-b) = var(a) + var(b)$.

The parameterized locus contains a set of ordered locus points, going from
the blue end to the red end.  Each locus point is described by a 
$(u'-g', g'-r', r'-i')$ position, along with an axial unit vector, 
$\hat{k}$, which points to the next locus point; a position angle, 
$\theta$; and a major and minor axis, $a_l, a_m$, which describe
an elliptical cylinder.  Objects which lie within this cylinder, and which
are between the planes which are equidistant from the current point and the
two adjacent points, are considered to be within the stellar locus.  In this
way we build up a tube of elliptical cross sections which can wind around in
multicolor space.

In some cases, an object might not be measured in all filters, resulting
in incomplete knowledge of the position in color space.  For example, if an
object was detected in $g'$, but was below the detection threshold in $u'$,
then we know that $u'-g' > u'_{lim} - g'$.  If the u' measurement was missed
due to a defect in the image, then nothing would be known about the
$u'-g'$ color.  In general, we either know the color and variances of an
object; or we know a limiting color, and the variance is set to a flag
which distinguishes between an upper limit, a lower limit, or no knowledge
of that color.  In this case, the algorithm tries to determine if there is any
allowed value of the unknown colors which places the object within the stellar
locus.  If so, then the object is considered to be consistent with a galactic
star and is not flagged as a candidate.  We require that at least one color 
out of three be known.

The ends of the stellar locus are described by half ellipsoids which fit onto
the ends of the tube.  The center of the ellipsoids are constrained to lie on
the axis of the stellar locus.  A single input parameter determines the position
of each ellipsoid center along this axis ($k_{blue}$ for the blue end and
$k_{red}$ for the red end).  The axes of the ellipsoids are aligned 
with the axes of the elliptical cylinders associated with the first (and last) 
locus points.  Two of the axis widths are given by the major and minor axes of 
the closest locus point interior to the locus.  The width along the locus,
$a_k$, is a separate input parameter for each end.  In the limit that $a_k = 0$,
the locus ends with a plane rather than an ellipsoid.

The input locus parameterization describes the extent of the locus, not including
any photometric errors in the measurements of stars.
An input parameter, $N_{\sigma_{err}}$, determines how many sigma of the errors
in an individual point should be used to ``convolve" with the parameterized
stellar locus.  Before determining whether an object is consistent with being
in the stellar locus, the widths of the locus and the ending ellipsoids are
increased by adding the locus widths with the estimated errors in the axial
directions in quadrature.  This is a correct treatment if the density across
the stellar locus and the error statistics are both Gaussian.
Any limiting colors must have been previously adjusted to include any uncertainty 
in the knowledge of the limit, since variances along an axis with a limit will 
be assumed to be zero.

\section{Finding the closest locus point}

The first step is to determine which locus point is closest to the data point.
We do this by looking at the Euclidean distance in color space, not taking into
account any variation in the width of the locus, to simplify the computation.
In the case that one or two of the colors are not known uniquely, the distance
to each locus point is calculated after first moving the data point as close as
is allowed by the color limits.

If two locus points are equally distant from the datapoint, the redder one is
chosen.  If the closest locus point is not between the centers of the ending
ellipsoids, then the closest locus point is reassigned to be the first one
interior to the ellipsoid centers.

Once the closest locus point is determined, we do not consider any other
locus point.  The convolution is done only on the cylinder associated with this
one locus point, and the object is considered to be consistent with a star if
the color is within that cylinder and interior to the ellipsoids at the ends
of the locus.

\section{Handling errors in individual objects}

First, we estimate the error ellipse in the $(l, m)$ plane.  The variance, covariance
matrix for this plane is given by:
\[V_{ab} = \sum_{x} \sum_{y} \frac{\partial a}{\partial x} \frac{\partial b}{\partial y} S_{xy},\]
where ${\bf V}$ is the (two dimensional) covariance matrix in the $l,m$ plane,
${\bf S}$ is the (three dimensional) covariance matrix in color space, and the
sums are over the three coordinates: $u'-g', g'-r',$ and $r'-i'$.
The error ellipse is given by:
\[{\bf \vec{l}^T V^{-1} \vec{l} } = N_{\sigma_{err}}^2,\]
where $\vec{l} = l\hat{l} + m\hat{m}$.
Solving this, the ellipse parameters are:
\[a_{err}^2 = N_{\sigma_{err}}^2 \frac{V_{ll} + V_{mm} \pm \sqrt{(V_{ll}-V_{mm})^2 + 4 V_{lm}^2}}{2},\]
\[ \tan \theta_{err} = \left\{ \begin{array}{ll}
                0       e V_{lm} = 0, V_{ll}>V_{mm} \\
                \infty  & V_{lm} = 0, V_{mm}>V_{ll} \\
                \frac{-(V_{ll}-V_{mm})+\sqrt{(V_{ll}-V_{mm})^{2}+4V_{lm}V_{lm}}}{2V_{lm}}       & V_{lm} \neq 0
                \end{array}
        \right. \]

Now that we have the parameters for the error ellipse, we will adjust the stellar locus to reflect this information by 'convolving' the error ellipse and the cross-sectional locus ellipse. When convolving the two ellipses, we will make two important assumptions. The error distribution around each data point is gaussian, and the distribution of stars within the stellar locus is gaussian. These assumptions allow us to represent each of the ellipses as a Gaussian function, and, as such, we will convolve them.
If one convolves two bivariate Gaussians, G1 and G2 (with elliptical cross sections), one obtains a
third bivariate Gaussian with elliptical cross section. Here we will convolve the parameterized locus ellipse for the stellar locus in the l-m plane with the rotated error ellipse, also in the l-m plane. By convolving two Gaussians with these ellipses as their cross sections, we will discover a third ellipse which will allow us to determine if the data point is within the locus, or if it is a viable candidate for target selection.

\[G_{convolved} = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{G1(l-l', m-m')G2(l', m')}dl'dm'\]
\[G1 = e^(-(\frac{l^2}{2a_{l}^2} + \frac{m^2}{2a_{m}^2}))\]
\[G2 = e^(-(\frac{(l\cos \theta_{err} + m\sin \theta_{err})^2}{2a_{err maj}^2} + \frac{(-l\sin \theta_{err} + m\cos \theta_{err})^2}{2a_{err min}^2})\]

In this convolution, G1 represents the parametrized locus and G2 represented the error ellipse rotated to the l-m plane.
After the convolution we are left with a Gaussian which is given by:
\[K_{1}K_{2}e^(\frac{-l^2}{2a_{l}^2} -  \frac{m^2}{2a_{m}^2} + \frac{l^2}{4Aa_{l}^4} + \frac{T^2}{4R})\]

Where:
\[A = \frac{ \cos^2 \theta_{err}^2}{2a_{err maj}^2} + \frac{\sin^2\theta_{err}^2}{2a_{err min}^2} + \frac{1}{2a_{l}^2}\]\
\[T = \frac{2l \sin \theta_{err} \cos \theta_{err} (\frac{1}{a_{err maj}^2}-\frac{1}{a_{err min}^2})}{4Aa_{l}^2}\]\
\[R=\frac{\sin^2 \theta_{err}}{2a_{err min}^2} + \frac{\cos^2 \theta_{err}}{2a_{err maj}^2} + \frac{1}{a_{l}^2} + \frac{\cos^2 \theta_{err} \sin^2 \theta_{err} (\frac{1}{a_{err maj}} - \frac{1}{a_{err min}})}{4A}\]\

Here $K_{1}$ and $K_{2}$ are products of the integration, while the various terms within the exponent come from completing the square to facilitate the integration. Since we are only concerned with the elliptical cross section contained in the exponent, we will ignore the rest of the Gaussian and concentrate solely on the exponent. After considerable algebra the convolved ellipse looks like the following:
\[\alpha_{o}l^2-\beta_{o}lm+\gamma_{o}m^2=1\]
Where:
\[\alpha_{o} = \frac{\frac{4R}{4Aa_{l}^4} - \frac{4R}{2a_{l}^2} + \frac{\sin^2\theta_{err}\cos^2\theta_{err}(\frac{1}{a_{err maj}^2} - \frac{1}{a_{err min}^2})^2}{Aa_{l}^2}}{4R}\]\
\[\beta_{o} = \frac{(\frac{\sin\theta_{err}\cos\theta_{err}(\frac{1}{a_{err maj}^2} - \frac{1}{a_{err min}^2})}{4Aa_{l}^2})}{4R}\]\
\[\gamma_{o} = \frac{(\frac{1}{a_{m}^4} - \frac{4R}{2a_{m}^2})}{4R} \]\
 

And, after a great deal of simplification, the convolved ellipse is given by:
\[\alpha l^2 - 2 \beta lm + \gamma m^2 = d,\]
where
\[\begin{array}{l}
d = a_{err maj}^2 a_{err min}^2 + (a_l^2 a_{err maj}^2 + a_m^2 a_{err min}^2) \sin^2 \theta_{err} + (a_{err maj}^2 a_m^2 + a_l^2 a_{err min}^2) \cos^2 \theta_{err} + a_l^2 a_m^2\\
\alpha = a_{err min}^2 \cos^2 \theta_{err} + a_{err maj}^2 \sin^2 \theta_{err} + a_m^2 \\
\beta = \sin \theta_{err} \cos \theta_{err} (a_{err maj}^2 - a_{err min}^2) \\
\gamma = a_{err min}^2 \sin^2 \theta_{err} + a_{err maj}^2 \cos^2 \theta_{err} + a_l^2
\end{array}.\]

The major and minor axes of the convolved ellipse are given by the eigenvalues of the matrix form of the ellipse and the angle of inclination is simply the angle of the eigenvector pointed in the direction of the major axis. 
The variance/covariance matrix for the convolved ellipse is:
\[V_{convolved} = \frac{1}{d}\left( \begin{array}{cc} \alpha & -\beta \\ -\beta & \gamma\end{array} \right)\\
\] 
The characterisitc equation of this matrix is given by
\[\lambda^2 - \frac{\alpha + \gamma}{d}\lambda +\frac{\alpha\gamma - \beta^2}{d^2}\]
Solving this for the eigenvalues we get $\lambda_{1}=\frac{1}{a+b}$ and $\lambda_{2}=\frac{1}{a-b}$. These eigenvalues are the values which diagonalize the matrix, and thus they give the coefficients to any point $(l,m)$ which happens to lie on either the major or minor axis. Thus the major and minor axes are given by:
\[a_{tot maj}^2 = a + b,\]
\[a_{tot min}^2 = a - b,\]
The angle of rotation of the ellipse is found easily from the angle of the eigenvector pointed towards the major axis:
\[\tan \theta_{tot} = \left\{ \begin{array}{ll}
		0			&\beta=0, \alpha=\gamma\\
		\infty			&\beta=0, \alpha\neq\gamma\\
		\frac{\alpha-\gamma+\sqrt{(\alpha-\gamma)^2+4\beta^2}}{2\beta} &\beta  \neq 0, \theta_{err} >= 0\\
		-\frac{\alpha-\gamma+\sqrt{(\alpha-\gamma)^2+4\beta^2}}{2\beta} &\beta \neq 0, \theta_{err} < 0
		\end{array}
	\right. ,\]
where
\begin{eqnarray}
a &=& \frac{(\alpha+\gamma)}{2} = \frac{a_{err maj}^2 + a_{err min}^2 + a_l^2 + a_m^2}{2}, \nonumber\\
b &=& \frac{\sqrt{(\alpha + \gamma)^2 - 4 (\alpha \gamma - \beta^2)}}{2} =  \nonumber\\
  & & \frac{\sqrt{(a_{err maj}^2-a_{err min}^2)^2 + 2 (a_{err maj}^2-a_{err min}^2)(a_l^2-a_m^2)(\cos^2 \theta_{err} - sin^2 \theta_{err}) + (a_l^2-a_m^2)^2}}{2}. \nonumber
\end{eqnarray}
Here, $a$ and $b$ are positive by definition. 
We need not worry about taking the square root of
a negative number when computing the angle, since $a-b \leq \alpha \leq a+b$ for all cases.
If $b=0$, then the resulting ellipse is circular, so we might as well make $\theta=0$.

We now can redefine the $\hat{l}, \hat{m}$ axes associated with the closest locus point
to describe this new elliptical cross section.  The value of $\theta_{tot}$ must be
added to the original angle describing
the locus ellipse, since we have done this with respect to the $\hat{l}$ axis.  When adding,
we must be sure to keep the new value of $\theta$ in the allowed range: $-\frac{\pi}{2} < \theta \leq \frac{\pi}{2}$.

In addition to increasing the width of the excluded locus, we increase $a_k$ on the ends of the locus
in a similar manner:
\[var(k) = \sum_x \sum_y  \frac{\partial k}{\partial x} \frac{\partial k}{\partial y} S_{xy},\]
\[a_{k tot} = \sqrt{a_k^2 + N^2_{\sigma_{err}} var(k) }, \]
where the k in this equation is not for the locus point which is closest to the datapoint, but instead 
for the locus point which is closest to the respective endpoints (subject to being between the two
endpoints).

In the remainder of this paper we will use $a_l, a_m, a_k,$ and $\theta$ to refer
to the derived $a_{tot maj}, a_{tot min}, a_{k tot}$ and $\theta + \theta_{tot}$.  Likewise,
the $\hat{l}, \hat{m}$ unit vectors are in the new coordinate system.

\section{Dealing with non-detections in some filters}

If not all of the colors of the datapoint are known, then we are free to find the allowed point
in color space that is most likely to be in the locus.  This is {\it not} in general the same as
finding the colors that are closest to the locus point in the Euclidean sense, since the
parameterized region is an elliptical cylinder, not a sphere centered at the locus point.
We instead wish to minimize the $r^*$ distance to the locus, where $r^*$ is given by:
\[r^* = \sqrt{\left( \frac{l}{a_l} \right)^2 + \left( \frac{m}{a_m} \right)^2},\]
where $l \equiv \Delta \vec{r} \cdot \hat{l}, m \equiv \Delta \vec{r} \cdot \hat{m}$,
$\Delta \vec{r} \equiv \vec{r} - \vec{r}_p$, and $\vec{r}_p$ is the position in color space
of the closest point.

First, let's solve the case where only the $u'-g'$ color, is unknown.  To make the equations easier
to read, I will use the notation $\Delta \vec{r} \equiv (\Delta x, \Delta y, \Delta z), l = l_x \hat{x}
+l_y \hat{y} + l_z \hat{z}, m = m_x \hat{x} + m_y \hat{y} + m_z \hat{z}$.  The solution is:
\[\Delta x = \left\{ \begin{array}{ll}
	-\frac{\{(a_m^2 l_x l_y + a_l^2 m_x m_y) \Delta y + (a_m^2 l_x l_z + a_l^2 m_x m_z) \Delta z \}}{a_m^2 l_x^2 + a_l^2 m_x^2}	& a_m^2 l_x^2 + a_l^2 m_x^2 \neq 0\\
	-\frac{l_y \Delta y + l_z \Delta z}{l_x}	& a_m^2 l_x^2 + a_l^2 m_x^2 = 0, l_x \neq 0\\
	-\frac{m_y \Delta y + m_z \Delta z}{m_x}	& a_m^2 l_x^2 + a_l^2 m_x^2 = 0, m_x \neq 0\\
	0	& a_m^2 l_x^2 + a_l^2 m_x^2 = 0, l_x = 0, m_x = 0
	\end{array}
	\right. .\]
The first case above is the result of the minimization process.  If the denominator is zero, but
$l_x \neq 0$, then it follows that $a_m = 0$.  Also, either $a_l = 0$ or $m_x = 0$.  In
the first case, we will only be able to find a color point in the locus if we can fortuitously set
$m=l=0$ by moving along the $x$ axis.  In the second case, we cannot affect the magnitude of
m.  Either case is optimized by setting $l = 0$.  If the denominator is zero, but instead
$l_x = 0$, then we still need either $a_l = 0$ or $m_x = 0$.  The first case we cannot affect
the distance along the major axis, so we set the distance along the minor axis to zero.  The second 
case gives us $\hat{x} = \hat{k}$, so we might as well set $\Delta x = 0.$

If the value of $x$ is completely unknown, then we use this calculated $\Delta x$ to determine whether
the datapoint is consistent with being a star.  If the value of the color $x$ is a limit, then we
must ask whether the computed value is within the limits.  If it is not, then we instead use the
limiting $x$ value to determine whether the datapoint is consistent with being a star.

Next, we tackle the case where both $u'-g'$ and $g'-r'$ are unknown.  When two of the colors are unknown,
one can find values that will put the point exactly on the line $\vec{r}-\vec{r}_p = \lambda \hat{k},$
where $\lambda$ is a free parameter.  We can solve for $\lambda$, and then the two unknown colors
using:
\begin{eqnarray}
\lambda &=& \left\{ \begin{array}{ll}
	\frac{\Delta z}{k_z},		& k_z \neq 0\\
	0,				& k_z = 0,\\
	\end{array}
	\right. \nonumber\\ 
\Delta y &=& \lambda k_y, \nonumber\\
\Delta x &=& \lambda k_x. \nonumber
\end{eqnarray}
If $k_z=0$, then the line down the center of the locus is in the $x, y$ plane, so there are many
values of $x$ and $y$ which will be on it.  We arbitrarily choose the one on the locus point,
since we know that this value is also within the $k$ limits of the stellar locus.

The existance of limiting values in $x$ and $y$ make this a little trickier.  First, we
calculate the value of $\Delta y$ which places the point exactly on the center of the locus.
In the case that this value is not consistent with the $y$ limits, we assign 
$\Delta y = y_{lim} - y_p$.  Since we now have only one missing color, we can use the
procedure outlined above (the case where only the $u'-g'$ color is unknown) to calculate
the optimal value of $\Delta x$.  If a limit was not reached, one can verify that the
computed $\Delta x$ will be the same as if we had used $\Delta x = \lambda k_x$.  If
this computed value of $\Delta x$ within the allowed limits, then we have done the best we can.
If it is not, then we set $\Delta x = x_{lim} - x_p$ and then compute the optimal value
of $\Delta y$ given $\Delta x$ and $\Delta z$.  The optimal value of $\Delta y$, given
$\Delta x$ and $\Delta z$, can be computed by reassigning the axes ($x \rightarrow y, y \rightarrow
z, z \rightarrow x$) in the equation that optimizes $\Delta x$ given $\Delta y$ and $\Delta z$.

The above equations and their permutations ($x,y,z \rightarrow y,z,x$ and $x,y,z \rightarrow z,x,y$)
are used to determine the values of $\Delta x$, $\Delta y$, and $\Delta z$ which are most likely
to be included in the stellar locus, given the input limits on these quantities.  Once we
have this coordinate in three dimensions, we can ask whether it is within the parameterized
locus.  We must have $r^* \leq 1$, where
\[r^* = \left\{ \begin{array}{ll}
	\sqrt{\left( \frac{l}{a_l} \right)^2 + \left( \frac{m}{a_m} \right)^2}	& a_l > 0, a_m > 0\\
	\frac{l}{a_l}					& a_l > 0, a_m=0, m = 0\\
	0									& l = m = 0\\
	\infty									& a_l = 0, l \neq 0\\
	\infty									& a_m = 0, m \neq 0
	\end{array}
	\right.  \]
The point must also be within the ellipsoids at the ends of the locus.  If $\Delta k <0$,
then we already know we are within the red end limit, since the locus point is guaranteed
to be within the locus.  Similarly, if $\Delta k >0$, we are guaranteed to be within the
blue end limit.  So, we need only check one of the ellipsoids.  We ignore the fact that the
locus may curve around in $x, y, z$ coordinates, and only look at the $l, m, k$ coordinates.
We are within the ends of the locus if $k_{blue\_limit}(l,m) < k < k_{red\_limit}(l,m)$,
where:
\[k_{blue\_limit} = k_{blue} - a_{k\ blue} \sqrt{1 - {r^*_{blue}}^2} \]
\[k_{red\_limit} = k_{red} + a_{k\ red} \sqrt{1 - {r^*_{red}}^2}. \]
Here, the values of $r^*$ are computed using $a_l$ and $a_m$ for the end ellipses, but
values of $l$ and $m$ calculated from the closest locus point.

\section{End conditions}

If the previously calculated optimal values of $\Delta x, \Delta y,$ and $\Delta z$ put us
within the locus, then we are done.  If we could not find values for the coordinates
that made $r^* \leq 1$, then we are done.  However, if the values failed only the end
conditions, then there is a chance we could still put the datapoint in the locus by
minimizing the distance to the center of the end ellipsoid:
\[r^{**} = \sqrt{\left( \frac{l'}{a'_l} \right)^2 + \left( \frac{m'}{a'_m} \right)^2 + \left( \frac{k' - k'_{end}}{a_k} \right)^2}\]
rather than the center of the cylinder associated with the closest locus point.  Here, $k_{end}$ is either
$k_{blue}$ or $k_{red}$, depending on whether $\Delta k<0$ or $\Delta k>0$.  The primes indicate that
the values are measured with respect to the position in color space of the locus point which is
closest to the center of the end ellipsoid.  The $a_l$ and $a_m$ widths are for the end ellipsoid, not
the cylinder associated with the closest locus point.

We now calculate for the datapoint a new optimal color which minimizes $r^{**}$.  If only $x$
is unknown, we compute $\Delta z'$ and $\Delta y'$ using the definition $\Delta \vec{r}' \equiv 
\vec{r} - \vec{r}_e \equiv (\Delta x', \Delta y', \Delta z'),$ and $\vec{r}_e$ is the position in
color space of the locus point closest to the center of the end ellipsoid (subject to
$k_{blue} < k_e < k_{red}$).  The value of $\Delta x'$ which minimizes $r^{**}$ is
\[\Delta x' = \left\{ \begin{array}{ll}
	- \frac{{a'}_m^2 a_k^2 (l'_y \Delta y' + l'_z \Delta z') l'_x + {a'}_l^2 a_k^2 (m'_y \Delta y' + m'_z \Delta z') m'_x + {a'}_l^2 {a'}_m^2 (k'_y \Delta y' + k'_z \Delta z' - k'_{end}) k'_x}{{a'}_m^2 a_k^2 {l'}_x^2 + {a'}_l^2 a_k^2 {m'}_x^2 + {a'}_l^2 {a'}_m^2 {k'}_x^2}
										& D_1 \neq 0 \\
	- \frac{a_k^2 (l'_y \Delta y' + l'_z \Delta z') l'_x + {a'}_l^2 (k'_y \Delta y' + k'_z \Delta z' - k'_{end}) k'_x}{a_k^2 {l'}_x^2 + {a'}_l^2 {k'}_x^2}
										& D_1 = 0, {a'}_m^2 = 0, m'_x = 0, {a'}_l^2 \neq 0 \\
	-\frac{m'_y \Delta y' + m'_z \Delta z'}{m'_x}				& D_1 = 0, {a'}_m^2 = 0, m'_x \neq 0, {a'}_l^2 \neq 0 \\
	- \frac{a_m^2 (l'_y \Delta y' + l'_z \Delta z') l'_x + {a'}_l^2 (m'_y \Delta y' + m'_z \Delta z') m'_x}{a_m^2 {l'}_x^2 + {a'}_l^2 {m'}_x^2}
										& D_1 = 0, {a'}_m^2 > 0 \\
	0									& D_1 = 0, {a'}_l^2 = 0, m'_x = 0, k'_x \neq 1 \\
	0									& D_1 = 0, l'_x = 1
	\end{array}
	\right. \]
where $D_1 \equiv {a'}_m^2 a_k^2 {l'}_x^2 + {a'}_l^2 a_k^2 {m'}_x^2 + {a'}_l^2 {a'}_m^2 {k'}_x^2$.
The first case is the formal result of the minimization.  The second case is for $a_m^2 = 0$ and $m'_x =0$;
 in this case we cannot affect the $m$ value, so we might as well find the best position within the
$l,k$ ellipse.  If instead we have $m'_x \neq 0$, we instead put the datapoint on the $m=0$ plane.
If $D_1=0$ but ${a'}_m^2>0$, then we must have $a_k^2=kx=0$, since ${a'}_l^2>={a'}_m^2$.  If ${a'}_l^2 = 0$,
then the data point is either on the line of the locus or not, no change in x will give different results than the previous calculated optimal $x$.  All cases where $D_1=0$ and $l'_x=1$ reduce to $\Delta x' = 0$, and this sweeps up all of the cases not covered by the other criteria.  Again, if this computed value for $\Delta x'$ violates a limit in $x$, then we reassign $x$ to the limit.

If both $x$ and $y$ are unknown, then we start by calculating the values which minimize
$r^{**}$.  We deal with limits in the identical way as we did when we were minimizing
$r^*$.  That is, we figure out the optimal value of $x$ and $y$ together, but only assign
the $y$ value.  Then we figure out the optimal value of $x$ given that $y$ value and the
fixed value of $z$.  This way we can deal with the limits in a sensible way.  The optimal 
value of $y$ is given by:
\[\Delta y' = \left\{ \begin{array}{ll}
	\frac{(k'_y k'_z a_k^2 + m'_y m'_z {a'}_m^2 + l'_y l'_z {a'}_l^2) \Delta z' - (l'_x m'_z {a'}_m^2 - m'_x l'_z {a'}_l^2) k'_{end}}{a_k^2 {k'}_z^2 + {a'}_m^2 {m'}_z^2 + {a'}_l^2 {l'}_z^2}
		& D_2 \equiv a_k^2 {k'}_z^2 + {a'}_m^2 {m'}_z^2 + {a'}_l^2 {l'}_z^2 \neq 0 \\
	k'_y k'_{end}					& D_2 = 0
	\end{array}
	\right. .\]
The first option is the formal minimization of $r^{**}$.  If $k'_z = 0$, then the center of
the ellipsoid is in the $x,y$ plane.  Therefore, we try to place the datapoint on that center
(or as close as we can get, given that $\Delta z'$ might not be zero).
If the denominator is zero, then
${a'}_m^2=0$.  This is true because we cannot have $k'_z=1, l'_z=0, m'_z=0$ due to the
definition of our system.  So, either ${a'}_l^2$ or ${a'}_m^2$ must be zero.  Either way,
${a'}_m^2=0$.  If $k'_z \neq 0$, we must have ${a'}_k^2 = 0$ as well.  Additionally,
either ${a'}_l^2= 0$ (the endpoint is the only hope) or $l'_z=0$ ($\hat{l}$ is in the
$x,y$ plane, so the endpoint is still the best choice, if we can get there).

As before, we then adjust $\Delta y'$ if it is outside the allowed limits, then use our
equations for only $\Delta x'$ missing to find the optimal $x'$.  If the computed
value of $x'$ is outside the limits, then we set $x'$ to the limit, and recompute
the optimal value of $y'$.

I cannot think of a case in which the new position is in the interior side of the
endplane at $k_{end}$.  So, all I have to do is figure out if the new position 
satisfies ${r'}^* \leq 1$ and is interior to the $k$ limits on the outer surface
of the ellipsoid.  If it is, then it is consistent with being a star, otherwise
it is not.

%\section{Discussion and Results}

%\acknowledgments

%\clearpage


%\begin{deluxetable}{rrrrrrrrrrrr}
%\footnotesize
%\tablecaption{Parameterization of the locus of $7^{th}$ magnitude stars in WBVR \label{table1}}
%\tablewidth{0pt}
%\tablehead{
%\colhead{} & \colhead{$k$} & \colhead{$N_p$} & \colhead{$(W-B)_p$} & \colhead{$(B-V)_p$} & \colhead{$(V-R)_p$} & \colhead{$k_{(W-B)_p}$} & \colhead{$k_{(B-V)_p}$} & \colhead{$k_{(V-R)_p}$} & \colhead{$a_p$} & \colhead{$b_p$} & \colhead{$\theta_p$}
%}
%\startdata
%1& 0.000& 383& $-$0.469& $-$0.073& $-$0.034& 0.991& 0.124& 0.060& 0.038& 0.012& $-$0.828\nl
%2& 0.268& 427& $-$0.203& $-$0.040& $-$0.018& 0.982& 0.171& 0.085& 0.037& 0.012& $-$0.799\nl
%3& 0.484& 659& 0.005& 0.010& 0.008& 0.876& 0.399& 0.271& 0.026& 0.013& $-$0.937\nl
%4& 0.643& 891& 0.111& 0.104& 0.080& 0.313& 0.732& 0.606& 0.043& 0.013& $-$1.408\nl
%5& 0.802& 779& 0.097& 0.223& 0.184& $-$0.332& 0.696& 0.636& 0.055& 0.013& 1.270\nl
%6& 0.967& 610& 0.007& 0.322& 0.280& $-$0.624& 0.574& 0.530& 0.056& 0.010& 1.102\nl
%7& 1.149& 734& $-$0.119& 0.421& 0.367& $-$0.366& 0.728& 0.579& 0.048& 0.009& 1.352\nl
%8& 1.294& 610& $-$0.102& 0.538& 0.451& 0.409& 0.755& 0.513& 0.053& 0.009& $-$1.194\nl
%9& 1.473& 264& 0.009& 0.656& 0.527& 0.736& 0.563& 0.377& 0.040& 0.009& $-$0.908\nl
%10& 1.728& 166& 0.215& 0.781& 0.613& 0.789& 0.521& 0.327& 0.056& 0.019& $-$1.006\nl
%11& 1.926& 239& 0.366& 0.892& 0.675& 0.807& 0.525& 0.271& 0.034& 0.010& $-$0.852\nl
%12& 2.052& 328& 0.475& 0.950& 0.701& 0.879& 0.433& 0.201& 0.027& 0.010& $-$0.731\nl
%13& 2.152& 365& 0.564& 0.990& 0.720& 0.884& 0.410& 0.225& 0.027& 0.011& $-$0.710\nl
%14& 2.246& 372& 0.646& 1.029& 0.744& 0.868& 0.421& 0.263& 0.029& 0.010& $-$0.653\nl
%15& 2.337& 329& 0.725& 1.068& 0.769& 0.887& 0.397& 0.234& 0.028& 0.011& $-$0.738\nl
%16& 2.447& 273& 0.825& 1.109& 0.792& 0.896& 0.384& 0.223& 0.029& 0.010& $-$0.741\nl
%17& 2.545& 275& 0.912& 1.148& 0.815& 0.865& 0.417& 0.278& 0.033& 0.013& $-$0.804\nl
%18& 2.661& 254& 1.009& 1.198& 0.851& 0.865& 0.411& 0.288& 0.041& 0.014& $-$0.758\nl
%19& 2.786& 240& 1.120& 1.247& 0.885& 0.864& 0.414& 0.288& 0.047& 0.014& $-$0.729\nl
%20& 2.940& 228& 1.250& 1.314& 0.931& 0.851& 0.419& 0.316& 0.052& 0.014& $-$0.735\nl
%21& 3.108& 179& 1.393& 1.382& 0.987& 0.836& 0.419& 0.355& 0.031& 0.013& $-$0.758\nl
%22& 3.278& 194& 1.533& 1.455& 1.051& 0.818& 0.429& 0.384& 0.034& 0.011& $-$0.845\nl
%23& 3.459& 194& 1.680& 1.532& 1.121& 0.794& 0.409& 0.450& 0.027& 0.015& $-$1.181\nl
%24& 3.646& 207& 1.824& 1.605& 1.216& 0.643& 0.396& 0.656& 0.027& 0.019& $-$1.051\tablebreak
%25& 3.827& 217& 1.912& 1.675& 1.358& 0.187& 0.253& 0.949& 0.052& 0.019& $-$0.178\nl
%26& 4.008& 147& 1.887& 1.691& 1.537& $-$0.367& $-$0.036& 0.930& 0.096& 0.021& 0.581\nl
%27& 4.233& 95& 1.767& 1.661& 1.725& $-$0.593& $-$0.142& 0.792& 0.080& 0.026& 0.396\nl
%28& 4.519& 62& 1.585& 1.618& 1.940& $-$0.639& $-$0.149& 0.755& 0.069& 0.030& 0.320
%\enddata
%\tablecomments{
%This table contains one line for each of the 28 locus points along the stellar
%locus.  For each locus point we record: column 1 ($k$) - the distance in
%magnitudes along the stellar locus, starting from the first locus point;
%column 2 ($N_p$) - the number of sources associated with the locus point
%(not including the outliers which were rejected); columns 3, 4, 5 ($(W-B)_p, (B-V)_p,
%(V-R)_p$) - the $W-B, B-V, V-R$ position of the locus point; columns 6, 7, 8
%($k_{(W-B)_p}, k_{(B-V)_p}, k_{(V-R)_p}$) - the components of the unit vector
%$\hat{k_p}$ along the locus; and columns 9, 10, 11 ($a_p, b_p, \theta_p$) -
%the major axis, minor axis, and position angle of the ellipse fit to
%the cross section of the stellar locus.
%}
%\end{deluxetable}

%\clearpage
% Each reference has a \bibitem command to define the citation format
% to be placed in the text (in []) and the symbolic tag used for 
% cross referencing (in {}).
%
% See sample1.tex, or the AASTeX guide, for an alternative to the \cite-
% \bibitem command.

%\begin{thebibliography}{}
%\end{thebibliography}

%\clearpage

%\figcaption{ Algorithm Demonstration.  We show here the progress of the
%algorithm with each iteration.  The algorithm starts out (0 iterations)
%with a locus point at the position of each endpoint.  During the first
%iteration, a new locus point is added at the midpoint between the first
%two locus points, ellipses are fit to the cross section at each of the
%three locus points, and then the points are moved to the centroids of their 
%associated stars.  Since
%the ellipse fits were calculated before the locus points were moved, the
%calculated ellipse widths are large.  Therefore, on the second iteration
%the locus points are closer together than 3 times the larger axes of the
%ellipses, so no new locus points are added.  The standard deviations are now
%re-calculated, and the three locus points are again moved to the centroid
%of their associated stars.  Note that the set of associated stars for each point
%in the second iteration is different from the set of associated stars in
%the first iteration.  During each iteration up to iteration 100, the same
%three steps: add points, calculate standard deviations, move points - are repeated.  The
%last plot shows the positions of the locus points after outlier rejection
%has been run.  Because the outliers of this distribution are fairly
%symmetric, the locus points hardly moved.  However, outlier rejection
%decreased the size of the major axes of the ellipse fits at each locus point 
%by 17\% on average.
%\label{fig1}}

%\clearpage

%\plotone{slpfig1.ps}

\end{document}
