As far as I understand the Kepler problem tries to predict the motion of two bodies under the influence of gravity.
Here the trajectory of each particle is chosen to be:
$$\vec{r}_1 = \vec{R}+\dfrac{m_2}{M}\,\vec{r} \qquad \vec{r}_2 = \vec{R}-\dfrac{m_1}{M}\,\vec{r}$$
where $\vec{R} = \dfrac{m_1\,\vec{r}_1+m_2\,\vec{r}_2}{M}$ is the the common barycentre, $\vec{r} = \vec{r}_1 -\vec{r}_2$ the relative Distance and $M = m_1+m_2$ the total mass.
Now on the way to the solution some simplifications have been made. For instance motion is just locked to the $x$–$y$-plane. Subsequently angular momentum $\vec{L}$ is constant and just non-zero in the $z$-combonent: $L_z = x\,\dot{y}-\dot{y}\,x $. Energy is defined as usual: $E = \dfrac{1}{2}\,\mu\vec{\dot{r}}^2+\dfrac{1}{2}\,M\,\vec{\dot{R}}^2+V$, $\quad\mu = \dfrac{m_1\,m_2}{M}$ being the reduced mass and $V$ the potential. Respectively using polar coordinates: $x = r\,\cos\varphi \quad y = r\,\sin\varphi$ the equations simplifies further more: $L_z
= \mu\,r^2\,\dot{\varphi} \quad E_r = \dfrac{1}{2}\,\mu\dot{r}^2+\dfrac{{L_z}^2}{2\,\mu\,r^2}+V(r)$.
At this point using differentials the equation of motion becomes:
$$\varphi-\varphi_0 = \int_{r_0}^r \dfrac{L_z}{r^2\,\sqrt{(E_r-V(r))\,2\,\mu-\dfrac{{L_z}^2}{r^2}}}\,\mathrm{dr}$$
During the Kepler problem the potential is predefined: $V(r) =
-\dfrac{\gamma}{r}$. In this case the integral can be solved and rearranged:
$$r(\varphi) = \dfrac{\rho}{1+\varepsilon\cos{\varphi}}\quad \text{where} \quad \varepsilon = \sqrt{\dfrac{2\,E_r\,{L_z}^2}{\gamma^2\,\mu}+1},\quad \rho = \dfrac{{L_z}^2}{\mu\,\gamma}$$
Now here is what the actual question is. How do I plot this thing? Apparently using the equations above it follows: $$\vec{r_1} = \left(\begin{array}{c}x_1 \\y_1\end{array}\right) = R+\dfrac{m_2}{M}\,\left(\begin{array}{c}r(\varphi)\,\cos\varphi \\r(\varphi)\,\sin\varphi\end{array}\right)\qquad \text{$\vec{r}_2$ by analogy}$$
But how about initial conditions? For example if the trajectory of one object moving is set to be $\vec{r}_a = \left(\begin{array}{c}a_1+t^2 \\ a_2+t\end{array}\right)$ and the other $\vec{r}_b = \left(\begin{array}{c}b_1-2\,t \\ b_2-1\end{array}\right)$ it seems not matching the formula for $R_1$, $L_z$ and $E_r$ thus $\vec{r}_1$ and $\vec{r}_2$ the way it was clarified earlier. Probably because time is now important. In general I do not understand how $\quad \vec{r}_1, \quad \vec{r}_2\quad $ are calculated at all, because therefore $\vec{R}$ needs to be determined, but it depends on $\vec{r}_{1,2}$ itself.
Best Answer
Two body equation of motions
\begin{align*} & \ddot{r} -r{\dot\varphi }^{2}+{\frac { \left( m+M \right) G}{{r}^{2}}}=0\\ &\ddot\varphi +{2\, \frac {\dot\varphi \,{\dot r}}{r}}= 0 \end{align*}
solve those differential equations with the initial conditions \begin{align*} &\textbf{vis-viva -equation}\\\\ &v=\sqrt{G\,(M+m)\,\left(\frac{2}{r}-\frac{1}{a}\right)}\\ &\text{with}\\ &r=\frac{a\,(1-e^2)}{1+e\,\cos(\varphi)}\\ &\text{the value of r at the perigee}\\ &r_0=r(\varphi=0)=\frac{a\,(1-e^2)}{1+e}=a\,(1+e)\\ &v_0=v(r=r_0)=\sqrt {G \left( M+m \right) \left( 2\,{\frac {1+e}{a \left( 1-{e}^{2} \right) }}-\frac{1}{a} \right) } \\\\ &\text{hence}\\ &r(0)=r_{0}\,,\dot r(0)=0\\\\ &\varphi(0)=0\,,\dot \varphi (0)=\frac{v_0}{r_0} \end{align*}
you obtain the position of the bodies
\begin{align*} &x_1=\frac{\mu}{m}\,r(t)\\ &x_2=-\frac{\mu}{m}\,r(t)\\ &\text{with}\\ &\mu=\frac{m\,M}{m+M} \end{align*}