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 with planetSpin() during each gravity evaluation and advances the orientation with advanceDcm() to that evaluation epoch. 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.
   • Midpoint extrapolation starts only after the observed spacecraft write interval matches the module interval and all spacecraft messages were written at the previous module update. During startup, or whenever midpoint extrapolation is disabled but all spacecraft messages share that previous-update epoch, spacecraft states remain as written and planet states are moved to their epoch. Otherwise, the states remain as written; unsupported mixed message epochs do not guarantee consistent relative geometry. A detected scheduling mismatch produces a one-time warning. No spacecraft messages means no extrapolation and no mismatch warning.
   • 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.
   • A planet message that only sets position and velocity leaves J20002Pfix at zero. applyPlanet() treats an all-zero orientation as the identity on every path, with or without extrapolation, so all the environment modules (AtmosphereBase, WindBase, MagneticFieldBase, Eclipse, SolarFlux) see a valid rotation. gravityEffector applies the same rule.
   • When extrapolation is enabled, time-dependent models (MSIS local solar time, WMM, wind) are evaluated at the extrapolated epoch rather than at the current time. With extrapolation off, they behave as before.

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 ellipsoidal-altitude, apparent-solar-time and Ap-history options in items 3 and 6 are opt-in. The UTC/time-zone correction in item 4 applies unconditionally.
• 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 for 17 comparison cases are on the documentation page: thirteen cases run for 30 days, while leo_egm70, leo_box_fixed, leo_box_spin and geo_box_spin run for seven days.
• 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 environment-module update and held between updates, including through the spacecraft integrator's internal stages. 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]

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

[schaubh]

Thanks, we are getting close now. I had two follow up comments and one new comment to consider.

[carlo98]

The latest commit addresses your comments and a few more issues. Forced push to rebase and solve a conflict in bskKnownIssues.rst.

[carlo98]

Removed caching in gravityEffector previously added by this branch as it would have led to rare time savings at the cost of additional operations for most of the updates. Moved the shared extrapolation logic in stateExtrapolation.h to avoid inconsistencies.

[schaubh]

[P3] Update the PR description after removing the gravity spin cache

Item 1 still says that gravityEffector derives planet angular velocity “once per planet message” and only applies the rotation at each integrator stage. The latest commit removes that cache: GravBodyData::computeGravityInertial() now calls planetSpin() during each gravity evaluation.

Please replace that sentence with: “gravityEffector derives the planet angular velocity with planetSpin() during each gravity evaluation and advances the orientation with advanceDcm() to that evaluation epoch.”

[schaubh]

Ok, I think we are there. Thanks for the contribution.