Non-Terran satellite support?

[ghost]

Does skyfield provide position data for non-Terran satellites, in particular the 4 large Jovian satellites, Titan, and other non-Terran satellites "easily" visible from Earth?

[brandon-rhodes]

No, I have not yet run across any good high-precision data sets predicting their positions. I want something like “DE421 but for satellites” — do you happen to know of any?

[ghost]

I assume you are familiar with:

http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/

or is that too inaccurate?

[ghost]

Brandon, I was wondering if you'd figured out how to read the files in http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/

I think I just figured it out, but, if you've already done it, no point in posting my findings.

[brandon-rhodes]

No, I was too busy in August to tackle this! Please do feel free to paste in any code snippets or adventures that you had in deciphering and using this data. Thanks!

[ghost]

Brandon, I tried to writeup a document explaining this, but I'm not happy with it: it's too long and verbose and implies finding satellite positions is difficult, when it's really not.

Despite this, it should provide enough information for programmers. The "living document" is at:

https://github.com/barrycarter/bcapps/blob/master/ASTRO/README.bsp

and I've cut and pasted the current version below:

This informal document explains how to decode the satellite SPICE/SPK .bsp files found at:

http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/

which are essentially sets of Chebyshev coefficients.

Intended audience: you should know what the following terms mean:

• DE430
• SPICE
• Chebyshev polynomials

Summary of steps:

• Convert .bsp to .xsp
• Convert DAF numbers to decimals
• Profit!
• Understand file format

Step 1: Using sat365 as an example, run "toxfr" on sat365.bsp. You can find toxfr at:

ftp://naif.jpl.nasa.gov/pub/naif/toolkit/

You'll have to dig a bit to find your platform/language. I got mine at:

ftp://naif.jpl.nasa.gov/pub/naif/toolkit/C/PC_Linux_GCC_32bit/executables/

Note: Also download all other files in the toolkit, all are useful.

Step 2: "toxfrp sat365.bsp" yields "sat365.xsp". Looking in sat365.xsp, we see quoted numbers like:

'-E6BEBE0789D848^-2'

These numbers are in IEEE-754 hex format and used by Fortran. In base 16, the number above is:

-0.E6BEBE0789D848*(10^-2)

where "10" means 16, since we're in base 16.

To convert do decimal, we first convert "E6BEBE0789D848" to decimal, yielding 64948968022988872 [1].

"E6BEBE0789D848" is 14 characters long, so we divide by 16^14, yielding (to 30 digits):

0.901347996559565989294071641780

"^-2" means we multiply by 16^-2, yielding (30 digits):

0.00352089061156080464567996735070

Of course, we could've combined the two previous steps and just done:

hex("E6BEBE0789D848")/16^16

to get the same thing. More efficient methods almost certainly exist.

bc-ieee754.pl in this directory (https://github.com/barrycarter/bcapps/tree/master/ASTRO) will also do the conversion.

Step 3: Looking at the last few lines of sat365.xsp, we see things like:

Master Interval: 18.0000 Days

Name Number GM NDIV NDEG Model Mimas 601 2.503524000000000E+00 24 15 SATORBINT Enceladus 602 7.211454165826105E+00 24 12 SATORBINT

Taking Enceladus as an example:

• Since Enceladus is the 2nd satellite listed, it's data will follow "BEGIN_ARRAY 2 7792004" in the xsp file (the 7792004 is the length of the array, so the next 7792004 lines are elements of the array) [2]
• NDIV=24 means each set of coefficients covers 18/24 days (18 = master interval).
• NDEG=12 means each coordinate is a 12th degree polynomial, so there will be 13 (not 12) Chebyshev coefficients for each coordinate.

IMPORTANT: The NDEG convention used here is different from the convention used in GROUP 1050 of header.430_572

Step 4: Looking at the first few lines after "BEGIN_ARRAY 2 7792004" after decimal conversion ("sat365.xsp.dec" if you're using bc-ieee754.pl), we see:

BEGIN_ARRAY 2 7792004 'SAT365 ' '-3155716758.81603' '3155716867.18389' '25A' '6' '1' '3' 1024 '-3155727600' '32400' '31713.2412523424'

"-3155727600" is the number of seconds after the Julian day 2451545 (2000-01-01T12:00:00), the same as the "Equinox Reference Julian Date".

IMPORTANT: Ignore the "File Epoch Julian Date": it is not useful to us

The 32400 means the first set of coefficients are good for -3155727600 seconds after the epoch (ie, 1899-12-31T21:00:00) plus or minus 32400 seconds.

TODO: Why isn't the epoch date more prominently mentioned?

The "31713.2412523424" above is the first Chebyshev coefficient itself.

Step 5: Since NDEG=12 for Enceladus, each coordinate has 13 (not 12) Chebyshev coefficients. The coordinates are for X, Y, Z, VX, VY, and VZ [3].

Thus, each set has 13 coefficients times 6 coordinates for a total of 78 lines. If we count "BEGIN_ARRAY 2 7792004" as line 1:

Lines 12-24 are Chebyshev coefficients for X Lines 25-37 are Chebyshev coefficients for Y Lines 38-50 are Chebyshev coefficients for Z Lines 51-63 are Chebyshev coefficients for VX Lines 64-76 are Chebyshev coefficients for VY Lines 77-89 are Chebyshev coefficients for VZ

Note the 90th line after "BEGIN_ARRAY 2 7792004" is "-3155662800", indicating the start of the next set of coefficients.

Important notes:

• Ignore any unquoted numbers (most often "1024") in xsp files: these are NOT array entries [4]
• To see .bsp info without converting to .xsp, use "commnt -r sat365.bsp" and "brief -c sat365.bsp"
• Satellite positions are given relative to the barycenter of the satellite's planetary system, which means:
  • Positions are NOT given relative to the solar system barycenter.
  • Positions are NOT given relative to the planet itself. The planet's barycenter is the center of mass of the planet and its satellites.
  • The position of the planet itself (relative to it's own barycenter) is provided, in this case as the 13th array.

Spotchecking the answer (this section is incomplete!!!):

• To check via email, send the following email (for example) to horizons@ssd.jpl.nasa.gov:

!$$SOF COMMAND= '602' CENTER= '500@6' MAKE_EPHEM= 'YES' TABLE_TYPE= 'VECTORS' START_TIME= '2014-09-02' STOP_TIME= '2014-09-12' STEP_SIZE= '1 h' OUT_UNITS= 'KM-S' VECT_TABLE= '3' REF_PLANE= 'FRAME' REF_SYSTEM= 'J2000' VECT_CORR= 'NONE' VEC_LABELS= 'NO' CSV_FORMAT= 'YES' OBJ_DATA= 'YES' !$$EOF

• To check via the Horizons web interface, use these settings (for example):

Ephemeris Type [change] : VECTORS Target Body [change] : Enceladus (SII) [602] Coordinate Origin [change] : Saturn System Barycenter [500@6] Time Span [change] : Start=2014-09-02, Stop=2014-09-12, Step=1 h Table Settings [change] : output units=KM-S; reference plane=FRAME; CSV format=YES Display/Output [change] : plain text

IMPORTANT NOTES (if using web interface):

* The coordinate origin must be "Saturn System Barycenter

[500@6]", NOT Saturn itself.

* Under Table Settings, the reference plane must be "FRAME".

TODO: add more steps above

TODO: spellcheck

TODO: figure out what the GM heading means (although it appears to be unimportant)

[1] Some languages/operating systems may choke on this. My Perl non-fatally warns "Integer overflow in hexadecimal number", but yields "6.49489680229889e+16", an approximation to the actual value. You can probably use clever bit-shifting to avoid this.

[2] This isn't strictly true, but I explain more later.

[3] The coefficients for VX, VY, VZ are redundant and can be obtained by differentiating the polynomials for X, Y, and Z.

[4] These numbers probably had legacy use at some point, indicating, for example, that the next 1024 lines were array entries (useful to know when machines had less memory).

[brandon-rhodes]

@barrycarter — thank you very much for putting together this information! These ephemerides are important resources that I think Skyfield should be able to process natively, maybe as soon as the next version.

Having looked over some of the documentation of the .bsp format, I think, happily, that we can skip the step of turning the ephemeris into text — which, after all, would then have the unhappy result that Python would have to parse the decimal numbers to rebuild the binary floating-point numbers again. Instead, I think that we can memory-map a .bsp file that the user is interested in, and directly view the data as a NumPy array using the ndarray() type’s ability to understand a raw array of double-precision IEEE floats, which is exactly what a .bsp file happens to be made of.

A quick experiment in reading the binary file format has started producing coefficients given only about two hours of work on my part, so now the final challenge will be hooking the numbers up to the polynomial logic that I already have inside of jplephem. I expect results within a few days, and will post back here with the results.

Thanks again for prodding me to move Skyfield in this direction!

[brandon-rhodes]

Oh — and, for reference, the way to read a .bsp file is described in two documents, that should be read in this order:

http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/daf.html http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/spk.html

[brandon-rhodes]

After 20 hours of work on a new version, I have released jplephem 2.0 with full support for SPK kernels:

https://pypi.python.org/pypi/jplephem

Now to figure out how to get this integrated into Skyfield. Stay tuned!

[ghost]

Brandon, I 100% agree that the binary version is better, but I couldn't figure it out, even after reading the references above. I'm sure it gives the same information as the text format, but I couldn't figure out how to map the bytes in the bsp file to the information I wanted.

The references mostly mention looking at the FORTRAN code, which I didn't really understand.

Is there any chance you could create a brief human-readable text version of how BSP files store their information? I'd like to create a Perl subroutine that reads them in native binary format, instead of the subroutine I have now which reads them in text format.

[brandon-rhodes]

If you can read Python code, here is how to read the underlying binary format:

https://github.com/brandon-rhodes/python-jplephem/blob/master/jplephem/daf.py

and then how to read the .spk records that is built on top of it:

https://github.com/brandon-rhodes/python-jplephem/blob/master/jplephem/spk.py

[brandon-rhodes]

And, Skyfield now ships with native support for SPK files! You can learn the positions of the Jovian satellites by asking Skyfield to download the (enormous) jup310.bsp file:

eph = api.load('jup310.bsp')
astrometric = eph['earth'].observe(eph['callisto']).at(tdb=2471184.5)

The positions seem to agree with HORIZONS to within 1mm, so I am going to close this issue, but please let me know if you have any further ideas for Skyfield. Thanks!