Unexpected eccentricity in moon orbit simulation

Question

I've used Mathematica's NDSolve function to calculate the orbit of the Moon around the Earth.

I used the following initial positions (perigee):

$\vec{r}_{Earth}=\begin{pmatrix}-\frac{m_{Moon}}{m_{Earth}}362.6\cdot10^{6}\\ 0\\ 0 \end{pmatrix}$

$\vec{r}_{Moon}=\begin{pmatrix}362.6\cdot10^{6}\\ 0\\ 0 \end{pmatrix}$

Where the origin is the barycentre

I then calculated the velocity vectors to be:

$\vec{v}_{Moon}=\begin{pmatrix}0\\ \sqrt{G\left(m_{Earth}+m_{Moon}\right)\left(\frac{2}{|r_{Moon}|}-\frac{1}{a}\right)}\\ 0 \end{pmatrix}$

$\vec{v}_{Earth}=\begin{pmatrix}0\\ -\frac{m_{Moon}}{m_{Earth}}\sqrt{G\left(m_{earth}+m_{Moon}\right)\left(\frac{2}{|r_{moon}|}-\frac{1}{a}\right)}\\ 0 \end{pmatrix}$

Here I used the preservation of momentum to find the Earth's wobble's velocity (I'm not 100% sure about this, but it seems to work).

Where

$a=384.399\cdot10^{6}$

the semi-major axis.

I then used Mathematica to plot the orbit of the Earth and the Moon around the barycentre.

Mathematica file here. The notebook is a generic two-body code hence the sometimes circuitous code.

Here's my problem (image):

http://i16.photobucket.com/albums/b21/ApocalypseVolcano/Aplots.png

The apogee is nowhere near the distance from the barycentre as it should be.

I get that other bodies in the solar system influence the orbit of the Earth and the Moon, but is that what's solely responsible?

Also, if I were to have the Earth orbit the Sun, can I use the combined mass of the Earth and the Moon and merely have this system's barycentre orbit the Sun?

The notebook works well for earlier tries with a circular orbit.

Answer 1

You have four big problems and two small problems. The big problems are that you are initializing the Earth's and Moon's initial position and velocity incorrectly. The initial distance between the Earth and Moon is off by a factor of 1.0123, as is the initial relative velocity. The small problems are (1) an incorrect value for the Earth-Moon semi-major axis and (2) your use of $G$, $m_{moon}$, and $m_{earth}$.

The big problems were the key cause of your larger than expected apogee distance.

The small problems: You should want to fix those as well.

Issue #1: You are using 384,399 km as the length of the semi-major axis of the the orbit of the Moon. That is incorrect. That value is the inverse sine parallax of the Moon, the inverse of the mean value of the inverse of the distance. A better value for the semi-major axis length is 385,000 km, which is the mean distance between the Earth and Moon. An even better value is 384,748 km, from Chapront-Touzé, M., & Chapront, J. (1983). The lunar ephemeris ELP 2000. Astronomy and Astrophysics, 124, 50-62.

Issue #2: You are using the product $G (m_{earth}+m_{moon})$. You have obliterated your accuracy when you do that. Solar system astronomers instead use what are called the "standard gravitational parameters" to describe the masses of the Sun, the planets, and our Moon. Conceptually, this is just the product $\mu_{body} = G m_{body}$. There's a big difference, though, between using $\mu_{body}$ and $G m_{body}$. Scientists know many of those gravitational parameters to six places of accuracy or more. On the other hand, the gravitational constant G is known a paltry four places. Even more importantly, when you use G and mass you are almost certainly using values that are inconsistent with one another. Use the standard gravitational parameters. You can find a list of them at this wikipedia article.

---

So what should you do?

Denoting

• $r_p=362600\,\text{km}$ as the distance between the Earth and Moon at perigee,
• $a= 384748\,\text{km}$ as the semi-major axis of the Earth and Moon about one another,
• $u_e=398600.4418\,\text{km}^3/\text{s}^2$ as the Earth's standard gravitational parameter,
• $u_m=4902.8000\,\text{km}^3/\text{s}^2$ as the Moon's standard gravitational parameter, and
• $v_p = \sqrt{(\mu_e+\mu_m)\left(\frac 2 {r_p} - \frac 1 a\right)}\,$ as the relative velocity between the Earth and Moon at perigee, per the vis-viva equation.

You need to emplace the Earth and Moon such that the distance and velocity between them are $r_p$ and $v_p$. Since you want the barycenter to be at the origin, one way to do this is

$$\begin{aligned} r_{moon} &= \phantom{-} \frac {r_p} {1+\mu_m/\mu_e} \hat x \\ v_{moon} &= \phantom{-} \frac {v_p} {1+\mu_m/\mu_e} \hat y \\ r_{earth} &= -\frac {r_p} {1+\mu_e/\mu_m} \hat x \\ v_{earth} &= -\frac {v_p} {1+\mu_e/\mu_m} \hat y \end{aligned}$$

Finally, you should compute your accelerations using the standard gravitational parameters rather than using G*M.