%\documentclass[panstarrs,psreport,12pt]{panstarrs}
\documentclass[panstarrs]{panstarrs}

% basic document variables
\title{Pan-STARRS Moving Object Processing System}
\subtitle{Trailed Objected Fitting}
\author{Joseph Masiero, Eugene Magnier, Paul Price, Robert Jedicke}
\shorttitle{Trailed Objects}
\group{Pan-STARRS MOPS development group}
\project{Pan-STARRS Moving Object Processing System}
\organization{Institute for Astronomy}
% \distribution{}
\version{01}
\docnumber{PSDC-xxx-xxx}
% note the use of the docnumber & version number:
% the complete PSDC document number is given by
% \thedocnumber-\theversion

\begin{document}
\maketitle

% -- Revision History --
\RevisionsStart
% version   Date         Description
01        & 2006.10.04 & First draft \\
\RevisionsEnd

% -- TBD / TBR Listing --
%\TBDsStart
% section % page Date        Description
%1.2.3.4   & 2    & xxx-xxx & missing a dingle on a dangle \\
%\TBDsEnd

\pagebreak
\tableofcontents

\pagebreak
\listoffigures

\pagebreak 
\pagenumbering{arabic}
\section{Overview}

Many of the Near-Earth Objects (NEOs) that will be imaged by Pan-STARRS will
have significant ($>1~$pixel) apparent motions over the $30~$second exposure
time of each image.  These objects, instead of being point sources, will be
imaged as trails moving across the CCD.  While these trails provide more
information as to the rate of motion of the object, they also present a
problem in extracting that information in an automated manner using the
PSF-fitting code that is being developed by the Pan-STARRS Image Processing
Pipeline (IPP).  In order to accurately determine orbits of NEOs that are
trailed, we need to determine their centers, length (and thereby rate of
motion), angle (and thereby their direction of motion to within an
uncertainty of a negative sign), and brightness (and thereby their magnitude
and size, assuming an albedo and calculating an orbit).  To do this, we have
developed an analytical description of a trailed asteroid, which the IPP's
{\tt psphot} code will use as its model for fitting non-point sources in the
Pan-STARRS difference images.  Because these images will have had the static
sky removed, the only this remaining will be point-source transients,
non-point source transients, cosmic rays, and other noise.  Once cosmic ray
removal routines have been run (or transients not seen in multiple images
removed, in the case of PS4), and PSF fitters run on the images, only trailed
asteroids and comets will remain.  Comet fitting routines are discussed in
PSDC-xxx-xxx (Kleyna, et al), and so we shall only concern ourselves with
trailed asteroids.

In order to test the effectiveness and efficiency of our trailed object
fitting routine, we have developed a synthetic trail generator, which can
create artificial trails at a variety of magnitudes, lengths, angles, and
positions.  We then applied our fitting model to these synthetic trails.
Results of these tests are described, as well as problems with the model they
have illuminated that need to be addressed.


\section{Synthetic Trail Generation}
We generated our synthetic trailed objects by inserting a large number of
individual ``sub-PSFs'' along the line of motion at $0.2~$pixel spacing to
build up a trail.  The ``sub-PSFs'' spacing size was chosen to ensure a
smooth final trail without bogging down the generation code due to excessive
numbers of iterations.  In principle, anything with motions greater than
$\sim0.5~$pixels per image should be detectable however in practice it is
difficult to distinguish between motions of this magnitude and slight
variations in object PSF.

Each ``sub-PSF'' is a Gaussian ellipsoid, where the major axis was matched to
the assumed-circular PSF measured for stars in the field, and the minor axis
had a value of $90\%$ that of the major axis.  This was done to model real
sources which may not necessarily have perfectly circular PSFs.  The
``sub-PSFs'' were then co-added to create a smooth trail.  Poissonian noise
was added to the counts on each pixel, and the trail was inserted into the
data image at a random location.  Each trail was given a randomly
distributed:

\begin{itemize}
\item location of trail center
\item 'stationary' $g'$ magnitude between $13$ and $20$~mag
\item angle of motion (from the pixel x-axis)
\item rate of motion between ($0~-~10$~deg/day)
\end{itemize}

Trail positions were generated without regard to chip gaps, bad columns,
stars, or other image features.  Fig~\ref{fig.strail} shows an example
of one of our trails, which had a center location of x$~=~108.3~$pixels,
y$~=~1270.2~$pixels, 'stationary' magnitude of $18.96~$mag, angle of motion of
$29.62^\circ$, and rate of motion of $6.04~$deg/day.


\section{Analytical Model}

We developed an analytical model to fit both synthetic and real trailed
objects with an LM $\chi^2$ minimizer (Magnier et al, 2006).  Although the most
appropriate model is a convolution of a PSF with a line of the desired length
and orientation, this method is incredibly computationally intensive, making
the run-time prohibitive, and finding analytical derivatives of such a
function that are easily implemented in a program quickly exceeds the scope
of this project.

We have chosen to represent our trailed objects with a piecewise function
consisting of a 1-D Gaussian ``pipe'' with 2-D Gaussian ``endcaps'', as
illustrated by a wire diagram in Fig~\ref{fig.wirestrail}.  First we define a new
coordinate system $\left(x',~y'\right)$, where $x'$ is along the trailed
objects line of motion, and $y'$ is perpendicular to the line of motion, and
the location $\left(0,~0\right)$ is the center of the trail.  Beyond a
half-length in both the $+x'$ and $-x'$ directions, we follow the fitted
stellar PSF exactly.  Between these two point we use a constant brightness
line with Gaussian falloff along the $y'$ axis, matching the FWHM of the
fitted stellar PSF in that direction.

To use the LM $\chi^2$ minimizer, we employed the following equations for the
model:


\begin{eqnarray}
      F_{(x_1,y_1)} & = & F_{peak}~\left(1+z+\frac{z^2}{2}+\frac{z^3}{6}\right)^{-1} + F_{sky}\\
      z & = & \frac{x_1^2}{2~\sigma_x^2} + \frac{y_1^2}{2~\sigma_y^2} + \sigma_{xy}~x_1~y_1
\end{eqnarray}
where $F$ is the flux at a given location; $F_{peak}$ is the peak flux;
$F_{sky}$ is the sky background flux; $\sigma_x$, $\sigma_y$, and
$\sigma_{xy}$ are the ellipsoid PSF shape parameters in the unrotated frame;
and $x_1$ and $y_1$ are the image pixel coordinates shifted so that (0,0) is
the center of the trail (i.e. $x1 = x - x0$ and $y1 = y - y0$), which can be
expressed as (when $l$ is the fitted length and $\theta$ is the fitted angle
with respect to the pixel x-axis):

\begin{itemize}
\item for $x' > \frac{l}{2}$:
\begin{eqnarray}
      x_1 & = & \left(x'-\frac{l}{2}\right)~\cos\theta - y'~\sin\theta\\
      y_1 & = & \left(x'-\frac{l}{2}\right)~\sin\theta + y'~\cos\theta
\end{eqnarray}
\item for $x' < -\frac{l}{2}$:
\begin{eqnarray}
      x_1 & = & \left(x'+\frac{l}{2}\right)~\cos\theta - y'~\sin\theta\\
      y_1 & = & \left(x'+\frac{l}{2}\right)~\sin\theta + y'~\cos\theta
\end{eqnarray}
\item for $\frac{l}{2} \ge x' \ge  -\frac{l}{2}$:
\begin{eqnarray}
      x_1 & = & -y'~\sin\theta\\
      y_1 & = & y'~\cos\theta
\end{eqnarray}
\end{itemize}

The LM minimizer also requires the derivative of the model with respect to
all variables.  The derivatives are:

\begin{eqnarray}
	\frac{dF}{dF_{sky}} & = & 1\\
        \frac{dF}{dF_{peak}} & = & \left(1+z+\frac{z^2}{2}+\frac{z^3}{6}\right)^{-1}
\end{eqnarray}

and letting $q = F_{peak}~\left(1+z+\frac{z^2}{2}+\frac{z^3}{6}\right)^{-2}~\left(1+z+\frac{z^2}{2}\right)$:

\begin{itemize}
\item for $x' > \frac{l}{2}$:
\begin{eqnarray}
	  \frac{dF}{dx_0} & = & q~\left(S_x^2~x_1 + S_{xy}~y_1\right)\\
	  \frac{dF}{dy_0} & = & q~\left(S_y^2~y_1 + S_{xy}~x_1\right)\\
	  \frac{dF}{dS_x} & = & -q~S_x~x_1^2\\
	  \frac{dF}{dS_y} & = & -q~S_y~y_1^2\\
	  \frac{dF}{dS_{xy}} & = & -q~x_1~y_1\\
	  \frac{dF}{dl} & = & q~\left(S_x^2~x_1~\frac{\cos\theta}{2} + S_y^2~y_1~\frac{\sin\theta}{2} + S_{xy}~x_1~\frac{\sin\theta}{.}+ S_{xy}~y_1~\frac{\cos\theta}{2}\right)\\
	  \frac{dF}{d\theta} & = & -q~\left(S_x^2~x_1~\frac{l}{2}~\sin\theta - S_y^2~y_1~\frac{l}{2}~\cos\theta - S_{xy}~x_1~\frac{l}{2}~\cos\theta + S_{xy}~y_1~\frac{l}{2}~\sin\theta\right)
\end{eqnarray}
\item for $x' < -\frac{l}{2}$:
\begin{eqnarray}
	  \frac{dF}{dx_0} & = & q~\left(S_x^2~x_1 + S_{xy}~y_1\right)\\
	  \frac{dF}{dy_0} & = & q~\left(S_y^2~y_1 + S_{xy}~x_1\right)\\
	  \frac{dF}{dS_x} & = & -q~S_x~x_1^2\\
	  \frac{dF}{dS_y} & = & -q~S_y~y_1^2\\
	  \frac{dF}{dS_{xy}} & = & -q~x_1~y_1\\
	  \frac{dF}{dl} & = & -q~\left(S_x^2~x_1~\frac{\cos\theta}{2} + S_y^2~y_1~\frac{\sin\theta}{2} + S_{xy}~x_1~\frac{\sin\theta}{2} + S_{xy}~y_1~\frac{\cos\theta}{2}\right)\\
	  \frac{dF}{d\theta} & = & -q~\left(-S_x^2~x_1~\frac{l}{2}~\sin\theta + S_y^2~y_1~\frac{l}{2}~\cos\theta + S_{xy}~x_1~\frac{l}{2}~\cos\theta - S_{xy}~y_1~\frac{l}{2}~\sin\theta\right)
\end{eqnarray}
\item for $\frac{l}{2} \ge x' \ge  -\frac{l}{2}$:
\begin{eqnarray}
	  \frac{dF}{dx_0} & = & q~\sin\theta~\left(S_x^2~\sin\theta~x_1 - S_y^2~y_1~\cos\theta - S_{xy}~x_1~\cos\theta + S_{xy}~y_1~\sin\theta\right)\\
	  \frac{dF}{dy_0} & = & q~\cos\theta~\left(-S_x^2~x_1~\sin\theta + S_y^2~\cos\theta~y_1 + S_{xy}~x_1~\cos\theta - S_{xy}~y_1~\sin\theta\right)\\
	  \frac{dF}{dS_x} & = & -q~S_x~x_1^2\\
	  \frac{dF}{dS_y} & = & -q~S_y~y_1^2\\
	  \frac{dF}{dS_{xy}} & = & -q~x_1~y_1\\
	  \frac{dF}{dl} & = & 0\\
	  \frac{dF}{d\theta} & = & -q~(S_x^2~x_1~\left(x'~\sin\theta - y'~\cos\theta\right) + S_y^2~y_1~\left(x'~\cos\theta - y'~\sin\theta\right) + \nonumber \\
	& & S_{xy}~x_1~\left(x'~\cos\theta - y'~\sin\theta\right) + S_{xy}~y_1~\left(x'~\sin\theta - y'~\cos\theta\right))
\end{eqnarray}
\end{itemize}

The variables of the equations were: $F_{peak}$, the
pixel-coordinate location of the trail center: $x_0$ and $y_0$, the
length $l$, and the angle with respect to the pixel x-axis $\theta$.

The $S_x = \frac{1}{\sigma_x}$, $S_y = \frac{1}{\sigma_y}$, and $S_{xy} =
\sigma_{xy}$ values will be fitted by the earlier PSF-fitting stage of the
image processing and passed to our code.  The IPP finds median values for the
PSF parameters from all the stellar-like objects in the image, and allows
these values to vary as a function of position in the field.  $F_{sky}$ is
likewise calculated during an earlier stage of image processing.

We have also developed a model that will co-fit a circular PSF to the trail
in addition to the other parameters, however when trying to co-fit the trail
parameters and the ellipsoid PSF parameters, we find that the number of
variables causes a steep degradation in fit quality.


\section{Trail Identification Efficiency and Fitting Accuracy}
To test the efficiency of our trail-fitting routines, we first inserted our
synthetic trails randomly into Canada-France-Hawaii Telescope (CFHT)/Megacam
data.  We then ran the IPP source determination routine {\tt psphot} on the
images, first fitting all PSF-like objects with a 2-dimensional third-order
approximation to a Gaussian.  All non-PSF objects in the field were assumed
to be trails and the above model was fitted to them.  Because there were
galaxies in our image, not everything that was passed to our routine was a
trail, and fitted center pixels were compared to generated center pixels to
recover trails from the data.  Fitted trails within 2 pixels of a generated
trail in the x and/or y directions were considered ``confirmed'' fits, while fitted
trailed within 10 pixels in x and/or y were considered ``potential'' fits.

Our overall fitting efficiency, comparing the number of generated trails to
the number of fitted trails (both ``confirmed'' and ``potential'', can be seen in
Figs~\ref{fig.maghist}\&\ref{fig.mageff}.  We find that our model is applied
to $\ge50\%$ of the generated trails down to a 'stationary' magnitude of
$g=17$, with a significant decrease at larger magnitudes.  We find the source
low efficiency measurement to be the code that passes objects to our
trail-fitting routine.  It seems that the synthetic trailed asteroids are
different enough from expected object shapes that the core {\tt psphot} code
cannot decide whether they are PSFs, extended sources, or bad objects
(e.g. cosmic rays, diffraction spikes), and so punts on some of them.

For those trails that were ``confirmed'' fits, we determine the RMS of the
scatter in pixel location to be $0.20~$pixels from the center of the
generated trail to the center of the fitted trail.  A scatter plot of the
difference in location centers is shown in Fig~\ref{fig.xy}, with the dotted
line describing the cut between the black ``confirmed'' points and the red
``potential'' points.  This color distinction will hold for the other plots
used to compared fitted and generated trails.  Looking at the comparison
between generated and fitted trail length in Fig~\ref{fig.length}, there are
two features that are cause for concern.  The first is the $\sim3.9~$pixel
offset between the apparent linear relation between the generated and fitted
lengths and the one-to-one relation for these same values.  The second is the
distinction between red and black points, especially at longer trail lengths,
with the red points being far more poorly fit than the black.  The first
feature shows a relation with the assumed PSF of the generated trails, where
a smaller PSF causes a smaller offset between the two lines.  Ideally this
can be modeled out, and a process to do this is being investigated.  After
artificially removing the offset in the generated-to-fitted comparison we
find an RMS scatter of $1.9~$pixels in the ``confirmed'' black points.  The
second feature is caused by trails that are only partially fit, e.g. those
broken by stars, bad columns, or noise.  This will show offsets in length and
pixel center, but should show the correct fitted angle, as shown below.  This
problem is much more difficult to deal with, as it involves developing a
smarter initial object determinator.

Looking at the relation between the fitted and generated angles of motion,
we see a much tighter correlation, as shown in Fig~\ref{fig.theta}.  There
are some points randomly scattered, but this is unrelated to whether the
fitted trails were labeled as ``confirmed'' or ``potential'', and in fact
many of the red ``potential'' points have reasonably accurate angular fits,
as predicted above.  Much of this scatter can be explained by trails with
very short lengths, and thus very ambiguous angles.  As trails become longer,
angle fitting accuracy improves, though some scatter does still remain.  We
find an RMS scatter in angle of motion of $6.8^\circ$ for the ``confirmed''
trail fits when looking across all trail lengths.  For trails with motions
$5\le~v~\le10~$deg/day, see Fig~\ref{fig.thetalong}, we find an RMS scatter
of $1.4^\circ$ for the same points.  There are a few points with very large
errors ($\sim180^\circ$) near the $-90^\circ$ and $+90^\circ$ extrema, but
this is to be expected as fits near these values are nearly identical.  There
also appears to be a large increase in the scatter of the fits from
one-to-one near angles of $0^\circ$ and $90^\circ$.  This effect is currently
unexplained, but may be a result of the trigonometric functions used in the
angle determination, and is currently being investigated.  

When comparing fitted magnitude to the generated magnitude, we first needed
to shift the fitted magnitude by the $26.28~$mag zero-point, as given in the
image header.  This value was calculated during the image processing step by
the MegaCam Elixir system.  A comparison of these two values can be seen in
Fig~\ref{fig.mag}.  There is significant scatter, but no seeming pattern,
other than that the ``potential'' trails show much dimmer fitted magnitudes
than the ``confirmed'' ones, as would be expected if they were only partial
fits of the trail length.  However, when we compare the difference in
generated and fitted magnitudes as a function of length, as seen in
Fig~\ref{fig.mloff}, we see an interesting, and worrying trend.  As the
generated trail length gets longer, the fitter calculates bigger (dimmer)
magnitudes for the objects than are actually generated.  Additionally, as the
trail length goes to zero, the fitted magnitudes are actually brighter than
the generated ones.  The latter problem may be due to incorrect trail
generation techniques, or could be due to the same effect that is causing the
underestimation at longer trail lengths, which is currently unexplained and
being investigated.

\section{Summary}

We have developed an analytical model of a trailed asteroid, to be
incorporated into the Pan-STARRS IPP.  This model will be fit to trailed
objects in the Pan-STARRS images to extract information about the position
and motion of nearby objects, which will aid in the orbit identification and
recovery of the object.  Application of our analytical trailed asteroid model
to synthetic trails has confirmed that our model and fitting routines are
reasonably accurate, as well as illuminated areas of that need to be
improved.

In particular, before this fitting algorithm can be released into the main
pipeline, a few issues need to be resolved, including:
\begin{itemize}
\item The PSF dependent offset between the generated and fitted lengths
\item The fitting of parts of trails instead of the full trail
\item The increase in error of the determined angle near $0^\circ$ and $90^\circ$
\item The difference between generated and fitted magnitudes, especially as a function of length
\item The passing of only $\sim50\%$ of the observed trailed objects to the fitting routine
\end{itemize}

Despite the areas that require investigation, the accuracy of the trailed
objects that are fit indicates that our model and this method of fitting will
be effective for determining the location and motion of NEOs and other
quickly moving objects.

\bibliographystyle{plain}
\begin{center}
{\bf Bibliography}
\end{center}
Bertin, E. \& Arnouts, S., 1996, A\&AS, 117, 393.

Press, W.~H., Teukolsky, S.~A., Vetterling, W.~T. \& Flannery, B.~P., 2002, Numerical Recipes, Cambridge University Press.

Magnier, E.~A., Price, P.~A., Lupton, R.~H. \& Hoblitt, J., 2006, Pan-STARRS PS-1 IPP Library Supplementary Design Requirements Specification, PSDC-430-007, http://project.pan-starrs.ifa.hawaii.edu/project/PSDC/PSDC-4xx/PSDC-430-012-10.pdf.

%\bibitem [Masiero et~al.(2006)]{SPD} 
%Masiero, J.~R., Jedicke, R., Price, P.~A., Kubica, J., \& Heasley, J.~N., 2006, in prep.


%\bibliography{panstarrs}


\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{trail.jpg} 
\end{center}
\protect\caption{An example of a synthetic trail generated by our code.  In this example, the pixel x-axis runs to the right of the figure, and pixel y-axis runs toward the top.  }
\label{fig.strail}
\end{figure}
\clearpage

\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=1]{wirestrail.jpg}  
\end{center}
\protect\caption{A wire diagram of the analytical model used to fit trailed objects.  The axes are labeled in the rotated system aligned with the object's motion.
}
\label{fig.wire}
\end{figure}
\clearpage

\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailHistmag.jpg}  
\end{center}
\protect\caption{Histogram of magnitudes of trails generated (red) and found (black).
}
\label{fig.maghist}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailHistmageff.jpg}  
\end{center}
\protect\caption{Efficiency (fitted/generated).
}
\label{fig.mageff}
\end{figure}
\clearpage


\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailCompxy.jpg}  
\end{center}
\protect\caption{Differences in x and y positions between the generated and fitted synthetic trails.  The red points are defined as fitted trails that were more than 2 pixels in x or y away from the generated trail center (anything falling outside the dotted-line box), and will be defined as such in the remaining plots.
}
\label{fig.xy}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailComplength.jpg}  
\end{center}
\protect\caption{Generated trail length compared to fitted length.  Note the linear correlation which is offset from a 1:1 relation by a constant $\sim3~$pixels.  Also note that the red points tend to be fit by too short a trail length.
}
\label{fig.length}
\end{figure}
\clearpage

\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailComptheta.jpg}  
\end{center}
\protect\caption{Generated theta (with respect to the pixel x axis) compared to fitted theta.  Note the lack of correlation between red points and goodness of fit.  Also note the apparent increase in scatter near 0 and 90 degrees.
}
\label{fig.theta}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailLongtheta.jpg}  
\end{center}
\protect\caption{Generated theta (with respect to the pixel x axis) compared to fitted theta for trails with lengths between $5$ and $10$ degrees per day.
}
\label{fig.thetalong}
\end{figure}
\clearpage

\newpage
\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailCompmag.jpg}  
\end{center}
\protect\caption{Generated magnitude compared to fitted magnitude, where the fitted magnitude has been shifted by $28.26$ magnitudes, which is the magnitude zero-point of the images used, as given by the MegaCam Elixir image processing system.
}
\label{fig.mag}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{strailCompmloff.jpg}  
\end{center}
\protect\caption{The change in magnitude (generated - fitted) compared to the generated length of the trail.  A negative value on the y-axis implies an object fitted to a fainter flux level than was generated.
}
\label{fig.mloff}
\end{figure}
\clearpage


\end{document}

