5
$\begingroup$

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), 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?

$\endgroup$
7
  • 1
    $\begingroup$ 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 ssd.jpl.nasa.gov/astro_par.html The value for the Sun is 132712440041.279419. $\endgroup$
    – PM 2Ring
    Commented Jul 23, 2022 at 18:36
  • 1
    $\begingroup$ 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. $\endgroup$
    – PM 2Ring
    Commented Jul 23, 2022 at 18:40
  • 1
    $\begingroup$ @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. $\endgroup$
    – Peter
    Commented Jul 25, 2022 at 16:54
  • 1
    $\begingroup$ 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? $\endgroup$
    – Prallax
    Commented Jul 25, 2022 at 20:09
  • 1
    $\begingroup$ Also, I assume the missing square in your seventh equation is a typo ;) $\endgroup$
    – Prallax
    Commented Jul 25, 2022 at 20:22

1 Answer 1

3
+100
$\begingroup$

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 or any other compiler of your choice.

Alternatively, you can just Try it online!

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

And here it is after about 1 Julian year

Mars position after 1 year

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, and, when expressed in units of Julian years and AU, 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.

$\endgroup$
7
  • $\begingroup$ +1! You had me at "but I wrote the same program in Fortran." :D $\endgroup$
    – pela
    Commented Jul 26, 2022 at 11:01
  • $\begingroup$ 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. $\endgroup$
    – Peter
    Commented Jul 26, 2022 at 15:18
  • $\begingroup$ @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 $\endgroup$
    – Prallax
    Commented Jul 26, 2022 at 16:34
  • $\begingroup$ 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? $\endgroup$
    – Peter
    Commented Jul 26, 2022 at 19:43
  • $\begingroup$ Also, your intitial Mars y values are a little off compared to the JPL values? $\endgroup$
    – Peter
    Commented Jul 26, 2022 at 20:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .