planet coordinate systems (moving observer to another planet) ?

[ifmihai]

I'm new to pyephem. I searched the docs, and it appears it's not possible, at least for now.

would it be possible to move the observer to another planet?

I'm interested right now mostly in calculating declination of venus (for example) as seen from jupiter (for example) Declination => find moments when venus crosses jupiter's equator

The problem is I'm not a physics guy, or an advanced math guy, so I try to find a library to be able to do this

Any help or insight on this?

[brandon-rhodes]

Unfortunately not. The underlying astronomy library beneath PyEphem only works for Earth. If you want to do computations from other planets, try the Skyfield library that I am writing to replace PyEphem. At the moment it can do half of what you need: it can tell you the direction in which Venus lies from Jupiter, through a call like:

ra, dec, distance = jupiter(jd).observe(venus).radec()

Unfortunately this right ascension and declination will be relative to Earth's equator (more or less — more strictly speaking, they will be ICRS coordinates) and not relative to Jupiter's equator. To predict where Venus's position would lie for a Jovian observer, we would have to find a model for Jupiter that predicts where its axis is pointing each day.

What is the end result you are trying to compute, if I may ask? That will let me figure out if Skyfield will soon provide you with a way there!

[ifmihai]

It seems sky field will rock! especially for folks like me, the non-experts. It's the only alternative so far, so I'm glad you mentioned it (I already starred it :P)

I have found a nasa link that can move observer to another planet. But it's cumbersome to get each pair of planets, calculate ephemeris, then glue all ephemeris together into a database, etc.

The end result I'm targeting, is to find the dates (in a given range) when any planet crosses other planet's equator. I want to test and investigate some affirmations made by William Foster about weather forecasting (I don't know how well known he is)

If I don't have a choice, I will grab slowly the ephemeris from nasa horizons, and compile this list of dates manually.

[ghost]

Why "slowly"?

On Sun, 3 Aug 2014, Michael wrote:

Date: Sun, 03 Aug 2014 01:08:36 -0700 From: Michael notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+i-39175826-fa9d326f6ed25bdc01090464a8c62f3a3128a1e0-312349@reply.gi thub.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Subject: Re: [pyephem] planet coordinate systems (moving observer to another planet) ? (#48)

It seems sky field will rock! especially for folks like me, the non-experts. It's the only alternative so far, so I'm glad you mentioned it (I already starred it :P)

I have found a nasa link that can move observer to another planet. But it's cumbersome to get each pair of planets, calculate ephemeris, then glue all ephemeris together into a database, etc.

The end result I'm targeting, is to find the dates (in a given range) when any planet crosses other planet's equator. I want to test and investigate some affirmations made by William Foster about weather forecasting (I don't know how well known he is)

If I don't have a choice, I will grab slowly the ephemeris from nasa horizons, and compile this list of dates manually.

---

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/48#issuecomment-50984715

[ifmihai]

Slowly in the sense that there are 45 pairs of planets? (45,if im correct)

And for each pair i have to generate ephemeris for a range of dates. Then glue all results. Then parse to search days when declination changed signs.

Slow process.

[ghost]

OK. I thought you would download the positions and tilts of each planet and then compute the RA/DEC for each pair, instead of downloading each pair itself.

It's 45 if there are 10 planets, so I'm guessing you're adding 2 to the currently recognized 8? (although Pluto will always be a planet in my heart!). Eris is the other one?

You might want to look at their mail interface and my:

https://github.com/barrycarter/bcapps/blob/master/bc-email-horizons.pl

for some scripting. Or, I'm bored and can script it for you (no charge) if you tell me the data you're looking for?

On Sun, 3 Aug 2014, Michael wrote:

Date: Sun, 03 Aug 2014 08:52:33 -0700 From: Michael notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+i-39175826-fa9d326f6ed25bdc01090464a8c62f3a3128a1e0-312349@reply.gi thub.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Cc: barrycarter github@barrycarter.info Subject: Re: [pyephem] planet coordinate systems (moving observer to another planet) ? (#48)

Slowly in the sense that there are 45 pairs of planets? (45,if im correct)

And for each pair i have to generate ephemeris for a range of dates. Then glue all results. Then parse to search days when declination changed signs.

Slow process.

---

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/48#issuecomment-50993966

[ifmihai]

@barrycarter well, I counted 10 planets, from the sun to pluto (pluto is a planet for me also :P) but indeed, planets crossing sun's equator are easy to calculate => -9 already and also, planets crossing earth's equator are easy again => -9 again so now I think there are 10 planets x 9 others - 9 - 9 = 72 pairs

I was aware of mail interface, but I didn't take it into consideration Now that you point it out, it might be a solution I will have to get into it deeper

well, if you are bored, and want to do something about it :P I guess a range from 1970 to 2014 will do actually make it 1970 to 2020

there are enough hurricanes or major weather reports in this range, i believe

then, the pairs, written as (observer location, target body) (mercury, sun) (mercury, venus) mercury, earth) etc

then (venus, sun) (venus, earth) (venus, mars) etc

then (mars, sun) (mars, mercury) (mars, venus) (mars, earth) etc

until pluto (pluto, sun) (pluto, mercury) (pluto, venus) etc

(but I guess, seeing from pluto, all inner planets are crossing its equator kind of the same time)

meanwhile I will get into mail possibility thank you for you boredom anyway :D

[ghost]

Michael, email or google talk me at carter.barry@gmail.com. This doesn't seem difficult, but I want to check a couple of details with you.

Brandon, if you prefer, we can keep the conversation here, but I wasn't sure it was relevant.

[ghost]

Just for fun, and as part of my never-doing-what-I-say campaign, I fired up stellarium and looked at the night sky from Jupiter.

As I sort of expected, Jupiter's axis, like Earth's is fixed with respect to the stars (ignoring precession and nutation).

Earth's axis, with a tilt of ~23.5 degrees, points (roughly) to Polaris. Jupiter's axis points (very roughly) to HIP 87179:

http://simbad.u-strasbg.fr/simbad/sim-id?Ident=HIP+87179

This star has a declination of 64N 46' on Earth, and 89N 32' on Jupiter. In other words, Jupiter's "north star" is must closer to the north ecliptic pole, since Jupiter's axial tilt of only 1.305 degrees is much smaller than our axial til of ~23.5 degrees.

With a little work, we could refine this information and use it to calculate RA/DEC from other planets. And by "we", I mean "Brandon" :)

[ghost]

It turns out there's a table of "pole positions":

http://en.wikipedia.org/wiki/Poles_of_astronomical_bodies

and individual planet/satellite pages also list this. On:

http://en.wikipedia.org/wiki/Jupiter

it's under "north pole right ascension" and "north pole declination". Assuming these are relatively fixed, combining this with a coordinate transform should solve the problem?

[drbitboy]

Have you consider JPLs NAIF/SPICE? Its Achilles heel for complex tasks is kernel (data) management and it is yet another steep learning curve to climb, but once you have it in your toolbox it is extremely powerful. For simple stuff like the original post in this thread I think there is a one- or few-call solution, and even the kernel management would be pretty easy.

It is what the planetary remote sensing community uses for science and operations, and there is a decent and functional, if informal and slightly incomplete, Python interface (https://github.com/rca/PySPICE).

[brandon-rhodes]

@barrycarter, I was about to ask you about that list of poles! I will see whether Skyfield can use those pole coordinates to determine where in Jupiter's sky another planet is, and then you could extent that solution out to other planets. If I come up with a solution I will post it here.

[drbitboy]

@ifmihai, NAIF/SPICE has a Geometry Finder suite of tools, one of which is for solving the exact problem from your original post:

http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfposc_c.html

The example section finds times when the Sun is in a particular geometry relative to the Earth's body-fixed frame. Replace SUN and EARTH with VENUS and JUPITER in the assignments, and Example 3 is your original request exactly.

A couple dozen C statements for variable declarations and setup, one call to actually solve the problem (determine times when Venus crosses Jupiter's equator), and a few dozen lines for output. Plus the kernel management which is not tough for this problem.

SPICE has a steep learning curve; let me know if you want to pursue this.

[brandon-rhodes]

I agree with @barrycarter that probably only a rotation or two is necessary, given the pole direction on the Wikipedia page, to get an approximate position for a Jovian coordinate system. But without a way to generate a correct coordinate to aim for, I have no idea whether I have rotated in the right direction or the wrong one for this first try using Skyfield:

# From Jupiter, observe the position of Venus

from skyfield.api import JulianDate, venus, jupiter
from skyfield.functions import rot_x, rot_z
from skyfield.positionlib import to_polar
from skyfield.units import Angle, Distance

jd = JulianDate(utc=(2014, 8, 5))
xyz = jupiter(jd).observe(venus).position.AU

# The resulting ICRS x,y,z coordinates need to be rotated to directions
# relative to Jupiter's equator and pole.  (These rotations are probably
# backwards.  Or maybe even just wrong.  Can someone generate some
# coordinates using another tool to double-check?)

jpole_ra = Angle(degrees=268.06)
jpole_dec = Angle(degrees=64.50)

W = (rot_z(-jpole_ra.radians)).dot(rot_x(-jpole_dec.radians))

xyz2 = (W).dot(xyz)
d, theta, phi = to_polar(xyz2)
print 'Jupiter declination?', Angle(radians=theta)
print 'Jupiter HA?', Angle(radians=phi, preference='hours')
print Distance(AU=d)

This all strikes me as a bit clunky, but I think it is because Skyfield has no idea yet about the coordinate systems of other planets — so we are having to do the rotations, and change to polar coordinates, out in the open whereas normally it would occur under the covers if you were asking about a traditional Earth-based coordinate system.

[ghost]

Brandon,

The resulting ICRS x,y,z coordinates need to be rotated to directions

relative to Jupiter's equator and pole. (These rotations are probably

backwards. Or maybe even just wrong. Can someone generate some

coordinates using another tool to double-check?)

Actually, I believe Horizons will do this (as the OP discovered). You can use the web interface, or send this email to horizons@ssd.jpl.nasa.gov with the subject "JOB":

!$$SOF COMMAND= '499' CENTER= '@799' OBJ_DATA= 'YES' MAKE_EPHEM= 'YES' TABLE_TYPE= 'OBS' START_TIME= '2014-01-01' STOP_TIME= '2015-01-01' STEP_SIZE= '1 d' QUANTITIES= '1' CAL_FORMAT= 'CAL' REF_SYSTEM= 'J2000' ANG_FORMAT= 'DEG' RANGE_UNITS= 'AU' APPARENT= 'AIRLESS' SOLAR_ELONG= '0,180' SUPPRESS_RANGE_RATE= 'YES' SKIP_DAYLT= 'NO' EXTRA_PREC= 'NO' R_T_S_ONLY= 'NO' CSV_FORMAT= 'YES' !$$EOF

and you will get back the attached (so my example above is for reference purposes only). This gives the position of Mars as seen from Uranus.

 Automated mail xmit by MAIL_REQUEST, PID=   3552 Sun Aug  3 15:58:40 2014
++++++++++++++++++++++++++++++++ (part 1 of 1)  +++++++++++++++++++++++++++++++

*******************************************************************************
 Revised: Sep 28, 2012                 Mars                             499 / 4

 GEOPHYSICAL DATA (updated 2009-May-26):
  Mean radius (km)      = 3389.9(2+-4)    Density (g cm^-3)     =  3.933(5+-4)
  Mass (10^23 kg )      =    6.4185       Flattening, f         =  1/154.409
  Volume (x10^10 km^3)  =   16.318        Semi-major axis       =  3397+-4
  Sidereal rot. period  =   24.622962 hr  Rot. Rate (x10^5 s)   =  7.088218
  Mean solar day        =    1.0274907 d  Polar gravity ms^-2   =  3.758
  Mom. of Inertia       =    0.366        Equ. gravity  ms^-2   =  3.71
  Core radius (km)      =  ~1700          Potential Love # k2   =  0.153 +-.017

  Grav spectral fact u  =   14 (x10^5)    Topo. spectral fact t = 96 (x10^5)
  Fig. offset (Rcf-Rcm) = 2.50+-0.07 km   Offset (lat./long.)   = 62d / 88d
  GM (km^3 s^-2)        = 42828.3         Equatorial Radius, Re = 3394.0 km 
  GM 1-sigma (km^3 s^-2)= +- 0.1          Mass ratio (Sun/Mars) = 3098708+-9

  Atmos. pressure (bar) =    0.0056       Max. angular diam.    =  17.9"
  Mean Temperature (K)  =  210            Visual mag. V(1,0)    =  -1.52
  Geometric albedo      =    0.150        Obliquity to orbit    =  25.19 deg
  Mean sidereal orb per =    1.88081578 y Orbit vel.  km/s      =  24.1309
  Mean sidereal orb per =  686.98 d       Escape vel. km/s      =   5.027
  Hill's sphere rad. Rp =  319.8          Mag. mom (gauss Rp^3) = < 1x10^-4 
*******************************************************************************

*******************************************************************************
Ephemeris / MAIL_REQUEST Sun Aug  3 15:58:40 2014 Pasadena, USA  / Horizons    
*******************************************************************************
Target body name: Mars (499)                      {source: MAR097}
Center body name: Uranus (799)                    {source: URA111}
Center-site name: BODYCENTRIC
*******************************************************************************
Start time      : A.D. 2014-Jan-01 00:00:00.0000 UT      
Stop  time      : A.D. 2015-Jan-01 00:00:00.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_MARS                        {East-longitude -}
Target radii    : 3396.2 x 3396.2 x 3376.2 km     {Equator, meridian, pole}    
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : IAU_URANUS                      {East-longitude +}
Center radii    : 25559.0 x 25559.0 x 24973.0 km  {Equator, meridian, pole}    
Target primary  : Sun
Vis. interferer : TITANIA (R_eq= 788.900) km      {source: URA111}
Rel. light bend : Sun, URANUS                     {source: URA111}
Rel. lght bnd GM: 1.3271E+11, 5.7940E+06 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : CAL 
EOP file        : eop.140801.p141023                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2014-AUG-01. PREDICTS-> 2014-OCT-22
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass     n.a.    , Daylight (NO )
Table cut-offs 2: Solar Elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table format    : Comma Separated Values (spreadsheet)
*******************************************************************************
 Date__(UT)__HR:MN, , ,R.A._(ICRF/J2000.0), DEC_(ICRF/J2000.0),
***************************************************************
$$SOE
 2014-Jan-01 00:00, , ,188.26780, -2.71827,
 2014-Jan-02 00:00, , ,188.30290, -2.73395,
 2014-Jan-03 00:00, , ,188.33813, -2.74968,
⋮
$$EOE
*******************************************************************************
Column meaning:

TIME

  Prior to 1962, times are UT1. Dates thereafter are UTC. Any 'b' symbol in
the 1st-column denotes a B.C. date. First-column blank (" ") denotes an A.D.
date. Calendar dates prior to 1582-Oct-15 are in the Julian calendar system.
Later calendar dates are in the Gregorian system.

  Time tags refer to the same instant throughout the universe, regardless of
where the observer is located.

  The dynamical Coordinate Time scale is used internally. It is equivalent to
the current IAU definition of "TDB". Conversion between CT and the selected
non-uniform UT output scale has not been determined for UTC times after the
next July or January 1st.  The last known leap-second is used over any future
interval.

  NOTE: "n.a." in output means quantity "not available" at the print-time.

 R.A._(J2000.0)_DEC. =
   J2000.0 astrometric right ascension and declination of target center.
Adjusted for light-time. Units: DEGREES

 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information: http://ssd.jpl.nasa.gov/
     Connect    : telnet://ssd.jpl.nasa.gov:6775  (via browser)
                  telnet ssd.jpl.nasa.gov 6775    (via command-line)
     Author     : Jon.Giorgini@jpl.nasa.gov

*******************************************************************************

[drbitboy]

NAIF/SPICE solution: 82 Venus crossings of the Jupiter equator between 1901 and 2100; one-second resolution; currently uses DE-430 and JUP-310 ephemerides and IAU Jupiter frame; trivially translatable to almost any other observer/target bodies. Took me an hour and a half to write; about thirty lines of code.

https://github.com/drbitboy/vj

[ifmihai]

@drbitboy SPICE is very nice indeed. Powerful even though a little steep learning curve. Thank you for pointing it out and for python example code. Lots to crunch, but its a great tool, the only one apparently for now.

@barrycarter if stellarium and spice can move observer to other planets, can we get the method out of their code? It seems the most straight forward way.

[brandon-rhodes]

@barrycarter Observing one planet from another is a single line of code in Skyfield. Here, for example, is the Skyfield code to generate the first position from your HORIZONS result:

from skyfield.api import uranus, mars
ra, dec, distance = uranus(utc=(2014, 1, 1)).observe(mars).radec()
print ra._degrees, dec.degrees

The problem is that these two coordinates are not a sufficient answer for this question, because they are ICRS coordinates which measure angles relative to the Earth’s equator and poles. The question is about whether we can express these same coordinates relative to the equator and poles of the planet from which we are observing. That is what, you will see, I was trying to do with the pair of rotations in my sample code above.

I did click through all of the HORIZONS features last night, and did not find any ability to use local coordinate systems like “the dynamical reference system of Jupiter” or anything — so my guess is that HORIZONS cannot do this calculation.

@drbitboy Thank you, I will take a look at these results!

[drbitboy]

Did you try the frame IAU_JUPITER with Horizons? I haven't been there in a while, I am just guessing.

Using SPICE code is probably an all or nothing proposition; either you put it under the hood of your app, or you encapsulate it in a web app (i.e. I think it's integrated into Horizons for some functions), or you don't use it.

The real beauty of SPICE is that, since the 1980s, it has separated the code (functions) from data:

S - spacecraft ephemeris; P - Planet (and any body or surface point) ephemerides and constants I - Instrument geometry; C - C-matrix i.e. pointing and attitude E - Events; implemented as a relational database (think SQLite)

It's an objected oriented toolkit originally written in FORTRAN. And the documentation is outstanding; here is one of my favorites:

C C Mathematical background for the solution---the lambda equation. C ================================================================ C ---------------------------------------------------------------- C C
C Here is the background and general outline of how this problem is C going to be solved. C
C We want to find, a point on the ellipsoid C
C
C X2 Y2 Z2 C ------ + ------ + ------ = 1 C A2 B2 C2 C
C that is closest to the input point POSITN. C
C If one cares about the gory details, we know that C such a point must exist because the C ellipsoid is a compact subset of Euclidean 3-space C and the distance function between the input point C and the ellipsoid is a continuous functions. C Since continuous functions on compact sets C actually achieve their minimums at some point of C the compact set, we are guaranteed that a closest C point exists.

[ifmihai]

@brandon-rhodes this is the web GUI interface to SPICE

I did a search for declination == 0, target, observer, reference = EARTH, JUPITER, IAU_JUPITER

it would be interesting to see if SPICE is different than horizon (I doubt it) If SPICE with IAU_JUPITER doesn't work also, then...

[drbitboy]

@brandon-rhodes Did you see a "Reference plane [eclip, frame, body ]" option in your clicking through the Horizons features last night?

I was able to coax the Horizons telnet interface (port 6775) to put Venus in the Jupiter frame (I think; the results eyeball-agree to within a minute of the SPICE result for the 2014 equator crossing):

https://github.com/drbitboy/vj/blob/master/horizons.md

I used python and the expect module to extract data from the Horizons telnet interface in PySPICE (see the test code), and I once wrote a VBA script to fill out an eXcel spreadsheet via the same route.

[brandon-rhodes]

@ifmihai — What were the results of your search — on what dates and times did it report a declination of zero?

@drbitboy — No. I just clicked through it again, and see only ICRF and FK5 as possible output reference frames.

[brandon-rhodes]

Here is the answer.

The Python code in the following notebook solves the original question for Jupiter and Venus. Given the other planets’ pole orientations as given in the Wikipedia article cited above, it could do so for any planet pairing.

It took me a few tries to get the rotations in the right direction; I think they are correct now. You will see that the zero crossings now agree to very good approximation to the SPICE predictions @drbitboy offered above:

http://nbviewer.ipython.org/gist/brandon-rhodes/60503689abeef7fc099c

Thank you VERY much, @drbitboy — without your having solved the problem in SPICE, I would have had no way to find and verify the correct rotation directions to convert into Jovian coordinates!

I suspect that the SPICE solution is superior since it is presumably approximating the Jovian pole progression over the eons, and is thus more accurate than the rough point-in-time figures given in the Wikipedia.

In the future the operation might be even simpler in Skyfield, without requiring these extra imports and manual rotations, but for now this is the best approach.

[brandon-rhodes]

And, no, I had not heard of William Foster before — PyEphem users doing astrology before have, I gather, been using systems that only use the relative positions of other planets to Earth rather than positing that planet-to-planet relationships not involving the Earth could influence events here.

[drbitboy]

@brandon-rhodes

On Horizons, you need to choose vectors for [Ephemeris type], then in the [Table settings](the scalpels go on the left, with the pitchforks. Igor! Igor!) you can select the option for body-fixed coordinates.

[brandon-rhodes]

@drbitboy — And there I had gone, clicking through all of the options under "Ephemeris Type: OBSERVER" thinking that I had exhausted the possibilities! Thank you for the detailed diagram. The HORIZONS web interface can indeed produce this result.

[drbitboy]

This is like the experienced coder sitting behind the newbie who just read the emacs or vi manual, and the newbie does something and the coder says "Hey! How did you do that?"

There are just too many options on this blasted interwebnet thingy.

[ifmihai]

@brandon-rhodes Here is the list for venus, jupiter, with reference IAU_JUPITER, from 1950 to 2020.

This is copy paste, columns being: nr, Start TimeStop TimeDuration (secs)

1 1950-01-25 10:11:01.502326 UTC 1950-01-25 10:11:01.502326 UTC 0.00000000E+00 2 1950-06-09 16:44:03.385450 UTC 1950-06-09 16:44:03.385450 UTC 0.00000000E+00 3 1950-07-07 12:26:40.785220 UTC 1950-07-07 12:26:40.785220 UTC 0.00000000E+00 4 1955-05-29 23:56:05.298093 UTC 1955-05-29 23:56:05.298093 UTC 0.00000000E+00 5 1955-07-01 03:20:29.997507 UTC 1955-07-01 03:20:29.997507 UTC 0.00000000E+00 6 1955-11-12 08:23:25.496784 UTC 1955-11-12 08:23:25.496784 UTC 0.00000000E+00 7 1961-10-18 12:19:02.930748 UTC 1961-10-18 12:19:02.930748 UTC 0.00000000E+00 8 1962-01-17 04:50:43.869918 UTC 1962-01-17 04:50:43.869918 UTC 0.00000000E+00 9 1962-04-09 08:24:24.240686 UTC 1962-04-09 08:24:24.240686 UTC 0.00000000E+00 10 1967-08-04 02:41:24.750592 UTC 1967-08-04 02:41:24.750592 UTC 0.00000000E+00 11 1973-07-17 10:35:59.737340 UTC 1973-07-17 10:35:59.737340 UTC 0.00000000E+00 12 1973-08-29 21:00:02.450309 UTC 1973-08-29 21:00:02.450309 UTC 0.00000000E+00 13 1974-01-01 04:33:28.179947 UTC 1974-01-01 04:33:28.179947 UTC 0.00000000E+00 14 1979-04-26 05:43:53.084149 UTC 1979-04-26 05:43:53.084149 UTC 0.00000000E+00 15 1979-08-06 12:19:57.630937 UTC 1979-08-06 12:19:57.630937 UTC 0.00000000E+00 16 1979-10-17 10:02:58.720847 UTC 1979-10-17 10:02:58.720847 UTC 0.00000000E+00 17 1985-09-23 16:40:28.942229 UTC 1985-09-23 16:40:28.942229 UTC 0.00000000E+00 18 1991-01-19 16:46:16.732645 UTC 1991-01-19 16:46:16.732645 UTC 0.00000000E+00 19 1991-03-22 16:07:08.698427 UTC 1991-03-22 16:07:08.698427 UTC 0.00000000E+00 20 1991-07-11 14:57:21.394607 UTC 1991-07-11 14:57:21.394607 UTC 0.00000000E+00 21 1997-06-16 03:41:34.577125 UTC 1997-06-16 03:41:34.577125 UTC 0.00000000E+00 22 1997-10-05 23:05:35.398998 UTC 1997-10-05 23:05:35.398998 UTC 0.00000000E+00 23 1997-12-03 19:47:42.964686 UTC 1997-12-03 19:47:42.964686 UTC 0.00000000E+00 24 2003-04-02 14:41:04.377851 UTC 2003-04-02 14:41:04.377851 UTC 0.00000000E+00 25 2009-03-11 02:10:09.145514 UTC 2009-03-11 02:10:09.145514 UTC 0.00000000E+00 26 2009-05-20 01:45:02.436037 UTC 2009-05-20 01:45:02.436037 UTC 0.00000000E+00 27 2009-08-30 04:49:55.770999 UTC 2009-08-30 04:49:55.770999 UTC 0.00000000E+00 28 2014-12-23 05:05:26.629921 UTC 2014-12-23 05:05:26.629921 UTC 0.00000000E+00 29 2015-04-25 21:32:45.501325 UTC 2015-04-25 21:32:45.501325 UTC 0.00000000E+00 30 2015-06-11 12:29:53.269830 UTC 2015-06-11 12:29:53.269830 UTC 0.00000000E+00

@brandon-rhodes I doubt Foster was into astrology. When a men refuses to give away his method, I believe he is beyond (usual crap) astrology. Personally, I don't care, I want to test as I was intrigued by the person and his ways.

@drbitboy emacs/vi image was funny :P

[brandon-rhodes]

@ifmihai I am glad you now have the list of events you were looking for, hopefully together with a good way to produce them in the future for the planets you are interested in! And I am happy to see that your list includes the third zero-crossing in late 2009 that SPICE missed, but that is clearly visible in the graph from my notebook:

http://nbviewer.ipython.org/gist/brandon-rhodes/60503689abeef7fc099c

Since PyEphem itself will never be capable of a calculation like this because of the limitations of the underlying C library, I am going to go ahead and close this issue. Good luck!

[ghost]

Here's a solution to the OP's question that I'm particularly proud of because it is not instructive or useful for anything else in any way.

Example: find when Mercury crosses Jupiter's equator.

Wikipedia tells us Jupiter's north pole points to 268.057 degrees right ascension and 64.496 degrees declination in the ICRS plane.

Converting this to rectangular coordinates (with an arbitrary radius of 1) gives us the vector:

[-0.0145987218963993, -0.430326550289791, 0.902555226805917]

Jupiter's equator is the plane perpendicular to this vector going through the center of Jupiter.

Thus, when Jupiter is at [jx, jy, jz], the plane of Jupiter's equator is defined by:

-0.0145987218963993(x-jx) + -0.430326550289791(y-jy) + 0.902555226805917*(z-jz) = 0

Thus, if Mercury is at [mx, my, mz], we are checking if:

-0.0145987218963993(mx-jx) + -0.430326550289791(my-jy) + 0.902555226805917*(mz-jz) = 0

I tested using this program:

for t in range (2456876,2456876+1000): mx,my,mz = eph.position('mercury', t) jx,jy,jz = eph.position('jupiter', t) eq = -0.0145987218963993(jx-mx)-0.430326550289791(jy-my)+0.90255522680591\ 7*(jz-mz) print "ALPHA",t,eq

and looking at the results visually. Two crossings I found:

ALPHA 2457004 [-392992.70703286] ALPHA 2457005 [ 18670.08156499]

ALPHA 2457040 [ 472299.13722086] ALPHA 2457041 [-38209.33185893]

which correspond to "2014-12-12 12:00:00" and "2015-01-17 12:00:00". Stellarium confirms these are accurate.

Of course, you should really be using a binary search, further test this code, and improve it many other ways before using it "in production".

[ghost]

@brandon-rhodes It turns out a rotation is not sufficient because there are an infinite number of rotations that will map a given RA/DEC to the north pole. However, any of these rotations will yield the same declination. So, you can use a rotation to find the declination, which technically answers the OP's question :)

[ghost]

Will skyfield include IAU_JUPITER and other planetary reference frames?

[ifmihai]

@barrycarter sounds good (but i dont understand the details). But +1 for IAU_planet reference systems. As these reference systems i used with webgeocalc and they gave(apparently) accurate results.

[brandon-rhodes]

I have opened an issue on Skyfield and will add that feature — thanks for the idea!

[drbitboy]

Brandon, going back over old items; I may get an opportunity to use SPICE for astrological calculation.

What crossing in 2009 did SPICE miss? I set the step size to 10d, so it certainly is possible, but I see three in my results here.

[brandon-rhodes]

Comparing your list with the short list in my code:

http://nbviewer.jupyter.org/gist/brandon-rhodes/60503689abeef7fc099c

it's the 2009-08 crossing that's not marked on my graph as a SPICE-detected crossing. I can't now decide, all these years later, whether I was cutting and pasting from a SPICE result set of my own (or someone else's besides yours), or whether I made a cut and paste error that left 2009-08 out. So unless we find new evidence, let's assume that the error was mine, and that SPICE in fact found that crossing!