From ca4807940a1d00a4839c939329d207928b43cd68 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Mon, 5 Oct 2026 18:47:22 -0700 Subject: [PATCH] Convert observer geocentric offset from km to AU only once in geo_to_ssbary. Fixes OSSOS/liborbfit#19. Co-authored-by: Cursor --- orbfit/ephem_earth.c | 7 ++++--- test/test_orbit.py | 21 +++++++++++++++++++-- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index 84598c2..0c760bb 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -803,14 +803,15 @@ void geo_to_ssbary(double jd, double *x, double *y, double *z) geocenter_ssbary(jd, xgeo); /* fprintf(stderr, "XGEO: %lf %lf %lf\n", xgeo[0], xgeo[1], xgeo[2]); fprintf(stderr, "OBS: %lf %lf %lf (km)\n", *x, *y, *z); */ + /* R1.AU is in km */ *x /= R1.AU; *y /= R1.AU; *z /= R1.AU; /* fprintf(stderr, "OBS: %lf %lf %lf (au)\n", *x, *y, *z); */ /* Add in to give observatory coords in ICRS */ - *x = xgeo[0] + *x/R1.AU; - *y = xgeo[1] + *y/R1.AU; - *z = xgeo[2] + *z/R1.AU; + *x += xgeo[0]; + *y += xgeo[1]; + *z += xgeo[2]; /* fprintf(stderr, "PRO: %lf %lf %lf\n", *x, *y, *z); */ return ; } diff --git a/test/test_orbit.py b/test/test_orbit.py index 609d6bb..0e580d4 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -22,7 +22,9 @@ def test_parser(self): def test_fitradec(self): orbit = mp_ephem.BKOrbit(None, self.mpc_filename) orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.a.to('au').value, 17.08, 0) + self.assertAlmostEqual(orbit.a.to('au').value, 91.82, 0) + # JPL SBDB gives i=3.32 deg for 2020 KQ11 + self.assertAlmostEqual(orbit.inc.to('degree').value, 3.3, delta=0.5) class SimonFormat(unittest.TestCase): @@ -36,7 +38,7 @@ def test_parser(self): def test_fitradec(self): orbit = mp_ephem.BKOrbit(None, self.mpc_filename) orbit.predict(orbit.observations[0].date) - self.assertAlmostEqual(orbit.a.to('au').value, 43.8015219, 3) + self.assertAlmostEqual(orbit.a.to('au').value, 43.68817, 3) class OrbitFit(unittest.TestCase): @@ -152,6 +154,21 @@ def test_alphanumeric_obscode(self): self.assertAlmostEqual(orbit.a.to(units.AU).value, self.orbit.a.to(units.AU).value, 3) self.assertAlmostEqual(orbit.e.value, self.orbit.e.value, 4) + def test_geocentric_location_matches_obscode(self): + """ + Supplying the geocentric position (km) of Mauna Kea must give the same orbit as obscode 568. + """ + mko = coordinates.EarthLocation.from_geodetic(-(155 + 28/60. + 19.9/3600.) * units.degree, + (19 + 49/60. + 33.9/3600.) * units.degree, + 4207.24 * units.m) + observations = [obs for obs in self.observations if not obs.null_observation] + for obs in observations: + position, _ = mko.get_gcrs_posvel(obs.date) + obs.location = ObserverLocation(*position.xyz.to(units.km).value) + orbit = mp_ephem.BKOrbit(observations) + self.assertAlmostEqual(orbit.a.to(units.AU).value, self.orbit.a.to(units.AU).value, 3) + self.assertAlmostEqual(orbit.e.value, self.orbit.e.value, 4) + def test_predict_space_observatory(self): """ A space observatory without a supplied position is predicted from the geocenter.