<p>After writing <a href="https://astronomy.stackexchange.com/a/49678/7982">this answer to <em>JPL Horizons - &quot;highly accurate measurements of planetary positions&quot; - how do they do it?</em></a> which draws from Park, Folkner, Williams &amp; Boggs 2021 <a href="https://iopscience.iop.org/article/10.3847/1538-3881/abd414" rel="nofollow noreferrer">The JPL Planetary and Lunar Ephemerides DE440 and DE441</a> I'm left wondering how the heck they actually do the numerical integration.</p>
<p><strong>Questions:</strong></p>
<ol>
<li>Do they use a fixed time step for the whole solar system, or is there higher time granularity (whatever that might mean) within say the Jovian moon system than for interactions between planetary barycenters that never get near each other or have rapidly changing relative accelerations.</li>
<li>Are two body interactions calculated hierarchically (e.g. Jovian system as a whole on Saturnian system as a whole), or are all <span class="math-container">$n(n-1)$</span> interactions explicitly, or is it somewhere in-between?</li>
<li>What integrator do they use? Do they use more than one?</li>
<li>Order of magnitude how many time steps do they use per year of simulation? Yes higher order integrators can have bigger steps and several evaluations at different times within the step, but I'm just trying to understand if it's one step per minute or millisecond.</li>
</ol>
<p>Since computers keep getting faster, what was done for the early ephemerides several decades ago may be different than what's done today, or it may not be. I can imagine some tricks to reduce computation time are no longer needed, but there might be hesitancy to switch to simpler algorithms to simply avoid breaking things.</p>

CC BY-SA 4.0
2022-06-26T00:11:40.103Z
Yes there are several parts to this question, but I'm confident that they will all likely be answerable from the same source or two, and/or an answer author familiar with one of the issues will also be familiar with the rest. So in this particular case it's best to keep these together.
uhoh
CC BY-SA 4.0
2022-06-26T01:30:02.870Z
I think it's even more involved than your question implies.  The computations are fit to actual observations, so it's not as simple as setting a time step and letting it run, there's a lot of back and forth between what the known physics says it should be and what was actually observed.  From conversations I've read from people who would likely know, I get the impression that the hardware and software aren't just non-public, they're classified.  If you don't get an answer here, you can always just e-mail JPL directly and ask.  Just remember to update us on the answer.
Greg Miller
CC BY-SA 4.0
2022-06-26T01:43:26.430Z
I'm quite confident this question cannot be answered. JPL intentionally does not reveal the details of their integration. The paper describing the most recent Development Ephemeris says repeatedly that the solar system is "numerically integrated" with no details. Previous releases have given a bit more detail, along the lines of a "variable step-size, variable order Adams family integrator." That's not saying much as there are hundreds of variants of Adams integrators, and getting them to be variable step-size & variable order is more than a bit (pun intended) creepy and kooky.
David Hammen
CC BY-SA 4.0
2022-06-26T01:44:10.977Z
@GregMiller I talk about fitting in my linked answer, and of course essentially all good integrators will provide dynamically updated time step sizes (when I code RK45 I include that function, and [when I integrate the solar system using canned algorithms](https://space.stackexchange.com/a/23409/12102) I'm aware of it as well], but I think that this will turn out to be just about as complicated *as I think it is* and I don't see any implications like that in my question.
uhoh
CC BY-SA 4.0
2022-06-26T01:45:09.557Z
@DavidHammen then I'm confident that it can not be answered *to your satisfaction*, but being the OP my standards will differ from yours. Part of the price of not asking questions is not being able to decide which answers are accepted. So answers of the form "Adams-like" and/or "order of one minute" would do it for me :-) Surely there must be something written about hierarchies and local barycenters somewhere.
uhoh
CC BY-SA 4.0
2022-06-26T01:53:55.210Z
I suspect this lack of publications on the exact technique they are using is intentional. Their competition (there are only two organizations with JPL's stature with regard to solar system ephemerides, the Observatory of Paris's IMCCE (Institute of Celestial Mechanics and Ephemeris Calculations; the acronym is French), and the  Institute of Applied Astronomy of the Russian Academy of Sciences) do much the same. Not releasing details on the technique that lies at the heart of the process of developing an ephemeris means that these three organizations can share observations and compare results.
David Hammen
CC BY-SA 4.0
2022-06-26T01:58:19.583Z
@DavidHammen yes that makes a lot of sense
uhoh
CC BY-SA 4.0
2022-06-26T02:09:30.483Z
@uhoh "One minute step size"?? One might as well use RK4 for this small of a step size. Or maybe even symplectic Euler. One of the key advantages of a linear multi-step integrator (aka Adams family) is that this can enable very large step sizes. I would be very surprised if the step size for Venus orbit (for example) was much smaller than one day per step. I know of single-step / multi-stage integrators (e.g., Runge-Kutta integrators) that take steps that are larger than one orbit. They have lots of internal stages (e.g., RK4 has four stage per step).
David Hammen
CC BY-SA 4.0
2022-06-26T02:12:13.337Z
@DavidHammen you might be "leaking" answers! Since the question is about the DE's and since some Jovian moons for example have periods of only  ̶a̶ ̶f̶e̶w̶ ̶d̶a̶y̶s̶ *7 hours*, I'm starting to get the feeling that there can indeed be different time step sizes in different parts of a hierarchy.
uhoh
CC BY-SA 4.0
2022-06-26T02:15:24.520Z
@uhoh The only moon addressed by the Development Ephemerides is the Earth's Moon. All other planetary systems are modeled as a point mass at the planetary system's barycenter. Modeling the Moon is critical as the Moon's mass is about 0.0123 times the Earth's mass. (That's a handy and rather precise number. There are several zeros after 0.0123 before the next non-zero number arises.) I suppose they should be modeling the Pluto system similarly, but they don't. The Pluto system has very little affect on the orbits of the planets.
David Hammen
CC BY-SA 4.0
2022-06-26T02:21:19.040Z
@DavidHammen Curiouser and curiouser! So while modeling the motion of the major moons of Jupiter and Saturn may be necessary to get from raw delay-Doppler data from spacecraft to those systems' barycenters, those calculations are done separately; *outside of the DE integration?* Now that separation sounds kind-of difficult to me in practice since the Earth's motion is an intimate part of the delay-Doppler data.
uhoh
CC BY-SA 4.0
2022-06-26T02:25:29.070Z
@uhoh Correct. The key concern of a Development Ephemeris is to model the solar system. Modeling a planetary system would get in the way of that goal. The Moon is the one exception because it is only moon of a planet (Pluto is not a planet) that comes anywhere close in mass to that of the primary.
David Hammen
CC BY-SA 4.0
2022-06-26T02:36:22.710Z
@DavidHammen I'm not talking about concerns or displayed output, in this these recent comments I am talking about mathematical necessity. I don't see how they can be separated in practice. For the periods of time that the delay-Doppler data is used, I think they must have the motion of at least some of the moons appear explicitly in the integration. It doesn't mean those trajectories are good for ephemerical purposes nor that they'd be shared, but I don't see how one can leave them out of the calculation itself.
uhoh
CC BY-SA 4.0
2022-06-26T02:40:26.267Z
@DavidHammen I'm not saying what is or isn't so, I'm only saying "...I don't see how..." again the luxury of a question-asker :-)
uhoh
CC BY-SA 4.0
2022-06-26T03:16:40.530Z
@uhoh The error that results from treating the Jupiter system as a point mass located at the Jupiter system barycenter versus a large planet and a bunch of moons on the rest of the solar system is deeply into the noise. This means that a solar system ephemeris model of the planets can get away with treating the Jupiter system as a point mass located at the Jupiter system barycenter.
David Hammen
CC BY-SA 4.0
2022-06-26T03:17:54.683Z
With regard to observing Juno, yes, it's important to model the behavior within the Jupiter system as well as the behavior of the Jupiter system as a whole. The problems are separable, at least with regard to impact on accuracy.
David Hammen
CC BY-SA 4.0
2022-06-26T03:35:47.957Z
@DavidHammen "The problems are separable..." so perhaps off to the side somewhere there can be an integration of a spacecraft in the Jovian system only, and observations of the satellites' positions relative to Jupiter as well as the stat and ending points of Juno's trajectory can all be matched via fitting masses and orbital elements. That gives Juno's position vs time relative to the Jovian barycenter. With that one walks across the hall and gives it to the DE folks who incorporate it into their fitting procedures which includes matching to terrestrial delay-Doppler measurements?
uhoh
CC BY-SA 4.0
2022-06-26T13:19:52.023Z
The great thing about Kepler's theory is that, unlike numerical integration, it doesn't accumulate numerical error. So, all orbit calculation methods I know of make as much use as possible of Keplerian theory. To the extent that they are integrating, they integrate perturbations to Keplerian elements rather than mindlessly integrating the state vector.
John Doty
CC BY-SA 4.0
2022-06-26T15:11:18.740Z
http://ilrs.gsfc.nasa.gov/docs/2014/196C.pdf is an older version of the documentation, and MAY contain more or different detail than the current documentation. You might also look at VSOP2013 which gets similar results, but is more open about how they do things.
Barry Carter
CC BY-SA 4.0
2022-06-26T19:04:58.440Z
@JohnDoty certainly that sounds like a good idea in principle; is there a citable example showing it might actually be used for the calculation of JPL DE's? Actually I'd like to see how that's used in any n>2 system. If you like, I can ask a new question to provide more space for an answer.
uhoh
CC BY-SA 4.0
2022-06-26T19:11:10.733Z
@barrycarter that's a really good point! Each DE release's paper is unique and I wouldn't be surprised if as each improvement develops it's mentioned once or twice but not repeated going forward. I think the page counts of some are much longer than the current DE 440/441 I linked to (15 pages) the DE 430/431 you link to has 81 pages for example, though much of that is details on measured data used in the fitting.
uhoh
CC BY-SA 4.0
2022-06-27T12:30:38.367Z
I actually noticed the page size difference too and even converted to text (also longer), but, after looking at both, I wondered if the page size difference was due to differences in formatting.
Barry Carter
<mathematics><ephemeris><nasa><algorithm>
How the heck does JPL actually do the numerical integration for their development ephemerides? (general questions)
uhoh
# How the heck does JPL actually do the numerical integration for their development ephemerides? (general questions)

(Asked by: username_0 on 2022-06-25)

After writing [this answer to *JPL Horizons - "highly accurate measurements of planetary positions" - how do they do it?*](https://astronomy.stackexchange.com/a/49678/7982) which draws from Park, Folkner, Williams & Boggs 2021 [The JPL Planetary and Lunar Ephemerides DE440 and DE441](https://iopscience.iop.org/article/10.3847/1538-3881/abd414) I'm left wondering how the heck they actually do the numerical integration.

**Questions:**

1. Do they use a fixed time step for the whole solar system, or is there higher time granularity (whatever that might mean) within say the Jovian moon system than for interactions between planetary barycenters that never get near each other or have rapidly changing relative accelerations.
2. Are two body interactions calculated hierarchically (e.g. Jovian system as a whole on Saturnian system as a whole), or are all $n(n-1)$ interactions explicitly, or is it somewhere in-between?
3. What integrator do they use? Do they use more than one?
4. Order of magnitude how many time steps do they use per year of simulation? Yes higher order integrators can have bigger steps and several evaluations at different times within the step, but I'm just trying to understand if it's one step per minute or millisecond.

Since computers keep getting faster, what was done for the early ephemerides several decades ago may be different than what's done today, or it may not be. I can imagine some tricks to reduce computation time are no longer needed, but there might be hesitancy to switch to simpler algorithms to simply avoid breaking things.

**Tags:** mathematics, ephemeris, nasa, algorithm
> **Comments:**
> **username_0**: Yes there are several parts to this question, but I'm confident that they will all likely be answerable from the same source or two, and/or an answer author familiar with one of the issues will also be familiar with the rest. So in this particular case it's best to keep these together.  
> **username_1**: I think it's even more involved than your question implies.  The computations are fit to actual observations, so it's not as simple as setting a time step and letting it run, there's a lot of back and forth between what the known physics says it should be and what was actually observed.  From conversations I've read from people who would likely know, I get the impression that the hardware and software aren't just non-public, they're classified.  If you don't get an answer here, you can always just e-mail JPL directly and ask.  Just remember to update us on the answer.  
> **username_2**: I'm quite confident this question cannot be answered. JPL intentionally does not reveal the details of their integration. The paper describing the most recent Development Ephemeris says repeatedly that the solar system is "numerically integrated" with no details. Previous releases have given a bit more detail, along the lines of a "variable step-size, variable order Adams family integrator." That's not saying much as there are hundreds of variants of Adams integrators, and getting them to be variable step-size & variable order is more than a bit (pun intended) creepy and kooky.  
> **username_0**: @GregMiller I talk about fitting in my linked answer, and of course essentially all good integrators will provide dynamically updated time step sizes (when I code RK45 I include that function, and [when I integrate the solar system using canned algorithms](https://space.stackexchange.com/a/23409/12102) I'm aware of it as well], but I think that this will turn out to be just about as complicated *as I think it is* and I don't see any implications like that in my question.  
> **username_0**: @DavidHammen then I'm confident that it can not be answered *to your satisfaction*, but being the OP my standards will differ from yours. Part of the price of not asking questions is not being able to decide which answers are accepted. So answers of the form "Adams-like" and/or "order of one minute" would do it for me :-) Surely there must be something written about hierarchies and local barycenters somewhere.  
> **username_2**: I suspect this lack of publications on the exact technique they are using is intentional. Their competition (there are only two organizations with JPL's stature with regard to solar system ephemerides, the Observatory of Paris's IMCCE (Institute of Celestial Mechanics and Ephemeris Calculations; the acronym is French), and the  Institute of Applied Astronomy of the Russian Academy of Sciences) do much the same. Not releasing details on the technique that lies at the heart of the process of developing an ephemeris means that these three organizations can share observations and compare results.  
> **username_0**: @DavidHammen yes that makes a lot of sense  
> **username_2**: @uhoh "One minute step size"?? One might as well use RK4 for this small of a step size. Or maybe even symplectic Euler. One of the key advantages of a linear multi-step integrator (aka Adams family) is that this can enable very large step sizes. I would be very surprised if the step size for Venus orbit (for example) was much smaller than one day per step. I know of single-step / multi-stage integrators (e.g., Runge-Kutta integrators) that take steps that are larger than one orbit. They have lots of internal stages (e.g., RK4 has four stage per step).  
> **username_0**: @DavidHammen you might be "leaking" answers! Since the question is about the DE's and since some Jovian moons for example have periods of only  ̶a̶ ̶f̶e̶w̶ ̶d̶a̶y̶s̶ *7 hours*, I'm starting to get the feeling that there can indeed be different time step sizes in different parts of a hierarchy.  
> **username_2**: @uhoh The only moon addressed by the Development Ephemerides is the Earth's Moon. All other planetary systems are modeled as a point mass at the planetary system's barycenter. Modeling the Moon is critical as the Moon's mass is about 0.0123 times the Earth's mass. (That's a handy and rather precise number. There are several zeros after 0.0123 before the next non-zero number arises.) I suppose they should be modeling the Pluto system similarly, but they don't. The Pluto system has very little affect on the orbits of the planets.  
> **username_0**: @DavidHammen Curiouser and curiouser! So while modeling the motion of the major moons of Jupiter and Saturn may be necessary to get from raw delay-Doppler data from spacecraft to those systems' barycenters, those calculations are done separately; *outside of the DE integration?* Now that separation sounds kind-of difficult to me in practice since the Earth's motion is an intimate part of the delay-Doppler data.  
> **username_2**: @uhoh Correct. The key concern of a Development Ephemeris is to model the solar system. Modeling a planetary system would get in the way of that goal. The Moon is the one exception because it is only moon of a planet (Pluto is not a planet) that comes anywhere close in mass to that of the primary.  
> **username_0**: @DavidHammen I'm not talking about concerns or displayed output, in this these recent comments I am talking about mathematical necessity. I don't see how they can be separated in practice. For the periods of time that the delay-Doppler data is used, I think they must have the motion of at least some of the moons appear explicitly in the integration. It doesn't mean those trajectories are good for ephemerical purposes nor that they'd be shared, but I don't see how one can leave them out of the calculation itself.  
> **username_0**: @DavidHammen I'm not saying what is or isn't so, I'm only saying "...I don't see how..." again the luxury of a question-asker :-)  
> **username_2**: @uhoh The error that results from treating the Jupiter system as a point mass located at the Jupiter system barycenter versus a large planet and a bunch of moons on the rest of the solar system is deeply into the noise. This means that a solar system ephemeris model of the planets can get away with treating the Jupiter system as a point mass located at the Jupiter system barycenter.  
> **username_2**: With regard to observing Juno, yes, it's important to model the behavior within the Jupiter system as well as the behavior of the Jupiter system as a whole. The problems are separable, at least with regard to impact on accuracy.  
> **username_0**: @DavidHammen "The problems are separable..." so perhaps off to the side somewhere there can be an integration of a spacecraft in the Jovian system only, and observations of the satellites' positions relative to Jupiter as well as the stat and ending points of Juno's trajectory can all be matched via fitting masses and orbital elements. That gives Juno's position vs time relative to the Jovian barycenter. With that one walks across the hall and gives it to the DE folks who incorporate it into their fitting procedures which includes matching to terrestrial delay-Doppler measurements?  
> **username_3**: The great thing about Kepler's theory is that, unlike numerical integration, it doesn't accumulate numerical error. So, all orbit calculation methods I know of make as much use as possible of Keplerian theory. To the extent that they are integrating, they integrate perturbations to Keplerian elements rather than mindlessly integrating the state vector.  
> **username_4**: http://ilrs.gsfc.nasa.gov/docs/2014/196C.pdf is an older version of the documentation, and MAY contain more or different detail than the current documentation. You might also look at VSOP2013 which gets similar results, but is more open about how they do things.  
> **username_0**: @JohnDoty certainly that sounds like a good idea in principle; is there a citable example showing it might actually be used for the calculation of JPL DE's? Actually I'd like to see how that's used in any n>2 system. If you like, I can ask a new question to provide more space for an answer.  
> **username_0**: @barrycarter that's a really good point! Each DE release's paper is unique and I wouldn't be surprised if as each improvement develops it's mentioned once or twice but not repeated going forward. I think the page counts of some are much longer than the current DE 440/441 I linked to (15 pages) the DE 430/431 you link to has 81 pages for example, though much of that is details on measured data used in the fitting.  
> **username_4**: I actually noticed the page size difference too and even converted to text (also longer), but, after looking at both, I wondered if the page size difference was due to differences in formatting.  

---

## 0 answers

