Estimating close approach between all asteroids themselves

I am currently writing a C++ program to show asteroids in 3D, and find close approaches or collisions. I got my orbital elements from JPL https://ssd.jpl.nasa.gov/sbdb_query.cgi. So far so good, with over a million asteroids drawn in 5ms with my old Quadro GPU

But the problem is the variation of the elements over time. For example, CERES:

epoch,a,e,i,om,w,ma
2459200.5,2.766089105818,.07816842657453,10.58789954719,80.27235841368,73.72488984426,205.5454154582
2459310.5,2.765760313090,.07831877879848,10.58807660401,80.26860808947,73.73699886586,229.1146825391


After only 4 months, the eccentricity and all other elements have changed as seen on the second line.

How can I compute the change in the parameters without knowing their derivatives?

Or where to get those derivatives?

I searched all JPL small bodies site and their telnet access or email, but could not find a way to download orbital elements with derivatives or a way to compute the change.

I am aware that the difference might be just a pixel on the screen but it represents hundreds of thousands of kilometres. What I want to do is study the close approach between the asteroids themselves and eventually near collision. I already did that for artificial satellites using SGP4 propagator using OpenCL running on GPU video card. I can propagate 20,000 satellites (include debris) in a few milliseconds and get results exact to 1 kilometre. Compare with Celestrak Socrate. A prediction of a 1,000-kilometre approach must be possible.

Does anybody know how?

• few quick questions: 1) if you are drawing a million asteroids will anybody know that a few are a half-pixel off? 2) wait you downloaded orbital elements for a million asteroids from JPL? 3) what are the uncertainties on those osculating orbital elements? 4) if you had all their derivatives, are you confident you would know how to use them to propagate correctly? 5) could there be a reason that osculating elements don't come with derivatives? – uhoh Apr 7 at 0:57
• Stupid answer from me: Grab every single object in the solar system at epoch J2000 and run a huge nbody simulation over it. If you compare that to your given parameters and it is reasonably accurate, keep doing it. But remember that a butterfly flapping its wings can cause a hurricane! – slowerthanstopped Apr 7 at 1:53
• Salut François! The equations describing the precession of nodes and the precession of apsides are quite complex. I doubt you’d want to perform them for millions of bodies! As @uhoh pointed out, nobody will notice that your asteroids are “a half-pixel off” (if that much!). Unless you want to go into the distant past or distant future, I see no use for that. – Pierre Paquette Apr 7 at 1:57
• @FrancoisBilodeau that's exactly the kind of information that you should put in the question to begin with. The more you can share about the type of solution you need, the better that answers can address your specific need. You emphasized fast screen painting in the question and said nothing about accuracy or close approach detection. If that's what you are after, great! But please emphasize that clearly by editing your question post accordingly. Many/most users won't dig down into the comments before writing an answer. Thanks! – uhoh Apr 7 at 18:44
• In this case you really might want to consider switching to propagating state vectors using a gravity model of several solar system bodies, nor simply orbital elements. – uhoh Apr 7 at 18:45

Thanks folks for helping. I knew about SPK files but from Horizon's telnet inteface there is a limit of 200 bodies per request. Making over 5000 request might be possible but not very productive. I also looked at the DE421 file from ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/ but it looks like the corrections for the 8+1 planets

here is the output from BRIEF, a SPICE utility: https://naif.jpl.nasa.gov/naif/utilities.html

BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066
Summary for: de421.bsp
Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
JUPITER BARYCENTER (5)  SUN (10)                MARS (499)
Start of Interval (ET)              End of Interval (ET)
-----------------------------       -----------------------------
1899 JUL 29 00:00:00.000            2053 OCT 09 00:00:00.000


I can always do a statistical estimate, any close approach being false but the overall count being close to reality. Too bad to drop the project, I was getting nices pictures from it. The coloring is done against eccentricity, red=0

• 1. It would be a shame to drop this project. Maybe contact Horizons, explain what you want to do, and see if they have any suggestions. They may be able to supply the SPK files in bulk (depending on what time span you want to cover), or have some other practical suggestions. 2. Since you're investigating near encounters, I guess even the SPK files aren't adequate when two asteroids are close enough to have significant gravitational affect on each other, unless one or both are in the set of 343 asteroids included in the dynamical model. – PM 2Ring Apr 8 at 13:43
• well, thanks PM, you can give it a ride at geomaitre.com/geomaitre You must have a graphic card that support opengl 4.5 most modern one do Good luck – Francois Bilodeau Apr 8 at 15:27

This is a common problem with ephemeris or TLEs. They change over time and one might want to know the values after the times that are given.

Most folks use an interpolation scheme to estimate values at a certain time. I would suggest starting with linear interpolation. If this doesn't give you good enough fidelity, you could move on to a more complicated interpolation scheme. Industry often uses a cubic spline.

For your example, using linear interpolation for the Ceres data, the epoch values are $$t_1 = 2459200.5$$ and $$t_2 = 2459310.5$$ (in days) with semi-major axes $$a_1 = 2.766089105818$$ and $$a_2 = 2.765760313090$$ (in AU). The linear interpolation estimate of the change in semi-major axis in AUs per day is $$\dfrac{da}{dt} \approx \dfrac{a_2-a_1}{t_2-t_1} \approx -2.98902 \times 10^{-6}$$. If you want to know the estimate for the semi-major axis at $$t_2+25$$ days, the answer is $$a_2+25\dfrac{da}{dt} = 2.76568558759$$.

If you have the elements for several times, you can compare the estimates for pairs of times. If you get similar results, then a linear interpolation will be fairly accurate.

• How do "most folks use an interpolation scheme" to actually propagate a Keplerian orbit? Can you show an example or cite a source? I've never seen it done and don't see how at all. "This is how to interpolate" seems to miss the point of the question. – uhoh Apr 7 at 18:48
• You may want to look at ast343de430.bsp in ssd.jpl.nasa.gov/pub/eph/small_bodies/asteroids_de430 and other files in the ssd.jpl.nasa.gov/pub/eph/small_bodies directory, but these will just give you a few thousand more, not the 1 million you want. – Guest Apr 12 at 16:52