Add GMAT/Orekit accuracy comparison

[carlo98]

• Tickets addressed: docs: Add documentation comparing Basilisk's accuracy to other tools #1456
• Review: By file
• Merge strategy: Merge (no squash)

Description

Adds the documentation page requested in #1456, Support/User/accuracyComparison.rst, which compares Basilisk orbit propagation with GMAT and Orekit, together with the scripts and a guide to reproduce it. Seventeen cases add the perturbations one at a time (two-body, zonal and tesseral gravity up to degree and order 70, Sun and Moon, SRP, drag), so a disagreement can be traced to a single effect. They cover LEO, sun-synchronous, Molniya, GTO and GEO orbits, and faceted drag and SRP of a box spacecraft with a fixed and with a spinning attitude (compared with Orekit). All tools propagate in the ICRF axes (GCRF in Orekit, EarthICRF in GMAT, J2000 of the DE430 kernels in Basilisk), and the drag cases use a single-scale exponential atmosphere, which all three tools can model identically.

The comparison exposed a few Basilisk behaviours that are touched by this PR:

1. Planet orientation in gravityEffector (gravityEffector.cpp): the orientation of a SPICE-driven planet was advanced with DCM + DCM_dot * dt, which is not a rotation. It scales the evaluated field by a relative error of order (omega dt)^2 and distorts even a point-mass field. It is now advanced as a rotation about the planet angular velocity, with the new planetSpin() and advanceDcm() of stateExtrapolation.h, which are shared with the environment modules. gravityEffector derives the planet angular velocity once per planet message (planetSpin(), cached in GravBodyData) and only applies the rotation (advanceDcm()) at each integrator stage. advanceDcmDot() advances J20002Pfix_dot consistently with the advanced orientation, both in extrapolatePlanetStateToEpoch() (so the angular velocity reconstructed by WindBase from the extrapolated planet state is unchanged) and in gravityEffector, which now stores the advanced J20002Pfix and J20002Pfix_dot of the same epoch. These two state properties are stored as [PN] and [PN_dot], as documented in gravityEffector.h; they were the transposed [NP] after the first update, while their initial value was [PN].

2. One-step lag of the environment modules (new architecture/utilities/stateExtrapolation.h; AtmosphereBase, WindBase, MagneticFieldBase, Eclipse, SolarFlux): these modules run before the spacecraft, so they read its state from the previous step. They can now use the position extrapolated to the middle of the interval the next spacecraft update integrates. The extrapolation is off by default and is enabled with setExtrapolateScStateToStepMidpoint(True).

   • The planet (or Sun) state is moved to the same epoch (applyPlanet(), extrapolatePlanetStateToEpoch()), so a translating planet does not bias the relative position. This matters for the Sun in SolarFlux and for planet-fixed frames in AtmosphereBase.
   • The extrapolation is applied only if the spacecraft message was written at the previous module update and two successive spacecraft write times were observed with an interval equal to the module update interval (the first updates are not extrapolated). ScStateExtrapolation::prepare() (with a reusable write-time buffer, writeTimesBuffer()) makes one decision per module update for all spacecraft and the shared planets: if there is no spacecraft, or any spacecraft message was not written at the previous module update (stale, or a spacecraft that runs slower or faster than the module, as seen from its write interval), nothing is extrapolated and a warning is logged once (mismatchDetected()). This keeps every spacecraft and the planets at one epoch. A task period mismatch is not always detectable from the write times: a message written at the current module update (for example by a faster spacecraft that runs before the module) is used as written and gives no warning, but a message that is not the output of the previous module update is never extrapolated.
   • The requirement is on the task period, not on the integrator: the internal substeps of a variable-step integrator do not matter.
   • WindBase reads and extrapolates the planet state only when both the spacecraft and the planet messages were written.

3. MSIS altitude and solar time (AtmosphereBase, MsisAtmosphere), both opt-in so existing results are unchanged: setPlanetPolarRadius() gives the geodetic altitude (also used by ExponentialAtmosphere in the comparison); setUseApparentSolarTime() adds the equation of time to the local solar time (epochs before 1970 raise a BSK_ERROR).

4. Time zone dependence of the environment modules (new architecture/utilities/utcTime.h; AtmosphereBase, MsisAtmosphere, WindBase, MagneticFieldBase, MagneticFieldWMM): the UTC epoch was normalized with mktime, which uses the time zone of the computer. On a computer in a time zone with daylight saving time, the local solar time of MsisAtmosphere was one hour wrong once a simulation crossed a transition. The epoch is now normalized as UTC with timegm (_mkgmtime on Windows). Simulations that did not cross a transition are unchanged.

5. Solar flux of FacetSRPDynamicEffector (facetSRPDynamicEffector.cpp): the module used 1368 W/m^2 at 1 AU while RadiationPressure uses SOLAR_FLUX_EARTH (1361 W/m^2). Both now use SOLAR_FLUX_EARTH (1361 W/m^2), and the module takes the astronomical unit and the speed of light from astroConstants.h (AU2M, SPEED_LIGHT) instead of local constants. The two Python tests that recompute the force use the same values.

6. MSIS 3-hour Ap history (MsisAtmosphere), opt-in so existing results are unchanged: setUseApHistory() sets switch 9 of NRLMSISE-00 to -1, so the model uses the 3-hour Ap array that the module already builds from messages 1 to 20 instead of the daily Ap of message 0. The default is the daily Ap, as before.

7. PCPF2LLA() at a pole: the altitude was wrong for a position exactly at a pole of an oblate planet, where the equatorial component is zero.

test_radiationPressureIntegratedTest failed after the gravity-orientation fix because its stored truth contained the old orientation error. The truth is replaced with an independent scipy integration.

Notes for reviewers:

• The drag comparison uses an exponential atmosphere, which all three tools can model with the same two parameters. The MSIS options (items 3, 4 and 6) are opt-in features of MsisAtmosphere with their own unit tests.
• The GMAT and Orekit ephemerides are not shipped. The scripts that generate them and a "Reproducing the Results" section on the page are included; the tables and figures on the page come from a local run with GMAT R2026a and Orekit 13.1.
• The exponential atmosphere is defined by its native parameters (density at zero altitude and scale height), so Basilisk, GMAT and Orekit receive the same numbers without a conversion formula. Before any propagation with drag, compare_with_basilisk.py compares the altitude and density that each tool computes at nine Earth-fixed probe points (density_probe in cases.json, files <tool>_density_probe.csv written by the generators, with a manifest entry) against Basilisk, and stops if the altitude definition or the density differ beyond the tolerances there.
• Each generated reference writes a manifest with the case configuration, frame, tool versions, variant and data checksums. compare_with_basilisk.py validates it before running Basilisk and rejects a stale or mismatched reference with a message naming the field.
• The comparison lives under benchmarks/ and is not an automated test. The rotating-Earth cases use the high-precision earth_000101_260711_260415.bpc and earth_assoc_itrf93.tf kernels, which are not in the support-data registry or this repository; the script looks for them with --kernel-dir, in benchmarks/accuracyComparison/data/spice, or in the support-data cache.
• The Orekit reference uses a spherical occulting body for the eclipse and a 10 s (SRP) or 2 s (box) maximum step. With its defaults it differs from Basilisk and GMAT by 24 to 178 m. The page explains both effects, and --oblate-shadow and --max-step select the Orekit default (such references get a variant label and are rejected by a default comparison unless --orekit-variant is given).
• Other environment modules that only produce measurements (ground/strip location, sensors, albedo/earth radiation, Denton flux, charging) are unchanged.

Verification

• New and re-baselined unit tests: gravity of a rotating planet (orthonormal orientation, and stored orientation and rate of the same epoch), state extrapolation (C++ and Python step-lag tests for the exponential atmosphere, wind, magnetic field, eclipse and solar flux), utcTime, polar radius, MSIS apparent solar time and Ap history.
• Tests of the benchmark manifest and density probe validation (benchmarks/tests/test_accuracy_comparison_manifest.py), and two-spacecraft tests of the all-or-nothing extrapolation (C++ and Python), and tests of messages written at the module update by a faster or slower spacecraft (C++ and Python).
• Results of the 17 comparison cases over 30 days are on the documentation page. compare_with_basilisk.py enables the step extrapolation on the atmosphere, the wind and the eclipse modules.
• Constants in the new and modified unit tests are consistent with astroConstants.h.

Documentation

• New Support/User/accuracyComparison.rst (linked from Support/User.rst and the validation bullet in index.rst) with static tables, SVG figures, and a "Reproducing the Results" guide.
• Updated module pages: gravityEffector, atmosphereBase, windBase, magneticFieldBase, eclipse, solarFlux, msisAtmosphere, facetSRPDynamicEffector.
• Release-note snippet.
• bskKnownIssues.rst entries.

Future work

1. The atmosphere density is evaluated once per integration step, so the drag effector sees a density that is constant over the step. Option 1: let the drag effector query the atmosphere model at each integrator stage, or interpolate the density over the step. Querying at each stage would also remove the same-task-period requirement of the state extrapolation in stateExtrapolation.h. Option 2: handle a spacecraft period that differs from the module period instead of disabling the extrapolation. ScStateExtrapolation already observes the spacecraft period from the message write times, so the position could be extrapolated to the middle of the next spacecraft integration interval.
2. Apply the same handling to the modules left out here (planetRadiationBase/albedo/earth radiation, dentonFluxModel, charging).
3. Decide on defaults: the ellipsoidal altitude and the apparent solar time are opt-in; consider defaulting them after a deprecation cycle, and compute the solar time from the Sun position directly.
4. Add an NRLMSISE-00 comparison against GMAT and Orekit with a common implementation of the inputs (solar time, altitude, Ap), and with observed space weather.
5. Add the high-precision Earth orientation kernels (ITRF93) to the support-data registry so the comparison runs without manual setup.
6. Optional CI smoke test for compare_with_basilisk.py (only the manifest and probe validation are tested now).

[carlo98]

@ReeceHumphreys @suhaslord, I Started working on a comparison, if you want to take a look. It is leading to several changes though, it could be better to split them in several PRs, what do you think @schaubh ?

On another note, It would be great if someone with experience on GMAT and Orekit could review their setups.

Leaving this in draft as it could still change significantly and probably isn't at the top of the priorities.

[schaubh]

Requesting changes to address the four inline findings below: consistent epochs for midpoint extrapolation, aligned inertial frames in the comparison, complete rate-mismatch detection, and reference-file provenance. Reviewed commit 7208276.

[schaubh]

[P2] Midpoint extrapolation mixes spacecraft and planet epochs

The new code advances the spacecraft to the interval midpoint, then updateRelativePos() subtracts the planet's unadjusted position. This is incorrect when the planet translates in the simulation frame. I reproduced it with a spacecraft and planet moving together at 30 km/s: a 10-second interval introduced a 150 km altitude error despite constant physical separation. Planet orientation and time-dependent environment inputs also remain at the current epoch.

Please evaluate the relevant inputs at one consistent epoch, or explicitly restrict the supported configuration, and add a regression test with a translating planet. The same timing assumption should be checked in the other environment modules adopting this helper.

[schaubh]

[P2] The comparison does not align the inertial frames

Orekit initializes and exports states in EME2000. Basilisk uses SPICE's default J2000, which, for the selected DE430 and high-precision Earth kernels, means ICRF axes. The scripts assign identical numerical initial coordinates and subtract outputs without accounting for the rotation between these frames. That matters when interpreting meter-level differences.

This conclusion follows from the code and the documented frame definitions; I have not measured its contribution to each published result. Please use a consistent frame for initial conditions, environment data, attitudes, and output across the tools, then regenerate the affected comparisons.

References: NAIF's ICRF versus J2000 explanation and Orekit 13.1's EME2000 frame-bias implementation.

[schaubh]

[P2] The rate-mismatch warning misses a common mismatch

This condition only detects messages older than the previous environment update. When the spacecraft runs faster than the environment, its message is newer, so the warning never fires. I reproduced this with spacecraft updates every second and environment updates every ten seconds: at environment times 10, 20, and 30 s, spacecraft messages written at 9, 19, and 29 s are extrapolated without a warning.

The option therefore continues extrapolating silently outside its documented scheduling assumptions. Please add coverage for both directions of rate mismatch and either disable unsupported extrapolation or report it reliably.

[schaubh]

[P2] Reference files cannot establish which experiment they represent

The validation checks shape, length, and elapsed sample times, but the reference CSVs carry no epoch, configuration hash, frame, tool version, or data provenance. For example, generating Orekit references with --oblate-shadow overwrites the same filenames, and a later default comparison accepts them without identifying the changed physics. Changing an epoch in cases.json is similarly undetectable because elapsed sample times still match.

Please add a reference manifest that records the effective case configuration, generator options, frame, tool versions, and external-data identifiers/checksums, and validate it before running Basilisk. Intentional alternative configurations should be identified explicitly in the results rather than silently appearing as the default reference.

[schaubh]

[P2] Distinguish task periods from adaptive integrator substeps

This warning conflates task periods with adaptive integrator substeps. It says variable-step integrators invalidate extrapolation, yet the comparison itself uses adaptive RKF78. Internal adaptive substeps still cover the requested spacecraft update interval.

Please explain the scheduling requirement precisely and record Basilisk's integration tolerances in the accuracy-comparison documentation.

[schaubh]

Posted some feedback. Thanks for taking a stab at this.

I would not include the "Run Time" section of the benchmark comparison as it really is not apples-to-apples comparison. The tools have different features as you point out. If BSK has a 1Hz task rate, we force updates every 1s even if larger steps would be possible from an orbital mechanics point of view.

I would add a note at the top that these comparison are meant as a rough comparison. These are not auto-generated, so this page will become out of date with new releases of all tools!

You try hard to compare with MSIS drag included. There are a lot of factors here how this model is setup. Are you sure this is apples-and-apples at this point? Would it make more sense to just include a basic exponential density model to avoid these issues?

[carlo98]

Thank you! Yes, you are right. It would be great to compare everything, but it is better to start with the basics. I changed to exponential density and worked on your comments. I'll push a few changes in a bit