From e9bced14679791104273e78a09e60deafab3ecca Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Sat, 26 Jul 2025 09:47:05 -0700 Subject: [PATCH 1/8] Add likelihood from CLASSY formated comment lines. --- src/mp_ephem/ephem.py | 10 ++++++++-- src/mp_ephem/time_mpc.py | 2 +- test/data/classy.mpc | 18 ++++++++++++++++++ test/test_orbit.py | 16 ++++++++++++++++ 4 files changed, 43 insertions(+), 3 deletions(-) create mode 100644 test/data/classy.mpc diff --git a/src/mp_ephem/ephem.py b/src/mp_ephem/ephem.py index 15867b8..94529d2 100644 --- a/src/mp_ephem/ephem.py +++ b/src/mp_ephem/ephem.py @@ -602,6 +602,9 @@ def from_string(cls, input_line): # noqa: C901 mpc_format = '0s5s7s1s1s1s17s12s12s9x5s1s6x3s' ossos_format1 = '1s4s7s1s1s1s17s12s12s9x5s1s6x3s' ossos_format2' = '1s0s11s1s1s1s17s12s12s9x5s1s6x3s' + classy_format = '1s2s11s1s8s3s1s2s8s1s8s1s4s1s1s1s5s1s4s1s4s' + 112123456789011123456781231112123456789912345 +'#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G if those fail then tries Alex Parker's .ast format. @@ -1076,6 +1079,8 @@ def __init__(self, version, frame, source_name, photometry_note, mpc_note, x, y, self._plate_uncertainty = None self._astrometric_level = 0 self.likelihood = likelihood + if self.likelihood is not None and not isinstance(self.likelihood, int): + self.likelihood = int(self.likelihood) try: self.plate_uncertainty = plate_uncertainty self.astrometric_level = astrometric_level @@ -1115,11 +1120,12 @@ def from_string(cls, comment): # noqa: C901 if len(values) > 1: comment_string = values[1].lstrip(' ') # O 1631355p21 O13AE2O Z 1632.20 1102.70 0.21 3 ----- ---- % Apcor failure. + # '#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G ossos_comment_format = '1s1x12s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' old_ossos_comment_format = '1s1x10s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' bad_ossos_comment_format = '1s3x10s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' orig_ossos_comment_format = '1s1x10s1x8s1x1s1s1x6s1x6s1x5s1x4s1x4s1x' - classy_comment_format = '1s1x10s1x8s1x1s1s1x6s1x6s1x5s1x4s1x4s1x4s1x' + classy_comment_format = '1s1x12s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x4s1x' for struct_ in [ossos_comment_format, old_ossos_comment_format, bad_ossos_comment_format, orig_ossos_comment_format, classy_comment_format]: @@ -1335,7 +1341,7 @@ def __str__(self): comm += self.to_str('{:5.2f}', self.mag, "-" * 5) comm += self.to_str('{:4.2f}', self.mag_uncertainty, "-" * 4) if self.likelihood is not None: - comm += self.to_str('{:4d}', int(numpy.floor(self.likelihood)), '-'*4) + comm += self.to_str('{:4d}', int(numpy.floor(float(self.likelihood))), '-'*4) comm += ' % {}'.format(self.comment) return comm diff --git a/src/mp_ephem/time_mpc.py b/src/mp_ephem/time_mpc.py index 077dfc6..562d3c3 100644 --- a/src/mp_ephem/time_mpc.py +++ b/src/mp_ephem/time_mpc.py @@ -110,7 +110,7 @@ def str_kwargs(self): """ scale = self.scale.upper().encode('ascii'), iys, ims, ids, ihmsfs = d2dtf(scale, self.precision, - self.jd1, self.jd2_filled) + self.jd1, self.jd2) # Get the str_fmt element of the first allowed output subformat _, _, str_fmt = self._select_subfmts(self.out_subfmt)[0] diff --git a/test/data/classy.mpc b/test/data/classy.mpc new file mode 100644 index 0000000..856b53e --- /dev/null +++ b/test/data/classy.mpc @@ -0,0 +1,18 @@ +#C AS222083100 FbAS2001 Z 1392.00 2932.00 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.35746 22 24 20.536-11 34 27.54 568 +#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.41996 22 24 20.246-11 34 28.99 568 +#C AS222083100 FbAS2001 Z 1345.95 2916.27 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.48246 22 24 19.956-11 34 30.45 568 +#C AS222090100 FbAS2001 Z 1374.00 2930.00 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.37455 22 24 16.049-11 34 52.70 568 +#C AS222090100 FbAS2001 Z 1351.30 2922.22 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.43705 22 24 15.763-11 34 54.14 568 +#C AS222090100 FbAS2001 Z 1328.59 2914.43 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.49955 22 24 15.478-11 34 55.58 568 +#C AS222090200 FbAS2001 Z 1374.00 2925.00 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.34810 22 24 11.737-11 35 16.36 568 +#C AS222090200 FbAS2001 Z 1352.35 2915.59 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.41060 22 24 11.464-11 35 18.10 568 +#C AS222090200 FbAS2001 Z 1330.70 2906.19 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.47310 22 24 11.192-11 35 19.84 568 \ No newline at end of file diff --git a/test/test_orbit.py b/test/test_orbit.py index 2c5f525..40e9b68 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -144,3 +144,19 @@ def test_tnodb_discovery_flags(self): orbit = mp_ephem.BKOrbit(None, ast_filename=os.path.join(__PATH__,'data/o4h29.ast')) for observation in orbit.observations: self.assertTrue(observation.discovery) + +class CLASSYFORM(unittest.TestCase): + + def setUp(self): + self.mpc_filename = os.path.join(__PATH__, 'data/classy.mpc') + + def test_parser(self): + obs = mp_ephem.EphemerisReader().read(self.mpc_filename) + self.assertIsInstance(obs[0], mp_ephem.Observation) + self.assertIsInstance(obs[0].comment, mp_ephem.OSSOSComment) + self.assertAlmostEqual(obs[0].comment.likelihood, 12) + + 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, 44.48, 2) # Example value, adjust as needed \ No newline at end of file From 317ce41163f7d5d771c379290f1a03e5876cfca5 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Sat, 26 Jul 2025 09:47:05 -0700 Subject: [PATCH 2/8] Add likelihood from CLASSY formated comment lines. --- src/mp_ephem/ephem.py | 10 ++++++++-- src/mp_ephem/time_mpc.py | 2 +- test/data/classy.mpc | 18 ++++++++++++++++++ test/test_orbit.py | 16 ++++++++++++++++ test/test_query.py | 2 +- 5 files changed, 44 insertions(+), 4 deletions(-) create mode 100644 test/data/classy.mpc diff --git a/src/mp_ephem/ephem.py b/src/mp_ephem/ephem.py index 15867b8..94529d2 100644 --- a/src/mp_ephem/ephem.py +++ b/src/mp_ephem/ephem.py @@ -602,6 +602,9 @@ def from_string(cls, input_line): # noqa: C901 mpc_format = '0s5s7s1s1s1s17s12s12s9x5s1s6x3s' ossos_format1 = '1s4s7s1s1s1s17s12s12s9x5s1s6x3s' ossos_format2' = '1s0s11s1s1s1s17s12s12s9x5s1s6x3s' + classy_format = '1s2s11s1s8s3s1s2s8s1s8s1s4s1s1s1s5s1s4s1s4s' + 112123456789011123456781231112123456789912345 +'#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G if those fail then tries Alex Parker's .ast format. @@ -1076,6 +1079,8 @@ def __init__(self, version, frame, source_name, photometry_note, mpc_note, x, y, self._plate_uncertainty = None self._astrometric_level = 0 self.likelihood = likelihood + if self.likelihood is not None and not isinstance(self.likelihood, int): + self.likelihood = int(self.likelihood) try: self.plate_uncertainty = plate_uncertainty self.astrometric_level = astrometric_level @@ -1115,11 +1120,12 @@ def from_string(cls, comment): # noqa: C901 if len(values) > 1: comment_string = values[1].lstrip(' ') # O 1631355p21 O13AE2O Z 1632.20 1102.70 0.21 3 ----- ---- % Apcor failure. + # '#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G ossos_comment_format = '1s1x12s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' old_ossos_comment_format = '1s1x10s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' bad_ossos_comment_format = '1s3x10s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x' orig_ossos_comment_format = '1s1x10s1x8s1x1s1s1x6s1x6s1x5s1x4s1x4s1x' - classy_comment_format = '1s1x10s1x8s1x1s1s1x6s1x6s1x5s1x4s1x4s1x4s1x' + classy_comment_format = '1s1x12s1x11s1x1s1s1x7s1x7s1x4s1x1s1x5s1x4s1x4s1x' for struct_ in [ossos_comment_format, old_ossos_comment_format, bad_ossos_comment_format, orig_ossos_comment_format, classy_comment_format]: @@ -1335,7 +1341,7 @@ def __str__(self): comm += self.to_str('{:5.2f}', self.mag, "-" * 5) comm += self.to_str('{:4.2f}', self.mag_uncertainty, "-" * 4) if self.likelihood is not None: - comm += self.to_str('{:4d}', int(numpy.floor(self.likelihood)), '-'*4) + comm += self.to_str('{:4d}', int(numpy.floor(float(self.likelihood))), '-'*4) comm += ' % {}'.format(self.comment) return comm diff --git a/src/mp_ephem/time_mpc.py b/src/mp_ephem/time_mpc.py index 077dfc6..562d3c3 100644 --- a/src/mp_ephem/time_mpc.py +++ b/src/mp_ephem/time_mpc.py @@ -110,7 +110,7 @@ def str_kwargs(self): """ scale = self.scale.upper().encode('ascii'), iys, ims, ids, ihmsfs = d2dtf(scale, self.precision, - self.jd1, self.jd2_filled) + self.jd1, self.jd2) # Get the str_fmt element of the first allowed output subformat _, _, str_fmt = self._select_subfmts(self.out_subfmt)[0] diff --git a/test/data/classy.mpc b/test/data/classy.mpc new file mode 100644 index 0000000..856b53e --- /dev/null +++ b/test/data/classy.mpc @@ -0,0 +1,18 @@ +#C AS222083100 FbAS2001 Z 1392.00 2932.00 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.35746 22 24 20.536-11 34 27.54 568 +#C AS222083100 FbAS2001 Z 1368.97 2924.14 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.41996 22 24 20.246-11 34 28.99 568 +#C AS222083100 FbAS2001 Z 1345.95 2916.27 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 08 31.48246 22 24 19.956-11 34 30.45 568 +#C AS222090100 FbAS2001 Z 1374.00 2930.00 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.37455 22 24 16.049-11 34 52.70 568 +#C AS222090100 FbAS2001 Z 1351.30 2922.22 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.43705 22 24 15.763-11 34 54.14 568 +#C AS222090100 FbAS2001 Z 1328.59 2914.43 0.20 0 ----- ---- 12 % G + FbAS2001 C2022 09 01.49955 22 24 15.478-11 34 55.58 568 +#C AS222090200 FbAS2001 Z 1374.00 2925.00 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.34810 22 24 11.737-11 35 16.36 568 +#C AS222090200 FbAS2001 Z 1352.35 2915.59 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.41060 22 24 11.464-11 35 18.10 568 +#C AS222090200 FbAS2001 Z 1330.70 2906.19 0.20 0 ----- ---- 13 % G + FbAS2001 C2022 09 02.47310 22 24 11.192-11 35 19.84 568 \ No newline at end of file diff --git a/test/test_orbit.py b/test/test_orbit.py index 2c5f525..40e9b68 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -144,3 +144,19 @@ def test_tnodb_discovery_flags(self): orbit = mp_ephem.BKOrbit(None, ast_filename=os.path.join(__PATH__,'data/o4h29.ast')) for observation in orbit.observations: self.assertTrue(observation.discovery) + +class CLASSYFORM(unittest.TestCase): + + def setUp(self): + self.mpc_filename = os.path.join(__PATH__, 'data/classy.mpc') + + def test_parser(self): + obs = mp_ephem.EphemerisReader().read(self.mpc_filename) + self.assertIsInstance(obs[0], mp_ephem.Observation) + self.assertIsInstance(obs[0].comment, mp_ephem.OSSOSComment) + self.assertAlmostEqual(obs[0].comment.likelihood, 12) + + 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, 44.48, 2) # Example value, adjust as needed \ No newline at end of file diff --git a/test/test_query.py b/test/test_query.py index 9508252..f6d74b9 100644 --- a/test/test_query.py +++ b/test/test_query.py @@ -7,4 +7,4 @@ class TestQuery(TestCase): def test_target(self): body = horizons.Body('2002 MS4') self.assertIsInstance(body, horizons.Body) - self.assertAlmostEqual(body.mag, 20.4, places=1) + self.assertAlmostEqual(body.mag, 20.3, places=1) From ddf3e423035986e9d64b8419b1f5ed5099016d2c Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Tue, 14 Apr 2026 09:21:18 -0700 Subject: [PATCH 3/8] Increase version number to track jd modification --- src/mp_ephem/__version__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mp_ephem/__version__.py b/src/mp_ephem/__version__.py index 252b2d3..ef3d3f2 100644 --- a/src/mp_ephem/__version__.py +++ b/src/mp_ephem/__version__.py @@ -1 +1 @@ -version = "0.13.0" +version = "0.14.0" From 1eae694b7c851698a1371da47a58ad2ef2264ded Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Tue, 25 Aug 2026 14:08:27 -0700 Subject: [PATCH 4/8] become compatiable with astropy 6 and higher --- setup.py | 2 +- src/mp_ephem/time_mpc.py | 197 +++++++++++++++------------------------ 2 files changed, 76 insertions(+), 123 deletions(-) diff --git a/setup.py b/setup.py index ec869a5..512734f 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ import sys from setuptools import setup, find_packages, Extension -version = "0.13.0" +version = "0.14.0" dependencies = ['astropy', 'pyerfa', 'numpy', 'requests'] diff --git a/src/mp_ephem/time_mpc.py b/src/mp_ephem/time_mpc.py index 562d3c3..ed13c0e 100644 --- a/src/mp_ephem/time_mpc.py +++ b/src/mp_ephem/time_mpc.py @@ -1,134 +1,87 @@ -import logging -from datetime import datetime -import numpy -import re -import time -from erfa import d2dtf, dtf2d -from astropy.time import TimeString - - -class TimeMPC(TimeString): - """ - Override the TimeString class to convert from MPC format string to astropy.time.Time object. +"""MPC calendar-date time format for :class:`astropy.time.Time`. - usage: +Registers format ``mpc`` so strings like ``YYYY MM DD.ffffff`` (where the +fraction is a fraction of a day) can initialize and format ``Time`` objects. - from time_mpc import TimeMPC - from astropy.time.core import Time +Subclassing :class:`~astropy.time.TimeString` is enough for Astropy 6+: the +base class owns JD conversion, masking, and output. This module only teaches +it that digits after the decimal are a **day** fraction (not a second +fraction, which is the ``TimeString`` default). +""" - t = Time('2000 01 01.00001', format='mpc', scale='utc') - - str(t) == '2000 01 01.000001' - """ - name = 'mpc' - subfmts = (('mpc', '%Y %m %d', "{year:4d} {mon:02d} {day:02d}.{fracday:s}"),) +import re - def __init__(self, val1, val2, scale, precision=6, - in_subfmt=None, out_subfmt=None, from_jd=False): - super(TimeMPC, self).__init__(val1=val1, - val2=val2, - scale=scale, - precision=precision, - in_subfmt=in_subfmt, - out_subfmt=out_subfmt, - from_jd=from_jd) - self.precision = precision +from astropy.time import TimeString - def set_jds(self, val1, val2): - """Parse the time strings contained in val1 and set jd1, jd2""" - # Select subformats based on current self.in_subfmt - subfmts = self._select_subfmts(self.in_subfmt) +# Importing this module registers the format via TimeString's metaclass. +__all__ = ["TimeMPC"] - iterator = numpy.nditer([val1, None, None, None, None, None, None], - op_dtypes=[val1.dtype] + 5 * [numpy.intc] + [numpy.double]) - for val, iy, im, i_day, ihr, i_min, d_sec in iterator: - iy[...], im[...], i_day[...], ihr[...], i_min[...], d_sec[...] = ( - self.parse_string(val.item(), subfmts)) +class TimeMPC(TimeString): + """ + Minor Planet Center calendar date: ``YYYY MM DD.ffffff``. + + The fractional part is a fraction of a **day** (MPC convention), not a + fraction of a second. + + Examples + -------- + >>> from astropy.time import Time # doctest: +SKIP + >>> import mp_ephem.time_mpc # registers format='mpc' # doctest: +SKIP + >>> t = Time('2000 01 01.00001', format='mpc', scale='utc', precision=5) + >>> t.iso + '2000-01-01 00:00:00.86400' + >>> t.mpc + '2000 01 01.00001' + """ - self.jd1, self.jd2 = dtf2d( - self.scale.upper().encode('utf8'), *iterator.operands[1:]) + name = "mpc" + subfmts = ( + ( + "mpc", + re.compile( + r"^\s*(?P\d{4})\s+(?P\d{1,2})\s+(?P\d{1,2})" + r"(?:\.(?P\d+))?\s*$" + ), + "{year:4d} {mon:02d} {day:02d}.{fracday:s}", + ), + ) def parse_string(self, timestr, subfmts): - """ - Read time from a single string, using a set of possible formats. - - :param timestr: A string representing the time. eg: '2001 01 01.00001' - :param subfmts: which format is the string in, eg: 'mpc' - """ - # Datetime components required for conversion to JD by ERFA, along - # with the default values. - components = ('year', 'mon', 'mday') - defaults = (None, 1, 1, 0) - # Assume that anything following "." on the right side is a - # floating fraction of a second. - try: - idot = timestr.rindex('.') - except Exception as ex: - logging.debug(ex) - fracday = 0.0 - else: - timestr, fracday = timestr[:idot], timestr[idot:] - fracday = float(fracday) - - for _, strptime_fmt_or_regex, _ in subfmts: - if isinstance(strptime_fmt_or_regex, str): - try: - tm = time.strptime(timestr, strptime_fmt_or_regex) - except ValueError as ex: - logging.debug(str(ex)) - continue - else: - vals = [getattr(tm, 'tm_' + component) - for component in components] - vals.append(tm.tm_hour + int(24 * fracday)) - vals.append(tm.tm_min + int(60 * (24 * fracday - vals[-1]))) - vals.append(tm.tm_sec + 60 * (60 * (24 * fracday - vals[-2]) - vals[-1])) - else: - tm = re.match(strptime_fmt_or_regex, timestr) - if tm is None: - continue - tm = tm.groupdict() - vals = [int(tm.get(component, default)) for component, default - in zip(components, defaults)] - - hrprt = int(24 * fracday) - vals.append(hrprt) - mnprt = int(60 * (24 * fracday - hrprt)) - vals.append(mnprt) - scprt = 60 * (60 * (24 * fracday - hrprt) - mnprt) - vals.append(scprt) - return vals - else: - raise ValueError('Time {0} does not match {1} format' - .format(timestr, self.name)) + """Parse one MPC date string into ``(year, mon, mday, hour, min, sec)``.""" + for _, parser, _ in subfmts: + match = parser.match(timestr) + if match is None: + continue + groups = match.groupdict() + year = int(groups["year"]) + mon = int(groups["mon"]) + mday = int(groups["mday"]) + frac_digits = groups["fracday"] or "" + fracday = ( + int(frac_digits) / (10 ** len(frac_digits)) if frac_digits else 0.0 + ) + seconds = fracday * 86400.0 + hour = int(seconds // 3600) + minute = int((seconds % 3600) // 60) + sec = seconds - hour * 3600 - minute * 60 + return [year, mon, mday, hour, minute, sec] + raise ValueError(f"Time {timestr} does not match {self.name} format") def str_kwargs(self): - """ - Generator that yields a dict of values corresponding to the - calendar date and time for the internal JD values. - """ - scale = self.scale.upper().encode('ascii'), - iys, ims, ids, ihmsfs = d2dtf(scale, self.precision, - self.jd1, self.jd2) - - # Get the str_fmt element of the first allowed output subformat - _, _, str_fmt = self._select_subfmts(self.out_subfmt)[0] - - yday = None - has_yday = True if '{yday:' in str_fmt else False - - ihrs = ihmsfs['h'] - imins = ihmsfs['m'] - isecs = ihmsfs['s'] - ifracs = ihmsfs['f'] - for iy, im, id, ihr, imin, isec, ifracsec in numpy.nditer( - [iys, ims, ids, ihrs, imins, isecs, ifracs]): - if has_yday: - yday = datetime(iy, im, id).timetuple().tm_yday - - fracday = (((((ifracsec / 10**self.precision + isec) / 60.0 + imin) / 60.0) + ihr) / 24.0) - fracday = f'{fracday:f}'.split('.')[1][0:self.precision] - yield {'year': int(iy), 'mon': int(im), 'day': int(id), - 'hour': int(ihr), 'min': int(imin), 'sec': int(isec), - 'fracsec': int(ifracsec), 'yday': yday, 'fracday': fracday} + """Yield format kwargs, adding MPC ``fracday`` from h/m/s/fracsec.""" + for kwargs in super().str_kwargs(): + total_sec = ( + (kwargs["hour"] * 60 + kwargs["min"]) * 60 + + kwargs["sec"] + + kwargs["fracsec"] / (10 ** self.precision) + ) + if self.precision > 0: + frac_int = int(round(total_sec / 86400.0 * (10 ** self.precision))) + modulus = 10 ** self.precision + if frac_int >= modulus: + frac_int -= modulus + kwargs["fracday"] = f"{frac_int:0{self.precision}d}" + else: + kwargs["fracday"] = "" + yield kwargs From d282287a4588ae22f877596c6e65834499cdecae Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Tue, 25 Aug 2026 14:18:57 -0700 Subject: [PATCH 5/8] fix the time of test --- test/test_query.py | 36 ++++++++++++++++++++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/test/test_query.py b/test/test_query.py index f6d74b9..e3b896d 100644 --- a/test/test_query.py +++ b/test/test_query.py @@ -1,10 +1,42 @@ from unittest import TestCase + +from astropy import units +from astropy.time import Time + from mp_ephem import horizons +# Minimal Horizons CSV response captured for 2002 MS4 on 2020-01-01/02. +# Keeps this test independent of live JPL queries and Time.now(). +HORIZONS_FIXTURE = [ + b"****", + b" A= 41.96777612227916", + b" Date__(UT)__HR:MN:SS, Date_________JDUT, , , R.A._(ICRF), DEC_(ICRF)," + b" dRA*cosD, d(DEC)/dt, RA_3sigma, DEC_3sigma, APmag, S-brt," + b" SMAA_3sig, SMIA_3sig, Theta, Area_3sig, r, rdot," + b" delta, deldot, phi, PAB-LON, PAB-LAT, S-O-T,/r,", + b"***", + b"$$SOE", + b" 2020-Jan-01 00:00:00, 2458849.500000000,*,m, 281.94116, -6.97964," + b" 3.531514, 0.329693, 0.005, 0.003, 20.441, n.a.," + b" 0.00502, 0.00261, -7.980, 0.0000413, 46.55184511243, -0.3665496," + b" 47.4954884093131, 1.0730811, 0.3365, 282.3613, 16.1159, 16.1532,/T,", + b" 2020-Jan-02 00:00:00, 2458850.500000000,*,m, 281.96516, -6.97745," + b" 3.532164, 0.345306, 0.005, 0.003, 20.440, n.a.," + b" 0.00506, 0.00261, -7.947, 0.0000416, 46.55163339978, -0.3665946," + b" 47.4958693813771, 0.5839054, 0.3339, 282.3753, 16.1162, 16.0271,/T,", + b"$$EOE", + b"# obs: 496 (1954-2023)", +] + class TestQuery(TestCase): def test_target(self): - body = horizons.Body('2002 MS4') + start = Time("2020-01-01", scale="utc") + stop = start + 1 * units.day + body = horizons.Body("2002 MS4", start_time=start, stop_time=stop) + body._data = HORIZONS_FIXTURE + body.predict(start) + self.assertIsInstance(body, horizons.Body) - self.assertAlmostEqual(body.mag, 20.3, places=1) + self.assertAlmostEqual(body.mag, 20.441, places=2) From 26b666c343a3d76e2636503c17213843890a404b Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Tue, 25 Aug 2026 14:25:13 -0700 Subject: [PATCH 6/8] Manual publish --- .github/workflows/pypi-publish.yml | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/.github/workflows/pypi-publish.yml b/.github/workflows/pypi-publish.yml index c65d728..debecc4 100644 --- a/.github/workflows/pypi-publish.yml +++ b/.github/workflows/pypi-publish.yml @@ -1,10 +1,8 @@ -# This workflow will upload a Python Package using Twine when a release is created +# Build and upload the package to PyPI. Run manually from the Actions tab. name: Upload Package to PyPI -on: - push: - branches: - - master +on: + workflow_dispatch: jobs: build: From 922b3d721014a5de6db5d1f83947aae199c2d999 Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Tue, 25 Aug 2026 14:27:44 -0700 Subject: [PATCH 7/8] Upgrade PyPI workflow artifact actions to v4. GitHub deprecated upload-artifact v3; keep upload and download on matching v4 versions. Co-authored-by: Cursor --- .github/workflows/pypi-publish.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/pypi-publish.yml b/.github/workflows/pypi-publish.yml index debecc4..e5cb88f 100644 --- a/.github/workflows/pypi-publish.yml +++ b/.github/workflows/pypi-publish.yml @@ -19,7 +19,7 @@ jobs: - name: Build for upload run: python3 -m build --sdist - name: Store the distribution packages - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: python-package-distributions path: dist/ @@ -36,8 +36,8 @@ jobs: id-token: write steps: - name: Download all the dists - uses: actions/download-artifact@v3 - with: + uses: actions/download-artifact@v4 + with: name: python-package-distributions path: dist/ - name: Publish to PyPI From 657c36b0bc2cc5a96dc0f14b31e3efa30642afdc Mon Sep 17 00:00:00 2001 From: JJ Kavelaars Date: Sat, 3 Oct 2026 09:37:10 -0700 Subject: [PATCH 8/8] Evaluate DE405 at TT using the leap-second offset from UTC. The fixed 64 s UT-to-TDB shift matched 1999 only; TT-UTC is now 69.184 s and changes when leap seconds are inserted. Co-authored-by: Cursor --- orbfit/ephem_earth.c | 120 +++++++++++++++++++++++++++++++++++++-- orbfit/orbfit.h | 7 ++- orbfit/orbfit1.c | 3 + src/mp_ephem/bk_orbit.py | 13 +++-- test/data/o3o08.abg | 8 +-- test/test_orbit.py | 4 +- test/test_time.py | 49 ++++++++++++++++ 7 files changed, 185 insertions(+), 19 deletions(-) diff --git a/orbfit/ephem_earth.c b/orbfit/ephem_earth.c index bc2aabc..f5aac8c 100644 --- a/orbfit/ephem_earth.c +++ b/orbfit/ephem_earth.c @@ -7,12 +7,17 @@ static char vcid[] = "$Id: ephem_earth.c,v 1.1 2006/11/22 20:31:50 observe Exp $ *** returning the location of earth geocenter relative to SSBARY. Also *** have eliminated the nutation & libration routines. *** 8/9/99 gmb -*** add 64s difference between UT (usually used in position tables) and TDB -*** (used in this ephemeris). +*** +*** Callers pass UTC Julian Dates (MPC times). DE405 is argumented in +*** Teph, which differs from TT by under 2 ms, so the interpolators +*** evaluate the ephemeris at TT. TT - UTC = 32.184 s + (TAI - UTC). +*** TAI - UTC is the leap-second offset and is not constant: the old +*** fixed 64 s was the 1999 value (actually 64.184 s) and is 5.184 s +*** short for dates since the 2017 January 1 leap second. +*** Local sidereal time is still computed from the UTC date. ***/ #include "orbfit.h" #include -#define TDBOFFSET (64./86400.) #define FNAMESIZE 256 #define BUFFSIZE 512 /* structures used to keep observatory information: */ @@ -275,6 +280,109 @@ int Initialize_Ephemeris() } } +/**==========================================================================**/ +/** utc_jd_to_tt **/ +/** **/ +/** TT = TAI + 32.184 s, and TAI - UTC is the IERS leap-second offset. **/ +/** Before 1972 the offset also drifts (USNO tai-utc.dat). After the **/ +/** last tabulated leap second the latest offset is held. IERS Bulletin **/ +/** C 72 (2026 July 6) keeps TAI - UTC = 37 s through 2026 December; a **/ +/** newly announced leap second has to be added to the table below. **/ +/** Dates before 1960 (when UTC was not defined) use the 1960 Jan 1 **/ +/** expression. **/ +/** **/ +/**==========================================================================**/ + +/* TAI - UTC (seconds) at a UTC Julian Date. The fraction of jd_utc is the + * fraction of that UTC day, including a smeared leap second. */ +static double tai_minus_utc(double jd_utc) +{ + /* jd0: UTC JD when this expression starts. + dat0, mjd0, rate: TAI-UTC = dat0 + (MJD - mjd0) * rate seconds. + Values are the IERS/USNO tai-utc series. 1960-01-01 continues the + 1961 rate back to the start of UTC. */ + static const struct { + double jd0, dat0, mjd0, rate; + } leap[] = { + { 2436934.5, 1.4178180, 37300.0, 0.0012960 }, /* 1960-01-01 */ + { 2437300.5, 1.4228180, 37300.0, 0.0012960 }, /* 1961-01-01 */ + { 2437512.5, 1.3728180, 37300.0, 0.0012960 }, /* 1961-08-01 */ + { 2437665.5, 1.8458580, 37665.0, 0.0011232 }, /* 1962-01-01 */ + { 2438334.5, 1.9458580, 37665.0, 0.0011232 }, /* 1963-11-01 */ + { 2438395.5, 3.2401300, 38761.0, 0.0012960 }, /* 1964-01-01 */ + { 2438486.5, 3.3401300, 38761.0, 0.0012960 }, /* 1964-04-01 */ + { 2438639.5, 3.4401300, 38761.0, 0.0012960 }, /* 1964-09-01 */ + { 2438761.5, 3.5401300, 38761.0, 0.0012960 }, /* 1965-01-01 */ + { 2438820.5, 3.6401300, 38761.0, 0.0012960 }, /* 1965-03-01 */ + { 2438942.5, 3.7401300, 38761.0, 0.0012960 }, /* 1965-07-01 */ + { 2439004.5, 3.8401300, 38761.0, 0.0012960 }, /* 1965-09-01 */ + { 2439126.5, 4.3131700, 39126.0, 0.0025920 }, /* 1966-01-01 */ + { 2439887.5, 4.2131700, 39126.0, 0.0025920 }, /* 1968-02-01 */ + { 2441317.5, 10.0, 0.0, 0.0 }, /* 1972-01-01 */ + { 2441499.5, 11.0, 0.0, 0.0 }, /* 1972-07-01 */ + { 2441683.5, 12.0, 0.0, 0.0 }, /* 1973-01-01 */ + { 2442048.5, 13.0, 0.0, 0.0 }, /* 1974-01-01 */ + { 2442413.5, 14.0, 0.0, 0.0 }, /* 1975-01-01 */ + { 2442778.5, 15.0, 0.0, 0.0 }, /* 1976-01-01 */ + { 2443144.5, 16.0, 0.0, 0.0 }, /* 1977-01-01 */ + { 2443509.5, 17.0, 0.0, 0.0 }, /* 1978-01-01 */ + { 2443874.5, 18.0, 0.0, 0.0 }, /* 1979-01-01 */ + { 2444239.5, 19.0, 0.0, 0.0 }, /* 1980-01-01 */ + { 2444786.5, 20.0, 0.0, 0.0 }, /* 1981-07-01 */ + { 2445151.5, 21.0, 0.0, 0.0 }, /* 1982-07-01 */ + { 2445516.5, 22.0, 0.0, 0.0 }, /* 1983-07-01 */ + { 2446247.5, 23.0, 0.0, 0.0 }, /* 1985-07-01 */ + { 2447161.5, 24.0, 0.0, 0.0 }, /* 1988-01-01 */ + { 2447892.5, 25.0, 0.0, 0.0 }, /* 1990-01-01 */ + { 2448257.5, 26.0, 0.0, 0.0 }, /* 1991-01-01 */ + { 2448804.5, 27.0, 0.0, 0.0 }, /* 1992-07-01 */ + { 2449169.5, 28.0, 0.0, 0.0 }, /* 1993-07-01 */ + { 2449534.5, 29.0, 0.0, 0.0 }, /* 1994-07-01 */ + { 2450083.5, 30.0, 0.0, 0.0 }, /* 1996-01-01 */ + { 2450630.5, 31.0, 0.0, 0.0 }, /* 1997-07-01 */ + { 2451179.5, 32.0, 0.0, 0.0 }, /* 1999-01-01 */ + { 2453736.5, 33.0, 0.0, 0.0 }, /* 2006-01-01 */ + { 2454832.5, 34.0, 0.0, 0.0 }, /* 2009-01-01 */ + { 2456109.5, 35.0, 0.0, 0.0 }, /* 2012-07-01 */ + { 2457204.5, 36.0, 0.0, 0.0 }, /* 2015-07-01 */ + { 2457754.5, 37.0, 0.0, 0.0 } /* 2017-01-01 */ + }; + const int nleap = (int)(sizeof leap / sizeof leap[0]); + int i; + double mjd; + + if (jd_utc < leap[0].jd0) + i = 0; + else { + for (i = nleap - 1; i > 0; i--) + if (jd_utc >= leap[i].jd0) break; + } + + mjd = jd_utc - 2400000.5; + return leap[i].dat0 + (mjd - leap[i].mjd0) * leap[i].rate; +} + +double utc_jd_to_tt(double jd_utc) +{ + /* UTC JD runs one calendar day per JD day, so a leap-second day is + * squeezed into 86400 JD seconds. Undo that, then add TAI-UTC at 0h + * and the 32.184 s TT-TAI offset. Same construction as SOFA utctai. */ + double jd0, fd, dat0, dat12, dat24, dlod, dleap; + const double tt_minus_tai = 32.184; + const double sec_per_day = 86400.0; + + jd0 = floor(jd_utc - 0.5) + 0.5; + fd = jd_utc - jd0; + dat0 = tai_minus_utc(jd0); + dat12 = tai_minus_utc(jd0 + 0.5); + dat24 = tai_minus_utc(jd0 + 1.0); + dlod = 2.0 * (dat12 - dat0); + dleap = dat24 - (dat0 + dlod); + fd *= (sec_per_day + dleap) / sec_per_day; + fd *= (sec_per_day + dlod) / sec_per_day; + return jd0 + fd + (dat0 + tt_minus_tai) / sec_per_day; +} + /**==========================================================================**/ /** Interpolate_Position **/ /** **/ @@ -284,7 +392,7 @@ int Initialize_Ephemeris() /** the function Read_Coefficients (when necessary). **/ /** **/ /** Inputs: **/ -/** Time -- Time for which position is desired (Julian Date). **/ +/** Time -- UTC Julian Date. Converted to TT before interpolation. **/ /** Target -- Solar system body for which position is desired. **/ /** Position -- Pointer to external array to receive the position. **/ /** **/ @@ -298,7 +406,7 @@ void Interpolate_Position( double Time , int Target , double Position[3] ) int i , j; long int C , G , N , offset = 0; - Time += TDBOFFSET; /****!!!!! Account for difference between TDB and UT**/ + Time = utc_jd_to_tt(Time); /*--------------------------------------------------------------------------*/ /* This function doesn't "do" nutations or librations. */ /*--------------------------------------------------------------------------*/ @@ -432,7 +540,7 @@ void Interpolate_State(double Time, int i , j; long int C , G , N , offset = 0; - Time += TDBOFFSET; /****!!!!! Account for difference between TDB and UT**/ + Time = utc_jd_to_tt(Time); /*--------------------------------------------------------------------------*/ /* This function doesn't "do" nutations or librations. */ diff --git a/orbfit/orbfit.h b/orbfit/orbfit.h index 6155d18..1f8d9a6 100644 --- a/orbfit/orbfit.h +++ b/orbfit/orbfit.h @@ -131,7 +131,12 @@ void earth3d(double t, /* time is in years here */ int obscode, double *x, double *y, double *z); -/*ICRS vector from SSBary to observatory location */ +/* UTC Julian Date to TT Julian Date. DE405 is evaluated at TT. */ +double +utc_jd_to_tt(double jd_utc); + +/*ICRS vector from SSBary to observatory location. + * jd is a UTC Julian Date. */ void earth_ssbary(double jd, int obscode, diff --git a/orbfit/orbfit1.c b/orbfit/orbfit1.c index 6f1a915..f811269 100644 --- a/orbfit/orbfit1.c +++ b/orbfit/orbfit1.c @@ -266,6 +266,9 @@ int scan_observation(char *inbuff, OBSERVATION *obs, OBSERVATION *previous) obs->thetax = DTOR * hmsdeg(rastring); obs->thetay = DTOR * dmsdeg(decstring); obs->dthetax = obs->dthetay = obs->dthetay * ARCSEC; + /* No site code on this line. Keep the frame origin at the + * barycenter; mpc3d places the supplied geocentric vector. */ + obs->obscode = OBSCODE_SSBARY; return(-2); } } diff --git a/src/mp_ephem/bk_orbit.py b/src/mp_ephem/bk_orbit.py index 42ff32a..c23f39a 100644 --- a/src/mp_ephem/bk_orbit.py +++ b/src/mp_ephem/bk_orbit.py @@ -152,10 +152,11 @@ def _fit_radec(self, randomize=False): # noqa: C901 dec = coordinate.dec.to_string(unit=units.degree, decimal=False, sep=":", precision=1, alwayssign=True, pad=True) + # UTC JD. The ephemeris converts this to TT; do not pass TT here. if obs.location is None: - _mpc_file.write("{} {} {} {} {}\n".format(obs.date.jd, ra, dec, res, int(obs.observatory_code))) + _mpc_file.write("{} {} {} {} {}\n".format(obs.date.utc.jd, ra, dec, res, int(obs.observatory_code))) else: - _mpc_file.write("{} {} {} {} {} {} {}\n".format(obs.date.jd, ra, dec, res, + _mpc_file.write("{} {} {} {} {} {} {}\n".format(obs.date.utc.jd, ra, dec, res, obs.location.x, obs.location.y, obs.location.z)) @@ -541,7 +542,7 @@ def predict_helio(self, date, obs_code=568, abg_file=None, minimum_delta=None): if -minimum_delta < self.time - _date < minimum_delta: return - jd = ctypes.c_double(_date.jd) + jd = ctypes.c_double(_date.utc.jd) if abg_file is None: abg_file = tempfile.NamedTemporaryFile(suffix='.abg') abg_file.write(bytes(self.abg, 'utf-8')) @@ -586,7 +587,7 @@ def predict(self, date, obs_code=568, abg_file=None, minimum_delta=None): if -minimum_delta < self.time - _date < minimum_delta: return - jd = ctypes.c_double(_date.jd) + jd = ctypes.c_double(_date.utc.jd) if abg_file is None: abg_file = tempfile.NamedTemporaryFile(suffix='.abg') abg_file.write(bytes(self.abg, 'utf-8')) @@ -636,8 +637,8 @@ def rate_of_motion(self, date=None): date = Time(date, scale='utc') try: - start_date = date.jd - end_date = date.jd + 1 + start_date = date.utc.jd + end_date = date.utc.jd + 1 except Exception as e: logging.error(str(e)) return None diff --git a/test/data/o3o08.abg b/test/data/o3o08.abg index 9ade5a7..da5aeb0 100644 --- a/test/data/o3o08.abg +++ b/test/data/o3o08.abg @@ -1,11 +1,11 @@ # Exact a, adot, b, bdot, g, gdot: - 0.00001882 0.03844611 -0.00000100 -0.00319662 0.03237759 0.00873942 + 0.00001882 0.03844611 -0.00000100 -0.00319661 0.03237758 0.00873920 # Covariance matrix: - 2.3263e-13 2.9633e-13 -9.4644e-16 -4.9622e-14 3.8138e-13 1.4037e-11 + 2.3263e-13 2.9633e-13 -9.4643e-16 -4.9622e-14 3.8138e-13 1.4037e-11 2.9633e-13 1.5949e-12 1.9597e-15 -1.9623e-13 1.1534e-12 5.4228e-11 --9.4644e-16 1.9597e-15 9.3316e-14 -9.2562e-14 -1.3215e-15 1.1610e-13 +-9.4643e-16 1.9597e-15 9.3316e-14 -9.2562e-14 -1.3215e-15 1.1610e-13 -4.9622e-14 -1.9623e-13 -9.2562e-14 2.1178e-13 -1.5176e-13 -7.6622e-12 3.8138e-13 1.1534e-12 -1.3215e-15 -1.5176e-13 1.0685e-12 4.2575e-11 1.4037e-11 5.4228e-11 1.1610e-13 -7.6622e-12 4.2575e-11 2.1086e-09 # lat0 lon0 xBary yBary zBary JD0 - 6.8056487 -120.3801450 0.6426374 0.0913120 -0.7649650 2456392.051150 + 6.8056487 -120.3801450 0.6426369 0.0913121 -0.7649654 2456392.051150 diff --git a/test/test_orbit.py b/test/test_orbit.py index 40e9b68..e134863 100644 --- a/test/test_orbit.py +++ b/test/test_orbit.py @@ -36,7 +36,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.8026798, 3) + self.assertAlmostEqual(orbit.a.to('au').value, 43.8015219, 3) class OrbitFit(unittest.TestCase): @@ -68,7 +68,7 @@ def test_orbit(self): self.assertAlmostEqual(self.orbit.inc.to(units.degree).value, 8.05, 2) self.assertAlmostEqual(self.orbit.Node.to(units.degree).value, 113.85, 2) self.assertAlmostEqual(self.orbit.om.to(units.degree).value, 66.24, 2) - self.assertAlmostEqual(self.orbit.T.to(units.day).value, 2447884.5070, 3) + self.assertAlmostEqual(self.orbit.T.to(units.day).value, 2447884.5762, 3) self.assertAlmostEqual(self.orbit.epoch.jd, 2456392.05115, 4) def test_summarize(self): diff --git a/test/test_time.py b/test/test_time.py index d7f5951..46c677b 100644 --- a/test/test_time.py +++ b/test/test_time.py @@ -1,4 +1,8 @@ +import ctypes +import glob +import os import unittest + import mp_ephem from astropy.time import Time @@ -11,3 +15,48 @@ def test_precision(self): this_time = Time(time_str, scale='utc', precision=time_precision) # print(this_time.isot) self.assertEqual(str(this_time), time_str) + + +class UtcToTtTest(unittest.TestCase): + """DE405 is evaluated at TT. The UTC offset must follow leap seconds.""" + + @classmethod + def setUpClass(cls): + lib = glob.glob(os.path.join(os.path.dirname(mp_ephem.__file__), 'orbfit*.so'))[0] + cls.orbfit = ctypes.CDLL(lib) + cls.orbfit.utc_jd_to_tt.restype = ctypes.c_double + cls.orbfit.utc_jd_to_tt.argtypes = [ctypes.c_double] + + def assert_matches_astropy(self, utc_string): + utc = Time(utc_string, scale='utc') + got = self.orbfit.utc_jd_to_tt(utc.jd) + # A tenth of a millisecond is far inside the ephemeris time tolerance. + self.assertLess(abs(got - utc.tt.jd) * 86400.0, 1.0e-4, + msg='{}: C TT {} vs astropy {}'.format(utc_string, got, utc.tt.jd)) + + def test_leap_second_boundaries(self): + for stamp in ( + '1960-01-01', + '1965-06-15', + '1971-12-31', + '1972-01-01', + '1961-07-31 18:00:00', + '1998-12-31 12:00:00', + '1998-12-31 23:59:59', + '1999-01-01', + '2005-12-31 23:59:59', + '2006-01-01', + '2016-12-31 12:00:00', + '2016-12-31 23:59:59', + '2016-12-31 23:59:60', + '2017-01-01', + '2013-04-03 15:06:08', + '2026-10-03', + ): + self.assert_matches_astropy(stamp) + + def test_modern_offset_is_not_the_1999_constant(self): + utc = Time('2026-10-03', scale='utc') + offset_seconds = (self.orbfit.utc_jd_to_tt(utc.jd) - utc.jd) * 86400.0 + self.assertAlmostEqual(offset_seconds, 69.184, places=3) + self.assertGreater(abs(offset_seconds - 64.0), 5.0)