brandon-rhodes / python-jplephem

Python version of NASA DE4xx ephemerides, the basis for the Astronomical Alamanac
MIT License
108 stars 28 forks source link

Unable to excerpt kernel with many segments #54

Closed miccoli closed 1 year ago

miccoli commented 1 year ago

I found this problem while trying to excerpt kernels over a time period that overlaps more than one segment for the same body, see #20, but I suspect thath the problem is more general.

Let me show the bug with an example using de441.bsp from https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp.

Structure of `de441.bsp` ``` $ python -m jplephem spk de441.bsp 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) ``` ``` $ python -m jplephem daf de441.bsp 1 DE-0441LE-0441 -479654827200.0 -960120000.0 299 2 1 2 8065 8076 2 DE-0441LE-0441 -479654827200.0 -960120000.0 199 1 1 2 8077 8088 3 DE-0441LE-0441 -479654827200.0 -960120000.0 399 3 1 2 8089 56797684 4 DE-0441LE-0441 -479654827200.0 -960120000.0 301 3 1 2 56797685 113587280 5 DE-0441LE-0441 -479654827200.0 -960120000.0 10 0 1 2 113587281 125707014 6 DE-0441LE-0441 -479654827200.0 -960120000.0 9 0 1 2 125707015 129169798 7 DE-0441LE-0441 -479654827200.0 -960120000.0 8 0 1 2 129169799 132632582 8 DE-0441LE-0441 -479654827200.0 -960120000.0 7 0 1 2 132632583 136095366 9 DE-0441LE-0441 -479654827200.0 -960120000.0 6 0 1 2 136095367 140077567 10 DE-0441LE-0441 -479654827200.0 -960120000.0 5 0 1 2 140077568 144579185 11 DE-0441LE-0441 -479654827200.0 -960120000.0 4 0 1 2 144579186 150639054 12 DE-0441LE-0441 -479654827200.0 -960120000.0 3 0 1 2 150639055 164836456 13 DE-0441LE-0441 -479654827200.0 -960120000.0 2 0 1 2 164836457 175917356 14 DE-0441LE-0441 -479654827200.0 -960120000.0 1 0 1 2 175917357 206389824 15 DE-0441LE-0441 -960120000.0 479387937600.0 299 2 1 2 206389825 206389836 16 DE-0441LE-0441 -960120000.0 479387937600.0 199 1 1 2 206389837 206389848 17 DE-0441LE-0441 -960120000.0 479387937600.0 399 3 1 2 206389849 263375588 18 DE-0441LE-0441 -960120000.0 479387937600.0 301 3 1 2 263375589 320361328 19 DE-0441LE-0441 -960120000.0 479387937600.0 10 0 1 2 320361329 332522922 20 DE-0441LE-0441 -960120000.0 479387937600.0 9 0 1 2 332522923 335997666 21 DE-0441LE-0441 -960120000.0 479387937600.0 8 0 1 2 335997667 339472410 22 DE-0441LE-0441 -960120000.0 479387937600.0 7 0 1 2 339472411 342947154 23 DE-0441LE-0441 -960120000.0 479387937600.0 6 0 1 2 342947155 346943109 24 DE-0441LE-0441 -960120000.0 479387937600.0 5 0 1 2 346943110 351460275 25 DE-0441LE-0441 -960120000.0 479387937600.0 4 0 1 2 351460276 357541074 26 DE-0441LE-0441 -960120000.0 479387937600.0 3 0 1 2 357541377 371787814 27 DE-0441LE-0441 -960120000.0 479387937600.0 2 0 1 2 371787815 382906986 28 DE-0441LE-0441 -960120000.0 479387937600.0 1 0 1 2 382906987 413484702 ```

In this file 1969-07-30T00:00 is the cutoff between the two segments. If I ask for the 1969/7/29 1969/7/31 time span I get

$ python -m jplephem excerpt 1969/7/29 1969/7/31 de441.bsp de441_partial.bsp 
Date 1969/7/29  = JD 2440432
Date 1969/7/31  = JD 2440434

'de441_partial.bsp' written successfully with the following contents

File type DAF/SPK and format LTL-IEEE with 3 segments:
2440432.50..2440448.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
2440432.50..2440448.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
2440432.50..2440440.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
$ python -m jplephem daf de441_partial.bsp 
 1 XE-0441LE-0441 -960120000.0 -958737600.0 3 0 1 2 9217 9261
 2 XE-0441LE-0441 -960120000.0 -958737600.0 2 0 1 2 9262 9297
 3 XE-0441LE-0441 -960120000.0 -959428800.0 1 0 1 2 9298 9345

28 segments are expected, but only the last 3 segments are present in the summary.

My conjecture (I had no time to do a complete debugging) is that the summary gets overwritten in DAF.add_array. See the following examples.

Excerpetd file with 24 segments is OK:

 $ python -m jplephem excerpt --target 299,199,399,301,10,9,8,7,6,5,4,3 1969/7/29 1969/7/31 de441.bsp de441_partial.bsp 
Date 1969/7/29  = JD 2440432
Date 1969/7/31  = JD 2440434

'de441_partial.bsp' written successfully with the following contents

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

Excerpted file with 26 expected segments is truncated:

$ python -m jplephem excerpt --target 299,199,399,301,10,9,8,7,6,5,4,3,2 1969/7/29 1969/7/31 de441.bsp de441_partial.bsp 
Date 1969/7/29  = JD 2440432
Date 1969/7/31  = JD 2440434

'de441_partial.bsp' written successfully with the following contents

File type DAF/SPK and format LTL-IEEE with 1 segments:
2440432.50..2440448.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)

Note: all examples run against commit 785c155c29 on MacOS.

sgbmzm commented 1 year ago

I also have the same problem. I am trying to create a file from 1.1.1. until 4.1.4001 Don't know why this is happening. I used to be able to create a file up to the year 3000 but now I can't in any way create a file up to the year 4001

brandon-rhodes commented 1 year ago

@miccoli @sgbmzm — Alas that I didn't quickly have time to return to this issue from April! But I've just dedicated a bit of my Saturday evening to it, and I think I found the problem:

3f6b678a04ebca33e1753400afc9071eabeaba54

Each time a summary record got full and a new one needed to be created, the code was zeroing out the summary count of the previous record, instead of leaving it at its current count — making the DAF format ignore all the summaries it contained. With the fix, I wait and don't reset my variable back to zero until I turn my attention to the new summary block I'm building.

To install the development version of this project and test this out, run:

pip install -U https://github.com/brandon-rhodes/python-jplephem/archive/master.zip

Let me know if this works, and I can release a new version!

sgbmzm commented 1 year ago

Thanks. I replaced only the DAF file with the corrected one. now i get File type DAF/SPK and format LTL-IEEE with 28 segments: (I really don't understand why there are 28 and not 14. It looks like 14 segments appear twice and together it's 28 segments)

I thought it was correct but when I try to use it on year 200 I get an error even though it should cover from year 1

skyfield.errors.EphemerisRangeError: ephemeris segment only covers dates 1969-07-29 23:59:18Z through 4001-01-05 23:58:51Z UT

So apparently there is still a problem with the final summary

python -m jplephem excerpt -- 1/1/1 4001/1/4 de441.bsp excerpt441.bsp
Date 1/1/1      = JD 1721426
Date 4001/1/4   = JD 3182399

'excerpt441.bsp' written successfully with the following contents

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)
1721424.50..2440432.50  Type 2  Earth Barycenter (3) -> Earth (399)
1721424.50..2440432.50  Type 2  Earth Barycenter (3) -> Moon (301)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Sun (10)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
1721424.50..2440432.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
2440400.50..8000016.50  Type 2  Venus Barycenter (2) -> Venus (299)
2440400.50..8000016.50  Type 2  Mercury Barycenter (1) -> Mercury (199)
2440432.50..3182400.50  Type 2  Earth Barycenter (3) -> Earth (399)
2440432.50..3182400.50  Type 2  Earth Barycenter (3) -> Moon (301)
2440432.50..3182400.50  Type 2  Solar System Barycenter (0) -> Sun (10)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Pluto Barycenter (9)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Neptune Barycenter (8)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Uranus Barycenter (7)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Saturn Barycenter (6)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Jupiter Barycenter (5)
2440432.50..3182416.50  Type 2  Solar System Barycenter (0) -> Mars Barycenter (4)
2440432.50..3182400.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)
2440432.50..3182400.50  Type 2  Solar System Barycenter (0) -> Venus Barycenter (2)
2440432.50..3182400.50  Type 2  Solar System Barycenter (0) -> Mercury Barycenter (1)
sgbmzm commented 1 year ago

Although from the size of the file (425 MB) it seems that it does contain the years 1 to 1969, so maybe it's just not updated that these years were included, which is why I get an error?

brandon-rhodes commented 1 year ago

@sgbmzm —

(I really don't understand why there are 28 and not 14. It looks like 14 segments appear twice and together it's 28 segments)

Do you see any segments with exactly the same dates? Or do 14 segments have early dates and 14 segments have later dates?

So apparently there is still a problem with the final summary

Happily, that's not a jplephem excerpt problem, but a Skyfield problem: I haven't yet had time to sit down and teach it about ephemeris files that have 2 segments for each planet instead of having 1 segment for each planet. Here's the issue:

https://github.com/skyfielders/python-skyfield/issues/691

Maybe I should add something to the documentation about how to work around it? Or maybe that time would be better spent working on logic for handling n segments as though they were a single segment.

sgbmzm commented 1 year ago

When I loaded it using skyfield you can see the times of each segment. As I copied below. Yes, it seems to me that there is still a problem with jplephem because I asked to receive information from the year 1 to the year 4001 but there are sections for example VENUS and MERCURY that I received starting from the year minus 13200 to plus 17191 even though I did not ask for most of these years.

[Admittedly, I received it in two parts: the first part until 1969 and the second part from 1969, which is also something I did not ask for. But this is because in the large original file de441.eph the information is divided into 14 sections multiplied by 2, a total of 28 sections, one up to 1969 and the other starting from 1969 (I don't know why they didn't summarize everything into 14 large sections)]

SPICE kernel file 'excerpt441.bsp' has 28 segments
  JD -3100015.50 - JD 2440432.50  (-13200-05-06 through 1969-07-29)
      2 -> 299  VENUS BARYCENTER -> VENUS
      1 -> 199  MERCURY BARYCENTER -> MERCURY
  JD 1721424.50 - JD 2440432.50  (0-12-30 through 1969-07-29)
      3 -> 399  EARTH BARYCENTER -> EARTH
      3 -> 301  EARTH BARYCENTER -> MOON
      0 -> 10   SOLAR SYSTEM BARYCENTER -> SUN
      0 -> 9    SOLAR SYSTEM BARYCENTER -> PLUTO BARYCENTER
      0 -> 8    SOLAR SYSTEM BARYCENTER -> NEPTUNE BARYCENTER
      0 -> 7    SOLAR SYSTEM BARYCENTER -> URANUS BARYCENTER
      0 -> 6    SOLAR SYSTEM BARYCENTER -> SATURN BARYCENTER
      0 -> 5    SOLAR SYSTEM BARYCENTER -> JUPITER BARYCENTER
      0 -> 4    SOLAR SYSTEM BARYCENTER -> MARS BARYCENTER
      0 -> 3    SOLAR SYSTEM BARYCENTER -> EARTH BARYCENTER
      0 -> 2    SOLAR SYSTEM BARYCENTER -> VENUS BARYCENTER
      0 -> 1    SOLAR SYSTEM BARYCENTER -> MERCURY BARYCENTER
  JD 2440400.50 - JD 8000016.50  (1969-06-27 through 17191-03-14)
      2 -> 299  VENUS BARYCENTER -> VENUS
      1 -> 199  MERCURY BARYCENTER -> MERCURY
  JD 2440432.50 - JD 3182400.50  (1969-07-29 through 4001-01-05)
      3 -> 399  EARTH BARYCENTER -> EARTH
      3 -> 301  EARTH BARYCENTER -> MOON
      0 -> 10   SOLAR SYSTEM BARYCENTER -> SUN
  JD 2440432.50 - JD 3182416.50  (1969-07-29 through 4001-01-21)
      0 -> 9    SOLAR SYSTEM BARYCENTER -> PLUTO BARYCENTER
      0 -> 8    SOLAR SYSTEM BARYCENTER -> NEPTUNE BARYCENTER
      0 -> 7    SOLAR SYSTEM BARYCENTER -> URANUS BARYCENTER
      0 -> 6    SOLAR SYSTEM BARYCENTER -> SATURN BARYCENTER
      0 -> 5    SOLAR SYSTEM BARYCENTER -> JUPITER BARYCENTER
      0 -> 4    SOLAR SYSTEM BARYCENTER -> MARS BARYCENTER
  JD 2440432.50 - JD 3182400.50  (1969-07-29 through 4001-01-05)
      0 -> 3    SOLAR SYSTEM BARYCENTER -> EARTH BARYCENTER
      0 -> 2    SOLAR SYSTEM BARYCENTER -> VENUS BARYCENTER
      0 -> 1    SOLAR SYSTEM BARYCENTER -> MERCURY BARYCENTER
brandon-rhodes commented 1 year ago

Yes, it seems to me that there is still a problem with jplephem because I asked to receive information from the year 1 to the year 4001 but there are sections for example VENUS and MERCURY that I received starting from the year minus 13200 to plus 17191 even though I did not ask for most of these years.

That's because those segments have only a single coefficient for the whole time period: a zero vector! That's because Venus and Mercury have no moons, and so the planet itself is the barycenter of its system.

Let's write a quick routine to unpack and display a Type 2 segment. It took me maybe an hour-ish to write this, and it should help with future debugging and questions; maybe I should put it inside a future jplephem release somewhere:

import numpy as np
from sys import argv
from jplephem.names import target_names
from jplephem.spk import SPK, titlecase

_component_names = {1: 'x', 2: 'x,y', 3: 'x,y,z'}

def _compute_calendar_date(jd_integer, julian_before=None):
    """Convert Julian day ``jd_integer`` into a calendar (year, month, day).

    Uses the proleptic Gregorian calendar unless ``julian_before`` is
    set to a specific Julian day, in which case the Julian calendar is
    used for dates older than that.

    """
    use_gregorian = (julian_before is None) or (jd_integer >= julian_before)

    # See the Explanatory Supplement to the Astronomical Almanac 15.11.
    f = jd_integer + 1401
    f += use_gregorian * ((4 * jd_integer + 274277) // 146097 * 3 // 4 - 38)
    e = 4 * f + 3
    g = e % 1461 // 4
    h = 5 * g + 2
    day = h % 153 // 5 + 1
    month = (h // 153 + 2) % 12 + 1
    year = e // 1461 - 4716 + (12 + 2 - month) // 12
    return year, month, day

def _format_date(jd):
    year, month, day = _compute_calendar_date(int(jd))
    return '{:4}-{:02}-{:02}'.format(year, month, day)

def print_type2_segment(segment):
    initial_epoch, interval_length, coefficients = segment.load_array()
    component_count, n, coefficient_count = coefficients.shape

    center = titlecase(target_names.get(segment.center, 'Unknown center'))
    target = titlecase(target_names.get(segment.target, 'Unknown target'))

    start_date = _format_date(segment.start_jd)
    end_date = _format_date(segment.end_jd)
    print(f'Date: {start_date} to {end_date}')
    print(f'Center: {center} ({segment.center})')
    print(f'Target: {target} ({segment.target})')

    c = _component_names.get(component_count, str(component_count))
    s = 8 * (segment.end_i - segment.start_i)
    print(f'Components: {c}  '
          f'Coefficient count: {coefficient_count}  '
          f'Size: {s:,} bytes')

    indent = '\n' + ' ' * 14
    if n > 5:
        ii = [0, 1, 2, '   ...', n-3, n-2, n-1]
    else:
        ii = range(n)

    with np.printoptions(precision=1, linewidth=82 - len(indent)):
        for i in ii:
            if isinstance(i, str):
                print(i)
                continue
            d = initial_epoch + i * interval_length
            nums = str(coefficients[:,i]).replace('\n', indent)
            print(_format_date(d).rjust(13), nums)

        i += 1
        d = initial_epoch + i * interval_length
        print(_format_date(d).rjust(13))

if __name__ == '__main__':
    kernel = SPK.open(argv[1])  # 'de441.bsp' or 'de441_partial.bsp'
    segment_index = int(argv[2])  # like '2' or '0'
    print_type2_segment(kernel.segments[segment_index])

I saved that to a file tmp54.py. Okay! Let's look at some segments in the original ephemeris. Here's a normal one:

python tmp54.py de441.bsp 2
Date: -13200-05-06 to 1969-07-29
Center: Earth Barycenter (3)
Target: Earth (399)
Components: x,y,z  Coefficient count: 13  Size: 454,316,760 bytes
 -13200-05-06 [[-3.4e+03 -1.4e+03  1.6e+02  1.0e+01 -5.3e-01 -5.9e-03  4.6e-05
                -9.3e-05  1.7e-06  1.6e-07 -6.0e-09  8.1e-10  2.9e-11]
               [ 3.1e+03 -1.3e+03 -1.5e+02  1.0e+01  4.2e-01 -1.4e-02  8.1e-04
                -2.1e-05 -4.9e-06  8.4e-08 -2.6e-09 -4.4e-10  6.0e-11]
               [ 1.1e+03 -5.7e+02 -4.9e+01  4.4e+00  1.4e-01 -5.7e-03  3.2e-04
                -1.3e-05 -1.9e-06  4.1e-08 -1.3e-09 -1.3e-10  2.5e-11]]
 -13200-05-10 [[-4.6e+03  2.3e+02  2.3e+02  4.3e-01 -7.9e-01 -2.4e-02 -8.0e-04
                 8.0e-05  1.1e-05  3.9e-07 -9.5e-09 -3.1e-09 -2.3e-10]
               [-1.9e+02 -1.9e+03  9.3e+00  1.5e+01  2.1e-01 -1.7e-02 -1.6e-03
                -1.3e-04 -6.9e-08  5.0e-07  4.2e-08  1.5e-09 -1.1e-10]
               [-2.9e+02 -7.2e+02  1.4e+01  6.0e+00  4.5e-02 -7.9e-03 -6.8e-04
                -4.9e-05  4.8e-07  2.2e-07  1.6e-08  4.5e-10 -5.5e-11]]
 -13200-05-14 [[-2.5e+03  1.8e+03  1.4e+02 -1.6e+01 -1.1e+00  2.4e-02  5.7e-03
                 2.1e-04 -2.4e-05 -2.6e-06  2.7e-08  2.0e-08  8.8e-10]
               [-3.3e+03 -1.1e+03  1.9e+02  1.3e+01 -7.2e-01 -7.4e-02 -7.6e-04
                 3.6e-04  2.5e-05 -1.1e-06 -2.2e-07 -4.3e-09  1.4e-09]
               [-1.4e+03 -3.3e+02  8.1e+01  4.2e+00 -3.3e-01 -2.8e-02 -3.3e-05
                 1.5e-04  8.5e-06 -5.5e-07 -8.6e-08 -7.6e-10  5.8e-10]]
   ...
   1969-07-17 [[ 4.5e+03 -2.7e+01 -2.3e+02 -2.4e+00  8.5e-01  2.5e-02  2.4e-04
                -1.9e-05 -6.3e-06 -5.2e-07 -1.1e-08  1.1e-09  1.4e-10]
               [-1.9e+02  1.8e+03  7.0e+00 -1.5e+01 -2.6e-01  2.2e-02  1.1e-03
                 7.9e-05  3.6e-06 -1.6e-07 -3.0e-08 -1.7e-09 -4.1e-11]
               [-4.5e+01  9.8e+02  1.0e+00 -8.1e+00 -1.3e-01  1.2e-02  6.1e-04
                 4.3e-05  1.9e-06 -9.4e-08 -1.7e-08 -9.1e-10 -2.0e-11]]
   1969-07-21 [[ 2.7e+03 -1.7e+03 -1.6e+02  1.5e+01  1.2e+00 -1.3e-02 -5.1e-03
                -3.3e-04  5.5e-06  2.7e-06  1.6e-07 -6.3e-09 -1.7e-09]
               [ 2.9e+03  1.1e+03 -1.8e+02 -1.3e+01  6.3e-01  7.2e-02  1.8e-03
                -2.2e-04 -2.9e-05 -7.6e-07  1.5e-07  1.7e-08  2.4e-10]
               [ 1.6e+03  5.8e+02 -9.7e+01 -7.1e+00  3.6e-01  3.9e-02  9.0e-04
                -1.3e-04 -1.6e-05 -3.8e-07  8.2e-08  9.0e-09  1.1e-10]]
   1969-07-25 [[-1.3e+03 -2.1e+03  9.0e+01  2.3e+01 -5.4e-01 -1.3e-01  1.9e-03
                 8.9e-04 -2.4e-06 -7.2e-06 -8.0e-08  6.3e-08  1.6e-09]
               [ 3.4e+03 -6.5e+02 -2.3e+02  6.8e+00  1.6e+00 -2.9e-02 -9.2e-03
                 7.8e-05  6.9e-05  3.2e-07 -5.9e-07 -1.1e-08  5.3e-09]
               [ 1.8e+03 -3.8e+02 -1.2e+02  4.0e+00  8.4e-01 -1.7e-02 -5.0e-03
                 5.4e-05  3.8e-05  8.3e-08 -3.2e-07 -5.3e-09  2.9e-09]]
   1969-07-29

You can see that it records positions for the Earth (relative to the Earth-Moon barycenter) that are spaced 4 days apart. That's a huge amount of data, of course, that takes many bytes to store.

What happens when jplephem excerpts it? Well, jplephem doesn't ever try to do the math to split one of these 4-day segments into smaller pieces; instead, given the dates you want, it either takes a whole 4-day piece or it discards it. So, for example, let's look at that same segment but in the excerpt:

python tmp54.py de441_partial.bsp 2
Date: 1969-07-25 to 1969-07-29
Center: Earth Barycenter (3)
Target: Earth (399)
Components: x,y,z  Coefficient count: 13  Size: 352 bytes
   1969-07-25 [[-1.3e+03 -2.1e+03  9.0e+01  2.3e+01 -5.4e-01 -1.3e-01  1.9e-03
                 8.9e-04 -2.4e-06 -7.2e-06 -8.0e-08  6.3e-08  1.6e-09]
               [ 3.4e+03 -6.5e+02 -2.3e+02  6.8e+00  1.6e+00 -2.9e-02 -9.2e-03
                 7.8e-05  6.9e-05  3.2e-07 -5.9e-07 -1.1e-08  5.3e-09]
               [ 1.8e+03 -3.8e+02 -1.2e+02  4.0e+00  8.4e-01 -1.7e-02 -5.0e-03
                 5.4e-05  3.8e-05  8.3e-08 -3.2e-07 -5.3e-09  2.9e-09]]
   1969-07-29

You asked for a starting date of 1969/7/29, but that fell at the end of a coefficient interval, so jplephem grabbed the whole 4-day block for you, going all the way back to 7/25.

Now let's look at the segment that is more confusing to you. In the original ephemeris:

python tmp54.py de441.bsp 0
Date: -13200-05-06 to 1969-07-29
Center: Venus Barycenter (2)
Target: Venus (299)
Components: x,y,z  Coefficient count: 2  Size: 88 bytes
 -13200-05-06 [[0. 0.]
               [0. 0.]
               [0. 0.]]
   1969-07-29

Wow! It's so little! Thirteen thousands years of "zero" for x, y, and z takes only 88 bytes.

When you asked for an excerpt, jplephem had only one choice: either keep this one coefficient interval, or discard it. Since it contained a date you asked for (1969/7/29), it kept it!

python tmp54.py de441_partial.bsp 0
Date: -13200-05-06 to 1969-07-29
Center: Venus Barycenter (2)
Target: Venus (299)
Components: x,y,z  Coefficient count: 2  Size: 88 bytes
 -13200-05-06 [[0. 0.]
               [0. 0.]
               [0. 0.]]
   1969-07-29

I suppose I should add something to the documentation about this, for folks who aren't familiar with SPK kernels and how NASA constructs them.

Also — hmm. It looks like jplephem is being a bit cautious: if you ask about date X, and date X is at the very end of a segment, then out of an abundance of caution, jplephem grabs that segment. But that is preventing you from getting as clean a result as you would like when, in a crazy multi-segment file like DE441, you ask for the specific date that ends one "era" and starts the next. Maybe I should tighten jplephem's logic up and make it a bit more daring?

Anyway: does the above output clear up your confusion?

sgbmzm commented 1 year ago

Wow! thank you very much! So now it's clear. I understand that there is no benefit in: 299 VENUS BARYCENTER -> VENUS and also in: 199 MERCURY BARYCENTER -> MERCURY

Because I just need to use only: 2 SOLAR SYSTEM BARYCENTER -> VENUS BARYCENTER 1 SOLAR SYSTEM BARYCENTER -> MERCURY BARYCENTER

brandon-rhodes commented 1 year ago

I understand that there is no benefit in: 299 VENUS BARYCENTER -> VENUS…

Yes, that's correct: you can just use the barycenter. Some folks avoid doing that — maybe they would find it confusing to build a list of targets that included some barycenters and some planets, or maybe it would just seem odd to them to be running a mission that orbits the Venus barycenter rather than just Venus itself in their code — but unless the use of the barycenter makes your code harder to read, it will make things slightly simpler and faster to just deal with one vector rather than two!

miccoli commented 1 year ago

@brandon it seems to mee that there has been a regression here: if I run from 9f06060 I get

python -m jplephem excerpt --target 3 1969/7/29 1969/7/31 de441.bsp de441_partial.bsp
Date 1969/7/29  = JD 2440432.5
Date 1969/7/31  = JD 2440434.5

'de441_partial.bsp' written successfully with the following contents

File type DAF/SPK and format LTL-IEEE with 1 segments:
2440432.50..2440448.50  Type 2  Solar System Barycenter (0) -> Earth Barycenter (3)

but 1969-07-29 00:00:00 UT is 2440431.50 JD, so I would expect two segments, not just one.

brandon-rhodes commented 1 year ago

@miccoli — Well, drat; thank you for pointing out that I added 0.5 when I should have subtracted!

Instead of tossing up a quick fix, I really should write a test for the excerpter, so that things like this stay fixed going forward, since now folks besides just myself seem to be using it! Let's see if by the end of the day I can have a fix up with a test to keep it working. In which case I'll signal you so you can try your test over again!

brandon-rhodes commented 1 year ago

@miccoli — There! I've added a test and made the fix of that off-by-one bug in the JD computation. Could you try again and let me know if the regression now looks fixed to you?

miccoli commented 1 year ago

For me it's OK, thanks.

brandon-rhodes commented 1 year ago

@miccoli — Great, I'll go ahead and plan a new release today or tomorrow. Thanks for pointing out the issue!