skyfielders / python-skyfield

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

Support SPICE ephemeris with multi-segments #975

Open Tontyna opened 3 months ago

Tontyna commented 3 months ago

resolves #691

MultiKernel can load multiple SpiceKernels and combines their segments. When a target's position is queried, the appropriate SPICESegments are picked from all segments available in all loaded SPK files.

Per definition that includes support of ephemerides that have more than 1 segment per planet.

Tontyna commented 3 months ago

SPK files for celestial small bodies tend to be multi-segmented. But ... they probably consist of type 21 segments which aren't supported in Skyfield (yet).

Tontyna commented 3 months ago

Test multi-segmented DE441

from skyfield.api import load
from skyfield.jpl_multikernel import MultiKernel

ts = load.timescale()
kernel = MultiKernel()
# it's possible to create ` MultiSegments`  in advance
earth = kernel['earth']
mars = kernel['mars']   # thanks to `MultiSegment._fallback_target`  this will work

kernel.load_kernel('de441.bsp')

# a time covered by first segment
t = ts.tdb_jd(1440432)
barycentric = earth.at(t)

# a time covered by second segment 
t = ts.tdb_jd(3440432)
barycentric = earth.at(t)

## critical tdb: obeserve() requires mars-vectors from both segments
t = ts.tdb_jd(2440432.501)
barycentric = earth.at(t)
astrometric = barycentric.observe(mars)
Tontyna commented 3 months ago

The fact that MultiSegment.center is detected on the fly instead of being fixed right from the beginning has its drawbacks.

E.g. when a MultiSegment is involved in a VectorSum (by addition or subtraction) Skyfield checks and combines the VectorFuntions' centers and targets. No problem with the target, but the VectorSum.center determines the ICRF (sub)class used in subsequent calculations.

Without major changes to existing code, the only solution I can think of, is:

  1. assume that the center for the target will always be the same in the time range covered by a given set of SPK files
  2. call MultiSegment.at() to initialize its center

For example topocentric coordinates:

from skyfield.api import load, N, W, wgs84
from skyfield.jpl_multikernel import MultiKernel

ts = load.timescale()
t = ts.now()

eph = MultiKernel()
eph.load_kernel('de441.bsp')

earth = eph['earth']
mars = eph['mars']

# !important: initialize earth's center before adding
earth.at(t)

boston = earth + wgs84.latlon(42.3583 * N, 71.0636 * W)
astrometric = boston.at(t).observe(mars)
alt, az, d = astrometric.apparent().altaz()

print(alt)
print(az)
brandon-rhodes commented 2 months ago

@Tontyna — I'm not sure I'll have time this month to respond in detail, but I at least want to say, thanks for your interest in getting Skyfield set up to handle ephemerides that have two or more segments per target; it's definitely a use case I want to see added to Skyfield this year. I'm not sure I'll want to merge a solution that uses a for loop over time, though, as that's something I've been careful to avoid so far? The questions you raise about how the classes should be designed is a very important one, though, regardless of what kind of loop powers the computation under the hood; I plan to read over your solution carefully as I think through what direction the classes should take to add this ability to Skyfield.

Tontyna commented 2 months ago

Of course, looping over time in a for-loop, isn't great. But what else is there to find the segment suitable for the requested time? Of course, a pick-the-right-segment_for(t) could be built into the SPICESegment, but somewhere somebody has to pick. And neither the plain SPICESegment nor the single-file SpiceKernel know how to build the heliocentric chain when the required data stem form multiple SPK files.

BTW.: My major intention for this MultiKernel wasn't the two-segmented DE441. My pet project is an astrological application, so I wanted a class that automatically builds the proper heliocentric VectorFunction chains for various asteroids (some of them in multi-planet-multi-segmented SPK files, some of them of Typ 21). Covering a time span as big as possible.

Tontyna commented 2 months ago

Anyway, there is no need to merge this single jpl_multikernel.py file into skyfield. If required, a plain local import will do.