Issues when building frames using DE440

[AntoineRichard]

Hi there,

I'm running into issues using the ephemeris DE440 and observing planets from the moon.
I think there is a fair chance that I'm just mentally challenged and missed something, but on the account I didn't here is how to reproduce the issue:

The files you will need:

wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/moon_de440_220930.tf
wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/fk/satellites/a_old_versions/moon_de440_200625.tf
wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/pck00011.tpc
wget https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de440_200625.bpc

The code (it's the same as the demo code for that):

from skyfield.api import load, PlanetaryConstants

ts = load.timescale()
t = ts.utc(2019, 12, 20, 11, 5)

eph = load('de440.bsp')
earth, moon = eph['earth'], eph['moon']

pc = PlanetaryConstants()
pc.read_text(load('moon_de440_220930.tf'))
pc.read_text(load('pck00011.tpc'))
pc.read_binary(load('moon_pa_de440_200625.bpc'))

frame = pc.build_frame_named('MOON_ME_DE440_ME421')
aristarchus = moon + pc.build_latlon_degrees(frame, 26.3, -46.8)

apparent = aristarchus.at(t).observe(earth).apparent()
ra, dec, distance = apparent.radec(epoch='date')
print(ra)
print(dec)

alt, az, distance = apparent.altaz()
print(alt, 'degrees above the horizon')
print(az, 'degrees around the horizon from north')

The output:

Traceback (most recent call last):
  File "/home/antoine/Documents/OmniLRS/src/stellar/stelar_engine.py", line 19, in <module>
    apparent = aristarchus.at(t).observe(earth).apparent()
  File "/home/antoine/.local/lib/python3.10/site-packages/skyfield/vectorlib.py", line 91, in at
    p, v, gcrs_position, message = self._at(t)
  File "/home/antoine/.local/lib/python3.10/site-packages/skyfield/vectorlib.py", line 215, in _at
    p2, v2, _, message = vf._at(t)
  File "/home/antoine/.local/lib/python3.10/site-packages/skyfield/planetarylib.py", line 232, in _at
    R, dRdt = self._frame.rotation_and_rate_at(t)
  File "/home/antoine/.local/lib/python3.10/site-packages/skyfield/planetarylib.py", line 169, in rotation_and_rate_at
    components, rates = self._segment.compute(t.whole, t.tdb_fraction, True)
  File "/home/antoine/.local/lib/python3.10/site-packages/jplephem/pck.py", line 144, in compute
    raise ValueError(
ValueError: segment only covers dates 2426-02-16 through 2650-01-25

I double-checked the time range for moon_pa_de440_200625.bpc and it should be fine. I'm not sure what's happening.
I also tried with moon_de440_200625.tf which resulted in the same outcome.

Thanks for all your work on this lib!

Best,

Antoine

[brandon-rhodes]

@AntoineRichard — I suspect that these ephemeris files, instead of having one segment per target, have multiple segments for the same target. Here's an open issue that might provide background:

#691

The ability to automatically handle these new cases is on my to-do list!

[AntoineRichard]

Thanks for clarifying that, I may stick with DE421 for the time being then!

[Tontyna]

It's not DE440, but the multi-segmented BPC file

moon_pa_de440_200625.bpc contains 2 segments for the moon:

segment #1: 1549-12-31 00:00:00.00 through 2426-02-16 00:00:00.00
segment #2: 2426-02-16 00:00:00.00 through 2650-01-25 00:00:00.00

Sadly PlanetaryConstants.read_binary() collects only the last one in its _segment_map

Consequently this will raise the ValueError
t = ts.utc(2019, 12, 20)

No Problem with this one (supposing your DExxx covers that date)
t = ts.utc(2500, 12, 20)

ugly workaround

Force the moon's segment to be the one you want.
Of course the body-number and the list-number depend on the bcp and tf files you already loaded and the framename you select.

With DE440 being the planetary ephemeris this will do:

pc = PlanetaryConstants()
pc.read_text(load('pck00011.tpc'))
pc.read_text(load('moon_de440_200625.tf')) # der zugehörige  NAIF lunar frame kernel
pc.read_binary(load('moon_pa_de440_200625.bpc'))

# ugly workaround
pc._segment_map[31008] = pc._segment_list[0]

frame = pc.build_frame_named('MOON_ME_DE440_ME421')

aristarchus = moon + pc.build_latlon_degrees(frame, 26.3, -46.8)

t = ts.utc(2019, 12, 20, 11, 5)

apparent = aristarchus.at(t).observe(earth).apparent()
ra, dec, distance = apparent.radec(epoch='date')

[Tontyna]

BTW: _segment_list was introduced in skyfield version 1.49

See #952