Spherical harmonic gravity #8

New Issue

2026-09-30T11:08:03-04:00

connor commented

2026-09-30 11:08:03 -04:00

Arbitrary degree and order, with J2 to J4 as the first case to get right. Validated against closed-form secular nodal regression and apsidal rotation rates.

M1.

Arbitrary degree and order, with J2 to J4 as the first case to get right. Validated against closed-form secular nodal regression and apsidal rotation rates. M1.

connor added this to the M1 milestone 2026-09-30 11:24:56 -04:00

connor added the size/L label 2026-09-30 11:24:56 -04:00

connor referenced this issue

2026-09-30 16:10:11 -04:00

Planetary bodies instead of bare mu #21

connor referenced this issue

2026-10-01 15:32:20 -04:00

Orekit fixture generation #16

connor commented

2026-10-01 16:14:28 -04:00

Done. src/gravity holds the field, src/dynamics/models/gravity.rs wraps it as a model, and Dynamics::Gravity { field, degree, order } is the spec.

Gottlieb's form of Pines' method: the expansion is written in the direction cosines of the position, so nothing divides by cos(latitude) and the poles are ordinary points. Arbitrary degree and order up to the coefficients, which are EGM96 to 360, committed as Rust source with their provenance and regenerable by scripts/generate-gravity-field.py.

Performance, field alone, no frame rotation:

degree	0	2	4	8	16	32	70	180	360
	7.0 ns	24.8 ns	59.9 ns	153 ns	458 ns	1.54 us	6.76 us	43.3 us	168 us

That is 2.6 ns per coefficient at both degree 70 and degree 360, so cost is linear in coefficient count with nothing superlinear hiding in it. Three things got it there: the normalization constants are a table built once per process rather than four square roots per coefficient per evaluation; only three degrees of the Legendre table are ever live, so the working set is 546 values at degree 360 rather than 131044; and the table is only as wide as the order, which makes a zonal degree-70 field 412 ns against the full field's 6.76 us, a 16x saving.

Validation. Nine tests on the field, three on the model.

Closed-form J2 nodal regression and apsidal rotation, agreeing to 3.8e-3 and 6.2e-3 with the tolerance declared at 1e-2. The gap is systematic rather than noise — it grows with arc length — and has two causes, both physics: the closed form is first order and describes mean elements where the fit is on osculating ones, and the field's symmetry axis is the true pole where the inclination is measured in J2000, which differ by about 0.15 degrees of accumulated precession at this epoch. The tight comparison belongs to #17.
An independently written potential, differentiated numerically, agreeing to better than 1e-8 on the perturbation across five positions and four truncations.
The poles, degree zero reproducing the point mass to 1e-15, degree one contributing nothing, truncation converging, and the autodiff path against central differences.

Two things worth carrying forward. The expansion replaces TwoBody rather than adding to it, since the point mass is its degree-zero term and composing both counts the monopole twice. And the field carries its own mu and reference radius: EGM96's radius is 6378136.3 m where Body::EARTH's is WGS84's 6378137.0 m, and those are different quantities that happen to look alike.

examples/full.rs exercises it end to end, and reproduces the 5.0 deg/day nodal regression that closed-form theory predicts for this orbit.

Done. `src/gravity` holds the field, `src/dynamics/models/gravity.rs` wraps it as a model, and `Dynamics::Gravity { field, degree, order }` is the spec. Gottlieb's form of Pines' method: the expansion is written in the direction cosines of the position, so nothing divides by `cos(latitude)` and the poles are ordinary points. Arbitrary degree and order up to the coefficients, which are EGM96 to 360, committed as Rust source with their provenance and regenerable by `scripts/generate-gravity-field.py`. **Performance**, field alone, no frame rotation: | degree | 0 | 2 | 4 | 8 | 16 | 32 | 70 | 180 | 360 | |---|---|---|---|---|---|---|---|---|---| | | 7.0 ns | 24.8 ns | 59.9 ns | 153 ns | 458 ns | 1.54 us | 6.76 us | 43.3 us | 168 us | That is 2.6 ns per coefficient at both degree 70 and degree 360, so cost is linear in coefficient count with nothing superlinear hiding in it. Three things got it there: the normalization constants are a table built once per process rather than four square roots per coefficient per evaluation; only three degrees of the Legendre table are ever live, so the working set is 546 values at degree 360 rather than 131044; and the table is only as wide as the order, which makes a zonal degree-70 field 412 ns against the full field's 6.76 us, a 16x saving. **Validation.** Nine tests on the field, three on the model. - Closed-form J2 nodal regression and apsidal rotation, agreeing to 3.8e-3 and 6.2e-3 with the tolerance declared at 1e-2. The gap is systematic rather than noise — it grows with arc length — and has two causes, both physics: the closed form is first order and describes mean elements where the fit is on osculating ones, and the field's symmetry axis is the true pole where the inclination is measured in J2000, which differ by about 0.15 degrees of accumulated precession at this epoch. The tight comparison belongs to #17. - An independently written potential, differentiated numerically, agreeing to better than 1e-8 on the perturbation across five positions and four truncations. - The poles, degree zero reproducing the point mass to 1e-15, degree one contributing nothing, truncation converging, and the autodiff path against central differences. **Two things worth carrying forward.** The expansion replaces `TwoBody` rather than adding to it, since the point mass is its degree-zero term and composing both counts the monopole twice. And the field carries its own `mu` and reference radius: EGM96's radius is 6378136.3 m where `Body::EARTH`'s is WGS84's 6378137.0 m, and those are different quantities that happen to look alike. `examples/full.rs` exercises it end to end, and reproduces the 5.0 deg/day nodal regression that closed-form theory predicts for this orbit.

connor closed this issue

2026-10-01 16:14:28 -04:00

Sign in to join this conversation.

1 Participants

Notifications

Due Date

No due date set.

Dependencies

No dependencies set.

Reference: connor/astrokit#8