From ea9b2b1c25144bb6db30cb2597d709e937e9dccc Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 08:45:52 -0400 Subject: [PATCH 1/3] Choose the Gauss triplet by time span, not first/middle/last (#509) Gauss truncates the Lagrange f and g series, so the interval has to be short against the orbital period. The triplet was the first, middle and last observation of seq[0], and seq[0] is by construction the longest-span chunk, so on a long arc it was as wide as it could be. The objects that fail to converge cold span a median 53 degrees of mean anomaly against 29 for those that fit. Pick the triplet whose outer span is nearest a target instead. The target is in mean anomaly, which needs a period, which needs the orbit being sought -- so it is converted to days with an ASSUMED a = 2.5 au rather than the object's own, about 60 days. That keeps it a prior and not an oracle; using each object's published a is slightly worse. A balance guard requires each sub-interval to be at least 10% of the outer span. Without it a repeated epoch can put the middle observation on an endpoint, which gives a zero-length interval and no usable root. Where no triplet qualifies the old first/middle/last choice is used, so behaviour is never worse than before. Measured cold on objects from the MPC catalog: regression, objects the old selection already fit 398/400 (99.50%) residue, objects it did not 39/46 (84.8%) Net 39 gained against 2 lost. The orbits are the same ones -- median |da/a| against the old selection is 1.7e-4 -- so this is a convergence improvement rather than a different answer. Closes #509. --- src/layup/iod.py | 91 ++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 85 insertions(+), 6 deletions(-) diff --git a/src/layup/iod.py b/src/layup/iod.py index 46210857..980d72f0 100644 --- a/src/layup/iod.py +++ b/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,97 @@ 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] + 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 From 77d6597eb650fe120a8e85a163446117ce0dd1ea Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 13:34:00 -0400 Subject: [PATCH 2/3] Fall back to first/middle/last when the arc is shorter than the target The span selection exists to shorten over-wide triplets on long arcs. When the whole segment already fits inside the ~60 day target there is nothing to shorten -- the widest triplet is the best available, which is what first/middle/last takes anyway -- so the only thing the selection could change is 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 the 60 day target, so both selections take the same outer pair; moving the middle point alone shifted the fit epoch by 8 days and cost 1% in position and 1.5% in velocity against JPL Horizons, on a weakly-constrained hyperbolic orbit. Caught by tests/layup/test_3i_atlas_validation.py, which now passes. The long-arc behaviour is unchanged, so the measured results stand: the regression sample is objects the old selection already fit, and the residue objects are long-arc by construction. --- src/layup/iod.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/layup/iod.py b/src/layup/iod.py index 980d72f0..f7b73709 100644 --- a/src/layup/iod.py +++ b/src/layup/iod.py @@ -135,6 +135,16 @@ def _select_gauss_triplet(epochs, idx0, target_days=None): 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 From 4f558d4cac54ae43ebfbc6781c0bf8d91c3c5d11 Mon Sep 17 00:00:00 2001 From: matthewholman Date: Thu, 3 Sep 2026 13:40:21 -0400 Subject: [PATCH 3/3] Regenerate the two JPL references whose epochs moved, and apply black The fit epoch is the middle observation of the selected triplet, so choosing the triplet by time span moves it. Two of the three real-data objects shift: 119839 by 87 days and 742428 by 11 days. 609631 and 3I/ATLAS are unchanged once the short-arc guard is in. Both references were re-queried from JPL Horizons at the new epochs, at the full precision the epoch assertion needs, using the recipe recorded in the file's own _comment_. The fits were checked against the new states BEFORE they were stored: 119839 agrees to 3.8e-7 and 742428 to 3.9e-7 in relative position. So this is an epoch change, not a loosened comparison -- which matters, because regenerating a reference is exactly how a real regression gets hidden. That is not hypothetical here: before the short-arc guard, 3I/ATLAS was off by 1% in position, and regenerating its reference would have made that pass silently. --- src/layup/iod.py | 4 +- tests/data/jpl_reference_states.json | 78 ++++++++++++++-------------- 2 files changed, 42 insertions(+), 40 deletions(-) diff --git a/src/layup/iod.py b/src/layup/iod.py index f7b73709..96070746 100644 --- a/src/layup/iod.py +++ b/src/layup/iod.py @@ -192,7 +192,9 @@ def gauss_iod(observations, seq): 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) + solns = gauss( + GMtotal, observations[idx0_], observations[idx1], observations[idx2], 0.0001, SPEED_OF_LIGHT + ) return solns diff --git a/tests/data/jpl_reference_states.json b/tests/data/jpl_reference_states.json index 68e93206..2b990601 100644 --- a/tests/data/jpl_reference_states.json +++ b/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=, 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=, 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 + ] } -} + } +} \ No newline at end of file