# How to actually plot the Kepler Problem?

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.

• What is the problem? Can't you olot plot r$_i$ as a function of $\phi$, which a function of t? Nov 7, 2021 at 11:41

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*}

I'm not sure I understand the question... why are you trying to enforce that $$\vec{r}_1$$ and $$\vec{r}_2$$ are these arbitrary functions of time? You already have expressions for $$\vec{r}_1$$ and $$\vec{r}_2$$, and your initial conditions should enforce what the conserved quantities are.

If you specifically would like to plot the time-dependence, you should take the formula:

$$E_r=\frac{1}{2}\mu\dot{r}^2+\frac{L_z^2}{2\mu r^2}+V(r)$$

You can then manipulate this to

$$\frac{dr}{dt}=\sqrt{2\left(\frac{E_r}{\mu}-\frac{L_z^2}{2\mu^2r^2}-\frac{V(r)}{\mu}\right)}$$

Then the time can be related to the function $$r(\varphi)$$ (and thus $$\varphi$$) by the following integral, which can be done numerically:

$$t=\int dt=\int\frac{dr'}{\sqrt{2\left(\frac{E_r}{\mu}-\frac{L_z^2}{2\mu^2r'^2}-\frac{V(r')}{\mu}\right)}}$$