skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.41k stars 212 forks source link

Planetary Reference Frame ignores multiple segments for one body #952

Open brandon-rhodes opened 6 months ago

brandon-rhodes commented 6 months ago

Discussed in https://github.com/skyfielders/python-skyfield/discussions/950

Originally posted by **YingChangJ** March 22, 2024 Hello, I am playing with moon Planetary Reference Frames function. The demo code in documentation, ``` from skyfield.api import PlanetaryConstants, load ts = load.timescale() t = ts.utc(2019, 12, 20, 11, 5) eph = load('de421.bsp') earth, moon = eph['earth'], eph['moon'] pc = PlanetaryConstants() pc.read_text(load('moon_080317.tf')) pc.read_text(load('pck00008.tpc')) pc.read_binary(load('moon_pa_de421_1900-2050.bpc')) frame = pc.build_frame_named('MOON_ME_DE421') aristarchus = moon + pc.build_latlon_degrees(frame, 26.3, -46.8) apparent = earth.at(t).observe(aristarchus).apparent() ra, dec, distance = apparent.radec(epoch='date') print(ra) print(dec) ``` This works great, and then I looked quickly at the readme in jpl website: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/aareadme.txt Binary Lunar PCKs has a 440 version with an approximate time coverage of 1550-2650. So I made some change to the demo code. ``` from skyfield.api import PlanetaryConstants, load 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 = earth.at(t).observe(aristarchus).apparent() ra, dec, distance = apparent.radec(epoch='date') print(ra) print(dec) ``` there is an error: `segment only covers dates 13447252800.0 through 20514081600.0` which is 2396-2-16 to 2620-1-25. (it should be from 1550 to 2650) Btw, a quick check of the range ``` print(ephe) SPICE kernel file 'de440.bsp' has 14 segments JD 2287184.50 - JD 2688976.50 (1549-12-30 through 2650-01-24) ``` And the time range of `moon_pa_de440_200625.bpc` is claimed in https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/moon_pa_de440_200625.cmt Thanks a lot YCH
brandon-rhodes commented 6 months ago

This is an interesting new case. Previously, all the PCK files that I had encountered had one segment per body, and so Skyfield was written with that simple case in mind. But now it looks like PCK files are being issued which have multiple segments for the same body. The commit that I've just made (which should appear right above this comment) lets us for the first time actually see what segments the object has been asked to load:

print(pc)
PlanetaryConstants
  559 key-values in .variables dict
  2 segments loaded
    2287184.50..2607184.50 frame=1  Unknown Body (31008)
    2607184.50..2688976.50 frame=1  Unknown Body (31008)

Since Skyfield stores the segments in a dictionary {object_id: segment}, the second segment is tossing the first one out of the dictionary.

NumPy doesn't provide any way to transparently treat two arrays in memory as a single array, so it's not clear how to move forward here. I will draw up some possible approaches so they can be compared.