2

I'm modelling planetary orbits using a spreadsheet and Feynman's numerical method given here. I have put together a spreadsheet (10,000 rows) to model Mars' orbit in three $\left(x,y,z\right)$ dimensions. I'm using SI units with a time interval of 3600 seconds. I'm using initial state vector data from JPL Horizons. I believe Feynman's method is called leapfrog integration. I've also tried the Runge-Kutta method, but with no significant improvement in accuracy.

In terms of the orbit cycle, does it matter when I take the initial JPL Horizons data? Aphelion or perihelion, for example, or some other point in the orbit?

EDIT. Just to be clear what I'm asking. The first row of my spreadsheet uses state vector data $\left(\mathbf{x},\mathbf{y},\mathbf{z},\mathbf{v_{x}},\mathbf{v_{y}},\mathbf{v_{z}}\right)$ from JPL Horizons for an arbitrary date and time. I've chosen 2000-Jan-01 12:00:00.0000 TDB. Will I get more accurate results if I use a date and time for a specific point in Mars' orbit? Aphelion or perihelion, for example? If the answer is yes, an explanation would be useful.

Peter
  • 1,033
  • 6
  • 17
  • 1
    Don't use SI units! They cause all sorts of floating point issues. AU, Solar mass, year are much better set of units for solar system modelling than metre, kg, second. – James K Jul 18 '22 at 07:21
  • As I note in https://astronomy.stackexchange.com/questions/2416/where-can-i-find-a-set-of-data-of-the-initial-conditions-of-our-solar-system/13491#13491 you might want to use the same initial conditions JPL does. – Barry Carter Jul 18 '22 at 13:38
  • Yes, that's Leapfrog. There are a few variants of Leapfrog. I illustrate a couple of them here: https://astronomy.stackexchange.com/a/48477/16685 – PM 2Ring Jul 18 '22 at 13:49
  • Numerical integration can use any starting point. Analytical orbit calculations tend to use perihelion, because that's where the mean, eccentric, and true anomalies are measured from. However, there are practical issues in determining the exact moment of Earth's perihelion, as I mention here: https://astronomy.stackexchange.com/a/49605/16685 – PM 2Ring Jul 18 '22 at 13:53
  • @PM2Ring - I was surprised that the Runge–Kutta method wasn't (for me) much of an improvement on leapfrog. It took me ages to realise that comparing the radius of my predicted position and JPL's isn't very useful. If the Earth was a perfect sphere and Tom was standing on a table next to me our radii would be slightly different. If Tom travelled to the opposite side of the planet our radii would be the same though our positions on the Earth's surface would be very different. – Peter Jul 18 '22 at 15:51
  • Leapfrog is symplectic, Runge-Kutta isn't. https://en.wikipedia.org/wiki/Symplectic_integrator That gives Leapfrog an advantage for simple orbit integrations like this: your orbit conserves energy, so it won't spiral in or out due to accumulated integration errors. 4th order Runge-Kutta can be better than plain Leapfrog, especially with highly eccentric orbits, but 4th order Leapfrog is even better. ;) – PM 2Ring Jul 18 '22 at 16:11
  • Many useful suggestions in the comments :) but I am kind of missing the point of the question. What is it asking exactly? – Prallax Jul 18 '22 at 20:12
  • 1
    @Prallax - I've edited my question to make it clearer what I'm asking. – Peter Jul 19 '22 at 06:07
  • 1
    If I were a pathological nitpicker, I'd point out that coordinates not from the "epoch" are given as polynomial approximations to the real solution and are therefore trivially less "accurate". Of course this "inaccuracy" is vastly smaller than the inaccuracy of the positions themselves. – Barry Carter Jul 19 '22 at 15:07

1 Answers1

3

You can begin with any starting point. The numerical integration process doesn't "know" anything about ellipses, perhelion, or what the orbit is supposed to look like. The same iterative process is applied at all points on the orbit, so you can start anywhere.

It is a good idea to use astronomical units (AU, Solar mass, years) rather than SI units. With SI units, the intermediate calculations can involve very large and very small numbers and this can result in a loss of accuracy.

A symplectic integrator should do well enough, but note that the NASA values include all kinds of perturbations (by other planets, by the non-spherical sun etc) and a two-body model isn't going to match perfectly. Also remember that the sun is also moving, so if your x-y-z are relative to the centre of the sun, that is an issue. To deal with it, you would need to include other planets and the sun in your model.

For speed you probably will want to move to a "real" programming language. Spreadsheets are optimised for business calculations, rather than solving differential equations.

James K
  • 120,702
  • 5
  • 298
  • 423
  • Thanks for that. JPL offers the following choice of state vector units: km and seconds; AU and days; km and days. However, they give GM for the Sun in units of $\mathrm{km}^{3}\mathrm{s}^{-2}$. I've tried using this value of GM in the spreadsheet with distances in km, but the results come out exactly the same as with using metres and seconds. I haven't been able to find a figure for GM in AU and some other time units (days or years?). – Peter Jul 19 '22 at 07:57
  • 2
    in units of "sun mass, AU and years" G = GM = 4pi^2 – James K Jul 19 '22 at 08:12
  • Well, I never saw that coming. What an amazing result. – Peter Jul 19 '22 at 09:36
  • 1
    @Peter What James said about GM is true, but it depends on how you define "AU" & "year". See https://en.wikipedia.org/wiki/Gaussian_year & links therein. – PM 2Ring Jul 23 '22 at 18:23
  • Yes, there is a slight of hand, If you define "Year" as the time that a body in circular 1 AU orbit around a mass of the sun, you can see how the various units are related, and so the G=4pi^2 result becomes somewhat less surprising. Astronomers conventionally take 1 year = 365.25 days (exactly) for which G is only very close to 4pi^2 – James K Jul 23 '22 at 19:01