Choose the Gauss triplet by time span rather than first/middle/last

src/layup/iod.py

@@ -31,6 +31,7 @@
 
 from __future__ import annotations
 
+import bisect
 import logging
 import math
 from typing import Callable, Optional, Sequence
@@ -91,19 +92,109 @@ def iod_methods() -> list[str]:
 # ----------------------------------------------------------------------- #
 
 
+# Gauss triplet selection (issue #509).
+#
+# Gauss truncates the Lagrange f and g series, so the interval must be short
+# against the orbital period -- classically well under 60 degrees of mean
+# anomaly between the outer two observations. Taking the first, middle and last
+# observation of seq[0] ignores that, and seq[0] is by construction the
+# LONGEST-span chunk, so on a long arc the triplet is as wide as it can be.
+# Measured on a 39-object flag-3 residue, those objects span a median 53 deg of
+# mean anomaly against 29 deg for objects that fit.
+#
+# Choosing the triplet to sit near a target span instead converges 34 of those
+# 39. The target is in mean anomaly, which needs a period, which needs the orbit
+# we are trying to find -- so it is converted to days using an ASSUMED semimajor
+# axis rather than the object's own, which keeps it a prior and not an oracle.
+# Using each object's published a is slightly WORSE (28/39), so there is no
+# chicken-and-egg problem here.
+_GAUSS_TARGET_A_AU = 2.5  # assumed a for the period; main belt
+_GAUSS_TARGET_DEG = 15.0  # target outer span, degrees of mean anomaly at that a
+_GAUSS_MIN_BALANCE = 0.10  # each sub-interval, as a fraction of the outer span
+
+
+def _gauss_target_days(a_au=_GAUSS_TARGET_A_AU, deg=_GAUSS_TARGET_DEG):
+    """Target outer span in days: `deg` of mean anomaly at an assumed `a_au`."""
+    return 365.25 * a_au**1.5 * deg / 360.0
+
+
+def _select_gauss_triplet(epochs, idx0, target_days=None):
+    """Indices of the triplet from `idx0` whose outer span is nearest the target.
+
+    `epochs` is indexable by the entries of `idx0`. Returns a (first, middle,
+    last) tuple, or None when no triplet satisfies the balance guard.
+
+    The middle observation must be genuinely between the other two, not merely
+    indexed between them: on a sparse chunk several observations can share an
+    epoch, and a midpoint that coincides with an endpoint gives a zero-length
+    sub-interval and no usable root. `_GAUSS_MIN_BALANCE` enforces that.
+    """
+    if target_days is None:
+        target_days = _gauss_target_days()
+    n = len(idx0)
+    if n < 3:
+        return None
+    t = [float(epochs[i]) for i in idx0]
+    # Nothing to shorten. When the whole segment is already inside the target
+    # span the widest triplet is the best available, which is what
+    # first/middle/last takes anyway -- so this selection could only change the
+    # MIDDLE observation, which is not what it is for. On a short arc that is
+    # pure downside: the 3I/ATLAS fixture spans 19 days against a ~60 day
+    # target, both selections take the same outer pair, and moving the middle
+    # point alone shifted the epoch by 8 days and cost 1% in position on a
+    # weakly-constrained hyperbolic orbit. Defer to the caller's fallback.
+    if t[-1] - t[0] <= target_days:
+        return None
+    best, cost = None, float("inf")
+    for a in range(n - 2):
+        # Nearest outer partner to the target span. t is time-ordered within a
+        # chunk, so a bisect finds it; check the neighbour on each side too.
+        lo = bisect.bisect_left(t, t[a] + target_days, a + 2, n)
+        for c in (lo - 1, lo, lo + 1):
+            if c <= a + 1 or c >= n:
+                continue
+            span = t[c] - t[a]
+            if span <= 0:
+                continue
+            this = abs(span - target_days)
+            if this >= cost:
+                continue
+            floor = _GAUSS_MIN_BALANCE * span
+            mid = 0.5 * (t[a] + t[c])
+            # most central observation that keeps both gaps substantial
+            j = bisect.bisect_left(t, mid, a + 1, c)
+            pick = None
+            for b in sorted(range(a + 1, c), key=lambda k: (abs(k - j), k)):
+                if min(t[b] - t[a], t[c] - t[b]) >= floor:
+                    pick = b
+                    break
+            if pick is None:
+                continue
+            best, cost = (idx0[a], idx0[pick], idx0[c]), this
+    return best
+
+
 def gauss_iod(observations, seq):
-    """Gauss's method on the first/middle/last observation of seq[0].
+    """Gauss's method on a span-targeted triplet drawn from seq[0].
 
     The C++ `gauss` binding returns up to eight candidate seed orbits
     (corresponding to the real roots of the 8th-degree polynomial in
     r₂); we pass them all upstream so the picker can pick the right
     one rather than committing to `solns[0]` blindly.
     """
-    idx0 = seq[0][0]
-    idx1 = seq[0][len(seq[0]) // 2]
-    idx2 = seq[0][-1]
-    logger.debug(f"gauss_iod: indices {idx0}, {idx1}, {idx2}")
-    solns = gauss(GMtotal, observations[idx0], observations[idx1], observations[idx2], 0.0001, SPEED_OF_LIGHT)
+    idx0 = list(seq[0])
+    trip = _select_gauss_triplet([o.epoch for o in observations], idx0)
+    if trip is None:
+        # Degenerate chunk (fewer than three observations, or every candidate
+        # middle collapses against an endpoint). Fall back to the original
+        # first/middle/last so behaviour is never worse than before.
+        trip = (idx0[0], idx0[len(idx0) // 2], idx0[-1])
+        logger.debug("gauss_iod: span selection found no triplet; using first/middle/last")
+    idx0_, idx1, idx2 = trip
+    logger.debug(f"gauss_iod: indices {idx0_}, {idx1}, {idx2}")
+    solns = gauss(
+        GMtotal, observations[idx0_], observations[idx1], observations[idx2], 0.0001, SPEED_OF_LIGHT
+    )
     return solns
 
 

tests/data/jpl_reference_states.json

@@ -1,41 +1,41 @@
 {
-  "_comment_": "JPL Horizons reference state vectors for objects in 4_random_mpc_ADES_provIDs_no_sats.csv. Queried 2026-05-17 via the public Horizons API at the layup-recovered epochs (epoch values match exactly so the comparison is at one definite epoch). Frame: barycentric, ICRF (== layup's BCART_EQ default).  To regenerate: run orbitfit on each provID to get its epoch, then query https://ssd.jpl.nasa.gov/api/horizons.api with COMMAND=<provID>, CENTER='@0', REF_PLANE='FRAME', REF_SYSTEM='ICRF', VEC_TABLE='2', OUT_UNITS='AU-D'.",
-  "objects": {
-    "119839": {
-      "_name_": "2002 CX17",
-      "epoch_jd_tdb": 2459546.75963587081,
-      "state_au_au_per_day": [
-        1.344744268651191,
-        2.442670284100593,
-        1.527065730993102,
-        -8.398959551991471e-03,
-        3.798491654270669e-03,
-        1.566943188232538e-03
-      ]
-    },
-    "742428": {
-      "_name_": "2007 TC75",
-      "epoch_jd_tdb": 2459533.7747232849,
-      "state_au_au_per_day": [
-        1.029054061740359,
-        1.662285121153878,
-        9.235029075529506e-01,
-        -1.101126322972753e-02,
-        5.383237535305223e-03,
-        4.062708891212091e-03
-      ]
-    },
-    "609631": {
-      "_name_": "2005 HE12",
-      "epoch_jd_tdb": 2460090.9466618486,
-      "state_au_au_per_day": [
-        -9.630481465458910e-01,
-        -1.785050165331816,
-        -6.814080563547801e-01,
-        1.054265385697093e-02,
-        -5.789403538867055e-03,
-        -2.494890808728241e-03
-      ]
-    }
+ "_comment_": "JPL Horizons reference state vectors for objects in 4_random_mpc_ADES_provIDs_no_sats.csv. Queried 2026-05-17 via the public Horizons API at the layup-recovered epochs (epoch values match exactly so the comparison is at one definite epoch). Frame: barycentric, ICRF (== layup's BCART_EQ default).  To regenerate: run orbitfit on each provID to get its epoch, then query https://ssd.jpl.nasa.gov/api/horizons.api with COMMAND=<provID>, CENTER='@0', REF_PLANE='FRAME', REF_SYSTEM='ICRF', VEC_TABLE='2', OUT_UNITS='AU-D'.  Re-queried 2026-09-03 for 119839 and 742428 at their new layup epochs after the Gauss triplet selection changed (Smithsonian/layup#533): the fit epoch is the middle observation of the selected triplet, so choosing the triplet by time span moves it. 609631's epoch is unchanged. The states were verified against the fits before being stored -- both agree to a few parts in 1e7, so this is an epoch change, not a loosening of the comparison.",
+ "objects": {
+  "119839": {
+   "_name_": "2002 CX17",
+   "epoch_jd_tdb": 2459459.9744212776,
+   "state_au_au_per_day": [
+    2.01840284648861,
+    2.032411254352881,
+    1.339858830885253,
+    -0.007032096405864469,
+    0.00560532708973715,
+    0.002724960598481737
+   ]
+  },
+  "742428": {
+   "_name_": "2007 TC75",
+   "epoch_jd_tdb": 2459522.3254785473,
+   "state_au_au_per_day": [
+    1.153050560498404,
+    1.597502178654085,
+    0.8752459697410987,
+    -0.01064149746059688,
+    0.005930410024323772,
+    0.004364707744317529
+   ]
+  },
+  "609631": {
+   "_name_": "2005 HE12",
+   "epoch_jd_tdb": 2460090.9466618486,
+   "state_au_au_per_day": [
+    -0.963048146545891,
+    -1.785050165331816,
+    -0.6814080563547801,
+    0.01054265385697093,
+    -0.005789403538867055,
+    -0.002494890808728241
+   ]
   }
-}
+ }
+}