%\documentclass[preprint,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}

% Some other (several out of many) possibilities
\documentclass[preprint,aps]{revtex4}
%\documentclass[preprint,aps,draft]{revtex4}
%\documentclass[prb]{revtex4}% Physical Review B

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math

%\nofiles

\begin{document}

\preprint{APS}

\title{Determination of the Orbit of Mars\\ Using Kepler's Triangulation Technique}  

\author{N. Copernicus}
\affiliation{University of Central Arkansas}
\email{copernicus@uca.edu}


\date{\today}

\begin{abstract}


The eccentricity ($\epsilon$), semi-major axis ($a$), and perihelion position ($\omega$) for the orbit of Mars was determined 
using the same triangulation method developed by J. Kepler.
Data analyzed include T. Brahe's observations from the 1500's and Harvard Observatory data from
the 1931 through 1950.
The best fit orbital parameter values found were $\epsilon=0.0985 \pm 0.0034$, $a=1.529\pm0.004$ AU, and $\omega=326^\circ  \pm  2^\circ$.
The percent differences between these values and current modern values were $\%\Delta \epsilon=5\%$, 
$\%\Delta a=0.3\%$, and $\%\Delta\omega=3\%$.



\end{abstract}

\maketitle

\section{\label{sec:level1}INTRODUCTION}


Around 1600 Johannes Kepler developed a geometrical technique to determine 
the position of Mars in it's orbit from observations from the Earth \cite{kepler}.
If Mars can be observed from two different positions when it is at a particular point
in its orbit, then one can triangulate the location of Mars (see Fig. 1).    
Since both the Earth and Mars are revolving around the Sun, one cannot pick random
observation dates to do the triangulation, since Mars will likely not be in the same position
in its orbit on those dates.
By 1600 Kepler knew that the orbital period of Mars was approximately 687 days 
(current accepted value is about 686.980 days).
This means Mars will be at the same location in its orbit every 687 days, but the Earth 
will in a different location in its orbit, given the orbital period of the Earth is 365.256 days.
Therefore, one can triangulate the position of Mars by analyzing the apparent position of 
Mars on the sky for dates separated by some integer multiple integer of 687 days.
Given that Kepler had access to the Martian position data collected by Tycho Brahe \cite{kepler},
 which covered a couple of decades, he was able to carry out this triangulation 
 technique for several positions of Mars, and thus trace out the orbit of Mars.
 This lead to Kepler's first law of planetary motion, which states that the orbits of planets
 in general are elliptical with the Sun located at one of the foci of the ellipse.  
 
 Later Isaac Newton was able to derive from first principles of motion and gravity the
 trajectory of objects moving under the influence of a gravitational force\cite{newton}.
The position in polar coordinates is given by:
\begin{equation}
r(\theta)={{\alpha}\over{1+\epsilon cos(\theta-\omega)}}
\end{equation}
In Eq. (1) $r(\theta)$ is the distance from the focus, $\epsilon$ is the eccentricity, 2$\alpha$ is a 
parameter termed the latus rectum which along with $\epsilon$ gives the semi-major axis ($a$) of the orbit,
$\theta$ is the angular position, and $\omega$ is the angular position of the perihelion point of the
orbit.

In this project we analyzed some of the original Brahe data that Kepler also used, along with Martian
apparent position data from Harvard photographic taken in the 1930's and 1940's, to determine $\epsilon,
\alpha, \omega$, and $a$ for the orbit of Mars.

 
 
\section{PROCEDURES}

\subsection{Data Used}

The data needed for the triangulation method are the heliocentric longitude of the Earth and simultaneous geocentric
longitude of Mars.
The zero point of the longitudes is the direction to the vernal equinox, and the longitude increases in the counterclockwise 
direction as view from above the north pole of the solar system.
In addition their needs to be measurements separated by 687 days.

Table I shows five pairs of dates and the longitudinal angles that are from Tycho Brahe's observations that Kepler \cite{kepler} published
and used in his original analysis.
The measurements are given to the nearest arcminute, therefore we assumed that the uncertainty on the geocentric longitude of Mars
measurements to be $\pm$1 arcminute.

Table II shows the measurements off Harvard College Observatory photographic plates on eight pairs of dates\cite{harvard} between
1931 and 1950.
The photographic plates where taken with a wide angle camera with a field of view of about $35^\circ$.
Superimposing a grid geocentric coordinates along the ecliptic plane on the plates, the position of Mars was measured with an uncertainty of $\pm1^\circ$.


\subsection{Triangulation}

We made a few simplifying assumptions before doing the analysis.  
First, we assumed that the orbit of the Earth is perfectly circular, in fact it has an eccentricity of about 0.0167.
Second, we assumed that the orbits of Earth and Mars are exactly co-planer, but the orbits are inclined by 1.85$^\circ$.
For error propagation, it was was assumed that the position of the Earth was known exactly and that the geocentric
longitudes of Mars were the only values with uncertainties.  
The orbital radius of the Earth is, by definition, 1 AU, about $1.496\times10^{8}$km, and was used as the unit of length
in the analysis.

The orbital radius of the Earth, the heliocentric longitude of the Earth ($\theta$), and the geocentric longitude of Mars ($\phi$) 
on a particular date were used to determine the equation of a line that passed through Mars and Earth on that date.
The slope of the line ($m$) was found from $\phi$:
\begin{equation}
m=tan(\phi)
\end{equation}
A point on that line is the position of the Earth.
\begin{equation}
(cos\theta,sin\theta)
\end{equation}
Eq. (3) gives the position of the Earth in cartesian coordinates.
Putting Eqs. (2) and (3) together we found the y-intercept for the line.
The resulting equation for the line is given by:
\begin{equation}
y(x)=mx+(sin\theta-mcos\theta)
\end{equation}

For the pairs of dates separated by 687 days we were able to find using Eq. (4) the equations of the two
lines whose intersection was the position of Mars.
An example of the intersecting lines is shown in Fig. 1.
We solved the system of the two linear equations for their intersection giving us the position of Mars
in cartesian coordinates ($x_{Mars},y_{Mars}$).
\begin{equation}
x_{Mars}={{(sin\theta_2 -sin\theta_1)+(m_1 cos\theta_1-m_2 cos\theta_2)}\over{m_1-m_2}}
\end{equation}
\begin{equation}
y_{Mars}=m_1 x_{Mars}+(sin\theta_1 -m_1 cos\theta_1)
\end{equation}
The subscripts 1 and 2 in Eqs. (5) and (6) indicate the different dates of the pair of observations.
Converting to heliocentric polar coordinates the position of Mars is given by:
\begin{equation}
r_{Mars}= \sqrt{x_{Mars}^2+y_{Mars}^2}
\end{equation}

\begin{equation}
\theta_{Mars}= arctan({y_{Mars}\over{x_{Mars}}})
\end{equation}


\subsection{Orbit Fit}

Once several positions of Mars were determined, then Eq. (1) was fit to those positions.
The fitting was done with the plotting and fitting software {\tt gnuplot}.
It uses an implementation of the nonlinear least-squares Marquart-Levenberg algorithm \cite{marquardt, gnuplot} to
fit functions.
The quality of the fit is given by the reduced chi-square ($\chi^2_\nu$), which is the sum-of-squares
of the residuals between the fit and the data divided by the number of degrees of freedom ($\nu$).
In this case, there are three adjustable parameters ($\alpha$, $\epsilon$, and $\omega$) and therefore $\nu$=3.
The asymptotic standard errors \cite{gnuplot} of the parameters are also given by {\tt gnuplot}.
Strictly speaking, these are not the true one standard deviation confidence levels for the parameters,
but give a qualitative but over-optimistic estimate of the parameter uncertainties.



\section{RESULTS}

Fig. 1 shows an example of the triangulation using the first pair of observations from Table I.  
The error bars shown in the inset also shows the x and y error for the position of Mars which
results from the uncertainty in the geocentric longitudes.  
Table III lists all thirteen Mars positions in polar coordinates with the resulting uncertainties that
result when propagating the geocentric longitude errors.

Fig. 2 shows the Mars positions and the fitted ellipse.
The $\chi^2_\nu$ of the fit is 0.00014.
The best fit parameter values are $\alpha=1.514 \pm 0.004$ AU,
$\epsilon=0.0985 \pm 0.0034$, and $\omega=326^\circ  \pm  2^\circ$.
Given that the semi-major axis $a$ is related to $\alpha$ and $\epsilon$ by
\begin{equation}
a={ {\alpha}\over{1-\epsilon^2}}
\end{equation}
and gives $a=1.529\pm0.004$ AU.

 

\section{DISCUSSION}

A single agreed upon set of values for the orbital parameters of the planets does not exist.
This is partly due to the fact that those parameters are not fundamental physical constants
like the speed of light or planck's constant that presumably do not change over time.
Another reason for this is that continued research is still refining the orbital parameters for the planets.
One repository for these parameters that is being maintained with the latest values is the Planetary Fact Sheet
by D. R. Williams \cite{williams}
The orbital parameters for Mars listed there are  $\epsilon_{Mars}=0.09341233$, $a_{Mars}=1.52366231$, $\omega_{Mars}$=336.04084.
Taking these to be the current accepted values, the percent difference between our values and these are $\%\Delta \epsilon=5\%$, 
$\%\Delta a=0.3\%$, and $\%\Delta\omega=3\%$.

The current value of $\epsilon_{Mars}$ is smaller than the value we derived and outside of the optimistic uncertainty.
This means the fit resulted in a more elliptical shape than the true value.

The current value of  $a_{Mars}$ is larger than the value we derived and is also outside the optimistic uncertainty.
This means the fitted orbit is a bit smaller than the accepted orbit.

Finally the orientation of the orbit relative to the vernal equinox $\omega_{Mars}$ is larger than the value we derived from the fit.
This means that the fitted orbit has a perihelion point closer to the vernal equinox than the true orbit.

The results could be improved to a certain degree by using more Mars points than the thirteen used in this analysis.
Thirteen data points may have given slightly better results if they were more uniformly distributed around the orbit.
One can see from Fig. 2 that seven of the points are clustered between $\theta=140^\circ$ and $190^\circ$ where as the remaining
six points are outside of that 50$^\circ$ section.
A higher order correction would have been to include the eccentricity of the Earth's orbit, but this would only have been worth doing
if a greater number of more uniformly distributed points were used for the analysis.
A still higher order correction would be to make the corrections for the inclination of the Martian orbit, but again, not worth doing given
this particular data set.

The triangulation technique ingeniously developed by Kepler is no longer the method used to determine the orbits for objects in the solar system.
The downside to Kepler's technique is having to know the orbital period of the object, and then obtaining observations separated by
some integer number of orbital periods.
This is impractical for determining the orbits of asteroids and comets where none of those things may be possible.
Today two or three observations at random intervals are all that are necessary.
Given accurate date and time of observations are known and accurate and precise equatorial coordinates are measured,
the the orbit can be determined with a minimum of two or three observations  using the numerical multi-point boundary-value problem and
quasilinearization method \cite{bellman}.
This numerical technique forms the basis for the software developed to determine the orbits solar system bodies.

\section{CONCLUSIONS}

We have determined the basic orbital parameters of Mars using the same technique and data that Kepler
employed to establish his first law of planetary motion, supplemented with Harvard data from the 1930's and 1940's.
We were able to reproduce Kepler's results that the orbit of Mars is elliptical with the Sun at one of the foci of the ellipse, 
and we were able to get within a few percent or less the values for the eccentricity (5\%), semi-major axis length (0.3\%), and perihelion point position (3\%).  
The true values of these parameters are not within the optimistic-error bars for the parameter values, but a probably close to or
within the actual standard deviation of the parameter values that would result from a Monte Carlo mapping of parameter space.

In summary, we have shown that Kepler's method does indeed work, and can get within a few percent of the current values using
only a few data points.
But, we did not achieve the same accuracy and precision as the current modern values.


\begin{acknowledgments}
I wish to thank the UCA Department of Physics and Astronomy for
the use of their facilities to carry out this investigation.
\end{acknowledgments}



\begin{thebibliography}{00}

\bibitem{kepler} J. Kepler, {\it Astronomia Nova} (1609).
\bibitem{newton} I. Newton, {\it Philosophi� Naturalis Principia Mathematica } (1687).
\bibitem{harvard} \url{http://homepage.mac.com/dtrapp/ePhysics.f/labII_7.html}.
\bibitem{marquardt} D. W. Marquardt, Journal of the Society for Industrial and Applied Mathematics {\bf11}, 431, (1963).
\bibitem{gnuplot} \url{http://www.gnuplot.info/docs/gnuplot.pdf}.
\bibitem{williams} \url{http://nssdc.gsfc.nasa.gov/planetary/planetfact.html}.
\bibitem{bellman} R. Bellman, H. Kagiwada, and R. Kalaba, in {\it Procedings of the National Academy of Sciences of
the United States of America} {\bf 48}, 1327, (1962).
\end{thebibliography}

\pagebreak

FIG. 1. Shows an example of the resulting triangulation for the first pair of observations in Table 1.
  The inset shows the error bars on the Mars position assuming a one arc minute uncertainty on the geocentric longitude.

FIG. 2. Shows the triangulated positions of Mars and the elliptical fit.  The radial grid is spaced at 0.5 AU intervals from the Sun, 
and the angular intervals are 30$^\circ$.
The major axis of the fitted ellipse is also shown.

\pagebreak

\begin{table}
\caption{Tycho's observations of Mars from Kepler's Astronomia nova.}
\begin{ruledtabular}
\begin{tabular}{lcc}
Date & $\theta$ (Heliocentric Long. of Earth) & $\phi$ (Geocentric Long. of Mars) \\
\hline
1585 Feb. 17 & 159$^\circ$ 23$'$ & 135$^\circ$ 12$'$\\
1587 Jan. 5 & 115$^\circ$ 21$'$ & 182$^\circ$ 08$'$\\
\  &  \  &  \  \\
1591 Sep. 19 & 5$^\circ$ 47$'$ & 284$^\circ$ 18$'$\\
1583 Aug. 6 & 323$^\circ$ 26$'$ & 346$^\circ$ 56$'$\\
\  &  \  &  \  \\
1593 Dec. 7 & 85$^\circ$ 53$'$ & 3$^\circ$ 04$'$\\
1595 Oct. 25 & 41$^\circ$ 42$'$ & 49$^\circ$ 42$'$\\
\  &  \  &  \  \\
1587 Mar. 28 & 196$^\circ$ 50$'$ & 168$^\circ$ 12$'$\\
1589 Feb. 12 & 153$^\circ$ 42$'$ & 218$^\circ$ 48$'$\\
\  &  \  &  \  \\
1585 Mar. 10 & 179$^\circ$ 41$'$ & 131$^\circ$ 48$'$\\
1587 Jan. 26 & 136$^\circ$ 06$'$ & 184$^\circ$ 42$'$\\
\end{tabular}
\end{ruledtabular}
\end{table}

\pagebreak

\begin{table}
\caption{Harvard Mars observations.}
\begin{ruledtabular}
\begin{tabular}{lcc}
Date & $\theta$ (Heliocentric Long. of Earth) & $\phi$ (Geocentric Long. of Mars )\\
\hline
1931 Mar. 21 & 180$^\circ$ & 118$^\circ$ \\
1933 Feb. 5 & 136$^\circ$ & 168$^\circ$\\
\  &  \  &  \  \\
1933 Apr. 20 & 210$^\circ$ & 151$^\circ$ \\
1935 Mar. 8 & 167$^\circ$ & 204$^\circ$\\
\  &  \  &  \  \\
1935 May 26 & 245$^\circ$ & 187$^\circ$ \\
1937 Apr. 12 & 202$^\circ$ & 246$^\circ$\\
\  &  \  &  \  \\
1939 Sep. 26 & 3$^\circ$ & 301$^\circ$ \\
1941 Aug. 13 & 320$^\circ$ & 20$^\circ$\\
\  &  \  &  \  \\
1941 Nov. 22 & 60$^\circ$ & 12$^\circ$ \\
1943 Oct. 10 & 17$^\circ$ & 80$^\circ$\\
\  &  \  &  \  \\
1944 Jan. 21 & 121$^\circ$ & 66$^\circ$ \\
1945 Dec. 8 & 76$^\circ$ & 123$^\circ$\\
\  &  \  &  \  \\
1946 Mar. 19 & 178$^\circ$ & 108$^\circ$ \\
1948 Feb. 4 & 135$^\circ$ & 153$^\circ$\\
\  &  \  &  \  \\
1948 Apr. 4 & 195$^\circ$ & 138$^\circ$ \\
1950 Feb. 20 & 152$^\circ$ & 191$^\circ$\\
\end{tabular}
\end{ruledtabular}
\end{table}


\pagebreak


\begin{table}
\caption{Mars positions.}
\begin{ruledtabular}
\begin{tabular}{cc}
r (AU) (Sun-Mars Distance)& $\theta (^\circ)$ (Heliocentric Longitude) \\
\hline
1.691$\pm$0.002 & 149.22$\pm$0.07\\
1.379$\pm$0.011 & 149.22$\pm$0.80\\
1.503$\pm$0.001 & 149.22$\pm$0.03\\
1.639$\pm$0.001 & 149.22$\pm$0.02\\
1.675$\pm$0.002 & 149.22$\pm$0.08\\
1.68$\pm$0.10 & 149.6$\pm$5.3\\
1.64$\pm$0.04 & 182.5$\pm$1.0\\
1.57$\pm$0.04 & 219.7$\pm$1.3\\
1.38$\pm$0.05 & 341.0$\pm$4.5\\
1.46$\pm$0.08 & 42.5$\pm$2.9\\
1.63$\pm$0.01 & 96.3$\pm$0.4\\
1.67$\pm$0.29 & 142.3$\pm$2.5\\
1.65$\pm$0.04 & 168.6$\pm$1.9\\

\end{tabular}
\end{ruledtabular}
\end{table}




\end{document}

