Support ephemerides (like DE441) that have 2 segments per planet

[blucobalto]

Hi, sorry if I come back to the problem but I have a doubt. de441.bsp: range date is wrong https://github.com/skyfielders/python-skyfield/issues/682

If i use the files

• de441_part-1.bsp (1.6GB+)
• de441_part-2.bsp (1.6GB+)

downloaded from https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/ Skyfield has no problems.

If i use the file

• de441.bsp (3.2GB+)

downloaded from ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/ Skyfield returns this error:

Exception has occurred: EphemerisRangeError
ephemeris segment only covers dates 1969-07-29 23:59:18Z through 17191-03-14 23:58:51Z UT

But shouldn't it be the same?

Thanks a lot, Blu

[brandon-rhodes]

It's possible that the combined file still contains two separate segments for each planet, and maybe Skyfield only recognizes one of them instead of knowing how to switch back and forth depending on date. Try this:

python -m jplephem spk de441.bsp

[blucobalto]

I wouldn't want to be wrong, I understand that where the first segment ends, the second begins

-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
2440432.50..8000016.50  Type 2  Venus Barycenter (2) -> Venus (299)

This is the result:

File type DAF/SPK and format LTL-IEEE with 28 segments:
-3100015.50..2440432.50  Type 2  Venus Barycenter (2) -> Venus (299)
-3100015.50..2440432.50  Type 2  Mercury Barycenter (1) -> Mercury (199)
-3100015.50..2440432.50  Type 2  Earth Barycenter (3) -> Earth (399)
-3100015.50..2440432.50  Type 2  Earth Barycenter (3) -> Moon (301)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Sun (10)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
-3100015.50..2440432.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
2440432.50..8000016.50  Type 2  Venus Barycenter (2) -> Venus (299)
2440432.50..8000016.50  Type 2  Mercury Barycenter (1) -> Mercury (199)
2440432.50..8000016.50  Type 2  Earth Barycenter (3) -> Earth (399)
2440432.50..8000016.50  Type 2  Earth Barycenter (3) -> Moon (301)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Sun (10)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
2440432.50..8000016.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)

[brandon-rhodes]

Just as we suspected: the file contains two separate segments for each planet, and Skyfield expects only one segment per planet. For example, for Venus:

-3100015.50..2440432.50  Type 2  Venus Barycenter (2) -> Venus (299)
2440432.50..8000016.50  Type 2  Venus Barycenter (2) -> Venus (299)

I think there is an open issue somewhere involving this. I will make position generation slower, because instead of a simple vector operation, Skyfield will have to build a True/False mask over each array of dates to determine which dates need to be computed using the first array and which need to use the second array.

I'm going to update the task description. I'm the only very active contributor at the moment, and had a few other Skyfield improvements planned for the near-term — in the meantime is there another ephemeris that could support your project? How did you wind up selecting DE441 in particular? (Not meaning to criticize your choice of ephemeris! But more information will help me understand the level of urgency here.)

[blucobalto]

Hi, first of all thank you for your work on Skyfield. I chose the de441 for the rather large date range. I have no real urgency. I can continue working with the two separate files.

If there are urgencies, think about giving priority to those. There is no problem.

Thank you very much for your availability, Blu

P.S. sorry for my English. I hope I was clear

[brandon-rhodes]

P.S. sorry for my English. I hope I was clear

Yes, very clear, everything was easy to understand. Thanks!

[tzengyuxio]

Same problem here. I use workaround like this(not actual code, just for concept):

eph411a = load("de441_part-1.bsp")
eph411b = load("de441_part-2.bsp")

JD_SEPERATE = 2440432.50 

def find_discrete(start_time, end_time):
    if start_time < end_time < JD_SEPERATE:
        return almanac.find_discrete(start_time, end_time, almanac.seasons(eph411a))
    if JD_SEPERATE < start_time < end_time:
        return almanac.find_discrete(start_time, end_time, almanac.seasons(eph411b))
    if start_time < JD_SEPERATE and end_time > JD_SEPERATE:
        ta, ya = almanac.find_discrete(start_time, JD_SEPERATE, almanac.seasons(eph411a))
        tb, yb = almanac.find_discrete(JD_SEPERATE, end_time, almanac.seasons(eph411b)) 
        return ta + tb, ya + yb

[psbaltar]

I ran into the same issue with plu058.bsp. Unfortunately, it seems to be the only ephemeris for the Pluto and its moons. While it would be nice to be able to use it, my project is fine with PLUTO BARYCENTER from another ephemeris. So no added urgency from me, just an FYI.

In a previous comment, it was mentioned that adding a date check to choose the segment would slow things down. What an option to choose the segment? i.e. instead of automatically choosing which segment to use based on the date, have the user tell it which to use. That way, there's no performance impact if it isn't actually needed. In my specific case, I'd be fine with the 2015-2099 segment, but Skyfield is automatically picking the 1900-2015 segment. So plu058 is only usable from 1900-2015 even though there's more data in it.

https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/plu058.bsp

SPICE kernel file 'plu058.bsp' has 16 segments
  JD 2457357.50 - JD 2488065.50  (2015-11-30 through 2099-12-27)
      9 -> 901  PLUTO BARYCENTER -> CHARON
      9 -> 902  PLUTO BARYCENTER -> NIX
      9 -> 903  PLUTO BARYCENTER -> HYDRA
      9 -> 904  PLUTO BARYCENTER -> KERBEROS
      9 -> 905  PLUTO BARYCENTER -> STYX
      9 -> 999  PLUTO BARYCENTER -> PLUTO
  JD 2415021.50 - JD 2457357.50  (1900-01-01 through 2015-11-30)
      9 -> 901  PLUTO BARYCENTER -> CHARON
      9 -> 902  PLUTO BARYCENTER -> NIX
      9 -> 903  PLUTO BARYCENTER -> HYDRA
      9 -> 904  PLUTO BARYCENTER -> KERBEROS
      9 -> 905  PLUTO BARYCENTER -> STYX
      9 -> 999  PLUTO BARYCENTER -> PLUTO
  JD 2415021.50 - JD 2488065.50  (1900-01-01 through 2099-12-27)
      3 -> 399  EARTH BARYCENTER -> EARTH
      0 -> 10   SOLAR SYSTEM BARYCENTER -> SUN
      0 -> 9    SOLAR SYSTEM BARYCENTER -> PLUTO BARYCENTER
      0 -> 3    SOLAR SYSTEM BARYCENTER -> EARTH BARYCENTER

[brandon-rhodes]

@psbaltar — One approach would be to throw out the segments that end in 2015 after loading the ephemeris. If you are loading that exact ephemeris, just toss out segments 6 through 11:

from skyfield.api import load

eph = load('./plu058.bsp')
eph.segments = eph.segments[0:6] + eph.segments[12:]

Or be fancier and filter them by date:

ts = load.timescale()
eph = load('./plu058.bsp')
eph.segments = [segment for segment in eph.segments
                if segment.time_range(ts)[1].utc.year != 2015]

After either of the above steps, I think you'll find that charon = eph['Charon'] gets the more recent of the two segments.

Unfortunately you will not find that either of the above maneuvers changes what you see when you print(eph) because that routine steps through the contents of the file itself, rather than the .segments list derived from it. Maybe the next version of Skyfield should adjust that, so it steps through the .segments list instead?