% BB Feb 3 2010.

%\documentstyle[prd,aps,epsfig,floats]{revtex}
\documentclass[twocolumn,showpacs,amsmath,amssymb,prd,floats,floatfix]{revtex4}
%\documentstyle[prd,aps,epsfig,twocolumn,floats]{revtex}
%\documentstyle[prd,aps,multicol,epsfig]{revtex}
%\renewcommand{\narrowtext}{\begin{multicols}{2}\global\columnwidth20.5pc}
%\renewcommand{\widetext}{\end{multicols}\global\columnwidth42.5pc}

%\usepackage{amssymb,amsmath}
\usepackage{epsfig}

\begin{document}
%\wideabs{
\title{ Kepler Problem -- PHY 4241 Notes }

\author{Bernd A. Berg$^{\rm \,a}$}

\affiliation{
$^{\rm \,a)}$ Department of Physics, Florida State University,
  Tallahassee, FL 32306-4350 }
% (E-mail: berg@csit.fsu.edu)\\

% \date{Jan 7, 2007; revised May 9, 2007}
\date{\today }

\begin{abstract}
Kepler problem as discussed in class.
\end{abstract}

\maketitle

\section{Solution of the Kepler Problem for Numerical Implementation}

The purpose of this section is to cast the analytical solution of the 
Kepler problem into a form, which allows for easy implementation into 
a computer program (here done in Fortran). The final solutions are 
translated back to the initially given inertial system. There are 
many treatments of the Kepler problem, for instance \cite{LL,MT}. 
We follow to some extent by Landau and Lifschiz. 

\subsection{Central potential problem}

Let $t$ be the time and the masses be $m_1$ and $m_2$. In a frame 
$\Sigma$ defined by Cartesian orthonormal vectors $\hat{e}_1,\,
\hat{e}_2,\,\hat{e}_3$ positions and velocities of the masses 
$m_i$, ($i=1,2$) are given by:
\begin{equation} \label{sigma}
  \vec{r}_i(t) = \left(\begin{matrix} x_i^1(t)\cr x_i^2(t)\cr x_i^3(t)
  \end{matrix} \right)~~{\rm and}~~ \vec{v}_i(t) =
  \left(\begin{matrix} \dot{x}_i^1(t)\cr \dot{x}_i^2(t)\cr 
  \dot{x}_i^3(t) \end{matrix} \right)
\end{equation}
where the dot denotes the time derivative. For positions and 
velocities the initial conditions at time $t_0$ are
\begin{equation} \label{initial}
  \vec{r}_i(t_0)~~{\rm and}~~\vec{v}_i(t_0)\,,~~i=1,2
\end{equation}
in the frame $\Sigma$. The computer program will calculate 
$\vec{r}_i(t)$ and $\vec{v}_i(t)$, $i=1,2$ for any desired time $t$.

We denote the total mass by
\begin{equation} \label{M}
  M=m_1+m_2\,,
\end{equation}
the Center of Mass (CM) position and velocity are defined by
\begin{eqnarray} \label{rcm} \vec{r}_{\rm cm}(t)
  &=& \frac{m_1\,\vec{r}_1(t)+m_2\,\vec{r}_2(t)}{M}\,,
\\  \label{vcm} \vec{v}_{\rm cm}(t)
  &=&\frac{m_1\,\vec{v}_1(t)+m_2\,\vec{v}_2(t)}{M}\,.
\end{eqnarray}
Due to momentum conservation one finds 
\begin{eqnarray} \label{vCM}
  \vec{v}_{\rm cm}(t) &=& \vec{v}_{\rm cm}(t_0)~~{\rm and}\\ 
  \vec{r}_{\rm cm}(t) &=& \vec{r}_{\rm cm}(t_0) + \label{rCM}
  \vec{v}_{\rm cm}(t_0)\,\triangle t\,,~~\triangle t = (t-t_0)\,.
\end{eqnarray}
Defining the difference coordinates and velocities by
\begin{eqnarray} \label{vec12} 
  \vec{r}_{12}(t) &=& \vec{r}_1(t)-\vec{r}_2(t)\,,
  \\  \label{v12} 
  \vec{v}_{12}(t) &=& \vec{v}_1(t)-\vec{v}_2(t)\,,
\end{eqnarray}
the positions of the particles are
\begin{eqnarray} \label{r1} \vec{r}_1(t) 
  &=& \vec{r}_{\rm cm}(t) + \frac{m_2\,\vec{r}_{12}(t)}{M}\,,
  \\ \label{r2} \vec{r}_2(t) 
  &=& \vec{r}_{\rm cm}(t) - \frac{m_1\,\vec{r}_{12}(t)}{M}\,.
\end{eqnarray}
These equations show that the CM is located in between the particles 
on the straight line connecting them. Equations for the velocities 
are defined by taking the time derivatives.

So, the task has become to solve the equations of motion for the 
difference coordinates and to implement the solution numerically.
In the following we used the notation
\begin{equation} \label{r12}
  r_{12}=|\vec{r}_{12}(t)|\ .
\end{equation}
For a central potential $U(r_{12})$ the energy 
\begin{equation} \label{E}
  E=\frac{m_1\,\vec{v}_1(t)^2}{2}+\frac{m_2\,\vec{v}_2(t)^2}{2}
   + U(r_{12})
\end{equation}
is conserved, i.e., does not depend on the time $t$.

As the difference coordinate $\vec{r}_{12}(t)$ stays invariant under
transformation to the CM frame $\Sigma'$ defined by
\begin{equation} \label{CMframe}
  \vec{r}_{\rm cm}^{\,'}(t)=0\,. 
\end{equation}
We can calculate its time dependence in the CM system and still use 
Eq.~(\ref{r1}) and~(\ref{r2}) to find the time dependence of
the original coordinates. In the CM systems these equations
simplify to
\begin{equation} \label{r12cm}
  \vec{r}_1^{\,'}(t) =   \frac{m_2\,\vec{r}^{\,'}_{12}}{M}~~{\rm and}
  ~~\vec{r}^{\,'}_2(t) = - \frac{m_1\,\vec{r}^{\,'}_{12}}{M}\,.
\end{equation}
In the following we drop the primes and continue to work in the CM
frame (we shall use the primes for yet another frame soon).

Inserting the derivatives of equation (\ref{r12cm}) (without primes)
into the energy conservation (\ref{E}), we find
\begin{equation} \label{Ecm}
  E_{\rm cm} = \frac{\mu\,\vec{v}_{12}(t)^2}{2} + U(r_{12})
  ~~{\rm where}~~\mu = \frac{m_1\,m_2}{M}\,.
\end{equation}
$\mu$ is called reduced mass. To simplify the notation, we drop the 
subscript $12$: 
\begin{equation} \label{r}
  \vec{r}(t)=\vec{r}_{12}(t)~~{\rm and}~~\vec{v}(t)=\vec{v}_{12}(t)\,. 
\end{equation}
With $\vec{r}(t_0)$ and $\vec{v}(t_0)$ given, we want to find 
$\vec{r}(t)$ for the cental potential problem (\ref{Ecm}). Besides
momentum and energy the angular momentum is conserved:
\begin{equation} \label{Loriginal}
  \vec{L} = \vec{r}_1\times\vec{p}_1 + \vec{r}_2\times\vec{p}_2
\end{equation}
where $\vec{p}_i=m_i\vec{v}$, $i=1,2$ are the momenta of the masses. 
Using (\ref{r12cm}) we have
\begin{equation} \label{Lmu}
  \vec{L} = \mu\,\vec{r}(t)\times\vec{v}(t) =
            \mu\,\vec{r}(t_0)\times\vec{v}(t_0)\,.
\end{equation}
We deal with $\vec{L}=0$ later.
Assuming $\vec{L}\ne 0$, the motion takes place in the plane spanned by 
$\vec{r}(t_0)$ and $\vec{v}(t_0)$ and it is convenient to describe to
describe in in a new coordinate frame $\Sigma'$ defined by Cartesian
orthonormal vectors $\hat{e}'_1,\,\hat{e}'_2,\,\hat{e}'_3$ with
\begin{eqnarray} \label{ep1}
  \hat{e}'_1 &=& \frac{\vec{r}(t_0)}{|\vec{r}(t_0)|}\,,\\ \label{ep2}
  \hat{e}'_2 &=& \frac{\vec{v}(t_0)-[\vec{v}(t_0)\cdot\hat{e}'_1]\,
                 \hat{e}'_1}{|\vec{v}(t_0)-[\vec{v}(t_0)\cdot
                 \hat{e}'_1]\,\hat{e}'_1|} \,,\\ \label{ep3}
  \hat{e}'_3 &=& \hat{L}=\frac{\vec{L}}{L}~~{\rm with}~~L=|\vec{L}|\,.
\end{eqnarray}
By definition of the angular momentum this is a right-handed frame. 
In this frame the absolute value of the angular momentum reads
\begin{equation} \label{L}
  L = \mu\,r(t)^2\,\dot{\phi}^{\,'}(t)\,,
\end{equation}
where the azimuth angle $\phi'$ is defined by
\begin{eqnarray} \label{xp}
  x'(t) &=& r(t)\,\cos[\phi'(t)]\,, \\ \label{yp}
  y'(t) &=& r(t)\,\sin[\phi'(t)]\,.
\end{eqnarray}
Defining the effective potential by
\begin{equation} \label{Ueff}
  U_{\rm eff}(r) = U(r) + \frac{L^2}{2\mu\,r(t)^2}\,,
\end{equation}
where $L^2/\,[2\mu r(t)^2]$ is called centrifugal energy, the energy 
conservation (\ref{Ecm}) becomes
\begin{equation} \label{Eeff}
  E_{\rm cm} = \frac{\mu}{2}\,\dot{r}(t)^2 + U_{\rm eff}(r)\,.
\end{equation}
Separation of variables yields 
\begin{eqnarray} \nonumber
  dt &=& dr\,\sqrt{\frac{\mu}{2\left[E_{\rm cm}-U_{\rm eff}(r)\right]}}
  \ , \\ \label{tr}
  \triangle t &=& t-t_0\,=\,\int_{r(t_0)}^{r(t)} dr\,
  \sqrt{\frac{\mu}{2\left[E_{\rm cm}-U_{\rm eff}(r)\right]}}\ . 
\end{eqnarray}
To get $\phi'(t)$ we use~(\ref{L}):
\begin{eqnarray} \nonumber
  d\phi' &=& \frac{L\,dt}{\mu\,r^2} = \frac{L\,dr}{r^2\sqrt{2\mu
  \left[E_{\rm cm}-U_{\rm eff}(r)\right]}} \\ \label{phipt_int} 
  \triangle\phi' &=& \phi'(t)-\phi'(t_0)\\ \nonumber
  &=& \int_{r(t_0)}^{r(t)} \frac{L\,dr}{r^2 \sqrt{2\mu
  \left[E_{\rm cm}-U_{\rm eff}(r)\right]}}\,, 
\end{eqnarray}
where due to the choice of the $\Sigma'$ coordinate frame of
Eq.~(\ref{ep1}) to (\ref{ep3}) we have
\begin{equation} \label{phip0}
  \phi'(t_0) = 0\ .
\end{equation}

\subsection{Kepler problem}

We specialize now to the gravitational potential
\begin{equation} \label{Ug}
  U(r) = -\frac{\alpha}{r}~~{\rm with}~~\alpha=G\,m_1 m_2\,.
\end{equation}
The effective potential becomes
\begin{equation} \label{Ugeff}
  U_{\rm eff}(r) = -\frac{\alpha}{r} + \frac{L^2}{2\mu\,r(t)^2}\,.
\end{equation}
%
\begin{figure}[tb] \begin{center} % See subfolder FigKeplerUeff.
\epsfig{figure=fig/KeplerUeff.eps,width=\columnwidth} \vspace{-1mm}
\caption{Potential and effective potential for the Kepler problem. 
\label{fig_Ueff}} \end{center} \vspace{-3mm} \end{figure}
%
An example is shown in Fig.~\ref{fig_Ueff}. For $E_{\rm cm}\ge 0$ 
the distance between the masses escapes to infinity, while for 
$E_{\rm cm}<0$ it is confined between $r^{\min}$ and $r^{\max}$. 
There are no solutions for $E_{\rm cm}<U_{\rm eff}^{\min}$ with 
\begin{equation} \label{Ueffmin}
  U_{\rm eff}^{\min} = U_{\rm eff}(r^0_{\min}) = 
  -\frac{\alpha^2}{2} \frac{\mu}{L^2}\,,~~
  r^0_{\min} = \frac{L^2}{\alpha\mu}\ .
\end{equation}
At $r^0_{\min}$ the orbit is a circle. For (\ref{Ug}) the solution of 
the integral (\ref{phipt_int}) is elementary
\begin{eqnarray} \label{phipr} \nonumber
  \phi'(t) &=& \arccos\left( \frac{L/r(t)-\mu\,\alpha/L}
  {\sqrt{2\mu\,E_{\rm cm}+\mu^2\alpha^2/L^2} } \right) \\ 
  &-& \arccos\left( \frac{L/r(t_0)-\mu\,\alpha/L}
  {\sqrt{2\mu\,E_{\rm cm}+\mu^2\alpha^2/L^2} } \right)
\end{eqnarray}
We map on the coordinate conventions of the literature \cite{MT,LL}
in which the $r_{\min}$ position, called {\it perihelion} or {\it 
pericenter}, is at $\phi=0$. This defines the the azimuth angle
\begin{equation} \label{phi}
  \phi(t) = \arccos\left( \frac{L/r(t)-\mu\,\alpha/L}
  {\sqrt{2\mu\,E_{\rm cm}+\mu^2\alpha^2/L^2} } \right)\,,
\end{equation}
so that
\begin{equation} \label{phip_def}
  \phi'(t) = \phi(t) - \phi(t_0)
\end{equation}
holds for the previously introduced $\phi'(t)$. The angle $\phi(t)$ 
is defined with respect to a coordinate system, which we call 
$\Sigma_{\rm l}$ (l for literature). It is rotated with respect 
to $\Sigma'$ by $\phi(t_0)$ about the common $\hat{e}'_3$ axis.

We can write Eq.~(\ref{phi} as
\begin{equation} \label{cosphi}
  cos[\phi(t)] = \frac{L/r(t)-\mu\,\alpha/L}
  {(\mu\alpha/L)\sqrt{1+2 E_{\rm cm}L^2/(\mu\alpha^2)}}\,.
\end{equation}
With the definitions of the parameter $p$, $2p$ is called the {\it latus
rectum}, and the {\it eccentricity} $e$,
\begin{equation} \label{p_e}
  p = \frac{L^2}{\mu\alpha}~~{\rm and}~~
  e = \sqrt{1+\frac{2\,E_{\rm cm}\,L^2}{\mu\,\alpha^2}}\,,
\end{equation}
the orbit in $\Sigma_{\rm l}$ becomes
\begin{equation} \label{orbit}
  \frac{p}{r(t)} = 1 + e\,\cos[\phi(t)]~~{\rm or}~~
  r(t) = \frac{p}{1+e\,\cos[\phi(t)]}\ .
\end{equation}
This is a conic section with the coordinate origin in a focal point.
The shortest distance from the focal point is
\begin{equation} \label{rmin}
  r_{\min} = \frac{p}{1+e}\,,
\end{equation}
which corresponds to the pericenter. For $U^{\min}_{\rm eff}< E_{\rm cm} 
< 0$ the orbit is an ellipse and the largest distance from the focal 
point is
\begin{equation} \label{rmax}
  r_{\max} = \frac{p}{1-e}\,,
\end{equation}
which is corresponds to the {\it apocenter}. The {\it turning points}
$r_{\min}$ and $r_{\max}$ of the orbit are also called {\it apsides}. 
They are indicated in Fig.~\ref{fig_Ueff}. The large half-axis is
\begin{equation} \label{a}
  a=\frac{p}{1-e^2}=\frac{\alpha}{2|E_{\rm cm}|}
\end{equation}
and the small half-axis
\begin{equation} \label{b}
  b=\frac{p}{\sqrt{1-e^2}}=\frac{L}{\sqrt{2\mu|E_{\rm cm}|}}\ .
\end{equation}
The ellipse becomes a circle for $e=0$ ($E_{\rm cm}=U^{\min}_{\rm eff}$). 
The position and velocity vectors with be orthogonal $\vec{r}\cdot
\vec{v}=0$, and we have a relation between their magnitudes. With 
\begin{eqnarray} \label{r0}
  r_0 = r(t_0) &=& |\vec{r}(t_0)|\,,\\ \label{v0}
  v_0 = v(t_0) &=& |\vec{v}(t_0)|\,,
\end{eqnarray}
the relation $E_{\rm cm}=U_{\rm eff}^{\min}$ yields:
\begin{eqnarray} 
  \frac{\mu}{2}\,(v_0)^2 - \frac{\alpha}{r_0} = - \frac{\alpha^2}{2}
  \frac{1}{\mu\,(r_0)^2 (v_0)^2}\,,\\
  (v_0)^4 - \frac{2\alpha}{\mu\,r_0}\,(v_0)^2 + 
  \frac{\alpha^2}{\mu^2 (r_0)^2} = 0\,.
\end{eqnarray}
The argument of the $\pm\sqrt{}$ part of the solution turns out to be 
zero, so that we end up with the unique result
\begin{equation} \label{v0circ}
  v_0 = \sqrt{\frac{\alpha}{\mu\,r_0}}
\end{equation}
and $\vec{v}_0\cdot\vec{r}_0=0$.

For $E_{\rm cm}>0$ ($e>1$) the orbits are hyperbolic and escape to 
infinity for the two solutions for the equation $1+e\,
\cos(\phi_{\infty})=0$. For $E_{\rm cm}=0$ the orbit is parabolic 
and $\phi_{\infty}=\pi$.

\begin{figure}[tb] \begin{center} % Data created by running 
% phy4241orbits.f in KeplerExact. Plot with phy4241.plt.
\epsfig{figure=fig/Korbits.eps,width=\columnwidth} \vspace{-1mm}
\caption{Elliptical and hyperbolic orbits corresponding to the initial 
conditons of table~\ref{tab_Korbits}. \label{fig_Korbits}} \end{center} 
\vspace{-3mm} \end{figure}

\begin{table}[th]   % Generated by porbits.f in KeplerExact. Different
\centering \bigskip % numbers (here->there, e ellipse, h hyperbola): 
                    % 1->14 e, 2->8 h, 3->2 e, 4->12 h.
\begin{tabular}{|c|c|c|c|c|c|c|c|c|} \hline
 \multicolumn{1}{|c|}{ } & \multicolumn{2}{c|}{Masses}  & 
 \multicolumn{3}{c|}{Initial Positions}  &
 \multicolumn{3}{c|}{Initial Velocities}              \\ \hline
\#&$i$&$m_i$ &$x^1_{i,0}$ &$x^2_{i,0}$    &$x^3_{i,0}$
&$\dot{x}^1_{i,0}$&$\dot{x}^2_{i,0}$&$\dot{x}^3_{i,0}$\\ \hline
1&1&0.651& 0.585&-0.238&-0.755&-0.828&-0.865&-0.726\\ \hline
 &2&0.931&-0.096& 0.000& 0.357&-0.209& 0.107&-0.660\\ \hline
2&1&1.510& 0.460&-0.359&-0.234&-0.918&-0.941&-0.323\\ \hline
 &2&0.126&-0.066&-0.090&-0.809& 0.789& 0.788& 0.620\\ \hline
3&1&1.328&-0.125& 0.898& 0.194&-0.452& 0.172& 0.125\\ \hline
 &2&1.999&-0.449&-0.085&-0.454&-0.976&-0.990&-0.968\\ \hline
4&1&0.180& 0.204&-0.968&-0.753&-0.811&-0.632& 0.784\\ \hline
 &2&1.560&-0.889&-0.979& 0.854&-0.323&-0.774&-0.533\\ \hline
\end{tabular}
\caption{Examples of initial conditions corresponding to orbits
shown in Fig.~\ref{fig_Korbits}.}\label{tab_Korbits} \end{table} 
% \vspace*{0.2cm}

In Fig.~\ref{fig_Korbits} we give examples of orbits in the
$\Sigma_{\rm l}$ frame, which correspond to the initial conditions
and velocities compiled in table~\ref{tab_Korbits}. Crosses in the 
figure indicate the corresponding initial positions in the 
$\Sigma_{\rm l}$ frame. 
\relax

\begin{thebibliography}{19}

\bibitem{LL} L.D. Landau and E.M. Lifschiz, {\it Lehrbuch der
Theoretischen Physik I}, Akademie Verlag, Berlin 1962.

\bibitem{MT} J.B. Marion and S.T. Thornton {\it Classical Dynamics
of Particles and Systems}, Harcourt Brace \& Company, Orlando, FL 1995.

% \bibitem{paper_quote} L.D. McLerran and B. Svetitsky, Phys.\ Lett.\ B
% {\bf 98}, 195 (1981); % --198

\end{thebibliography}

\clearpage
\end{document}
