Fix propagate() dropping the orbit dimension for batch orbits at one shared time

[apos]

Current Status

This is a regression introduced in v1.51 by fb12074 ("Fix #959: Kepler orbits with
time arrays of len==1"), which replaced squeeze(position_prop) with a hand-computed
output_shape. That commit correctly fixed the len-1 time-array case, but broke the
case of several orbits propagated to one shared time. Confirmed still present on
master (as of this PR) and in the latest released version, 1.54.

Bug Description

_KeplerOrbit._at() / propagate() raises

ValueError: cannot reshape array of size 3*N into shape (3,)

whenever N>1 orbits that were built together as one batch _KeplerOrbit (e.g. via
skyfield.data.mpc._comet_orbits()) are propagated to a single shared Time (a 0-d
t1).

We hit this in PiFinder, an open-source
telescope-pointing device, in
_calc_comets_vectorized()
(python/PiFinder/comets.py, added in commit 194c9b5e / PR #470, "Vectorize comet
propagation to fix locked-on-target CPU hog"). The relevant call is exactly:

# python/PiFinder/comets.py, lines 254–255
kepler = mpc._comet_orbits(comets_df, sf_utils.ts, GM_SUN)
helio_pos = kepler._at(t)[0]

comets_df holds every active comet from the MPC comets file (currently ~965 rows);
t is a single Time for "now". This batches all comets into one _KeplerOrbit and
propagates them all to that one shared moment in a single call, replacing a per-comet
Python loop that pegged a CPU core continuously whenever PiFinder was locked on a
target. This exact call pattern — many orbits, one shared observation time — triggers
the bug on every call.

Notably, PiFinder pins skyfield==1.45 in its own requirements.txt, so stock
PiFinder never hits this (1.45 still used squeeze()). It surfaces for anyone who
upgrades past 1.51 — e.g. our downstream fork, which unpins skyfield for unrelated
NumPy 2.0 compatibility reasons and ends up on 1.54. PiFinder's own code already
wraps the vectorized call in a try/except that falls back to the old per-comet loop
on any exception, so this doesn't crash — but it silently reinstates the exact
CPU-hogging behavior PR #470 was written to eliminate, which is how we found it (a
"why is the UI laggy again" investigation).

Explanation of the Solution

output_shape = (3,) + t1.shape is computed at the very top of propagate(), before
position/velocity are normalized to their working shape. When t1 is 0-d (one
shared time), output_shape becomes (3,) — but the position/velocity arrays
computed later in the function actually have shape (3, N, 1) for N orbits. The final
in-place position_prop.shape = output_shape then fails, since the sizes (3*N vs.
3) don't match.

Solution

Recover the true orbit count N from position.shape[1] — already available at that
point, before the ndim==1 normalization runs — and only widen output_shape to
(3, N) for this batch/shared-time case. The existing single-orbit behavior (squeeze
to (3,)) is unchanged.

Added a regression test, test_kepler_shape_with_multiple_orbits_and_single_shared_time,
next to the existing test_kepler_shape_with_time_of_length_one. It builds a
multi-orbit _KeplerOrbit directly via _from_periapsis() (no network/data files
needed) and propagates it to one shared scalar time, plus checks the single-orbit case
still squeezes to (3,).

skyfield/tests/test_keplerlib.py and skyfield/tests/test_elementslib.py (the two
test modules exercising propagate()/_KeplerOrbit) pass in full, 186 further tests
across the rest of the suite that don't require network-fetched ephemeris/star data
also pass unchanged.