<p>I don't have enough information to determine what is wrong with your spreadsheet, but I wrote the same program in Fortran, trying to follow your description. My results seem to be quite accurate, considering that Leapfrog is only second order and that we are neglecting Saturn, the Earth and the other planets.
You can try and compare my results with yours, to understand what is the problem.</p>
<p>I chose to use Fortran because it is very readable, almost like writing math formulae directly, therefore you shouldn't have problems comparing my formulae with yours. Lines starting with an exclamation mark (!) are comments.</p>
<p>The program is the following.</p>
<pre><code>program ThreeBody
implicit none
! UNITS used:
! Length: AU
! Time: Julian year

! Declare all the variables
integer, parameter:: dp=kind(0.d0)
! 3D position vector of Mars, the Sun and Jupiter
real(dp) :: r_m(3), r_s(3), r_j(3)
! 3D velocity vector of Mars, the Sun and Jupiter
real(dp) :: v_m(3), v_s(3), v_j(3)
! 3D acceleration vector of Mars, the Sun and Jupiter
real(dp) :: a_m(3), a_s(3), a_j(3)
! Mass of Jupiter, mass of the Sun
real(dp) :: GM_j = 0.037685, GM_s = 39.478418  ! au^3 yr^-2 
! Integration time step
real(dp) :: eps = 0.000684 ! 6 hours in Julian years
! index of current step, unit of the output file
integer :: i, outputFile

! Initial conditions from JPL Horizons 2000-Jan-01 12:00:00.0000 TDB
! Mars
r_m = [1.383579466402647E+00, 1.621204117825275E-02, -3.426152579691657E-02] ! au
v_m = [6.768779495222499E-04, 1.517984118259482E-02, 3.015574279972501E-04]  ! au/day
v_m = v_m*365.25 ! au/yr

! Jupiter
r_j = [3.994040256427818E+00, 2.935779784974931E+00, -1.015789849456240E-01]
v_j = [-4.562948121087195E-03, 6.435847855854325E-03, 7.548691257292337E-05]
v_j = v_j*365.25

! The Sun
r_s = [-7.137179161607904E-03, -2.795997495567612E-03, 2.062985061910257E-04]
v_s = [5.378460339618783E-06, -7.406916207973516E-06, -9.434293432137470E-08]
v_s = v_s*365.25

! open output file
open(newunit=outputFile, file=&quot;output.csv&quot;, status=&quot;replace&quot;)
! Write header to file
write(outputFile, '(A)', advance=&quot;no&quot;) &quot;time,&quot;
write(outputFile, '(A)', advance=&quot;no&quot;) &quot;rx_m,ry_m,rz_m,rx_j,ry_j,rz_j,rx_s,ry_s,rz_s,&quot;
write(outputFile, '(A)', advance=&quot;no&quot;) &quot;vx_m,vy_m,vz_m,vx_j,vy_j,vz_j,vx_s,vy_s,vz_s,&quot;
write(outputFile, '(A)')               &quot;ax_m,ay_m,az_m,ax_j,ay_j,az_j,ax_s,ay_s,az_s&quot;

! Main loop. Here we repeat the Leapfrog step 10000 times
do i = 1, 10000
    ! Calculate the acceleration
    ! Ignoring the effect of Mars on the Sun and on Jupiter
    a_m = -GM_j * (r_m - r_j)/norm2(r_m - r_j)**3 - GM_s * (r_m - r_s)/norm2(r_m - r_s)**3
    a_j = -GM_s * (r_j - r_s)/norm2(r_j - r_s)**3
    a_s = -GM_j * (r_s - r_j)/norm2(r_s - r_j)**3

    ! Calculate the velocity at the 1/2 step
    if (i == 1) then
        ! If we are at the first step, use eps/2 instead of eps
        v_m = v_m + eps/2 * a_m
        v_j = v_j + eps/2 * a_j
        v_s = v_s + eps/2 * a_s
    else
        v_m = v_m + eps * a_m
        v_j = v_j + eps * a_j
        v_s = v_s + eps * a_s
    end if
    
    ! Calculate the position
    r_m = r_m + eps * v_m
    r_j = r_j + eps * v_j
    r_s = r_s + eps * v_s

    ! Write the output to a csv file
    write(outputFile, '(*(G0.6,:,&quot;,&quot;))') i*eps, r_m, r_j, r_s, v_m, v_j, v_s, a_m, a_j, a_s
    
end do

close(outputFile)

end program
</code></pre>
<p>To run it, you can copy and paste it into a text file called <code>ThreeBody.f</code>. If you are on Linux, you just have to open a terminal in the directory where you have saved the file and run</p>
<pre><code>gfortran ThreeBody.f &amp;&amp; ./a.out
</code></pre>
<p>If you are on Windows or Mac, you will probably first need to <a href="https://fortran-lang.org/learn/os_setup/install_gfortran" rel="nofollow noreferrer">install gfortran</a> or any other compiler of your choice.</p>
<p><strong>Alternatively, you can just <a href="https://tio.run/##nZZbr9M4EIDf@ytMX2hLGnyJ7bgSD7CgBQSrlThoHxAceRu3J900iZw2UP782ZnEbdOyLIiqSuvJ5Jvx3OJV5XfelvP1qv9zf1/7au3tltzceeeeVdlhROCTb@siX@Y7Ulal6yQPyPs/Xt28I/vGZYsgeePK9e5uQZ6@D4KbfOsW5OCsHwXJc7csrHfEFgXZ3TnSWp/bvwvX9GbKnVs7H5Hagg9u5/xiQbL6yT95mU1onNFpp@adLSZZPR3oEVCsc/KEiJglTBqupJCp0UZwkSaKq0SI4IJ4TuqqyXd5VZLWLXeVJ9WKvLW@iTqX3u1LYsuMvN7XOZAvLKIZf7udCLDtb5vwu4HfM7x1RQWhOvwSvA3wNsDbK7hdLl3hvP1l720wYIMBOzTw1jYN4sLDEdkGQSB/Q/v97e0Ggk5jKrRKZYSCBrNg4kSnCUuRavefBDn4T3NOgp1XmOewiR3UCGl2rv4G7uqmZ1Oq0gQeU@Su2vsGygQ8LHJbdqXVBCjUiPuCzi733rtyd4aGskJmPjq5ABVgC7KsyqyrhYasfLUlr/98Q15WPv@KEg6m56@hPygjjC8ohS@6Q8nN82enmAUPoC7A3Q8sFqmQ2iRKJRTqTr94RGlEWKw44zRhTKdcci1fzCmPyFzECVdMcnhCGaakRvnHLmodte2pKob4aqAayTlPjAG1BKmSaQOBZgA1Scp7qogpk1InXBujuaQMtT/2uXic2cOADNeZUDLmsr97OLXqRQV1af4gYmMSCtuSCtgpS/u98djAjrWBjBudGMF68ZyhGzoF/0wiFU8ouME@BuMdcJ7EcAMcZ5zRVDODYRERUXEC7QslJGUqE8GDWMcySSFMEC1uuBAYLDkEwjXsphMex9Cwerv6/DDXMRMaYscUU1QbmvQW5jzWRhrchpRKK8Z7OY8p@JlKqphhsH/dRTRY7ogyFrB/RYUwiqU6FaChAKjjhGJmOVjRQjIV5NAhIuFGwO7Ak0RjcNIhEa7HvRxT8peHjJA7Z7OQl88omMwi8nDydPoQujlrbbl0T8ZlNZ6SMfZWNP4ZTf/ldhv5A16@4uXL7QaXG1xucNngssFl83PEFoktElsktkhskdgisUVii8T2e8QpufyMLRItEi0SLRItEi0SLRItEmHZjEen9oRpUVRVHZOXDl48nx3MmNrZXTfT3jhbQ9uvu1lBWNfaGLS@o7OK4CuFRXhndPTjAfnNFst9YSEVyBiO5IHSq3UJY6RcdzputYJZfZzUBMfeYFZX5UWz4cd2zTnv5uuMTHC2zLELp4@BuuUDwWwm4F83dweKzbVig4oD/Cbgw1Ob66c2//lUc@lUc@1Uc3Zq9N2AnV6QIQnsMT/P6m5er8gEAg@Rn6LCOaohsivMYneK6AGr3Df9uI/wPIIvDkDmJUhshkEHwQXjNPrIo6A8w4hf6YSBcqGzudIJjXqhc7blisb9n@Efm/2x0WuTUE/56rT8bhKOR6CTQv/68gPX2oFr/QvAD9xqB271Q9UPXAL3hgXQDy60W@139R7SVhFLlk0LuSvOIRoMgNnkdxqraBGNo/EUh0E@AzKetrbdkas7f@H5aNsdkroTE55mtt2RJrqIyegYmKwaXazDYff@/l8" rel="nofollow noreferrer">Try it online!</a></strong></p>
<p>The output of the program is a Comma Separated Value (CSV) file called <code>output.csv</code> that you can read as a spreadsheet with Excel or LibreOffice. It contains 28 columns:</p>
<ul>
<li>1 time [yr]</li>
<li>2-4 3D position of Mars [au]</li>
<li>5-7 3D position of Jupiter [au]</li>
<li>8-10 3D position of the Sun [au]</li>
<li>11-19 same thing, but for the velocity [au/yr]</li>
<li>20-28 same thing, but for the acceleration [au/yr^2]</li>
</ul>
<p>Just for reference, here is the 3D position of Mars for the first 10 steps of integration</p>
<p><a href="https://i.sstatic.net/2IzDh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/2IzDh.png" alt="first 10 steps of integration" /></a></p>
<p>And here it is after about 1 Julian year</p>
<p><a href="https://i.sstatic.net/BUXnc.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/BUXnc.png" alt="Mars position after 1 year" /></a></p>
<h3>Remarks on the accuracy of the result:</h3>
<p>I don't know what accuracy you were expecting. In the comments I read that you considered <span class="math-container">$GM = 4\pi^2$</span> to be incorrect for your required level of accuracy. It best measured value can be found in the Wikipedia page for the <a href="https://en.wikipedia.org/wiki/Standard_gravitational_parameter" rel="nofollow noreferrer">Standard gravitational parameter</a>, and, when expressed in units of <a href="https://en.wikipedia.org/wiki/Julian_year_%28astronomy%29" rel="nofollow noreferrer">Julian years</a> and <a href="https://en.wikipedia.org/wiki/Astronomical_unit" rel="nofollow noreferrer">AU</a>, it is equal to <span class="math-container">$4\pi^2$</span> to at least one part per 10000.</p>
<p>If your required level of accuracy is smaller than this, then you have to be careful that 6 hours are not <span class="math-container">$0.000684$</span> Julian years, but (with some more digits) <span class="math-container">$0.00068446269$</span> years. A better value for the gravitational parameter of the Sun would be <span class="math-container">$GM_\odot = 39.476926414 au^3 yr^{-2}$</span>, and also the one of Jupiter could gain some more digits <span class="math-container">$GM_J = 0.037684447 au^3 yr^{-2}$</span>. But, in any case, if you really need to know the position of Mars so precisely, then you probably need to include at least Saturn, the Earth and probably Venus in your calculation. They could have a far greater influence on the result, than the 6th decimal digit of the gravitational parameter.</p>

CC BY-SA 4.0
2022-07-26T11:01:08.650Z
+1! You had me at "_but I wrote the same program in Fortran_." :D
pela
CC BY-SA 4.0
2022-07-26T15:18:25.113Z
Thanks. A very detailed response that I'm still trying to get my head around. I also can't quite believe that you were able to write an equivalent program in hours (maybe less?) when it's taken me weeks/months to put together my own dodgy spreadsheet :-) I really do need to learn some coding.
Peter
CC BY-SA 4.0
2022-07-26T16:34:06.543Z
@Peter That's the advantage of writing code in a programming language versus a spreadsheet I suppose :D. I also did it because you say you are new to coding, and starting can be difficult if you don't have a guide. I encourage you to take this template and extend it to suit your needs. Please feel free to ask if you don't understand some of the code or have trouble compiling/running
Prallax
CC BY-SA 4.0
2022-07-26T19:43:54.930Z
I notice you've used the solar system barycentre as coordinate centre instead of (as I have) the Sun (body center). That's more accurate, I assume?
Peter
CC BY-SA 4.0
2022-07-26T20:01:28.167Z
Also, your intitial Mars y values are a little off compared to the JPL values?
Peter
CC BY-SA 4.0
2022-07-26T21:00:07.763Z
@Peter The coordinates centered on the Sun are not an inertial system, therefore I would have had to take into account apparent forces. I chose to simplify my life and go for the barycenter instead :) I don't know why you see the y values off compared to JPL. It might be because I used the planet center (499) instead of barycenter (4). If that's the case, it's a sloppiness on my part, thank you for pointing it out.
Prallax
CC BY-SA 4.0
2022-07-27T10:05:34.507Z
My head is spinning with all this but I think what I've done is reversed the sign of your a_s (acceleration of Sun towards Jupiter) and considered it as an additional acceleration of Jupiter towards the Sun. The total acceleration of Mars then equals the acceleration of Mars towards the Sun, plus acceleration of Mars towards Jupiter, minus your a_s. I hope I've got that right - it seems to work, anyway. I've also just noticed that your initial y value for Mars should be negative. I've posted an “answer” listing some of my results. Many thanks.
Peter
CC BY-SA 4.0
2022-07-26T09:48:11.217Z
2022-07-26T12:57:22.577Z
2022-07-26T12:57:22.577Z
Prallax
<p>Apologies for the length of this question.</p>
<p>I'm modelling the motion of a single planet around the Sun using the equations:</p>
<p><span class="math-container">$$dv_{x_{i}}/dt=a(t)=-GM_{SUN}x_{i}/r^{3},$$</span></p>
<p>where <span class="math-container">$x_{1}=x$</span>, <span class="math-container">$x_{2}=y$</span> and <span class="math-container">$x_{3}=z$</span>. Also,</p>
<p><span class="math-container">$$r=\sqrt{x^{2}+y^{2}+z^{2}}.$$</span></p>
<p>I'm using the leapfrog numerical method of solving these equations (from Feynman, <a href="https://www.feynmanlectures.caltech.edu/I_09.html" rel="nofollow noreferrer">here</a>), where:</p>
<p><span class="math-container">$$x\left(t+\epsilon\right)=x\left(t\right)+\epsilon v\left(t+\epsilon/2\right)$$</span></p>
<p><span class="math-container">$$v\left(t+\epsilon/2\right)=v\left(t-\epsilon/2\right)+\epsilon a\left(t\right)$$</span></p>
<p><span class="math-container">$$v\left(\epsilon/2\right)=v\left(0\right)+\left(\epsilon/2\right)a\left(0\right).$$</span></p>
<p>Unable to code, I've put together a spreadsheet to model the orbit of Mars and (separately) Jupiter. I'm using units of au and Julian years. I'm using a time interval of <span class="math-container">$\epsilon = 6 \text{ hours} = 0.000684 \text{ yr}$</span>. I've used (first row) initial state vector data <span class="math-container">$\left(x,y,z,v_{x},v_{y},v_{z}\right)$</span> from JPL Horizons, with an arbitrary start date of 2000-Jan-01 12:00:00.0000 TDB. <span class="math-container">$GM_{SUN} = 39.478418\,\text{au}^{3}\text{yr}^{-2}$</span>. The results are pretty accurate for both Mars and Jupiter.</p>
<p>So far so good. The problem starts when I try to factor in Jupiter's gravitational effect on Mars' orbit by using the equations:</p>
<p><span class="math-container">$$\frac{dv_{ix}}{dt}=\sum_{j=1}^{N}-\frac{GM_{j}\left(x_{i}-x_{j}\right)}{r_{ij}^{3}}$$</span></p>
<p>and so on for <span class="math-container">$y$</span> and <span class="math-container">$z$</span>. Here, the index <span class="math-container">$i$</span> refers to Mars, and the index <span class="math-container">$j$</span> refers to the Sun and Jupiter. Also,</p>
<p><span class="math-container">$$r=\sqrt{\left(x_{i}-x_{j}\right)^{2}+\left(y_{i}-y_{j}\right)^{2}+\left(z_{i}-z_{j}\right)}.$$</span></p>
<p>To keep things simple, I'm ignoring the gravitational effect of Jupiter on the Sun. I find the <span class="math-container">$x$</span> component of Mars' acceleration toward Jupiter at <span class="math-container">$t\left(0\right)$</span>, which I'll call <span class="math-container">$MJa_{x}\left(0\right)$</span> using</p>
<p><span class="math-container">$$MJa_{x}\left(0\right)=-GM_{JUPITER}\frac{\left(x_{MARS}\left(0\right)-x_{JUPITER}\left(0\right)\right)}{\left(r_{MARS-JUPITER}\left(0\right)\right)^{3}},$$</span></p>
<p>where <span class="math-container">$x_{MARS}\left(0\right)$</span> and <span class="math-container">$x_{JUPITER}\left(0\right)$</span> values are from JPL Horizons and <span class="math-container">$GM_{JUPITER}$</span> equals <span class="math-container">$0.037685\,au^{3}yr^{-2}$</span>. And then, summing the accelerations, find</p>
<p><span class="math-container">$$v_{x}\left(\epsilon/2\right)=v_{x}\left(0\right)+\left(\epsilon/2\right)\left(a_{x}\left(0\right)+MJa_{x}\right(0)),$$</span></p>
<p>where <span class="math-container">$v_{x}\left(0\right)$</span> is the <span class="math-container">$x$</span>-component of Mars' velocity vector wrt the Sun (obtained from JPL), and <span class="math-container">$a_{x}\left(0\right)$</span> is the <span class="math-container">$x$</span>-component of Mars' acceleration vector solely due to the Sun's gravitational attraction. The corrected position of Mars at <span class="math-container">$t\left(1\right)$</span> is then given by</p>
<p><span class="math-container">$$x_{MARS}\left(1\right)=x_{MARS}\left(0\right)+\epsilon v_{x}\left(\epsilon/2\right).$$</span></p>
<p>I then find the next value of <span class="math-container">$x_{MARS}\left(1\right)-x_{JUPITER}\left(1\right)$</span> using my previously calculated <span class="math-container">$x_{JUPITER}\left(1\right)$</span> value (from the Jupiter orbit spreadsheet). Next, I repeat the above to find <span class="math-container">$MJa_{x}\left(1\right)$</span>, the <span class="math-container">$x$</span> component of Mars' acceleration toward Jupiter at <span class="math-container">$t\left(1\right)$</span>. And so on.</p>
<p>I was rather pleased with all this except for the fact that it doesn't work, i.e. the results are consistently less accurate than the original Sun-Mars model.</p>
<p>Any idea where I'm going wrong? Thanks.</p>
<p>EDIT - I've now amended my model by using <span class="math-container">$x_{MARS}-x_{JUPITER}$</span> values calculated directly from my existing Mars and Jupiter orbit spreadsheets (ie simply subtract <span class="math-container">$x_{JUPITER}$</span> from <span class="math-container">$x_{MARS}$</span>. With those positions known, I then don't need to worry about calculating velocities as I can then use <span class="math-container">$x_{MARS}-x_{JUPITER}$</span> to directly find the accelerations <span class="math-container">$MJa_{x}$</span>, etc. I've also taken into account the gravitational effect of Jupiter on the Sun as implied in the above <span class="math-container">$\frac{dv_{ix}}{dt}$</span> equation. The results are ambiguous at best. Am I doing this numerical integration correctly?</p>

CC BY-SA 4.0
2022-07-23T18:36:07.253Z
Your GM values look a little bit off, but I guess they're ok if they give acceptable results when you compute the orbits of Mars & Jupiter separately. The GM values from DE440 for the planets are given on https://ssd.jpl.nasa.gov/astro_par.html The value for the Sun is 132712440041.279419.
PM 2Ring
CC BY-SA 4.0
2022-07-23T18:40:13.347Z
Your initial vector data should be for the planet barycenters, not the planets. That is, body IDs 3 & 4, not 399 & 499. That won't matter much for Mars, but it does for Jupiter. And maybe you *do* need to consider the force of Jupiter on the Sun, since the Sun-Jupiter barycenter distance is roughly equal to the Sun's radius.
PM 2Ring
CC BY-SA 4.0
2022-07-25T16:54:43.627Z
@PM2Ring - Yes, I assumed, incorrectly, that the Sun $GM =4\pi^{2}$, which I've now corrected. I tried using Jupiter barycentre values but with no obvious increase in accuracy.
Peter
CC BY-SA 4.0
2022-07-25T20:09:47.040Z
Hi, what do mean you when you say that the results are consistently less accurate? Are you talking about a few percent, or is it totally off? Could you show some plots?
Prallax
CC BY-SA 4.0
2022-07-25T20:22:20.803Z
Also, I assume the missing square in your seventh equation is a typo ;)
Prallax
CC BY-SA 4.0
2022-07-25T22:14:43.307Z
I admire you setting this stuff up in a spreadheet, hard work! You wrote that you don't code, but doing that's coding big-time, and spreadsheet language is so awkward! The problem: no mention of barycentric coordinates, so I assume you're working heliocentrically. Did you make sure to pair out the direct and indirect attractions? I.e., perturbing accel. of Mars due to Jupiter = accel. of Mars due to Jupiter MINUS accel. of Sun due to Jupiter. Omit the 2nd term, you get wrong results. If working barycentrically, the accels of the sun due to planets define its barycentric motion. Good luck!
terry-s
CC BY-SA 4.0
2022-07-25T22:24:26.970Z
Another point: are you setting up checksums to see what's going wrong? Examples, at each step, vector sum of every paired  interaction (EITHER: force of A on B + force of B on A, OR mass A x accel of A due to B + mass B x accel of B due to A) should be zero, if not there's an error.  Sum of all the orbital angular momentum vectors should be a constant vector. Total energy should sum to a constant. If working barycentrically, the center of mass of the whole system should stay constant. Hope that helps.
terry-s
<orbital-elements><n-body-simulations><ephemeris><modeling>
Orbital modelling, three bodies, leapfrog numerical integration
Peter
# Orbital modelling, three bodies, leapfrog numerical integration

(Asked by: username_0 on 2022-07-23)

Apologies for the length of this question.

I'm modelling the motion of a single planet around the Sun using the equations:

$$dv\_{x\_{i}}/dt=a(t)=-GM\_{SUN}x\_{i}/r^{3},$$

where $x\_{1}=x$, $x\_{2}=y$ and $x\_{3}=z$. Also,

$$r=\sqrt{x^{2}+y^{2}+z^{2}}.$$

I'm using the leapfrog numerical method of solving these equations (from Feynman, [here](https://www.feynmanlectures.caltech.edu/I_09.html)), where:

$$x\left(t+\epsilon\right)=x\left(t\right)+\epsilon v\left(t+\epsilon/2\right)$$

$$v\left(t+\epsilon/2\right)=v\left(t-\epsilon/2\right)+\epsilon a\left(t\right)$$

$$v\left(\epsilon/2\right)=v\left(0\right)+\left(\epsilon/2\right)a\left(0\right).$$

Unable to code, I've put together a spreadsheet to model the orbit of Mars and (separately) Jupiter. I'm using units of au and Julian years. I'm using a time interval of $\epsilon = 6 \text{ hours} = 0.000684 \text{ yr}$. I've used (first row) initial state vector data $\left(x,y,z,v\_{x},v\_{y},v\_{z}\right)$ from JPL Horizons, with an arbitrary start date of 2000-Jan-01 12:00:00.0000 TDB. $GM\_{SUN} = 39.478418\,\text{au}^{3}\text{yr}^{-2}$. The results are pretty accurate for both Mars and Jupiter.

So far so good. The problem starts when I try to factor in Jupiter's gravitational effect on Mars' orbit by using the equations:

$$\frac{dv\_{ix}}{dt}=\sum\_{j=1}^{N}-\frac{GM\_{j}\left(x\_{i}-x\_{j}\right)}{r\_{ij}^{3}}$$

and so on for $y$ and $z$. Here, the index $i$ refers to Mars, and the index $j$ refers to the Sun and Jupiter. Also,

$$r=\sqrt{\left(x\_{i}-x\_{j}\right)^{2}+\left(y\_{i}-y\_{j}\right)^{2}+\left(z\_{i}-z\_{j}\right)}.$$

To keep things simple, I'm ignoring the gravitational effect of Jupiter on the Sun. I find the $x$ component of Mars' acceleration toward Jupiter at $t\left(0\right)$, which I'll call $MJa\_{x}\left(0\right)$ using

$$MJa\_{x}\left(0\right)=-GM\_{JUPITER}\frac{\left(x\_{MARS}\left(0\right)-x\_{JUPITER}\left(0\right)\right)}{\left(r\_{MARS-JUPITER}\left(0\right)\right)^{3}},$$

where $x\_{MARS}\left(0\right)$ and $x\_{JUPITER}\left(0\right)$ values are from JPL Horizons and $GM\_{JUPITER}$ equals $0.037685\,au^{3}yr^{-2}$. And then, summing the accelerations, find

$$v\_{x}\left(\epsilon/2\right)=v\_{x}\left(0\right)+\left(\epsilon/2\right)\left(a\_{x}\left(0\right)+MJa\_{x}\right(0)),$$

where $v\_{x}\left(0\right)$ is the $x$-component of Mars' velocity vector wrt the Sun (obtained from JPL), and $a\_{x}\left(0\right)$ is the $x$-component of Mars' acceleration vector solely due to the Sun's gravitational attraction. The corrected position of Mars at $t\left(1\right)$ is then given by

$$x\_{MARS}\left(1\right)=x\_{MARS}\left(0\right)+\epsilon v\_{x}\left(\epsilon/2\right).$$

I then find the next value of $x\_{MARS}\left(1\right)-x\_{JUPITER}\left(1\right)$ using my previously calculated $x\_{JUPITER}\left(1\right)$ value (from the Jupiter orbit spreadsheet). Next, I repeat the above to find $MJa\_{x}\left(1\right)$, the $x$ component of Mars' acceleration toward Jupiter at $t\left(1\right)$. And so on.

I was rather pleased with all this except for the fact that it doesn't work, i.e. the results are consistently less accurate than the original Sun-Mars model.

Any idea where I'm going wrong? Thanks.

EDIT - I've now amended my model by using $x\_{MARS}-x\_{JUPITER}$ values calculated directly from my existing Mars and Jupiter orbit spreadsheets (ie simply subtract $x\_{JUPITER}$ from $x\_{MARS}$. With those positions known, I then don't need to worry about calculating velocities as I can then use $x\_{MARS}-x\_{JUPITER}$ to directly find the accelerations $MJa\_{x}$, etc. I've also taken into account the gravitational effect of Jupiter on the Sun as implied in the above $\frac{dv\_{ix}}{dt}$ equation. The results are ambiguous at best. Am I doing this numerical integration correctly?

**Tags:** orbital-elements, n-body-simulations, ephemeris, modeling

> **Comments:**
> **username_1**: Your GM values look a little bit off, but I guess they're ok if they give acceptable results when you compute the orbits of Mars & Jupiter separately. The GM values from DE440 for the planets are given on https://ssd.jpl.nasa.gov/astro_par.html The value for the Sun is 132712440041.279419.  
> **username_1**: Your initial vector data should be for the planet barycenters, not the planets. That is, body IDs 3 & 4, not 399 & 499. That won't matter much for Mars, but it does for Jupiter. And maybe you *do* need to consider the force of Jupiter on the Sun, since the Sun-Jupiter barycenter distance is roughly equal to the Sun's radius.  
> **username_0**: @username_1 - Yes, I assumed, incorrectly, that the Sun $GM =4\pi^{2}$, which I've now corrected. I tried using Jupiter barycentre values but with no obvious increase in accuracy.  
> **username_2**: Hi, what do mean you when you say that the results are consistently less accurate? Are you talking about a few percent, or is it totally off? Could you show some plots?  
> **username_2**: Also, I assume the missing square in your seventh equation is a typo ;)  
> **username_3**: I admire you setting this stuff up in a spreadheet, hard work! You wrote that you don't code, but doing that's coding big-time, and spreadsheet language is so awkward! The problem: no mention of barycentric coordinates, so I assume you're working heliocentrically. Did you make sure to pair out the direct and indirect attractions? I.e., perturbing accel. of Mars due to Jupiter = accel. of Mars due to Jupiter MINUS accel. of Sun due to Jupiter. Omit the 2nd term, you get wrong results. If working barycentrically, the accels of the sun due to planets define its barycentric motion. Good luck!  
> **username_3**: Another point: are you setting up checksums to see what's going wrong? Examples, at each step, vector sum of every paired  interaction (EITHER: force of A on B + force of B on A, OR mass A x accel of A due to B + mass B x accel of B due to A) should be zero, if not there's an error.  Sum of all the orbital angular momentum vectors should be a constant vector. Total energy should sum to a constant. If working barycentrically, the center of mass of the whole system should stay constant. Hope that helps.  

---

## 1 answers

### Answer by username_2 on 2022-07-26. Score: 3 | Accepted answer

I don't have enough information to determine what is wrong with your spreadsheet, but I wrote the same program in Fortran, trying to follow your description. My results seem to be quite accurate, considering that Leapfrog is only second order and that we are neglecting Saturn, the Earth and the other planets.
You can try and compare my results with yours, to understand what is the problem.

I chose to use Fortran because it is very readable, almost like writing math formulae directly, therefore you shouldn't have problems comparing my formulae with yours. Lines starting with an exclamation mark (!) are comments.

The program is the following.

```
program ThreeBody
implicit none
! UNITS used:
! Length: AU
! Time: Julian year

! Declare all the variables
integer, parameter:: dp=kind(0.d0)
! 3D position vector of Mars, the Sun and Jupiter
real(dp) :: r_m(3), r_s(3), r_j(3)
! 3D velocity vector of Mars, the Sun and Jupiter
real(dp) :: v_m(3), v_s(3), v_j(3)
! 3D acceleration vector of Mars, the Sun and Jupiter
real(dp) :: a_m(3), a_s(3), a_j(3)
! Mass of Jupiter, mass of the Sun
real(dp) :: GM_j = 0.037685, GM_s = 39.478418  ! au^3 yr^-2 
! Integration time step
real(dp) :: eps = 0.000684 ! 6 hours in Julian years
! index of current step, unit of the output file
integer :: i, outputFile

! Initial conditions from JPL Horizons 2000-Jan-01 12:00:00.0000 TDB
! Mars
r_m = [1.383579466402647E+00, 1.621204117825275E-02, -3.426152579691657E-02] ! au
v_m = [6.768779495222499E-04, 1.517984118259482E-02, 3.015574279972501E-04]  ! au/day
v_m = v_m*365.25 ! au/yr

! Jupiter
r_j = [3.994040256427818E+00, 2.935779784974931E+00, -1.015789849456240E-01]
v_j = [-4.562948121087195E-03, 6.435847855854325E-03, 7.548691257292337E-05]
v_j = v_j*365.25

! The Sun
r_s = [-7.137179161607904E-03, -2.795997495567612E-03, 2.062985061910257E-04]
v_s = [5.378460339618783E-06, -7.406916207973516E-06, -9.434293432137470E-08]
v_s = v_s*365.25

! open output file
open(newunit=outputFile, file="output.csv", status="replace")
! Write header to file
write(outputFile, '(A)', advance="no") "time,"
write(outputFile, '(A)', advance="no") "rx_m,ry_m,rz_m,rx_j,ry_j,rz_j,rx_s,ry_s,rz_s,"
write(outputFile, '(A)', advance="no") "vx_m,vy_m,vz_m,vx_j,vy_j,vz_j,vx_s,vy_s,vz_s,"
write(outputFile, '(A)')               "ax_m,ay_m,az_m,ax_j,ay_j,az_j,ax_s,ay_s,az_s"

! Main loop. Here we repeat the Leapfrog step 10000 times
do i = 1, 10000
    ! Calculate the acceleration
    ! Ignoring the effect of Mars on the Sun and on Jupiter
    a_m = -GM_j * (r_m - r_j)/norm2(r_m - r_j)**3 - GM_s * (r_m - r_s)/norm2(r_m - r_s)**3
    a_j = -GM_s * (r_j - r_s)/norm2(r_j - r_s)**3
    a_s = -GM_j * (r_s - r_j)/norm2(r_s - r_j)**3

    ! Calculate the velocity at the 1/2 step
    if (i == 1) then
        ! If we are at the first step, use eps/2 instead of eps
        v_m = v_m + eps/2 * a_m
        v_j = v_j + eps/2 * a_j
        v_s = v_s + eps/2 * a_s
    else
        v_m = v_m + eps * a_m
        v_j = v_j + eps * a_j
        v_s = v_s + eps * a_s
    end if
    
    ! Calculate the position
    r_m = r_m + eps * v_m
    r_j = r_j + eps * v_j
    r_s = r_s + eps * v_s

    ! Write the output to a csv file
    write(outputFile, '(*(G0.6,:,","))') i*eps, r_m, r_j, r_s, v_m, v_j, v_s, a_m, a_j, a_s
    
end do

close(outputFile)

end program

```

To run it, you can copy and paste it into a text file called `ThreeBody.f`. If you are on Linux, you just have to open a terminal in the directory where you have saved the file and run

```
gfortran ThreeBody.f && ./a.out

```

If you are on Windows or Mac, you will probably first need to [install gfortran](https://fortran-lang.org/learn/os_setup/install_gfortran) or any other compiler of your choice.

**Alternatively, you can just [Try it online!](https://tio.run/##nZZbr9M4EIDf@ytMX2hLGnyJ7bgSD7CgBQSrlThoHxAceRu3J900iZw2UP782ZnEbdOyLIiqSuvJ5Jvx3OJV5XfelvP1qv9zf1/7au3tltzceeeeVdlhROCTb@siX@Y7Ulal6yQPyPs/Xt28I/vGZYsgeePK9e5uQZ6@D4KbfOsW5OCsHwXJc7csrHfEFgXZ3TnSWp/bvwvX9GbKnVs7H5Hagg9u5/xiQbL6yT95mU1onNFpp@adLSZZPR3oEVCsc/KEiJglTBqupJCp0UZwkSaKq0SI4IJ4TuqqyXd5VZLWLXeVJ9WKvLW@iTqX3u1LYsuMvN7XOZAvLKIZf7udCLDtb5vwu4HfM7x1RQWhOvwSvA3wNsDbK7hdLl3hvP1l720wYIMBOzTw1jYN4sLDEdkGQSB/Q/v97e0Ggk5jKrRKZYSCBrNg4kSnCUuRavefBDn4T3NOgp1XmOewiR3UCGl2rv4G7uqmZ1Oq0gQeU@Su2vsGygQ8LHJbdqXVBCjUiPuCzi733rtyd4aGskJmPjq5ABVgC7KsyqyrhYasfLUlr/98Q15WPv@KEg6m56@hPygjjC8ohS@6Q8nN82enmAUPoC7A3Q8sFqmQ2iRKJRTqTr94RGlEWKw44zRhTKdcci1fzCmPyFzECVdMcnhCGaakRvnHLmodte2pKob4aqAayTlPjAG1BKmSaQOBZgA1Scp7qogpk1InXBujuaQMtT/2uXic2cOADNeZUDLmsr97OLXqRQV1af4gYmMSCtuSCtgpS/u98djAjrWBjBudGMF68ZyhGzoF/0wiFU8ouME@BuMdcJ7EcAMcZ5zRVDODYRERUXEC7QslJGUqE8GDWMcySSFMEC1uuBAYLDkEwjXsphMex9Cwerv6/DDXMRMaYscUU1QbmvQW5jzWRhrchpRKK8Z7OY8p@JlKqphhsH/dRTRY7ogyFrB/RYUwiqU6FaChAKjjhGJmOVjRQjIV5NAhIuFGwO7Ak0RjcNIhEa7HvRxT8peHjJA7Z7OQl88omMwi8nDydPoQujlrbbl0T8ZlNZ6SMfZWNP4ZTf/ldhv5A16@4uXL7QaXG1xucNngssFl83PEFoktElsktkhskdgisUVii8T2e8QpufyMLRItEi0SLRItEi0SLRItEmHZjEen9oRpUVRVHZOXDl48nx3MmNrZXTfT3jhbQ9uvu1lBWNfaGLS@o7OK4CuFRXhndPTjAfnNFst9YSEVyBiO5IHSq3UJY6RcdzputYJZfZzUBMfeYFZX5UWz4cd2zTnv5uuMTHC2zLELp4@BuuUDwWwm4F83dweKzbVig4oD/Cbgw1Ob66c2//lUc@lUc@1Uc3Zq9N2AnV6QIQnsMT/P6m5er8gEAg@Rn6LCOaohsivMYneK6AGr3Df9uI/wPIIvDkDmJUhshkEHwQXjNPrIo6A8w4hf6YSBcqGzudIJjXqhc7blisb9n@Efm/2x0WuTUE/56rT8bhKOR6CTQv/68gPX2oFr/QvAD9xqB271Q9UPXAL3hgXQDy60W@139R7SVhFLlk0LuSvOIRoMgNnkdxqraBGNo/EUh0E@AzKetrbdkas7f@H5aNsdkroTE55mtt2RJrqIyegYmKwaXazDYff@/l8)**

The output of the program is a Comma Separated Value (CSV) file called `output.csv` that you can read as a spreadsheet with Excel or LibreOffice. It contains 28 columns:

* 1 time [yr]
* 2-4 3D position of Mars [au]
* 5-7 3D position of Jupiter [au]
* 8-10 3D position of the Sun [au]
* 11-19 same thing, but for the velocity [au/yr]
* 20-28 same thing, but for the acceleration [au/yr^2]

Just for reference, here is the 3D position of Mars for the first 10 steps of integration

[![first 10 steps of integration](https://i.sstatic.net/2IzDh.png)](https://i.sstatic.net/2IzDh.png)

And here it is after about 1 Julian year

[![Mars position after 1 year](https://i.sstatic.net/BUXnc.png)](https://i.sstatic.net/BUXnc.png)

### Remarks on the accuracy of the result:

I don't know what accuracy you were expecting. In the comments I read that you considered $GM = 4\pi^2$ to be incorrect for your required level of accuracy. It best measured value can be found in the Wikipedia page for the [Standard gravitational parameter](https://en.wikipedia.org/wiki/Standard_gravitational_parameter), and, when expressed in units of [Julian years](https://en.wikipedia.org/wiki/Julian_year_%28astronomy%29) and [AU](https://en.wikipedia.org/wiki/Astronomical_unit), it is equal to $4\pi^2$ to at least one part per 10000.

If your required level of accuracy is smaller than this, then you have to be careful that 6 hours are not $0.000684$ Julian years, but (with some more digits) $0.00068446269$ years. A better value for the gravitational parameter of the Sun would be $GM\_\odot = 39.476926414 au^3 yr^{-2}$, and also the one of Jupiter could gain some more digits $GM\_J = 0.037684447 au^3 yr^{-2}$. But, in any case, if you really need to know the position of Mars so precisely, then you probably need to include at least Saturn, the Earth and probably Venus in your calculation. They could have a far greater influence on the result, than the 6th decimal digit of the gravitational parameter.

> **Comments:**
> **username_4**: +1! You had me at "_but I wrote the same program in Fortran_." :D  
> **username_0**: Thanks. A very detailed response that I'm still trying to get my head around. I also can't quite believe that you were able to write an equivalent program in hours (maybe less?) when it's taken me weeks/months to put together my own dodgy spreadsheet :-) I really do need to learn some coding.  
> **username_2**: @username_0 That's the advantage of writing code in a programming language versus a spreadsheet I suppose :D. I also did it because you say you are new to coding, and starting can be difficult if you don't have a guide. I encourage you to take this template and extend it to suit your needs. Please feel free to ask if you don't understand some of the code or have trouble compiling/running  
> **username_0**: I notice you've used the solar system barycentre as coordinate centre instead of (as I have) the Sun (body center). That's more accurate, I assume?  
> **username_0**: Also, your intitial Mars y values are a little off compared to the JPL values?  
> **username_2**: @username_0 The coordinates centered on the Sun are not an inertial system, therefore I would have had to take into account apparent forces. I chose to simplify my life and go for the barycenter instead :) I don't know why you see the y values off compared to JPL. It might be because I used the planet center (499) instead of barycenter (4). If that's the case, it's a sloppiness on my part, thank you for pointing it out.  
> **username_0**: My head is spinning with all this but I think what I've done is reversed the sign of your a_s (acceleration of Sun towards Jupiter) and considered it as an additional acceleration of Jupiter towards the Sun. The total acceleration of Mars then equals the acceleration of Mars towards the Sun, plus acceleration of Mars towards Jupiter, minus your a_s. I hope I've got that right - it seems to work, anyway. I've also just noticed that your initial y value for Mars should be negative. I've posted an “answer” listing some of my results. Many thanks.  

---

