Hm, this seems to be somewhat overcomplicated. Although the most beautiful way to understand the Kepler problem is analytical Hamiltonian mechanics with the exploitation of the full SO(4) symmetry of the problem, it's also quite easily possible to use only "elementary Newtonian mechanics". This is what the paper does, but it can be somewhat simplified by using modern vector notation.
We treat the simplified Kepler problem with a fixed center (i.e., the limit of a "sun" with infinite mass). The more general case of the interacting two-body problem is only slightly more cumbersome, because it can be mapped to the fixed-center problem by introducing center-mass and relative coordinates.
Putting the sun in the origin of a Cartesian coordinate system leads to the equation of motion
m \ddot{\vec{r}}=-\frac{\gamma m M}{r^2} \frac{\vec{r}}{r} \qquad (1).
Since we have a central force, the angular momentum is conserved. This follows immediately by taking the cross product with the radius vector:
m \vec{r} \times \ddot{\vec{r}}=m \frac{\mathrm{d}}{\mathrm{d} t} (\vec{r} \times \dot{\vec{r}})=0 \; \Rightarrow \; \vec{L}=m \vec{r} \times \dot{\vec{r}}=\text{const}.
Assuming that \vec{L} \neq 0 we can choose the coordinate system such that the z axis is along the direction of L. Since both \vec{r} and \dot{\vec{r}} are perpendicular to this fixed direction, the motion is indeed restricted to the xy plane, and we can introduce polar coordinates,
\vec{r}=r (\cos \theta,\sin \theta,0) \qquad (2)
and the velocity is
\dot{\vec{r}}=\dot{r} (\cos \theta,\sin \theta,0)+r \dot{\theta} (-\sin \theta,\cos \theta,0).\qquad (3)
It follows
\vec{L}=m \vec{r} \times \dot{\vec{r}}=m r \dot{\theta} \vec{e}_z.
Thus
L=m r \dot{\theta} = \text{const}. \qquad (4)
This is already Kepler's 2nd law: The radius vector between sun and planet sweeps out equal areas per unit time, i.e., the "area velocity" is constant. This follows directly from the meaning of the cross product.
To get the form of the orbit, we use the conservation of energy, which follows from taking the dot product of (1) with \vec{r}:
m \dot{\vec{r}} \cdot \ddot{\vec{r}}=\frac{m}{2} \frac{\mathrm{d}}{\mathrm{d} t} \dot{\vec{r}}^2=-\gamma m M \frac{\dot{\vec{r}} \cdot \vec{r}}{r^3}.
Using (2) and (3) gives
\frac{m}{2} \frac{\mathrm{d}}{\mathrm{d} t} \dot{\vec{r}}^2=-\gamma m M \frac{\dot{r}}{r^2}=-\gamma m M \frac{\mathrm{d}}{\mathrm{d} t} \left (\frac{1}{r} \right ).
This leads to
E=\frac{m}{2} \dot{\vec{r}}^2-\frac{\gamma m M}{r}=\text{const}.
Of course
V(r)=-\frac{\gamma m M}{r}
is the potential for the gravitational force and E the total energy.
Using (3) again leads to
E=\frac{m}{2} (\dot{r}^2 + r^2 \dot{\theta}^2)-\frac{\gamma m M}{r}.
With help of (4) this can be written as
E=\frac{m}{2} \dot{r}^2+\frac{L^2}{2 m r^2} -\frac{\gamma m M}{r}. \qquad (5)
This equation is not solvable in closed form, but it is possible to get the form of the orbit as a polar equation.
r=r(\theta).
To that end we use (4) again to write
\dot{r}=r'(\theta) \dot{\theta}=r'(\theta) \frac{L}{m r^2}.
Plugging this into (5) we find
E=\frac{L^2}{2m r^4} r'^2+\frac{L^2}{2m r^2}-\frac{\gamma m M}{r}. \qquad(6)
Substituting
r=\frac{1}{s} \; \Rightarrow \; s'=-\frac{r'}{r^2}
into (6) we get
E=\frac{L^2}{2m} s'^2 + \frac{L^2}{2m} s^2 - \gamma m M s. \qquad (7).
Taking the derivative of this equation wrt. [\itex]\theta[/itex] gives
s' \left (\frac{L^2}{m} s'' + \frac{L^2}{m} s -\gamma m M \right)=0.
Because generally s' \neq 0 the bracket must vanish, and the linear differential equation has the general solution
s=\frac{\gamma m^2 M}{L^2} + C_1 \cos \varphi + C_2 \sin \varphi. \qquad (8)
Demanding that the perihelion (point of closest approach of the planet to the sun) at \varphi=0 leads to C_2=0.
The solution
r(\theta)=\frac{p}{1+\epsilon \cos \theta}
describes a conical section with a focus at the origin. Using (8) in (7) gives for the parameters
p=\frac{L^2}{\gamma m^2 M}, \quad \epsilon=\sqrt{1+\frac{2 E L^2}{\gamma^2 m^3 M^2}}.
This shows that the orbit is an ellipse for E<0, a parabola for E=0, and a hyperbola for E>0. The bound motion of a planet around the sun is thus always an ellipse (or for the special case that \epsilon=0 a circle). This is Kepler's 1st Law.
