skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 214 forks source link

Comets and Minor Planets throw ValueError with a narrow date in find_risings() #959

Open Bernmeister opened 7 months ago

Bernmeister commented 7 months ago

Using a date range of 24 hours to find risings (presumably too find settings ) for comets and minor planets throws a ValueError. See the example below. To test, comment/uncomment the t0 / t1 pairs.

The "fix" was simple: increase the date range to 48 hours.

#!/usr/bin/env python3

from skyfield import almanac
from skyfield.api import N, S, E, W, load, wgs84
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import mpc

ts = load.timescale()

eph = load('de421.bsp')
sun = eph['Sun']
earth = eph['earth']

bluffton = wgs84.latlon(40.8939 * N, 83.8917 * W)
observer = eph['Earth'] + bluffton

# FAILS
t0 = ts.utc(2024, 4, 26, 4)
t1 = ts.utc(2024, 4, 27, 4)

# PASSES
# t0 = ts.utc(2024, 4, 26, 4)
# t1 = ts.utc(2024, 4, 28, 4)

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

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

halley = sun + mpc.comet_orbit(comets.loc['1P/Halley'], ts, GM_SUN)

# ValueError: non-broadcastable output operand with shape (3,1) doesn't match the broadcast shape (3,3)
t, y = almanac.find_risings(observer, halley, t0, t1)
print(t.utc_iso(' '))
print(y)

with load.open('MPCORB.excerpt.DAT') as f:
    minor_planets = mpc.load_mpcorb_dataframe(f)

minor_planets = minor_planets.set_index('designation', drop=False)

ceres = sun + mpc.mpcorb_orbit(minor_planets.loc['(1) Ceres'], ts, GM_SUN)

# ValueError: non-broadcastable output operand with shape (3,1) doesn't match the broadcast shape (3,3)
t, y = almanac.find_risings(observer, ceres, t0, t1)
print(t.utc_iso(' '))
print(y)
Tontyna commented 4 months ago

Reason that the narrow dates fail is that there is only one rising inbetween.

Real reason is that KeplerOrbit._at(), more precisely propagate() in keplerlip.py, sqeezes the second dimension from its results' shape when fed with a time shaped (1,). Instead of position an velocity with expected shapes of (3,1) they have a shape of (3,)