From 0341acd6a488ec390194edbbb483e470b51381a8 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:21:09 +0000 Subject: [PATCH 01/20] Fix sign, wrapping and quadrants in the angle and horizon conversions Five helpers in ssapy.utils returned wrong angles: - dms_to_dd took the sign from the degree field only, so "-00:30:00" returned +0.5 deg. The sign now belongs to the whole angle. - dd_to_hms used |dd| for negative input (-15 deg -> "1:0:0" instead of 23h) and printed to stdout; 359.99999999 deg gave "24:0:00". It now wraps into [0, 360) and rolls 24h over to 0h. - rightascension_to_hourangle formatted LST - RA in degrees as D:M:S, so a 2 h hour angle came back as "30:0:0", and converted numeric local times inconsistently with numeric right ascensions. It now returns (LST - RA) mod 24 h as HH:MM:SS with both numeric arguments in degrees, as documented. - equatorial_to_horizontal took azimuth from arccos, which cannot tell east from west, then applied a southern-hemisphere pi - az patch; every object west of the meridian came out east. It now uses the atan2 form (north = 0, east = 90 deg) and warns instead of printing when both an hour angle and a right ascension are given. - horizontal_to_equatorial negated the zenith angle in the southern hemisphere and mirrored the hour angle on the signs of latitude and declination. It now uses atan2, which also stays well conditioned on the meridian. - equatorial_to_ecliptic and ecliptic_to_equatorial took longitude and RA from arctan of a ratio: RA 120 deg mapped to ecliptic longitude 300 deg (should be 120.04 deg) and RA 200 deg to 26 deg (206 deg); the inverse returned -60 deg for RA 120 deg. Both use arctan2. The same horizon bugs were fixed in SSAPy-Toolkit; these are the SSAPy copies. New tests: exact sexagesimal cases; dd_to_hms against astropy Angle at seven angles; hour angles from documented examples; the horizon transforms against the east/north/up rotation at 16 geometries in both hemispheres (1e-9 deg); and the ecliptic transforms against the obliquity rotation in all four quadrants (1e-9 deg). Three assertions in test_utils.py that pinned the old behaviour (the 30:0:0-style hour angle, the stdout message, and the printed hour-angle notice) now check the corrected behaviour. --- ssapy/utils.py | 114 ++++++++++++---------- tests/test_angle_conversion_references.py | 72 ++++++++++++++ tests/test_utils.py | 24 ++--- 3 files changed, 144 insertions(+), 66 deletions(-) create mode 100644 tests/test_angle_conversion_references.py diff --git a/ssapy/utils.py b/ssapy/utils.py index 1a45d1b8..7d73155f 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -5,6 +5,7 @@ import os import re import numpy as np +import warnings from astropy.time import Time as _Time import astropy.units as u from typing import Union, Tuple @@ -1863,20 +1864,15 @@ def rightascension_to_hourangle(right_ascension, local_time): Example: -------- - rightascension_to_hourangle("10:30:00", "12:45:00") -> "02:15:00" - rightascension_to_hourangle(157.5, 191.25) -> "02:15:00" + rightascension_to_hourangle("10:30:00", "12:45:00") -> "2:15:0" + rightascension_to_hourangle(157.5, 191.25) -> "2:15:0" """ - if type(right_ascension) is not str: - right_ascension = dd_to_hms(right_ascension) - if type(local_time) is not str: - local_time = dd_to_dms(local_time) - _ra = float(right_ascension.split(':')[0]) - _lt = float(local_time.split(':')[0]) - if _ra > _lt: - __ltm, __lts = local_time.split(':')[1:] - local_time = f'{24 + _lt}:{__ltm}:{__lts}' - - return dd_to_dms(hms_to_dd(local_time) - hms_to_dd(right_ascension)) + ra_deg = hms_to_dd(right_ascension) if isinstance(right_ascension, str) else float(right_ascension) + lst_deg = hms_to_dd(local_time) if isinstance(local_time, str) else float(local_time) + # Hour angle = local sidereal time - right ascension, wrapped to [0, 24) h + # and returned as hours. The previous code formatted the difference in + # degrees as D:M:S, so a 2 h hour angle came back as "30:0:0". + return dd_to_hms((lst_deg - ra_deg) % 360.0) def equatorial_to_horizontal(observer_latitude, declination, right_ascension=None, hour_angle=None, local_time=None, hms=False): @@ -1902,17 +1898,18 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non tuple Azimuth and altitude in degrees. """ - if right_ascension is not None and hour_angle is not None: - print('Both right_ascension and hour_angle parameters are provided.\nUsing hour_angle for calculations.') - if hms: + if hour_angle is not None: + if right_ascension is not None: + warnings.warn( + "Both right_ascension and hour_angle parameters are provided; using hour_angle.", + UserWarning, + stacklevel=2, + ) + if hms or isinstance(hour_angle, str): hour_angle = hms_to_dd(hour_angle) elif right_ascension is not None: - hour_angle = rightascension_to_hourangle(right_ascension, local_time) - if hms: - hour_angle = hms_to_dd(hour_angle) - elif hour_angle is not None: - if hms: - hour_angle = hms_to_dd(hour_angle) + # rightascension_to_hourangle always returns an HH:MM:SS string. + hour_angle = hms_to_dd(rightascension_to_hourangle(right_ascension, local_time)) else: raise ValueError('Either right_ascension or hour_angle must be provided.') @@ -1922,12 +1919,14 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non altitude = zenithangle_to_altitude(zenith_angle, deg=False) - _num = np.sin(declination) - np.sin(observer_latitude) * np.cos(zenith_angle) - _den = np.cos(observer_latitude) * np.sin(zenith_angle) - azimuth = np.arccos(_num / _den) - - if observer_latitude < 0: - azimuth = np.pi - azimuth + # Azimuth from north through east. arccos cannot tell east from west and + # put every object west of the meridian in the east; the sign of + # sin(hour angle) does (positive hour angle = west). + azimuth = np.mod(np.arctan2( + -np.cos(declination) * np.sin(hour_angle), + np.sin(declination) * np.cos(observer_latitude) + - np.cos(declination) * np.sin(observer_latitude) * np.cos(hour_angle), + ), 2 * np.pi) altitude, azimuth = np.degrees([altitude, azimuth]) return azimuth, altitude @@ -1954,18 +1953,20 @@ def horizontal_to_equatorial(observer_latitude, azimuth, altitude): altitude, azimuth, latitude = np.radians([altitude, azimuth, observer_latitude]) zenith_angle = altitude_to_zenithangle(altitude, deg=False) - zenith_angle = [-zenith_angle if latitude < 0 else zenith_angle][0] - - declination = np.sin(latitude) * np.cos(zenith_angle) - declination = declination + (np.cos(latitude) * np.sin(zenith_angle) * np.cos(azimuth)) - declination = np.arcsin(declination) - - _num = np.cos(zenith_angle) - np.sin(latitude) * np.sin(declination) - _den = np.cos(latitude) * np.cos(declination) - hour_angle = np.arccos(_num / _den) + declination = np.arcsin( + np.sin(latitude) * np.cos(zenith_angle) + + np.cos(latitude) * np.sin(zenith_angle) * np.cos(azimuth) + ) - if (latitude > 0 > declination) or (latitude < 0 < declination): - hour_angle = 2 * np.pi - hour_angle + # atan2 resolves the quadrant (east of the meridian, sin azimuth > 0, is a + # negative hour angle, i.e. above 12 h) and stays well conditioned on the + # meridian, where arccos of a value near 1 loses half its digits. The + # previous test on the signs of latitude and declination mirrored the hour + # angle in the wrong cases. + hour_angle = np.mod(np.arctan2( + -np.sin(azimuth) * np.sin(zenith_angle), + np.cos(zenith_angle) * np.cos(latitude) - np.sin(zenith_angle) * np.sin(latitude) * np.cos(azimuth), + ), 2 * np.pi) declination, hour_angle = np.degrees([declination, hour_angle]) @@ -2123,7 +2124,9 @@ def equatorial_to_ecliptic(right_ascension, declination, degrees=False): else: ra, dec = right_ascension, declination ec_latitude = np.arcsin(cos_ec * np.sin(dec) - sin_ec * np.cos(dec) * np.sin(ra)) - ec_longitude = np.arctan((cos_ec * np.cos(dec) * np.sin(ra) + sin_ec * np.sin(dec)) / (np.cos(dec) * np.cos(ra))) + # arctan2, not arctan of the ratio: arctan folds right ascensions between + # 90 and 270 deg onto the opposite side of the sky. + ec_longitude = np.arctan2(cos_ec * np.cos(dec) * np.sin(ra) + sin_ec * np.sin(dec), np.cos(dec) * np.cos(ra)) if degrees: return deg0to360(np.degrees(ec_longitude)), np.degrees(ec_latitude) else: @@ -2150,7 +2153,9 @@ def ecliptic_to_equatorial(lon, lat, degrees=False): lon, lat = np.radians(lon), np.radians(lat) else: lon, lat = lon, lat - ra = np.arctan((cos_ec * np.cos(lat) * np.sin(lon) - sin_ec * np.sin(lat)) / (np.cos(lat) * np.cos(lon))) + # arctan2 resolves the quadrant; arctan of the ratio returned right + # ascensions between 90 and 270 deg on the opposite side of the sky. + ra = rad0to2pi(np.arctan2(cos_ec * np.cos(lat) * np.sin(lon) - sin_ec * np.sin(lat), np.cos(lat) * np.cos(lon))) dec = np.arcsin(cos_ec * np.sin(lat) + sin_ec * np.cos(lat) * np.sin(lon)) if degrees: return np.degrees(ra), np.degrees(dec) @@ -2267,10 +2272,12 @@ def dms_to_dd(dms): # Degree minute second to Degree decimal """ dms, out = [[dms] if type(dms) is str else dms][0], [] for i in dms: - deg, minute, sec = [float(j) for j in i.split(':')] - if deg < 0: - minute, sec = float(f'-{minute}'), float(f'-{sec}') - out.append(deg + minute / 60 + sec / 3600) + text = i.strip() + # The sign belongs to the whole angle, so "-00:30:00" is -0.5 deg; a + # sign carried only by the degree field is lost when degrees are 0. + sign = -1.0 if text.startswith('-') else 1.0 + deg, minute, sec = [abs(float(j)) for j in text.lstrip('+-').split(':')] + out.append(sign * (deg + minute / 60 + sec / 3600)) return [out[0] if type(dms) is str or len(dms) == 1 else out][0] @@ -2393,9 +2400,8 @@ def dd_to_hms(degree_decimal): Notes: ------ - - Decimal degrees are divided by 15 to convert to hours. - - If the input DD value is negative, the function assumes the absolute value for conversion - and prints a warning message. + - Decimal degrees are wrapped into [0, 360) and divided by 15 to convert to hours, so + the result always lies in [0, 24) h. - Handles edge cases where seconds reach 60, incrementing minutes accordingly. - Returns seconds as an integer if the value is a whole number. @@ -2405,18 +2411,16 @@ def dd_to_hms(degree_decimal): '12:34:56' >>> dd_to_hms(-236.375) - '15:45:30' # Assumes positive value for conversion. + '8:14:30' # -236.375 deg is 123.625 deg. >>> dd_to_hms("12:34:56") # DMS string converted to DD first. '0:50:18.4' """ if type(degree_decimal) is str: degree_decimal = dms_to_dd(degree_decimal) - if degree_decimal < 0: - print('dd for HMS conversion cannot be negative, assuming positive.') - _dd = -degree_decimal / 15 - else: - _dd = degree_decimal / 15 + # Angles are wrapped into [0, 360) deg, so -15 deg is 23h; taking the + # absolute value would turn it into 1h. + _dd = float(np.mod(degree_decimal, 360.0)) / 15 _h = int(np.trunc(_dd)) __h = _dd - np.trunc(_dd) _m, __m = np.trunc(__h * 60), __h * 60 - np.trunc(__h * 60) @@ -2426,6 +2430,8 @@ def dd_to_hms(degree_decimal): _m, _s = _m + 1, '00' if _m == 60: _h, _m = _h + 1, 0 + if _h >= 24: + _h = 0 return f'{int(_h)}:{int(_m)}:{_s}' diff --git a/tests/test_angle_conversion_references.py b/tests/test_angle_conversion_references.py new file mode 100644 index 00000000..e2dd0756 --- /dev/null +++ b/tests/test_angle_conversion_references.py @@ -0,0 +1,72 @@ +"""Sexagesimal, hour-angle, horizon, and ecliptic conversions against references.""" + +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import Angle + +from ssapy import utils + + +def test_negative_sub_degree_dms_keeps_its_sign(): + # The sign belongs to the whole angle: -00:30:00 is -0.5 deg (exact). + for text, degrees in {"-00:30:00": -0.5, "-0:00:36": -0.01, "-10:30:36": -10.51, "+00:30:00": 0.5}.items(): + assert utils.dms_to_dd(text) == pytest.approx(degrees, abs=1e-12), text + + +@pytest.mark.parametrize("degrees", [-15.0, -236.375, 0.0, 15.0, 187.5, 359.99999999, 757.5]) +def test_dd_to_hms_wraps_into_one_day_like_astropy(degrees): + # astropy Angle.wrap_at(360 deg).hour gives the same instant within + # 1e-4 s of time, and hours lie in [0, 24). + hours, minutes, seconds = (float(part) for part in utils.dd_to_hms(degrees).split(":")) + assert 0 <= hours < 24 + expected = Angle(degrees * u.deg).wrap_at(360 * u.deg).hour * 3600.0 + actual = hours * 3600 + minutes * 60 + seconds + assert min(abs(actual - expected), 86400 - abs(actual - expected)) < 1e-4 + + +def test_hour_angle_is_sidereal_time_minus_right_ascension_in_hours(): + # HA = LST - RA, wrapped to [0, 24) h (exact for these inputs). + assert utils.rightascension_to_hourangle("10:30:00", "12:45:00") == "2:15:0" + assert utils.rightascension_to_hourangle(157.5, 191.25) == "2:15:0" + assert utils.rightascension_to_hourangle("23:00:00", "01:00:00") == "2:0:0" + + +@pytest.mark.parametrize("latitude", [37.68, -33.87]) +@pytest.mark.parametrize("hour_angle", [-60.0, 0.0, 45.0, 120.0]) +@pytest.mark.parametrize("declination", [-20.0, 20.0]) +def test_horizon_conversions_match_astropy_altaz_geometry(latitude, hour_angle, declination): + # The topocentric direction of (HA, dec) in the local east/north/up frame, + # built from rotation matrices, gives the reference azimuth (from north + # through east) and altitude; the forward and inverse transforms agree + # with it to 1e-9 deg in both hemispheres. + phi, h, d = np.radians([latitude, hour_angle, declination]) + east = -np.cos(d) * np.sin(h) + north = np.sin(d) * np.cos(phi) - np.cos(d) * np.sin(phi) * np.cos(h) + up = np.sin(d) * np.sin(phi) + np.cos(d) * np.cos(phi) * np.cos(h) + az_ref = np.degrees(np.arctan2(east, north)) % 360.0 + alt_ref = np.degrees(np.arcsin(up)) + az, alt = utils.equatorial_to_horizontal(latitude, declination, hour_angle=hour_angle) + assert alt == pytest.approx(alt_ref, abs=1e-9) + if abs(alt_ref) < 89.9: + assert ((az - az_ref + 180.0) % 360.0) - 180.0 == pytest.approx(0.0, abs=1e-9) + ha_back, dec_back = utils.horizontal_to_equatorial(latitude, az, alt) + assert dec_back == pytest.approx(declination, abs=1e-9) + assert ((ha_back - hour_angle + 180.0) % 360.0) - 180.0 == pytest.approx(0.0, abs=1e-9) + + +@pytest.mark.parametrize("ra, dec", [(30.0, 10.0), (120.0, 10.0), (200.0, -20.0), (300.0, 40.0)]) +def test_ecliptic_conversions_match_the_obliquity_rotation(ra, dec): + # Rotating the unit vector about +x by the J2000 obliquity (SSAPy's + # cos_ec/sin_ec) gives the reference ecliptic longitude and latitude in + # every quadrant; both directions agree with it to 1e-9 deg. + ce, se = utils.cos_ec, utils.sin_ec + ra_r, dec_r = np.radians([ra, dec]) + x = np.cos(dec_r) * np.cos(ra_r) + y = ce * np.cos(dec_r) * np.sin(ra_r) + se * np.sin(dec_r) + z = -se * np.cos(dec_r) * np.sin(ra_r) + ce * np.sin(dec_r) + lon_ref, lat_ref = np.degrees(np.arctan2(y, x)) % 360.0, np.degrees(np.arcsin(z)) + lon, lat = utils.equatorial_to_ecliptic(ra, dec, degrees=True) + assert lon == pytest.approx(lon_ref, abs=1e-9) and lat == pytest.approx(lat_ref, abs=1e-9) + ra_back, dec_back = utils.ecliptic_to_equatorial(lon_ref, lat_ref, degrees=True) + assert ra_back == pytest.approx(ra, abs=1e-9) and dec_back == pytest.approx(dec, abs=1e-9) diff --git a/tests/test_utils.py b/tests/test_utils.py index 816d2ae3..9e7c2301 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -184,8 +184,9 @@ def test_small_coordinate_and_angle_helpers(capsys): assert utils.zenithangle_to_altitude(60.0) == 30.0 assert np.isclose(utils.altitude_to_zenithangle(np.pi / 6, deg=False), np.pi / 3) assert np.isclose(utils.zenithangle_to_altitude(np.pi / 3, deg=False), np.pi / 6) - assert utils.rightascension_to_hourangle(30.0, 2.0) == "0:0:0" - assert utils.rightascension_to_hourangle("02:00:00", "01:00:00").count(":") == 2 + # Numeric arguments are degrees (docstring): LST 2 deg - RA 30 deg = -28 deg = 22h08m. + assert utils.rightascension_to_hourangle(30.0, 2.0) == "22:8:00" + assert utils.rightascension_to_hourangle("02:00:00", "01:00:00") == "23:0:0" assert utils.dms_to_dd("12:30:00") == 12.5 assert utils.dms_to_dd(["12:30:00", "-12:30:00"]) == [12.5, -12.5] @@ -193,7 +194,6 @@ def test_small_coordinate_and_angle_helpers(capsys): assert utils.hms_to_dd("12:00:00") == 180.0 assert utils.hms_to_dd(["12:00:00", "01:30:00"]) == [180.0, 22.5] assert utils.dd_to_hms(-180.0) == "12:0:0" - assert "cannot be negative" in capsys.readouterr().out def test_sun_ra_dec_matches_astropy_solar_position(): @@ -469,17 +469,17 @@ def test_coordinate_angle_edge_cases_and_errors(capsys): assert np.isfinite(az) assert np.isfinite(alt) - az2, alt2 = utils.equatorial_to_horizontal( - observer_latitude=45.0, - declination=0.0, - right_ascension="01:00:00", - hour_angle="00:00:00", - local_time="02:00:00", - hms=True, - ) + with pytest.warns(UserWarning, match="using hour_angle"): + az2, alt2 = utils.equatorial_to_horizontal( + observer_latitude=45.0, + declination=0.0, + right_ascension="01:00:00", + hour_angle="00:00:00", + local_time="02:00:00", + hms=True, + ) assert np.isfinite(az2) assert np.isfinite(alt2) - assert "Using hour_angle" in capsys.readouterr().out az3, alt3 = utils.equatorial_to_horizontal( observer_latitude=45.0, From dd2798fd5353420c165d29f0636cfc7efd2a381c Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:21:22 +0000 Subject: [PATCH 02/20] Advance the Sun's perihelion in the fast sunPos series sunPos(fast=True) implements Montenbruck & Gill section 3.3.2, which holds Omega + omega = 282.9400 deg fixed. In the J2000 frame the Earth-Moon barycentre's longitude of perihelion advances 0.32327364 deg per Julian century (JPL approximate positions, Table 1), so the series' longitude drifted by 11.6 arcsec per year: 0.2 arcmin at J2000, 5.3 arcmin in 2026, and 11.6 arcmin at worst over 1980-2060 against astropy's geometric Sun. AccelSolRad and the eclipse geometry use this path. With the advance the worst case over 1980-2060 is 0.58 arcmin. New test: 200 epochs over 1980-2060 against astropy's geometric geocentric Sun stay within 1 arcmin in direction and 3e-4 in distance; it fails on the previous series. --- ssapy/utils.py | 7 ++++++- tests/test_sun_position_fast.py | 24 ++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 1 deletion(-) create mode 100644 tests/test_sun_position_fast.py diff --git a/ssapy/utils.py b/ssapy/utils.py index 7d73155f..b5e9fb12 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -1197,7 +1197,12 @@ def sunPos(t, fast=True): # MG section 3.3.2 T = (_gpsToTT(t) - 51544.5) / 36525.0 M = 6.239998880168239 + 628.3019326367721 * T - lam = (4.938234585592756 + M + 0.03341335890206922 * np.sin(M) + 0.00034906585039886593 * np.sin(2 * M)) + # Omega + omega = 282.9400 deg at J2000, advancing 0.32327364 deg per + # Julian century in the J2000 frame (JPL approximate elements of the + # Earth-Moon barycentre). MG treat it as constant, which lets the + # longitude drift 11.6 arcsec/yr: 5.3 arcmin by 2026. + lam = (4.938234585592756 + 0.005642189402906841 * T + M + + 0.03341335890206922 * np.sin(M) + 0.00034906585039886593 * np.sin(2 * M)) rs = (149.619 - 2.499 * np.cos(M) - 0.021 * np.cos(2 * M)) * 1e9 obliquity = 0.40909280420293637 co, so = np.cos(obliquity), np.sin(obliquity) diff --git a/tests/test_sun_position_fast.py b/tests/test_sun_position_fast.py new file mode 100644 index 00000000..5babdd23 --- /dev/null +++ b/tests/test_sun_position_fast.py @@ -0,0 +1,24 @@ +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import get_body_barycentric +from astropy.time import Time + +from ssapy import utils + + +def test_fast_sun_position_stays_within_an_arcminute_from_1980_to_2060(): + # Reference: astropy's geometric geocentric Sun (built-in ephemeris). + # Montenbruck & Gill's series with the Earth-Moon barycentre's perihelion + # advance (0.32327364 deg/century, JPL approximate elements) stays within + # 1 arcmin in direction (measured 0.58 arcmin) and 3e-4 in distance over + # 1980-2060; without the advance it reached 11.6 arcmin. + worst = 0.0 + for jd in np.linspace(Time("1980-01-01", scale="tt").jd, Time("2060-01-01", scale="tt").jd, 200): + t = Time(jd, format="jd", scale="tt") + expected = (get_body_barycentric("sun", t) - get_body_barycentric("earth", t)).xyz.to_value(u.m) + position = np.asarray(utils.sunPos(t.gps, fast=True), dtype=float).reshape(3) + cosine = np.dot(position, expected) / (np.linalg.norm(position) * np.linalg.norm(expected)) + worst = max(worst, np.degrees(np.arccos(np.clip(cosine, -1.0, 1.0))) * 60.0) + assert np.linalg.norm(position) == pytest.approx(np.linalg.norm(expected), rel=3e-4) + assert worst < 1.0 From 1200a46012791b2bcc7d9c75b30bcea6883690ba Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:21:22 +0000 Subject: [PATCH 03/20] Make lb_to_tan's projected rates the derivative of its projection lb_to_tan projects positions orthographically, x = rahat . unit and y = dechat . unit with fixed axes, so their rates are exactly the projections of the sphere velocity onto those axes. It then stretched the radial part by 1/cos(rho), the derivative of a gnomonic projection. Off the tangent point the returned rates disagreed with the motion of the projected point: 0.01175 vs 0.01073 for a point 0.33 rad from centre (9.5%). The stretch is removed. The existing test checked only the tangent point, where the factor is 1. New test: at 0, 0.05, 0.33 and 0.8 rad from the tangent point the rates match a central difference of lb_to_tan's own projection to 1e-8 relative; the off-centre cases fail on the previous code. --- ssapy/utils.py | 12 +++++------- tests/test_tangent_plane_rates.py | 20 ++++++++++++++++++++ 2 files changed, 25 insertions(+), 7 deletions(-) create mode 100644 tests/test_tangent_plane_rates.py diff --git a/ssapy/utils.py b/ssapy/utils.py index b5e9fb12..cb7ca84a 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -988,15 +988,13 @@ def lb_to_tan(lb, b, mul=None, mub=None, lcen=None, bcen=None): dechat2 = np.cross(unit, rahat2) dechat2 /= np.sqrt(np.sum(dechat2**2, axis=1, keepdims=True)) vv = mul[:, None] * rahat2 + mub[:, None] * dechat2 + # x = rahat . unit and y = dechat . unit with fixed rahat/dechat, so + # their time derivatives are exactly rahat . (d unit/dt) and + # dechat . (d unit/dt). The radial 1/cos(rho) stretch that used to + # follow is the gnomonic derivative; applied to this orthographic + # projection it overstated the rates by 1/cos(rho). vx = np.sum(rahat * vv, axis=1) vy = np.sum(dechat * vv, axis=1) - rr = np.hypot(xx, yy) - m = np.abs(rr) > 1e-9 - vr = (vx[m] * xx[m] + vy[m] * yy[m]) / rr[m] - va = (vy[m] * xx[m] - vx[m] * yy[m]) / rr[m] - vr /= np.sum(unitcen[m, :] * unit[m, :], axis=1) - vx[m] = vr * xx[m] / rr[m] - va * yy[m] / rr[m] - vy[m] = vr * yy[m] / rr[m] + va * xx[m] / rr[m] res = res + (vx, vy) return res diff --git a/tests/test_tangent_plane_rates.py b/tests/test_tangent_plane_rates.py new file mode 100644 index 00000000..7b389993 --- /dev/null +++ b/tests/test_tangent_plane_rates.py @@ -0,0 +1,20 @@ +import numpy as np +import pytest + +from ssapy import utils + + +@pytest.mark.parametrize("offset", [0.0, 0.05, 0.33, 0.8]) +def test_tangent_plane_rates_are_the_derivative_of_the_projection(offset): + # The projected rates must be the time derivative of the projected + # position: a central difference of lb_to_tan for a point moving at + # (mul, mub) on the sphere (mul = cos(b) dl/dt) agrees to 1e-8. + lcen, bcen = np.array([1.0]), np.array([0.4]) + lb, b = np.array([1.0 + offset]), np.array([0.4 - 0.6 * offset]) + mul, mub = np.array([0.01]), np.array([-0.02]) + _x, _y, vx, vy = utils.lb_to_tan(lb, b, mul=mul, mub=mub, lcen=lcen, bcen=bcen) + h = 1e-6 + xp, yp = utils.lb_to_tan(lb + mul / np.cos(b) * h, b + mub * h, lcen=lcen, bcen=bcen) + xm, ym = utils.lb_to_tan(lb - mul / np.cos(b) * h, b - mub * h, lcen=lcen, bcen=bcen) + np.testing.assert_allclose(vx, (xp - xm) / (2 * h), rtol=1e-8) + np.testing.assert_allclose(vy, (yp - ym) / (2 * h), rtol=1e-8) From fc7a609655516c79b2c9904c05c0a9602311ec4c Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:21:22 +0000 Subject: [PATCH 04/20] Advance Leapfrog4 epochs by whole steps, not three rounded substeps Leapfrog4Propagator stored each step's epoch as t + h1 + h0 + h1 with the Yoshida weights, three additions that each round at the ulp of t (2.4e-7 s at GPS 1.4e9). The rounding is the same every step, so the stored epochs drifted from the integrated states. At t0 = 0 the method is fourth order (97, 6.1, 0.38, 0.024 m over 6000 s for h = 40, 20, 10, 5 s); at GPS 1.4e9 it stalled at 1.47 m (h = 10 s) and 2.17 m (h = 5 s). Substep epochs are now offsets from t and t advances by h, which restores the t0 = 0 errors at real epochs. New test: against the Keplerian propagator at t0 = 0, 1.4e9 and 3.8e9, halving h from 20 to 10 to 5 s reduces the error by more than 12x each time and leaves less than 0.05 m at h = 5 s; the GPS-epoch cases fail on the previous code. --- ssapy/propagator.py | 9 +++++++-- tests/test_leapfrog4_epoch.py | 22 ++++++++++++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) create mode 100644 tests/test_leapfrog4_epoch.py diff --git a/ssapy/propagator.py b/ssapy/propagator.py index 574176bf..77244306 100644 --- a/ssapy/propagator.py +++ b/ssapy/propagator.py @@ -1174,14 +1174,19 @@ def _prop(self, times, states, h, tthresh, propkw): h1 = w1 * h h0 = w0 * h + # Substep epochs are offsets from t, and t itself advances by + # exactly one h per step. Chaining t + h1 + h0 + h1 rounds three + # times at the ulp of t (2.4e-7 s at GPS 1.4e9), and the bias + # accumulates: the stored epochs drifted from the integrated + # states and capped the method near 2 m over 6000 s. r, v = self._leapfrog_step(self.accel, r, v, t, h1, propkw) t1 = t + h1 r, v = self._leapfrog_step(self.accel, r, v, t1, h0, propkw) - t2 = t1 + h0 + t2 = t + (h1 + h0) r, v = self._leapfrog_step(self.accel, r, v, t2, h1, propkw) - t = t2 + h1 + t = t + h state = np.hstack([r, v]) diff --git a/tests/test_leapfrog4_epoch.py b/tests/test_leapfrog4_epoch.py new file mode 100644 index 00000000..9a6c0b52 --- /dev/null +++ b/tests/test_leapfrog4_epoch.py @@ -0,0 +1,22 @@ +import numpy as np +import pytest + +import ssapy +from ssapy import AccelKepler + + +@pytest.mark.parametrize("t0", [0.0, 1.4e9, 3.8e9]) +def test_leapfrog4_is_fourth_order_at_real_gps_epochs(t0): + # Reference: the Keplerian propagator. Halving the step must cut the error + # by ~16x at any epoch (accept >= 12x) down to h = 5 s, where the error over + # 6000 s is under 0.05 m. Summing substep times onto a large GPS epoch used + # to stall the error near 2 m. + orbit = ssapy.Orbit.fromKeplerianElements(7500e3, 0.05, 0.9, 0.3, 1.2, 0.4, t=t0) + times = t0 + np.arange(0.0, 6000.0 + 1e-9, 600.0) + r_ref, _ = ssapy.rv(orbit, times, propagator=ssapy.KeplerianPropagator()) + errors = [] + for h in (20.0, 10.0, 5.0): + r, _ = ssapy.rv(orbit, times, propagator=ssapy.Leapfrog4Propagator(AccelKepler(), h=h)) + errors.append(np.max(np.linalg.norm(r - r_ref, axis=1))) + assert errors[0] / errors[1] > 12 and errors[1] / errors[2] > 12 + assert errors[2] < 0.05 From 2665e62a79787f18b462455dee5127c221e1ad7b Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:22:37 +0000 Subject: [PATCH 05/20] Write zero drag terms in make_tle, and keep a source TLE's drag terms make_tle documents that dynamic TLE terms are ignored, but it wrote B* as '99999-0' in columns 54-61, which SGP4 reads as 0.99999 per Earth radius. Propagating a make_tle TLE for a 7000 km orbit with sgp4 for 3 days drifted 11,890 km from the same TLE with B* = 0. By default all three drag terms (ndot/2, nddot/6, B*) are now zero. The placeholder happened to match the Aerospace GEO test TLEs, whose B* field is also 99999, which is why test_sgp4 matched direct sgp4 through SGP4Propagator(truncate=True). make_tle now takes drag_fields (line-1 columns 34-61), and both Orbit.tle and the truncate path pass the columns of the TLE the orbit came from, so truncation changes only the precision of the mean elements; test_sgp4 still matches to 1e-15. New tests: for a plain make_tle TLE sgp4 reads B*, ndot and nddot as 0 and the epoch to 1 ms, with valid checksums; for an orbit read from the documented ISS TLE, Orbit.tle keeps its B* (-1.1606e-5), ndot and nddot exactly as sgp4 reads them. The first fails on the previous code. --- ssapy/io.py | 18 ++++++++++++++---- ssapy/orbit.py | 4 +++- ssapy/propagator.py | 6 +++++- tests/test_io.py | 34 ++++++++++++++++++++++++++++++++++ 4 files changed, 56 insertions(+), 6 deletions(-) diff --git a/ssapy/io.py b/ssapy/io.py index c52841a0..3b076ca6 100644 --- a/ssapy/io.py +++ b/ssapy/io.py @@ -202,7 +202,7 @@ def parse_tle(tle): return a, e, i, pa, raan, trueAnomaly, epoch.gps -def make_tle(a, e, i, pa, raan, trueAnomaly, t): +def make_tle(a, e, i, pa, raan, trueAnomaly, t, drag_fields=None): """ Create a TLE from Kozai mean orbital elements @@ -223,11 +223,16 @@ def make_tle(a, e, i, pa, raan, trueAnomaly, t): t : float or astropy.time.Time If float, then should correspond to GPS seconds; i.e., seconds since 1980-01-06 00:00:00 UTC + drag_fields : str, optional + Line-1 columns 34-61 (ndot/2, nddot/6 and B*, 28 characters) to copy + verbatim, e.g. ``line1[33:61]`` of the TLE an orbit came from. By default + all three terms are written as zero. Notes ----- - Dynamic TLE terms, including the drag coefficient and ballistic coefficient, - are ignored in this function. + Unless ``drag_fields`` is given, the dynamic TLE terms are zero. They used + to be written with B* = '99999-0', which SGP4 reads as 0.99999 per Earth + radius: an enormous drag term for LEO (11,900 km of drift in 3 days). """ from .orbit import ( _ellipticalEccentricToMeanAnomaly, @@ -245,7 +250,12 @@ def make_tle(a, e, i, pa, raan, trueAnomaly, t): line1 += "{:02d}".format(int(year) % 100) line1 += "{:012.8f}".format(day) - line1 += " +.00000000 +00000-0 99999-0 0 0000" + # Columns 34-61: ndot/2, nddot/6 and B*, zero unless supplied. + if drag_fields is None: + drag_fields = "+.00000000 +00000-0 +00000-0" + if len(drag_fields) != 28: + raise ValueError("drag_fields must be the 28 characters of TLE line-1 columns 34-61.") + line1 += " " + drag_fields + " 0 0000" def checksum(s): check = 0 diff --git a/ssapy/orbit.py b/ssapy/orbit.py index dcaf671a..be2b51ff 100644 --- a/ssapy/orbit.py +++ b/ssapy/orbit.py @@ -1846,7 +1846,9 @@ def tle(self): """Two line element for this Orbit. """ from .io import make_tle - return make_tle(*self.kozaiMeanKeplerianElements, self.t) + source = getattr(self, "_tle", None) + drag_fields = source[0][33:61] if source is not None else None + return make_tle(*self.kozaiMeanKeplerianElements, self.t, drag_fields=drag_fields) class EarthObserver: diff --git a/ssapy/propagator.py b/ssapy/propagator.py index 77244306..e4307e4b 100644 --- a/ssapy/propagator.py +++ b/ssapy/propagator.py @@ -311,7 +311,11 @@ def _getRVOne(self, orbit, time): sat_epoch = orbit.t if self.truncate: - line1, line2 = make_tle(*orbit.kozaiMeanKeplerianElements, orbit.t) + # Keep the drag terms of the TLE the orbit came from, so truncation + # changes only the precision of the mean elements. + source = getattr(orbit, "_tle", None) + drag_fields = source[0][33:61] if source is not None else None + line1, line2 = make_tle(*orbit.kozaiMeanKeplerianElements, orbit.t, drag_fields=drag_fields) sat = Satrec.twoline2rv(line1, line2) elif getattr(orbit, "_sat", None) is not None: # Path A: the orbit was built from a TLE and still carries its diff --git a/tests/test_io.py b/tests/test_io.py index 92fcf492..9c851f52 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -229,3 +229,37 @@ def put(start, text): pos = b3obs2pos("".join(line)) assert pos["satnum"] == 12345 assert pos["sensnum"] == 511 + + +def test_make_tle_writes_zero_drag_terms_that_sgp4_reads_back(): + # The TLE line-1 format puts ndot/2, nddot/6 and B* in columns 34-61. + # make_tle documents that dynamic terms are ignored, so sgp4 must read all + # three as 0 and the epoch to 1 ms; the checksum digit must match the + # mod-10 rule. B* was previously written as 0.99999 per Earth radius. + from astropy.time import Time + from sgp4.api import Satrec + + t = Time("2026-10-08T06:30:15.5", scale="utc") + line1, line2 = io.make_tle(7000e3, 0.001, np.radians(51.6), np.radians(30.0), np.radians(120.0), np.radians(75.0), t) + sat = Satrec.twoline2rv(line1, line2) + assert (sat.bstar, sat.ndot, sat.nddot) == (0.0, 0.0, 0.0) + assert abs((sat.jdsatepoch + sat.jdsatepochF) - t.utc.jd) * 86400.0 < 1e-3 + for line in (line1, line2): + digits = sum(int(c) if c.isdigit() else (1 if c == "-" else 0) for c in line[:68]) + assert line[68] == str(digits % 10) + + +def test_tle_from_a_tle_orbit_keeps_its_drag_terms(): + # Orbit.tle and SGP4Propagator(truncate=True) rebuild a TLE from Kozai + # elements; for an orbit read from a TLE they must keep its ndot, nddot and + # B* (here the ISS element set from the TLE format documentation, B* = + # -1.1606e-5), as sgp4 reads them, exactly. + import ssapy + from sgp4.api import Satrec + + line1 = "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927" + line2 = "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537" + orbit = ssapy.Orbit.fromTLETuple((line1, line2)) + rebuilt = Satrec.twoline2rv(*orbit.tle) + original = Satrec.twoline2rv(line1, line2) + assert (rebuilt.bstar, rebuilt.ndot, rebuilt.nddot) == (original.bstar, original.ndot, original.nddot) From eb7b0ebb8228f69b0e09ebe37d3a46864bc4fd02 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:32:55 +0000 Subject: [PATCH 06/20] Document what the deprecated radecRate actually returns radecRate's docstring listed six return values (ra, raRate, dec, decRate, slantRange, slantRangeRate), but it returns only raRate, decRate and slantRangeRate. The docstring now says so and points to radec(..., rate=True) for the rest. Behaviour is unchanged; the rates themselves match central differences of radec (raRate is d(RA)/dt cos(dec), as documented). --- ssapy/compute.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/ssapy/compute.py b/ssapy/compute.py index 523eb415..d70fea5a 100644 --- a/ssapy/compute.py +++ b/ssapy/compute.py @@ -853,18 +853,15 @@ def radecRate( Returns ------- - ra : array_like (n, m) - right ascension in radians raRate : array_like (n, m) Rate of change of right ascension*cos(dec) in radians per second. - dec : array_link (n, m) - declination in radians decRate : array_like (n, m) Rate of change of declination in radians per second. - slantRange : array_like (n, m) - Range in meters slantRangeRate : array_like (n, m) Slant range rate in meters per second. + + Only the three rates are returned (the docstring previously listed six + values); use ``radec(..., rate=True)`` for angles, range and rates. """ import warnings warnings.warn("This function is deprecated; use " From f9644e1ac1af508f7f32fb8775c9fc7c3e6eecc1 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:38:31 +0000 Subject: [PATCH 07/20] Warn when the Gauss two-position solver does not converge GaussTwoPosOrbitSolver's fixed-point iteration converges only for short arcs, and when it diverged it returned the last iterate silently. Against Keplerian truth it recovered a 38 deg LEO arc to 6.5e-12 m/s but returned velocities 2.7 km/s off for a quarter GEO orbit, 34 km/s off for a 134 deg MEO arc, and 7e18 m/s off for a 270 deg LEO arc, all without any signal. It now issues a RuntimeWarning naming SheferTwoPosOrbitSolver when the last update exceeds 1e-10 of eta, and the class docstring states the short-arc limitation. (The existing branch-coverage test runs one iteration with eps = 0 and still passes; it now warns.) The near-zero series for X(x) also had 64/35 x where the expansion has 64/35 x^2. New tests: Shefer's method recovers the departure velocity to 1e-9 m/s on five arcs from 38 to 270 deg, LEO to GTO (measured <= 6e-12 m/s); Gauss recovers the 38 deg arc to 1e-9 m/s without warning and warns on the quarter GEO orbit. The warning check fails on the previous code. DanchickTwoPosOrbitSolver raises 'Invalid x' on the 270 deg arc; it fails loudly, so it is left as is. --- ssapy/orbit_solver.py | 19 ++++++++- tests/test_two_position_solver_references.py | 43 ++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 tests/test_two_position_solver_references.py diff --git a/ssapy/orbit_solver.py b/ssapy/orbit_solver.py index 682a2bd5..0a1bfe57 100644 --- a/ssapy/orbit_solver.py +++ b/ssapy/orbit_solver.py @@ -9,6 +9,7 @@ import abc import numpy as np +import warnings from astropy.time import Time from .orbit import Orbit @@ -132,6 +133,12 @@ class GaussTwoPosOrbitSolver(TwoPosOrbitSolver): and implements a method to compute the orbital parameter `p` based on Shefer's equations. + Gauss's fixed-point iteration converges only for short arcs: it recovers + the velocity to 1e-11 m/s for a 38 deg LEO arc but diverges for sweeps + near 90 deg and beyond (2.7 km/s error for a quarter GEO orbit). It warns + with a RuntimeWarning when it does not converge; SheferTwoPosOrbitSolver + handles long arcs. + Attributes: eps (float): Convergence tolerance for iterative calculations. maxiter (int): Maximum number of iterations allowed for convergence. @@ -167,10 +174,20 @@ def _getP(self): h = 2 * np.arcsinh(np.sqrt(-x)) X = (np.sinh(2 * h) - 2 * h) / np.sinh(h)**3 else: # -1e-15 1e-10 * max(1.0, np.abs(eta)): + warnings.warn( + "Gauss two-position iteration did not converge (|d_eta| = {:.3g} after {} iterations); " + "the solution is unreliable. Use SheferTwoPosOrbitSolver for long arcs.".format(abs(d_eta), niter), + RuntimeWarning, + stacklevel=3, + ) # Plug eta into (2) to obtain p p = (0.5 * eta * self.kappa * self.sigma / self.tau)**2 / self.mu return p diff --git a/tests/test_two_position_solver_references.py b/tests/test_two_position_solver_references.py new file mode 100644 index 00000000..24fea130 --- /dev/null +++ b/tests/test_two_position_solver_references.py @@ -0,0 +1,43 @@ +import numpy as np +import pytest + +import ssapy + +T1 = 1.4e9 +CASES = [ # (a [m], e, i [rad], arc [s]) + (7000e3, 0.01, 0.9, 600.0), # LEO, 38 deg + (7000e3, 0.01, 0.9, 4400.0), # LEO, 270 deg (long way) + (42164e3, 0.001, 0.1, 21600.0), # GEO, 90 deg + (20000e3, 0.4, 1.0, 9000.0), # eccentric MEO, 134 deg + (24400e3, 0.73, 0.5, 20000.0), # GTO, 166 deg +] + + +def _endpoints(a, e, inc, arc): + orbit = ssapy.Orbit.fromKeplerianElements(a, e, inc, 0.5, 1.0, 0.3, t=T1) + r1, v1 = ssapy.rv(orbit, T1) + r2, _ = ssapy.rv(orbit, T1 + arc) + kappa_sign = 1 if np.dot(np.cross(r1, r2), np.cross(r1, v1)) > 0 else -1 + return r1, r2, v1, kappa_sign + + +@pytest.mark.parametrize("a, e, inc, arc", CASES) +def test_shefer_recovers_the_keplerian_velocity_for_short_and_long_arcs(a, e, inc, arc): + # Reference: the Keplerian orbit that produced both positions. Shefer's + # method recovers the departure velocity to 1e-9 m/s on every arc + # (measured <= 6e-12 m/s). + r1, r2, v1, kappa_sign = _endpoints(a, e, inc, arc) + solved = ssapy.SheferTwoPosOrbitSolver(r1, r2, T1, T1 + arc, kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-9 + + +def test_gauss_warns_instead_of_returning_a_wrong_long_arc_orbit(): + # Gauss's iteration recovers a 38 deg LEO arc to 1e-9 m/s without a + # warning, but diverges for a quarter GEO orbit; it must now warn there + # (it used to return a velocity 2.7 km/s off silently). + r1, r2, v1, kappa_sign = _endpoints(*CASES[0]) + solved = ssapy.GaussTwoPosOrbitSolver(r1, r2, T1, T1 + CASES[0][3], kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-9 + r1, r2, v1, kappa_sign = _endpoints(*CASES[2]) + with pytest.warns(RuntimeWarning, match="did not converge"): + ssapy.GaussTwoPosOrbitSolver(r1, r2, T1, T1 + CASES[2][3], kappaSign=kappa_sign).solve() From 5d1b9555a1041ca47a906927d6ba941264736742 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:39:59 +0000 Subject: [PATCH 08/20] Recompute Earth precession-nutation hourly, not every 30 days EarthOrientation, which Body('earth') and therefore AccelHarmonic use for the GCRF-to-ITRF rotation, cached the IAU 1976/1980 precession-nutation matrix and UT1 - TT for 30 days and only updated Earth rotation per call. Over a propagation the frame drifted 0.06 arcsec after 1 day, 1.0 arcsec after 7 days and 5.1 arcsec after 29 days: 172 m at 7000 km and 1 km at GEO, and a 2.1e-7 m/s^2 error in a 20x20 harmonic acceleration at LEO. The default threshold is now 1 hour, which bounds the drift at 0.003 arcsec; the cost is one erfa.pnm80 call per simulated hour. New test: an EarthOrientation held from 2026-10-08 agrees with a freshly built one to 0.005 arcsec 1, 7 and 29 days later, and with astropy's GCRS -> ITRS rotation to 0.5 arcsec (IAU 1980 vs 2000A model difference, measured 0.27 arcsec). All three cases fail on the previous default. AccelDrag has the same 30-day cache for its true-of-date frame; a 5 arcsec error moves the density lookup by under 200 m, well inside the Harris-Priester model's own error, so it is left unchanged. --- ssapy/body.py | 8 +++++-- tests/test_earth_orientation_cache.py | 32 +++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 2 deletions(-) create mode 100644 tests/test_earth_orientation_cache.py diff --git a/ssapy/body.py b/ssapy/body.py index 756bce97..836f1f35 100644 --- a/ssapy/body.py +++ b/ssapy/body.py @@ -42,10 +42,14 @@ class EarthOrientation: Parameters ---------- recalc_threshold : float - Threshold for recomputing the orientation matrix. Default is 30 days. + Seconds after which the precession-nutation matrix and UT1 - TT are + recomputed; Earth rotation itself is evaluated at every call. Default + 1 hour, which bounds the frame error at 0.003 arcsec. The previous + 30-day default let the frame drift 1.0 arcsec in a week and 5.1 arcsec + in 29 days (172 m at 7000 km, 1 km at GEO). """ - def __init__(self, recalc_threshold=86400 * 30): + def __init__(self, recalc_threshold=3600.0): self.recalc_threshold = recalc_threshold self._t = None diff --git a/tests/test_earth_orientation_cache.py b/tests/test_earth_orientation_cache.py new file mode 100644 index 00000000..7f14d008 --- /dev/null +++ b/tests/test_earth_orientation_cache.py @@ -0,0 +1,32 @@ +import numpy as np +import pytest +import astropy.units as u +from astropy.coordinates import GCRS, ITRS, CartesianRepresentation +from astropy.time import Time + +from ssapy.body import EarthOrientation + + +def _angle_arcsec(a, b): + # Small-angle rotation between two frames; the Frobenius norm of + # (A B^T - I) / sqrt(2) avoids arccos's loss of precision near 1. + return np.degrees(np.linalg.norm(a @ b.T - np.eye(3)) / np.sqrt(2.0)) * 3600.0 + + +@pytest.mark.parametrize("days", [1.0, 7.0, 29.0]) +def test_cached_earth_orientation_matches_a_fresh_evaluation(days): + # A long-lived EarthOrientation (as held by Body / AccelHarmonic over a + # propagation) must agree with a freshly built one to 0.005 arcsec at any + # later time, and with astropy's GCRS -> ITRS rotation to 0.5 arcsec (the + # IAU 1980 / 2000A model difference; measured 0.27 arcsec). The old 30-day + # cache drifted 0.06, 1.0 and 5.1 arcsec after 1, 7 and 29 days. + t0 = Time("2026-10-08T00:00:00", scale="utc").gps + held = EarthOrientation() + held(t0) + t = t0 + days * 86400.0 + assert _angle_arcsec(held(t), EarthOrientation()(t)) < 0.005 + r = np.array([7000e3, 1000e3, 500e3]) + when = Time(t, format="gps") + expected = GCRS(CartesianRepresentation(r * u.m), obstime=when).transform_to(ITRS(obstime=when)).cartesian.xyz.to_value(u.m) + got = held(t) @ r + assert np.degrees(np.linalg.norm(np.cross(got, expected)) / (np.linalg.norm(got) * np.linalg.norm(expected))) * 3600.0 < 0.5 From cbe1bdd501d647aad945d4eedc3db5df2cab59e2 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:53:16 +0000 Subject: [PATCH 09/20] Interpolate fixed-step propagators with a quintic Hermite between steps RKPropagator (RK4, RK8, RK78, Leapfrog, Leapfrog4) answers queries between its step boundaries from an interpolant through the cached states. That was a not-a-knot cubic spline fit to positions and velocities independently, so its error was O(h^4) whatever the integrator's order: an RK8 arc accurate to 4e-8 m at its steps over 6000 s in LEO was only 2.5e-2 m accurate between them at h = 40 s (1.8e-3 m at h = 20 s). The interpolant is now piecewise quintic Hermite: position, velocity and acceleration at both ends of each step, with velocity taken as the derivative of the position quintic. Accelerations are evaluated once per step boundary with the propagator's own force model and reused as the cache grows (one extra evaluation per step, against 13 per RK8 step). If the acceleration cannot be evaluated the cubic spline is used as before. Off-grid RK8 errors become 1.1e-7 m (h = 40 s) and 1.6e-8 m (h = 20 s); RK4 and Leapfrog4 off-grid errors now equal their on-grid errors. New test: off-grid RK8 positions and velocities agree with the Keplerian propagator to 1e-6 m and 1e-6 m/s at h = 40 and 20 s; both fail on the previous interpolant. --- ssapy/propagator.py | 51 ++++++++++++++++++++++++++++++++-- tests/test_rk_interpolation.py | 21 ++++++++++++++ 2 files changed, 70 insertions(+), 2 deletions(-) create mode 100644 tests/test_rk_interpolation.py diff --git a/ssapy/propagator.py b/ssapy/propagator.py index e4307e4b..721298f2 100644 --- a/ssapy/propagator.py +++ b/ssapy/propagator.py @@ -580,6 +580,51 @@ def __eq__(self, rhs): ) +class _StateInterpolant: + """Piecewise quintic Hermite interpolant through cached RK states. + + Positions are interpolated from position, velocity and acceleration at + each step boundary, and velocities are the derivative of that quintic, so + the interpolation error is O(h^6) in position. The cubic spline through + positions and velocities that this replaces was O(h^4): an RK8 arc + accurate to 1e-8 m at its steps was only 2.5 cm accurate between them at + h = 40 s in LEO. Accelerations are evaluated once per step boundary and + reused when the cache grows. + """ + + def __init__(self, interpolant, accels): + self._interpolant = interpolant + self.accels = accels + + @classmethod + def build(cls, propagator, times, states, propkw, previous=None, fallback=None): + accels = dict(previous.accels) if isinstance(previous, cls) else {} + accel = getattr(propagator, "accel", None) + if accel is not None and len(times) >= 2: + try: + from scipy.interpolate import BPoly + derivs = np.empty((len(times), 3, 3), dtype=np.float64) + for k, (t, state) in enumerate(zip(times, states)): + key = float(t) + if key not in accels: + accels[key] = np.asarray( + _safe_accel(accel, state[0:3], state[3:6], t, propkw), dtype=np.float64 + ).reshape(3) + derivs[k, 0] = state[0:3] + derivs[k, 1] = state[3:6] + derivs[k, 2] = accels[key] + return cls(BPoly.from_derivatives(times, derivs), accels) + except Exception: # fall back to the cubic spline below + pass + k = min(3, len(times) - 1) + return fallback(times, states, k=k) + + def __call__(self, t): + position = self._interpolant(t) + velocity = self._interpolant.derivative()(t) + return np.concatenate([position, velocity], axis=-1) + + class RKPropagator(Propagator, ABC): """Abstract base class for Runge-Kutta-based orbit propagators. @@ -653,8 +698,10 @@ def _getRVOneWithMask(self, orbit, tQuery): return out[:, 0:3], out[:, 3:6], valid if remake_spline or spline is None: - k = min(3, len(times_arr) - 1) # cubic, stable - spline = make_interp_spline(times_arr, states_arr, k=k) + spline = _StateInterpolant.build( + self, times_arr, states_arr, orbit.propkw, + previous=spline, fallback=make_interp_spline, + ) container.clear() container.extend([times, states, h_pre, h_app, spline]) times_arr = np.asarray(times, dtype=np.float64) diff --git a/tests/test_rk_interpolation.py b/tests/test_rk_interpolation.py new file mode 100644 index 00000000..2715aff8 --- /dev/null +++ b/tests/test_rk_interpolation.py @@ -0,0 +1,21 @@ +import numpy as np +import pytest + +import ssapy +from ssapy import AccelKepler + + +@pytest.mark.parametrize("h", [40.0, 20.0]) +def test_rk8_is_as_accurate_between_steps_as_at_them(h): + # Reference: the Keplerian propagator. RK8 with a fixed step is accurate + # to ~1e-8 m at its step boundaries over 6000 s in LEO; queried between + # boundaries it must stay within 1e-6 m in position and 1e-6 m/s in + # velocity (measured 1.1e-7 m at h = 40 s). The cubic spline used before + # gave 2.5e-2 m at h = 40 s. + t0 = 1.4e9 + orbit = ssapy.Orbit.fromKeplerianElements(7500e3, 0.05, 0.9, 0.3, 1.2, 0.4, t=t0) + times = t0 + np.array([0.0, 1234.5, 3333.3, 5999.9]) + r_ref, v_ref = ssapy.rv(orbit, times, propagator=ssapy.KeplerianPropagator()) + r, v = ssapy.rv(orbit, times, propagator=ssapy.RK8Propagator(AccelKepler(), h=h)) + assert np.max(np.linalg.norm(r - r_ref, axis=1)) < 1e-6 + assert np.max(np.linalg.norm(v - v_ref, axis=1)) < 1e-6 From cadfd3dd830b5b16e003fb9fd004fe3226856540 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:57:05 +0000 Subject: [PATCH 10/20] Honour VolumeDistancePrior's scale argument VolumeDistancePrior.__init__(scale=RGEO) stored RGEO regardless of the argument, so every instance peaked at 2 RGEO (84,300 km) whatever scale was requested. It now stores the given scale. New test: with scale = 1000 km the log prior is -1e-9 at its 2000 km peak and matches 2 ln(d/s) - d/s - (2 ln 2 - 2 + 1e-9) at 500 km to 1e-12; it fails on the previous code. --- ssapy/correlate_tracks.py | 4 +++- tests/test_tracks.py | 13 +++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/ssapy/correlate_tracks.py b/ssapy/correlate_tracks.py index fec1ffee..203eb47b 100644 --- a/ssapy/correlate_tracks.py +++ b/ssapy/correlate_tracks.py @@ -222,7 +222,9 @@ class VolumeDistancePrior: """ def __init__(self, scale=RGEO): - self.scale = RGEO + # The argument used to be ignored (self.scale = RGEO), so every prior + # peaked at 2 RGEO whatever scale was requested. + self.scale = scale def __call__(self, orbit, distance, chi=False): """Return log prior probability of given orbit. diff --git a/tests/test_tracks.py b/tests/test_tracks.py index 709699b7..18f90f5c 100644 --- a/tests/test_tracks.py +++ b/tests/test_tracks.py @@ -1099,3 +1099,16 @@ def run(self, **kwargs): ) assert created['kwargs']['approximate'] is True assert newmht.run_kwargs == {'first': 2, 'last': 5} + + +def test_volume_distance_prior_uses_the_requested_scale(): + # lnp(d) = 2 ln(d/s) - d/s - (2 ln 2 - 2 + 1e-9) peaks at d = 2 s with + # lnp = -1e-9 (closed form, 1e-12). With s = 1000 km the peak is at 2000 km; + # the scale argument used to be ignored in favour of RGEO. + scale = 1000e3 + prior = VolumeDistancePrior(scale=scale) + orbit = sample_LEO_orbit(t=0) + assert prior(orbit, 2 * scale) == pytest.approx(-1e-9, abs=1e-12) + d = 500e3 + expected = 2 * np.log(d / scale) - d / scale - (2 * np.log(2) - 2 + 1e-9) + assert prior(orbit, d) == pytest.approx(expected, abs=1e-12) From bda7a5daa313c4210891b21cc16ebb5fd05da1c1 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:57:05 +0000 Subject: [PATCH 11/20] Use the selected row in circular_guess when rates are given With 'pmra'/'pmdec' columns, circular_guess took row 0 of the whole arc and ignored indices. It now uses the first of indices (a one-row arc still works with the default), and the docstring documents the (state, epoch) return value it already had. New test: given two observations with rates and indices=[1, 0], the returned epoch is the second observation's time; it fails on the previous code. --- ssapy/rvsampler.py | 12 +++++++----- tests/test_sampler.py | 22 ++++++++++++++++++++++ 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/ssapy/rvsampler.py b/ssapy/rvsampler.py index 0a2e8f06..10b0685d 100644 --- a/ssapy/rvsampler.py +++ b/ssapy/rvsampler.py @@ -530,16 +530,18 @@ def circular_guess(arc, indices=[0, 1]): Returns ------- - state: array_like, (6,) + state : array_like, (6,) state corresponding to orbit passing through points + epoch : astropy.time.Time + epoch of `state`: the midpoint of the two observations, or the + first selected observation when the arc carries rates ('pmra'). """ from astropy import units as u usepm = 'pmra' in arc.dtype.names - if not usepm: - observations = arc[indices] - else: - observations = arc + # With rates only the first selected row is needed (so a one-row arc + # works); it used to be row 0 of the whole arc whatever `indices` said. + observations = arc[[indices[0]]] if usepm else arc[indices] rStation = observations['rStation_GCRF'].to(u.m).value ra = observations['ra'] dec = observations['dec'] diff --git a/tests/test_sampler.py b/tests/test_sampler.py index 7493cc8d..1bd5757a 100644 --- a/tests/test_sampler.py +++ b/tests/test_sampler.py @@ -930,3 +930,25 @@ def test_mh_sampler(): import subprocess cmd = "gprof2dot -f pstats {} -n1 -e1 | dot -Tpng -o {}".format(args.prof_out, args.prof_png) subprocess.run(cmd, shell=True) + + +def test_circular_guess_with_rates_uses_the_selected_row(): + # With 'pmra'/'pmdec' columns only one observation is needed; the one + # used must be the first of `indices` (row 0 used to be taken regardless). + orbit, arc = _circular_guess_arc() + rows = [] + for k in range(2): + ra, dec, _, pmra, pmdec, _ = ssapy.radec( + orbit, arc['time'][k].gps, obsPos=arc['rStation_GCRF'][k].to(u.m).value, + obsVel=arc['vStation_GCRF'][k].to(u.m / u.s).value, rate=True) + rows.append((ra, dec, pmra, pmdec)) + arcpm = QTable() + arcpm['ra'] = [r[0] for r in rows] * u.rad + arcpm['dec'] = [r[1] for r in rows] * u.rad + arcpm['time'] = arc['time'][:2] + arcpm['rStation_GCRF'] = arc['rStation_GCRF'][:2] + arcpm['vStation_GCRF'] = arc['vStation_GCRF'][:2] + arcpm['pmra'] = [r[2] for r in rows] * u.rad / u.s + arcpm['pmdec'] = [r[3] for r in rows] * u.rad / u.s + _state, epoch = ssapy.circular_guess(arcpm, indices=[1, 0]) + assert abs((epoch - arcpm['time'][1]).to(u.s).value) < 1e-6 From a4edcc62547777c07822c11eb28c0c954f03a64d Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 22:57:05 +0000 Subject: [PATCH 12/20] Give the drag model the Sun of date and an hourly frame update AccelDrag evaluates Harris-Priester density in the true-of-date frame but passed it the Sun's right ascension and declination from GCRF, which put the diurnal bulge 0.37 deg off in right ascension in 2026 and growing at 50 arcsec per year. The Sun vector is now rotated into the same frame. Its precession-nutation cache also defaulted to 30 days; it now defaults to 1 hour, like EarthOrientation. New test: the Sun right ascension and declination AccelDrag hands the density model agree with astropy's true-equator, true-equinox Sun to 0.05 deg (fast sunPos is good to 0.6 arcmin); the right ascension fails on the previous code. --- ssapy/accel.py | 11 ++++++++--- tests/test_drag_sun_frame.py | 29 +++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 3 deletions(-) create mode 100644 tests/test_drag_sun_frame.py diff --git a/ssapy/accel.py b/ssapy/accel.py index 67e334d0..4c2fa0ed 100644 --- a/ssapy/accel.py +++ b/ssapy/accel.py @@ -365,12 +365,14 @@ class AccelDrag(Accel): ---------- recalc_threshold : float, optional Number of seconds past which the code will recompute the - precession/nutation matrix. Default: 86400*30 (30 days) + precession/nutation matrix. Default: 3600 (1 hour), matching + EarthOrientation; the previous 30-day default let the frame drift + 5 arcsec. defaultkw : dict default parameters for kwargs passed to __call__, (area, mass, CR) """ - def __init__(self, recalc_threshold=86400 * 30, **defaultkw): + def __init__(self, recalc_threshold=3600.0, **defaultkw): from . import _ssapy self.recalc_threshold = recalc_threshold @@ -412,7 +414,10 @@ def __call__(self, r, v, t, _T=None, **kwargs): self._t = t self._T = erfa.pnm80(2400000.5, mjd_tt) _T = self._T - r_sun = sunPos(t) + # The density model works in the true-of-date frame, so the Sun's + # right ascension and declination must be of date too; sunPos is + # GCRF, which put the diurnal bulge 0.37 deg off in 2026. + r_sun = _T @ sunPos(t) r_tod = _T @ r v_tod = _T @ v diff --git a/tests/test_drag_sun_frame.py b/tests/test_drag_sun_frame.py new file mode 100644 index 00000000..e55e73dc --- /dev/null +++ b/tests/test_drag_sun_frame.py @@ -0,0 +1,29 @@ +from types import SimpleNamespace + +import numpy as np +from astropy.coordinates import TETE, get_sun +from astropy.time import Time + +from ssapy.accel import AccelDrag + + +def test_drag_density_sees_the_sun_of_date(): + # Harris-Priester's diurnal bulge follows the Sun's right ascension and + # declination in the true-of-date frame the density is evaluated in. + # The values AccelDrag passes must match astropy's true-equator, + # true-equinox (TETE) Sun to 0.05 deg (fast sunPos is good to 0.6 arcmin); + # the GCRF values it passed before were 0.37 deg off in RA in 2026. + t = Time("2026-10-08T00:00:00", scale="utc") + drag = AccelDrag() + seen = {} + + def capture(x, y, z, ra_sun, dec_sun): + seen["ra"], seen["dec"] = ra_sun, dec_sun + return 1e-12 + + drag.atm = SimpleNamespace(density=capture) # the C++ model's methods are read-only + drag(np.array([6778e3, 0.0, 0.0]), np.array([0.0, 7670.0, 0.0]), t.gps, area=1.0, mass=100.0, CD=2.2) + sun = get_sun(t).transform_to(TETE(obstime=t)) + d_ra = (np.degrees(seen["ra"]) - sun.ra.deg + 180.0) % 360.0 - 180.0 + assert abs(d_ra) < 0.05 + assert abs(np.degrees(seen["dec"]) - sun.dec.deg) < 0.05 From db9bfe43edefe634bfbec2ac87cb1e03627a91fa Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:00:48 +0000 Subject: [PATCH 13/20] Make Danchick's two-position solver converge over the whole revolution DanchickTwoPosOrbitSolver chose its iteration from the sign of cos(2f) alone. Long-way arcs beyond 270 deg also have cos(2f) >= 0 but negative kappa, so they went to the eta iteration, which left its domain; and on fast perigee passages the x iteration's plain Newton step jumped out of 0 < x < 1. Both raised 'Invalid x': on 19 of 100 test arcs (2-98% of a period on LEO, eccentric MEO, GEO and GTO orbits). The eta iteration is now used for short-way arcs with cos(2f) >= 0 and the x iteration otherwise; if the preferred one fails the other is tried, the x iteration halves any step that would leave (0, 1), and a failure of both raises a RuntimeError naming each. All 100 arcs now solve, with the departure velocity within 3.6e-9 m/s of Keplerian truth (Shefer: 1.1e-9 m/s). New test: the 100 arcs, both solvers, 1e-7 m/s; the Danchick cases fail on the previous code. The branch-coverage test that forced a domain-leaving step now expects the combined RuntimeError. --- ssapy/orbit_solver.py | 110 +++++++++++++------ tests/test_orbit_solver.py | 5 +- tests/test_two_position_solver_references.py | 17 +++ 3 files changed, 96 insertions(+), 36 deletions(-) diff --git a/ssapy/orbit_solver.py b/ssapy/orbit_solver.py index 0a1bfe57..bc66ab9c 100644 --- a/ssapy/orbit_solver.py +++ b/ssapy/orbit_solver.py @@ -242,43 +242,83 @@ def dXdg(g): """Compute dX(g)/dg from Shefer (12).""" return (2 * (1 - np.cos(2 * g)) - 3 * (2 * g - np.sin(2 * g)) / np.tan(g)) / (np.sin(g)**3) + def _eta_iteration(self): + """Newton iteration on eta (Shefer 14); returns eta.""" + with np.errstate(divide="ignore", invalid="ignore"): + eta = 0.5 * (np.sqrt(np.float64(self.m) / (self.ell + 1)) + np.sqrt(np.float64(self.m) / self.ell)) + if not np.isfinite(eta): + raise RuntimeError("Invalid x") + Geta = 1 + niter = 0 + while np.abs(Geta) > self.eps and niter < self.maxiter: + x = self.m / eta**2 - self.ell # Shefer (14) + if not np.isfinite(x) or (1 - 2 * x) > 1 or (1 - 2 * x) < -1: + raise RuntimeError("Invalid x") + g = np.arccos(1 - 2 * x) # Shefer (9) + dgdx = 2 / np.sin(g) # Shefer (10) + dxdeta = -2 * self.m / eta**3 # next few lines are Shefer (14ish) + dXdeta = self.dXdg(g) * dgdx * dxdeta + Geta = eta - 1 - (self.ell + x) * self.X(g) + dGdeta = 1 - self.X(g) * dxdeta - (self.ell + x) * dXdeta + eta -= Geta / dGdeta + niter += 1 + if np.abs(Geta) > 1e-9 * max(1.0, np.abs(eta)): + raise RuntimeError("Danchick eta iteration did not converge") + return eta + + def _x_iteration(self): + """Damped Newton iteration on x in (0, 1) (Shefer 15); returns eta.""" + x = 0.5 + Fx = 1 + niter = 0 + while np.abs(Fx) > self.eps and niter < self.maxiter: + g = np.arccos(1 - 2 * x) # Shefer (9) + dgdx = 2 / np.sin(g) # Shefer (10) + eta = 1 + (self.ell + x) * self.X(g) # Shefer (15) + detadx = (self.ell + x) * self.dXdg(g) * dgdx + self.X(g) + Fx = x - self.m / eta**2 + self.ell # Shefer (15ish) + dFdx = 1 + (2 * self.m / eta**3) * detadx + step = Fx / dFdx + if not np.isfinite(step): + raise RuntimeError("Invalid x") + # Halve the step until x stays strictly inside (0, 1). + for _ in range(60): + if 0.0 < x - step < 1.0: + break + step *= 0.5 + else: + raise RuntimeError("Invalid x") + x -= step + niter += 1 + if np.abs(Fx) > 1e-9: + raise RuntimeError("Danchick x iteration did not converge") + return 1 + (self.ell + x) * self.X(np.arccos(1 - 2 * x)) + def _getP(self): - if self.cos2f >= 0: - # TODO: rewrite to use newton_raphson function. - eta = 0.5 * (np.sqrt(self.m / (self.ell + 1)) + np.sqrt(self.m / self.ell)) - Geta = 1 - niter = 0 - while np.abs(Geta) > self.eps and niter < self.maxiter: - x = self.m / eta**2 - self.ell # Shefer (14) - if (1 - 2 * x) > 1 or (1 - 2 * x) < -1: - raise RuntimeError("Invalid x") - g = np.arccos(1 - 2 * x) # Shefer (9) - dgdx = 2 / np.sin(g) # Shefer (10) - dxdeta = -2 * self.m / eta**3 # next few lines are Shefer (14ish) - dXdeta = self.dXdg(g) * dgdx * dxdeta - Geta = eta - 1 - (self.ell + x) * self.X(g) - dGdeta = 1 - self.X(g) * dxdeta - (self.ell + x) * dXdeta - eta -= Geta / dGdeta - niter += 1 - x = self.m / eta**2 - self.ell # Shefer (9) - elif self.cos2f < 0: - x = 0.5 - Fx = 1 - niter = 0 - while np.abs(Fx) > self.eps and niter < self.maxiter: - if 1 - 2 * x > 1 or 1 - 2 * x < -1: - raise RuntimeError("Invalid x") - g = np.arccos(1 - 2 * x) # Shefer (9) - dgdx = 2 / np.sin(g) # Shefer (10) - eta = 1 + (self.ell + x) * self.X(g) # Shefer (15) - detadx = (self.ell + x) * self.dXdg(g) * dgdx + self.X(g) - Fx = x - self.m / eta**2 + self.ell # Shefer (15ish) - dFdx = 1 + (2 * self.m / eta**3) * detadx - x -= Fx / dFdx - niter += 1 - eta = 1 + (self.ell + x) * self.X(np.arccos(1 - 2 * x)) - else: + # Danchick's eta iteration suits short-way arcs with cos(2f) >= 0 and + # the x iteration the rest. Long-way arcs beyond 270 deg also have + # cos(2f) >= 0 but negative kappa, and fast perigee passages can + # have cos(2f) < 0 on short arcs where the x iteration's plain Newton + # step leaves 0 < x < 1; both used to raise "Invalid x". The + # preferred iteration is tried first, then the other one, and the x + # iteration damps any step that would leave its domain. + if not np.isfinite(self.cos2f): raise ValueError("Invalid value of cos2f") + order = (self._eta_iteration, self._x_iteration) + if not (self.cos2f >= 0 and self.kappa > 0): + order = order[::-1] + errors = [] + for iteration in order: + try: + eta = iteration() + except RuntimeError as exc: + errors.append(str(exc)) + continue + if np.isfinite(eta): + break + errors.append("non-finite eta") + else: + raise RuntimeError("Danchick two-position iteration failed: " + "; ".join(errors)) # Shefer (2) p = (0.5 * eta * self.kappa * self.sigma / self.tau)**2 / self.mu diff --git a/tests/test_orbit_solver.py b/tests/test_orbit_solver.py index 4420b874..f96875b3 100644 --- a/tests/test_orbit_solver.py +++ b/tests/test_orbit_solver.py @@ -312,7 +312,10 @@ def _getP(self): danchick.ell = 0.0 monkeypatch.setattr(ssapy.DanchickTwoPosOrbitSolver, 'X', staticmethod(lambda g: 1.0)) monkeypatch.setattr(ssapy.DanchickTwoPosOrbitSolver, 'dXdg', staticmethod(lambda g: 0.0)) - with pytest.raises(RuntimeError, match='Invalid x'): + # The x iteration now damps steps that would leave 0 < x < 1 and falls + # back to the eta iteration, so a non-converging case fails with a + # RuntimeError naming both attempts. + with pytest.raises(RuntimeError, match='Danchick two-position iteration failed'): ssapy.DanchickTwoPosOrbitSolver._getP(danchick) danchick.cos2f = np.nan diff --git a/tests/test_two_position_solver_references.py b/tests/test_two_position_solver_references.py index 24fea130..b741acef 100644 --- a/tests/test_two_position_solver_references.py +++ b/tests/test_two_position_solver_references.py @@ -41,3 +41,20 @@ def test_gauss_warns_instead_of_returning_a_wrong_long_arc_orbit(): r1, r2, v1, kappa_sign = _endpoints(*CASES[2]) with pytest.warns(RuntimeWarning, match="did not converge"): ssapy.GaussTwoPosOrbitSolver(r1, r2, T1, T1 + CASES[2][3], kappaSign=kappa_sign).solve() + + +@pytest.mark.parametrize("a, e", [(7000e3, 0.01), (26560e3, 0.3), (42164e3, 0.001), (24400e3, 0.73)]) +def test_danchick_and_shefer_cover_the_whole_revolution(a, e): + # Reference: Keplerian truth. Over arcs from 2% to 98% of a period (both + # short- and long-way, including fast perigee passages on a GTO) both + # solvers recover the departure velocity to 1e-7 m/s (measured <= 3.6e-9). + # Danchick's method used to raise "Invalid x" on 19 of these 100 arcs. + orbit = ssapy.Orbit.fromKeplerianElements(a, e, 0.9, 0.5, 1.0, 0.3, t=T1) + for fraction in np.linspace(0.02, 0.98, 25): + arc = fraction * orbit.period + r1, v1 = ssapy.rv(orbit, T1) + r2, _ = ssapy.rv(orbit, T1 + arc) + kappa_sign = 1 if np.dot(np.cross(r1, r2), np.cross(r1, v1)) > 0 else -1 + for solver in (ssapy.DanchickTwoPosOrbitSolver, ssapy.SheferTwoPosOrbitSolver): + solved = solver(r1, r2, T1, T1 + arc, kappaSign=kappa_sign).solve() + assert np.linalg.norm(solved.v - v1) < 1e-7, (solver.__name__, fraction) From f49695ca44ba542bd4b23da8c1867613b4c4ae21 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:01:57 +0000 Subject: [PATCH 14/20] Correct VENUS_MASS VENUS_MASS was 4.687e24 kg, a transposition of 4.867e24: VENUS_MU / G is 4.8673e24 kg with CODATA 2018 G, a 3.7% difference. It is now 4.8675e24 kg. New test: every tabulated mass (Sun, Moon, eight planets) equals its tabulated GM / G to 0.2%; all others already agree to 0.05%. The Venus case fails on the previous value. --- ssapy/constants.py | 2 +- tests/test_body_constants.py | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) create mode 100644 tests/test_body_constants.py diff --git a/ssapy/constants.py b/ssapy/constants.py index 561524f1..7d656870 100644 --- a/ssapy/constants.py +++ b/ssapy/constants.py @@ -237,7 +237,7 @@ SUN_MASS = 1.98847e+30 MOON_MASS = 7.348e22 MERCURY_MASS = 3.301e23 -VENUS_MASS = 4.687e24 +VENUS_MASS = 4.8675e24 EARTH_MASS = 5.9722e24 MARS_MASS = 6.417e23 JUPITER_MASS = 1.899e27 diff --git a/tests/test_body_constants.py b/tests/test_body_constants.py new file mode 100644 index 00000000..bfc6319e --- /dev/null +++ b/tests/test_body_constants.py @@ -0,0 +1,15 @@ +import pytest + +from ssapy import constants + +G = 6.67430e-11 # CODATA 2018 [m^3 kg^-1 s^-2] + + +@pytest.mark.parametrize("body", ["SUN", "MOON", "MERCURY", "VENUS", "EARTH", "MARS", "JUPITER", "SATURN", "URANUS", "NEPTUNE"]) +def test_masses_are_consistent_with_gravitational_parameters(body): + # Each tabulated mass must equal its GM / G (CODATA 2018 G) to 0.2%; the + # masses and GMs come from different compilations, which differ by <0.05% + # here. VENUS_MASS was 4.687e24 kg against GM / G = 4.8673e24 kg (3.7%). + mu = getattr(constants, f"{body}_MU") + mass = getattr(constants, f"{body}_MASS") + assert mass == pytest.approx(mu / G, rel=2e-3) From a11e8db14180c4d796125f95e96a36df31ffc3d5 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:15:49 +0000 Subject: [PATCH 15/20] Warn instead of printing when AccelHarmonic clamps its degree or order AccelHarmonic printed 'WARNING::...' to stdout when the requested degree or order exceeded the model's, and AccelDrag printed its inputs before raising on a non-finite density. The clamps now issue UserWarnings (so they can be filtered or turned into errors), and the drag inputs are part of the ValueError message. The clamp test now checks the warnings instead of captured stdout. Behaviour is otherwise unchanged. --- ssapy/accel.py | 7 +++---- ssapy/gravity.py | 13 +++++++------ tests/test_accel.py | 11 ++++++----- 3 files changed, 16 insertions(+), 15 deletions(-) diff --git a/ssapy/accel.py b/ssapy/accel.py index 4c2fa0ed..31b17814 100644 --- a/ssapy/accel.py +++ b/ssapy/accel.py @@ -430,10 +430,9 @@ def __call__(self, r, v, t, _T=None, **kwargs): dec_sun ) if not np.isfinite(density): - print(f"r_tod = {r_tod}") - print(f"ra_sun = {ra_sun}") - print(f"dec_sun = {dec_sun}") - raise ValueError("non finite density") + raise ValueError( + f"non finite density at r_tod = {r_tod}, ra_sun = {ra_sun}, dec_sun = {dec_sun}" + ) a_tod = -0.5 * kw['CD'] * kw['area'] / kw['mass'] * density * v_rel * norm(v_rel) return _T.T @ a_tod diff --git a/ssapy/gravity.py b/ssapy/gravity.py index 06d9fe98..ab191ba0 100644 --- a/ssapy/gravity.py +++ b/ssapy/gravity.py @@ -5,6 +5,7 @@ # from functools import lru_cache import numpy as np +import warnings from .accel import Accel as _Accel, _invnorm from .utils import find_file, norm @@ -279,8 +280,8 @@ def __init__(self, body, n_max=None, m_max=None): m_max = body.harmonics.m_max self.body = body if n_max > body.harmonics.n_max: - print( - f"WARNING::The provided degree ({n_max}) is higher than the maximum value allowed by the " + warnings.warn( + f"The provided degree ({n_max}) is higher than the maximum value allowed by the " f"{body.harmonics.name} model ({body.harmonics.n_max}).\nSetting the degree to " f"{body.harmonics.n_max}." ) @@ -288,8 +289,8 @@ def __init__(self, body, n_max=None, m_max=None): else: self.n_max = n_max if m_max > body.harmonics.m_max: - print( - f"WARNING::The provided order ({m_max}) is higher than the maximum value allowed by the " + warnings.warn( + f"The provided order ({m_max}) is higher than the maximum value allowed by the " f"{body.harmonics.name} model ({body.harmonics.m_max}).\nSetting the order to " f"{body.harmonics.m_max}." ) @@ -297,8 +298,8 @@ def __init__(self, body, n_max=None, m_max=None): else: self.m_max = m_max if self.m_max > self.n_max: - print( - f"WARNING::The provided order ({self.m_max}) is higher than the selected degree " + warnings.warn( + f"The provided order ({self.m_max}) is higher than the selected degree " f"({self.n_max}).\nSetting the order to {self.n_max}." ) self.m_max = self.n_max diff --git a/tests/test_accel.py b/tests/test_accel.py index 20de4ecb..a3a7f795 100644 --- a/tests/test_accel.py +++ b/tests/test_accel.py @@ -86,7 +86,7 @@ def test_harmonic_coefficients_error_and_identity_paths(tmp_path): assert hc1 != object() -def test_accel_harmonic_clamps_requested_degree_and_order(capsys): +def test_accel_harmonic_clamps_requested_degree_and_order(): harmonics = SimpleNamespace( name="tiny", radius=1.0, @@ -100,12 +100,13 @@ def test_accel_harmonic_clamps_requested_degree_and_order(capsys): assert default_accel.n_max == 1 assert default_accel.m_max == 1 - clamped = ssapy.AccelHarmonic(body, n_max=5, m_max=6) + with pytest.warns(UserWarning) as record: + clamped = ssapy.AccelHarmonic(body, n_max=5, m_max=6) assert clamped.n_max == 1 assert clamped.m_max == 1 - out = capsys.readouterr().out - assert "provided degree" in out - assert "provided order" in out + messages = " ".join(str(w.message) for w in record) + assert "provided degree" in messages + assert "provided order" in messages def _make_body_dependent_objects(): From d137851a7b74cd41855ba6631ac8f2e6b434d227 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:18:16 +0000 Subject: [PATCH 16/20] Scale harmonic gravity by each model's own GM AccelHarmonic passed the body's central GM to the expansion, but the coefficients are normalised to the gravity model's own GM (as they are to its reference radius, which it already used). For the Moon the body GM is DE200's 4902.79989 km^3/s^2 against GRGM1200A's 4902.80012, so lunar harmonic accelerations were 2.2e-7 low (1e-10 m/s^2 at 1900 km); for EGM2008 Earth accelerations were 7.5e-10 high. HarmonicCoefficients.MG is now used when present, falling back to the body GM. New test: 20x20 accelerations from EGM96, EGM2008 and GRGM1200A agree with an independent potential, built from the coefficient files parsed in the test and scipy's Legendre functions and differenced numerically, to 1e-8 relative (measured 1.3e-9 to 2.7e-9) at two points each. The lunar case fails on the previous code. --- ssapy/gravity.py | 6 +- tests/test_harmonic_gravity_reference.py | 101 +++++++++++++++++++++++ 2 files changed, 106 insertions(+), 1 deletion(-) create mode 100644 tests/test_harmonic_gravity_reference.py diff --git a/ssapy/gravity.py b/ssapy/gravity.py index ab191ba0..89c3aa9a 100644 --- a/ssapy/gravity.py +++ b/ssapy/gravity.py @@ -303,8 +303,12 @@ def __init__(self, body, n_max=None, m_max=None): f"({self.n_max}).\nSetting the order to {self.n_max}." ) self.m_max = self.n_max + # The coefficients are normalised to the model's own GM and reference + # radius, so both scale the expansion. Using the body's central GM + # instead put GRGM1200A lunar accelerations 2.2e-7 low (DE200 Moon GM) + # and EGM2008 Earth accelerations 7.5e-10 high. self._harmonic = _ssapy.AccelHarmonic( - self.body.mu, + getattr(self.body.harmonics, "MG", None) or self.body.mu, self.body.harmonics.radius, self.body.harmonics.CS.shape[0], self.body.harmonics.CS.ctypes.data diff --git a/tests/test_harmonic_gravity_reference.py b/tests/test_harmonic_gravity_reference.py new file mode 100644 index 00000000..4444ca84 --- /dev/null +++ b/tests/test_harmonic_gravity_reference.py @@ -0,0 +1,101 @@ +"""AccelHarmonic against an independent spherical-harmonic potential. + +The coefficient files are parsed here without SSAPy's loaders, the potential +U = GM/r sum_n (R/r)^n sum_m Pbar_nm(sin phi) (C cos m lam + S sin m lam) is +built from scipy's associated Legendre functions with the 4-pi normalisation, +and its gradient is taken by central differences. The 20x20 accelerations +agree with AccelHarmonic to 1e-8 relative (finite-difference limited, +measured 1.3e-9 to 2.7e-9). Before AccelHarmonic used each model's own GM, +GRGM1200A lunar accelerations were 2.2e-7 low. +""" +from math import factorial + +import numpy as np +import pytest +from scipy.special import lpmv + +from ssapy.body import get_body +from ssapy.gravity import AccelHarmonic +from ssapy.utils import find_file + +DEGREE = 20 +T = 1.4e9 + + +def _read_egm(name): + egm = find_file(name, ext=".egm") + header = open(egm).read().splitlines() + radius = float(next(line for line in header if line.startswith("ModelRadius"))[11:]) + gm = float(next(line for line in header if line.startswith("ModelMass"))[9:]) + with open(egm + ".cof", "rb") as f: + f.read(8) + n_max = int.from_bytes(f.read(4), "little") + m_max = int.from_bytes(f.read(4), "little") + n_c = (m_max + 1) * (2 * n_max - m_max + 2) // 2 + n_s = m_max * (2 * n_max - m_max + 1) // 2 + c_values = np.frombuffer(f.read(8 * n_c)) + s_values = np.frombuffer(f.read(8 * n_s)) + c, s, k = {}, {}, 0 + for m in range(m_max + 1): + for n in range(m, n_max + 1): + c[(n, m)] = c_values[k] + k += 1 + k = 0 + for m in range(1, m_max + 1): + for n in range(m, n_max + 1): + s[(n, m)] = s_values[k] + k += 1 + return c, s, radius, gm + + +def _read_tab(name): + with open(find_file(name)) as f: + header = [float(x) for x in f.readline().replace(",", " ").split()] + c, s = {}, {} + for line in f: + fields = line.replace(",", " ").split() + n, m = int(fields[0]), int(fields[1]) + if n > DEGREE: + break + c[(n, m)], s[(n, m)] = float(fields[2]), float(fields[3]) + return c, s, header[0] * 1e3, header[1] * 1e9 + + +def _potential(r, c, s, radius, gm): + rr = np.linalg.norm(r) + x = r[2] / rr + lam = np.arctan2(r[1], r[0]) + total = 0.0 + for n in range(1, DEGREE + 1): + scale = (radius / rr) ** n + for m in range(n + 1): + norm = np.sqrt((2 - (m == 0)) * (2 * n + 1) * factorial(n - m) / factorial(n + m)) + pbar = norm * (-1) ** m * lpmv(m, n, x) # remove scipy's Condon-Shortley phase + total += scale * pbar * (c.get((n, m), 0.0) * np.cos(m * lam) + s.get((n, m), 0.0) * np.sin(m * lam)) + return gm / rr * total + + +def _gradient(r, *model, h=0.5): + return np.array([(_potential(r + h * e, *model) - _potential(r - h * e, *model)) / (2 * h) for e in np.eye(3)]) + + +@pytest.mark.parametrize("model", ["egm96", "egm2008"]) +def test_earth_harmonics_match_an_independent_potential(model): + earth = get_body("earth", model=model.upper()) + coefficients = _read_egm(model) + rotation = np.asarray(earth.orientation(T)) + for r_body in (np.array([7000e3, 1000e3, 500e3]), np.array([-2000e3, 3000e3, 6000e3])): + expected = _gradient(r_body, *coefficients) + got = rotation @ AccelHarmonic(earth, DEGREE, DEGREE)(rotation.T @ r_body, np.zeros(3), T) + assert np.linalg.norm(got - expected) / np.linalg.norm(expected) < 1e-8 + + +def test_moon_harmonics_match_an_independent_potential(): + moon = get_body("moon") + coefficients = _read_tab("gggrx_1200a_sha.tab") + rotation = np.asarray(moon.orientation(T)) + moon_position = np.asarray(moon.position(T)).reshape(3) + for r_body in (np.array([1900e3, 300e3, -200e3]), np.array([-500e3, 800e3, 1700e3])): + expected = _gradient(r_body, *coefficients) + got = rotation @ AccelHarmonic(moon, DEGREE, DEGREE)(moon_position + rotation.T @ r_body, np.zeros(3), T) + assert np.linalg.norm(got - expected) / np.linalg.norm(expected) < 1e-8 From 5ebcb98091afb171e1da2c3f833e9f9f48b90f19 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:19:31 +0000 Subject: [PATCH 17/20] Resolve data files through find_file, skipping git LFS pointers MoonOrientation, MoonPosition, SunPosition and PlanetPosition opened ssapy/data/ directly, and find_file returned whatever file it found. In a clone without 'git lfs pull' those are 130-byte LFS pointers, so jplephem failed with an opaque error about a malformed kernel. find_file now skips LFS pointer files (git-lfs spec v1 header), then searches the optional llnl-ssapy-data package by file name, and if only a pointer exists raises a FileNotFoundError that says so and how to fetch the data. The four ephemeris providers use find_file. With the data present (any release install) nothing changes; once SSAPy-Data ships these files, SSAPy finds them without further code changes. New tests: a pointer-only data directory raises FileNotFoundError naming the LFS pointer; with a pointer locally and a real file of the same name in llnl-ssapy-data, that file is returned. Both fail on the previous code. --- ssapy/body.py | 20 +++++-------- ssapy/utils.py | 55 +++++++++++++++++++++++++++------- tests/test_data_file_lookup.py | 38 +++++++++++++++++++++++ 3 files changed, 91 insertions(+), 22 deletions(-) create mode 100644 tests/test_data_file_lookup.py diff --git a/ssapy/body.py b/ssapy/body.py index 836f1f35..44401948 100644 --- a/ssapy/body.py +++ b/ssapy/body.py @@ -82,11 +82,10 @@ class MoonOrientation(_KernelBacked): orientation matrix at a given time. """ def __init__(self): - import os from jplephem.pck import PCK - from . import datadir + from .utils import find_file - fn = os.path.join(datadir, "moon_pa_de440_200625.bpc") + fn = find_file("moon_pa_de440_200625.bpc") self.kernel = PCK.open(fn) def __call__(self, t, _E=None): @@ -133,11 +132,10 @@ class MoonPosition(_KernelBacked): position vector at a given time. """ def __init__(self): - import os from jplephem.spk import SPK - from . import datadir + from .utils import find_file - fn = os.path.join(datadir, "de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + fn = find_file("de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl self.kernel = SPK.open(fn) def __call__(self, t): @@ -164,11 +162,10 @@ class SunPosition(_KernelBacked): position vector at a given time. """ def __init__(self): - import os from jplephem.spk import SPK - from . import datadir + from .utils import find_file - fn = os.path.join(datadir, "de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + fn = find_file("de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl self.kernel = SPK.open(fn) def __call__(self, t): @@ -197,11 +194,10 @@ class PlanetPosition(_KernelBacked): position vector at a given time. """ def __init__(self, planet_index): - import os from jplephem.spk import SPK - from . import datadir + from .utils import find_file - fn = os.path.join(datadir, "de430.bsp") + fn = find_file("de430.bsp") self.kernel = SPK.open(fn) self.planet_index = planet_index diff --git a/ssapy/utils.py b/ssapy/utils.py index cb7ca84a..26791f33 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -31,22 +31,57 @@ def _emcee_version_before_3(version_string): import astropy._erfa as erfa +_LFS_POINTER_PREFIX = b"version https://git-lfs.github.com/spec/v1" + + +def _is_lfs_pointer(path): + """True if ``path`` is an un-fetched git LFS pointer rather than data.""" + try: + if os.path.getsize(path) > 1024: + return False + with open(path, "rb") as f: + return f.read(len(_LFS_POINTER_PREFIX)) == _LFS_POINTER_PREFIX + except OSError: + return False + + +def _ssapy_data_files(): + """Paths of the files shipped by the optional ``llnl-ssapy-data`` package.""" + try: + import ssapy_data + except ImportError: + return [] + try: + return [os.fspath(path) for path in ssapy_data.iter_data_files()] + except Exception: + return [] + + def find_file(filename, ext=None): """ Find a file in the current directory or the ssapy datadir. If ext is not None, also try appending ext to the filename. + + Un-fetched git LFS pointer files (a clone without ``git lfs pull``) are + skipped. If the file is not found locally, the ``llnl-ssapy-data`` + package is searched by file name when it is installed. """ - candidates = [ - filename, - os.path.join(datadir, filename), - ] - if ext is not None: - candidates.extend([ - filename + ext, - os.path.join(datadir, filename + ext), - ]) + names = [filename] if ext is None else [filename, filename + ext] + candidates = [] + for name in names: + candidates.extend([name, os.path.join(datadir, name)]) for candidate in candidates: - if os.path.isfile(candidate): + if os.path.isfile(candidate) and not _is_lfs_pointer(candidate): return candidate + basenames = {os.path.basename(name) for name in names} + for path in _ssapy_data_files(): + if os.path.basename(path) in basenames and os.path.isfile(path) and not _is_lfs_pointer(path): + return path + for candidate in candidates: + if os.path.isfile(candidate): + raise FileNotFoundError( + f"{candidate} is a git LFS pointer, not the data file; run `git lfs pull` " + "in the SSAPy clone or install a release of SSAPy." + ) raise FileNotFoundError(filename) diff --git a/tests/test_data_file_lookup.py b/tests/test_data_file_lookup.py new file mode 100644 index 00000000..e2d8f964 --- /dev/null +++ b/tests/test_data_file_lookup.py @@ -0,0 +1,38 @@ +import os + +import pytest + +from ssapy import utils + +POINTER = ( + "version https://git-lfs.github.com/spec/v1\n" + "oid sha256:0123456789abcdef0123456789abcdef0123456789abcdef0123456789abcdef\n" + "size 119741440\n" +) + + +def test_lfs_pointer_is_not_returned_as_data(tmp_path, monkeypatch): + # A clone without `git lfs pull` holds 130-byte pointer files in place of + # the kernels; find_file used to return them, and jplephem then failed on + # an unreadable "kernel". It must now raise a FileNotFoundError that says + # what the file is (git LFS spec v1 pointer format). + (tmp_path / "de430.bsp").write_text(POINTER) + monkeypatch.setattr(utils, "datadir", str(tmp_path)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: []) + with pytest.raises(FileNotFoundError, match="git LFS pointer"): + utils.find_file("de430.bsp") + + +def test_ssapy_data_is_searched_when_the_package_data_is_missing(tmp_path, monkeypatch): + # With the package copy only a pointer, a real file of the same name + # shipped by llnl-ssapy-data is returned instead. + local = tmp_path / "local" + shipped = tmp_path / "shipped" / "ephemerides" + local.mkdir() + shipped.mkdir(parents=True) + (local / "de430.bsp").write_text(POINTER) + real = shipped / "de430.bsp" + real.write_bytes(b"DAF/SPK " + b"\0" * 2048) + monkeypatch.setattr(utils, "datadir", str(local)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: [str(real)]) + assert os.path.samefile(utils.find_file("de430.bsp"), real) From c71ba64ba19add3c0bd8fdf388a34f877dcea2f8 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Thu, 8 Oct 2026 23:41:28 +0000 Subject: [PATCH 18/20] Load SSAPy's data from llnl-ssapy-data, with DE440, and drop it here SSAPy stored 21 data files (ephemerides, gravity models, textures; 308 MB) in ssapy/data through Git LFS. They now come from the required llnl-ssapy-data >= 0.2.0 package, which carries them under ssapy/, and the repository no longer uses Git LFS: - ssapy/data and .gitattributes are removed; MANIFEST.in, setup.py and pyproject.toml no longer package data, and the sdist drops from ~300 MB to 441 kB. - llnl-ssapy-data>=0.2.0 is a dependency (pyproject.toml, requirements.txt). - ssapy.datadir is the ssapy/ tree of llnl-ssapy-data, so find_file and every loader (and SSAPy-Toolkit's find_file calls for earth.png, moon.png and Earth_graphics) resolve there. A missing file names the package to install. - The planetary ephemeris moves from JPL DE430 to DE440's short kernel, de440s.bsp (1849-12-26 to 2150-01-22). SSAPy's de430.bsp (119.7 MB) exceeds GitHub's 100 MB file limit, and DE440 matches the DE440 lunar orientation kernel SSAPy already uses. Over 1975-2050, geocentric positions change by at most 10 m for the Moon and 0.5 km for the Sun (103 m and 2.9 km over 1900-2150). A full de440.bsp (1550-2650) in the working directory is used instead when present (body.PLANETARY_EPHEMERIS_FILES). - CI and the publish workflow no longer install or pull Git LFS, so CI now runs the data-dependent tests it used to skip on LFS pointers. - The tests' data checks use find_file instead of inspecting LFS pointers; skills.md describes the new layout. New tests: datadir is ssapy/ inside llnl-ssapy-data, the orientation, gravity and texture files resolve there, the ephemeris resolves to de440s.bsp, and a missing file's error names llnl-ssapy-data. Requires llnl-ssapy-data 0.2.0 to be published first; its wheel is 184 MB, above PyPI's default 100 MB limit, so that release needs a PyPI file-size increase. --- .gitattributes | 9 - .github/workflows/ci.yml | 14 +- .github/workflows/python-publish.yml | 4 +- MANIFEST.in | 1 - pyproject.toml | 3 +- requirements.txt | 1 + setup.py | 1 - skills.md | 12 +- ssapy/__init__.py | 15 +- ssapy/body.py | 30 +- .../Earth_graphics/ne_50m_ocean.README.html | 395 ------------------ .../Earth_graphics/ne_50m_ocean.VERSION.txt | 1 - ssapy/data/Earth_graphics/ne_50m_ocean.cpg | 1 - ssapy/data/Earth_graphics/ne_50m_ocean.dbf | 3 - ssapy/data/Earth_graphics/ne_50m_ocean.prj | 1 - ssapy/data/Earth_graphics/ne_50m_ocean.shp | 3 - ssapy/data/Earth_graphics/ne_50m_ocean.shx | 3 - ssapy/data/de430.bsp | 3 - ssapy/data/earth.png | 3 - ssapy/data/egm2008.egm | 3 - ssapy/data/egm2008.egm.cof | 3 - ssapy/data/egm84.egm | 3 - ssapy/data/egm84.egm.cof | 3 - ssapy/data/egm96.egm | 3 - ssapy/data/egm96.egm.cof | 3 - ssapy/data/gggrx_1200a_sha.lbl | 3 - ssapy/data/gggrx_1200a_sha.tab | 3 - ssapy/data/moon.png | 3 - ssapy/data/moon_pa_de440_200625.bpc | 3 - ssapy/data/wgs84.egm | 3 - ssapy/data/wgs84.egm.cof | 3 - ssapy/utils.py | 20 +- tests/test_accel.py | 30 +- tests/test_body.py | 3 +- tests/test_data_file_lookup.py | 32 +- tests/test_orbit.py | 15 +- tests/test_tle_drag.py | 32 +- 37 files changed, 135 insertions(+), 536 deletions(-) delete mode 100644 .gitattributes delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.README.html delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.cpg delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.dbf delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.prj delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.shp delete mode 100644 ssapy/data/Earth_graphics/ne_50m_ocean.shx delete mode 100644 ssapy/data/de430.bsp delete mode 100644 ssapy/data/earth.png delete mode 100644 ssapy/data/egm2008.egm delete mode 100644 ssapy/data/egm2008.egm.cof delete mode 100644 ssapy/data/egm84.egm delete mode 100644 ssapy/data/egm84.egm.cof delete mode 100644 ssapy/data/egm96.egm delete mode 100644 ssapy/data/egm96.egm.cof delete mode 100644 ssapy/data/gggrx_1200a_sha.lbl delete mode 100644 ssapy/data/gggrx_1200a_sha.tab delete mode 100644 ssapy/data/moon.png delete mode 100644 ssapy/data/moon_pa_de440_200625.bpc delete mode 100644 ssapy/data/wgs84.egm delete mode 100644 ssapy/data/wgs84.egm.cof diff --git a/.gitattributes b/.gitattributes deleted file mode 100644 index b94f62d2..00000000 --- a/.gitattributes +++ /dev/null @@ -1,9 +0,0 @@ -*.png filter=lfs diff=lfs merge=lfs -text -*.jpg filter=lfs diff=lfs merge=lfs -text -*.egm* filter=lfs diff=lfs merge=lfs -text -gggrx_1200a_sha* filter=lfs diff=lfs merge=lfs -text -*.bpc filter=lfs diff=lfs merge=lfs -text -*.bsp filter=lfs diff=lfs merge=lfs -text -*.dbf filter=lfs diff=lfs merge=lfs -text -*.shp filter=lfs diff=lfs merge=lfs -text -*.shx filter=lfs diff=lfs merge=lfs -text diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e9c9deea..c0d56216 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,7 +23,6 @@ concurrency: env: CI_PROJECT_DIR: ${{ github.workspace }} - GIT_LFS_SKIP_SMUDGE: 1 jobs: test: @@ -39,7 +38,6 @@ jobs: - uses: actions/checkout@v4 with: submodules: 'recursive' - lfs: false token: ${{ github.token }} - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v5 @@ -50,15 +48,11 @@ jobs: run: | # Install dependencies sudo apt-get update - sudo apt-get install build-essential git git-lfs python3 python3-setuptools python3-venv graphviz cmake + sudo apt-get install build-essential git python3 python3-setuptools python3-venv graphviz cmake # Initialize and update submodules git submodule update --init --recursive - # Install Git LFS and pull LFS objects - # git lfs install - # git lfs pull - # Upgrade pip and install Python dependencies python -m pip install --upgrade pip setuptools flake8 python -m pip install -r requirements.txt @@ -77,12 +71,6 @@ jobs: # Use the actions/setup-python toolchain and avoid Homebrew updates, # which can fail CI when the hosted runner's brew state is stale. python -m pip install --upgrade pip setuptools flake8 cmake - - # Make LFS usable if the runner already has it, but do not require - # Homebrew just for CI wheel/import checks. - if command -v git-lfs >/dev/null 2>&1; then - git lfs install --skip-repo - fi git config --global credential.helper "" # avoid Keychain prompts # Use the actions/setup-python toolchain consistently diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml index db431709..e951b39d 100644 --- a/.github/workflows/python-publish.yml +++ b/.github/workflows/python-publish.yml @@ -25,9 +25,7 @@ jobs: - name: Install dependencies run: | sudo apt-get update - sudo apt-get install build-essential git git-lfs cmake - git lfs install - git lfs pull + sudo apt-get install build-essential git cmake python -m pip install --upgrade pip build twine python -m pip install -r requirements.txt diff --git a/MANIFEST.in b/MANIFEST.in index ab5c37cf..b5ce3fd2 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,4 +11,3 @@ recursive-include pybind11/include/ * include pybind11/CMakeLists.txt include pybind11/LICENSE recursive-include pybind11/tools * -recursive-include ssapy/data * diff --git a/pyproject.toml b/pyproject.toml index 57448de6..72d4ba66 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ "ipython_genutils", "jplephem", "tqdm", + "llnl-ssapy-data>=0.2.0", ] [project.optional-dependencies] @@ -70,7 +71,7 @@ packages = { find = { include = ["ssapy", "ssapy.*"], namespaces = true } } include-package-data = true [tool.setuptools.package-data] -ssapy = ["*.cpp", "*.h", "data/*"] +ssapy = ["*.cpp", "*.h"] [tool.pytest.ini_options] testpaths = ["tests"] diff --git a/requirements.txt b/requirements.txt index cac06051..c9687c4c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,7 @@ graphviz numpy scipy jplephem +llnl-ssapy-data>=0.2.0 astropy pyerfa emcee diff --git a/setup.py b/setup.py index d8cca2af..5e3da68e 100644 --- a/setup.py +++ b/setup.py @@ -45,7 +45,6 @@ def build_extension(self, ext): ext_modules=[CMakeExtension("ssapy._ssapy")], cmdclass={"build_ext": CMakeBuild}, packages=find_namespace_packages(include=["ssapy", "ssapy.*"]), - package_data={'ssapy': ['data/*']}, zip_safe=False, include_package_data=True, ) diff --git a/skills.md b/skills.md index 19e8ea36..2c7f9cc0 100644 --- a/skills.md +++ b/skills.md @@ -115,11 +115,13 @@ r, v = rv(orbit, times, propagator=KeplerianPropagator()) ## Data and Binary Rules -SSAPy currently includes required package data under `ssapy/data`, and that data -is still part of the base package. Do not migrate or delete it unless the task -explicitly calls for a coordinated data migration. +SSAPy's data (ephemerides, gravity models, textures) comes from the required +`llnl-ssapy-data` package, under its `ssapy/` tree; the repository no longer +stores data or uses Git LFS. The planetary ephemeris is JPL DE440's short +kernel (`de440s.bsp`, 1849-2150); a full `de440.bsp` in the working directory +is used instead when present. -Use `ssapy.datadir` or `ssapy.utils.find_file()` when reading existing bundled +Use `ssapy.datadir` or `ssapy.utils.find_file()` when reading bundled resources. Some model loaders already resolve packaged gravity, SPICE, and Earth orientation files internally; prefer those public loaders over hard-coded paths. @@ -184,7 +186,7 @@ python -m twine check dist/* ## Common Pitfalls -- Do not assume all packaged resources are small; base SSAPy still carries existing required data. +- Do not add data files to this repository; add them to SSAPy-Data and release a new `llnl-ssapy-data`. - Do not add Git LFS as a default solution for new data in this repository. - Do not commit generated `_ssapy*.so`, `build/`, `dist/`, egg-info, docs build output, caches, or local environments. - Do not rely on unbuilt in-tree imports for code paths that require `ssapy._ssapy`. diff --git a/ssapy/__init__.py b/ssapy/__init__.py index 83d10fd0..a11aaa35 100644 --- a/ssapy/__init__.py +++ b/ssapy/__init__.py @@ -2,7 +2,20 @@ __version__ = "1.1.10" -datadir = os.path.join(os.path.dirname(__file__), "data") + + +def _datadir(): + """Directory holding SSAPy's data files: ssapy/ in llnl-ssapy-data.""" + try: + from importlib.resources import files + return os.fspath(files("ssapy_data") / "data" / "ssapy") + except (ImportError, TypeError): + # llnl-ssapy-data is a required dependency; without it find_file + # raises a FileNotFoundError that says how to install it. + return os.path.join(os.path.dirname(__file__), "data") + + +datadir = _datadir() from . import _ssapy from .orbit import Orbit, EarthObserver, OrbitalObserver diff --git a/ssapy/body.py b/ssapy/body.py index 44401948..c740b477 100644 --- a/ssapy/body.py +++ b/ssapy/body.py @@ -15,6 +15,27 @@ def _close_if_possible(obj): close() +#: Planetary ephemeris kernels, in order of preference. llnl-ssapy-data ships +#: JPL DE440's short kernel ``de440s.bsp`` (1849-12-26 to 2150-01-22), consistent +#: with the DE440 lunar orientation kernel; the full ``de440.bsp`` (1550-2650, +#: 114 MB) placed in the working directory is used instead when present. +PLANETARY_EPHEMERIS_FILES = ("de440.bsp", "de440s.bsp") + + +def _planetary_ephemeris_path(): + from .utils import find_file + for name in PLANETARY_EPHEMERIS_FILES: + try: + return find_file(name) + except FileNotFoundError: + continue + raise FileNotFoundError( + "No planetary ephemeris found ({}); install 'llnl-ssapy-data>=0.2.0'.".format( + ", ".join(PLANETARY_EPHEMERIS_FILES) + ) + ) + + class _KernelBacked: def close(self): kernel = getattr(self, "kernel", None) @@ -133,9 +154,8 @@ class MoonPosition(_KernelBacked): """ def __init__(self): from jplephem.spk import SPK - from .utils import find_file - fn = find_file("de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + fn = _planetary_ephemeris_path() self.kernel = SPK.open(fn) def __call__(self, t): @@ -163,9 +183,8 @@ class SunPosition(_KernelBacked): """ def __init__(self): from jplephem.spk import SPK - from .utils import find_file - fn = find_file("de430.bsp") # https://naif.jpl.nasa.gov/pub/naif/LUCY/kernels/spk/de430s.bsp.lbl + fn = _planetary_ephemeris_path() self.kernel = SPK.open(fn) def __call__(self, t): @@ -195,9 +214,8 @@ class PlanetPosition(_KernelBacked): """ def __init__(self, planet_index): from jplephem.spk import SPK - from .utils import find_file - fn = find_file("de430.bsp") + fn = _planetary_ephemeris_path() self.kernel = SPK.open(fn) self.planet_index = planet_index diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.README.html b/ssapy/data/Earth_graphics/ne_50m_ocean.README.html deleted file mode 100644 index 9f9b1de5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.README.html +++ /dev/null @@ -1,395 +0,0 @@ - - - - - - - - -Ocean | Natural Earth - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - - -
-
- - - « 1:50m Physical Vectors  - « Downloads  -
-

Ocean

-
-
-
water_thumb
-
Ocean polygon split into contiguous pieces.

-
- Download ocean (450.91 KB) version 4.0.0 -
-

-
-
About

-

Derives from 10m ocean coastlines. Ocean polygon is the same as Bathymetry 0 (sea level).

-

(below) Yucatan peninsula, Cuba, and Hispaniola.

-

water_banner

-

Issues

-

Caspian’s Garabogaz Aylagy and other “lagoons” might better belong in the lake theme than oceans.

-

Dependencies

-

This theme (and it’s inverse cousin, Land) are used to clip the raster imagery.

-

Version History

-
    -
  • - 4.0.0 -
    • -
  • - 2.0.0 -
    • -
  • - 1.3.0 -
    • -
  • - 1.1.0 -
    • -
  • - 1.0.0 -
    • -
- -

The master changelog is available on Github » -

- - -
- -
- - -
- - - - -
-
- - - - - - \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt b/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt deleted file mode 100644 index 5dbaad6d..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.VERSION.txt +++ /dev/null @@ -1 +0,0 @@ -4.1.0 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.cpg b/ssapy/data/Earth_graphics/ne_50m_ocean.cpg deleted file mode 100644 index 3ad133c0..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.cpg +++ /dev/null @@ -1 +0,0 @@ -UTF-8 \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.dbf b/ssapy/data/Earth_graphics/ne_50m_ocean.dbf deleted file mode 100644 index cb8ca214..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.dbf +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:162d515911b17e35d36e22c5757e75fd37bdd38db2bc0011b9c79888a28bb4f5 -size 176 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.prj b/ssapy/data/Earth_graphics/ne_50m_ocean.prj deleted file mode 100644 index a30c00a5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.prj +++ /dev/null @@ -1 +0,0 @@ -GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] \ No newline at end of file diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.shp b/ssapy/data/Earth_graphics/ne_50m_ocean.shp deleted file mode 100644 index 014215f5..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.shp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6ff1190e8430e493cf8eed99393f021b7562b600b9a7a11aa23ce45b6035d8b8 -size 978340 diff --git a/ssapy/data/Earth_graphics/ne_50m_ocean.shx b/ssapy/data/Earth_graphics/ne_50m_ocean.shx deleted file mode 100644 index 099819b2..00000000 --- a/ssapy/data/Earth_graphics/ne_50m_ocean.shx +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:449353f638399c48c4cfe59e8f945e5c3025a4e412e9bcd0eeba64ba137ae7c1 -size 108 diff --git a/ssapy/data/de430.bsp b/ssapy/data/de430.bsp deleted file mode 100644 index 533d4b74..00000000 --- a/ssapy/data/de430.bsp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6e1b277c5f07135a84950604b83e56b736be696a7f3560bcddb1d4aeb944fca1 -size 119741440 diff --git a/ssapy/data/earth.png b/ssapy/data/earth.png deleted file mode 100644 index 16d9db60..00000000 --- a/ssapy/data/earth.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:6e8a1296caef991dcdba963fdc72816e4095720d005e7af6688b060e5fb85aa4 -size 8825662 diff --git a/ssapy/data/egm2008.egm b/ssapy/data/egm2008.egm deleted file mode 100644 index 9442ced7..00000000 --- a/ssapy/data/egm2008.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:031ecae647986faadf3f3a0177ba07ffcb73a2504b084c2e1c82b97cd9e61c75 -size 790 diff --git a/ssapy/data/egm2008.egm.cof b/ssapy/data/egm2008.egm.cof deleted file mode 100644 index 4268c37d..00000000 --- a/ssapy/data/egm2008.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:ebed5b09e3da7aeae37ac49e7b8c803fda0a525462bba374f762d6345624d132 -size 75755304 diff --git a/ssapy/data/egm84.egm b/ssapy/data/egm84.egm deleted file mode 100644 index a6512e8b..00000000 --- a/ssapy/data/egm84.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:549e21625b71a0f0daec8477e0f4928338885ed2269ff70e351b15c3e6e0d43e -size 800 diff --git a/ssapy/data/egm84.egm.cof b/ssapy/data/egm84.egm.cof deleted file mode 100644 index 2728f31d..00000000 --- a/ssapy/data/egm84.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:f95fa2967091b9a0630ae7f34934446f17a02f73e5b648ea919db4174b047633 -size 262112 diff --git a/ssapy/data/egm96.egm b/ssapy/data/egm96.egm deleted file mode 100644 index d6f2e3b7..00000000 --- a/ssapy/data/egm96.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:cdf0896bb34c1f8ad7ea1c4126cf949fe5069ad9df17c2a66cf3d81da096970e -size 753 diff --git a/ssapy/data/egm96.egm.cof b/ssapy/data/egm96.egm.cof deleted file mode 100644 index 8323c719..00000000 --- a/ssapy/data/egm96.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:fd427a88a944e2df46d3a3e1aefdb17a68e892e3decdf1382f3ca14de9b6f2ae -size 2085160 diff --git a/ssapy/data/gggrx_1200a_sha.lbl b/ssapy/data/gggrx_1200a_sha.lbl deleted file mode 100644 index bb784ae3..00000000 --- a/ssapy/data/gggrx_1200a_sha.lbl +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8d46351c2e73db699d4669a98da1f47551086cbee321d8f9baa3fdb6b65d3d5f -size 10637 diff --git a/ssapy/data/gggrx_1200a_sha.tab b/ssapy/data/gggrx_1200a_sha.tab deleted file mode 100644 index 0d404d19..00000000 --- a/ssapy/data/gggrx_1200a_sha.tab +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:fa04c3dce9376948ad243f3df74144e2602f12d183ea4d179604ed0a79da7ded -size 88059844 diff --git a/ssapy/data/moon.png b/ssapy/data/moon.png deleted file mode 100644 index 754fbf09..00000000 --- a/ssapy/data/moon.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4d9e6eebd180262aefe7555a0231927e5142954144cfcf2c4a10876b7c218027 -size 13728894 diff --git a/ssapy/data/moon_pa_de440_200625.bpc b/ssapy/data/moon_pa_de440_200625.bpc deleted file mode 100644 index ed2f32ac..00000000 --- a/ssapy/data/moon_pa_de440_200625.bpc +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:60cd55aa401ea2ea97360636f567554bfe4e37bb829f901b4460a455dfaf783f -size 12863488 diff --git a/ssapy/data/wgs84.egm b/ssapy/data/wgs84.egm deleted file mode 100644 index c27532ed..00000000 --- a/ssapy/data/wgs84.egm +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:594ff1c9dc2363e50d663003deb9baca2e144d096d70fc9b7a6c1c86f5b781f0 -size 724 diff --git a/ssapy/data/wgs84.egm.cof b/ssapy/data/wgs84.egm.cof deleted file mode 100644 index 0734d247..00000000 --- a/ssapy/data/wgs84.egm.cof +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:aeedf849678826d171b731eb7a783ab5d7e47552413da551418d21445fc3a2ec -size 192 diff --git a/ssapy/utils.py b/ssapy/utils.py index 26791f33..0cb6d1a4 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -58,12 +58,13 @@ def _ssapy_data_files(): def find_file(filename, ext=None): - """ Find a file in the current directory or the ssapy datadir. If ext is - not None, also try appending ext to the filename. + """ Find a file in the current directory or the ssapy datadir (the + ``ssapy/`` tree of the ``llnl-ssapy-data`` package). If ext is not None, + also try appending ext to the filename. - Un-fetched git LFS pointer files (a clone without ``git lfs pull``) are - skipped. If the file is not found locally, the ``llnl-ssapy-data`` - package is searched by file name when it is installed. + Git LFS pointer files (left over from SSAPy versions that stored data + with Git LFS) are skipped, and the rest of ``llnl-ssapy-data`` is searched + by file name. """ names = [filename] if ext is None else [filename, filename + ext] candidates = [] @@ -79,10 +80,13 @@ def find_file(filename, ext=None): for candidate in candidates: if os.path.isfile(candidate): raise FileNotFoundError( - f"{candidate} is a git LFS pointer, not the data file; run `git lfs pull` " - "in the SSAPy clone or install a release of SSAPy." + f"{candidate} is a git LFS pointer, not the data file. SSAPy's data now " + "comes from the llnl-ssapy-data package: pip install 'llnl-ssapy-data>=0.2.0'." ) - raise FileNotFoundError(filename) + raise FileNotFoundError( + f"{filename} was not found in the working directory or in llnl-ssapy-data " + f"({datadir}); install it with pip install 'llnl-ssapy-data>=0.2.0'." + ) def _wrapToPi(angle): diff --git a/tests/test_accel.py b/tests/test_accel.py index a3a7f795..33b19e22 100644 --- a/tests/test_accel.py +++ b/tests/test_accel.py @@ -4,7 +4,6 @@ import sys import subprocess from types import SimpleNamespace -from pathlib import Path import pytest import erfa @@ -18,20 +17,29 @@ from ssapy.constants import EARTH_MU, EARTH_RADIUS -def _is_lfs_pointer(path): + +def _has_data(name): + # SSAPy's data comes from llnl-ssapy-data; skip only if it is missing. + try: + ssapy.utils.find_file(name) + return True + except FileNotFoundError: + return False + + +def _has_ephemeris(): try: - with open(path, "rb") as f: - first = f.readline() - return first.startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: return False -HAS_EGM84 = not _is_lfs_pointer(Path(ssapy.datadir) / "egm84.egm.cof") -HAS_DE430 = not _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp") -HAS_MOON_PA = not _is_lfs_pointer(Path(ssapy.datadir) / "moon_pa_de440_200625.bpc") +HAS_EGM84 = _has_data("egm84.egm.cof") +HAS_EPHEMERIS = _has_ephemeris() +HAS_MOON_PA = _has_data("moon_pa_de440_200625.bpc") -HAS_BODY_DATA = HAS_EGM84 and HAS_DE430 and HAS_MOON_PA +HAS_BODY_DATA = HAS_EGM84 and HAS_EPHEMERIS and HAS_MOON_PA def test_harmonic_coefficients_error_and_identity_paths(tmp_path): @@ -143,7 +151,7 @@ def _make_body_dependent_objects(): earth = moon = sun = earth_MG = moon_MG = sun_MG = None -@pytest.mark.skipif(not HAS_EGM84, reason="EGM84 gravity file unavailable (Git LFS pointer)") +@pytest.mark.skipif(not HAS_EGM84, reason="EGM84 gravity file unavailable (llnl-ssapy-data not installed)") @timer def test_MG_3_1(): """Exercise 3.1 from Montenbruck and Gill diff --git a/tests/test_body.py b/tests/test_body.py index 5069e5f5..8beabb10 100644 --- a/tests/test_body.py +++ b/tests/test_body.py @@ -1,3 +1,4 @@ +import os import numpy as np import pytest @@ -163,7 +164,7 @@ def fake_open(path): assert planet_position.kernel == "planet-kernel" assert planet_position.planet_index == 5 - assert opened[0].endswith("de430.bsp") + assert os.path.basename(opened[0]) in body.PLANETARY_EPHEMERIS_FILES def test_orientation_helpers(monkeypatch): diff --git a/tests/test_data_file_lookup.py b/tests/test_data_file_lookup.py index e2d8f964..1472528d 100644 --- a/tests/test_data_file_lookup.py +++ b/tests/test_data_file_lookup.py @@ -16,11 +16,11 @@ def test_lfs_pointer_is_not_returned_as_data(tmp_path, monkeypatch): # the kernels; find_file used to return them, and jplephem then failed on # an unreadable "kernel". It must now raise a FileNotFoundError that says # what the file is (git LFS spec v1 pointer format). - (tmp_path / "de430.bsp").write_text(POINTER) + (tmp_path / "de440s.bsp").write_text(POINTER) monkeypatch.setattr(utils, "datadir", str(tmp_path)) monkeypatch.setattr(utils, "_ssapy_data_files", lambda: []) - with pytest.raises(FileNotFoundError, match="git LFS pointer"): - utils.find_file("de430.bsp") + with pytest.raises(FileNotFoundError, match="git LFS pointer.*llnl-ssapy-data"): + utils.find_file("de440s.bsp") def test_ssapy_data_is_searched_when_the_package_data_is_missing(tmp_path, monkeypatch): @@ -30,9 +30,29 @@ def test_ssapy_data_is_searched_when_the_package_data_is_missing(tmp_path, monke shipped = tmp_path / "shipped" / "ephemerides" local.mkdir() shipped.mkdir(parents=True) - (local / "de430.bsp").write_text(POINTER) - real = shipped / "de430.bsp" + (local / "de440s.bsp").write_text(POINTER) + real = shipped / "de440s.bsp" real.write_bytes(b"DAF/SPK " + b"\0" * 2048) monkeypatch.setattr(utils, "datadir", str(local)) monkeypatch.setattr(utils, "_ssapy_data_files", lambda: [str(real)]) - assert os.path.samefile(utils.find_file("de430.bsp"), real) + assert os.path.samefile(utils.find_file("de440s.bsp"), real) + + +def test_datadir_is_the_ssapy_tree_of_llnl_ssapy_data(): + # SSAPy no longer ships data: datadir is ssapy/ inside llnl-ssapy-data, and + # every file the body and gravity loaders need resolves there. + import ssapy + import ssapy_data + from ssapy.body import _planetary_ephemeris_path + + assert os.path.samefile(ssapy.datadir, os.path.join(os.path.dirname(ssapy_data.__file__), "data", "ssapy")) + for name in ("moon_pa_de440_200625.bpc", "egm84.egm", "egm2008.egm.cof", "gggrx_1200a_sha.tab", "earth.png", "moon.png"): + assert os.path.dirname(utils.find_file(name)) == os.path.normpath(ssapy.datadir) + assert os.path.basename(_planetary_ephemeris_path()) == "de440s.bsp" + + +def test_missing_file_names_the_data_package(monkeypatch, tmp_path): + monkeypatch.setattr(utils, "datadir", str(tmp_path)) + monkeypatch.setattr(utils, "_ssapy_data_files", lambda: []) + with pytest.raises(FileNotFoundError, match="llnl-ssapy-data"): + utils.find_file("no_such_file.bsp") diff --git a/tests/test_orbit.py b/tests/test_orbit.py index ad1f6c47..c0b7c478 100755 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -18,13 +18,13 @@ from . import ssapy_test_helpers from .ssapy_test_helpers import timer, checkAngle, checkSphere, sample_orbit, sample_LEO_orbit, sample_GEO_orbit -from pathlib import Path -def _is_lfs_pointer(path): +def _has_ephemeris(): + # SSAPy's ephemeris comes from llnl-ssapy-data; skip only if it is missing. try: - with open(path, "rb") as f: - return f.readline().startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: return False try: @@ -1570,10 +1570,7 @@ def test_light_time_correction(): np.testing.assert_allclose(dircos, dc_linear, rtol=0, atol=1e-7) np.testing.assert_allclose(dircos, dc_exact, rtol=0, atol=1e-14) -@pytest.mark.skipif( - _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp"), - reason="de430.bsp is a Git LFS pointer, not real ephemeris data", -) +@pytest.mark.skipif(not _has_ephemeris(), reason="planetary ephemeris unavailable (llnl-ssapy-data not installed)") @pytest.mark.timeout(90) @timer def test_find_passes(): diff --git a/tests/test_tle_drag.py b/tests/test_tle_drag.py index 066c00e5..18e102ca 100644 --- a/tests/test_tle_drag.py +++ b/tests/test_tle_drag.py @@ -1,5 +1,4 @@ import warnings -from pathlib import Path import numpy as np import pytest @@ -16,22 +15,31 @@ _sgp4_reference_arc, fit_drag) -def _is_lfs_pointer(path): +# The numerical drag propagator needs the gravity and ephemeris data from +# llnl-ssapy-data; skip those tests if it is not installed. +def _has_data(name): + # SSAPy's data comes from llnl-ssapy-data; skip only if it is missing. try: - with open(path, "rb") as f: - return f.readline().startswith(b"version https://git-lfs.github.com/spec/v1") - except Exception: + ssapy.utils.find_file(name) + return True + except FileNotFoundError: return False -# The numerical drag propagator needs the (Git LFS) gravity + ephemeris data, -# which is absent in CI; skip those tests there, as tests/test_accel.py does. -HAS_EGM84 = not _is_lfs_pointer(Path(ssapy.datadir) / "egm84.egm.cof") -HAS_DE430 = not _is_lfs_pointer(Path(ssapy.datadir) / "de430.bsp") -HAS_MOON_PA = not _is_lfs_pointer(Path(ssapy.datadir) / "moon_pa_de440_200625.bpc") -HAS_BODY_DATA = HAS_EGM84 and HAS_DE430 and HAS_MOON_PA +def _has_ephemeris(): + try: + ssapy.body._planetary_ephemeris_path() + return True + except FileNotFoundError: + return False + + +HAS_EGM84 = _has_data("egm84.egm.cof") +HAS_EPHEMERIS = _has_ephemeris() +HAS_MOON_PA = _has_data("moon_pa_de440_200625.bpc") +HAS_BODY_DATA = HAS_EGM84 and HAS_EPHEMERIS and HAS_MOON_PA needs_body_data = pytest.mark.skipif( - not HAS_BODY_DATA, reason="body data unavailable (Git LFS pointer)") + not HAS_BODY_DATA, reason="body data unavailable (llnl-ssapy-data not installed)") ISS = ("1 25544U 98067A 24015.54791435 .00016717 00000-0 30074-3 0 9993", "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.49514637123456") From 460bc01d73f3b2ecd1e3ece6bcf378585f851147 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Fri, 9 Oct 2026 03:26:36 +0000 Subject: [PATCH 19/20] Select DE440 or DE430 and fetch long-span kernels when epochs need them llnl-ssapy-data 0.2.0 ships two short planetary kernels: DE440's de440s.bsp (1849-12-26 to 2150-01-22, the default) and a 1900-2150 excerpt of DE430 (de430_1900_2150.bsp) for reproducing SSAPy <= 1.1.10 results. The new ssapy.ephemeris module selects between them and covers epochs outside them: - set_planetary_ephemeris("de440" | "de430") or SSAPY_EPHEMERIS chooses the ephemeris for Sun, Moon and planet positions created afterwards. In 2026 the two differ by under 10 m for the Moon and 0.5 km for the Sun. - Each epoch is evaluated with the first kernel of the selected ephemeris that covers it: the shipped short kernel, then the full-span kernel of the same solution (de440.bsp or de430.bsp, 1550-2650, 120 MB), then for DE440 the two halves of DE441 (-13200 to 1969 and 1969 to 17191, 1.65 GB each). Shipped and full kernels of one solution hold the same records, so positions are continuous where SSAPy switches. - Kernels not found locally are downloaded from NAIF into ~/.cache/ssapy (SSAPY_DATA_CACHE), with an EphemerisDownloadWarning first, and kept only if their size and SHA-256 match the values pinned here (computed from NAIF's files on 2026-10-09; de430.bsp's equals SSAPy's LFS copy). Using DE441 for a DE440 request also raises an EphemerisRangeWarning, since DE441 is a different solution (no lunar core-mantle damping). - Without network access, or with SSAPY_EPHEMERIS_DOWNLOAD=0, the error names the kernel, its URL and where to put it; ssapy.ephemeris.fetch() prefetches on a networked node for jobs that have none. - Epochs inside the shipped kernel take a single-kernel fast path: 159 us per scalar Moon position against 135 us for direct jplephem. The lunar orientation kernel (moon_pa_de440_200625.bpc) covers 1549-12-31 to 2650-01-25 only, so lunar-frame quantities still fail outside that span. New tests (offline, with a local file standing in for NAIF): DE440 and DE430 Moon and Sun positions agree to 20 m and 1 km in 2026; the environment variable selects DE430; epochs inside the shipped span never download; an epoch outside downloads once with a warning, is cached, and matches direct evaluation of the long kernel to 1e-6 m; a different long-span solution warns; offline and disabled downloads fail with a clear error and leave no partial file; a corrupt download is rejected; epochs in 1700, 2500, 1000 and 3000 route to de440.bsp, de440.bsp, DE441 part 1 and DE441 part 2. An opt-in test (SSAPY_TEST_NETWORK=1) downloads de440.bsp from NAIF and evaluates the Moon in 1700; it passed here. The installation docs describe selection, downloads and prefetching; git-lfs is no longer a prerequisite. --- docs/source/api.rst | 1 + docs/source/installation.rst | 30 +- docs/source/tables/system_prerequisites.csv | 1 - skills.md | 8 +- ssapy/__init__.py | 1 + ssapy/body.py | 63 +++-- ssapy/ephemeris.py | 292 ++++++++++++++++++++ tests/test_body.py | 2 +- tests/test_ephemeris_selection.py | 174 ++++++++++++ 9 files changed, 534 insertions(+), 38 deletions(-) create mode 100644 ssapy/ephemeris.py create mode 100644 tests/test_ephemeris_selection.py diff --git a/docs/source/api.rst b/docs/source/api.rst index f847b50b..b7e0087f 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -20,6 +20,7 @@ and/or heavy dependencies. Below lists the packages available in the ``ssapy`` ssapy.constants ssapy.correlate_tracks ssapy.ellipsoid + ssapy.ephemeris ssapy.gravity ssapy.io ssapy.linker diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 15593cfe..a364e85b 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -19,7 +19,7 @@ These requirements can be easily installed on most modern macOS and Linux system .. code-block:: console apt update - apt install build-essential git git-lfs python3 python3-distutils python3-venv graphviz + apt install build-essential git python3 python3-distutils python3-venv graphviz .. tab:: RHEL @@ -27,14 +27,14 @@ These requirements can be easily installed on most modern macOS and Linux system dnf install epel-release dnf group install "Development Tools" - dnf install git git-lfs gcc-gfortran python3 python3-pip python3-setuptools graphviz + dnf install git gcc-gfortran python3 python3-pip python3-setuptools graphviz .. tab:: macOS Brew .. code-block:: console brew update - brew install gcc git git-lfs python3 graphviz + brew install gcc git python3 graphviz Installation ------------ @@ -45,6 +45,30 @@ As the package has been published on `PyPI =0.2.0'.".format( - ", ".join(PLANETARY_EPHEMERIS_FILES) - ) - ) +class _KernelBacked: + _ephemeris = None + def _positions(self, chain, mjd_tt): + """Sum of SPK segment chains at TT/TDB MJDs, in km, choosing a longer + kernel of the selected ephemeris for epochs the shipped one lacks.""" + if self._ephemeris is None: + from .ephemeris import PlanetaryEphemeris + self._ephemeris = PlanetaryEphemeris() + return self._ephemeris.compute(chain, 2400000.5, mjd_tt, self.kernel) -class _KernelBacked: def close(self): kernel = getattr(self, "kernel", None) if kernel is not None: @@ -151,11 +146,17 @@ def _Rz(alpha): class MoonPosition(_KernelBacked): """Position of moon in GCRF. This is a callable class that returns the position vector at a given time. + + Uses the planetary ephemeris selected with + ``ssapy.ephemeris.set_planetary_ephemeris`` (DE440 by default), switching + to a longer kernel for epochs outside the shipped one. """ def __init__(self): from jplephem.spk import SPK + from .ephemeris import PlanetaryEphemeris - fn = _planetary_ephemeris_path() + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) def __call__(self, t): @@ -172,8 +173,8 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - pos = self.kernel[3, 301].compute(2400000.5, mjd_tt) # Earth-moon barycenter -> moon - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon barycenter -> earth + # Earth-moon barycenter -> moon, minus Earth-moon barycenter -> earth + pos = self._positions(((3, 301, 1), (3, 399, -1)), mjd_tt) return pos * 1e3 @@ -183,8 +184,10 @@ class SunPosition(_KernelBacked): """ def __init__(self): from jplephem.spk import SPK + from .ephemeris import PlanetaryEphemeris - fn = _planetary_ephemeris_path() + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) def __call__(self, t): @@ -201,10 +204,8 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - - pos = self.kernel[0, 10].compute(2400000.5, mjd_tt) # SS bary -> sun - pos -= self.kernel[0, 3].compute(2400000.5, mjd_tt) # SS bary -> Earth-moon bary - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon bary -> Earth + # SS bary -> sun, minus SS bary -> Earth-moon bary, minus Earth-moon bary -> Earth + pos = self._positions(((0, 10, 1), (0, 3, -1), (3, 399, -1)), mjd_tt) return pos * 1e3 @@ -214,8 +215,10 @@ class PlanetPosition(_KernelBacked): """ def __init__(self, planet_index): from jplephem.spk import SPK + from .ephemeris import PlanetaryEphemeris - fn = _planetary_ephemeris_path() + self._ephemeris = PlanetaryEphemeris() + fn = _planetary_ephemeris_path(self._ephemeris.family) self.kernel = SPK.open(fn) self.planet_index = planet_index @@ -233,9 +236,9 @@ def __call__(self, t): Position vector at time t in meters. """ mjd_tt = _gpsToTT(t) - pos = self.kernel[0, self.planet_index].compute(2400000.5, mjd_tt) # SS bary -> Jupiter - pos -= self.kernel[0, 3].compute(2400000.5, mjd_tt) # SS bary -> Earth-moon bary - pos -= self.kernel[3, 399].compute(2400000.5, mjd_tt) # Earth-moon bary -> Earth + # SS bary -> planet barycenter, minus SS bary -> Earth-moon bary, minus Earth-moon bary -> Earth + chain = ((0, self.planet_index, 1), (0, 3, -1), (3, 399, -1)) + pos = self._positions(chain, mjd_tt) return pos * 1e3 diff --git a/ssapy/ephemeris.py b/ssapy/ephemeris.py new file mode 100644 index 00000000..ab5ec8b6 --- /dev/null +++ b/ssapy/ephemeris.py @@ -0,0 +1,292 @@ +""" +Planetary ephemeris selection and long-span kernels. + +SSAPy evaluates Sun, Moon and planet positions from JPL development +ephemerides (DE). ``llnl-ssapy-data`` ships two short kernels: + +========== ======================== ========================= +Ephemeris Shipped kernel Coverage +========== ======================== ========================= +``de440`` ``de440s.bsp`` 1849-12-26 to 2150-01-22 +``de430`` ``de430_1900_2150.bsp`` 1900-01-01 to 2150-01-01 +========== ======================== ========================= + +DE440 is the default. Select DE430, for example to reproduce results made +with SSAPy <= 1.1.10, with :func:`set_planetary_ephemeris` or the +``SSAPY_EPHEMERIS`` environment variable. + +For epochs outside the shipped kernel, SSAPy uses a longer kernel of the same +ephemeris (``de440.bsp`` or ``de430.bsp``, 1550-2650) and, for DE440 beyond +that, DE441 (-13200 to +17191). These are 120 MB to 1.65 GB, so they are not +shipped: SSAPy looks for them in the working directory, ``llnl-ssapy-data`` +and the cache directory, and otherwise downloads them from NAIF, checks their +SHA-256, and caches them, warning with :class:`EphemerisDownloadWarning` +before each download. Set ``SSAPY_EPHEMERIS_DOWNLOAD=0`` to forbid downloads +(e.g. on compute nodes without network access; prefetch on a login node with +:func:`fetch`) and ``SSAPY_DATA_CACHE`` to choose the cache directory. + +The shipped and long kernels of one ephemeris hold the same solution, so +positions are continuous where SSAPy switches between them. DE441 is a +different solution from DE440 (no lunar core-mantle damping; less accurate +in the current century), so using it raises :class:`EphemerisRangeWarning`. +""" +from __future__ import annotations + +import hashlib +import os +import shutil +import tempfile +import warnings +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +__all__ = [ + "EphemerisDownloadWarning", + "EphemerisRangeWarning", + "EphemerisUnavailableError", + "set_planetary_ephemeris", + "get_planetary_ephemeris", + "fetch", + "cache_dir", +] + +NAIF_PLANETS_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/" + + +class EphemerisDownloadWarning(UserWarning): + """A long-span ephemeris kernel is being downloaded.""" + + +class EphemerisRangeWarning(UserWarning): + """An epoch is evaluated with a different solution than the one selected.""" + + +class EphemerisUnavailableError(FileNotFoundError): + """No kernel of the selected ephemeris covers the requested epoch.""" + + +@dataclass(frozen=True) +class Kernel: + name: str + start_jd: float # TDB Julian dates covered by every segment SSAPy uses + end_jd: float + sha256: Optional[str] = None # None for kernels shipped in llnl-ssapy-data + nbytes: Optional[int] = None + solution: str = "" + + @property + def url(self): + return NAIF_PLANETS_URL + self.name + + def covers(self, jd): + return (jd >= self.start_jd) & (jd <= self.end_jd) + + +EPHEMERIDES = { + "de440": ( + Kernel("de440s.bsp", 2396752.5, 2506352.5, solution="DE440"), + Kernel("de440.bsp", 2287184.5, 2688976.5, + "a4ce9bf9b3282becc9f4b2ac3cebe03a2ae7599981aabd7265fd8482fff7c4b5", 119799808, "DE440"), + Kernel("de441_part-2.bsp", 2440400.5, 8000016.5, + "3abb17dae2d78dd34880377544aacb54892104a0d4462b322cb9f4454d4887f6", 1656830976, "DE441"), + Kernel("de441_part-1.bsp", -3100015.5, 2440432.5, + "13757827f5db41b835a24bbd637488636ce79a8ca754062fed17844f7d5b618e", 1651119104, "DE441"), + ), + "de430": ( + Kernel("de430_1900_2150.bsp", 2415020.5, 2506331.5, solution="DE430"), + Kernel("de430.bsp", 2287184.5, 2688976.5, + "6e1b277c5f07135a84950604b83e56b736be696a7f3560bcddb1d4aeb944fca1", 119741440, "DE430"), + ), +} + +_selected = None +_open_kernels = {} + + +def set_planetary_ephemeris(name): + """Select the planetary ephemeris, ``"de440"`` (default) or ``"de430"``. + + Affects Sun, Moon and planet positions created afterwards + (``get_body``, ``SunPosition``, ``MoonPosition``, ``PlanetPosition``). + """ + global _selected + key = str(name).lower() + if key not in EPHEMERIDES: + raise ValueError(f"Unknown planetary ephemeris {name!r}; choose from {sorted(EPHEMERIDES)}") + _selected = key + + +def get_planetary_ephemeris(): + """Name of the selected planetary ephemeris.""" + if _selected is not None: + return _selected + key = os.environ.get("SSAPY_EPHEMERIS", "de440").lower() + if key not in EPHEMERIDES: + raise ValueError(f"SSAPY_EPHEMERIS={key!r} is not one of {sorted(EPHEMERIDES)}") + return key + + +def cache_dir(): + """Directory where downloaded kernels are kept (``SSAPY_DATA_CACHE``).""" + return os.path.expanduser(os.environ.get("SSAPY_DATA_CACHE", os.path.join("~", ".cache", "ssapy"))) + + +def _downloads_allowed(): + return os.environ.get("SSAPY_EPHEMERIS_DOWNLOAD", "1").strip().lower() not in ("0", "false", "no", "off") + + +def _sha256(path): + digest = hashlib.sha256() + with open(path, "rb") as handle: + for chunk in iter(lambda: handle.read(1 << 20), b""): + digest.update(chunk) + return digest.hexdigest() + + +def _open_url(url): # separated so tests can replace the network + from urllib.request import urlopen + return urlopen(url, timeout=60) + + +def _download(kernel, destination): + os.makedirs(os.path.dirname(destination), exist_ok=True) + handle, partial = tempfile.mkstemp(dir=os.path.dirname(destination), suffix=".part") + try: + with os.fdopen(handle, "wb") as out, _open_url(kernel.url) as response: + shutil.copyfileobj(response, out, 1 << 20) + if os.path.getsize(partial) != kernel.nbytes or _sha256(partial) != kernel.sha256: + raise EphemerisUnavailableError( + f"{kernel.name} downloaded from {kernel.url} does not match its recorded SHA-256" + ) + os.replace(partial, destination) + finally: + if os.path.exists(partial): + os.remove(partial) + + +def _local_path(kernel): + from .utils import find_file + try: + return find_file(kernel.name) + except FileNotFoundError: + pass + cached = os.path.join(cache_dir(), kernel.name) + return cached if os.path.isfile(cached) else None + + +def fetch(name=None, include_de441=True): + """Download (if needed) the long-span kernels of an ephemeris and return + their paths. Use on a login node to fill ``SSAPY_DATA_CACHE`` before + running jobs without network access.""" + family = EPHEMERIDES[(name or get_planetary_ephemeris()).lower()] + return [ + _path(kernel, "requested with ssapy.ephemeris.fetch") + for kernel in family[1:] + if include_de441 or kernel.solution == family[0].solution + ] + + +def _path(kernel, reason): + path = _local_path(kernel) + if path is not None: + return path + if kernel.sha256 is None: + raise EphemerisUnavailableError( + f"{kernel.name} is missing; install 'llnl-ssapy-data>=0.2.0'." + ) + destination = os.path.join(cache_dir(), kernel.name) + manual = (f"Download {kernel.url} ({kernel.nbytes / 1e6:.0f} MB) into {cache_dir()} " + f"(or set SSAPY_DATA_CACHE), or place it in the working directory.") + if not _downloads_allowed(): + raise EphemerisUnavailableError(f"{kernel.name} is needed ({reason}) and downloads are disabled " + f"(SSAPY_EPHEMERIS_DOWNLOAD=0). {manual}") + warnings.warn( + f"Downloading {kernel.name} ({kernel.nbytes / 1e6:.0f} MB) from NAIF to {destination}: {reason}.", + EphemerisDownloadWarning, + stacklevel=6, + ) + try: + _download(kernel, destination) + except EphemerisUnavailableError: + raise + except Exception as exc: # no network, proxy, HTTP error, disk full + raise EphemerisUnavailableError(f"Could not download {kernel.name}: {exc}. {manual}") from exc + return destination + + +def _open(kernel, path): + from jplephem.spk import SPK + key = os.path.abspath(path) + if key not in _open_kernels: + _open_kernels[key] = SPK.open(path) + return _open_kernels[key] + + +def _jd_text(jd): + from astropy.time import Time + return Time(float(jd), format="jd", scale="tdb").iso[:10] + + +class PlanetaryEphemeris: + """Evaluate SPK segment chains, choosing a kernel per epoch. + + The shipped kernel of ``family`` (passed to :meth:`compute` as + ``primary`` once opened) is used for every epoch it covers; other epochs + go to the first longer kernel of the family that covers them. + """ + + def __init__(self, family=None): + self.family = (family or get_planetary_ephemeris()).lower() + self.kernels = EPHEMERIDES[self.family] + + def _kernel_for(self, kernel, reason, primary): + if kernel is self.kernels[0] and primary is not None: + return primary + return _open(kernel, _path(kernel, reason)) + + def shipped_path(self): + return _path(self.kernels[0], "shipped kernel") + + def compute(self, chain, jd1, jd2, primary=None): + """Sum ``sign * kernel[center, target].compute(jd1, jd2)`` over the + ``(center, target, sign)`` terms in ``chain``; returns km.""" + jd2 = np.asarray(jd2, dtype=float) + shipped = self.kernels[0] + if primary is not None and np.all(shipped.covers(jd1 + jd2)): + # Fast path (every call in normal use): one kernel, no masking. + total = 0.0 + for center, target, sign in chain: + total = total + sign * np.asarray(primary[center, target].compute(jd1, jd2)) + return total + scalar = jd2.ndim == 0 + jd2 = np.atleast_1d(jd2) + jd = jd1 + jd2 + out = np.zeros((3, jd2.size)) + pending = np.ones(jd2.size, dtype=bool) + for kernel in self.kernels: + mask = pending & kernel.covers(jd) + if not np.any(mask): + continue + reason = None + if kernel is not shipped: + reason = (f"epochs {_jd_text(jd[mask].min())} to {_jd_text(jd[mask].max())} are outside " + f"{shipped.name} ({_jd_text(shipped.start_jd)} to {_jd_text(shipped.end_jd)})") + spk = self._kernel_for(kernel, reason, primary) + if kernel.solution != self.kernels[0].solution: + warnings.warn( + f"Using {kernel.solution} ({kernel.name}) for {reason}; {kernel.solution} differs from " + f"{self.kernels[0].solution}, and the switch makes positions jump slightly at its boundary.", + EphemerisRangeWarning, + stacklevel=4, + ) + for center, target, sign in chain: + out[:, mask] += sign * np.asarray(spk[center, target].compute(jd1, jd2[mask])).reshape(3, -1) + pending &= ~mask + if np.any(pending): + raise EphemerisUnavailableError( + f"No {self.family.upper()} kernel covers {_jd_text(jd[pending].min())} to " + f"{_jd_text(jd[pending].max())}." + ) + return out[:, 0] if scalar else out diff --git a/tests/test_body.py b/tests/test_body.py index 8beabb10..96bed091 100644 --- a/tests/test_body.py +++ b/tests/test_body.py @@ -164,7 +164,7 @@ def fake_open(path): assert planet_position.kernel == "planet-kernel" assert planet_position.planet_index == 5 - assert os.path.basename(opened[0]) in body.PLANETARY_EPHEMERIS_FILES + assert os.path.basename(opened[0]) == "de440s.bsp" def test_orientation_helpers(monkeypatch): diff --git a/tests/test_ephemeris_selection.py b/tests/test_ephemeris_selection.py new file mode 100644 index 00000000..d4e1ce21 --- /dev/null +++ b/tests/test_ephemeris_selection.py @@ -0,0 +1,174 @@ +"""Planetary ephemeris selection, long-span downloads, and their failure modes. + +The download tests stand in a local file for NAIF, so they run offline: a +2000-2050 excerpt of the shipped de440s.bsp plays the shipped kernel and the +full de440s.bsp plays the long-span kernel "downloaded" for other epochs. +""" +import io +import os +import subprocess +import sys +import warnings + +import numpy as np +import pytest +from astropy.time import Time +from jplephem.spk import SPK + +import ssapy +from ssapy import ephemeris +from ssapy.body import MoonPosition, SunPosition +from ssapy.utils import find_file + + +@pytest.fixture(autouse=True) +def _reset_selection(monkeypatch, tmp_path): + monkeypatch.setattr(ephemeris, "_selected", None) + monkeypatch.delenv("SSAPY_EPHEMERIS", raising=False) + monkeypatch.delenv("SSAPY_EPHEMERIS_DOWNLOAD", raising=False) + monkeypatch.setenv("SSAPY_DATA_CACHE", str(tmp_path / "cache")) + monkeypatch.setattr(ephemeris, "_open_kernels", {}) + + +def _gps(iso): + # TT avoids ERFA's "dubious year" warnings for UTC far from the present. + return Time(iso, scale="tt").gps + + +def test_de440_and_de430_agree_to_their_published_difference(): + # DE440 is the default; DE430 is selectable. In 2026 the two solutions' + # geocentric Moon differs by <= 10 m and Sun by <= 0.5 km over 1975-2050 + # (measured from the NAIF kernels); hold them to 20 m and 1 km. + t = _gps("2026-10-08T00:00:00") + assert ephemeris.get_planetary_ephemeris() == "de440" + moon440, sun440 = MoonPosition()(t), SunPosition()(t) + ephemeris.set_planetary_ephemeris("de430") + moon430, sun430 = MoonPosition()(t), SunPosition()(t) + assert 0 < np.linalg.norm(moon440 - moon430) < 20.0 + assert np.linalg.norm(sun440 - sun430) < 1e3 + with pytest.raises(ValueError): + ephemeris.set_planetary_ephemeris("de405") + + +def test_environment_variable_selects_the_ephemeris(monkeypatch): + monkeypatch.setenv("SSAPY_EPHEMERIS", "DE430") + assert ephemeris.get_planetary_ephemeris() == "de430" + assert os.path.basename(ssapy.body._planetary_ephemeris_path()) == "de430_1900_2150.bsp" + + +def test_no_download_inside_the_shipped_span(monkeypatch): + def no_network(url): + raise AssertionError("no download expected") + + monkeypatch.setattr(ephemeris, "_open_url", no_network) + with warnings.catch_warnings(): + warnings.simplefilter("error") + MoonPosition()(_gps(["1900-01-01T00:00:00", "2026-10-08T00:00:00", "2149-12-31T00:00:00"])) + + +def _harness(monkeypatch, tmp_path, *, long_solution="DE440", served=None): + """Shipped kernel = 2000-2050 excerpt of de440s; long kernel = de440s.""" + shipped_full = find_file("de440s.bsp") + excerpt = tmp_path / "short_test.bsp" + subprocess.run([sys.executable, "-m", "jplephem", "excerpt", "2000/1/1", "2050/1/1", shipped_full, str(excerpt)], + check=True, stdout=subprocess.DEVNULL) + payload = open(shipped_full, "rb").read() if served is None else served + monkeypatch.chdir(tmp_path) + full = SPK.open(shipped_full) + start = max(seg.start_jd for seg in full.segments) + end = min(seg.end_jd for seg in full.segments) + family = ( + ephemeris.Kernel("short_test.bsp", 2451544.5, 2469807.5, solution="DE440"), + ephemeris.Kernel("long_test.bsp", start, end, ephemeris._sha256(shipped_full), + os.path.getsize(shipped_full), long_solution), + ) + monkeypatch.setitem(ephemeris.EPHEMERIDES, "de440", family) + calls = [] + + def fake_open_url(url): + calls.append(url) + return io.BytesIO(payload) + + monkeypatch.setattr(ephemeris, "_open_url", fake_open_url) + return full, calls + + +def _direct_moon(kernel, t): + mjd = ssapy.utils._gpsToTT(t) + return (kernel[3, 301].compute(2400000.5, mjd) - kernel[3, 399].compute(2400000.5, mjd)) * 1e3 + + +def test_epochs_outside_the_shipped_span_download_with_a_warning(monkeypatch, tmp_path): + # An epoch beyond the shipped span downloads the long kernel once, warns, + # caches it, and evaluates it exactly; epochs inside stay on the shipped + # kernel, and the two agree at the boundary (same solution). + full, calls = _harness(monkeypatch, tmp_path) + moon = MoonPosition() + times = _gps(["2026-10-08T00:00:00", "2100-06-01T00:00:00"]) + with pytest.warns(ephemeris.EphemerisDownloadWarning, match="long_test.bsp"): + got = moon(times) + assert len(calls) == 1 and calls[0].endswith("long_test.bsp") + assert os.path.isfile(tmp_path / "cache" / "long_test.bsp") + np.testing.assert_allclose(got[:, 1], _direct_moon(full, times[1]), rtol=0, atol=1e-6) + np.testing.assert_allclose(got[:, 0], _direct_moon(full, times[0]), rtol=0, atol=1e-3) + with warnings.catch_warnings(): + warnings.simplefilter("error") + moon(times) # cached: no second download or warning + assert len(calls) == 1 + + +def test_a_different_long_span_solution_warns(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path, long_solution="DE441") + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.warns(ephemeris.EphemerisRangeWarning, match="DE441"): + MoonPosition()(_gps("2100-06-01T00:00:00")) + + +@pytest.mark.parametrize("iso, kernel", [ + ("1700-01-01T00:00:00", "de440.bsp"), + ("2500-01-01T00:00:00", "de440.bsp"), + ("1000-01-01T00:00:00", "de441_part-1.bsp"), + ("3000-01-01T00:00:00", "de441_part-2.bsp"), +]) +def test_long_span_epochs_route_to_full_de440_then_de441(monkeypatch, iso, kernel): + # With downloads disabled the error names the kernel that would be fetched: + # full DE440 for 1550-2650 outside de440s, DE441 beyond (NAIF coverage). + monkeypatch.setenv("SSAPY_EPHEMERIS_DOWNLOAD", "0") + with pytest.raises(ephemeris.EphemerisUnavailableError, match=kernel.replace(".", r"\.")): + MoonPosition()(_gps(iso)) + + +def test_offline_or_disabled_downloads_fail_clearly(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path) + t = _gps("2100-06-01T00:00:00") + + def offline(url): + raise OSError("network unreachable") + + monkeypatch.setattr(ephemeris, "_open_url", offline) + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.raises(ephemeris.EphemerisUnavailableError, match="Could not download long_test.bsp"): + MoonPosition()(t) + assert not os.listdir(tmp_path / "cache") # no partial file left + + monkeypatch.setenv("SSAPY_EPHEMERIS_DOWNLOAD", "0") + with pytest.raises(ephemeris.EphemerisUnavailableError, match="downloads are disabled"): + MoonPosition()(t) + + +def test_a_corrupt_download_is_rejected(monkeypatch, tmp_path): + _harness(monkeypatch, tmp_path, served=b"not a kernel") + with pytest.warns(ephemeris.EphemerisDownloadWarning): + with pytest.raises(ephemeris.EphemerisUnavailableError, match="SHA-256"): + MoonPosition()(_gps("2100-06-01T00:00:00")) + assert not os.path.exists(tmp_path / "cache" / "long_test.bsp") + + +@pytest.mark.skipif(os.environ.get("SSAPY_TEST_NETWORK") != "1", reason="set SSAPY_TEST_NETWORK=1 to download from NAIF") +@pytest.mark.timeout(600) +def test_real_download_of_full_de440(): + # Downloads de440.bsp (119.8 MB) from NAIF, checks its SHA-256, and + # evaluates the Moon in 1700, outside de440s.bsp. + with pytest.warns(ephemeris.EphemerisDownloadWarning, match="de440.bsp"): + position = MoonPosition()(_gps("1700-01-01T00:00:00")) + assert 3.5e8 < np.linalg.norm(position) < 4.1e8 From d0a9cb719132e6c54d96a82e74e9cb78e4a45783 Mon Sep 17 00:00:00 2001 From: Travis Yeager Date: Fri, 9 Oct 2026 11:01:05 -0700 Subject: [PATCH 20/20] Use split SSATK data distributions --- README.rst | 3 +++ pyproject.toml | 17 ++++++++++++++++- 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 2d938873..f8a614c5 100644 --- a/README.rst +++ b/README.rst @@ -29,6 +29,9 @@ wrappers around coordinate conversions, Lambertian magnitude / brightness calculations, and related extensions, see the companion project `SSAPy-Toolkit `__. +Optional packaged datasets are available with ``llnl-ssapy[data]``. Install +``llnl-ssapy[all-data]`` to install every split ``ssatk-data-*`` distribution. + SSAPy includes: - Ability to define satellite parameters (area, mass, radiation and drag coefficients, etc.) diff --git a/pyproject.toml b/pyproject.toml index 72d4ba66..79fbd7aa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,10 +28,25 @@ dependencies = [ "ipython_genutils", "jplephem", "tqdm", - "llnl-ssapy-data>=0.2.0", ] [project.optional-dependencies] +data = [ + "ssatk-data-core==0.0.1", + "ssatk-data-gravity==0.0.1", + "ssatk-data-lunar==0.0.1", + "ssatk-data-lunar-gravity==0.0.1", +] + +all-data = [ + "ssatk-data-core==0.0.1", + "ssatk-data-gravity==0.0.1", + "ssatk-data-lunar==0.0.1", + "ssatk-data-lunar-gravity==0.0.1", + "ssatk-data-propulsion==0.0.1", + "ssatk-data-benchmarks==0.0.1", +] + dev = [ "pytest", "pytest-cov",