Summary: What's an efficient way of finding planetary conjuctions with Mathematica using the SPICE kernels, and without using AstronomicalData (which I believe is less efficient).

NASA publishes planetary positions (ICRF reference frame) as polynomial interpolations for the past and future 15,000 years or so:



  • Each component of Mercury's position (x,y,z) is described by a 14th degree polynomial that changes every 8 days. The first 13 derivatives are continuous at the change points, as one would expect from an interpolation.

  • Each component of Jupiter's position is described by an 8th degree polynomial that changes every 32 days (with 7 continuous derivatives).

  • In other words, the interpolation's polynomial degree and interval change with each planet.

More details for those interested:



(despite the header, the document above also describes the current DE430 format).

I've created Mathematica versions of all of these files (but only for the visible planets), and uploaded some of them to (ran out of room on dropbox; if anyone has space they want to donate for the rest, please let me know):


These files are fairly large (about 350M each), so I found a file with data for only 100 years (1950 to 2050, Julian days 2433264.5 to 2469808.5):


The Mathematica version is 38M:


[Note: I realize that applying N to the coefficients would speed things up, but even that is slow]

For reference, these files were created using:


The dump files define the variable pos, which contains the raw NASA coefficients in a slightly more useful format.

You can get from pos to the actual x,y,z position of a planet at a given time (in the ICRS J2000 reference frame) with this code:

(* A planets position *)

posxyz[jd_,planet_] := Module[{jd2,chunk,days,t},

   (* normalize to boundary *)
   jd2 = jd-33/2;

   (* days in a given chunk *)
   days = 32/info[planet][chunks];

   (* which chunk *)
   chunk = Floor[Mod[jd2,32]/days]+1;

   (* where in chunk *)
   t = Mod[jd2,days]/days*2-1;

   (* and Chebyshev *)

(* Chebyshev of a list at a variable *)

chebyshev[list_,t_] := Sum[list[[i]]*ChebyshevT[i-1,t],{i,1,Length[list]}]

(* Definition of info, not all lines are relevant to this question *)
info[earthmoon][chunks] = 2
info[earthmoon][name] = earthmoon
info[earthmoon][num] = 13
info[earthmoon][pos] = 231
info[jupiter][chunks] = 1
info[jupiter][name] = jupiter
info[jupiter][num] = 8
info[jupiter][pos] = 342
info[mars][chunks] = 1
info[mars][name] = mars
info[mars][num] = 11
info[mars][pos] = 309
info[mercury][chunks] = 4
info[mercury][name] = mercury
info[mercury][num] = 14
info[mercury][pos] = 3
info[moongeo][chunks] = 8
info[moongeo][name] = moongeo
info[moongeo][num] = 13
info[moongeo][pos] = 441
info[neptune][chunks] = 1
info[neptune][name] = neptune
info[neptune][num] = 6
info[neptune][pos] = 405
info[pluto][chunks] = 1
info[pluto][name] = pluto
info[pluto][num] = 6
info[pluto][pos] = 423
info[saturn][chunks] = 1
info[saturn][name] = saturn
info[saturn][num] = 7
info[saturn][pos] = 366
info[sun][chunks] = 2
info[sun][name] = sun
info[sun][num] = 11
info[sun][pos] = 753
info[uranus][chunks] = 1
info[uranus][name] = uranus
info[uranus][num] = 6
info[uranus][pos] = 387
info[venus][chunks] = 2
info[venus][name] = venus
info[venus][num] = 10
info[venus][pos] = 171
info[jend] = 2.4698085*^6
info[jstart] = 2.4332645*^6

This gives the coordinates from the solar system barycenter. Converting to Earth coordinates is fairly easy, as if finding the angular separation between two planets as viewed from earth:

(* the vector between earth and a planet *)
earthvector[jd_,planet_] := posxyz[jd,planet]-posxyz[jd,earthmoon];

(* angle between two planets, as viewed from earth *)
earthangle[jd_,p1_,p2_] :=  VectorAngle[earthvector[jd,p1],earthvector[jd,p2]];

and there are other helper functions defined in:


(earthmoon is actually the position of the Earth-Moon barycenter, which is about 1000 miles below the surface of the Earth, and thus a close enough approximation)

Note that earthangle is measuring the angle of two vectors, both of whom's components are polynomials.

My goal: for every pair of planets, efficiently find each minimal angular separation that is smaller that 6 degrees (but this number may change).

In other words, efficiently find local minima of earthangle[jd,p1,p2] for every pair {p1,p2}.

The function earthangle is not pretty. Here it is for Mercury and Venus for a two year period (x axis is Julian dates, y axis is angular separation in radians):

enter image description here

As the image shows, there are many minima, some less than 6 degrees, others more than 6 degrees.

Since I'm doing this for a 30,000 year period (although only 1,000 years at a time due to memory limitations) for every pair of visible planets (15 combinations total), efficiency is important.

Perhaps because of the odd way I define my functions, Plot yielded these errors when constructing the plot above:

In[11]:= Plot[earthangle[jd,mercury,venus],{jd,2457267,2457267+720}]

Thread::tdlen: Objects of unequal length in
                             49993882064   290654274823  8330274361
    0.943016 + {>} + {-(-----------), ------------, ----------,
                                68169         21727        581364
        12136582411     1832669651     13799713729     78664514
      -(-----------), -(----------), -(-----------), -(---------), >,
          457477         8824344        982691536      784737751
         70505        25484          2723            2911
      -----------, ------------, -------------, --------------} cannot be
      84669807291  541431163313  2431800831074  12152414738889

Thread::tdlen: Objects of unequal length in
                             747975063017   213754948321  34683805773
    0.943016 + {>} + {-(------------), ------------, -----------,
                                12367          216355        47894
      4599643706    1091110147     1417781194     74049527
      ----------, -(----------), -(----------), -(---------), >,
       2681989       3856746       267621909      118474419
           129019          27629          3831           3761
      -(------------), -------------, -------------, -------------} cannot be
        922677826651   2127178466734  6849055148965  5827735896081

Thread::tdlen: Objects of unequal length in
                             4547634418852     87664785731   9636715468
    0.943016 + {>} + {-(-------------), -(-----------), ----------,
                                141179           102014        25007
      4089570547    1158179044     200361409     58983665
      ----------, -(----------), -(---------), -(---------), >,
       1115415       8936817       145805420     182336278
           30261           6807           3252            4155
      -(------------), -------------, -------------, --------------} cannot be
        187421371331   3200555072527  5932331907649  77179719891904

General::stop: Further output of Thread::tdlen
     will be suppressed during this calculation.

                   Pi (-(--) + Re[jd])
                         2                       Pi Im[jd]
Plot::exclul: {Sin[-------------------] - 0, Sin[---------] - 0,
                           32                       32
                     33                                   33
         Pi Re[Mod[-(--) + jd, 32]]           Pi Im[Mod[-(--) + jd, 32]]
                     2                                    2
     Sin[--------------------------] - 0, Sin[--------------------------] - 0,
                     8                                    8
         Pi >                   Pi Im[jd]
     Sin[--------] - 0, >, Sin[---------] - 0, Im[>] - 0} must be a
            8                          16
     list of equalities or real-valued functions.

Out[11]= -Graphics-

This isn't a huge deal for Plot (especially since it actually produces a plot), but functions like (N)Minimize, FindMinimum, etc, choke on my functions and won't produce any results.

I suspect that Plot is trying to take symbolic derivatives of my functions, which doesn't work in this case.

My current approach is purely iterative:

  • Compute angular separations of each pair of planets daily.

  • Find minimal elements in the list.

  • Use the ternary method to find the intraday minimum separation.

I actually did this, and it worked, but was very slow, and I'm sure there's a better way.

Long-term, I'm also trying to find minimum angular separations for 3 or more planets. In this case, angular separation is defined as the maximum separation between any two planets in the grouping. However, this gets harder (and introduces Max into the equations), so I'm holding off on that for now.

Also long-term, I'm trying to find minimum angular separations between planets and fixed stars, although, in this case, I'm only interested in separations within 3 degrees. In the ICRF frame, stars are fixed vectors, so this shouldn't be that hard(?), but there are over 300 such stars, so I could be wrong.

I've asked several questions relating to this, but have been frustrated because it's hard to explain just part of what I'm doing. I'm hoping there is enough general interest here not to close this question.

For linkage purposes, this question relates to these other questions I've asked:

Techniques to find all local minima of black box function with n continuous derivatives?

What Method values are available for Plot?

Computing planet conjunctions with 2D circular orbits still hard?

Find *all* numerical solutions to cosine based equality

EDIT (to answer @Mark_Adler's concerns): my goal here is to use Mathematica's advanced numerical (and non-numerical) techniques to find a non-iterative solution.

For example, the angle between two planets can be computed using arccos, and the minimum separation can be found by setting the derivative to 0. Since the derivative of arccos isn't a trignometric function, and all vectors have polynomial components, we are ultimately solving a polynomial roots problem (I realize that the product of the Norm of the planetary distances [necessary to compute the arccos] is actually the square root of a polynomial, so I'm really saying the problem we're solving is "polynomialish")

I realize that solving the polynomial equation may take more time than the iterative method, but I still think it would be interesting to find and use, in part because it can be used on a broader range of problems.

  • 3
    $\begingroup$ The problem looks interesting, but I believe it's a bad match for this site. Perhaps someone could suggest another fora where you can share projects this size. $\endgroup$ Sep 1, 2015 at 16:54
  • 2
    $\begingroup$ I agree with @belisarius. The question will probably be closed as "too broad" unless you boil it down to much more specific issues. Some related issues have of course been raised in your earlier questions, and maybe you could just provide more feedback to the people who responded there, to elicit more answers if what you learned so far doesn't satisfy your needs. $\endgroup$
    – Jens
    Sep 1, 2015 at 16:59
  • 1
    $\begingroup$ It's a bit complicated, yes. But I can offer a tiny tip now: considering the degree of the polynomials you're dealing with, I would suggest looking into the Clenshaw recurrence to evaluate your Chebyshev series. $\endgroup$ Sep 1, 2015 at 17:24
  • 2
    $\begingroup$ It is very difficult to comment on code that is not shown. A guess is that you will need to restrict some functions to only operate on explicitly numeric input, e.g. earthangle[jd_?NumericQ,...] := ... $\endgroup$ Sep 1, 2015 at 17:37
  • $\begingroup$ I don't think you have clarified the question by making it even longer and introducing the undefined term "polynomialish." $\endgroup$
    – Jens
    Sep 3, 2015 at 18:39

1 Answer 1


You should use the SPICE toolkit (in C), which you can access from Mathematica with LibraryLink. The toolkit routines will handle all of the kernel access and coordinate transformations for you transparently and fast. In my opinion it would be both a waste of time and highly error prone to try to duplicate their functionality. Not to mention slower.

The SPICE routine gfsep_c() does exactly what you want:

Determine time intervals when the angular separation between the position vectors of two target bodies relative to an observer satisfies a numerical relationship.

Replying to the edited question: No.

These are vectors of polynomials of degree 5 to 13, which you subtract and dot product to get even higher degree polynomials, and normalize with the square root of more dot products on the bottom. The numerator of the derivative of the normalized dot product of the Earth-Venus and Earth-Mercury vectors is a polynomial of degree 63. Every eight days you have a different such polynomial. Most won't even have a zero over that range. The ones that do will require numeric root finding. And then you still haven't found where the normalized dot product is greater than $\cos\left(6^\circ\right)$, so you will need need to numerically search for those, which may take you to adjacent eight-day slices.

SPICE has efficient condition search algorithms, so you're probably not going to be able to implement anything significantly better. I just tried gfsep_c() on the DE431 kernel, and found about 64,000 Mercury/Venus conjunctions over 30,000 years, identifying the times going below and above 6° to an accuracy of one second, all in two minutes. With one function call.

  • $\begingroup$ Point taken, but note the parameter step in that routine: "Step size in seconds for finding angular separation events". In other words, it still appears to be an iterative search. $\endgroup$
    – user1722
    Sep 2, 2015 at 2:39
  • 2
    $\begingroup$ How can it not be? $\endgroup$
    – Mark Adler
    Sep 2, 2015 at 4:03
  • $\begingroup$ Please see edit to original question, thanks! $\endgroup$
    – user1722
    Sep 3, 2015 at 17:00