Skip to content

Cursor/utc to tt ephemeris #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ijiraq merged 9 commits into from
Oct 3, 2026
Merged

Cursor/utc to tt ephemeris #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
e9bced1
Add likelihood from CLASSY formated comment lines.
ijiraq Jul 26, 2025
317ce41
Add likelihood from CLASSY formated comment lines.
ijiraq Jul 26, 2025
c9fdd9c
Merge remote-tracking branch 'origin/master'
ijiraq Jul 26, 2025
ddf3e42
Increase version number to track jd modification
ijiraq Apr 14, 2026
1eae694
become compatiable with astropy 6 and higher
ijiraq Aug 25, 2026
d282287
fix the time of test
ijiraq Aug 25, 2026
26b666c
Manual publish
ijiraq Aug 25, 2026
922b3d7
Upgrade PyPI workflow artifact actions to v4.
ijiraq Aug 25, 2026
657c36b
Evaluate DE405 at TT using the leap-second offset from UTC.
ijiraq Oct 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 6 additions & 8 deletions .github/workflows/pypi-publish.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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:
Expand All @@ -21,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/
Expand All @@ -38,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
Expand Down
120 changes: 114 additions & 6 deletions orbfit/ephem_earth.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <string.h>
#define TDBOFFSET (64./86400.)
#define FNAMESIZE 256
#define BUFFSIZE 512
/* structures used to keep observatory information: */
Expand Down Expand Up @@ -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 **/
/** **/
Expand All @@ -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. **/
/** **/
Expand All @@ -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. */
/*--------------------------------------------------------------------------*/
Expand Down Expand Up @@ -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. */
Expand Down
7 changes: 6 additions & 1 deletion orbfit/orbfit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
3 changes: 3 additions & 0 deletions orbfit/orbfit1.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
Expand Down
2 changes: 1 addition & 1 deletion setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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']

Expand Down
2 changes: 1 addition & 1 deletion src/mp_ephem/__version__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1 @@
version = "0.13.0"
version = "0.14.0"
13 changes: 7 additions & 6 deletions src/mp_ephem/bk_orbit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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'))
Expand Down Expand Up @@ -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'))
Expand Down Expand Up @@ -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
Expand Down
10 changes: 8 additions & 2 deletions src/mp_ephem/ephem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]:
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading