Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
Bernmeister opened this issue Apr 28, 2024 · 1 comment

Comments

@Bernmeister
Copy link

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
Copy link

Tontyna commented Jul 23, 2024

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,)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants