Osculating elements returning periapsis time far in the past

[jurezakrajsek]

Calculating orbital elements with Osculating_elements_of a comet returns periapsis time far in the past, which in some cases would return an error due to python date-time limitations.

Comet C/2021 P4 (ATLAS) orbital elements used: CK21P040 2022 07 30.4017 1.080556 0.996483 175.8228 348.0963 56.3118 20211117 9.5 4.0 C/2021 P4 (ATLAS) MPEC 2021-UA9

Calculated periapsis time: print(elements.periapsis_time)

print(elements.periapsis_time.utc_datetime())
Traceback (most recent call last):
  File ".\update_jpl_orbit.py", line 218, in <module>
    main(sys.argv[1:])
  File ".\update_jpl_orbit.py", line 75, in main
    'periapsis_time_utc': ele.periapsis_time.utc_datetime(),
  File "C:\Python37\lib\site-packages\skyfield\timelib.py", line 505, in utc_datetime
    dt, leap_second = self.utc_datetime_and_leap_second()
  File "C:\Python37\lib\site-packages\skyfield\timelib.py", line 538, in utc_datetime_and_leap_second
    dt = datetime(year, month, day, hour, minute, second, micro, utc)
ValueError: year -3363 is out of range

My script:

def preiapsis_time():
    local = pytz.timezone('UTC')
    today = datetime.date.today()
    date_now = datetime.datetime(today.year, today.month, today.day)
    date_now = local.localize(date_now)
    ecliptic = inertial_frames['ECLIPJ2000']

    with load.open(mpc.COMET_URL) as f:
        comets = mpc.load_comets_dataframe(f)

    print(len(comets), 'comets loaded')
    comets = (comets.sort_values('reference')
              .groupby('designation', as_index=False).last()
              .set_index('designation', drop=False))

    # Sample lookups.
    row = comets.loc['C/2021 P4 (ATLAS)']

    ts = load.timescale()
    t1 = ts.utc(date_now)
    eph = load('de421.bsp')
    sun, earth = eph['sun'], eph['earth']

    comet = sun + mpc.comet_orbit(row, ts, GM_SUN)
    position = (comet - sun).at(t1)
    elements = osculating_elements_of(position, ecliptic)
    print(elements.periapsis_time)
    print(elements.periapsis_time.utc_datetime())

It seams the osculating_elements_of is searching only for the first periapsis time in the past, while it would be nice that it would return the closest periapsis time to the given date, even if it is in the near future (in case of comet 2021P4 i would expect a date 2022-07-30 9:38:26.88) as the comet is currently heading to the sun and this periapsis time is the most relevant. It could also return an array of calculated periapsis times (for periodic comets, the first one being the most relevant one and than a few in the past) where the user could then select which periapsis time to print out).

Clear skies Jure

[brandon-rhodes]

That's a good question. I have just spent some time in the periapsis calculation, and it would make it slower if it did something fancy like try to guess which periapsis would be most convenient for the user. Right now it does something that is faster to compute and conceptually simpler (if I'm reading it correctly): it chooses the most recent periapsis before the time .t of the position that it's been given.

1. We should document that choice, so that users know which periapsis to expect when they access that attribute.
2. The documentation can also go ahead and illustrate how to generate further periapsides by adding or subtracting multiples of .period_in_days.
3. I think the documentation should also warn people away from using Python datetimes, to avoid the problem you're running into. They simply aren't safe for many astronomical dates, which can range outside of their allowed limits. Can we talk about the limitations in Skyfield's Time class that led you to try using the underpowered and awkwardly constrained Python datetime class? Maybe we can improve the Time object to do what you needed from datetime.

[jurezakrajsek]

Hi Brandon,

Yes the periapsis time returned is always before the current time, which might not be the relevant one when looking at the given .t

Addition to my above script

    t2 = ts.tt_jd(elements.periapsis_time.ut1 + elements.period_in_days)
    print(t2.utc_datetime())

I will need to do some additional logic to find the closest periapsis time to the given date, so it should not be more than 1/2 of period_in_days from the given time. That would then be my starting point from which I can then calculate periapsis times back to the past as many times that would be needed for short period comets.

Perhaps it would be more convenient if the periapsis_time would return the closest periapsis to the given .t? So it would look in a time frame -1/2 period - +1/2 period and return a periapsis_time that falls in this time frame.

Clear skies

[brandon-rhodes]

Perhaps it would be more convenient if the periapsis_time would return the closest periapsis to the given .t? So it would look in a time frame -1/2 period - +1/2 period and return a periapsis_time that falls in this time frame.

My guess is that Skyfield should continue doing something very simple — returning the most recent periapsis before the epoch of the elements — instead of trying something more complex that a particular user might or might not want. Among other things, one problem with extrapolating to a future periapsis is that in some cases the object won't exist by then, as with some sun-grazing comets.

The commit shown just above this comment shows how I've expanded the documentation to try addressing this issue, by showing how periapsis dates work. Hopefully that gets folks started who run into a similar issue.

I'm going to close this in the hopes that the new example code helps folks like you who need to aim for a closer periapsis. But if you can think of ways that the documentation can be expanded further, simply comment or re-open, and we'll try to get more written up!