Repository navigation
Fix/bug review - #152
Merged
Merged
Fix/bug review#152
Conversation
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.
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.
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.
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.
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.
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).
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
MoonOrientation, MoonPosition, SunPosition and PlanetPosition opened ssapy/data/<file> 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 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.
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.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.