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:
ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/de431/
Notes/examples:
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:
https://github.com/barrycarter/bcapps/blob/master/ASTRO/format.txt
https://github.com/barrycarter/bcapps/blob/master/ASTRO/header.430_572
(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):
https://www.dropbox.com/sh/i5x7qdxaaer7to6/AADwEtmDpq-z2zGI-Swr3BrQa?dl=0
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):
https://github.com/barrycarter/bcapps/blob/master/ASTRO/ascp1950.430.bz2
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:
https://github.com/barrycarter/bcapps/blob/master/ASTRO/bc-dump-cheb.pl
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 *)
Table[chebyshev[pos[planet][Quotient[jd2,32]*32+33/2][[chunk]][[i]],t],
{i,1,3}]
];
(* 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:
https://github.com/barrycarter/bcapps/blob/master/bclib.m#L359
(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):
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
combined.
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
combined.
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
combined.
General::stop: Further output of Thread::tdlen
will be suppressed during this calculation.
33
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.
earthangle[jd_?NumericQ,...] := ...– Daniel Lichtblau Sep 01 '15 at 17:37