Unify propagator force model; cap gravity coefficient table#114
Merged
Conversation
Collapse the four near-identical force closures (RKV ydot, RODAS4 ydot_vec/jac_fn, Gauss-Jackson accel_fn) into a single inlined force_model() with a ForceEval mode selector. Each integrator closure is now a thin adapter over its own state layout. Because every call site passes a compile-time-constant mode, the partials/accel branches constant-fold away, reproducing the prior per-integrator specialization with no extra call or runtime dispatch. End states are bit-identical. Also gate atmospheric drag on norm_squared() against a precomputed squared limit instead of norm(), dropping a sqrt from every force eval. Cap the stored Earth-gravity coefficient table at degree 44 (MAX_COEFF_DIM) instead of the file's full resolution. The evaluator dispatches at degree <= 40 and the divisor tables are 44x44, so higher coefficients were never used; EGM96 (file degree 360) was holding a ~1 MB 361x361 DMatrix. Capping drops it to ~15 KB and shrinks the column stride 361 -> 44, keeping the strided S-coefficient reads in L1. Coefficients for n,m <= 43 are untouched, so results are bit-identical; ~5-15% faster accel/accel_and_partials on EGM96 plus removal of cache-miss timing variance. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…on-gravity-cap # Conflicts: # CHANGELOG.md
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
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.
Summary
Two internal refactor/perf changes to the high-precision propagator and Earth-gravity model. No behavior change — end states are bit-identical and the full Rust test suite (including the ESA SP3 GPS regression) passes unchanged.
1. Force model unified across all integrators
The RKV
ydot, RODAS4ydot_vec/jac_fn, and Gauss-Jacksonaccel_fneach carried a near-identical copy of the acceleration physics (gravity, Sun/Moon, solid tides, GR, SRP, drag, thrust) — ~200 lines duplicated four ways, where adding a force term meant editing every copy in lockstep.force_model()now holds the physics; each integrator closure is a thin adapter over its own state layout (Matrix<6, C>/Vector<f64, 6>/(r, v)).ForceEvalmode selector (Accel/AccelAndPartials/PartialsOnly) preserves each integrator's cost. Every call site passes a compile-time-constant mode andforce_modelis#[inline], so the branches constant-fold away — no extra call, no runtime dispatch. The RODAS4 Jacobian path skips the tide/GR/SRP/thrust terms (no partials), so it does not regress.norm_squared()against a precomputed squared limit instead ofnorm(), dropping a sqrt from every force eval.2. Earth gravity: coefficient table capped at the evaluation degree
The evaluator dispatches at degree ≤ 40 (divisor tables are 44×44), so higher coefficients were never used — yet the default EGM96 model (file degree 360) stored a 361×361
DMatrix(~1 MB).Gravity::parsenow caps the stored table at degree 44 (MAX_COEFF_DIM). Drops EGM96 to ~15 KB and shrinks the column stride 361 → 44, keeping the strided S-coefficient readscoeffs[(m-1, n)]resident in L1.accel/accel_and_partialson EGM96 and removes the cache-miss timing variance in the hot loop.Test plan
cargo testfull suite withSATKIT_DATA+ test vectors — zero failuresearthgravityunit tests 10/10cargo fmt --allcleanpip install -e .rebuild)🤖 Generated with Claude Code