connor/astrokit
1
0 Fork 0
Code Issues 20 Pull Requests Actions Packages Projects Releases Wiki Activity
Files
98f19b700a540a630b4a658355a770feb4f5dfd7
astrokit/scripts/generate-orekit-fixtures.sh
T

787 lines
40 KiB
Bash
Executable File
Raw Blame History

#!/bin/sh
# Generates tests/fixtures/orekit/*.json: reference trajectories from Orekit, for this
# crate's propagator to be compared against.
#
# Fetches Orekit and its Hipparchus dependencies from Maven Central at pinned digests,
# compiles a driver against them, runs it and throws the jars away. What is committed is
# the trajectories, which are numbers rather than anyone's code.
#
# Orekit is BSD-3-Clause, by CS GROUP and contributors. Hipparchus is Apache-2.0.
#
# Needs a JDK and network access.
#
# Usage: scripts/generate-orekit-fixtures.sh [--time]
set -e
OREKIT=13.1.8
HIPPARCHUS=4.0.3
MAVEN=https://repo1.maven.org/maven2
# Moves as Earth orientation is updated, so there is nothing fixed to pin. The digest of
# what was downloaded goes into every fixture this run writes.
DATA_URL=https://gitlab.orekit.org/orekit/orekit-data/-/archive/main/orekit-data-main.zip
# `--time` runs each scenario under the JVM and reports what it costs, instead of
# writing fixtures. Everything else about the run is identical, so the timings are of the
# very configurations the fixtures were generated from.
mode=fixtures
if [ "$1" = "--time" ]; then
mode=time
fi
root=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd)
out="$root/tests/fixtures/orekit"
work=$(mktemp -d)
trap 'rm -rf "$work"' EXIT
if ! command -v javac >/dev/null 2>&1; then
echo "a JDK is needed and javac is not on PATH" >&2
exit 1
fi
mkdir -p "$work/lib"
# The jars are immutable on Maven Central, so these digests are a lock rather than a
# record. A function rather than a loop over a list: a `while` fed by a pipe runs in a
# subshell, where the `exit` on a mismatch would not stop this script.
fetch() {
artifact=$2
name="$artifact-$3.jar"
curl -sS --fail --max-time 300 -L -o "$work/lib/$name" "$MAVEN/$1/$artifact/$3/$name"
got=$(sha256sum "$work/lib/$name" | cut -d' ' -f1)
if [ "$got" != "$4" ]; then
echo "digest mismatch for $name" >&2
echo " expected $4" >&2
echo " got $got" >&2
exit 1
fi
}
fetch org/orekit orekit "$OREKIT" \
5d035486d916a458d6a541866db99885a21d5bc0cbdb3785bb8f42ee8a70ddf4
fetch org/hipparchus hipparchus-core "$HIPPARCHUS" \
5af54dc3c21b8d85e1264157ee91d93fa67e68c825cfb38834dc2bbb9b0f6583
fetch org/hipparchus hipparchus-geometry "$HIPPARCHUS" \
19594d5106ec6a46951e859e2924bfcade02a834c06c56bd3d2130a74f0fbb7e
fetch org/hipparchus hipparchus-ode "$HIPPARCHUS" \
75f0f0b3fdcd99d13e9c394b7ae809cce5a9c22f1c14b1b90a5109ed951467fe
fetch org/hipparchus hipparchus-fitting "$HIPPARCHUS" \
3712c3b4f43d9c12bf7b7aad187764f4b3c09cc962504899d8951eb310b1e4be
fetch org/hipparchus hipparchus-optim "$HIPPARCHUS" \
4bc0dbb5c6d5fefda5db9d92141c6cbb807174a1d491a85df46bcc57a67b8992
fetch org/hipparchus hipparchus-filtering "$HIPPARCHUS" \
279eece43411f372bd8340355cb89101def0bf9b1b6824231b76181d1ec63a81
fetch org/hipparchus hipparchus-stat "$HIPPARCHUS" \
8b066e3060dc119efcce241fb0720a3b9e59c026d0f5fe9e70abfd5e22d3817b
echo "fetched orekit $OREKIT and hipparchus $HIPPARCHUS, digests verified"
curl -sS --fail --max-time 600 -L -o "$work/orekit-data.zip" "$DATA_URL"
data_digest=$(sha256sum "$work/orekit-data.zip" | cut -d' ' -f1)
unzip -q -o "$work/orekit-data.zip" -d "$work"
data="$work/orekit-data-main"
echo "fetched orekit-data, sha256 $data_digest"
# The gravity comparison is only worth anything if both sides read the same numbers, so
# Orekit is given the very file this crate's own coefficients are generated from, pinned
# at the same digest as scripts/generate-gravity-field.py. GMAT wraps the standard
# egm96_to360.ascii in its own record format, and Orekit's reader wants it plain.
EGM96_URL=https://raw.githubusercontent.com/nasa/GMAT/main/application/data/gravity/earth/EGM96.cof
EGM96_SHA256=4735c9301e20ee8db09b6d0882a24671b327cd35ec5740ea497c8eeb43a39807
curl -sS --fail --max-time 300 -L -o "$work/EGM96.cof" "$EGM96_URL"
got=$(sha256sum "$work/EGM96.cof" | cut -d' ' -f1)
if [ "$got" != "$EGM96_SHA256" ]; then
echo "digest mismatch for EGM96.cof" >&2
echo " expected $EGM96_SHA256" >&2
echo " got $got" >&2
exit 1
fi
mkdir -p "$work/gravity"
python3 - "$work/EGM96.cof" "$work/gravity/egm96_to360.ascii" <<'PY'
import re
import sys
# Fixed-width RECOEF records, degree in columns 6-10 and order in 11-13, the same fields
# scripts/generate-gravity-field.py reads. The cosine and sine terms are not separated
# when the sine is negative, so they are matched rather than split.
rows = []
for line in open(sys.argv[1]):
if not line.startswith("RECOEF"):
continue
degree, order = int(line[6:11]), int(line[11:14])
values = re.findall(r"[-+]?\d\.\d+E[-+]\d+", line[14:])
cosine = float(values[0])
sine = float(values[1]) if len(values) > 1 else 0.0
# Explicit separators: a negative value's sign would otherwise run into the column
# before it and the free-format reader would see one token where there are two.
rows.append(f"{degree:6d} {order:5d} {cosine: .12E} {sine: .12E} 0.0E+00 0.0E+00")
open(sys.argv[2], "w").write("\n".join(rows) + "\n")
print(f"converted {len(rows)} EGM96 coefficients for Orekit")
PY
cat > "$work/Generate.java" <<'JAVA'
import java.io.File;
import java.io.PrintWriter;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.orekit.bodies.CelestialBody;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.NewtonianAttraction;
import org.orekit.forces.gravity.ThirdBodyAttraction;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.models.earth.atmosphere.HarrisPriester;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.NRLMSISE00InputParameters;
import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherData;
import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
import org.orekit.forces.radiation.SolarRadiationPressure;
import org.orekit.forces.gravity.potential.EGMFormatReader;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.ITRFVersion;
import org.orekit.orbits.CartesianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.ExtendedPositionProvider;
import org.orekit.utils.PVCoordinates;
/** Writes reference arcs in this crate's fixture format. */
public class Generate {
/** CSSI's space weather, with the 81-day flux average recomputed rather than read.
*
* CSSI publishes that average rounded to a tenth, and this crate computes it from the
* daily values, so the two differ by about 0.028% of the density for no reason other
* than the rounding. The daily fluxes themselves are published in full, so the average
* is recovered by averaging them — the same 81 days, centered, that this crate uses.
*
* `getDailyFlux(d)` is the flux of the day before `d`, which is what the model wants,
* so the window is shifted by one day to cover the same calendar days.
*/
static class UnroundedWeather implements NRLMSISE00InputParameters {
private final CssiSpaceWeatherData inner;
/** One average per day, because that is how often it changes.
*
* Without this the eighty-one lookups run on every density evaluation, which is
* every stage of every step: measured at 279 ms of a 641 ms arc, so the timing
* would have been reporting this wrapper rather than Orekit.
*/
private final java.util.Map<Long, Double> cache = new java.util.HashMap<>();
UnroundedWeather(CssiSpaceWeatherData inner) {
this.inner = inner;
}
public AbsoluteDate getMinDate() {
return inner.getMinDate().shiftedBy(41 * 86400.0);
}
public AbsoluteDate getMaxDate() {
return inner.getMaxDate().shiftedBy(-41 * 86400.0);
}
public double getDailyFlux(AbsoluteDate date) {
return inner.getDailyFlux(date);
}
public double getAverageFlux(AbsoluteDate date) {
long day = (long) Math.floor(date.durationFrom(AbsoluteDate.MODIFIED_JULIAN_EPOCH)
/ 86400.0);
Double known = cache.get(day);
if (known != null) {
return known;
}
double sum = 0.0;
for (int offset = -39; offset <= 41; offset++) {
sum += inner.getDailyFlux(date.shiftedBy(offset * 86400.0));
}
double average = sum / 81.0;
cache.put(day, average);
return average;
}
public double[] getAp(AbsoluteDate date) {
return inner.getAp(date);
}
}
/** EGM96's, which is what this crate's Earth carries. */
static final double MU = 3.986004415e14;
/** Position accuracy the integrator is asked for, well inside what is compared.
*
* A nanometer, which is absurd for a result and is the point: the comparison has to
* measure the models and not either integrator. One micrometer was tried first and is
* enough for the smooth cases but not for radiation pressure, where the arc passes
* through the shadow boundary thirty-odd times a day — it left that fixture 0.4 m out,
* which is the size of the disagreement being looked for.
*/
static final double TARGET_METERS = 1.0e-9;
/** Matching this crate's constants module rather than Orekit's own. */
static final double IRRADIANCE_AT_1AU = 1361.0;
static final double LIGHT_SPEED = 299792458.0;
static final double ASTRONOMICAL_UNIT = 1.495978707e11;
static final double EARTH_RADIUS = 6378137.0;
/** The cannonball: mass, cross-section and radiation coefficient. */
static final double MASS = 500.0;
static final double AREA = 5.0;
static final double CR = 1.3;
static final double CD = 2.2;
static final double EARTH_FLATTENING = 1.0 / 298.257223563;
static final double BULGE_EXPONENT = 2.0;
static final double ARC_SECONDS = 86400.0;
static final double STEP_SECONDS = 60.0;
static String version;
static String dataDigest;
static String generated;
public static void main(String[] args) throws Exception {
File gravity = new File(args[0]);
File data = new File(args[1]);
File out = new File(args[2]);
version = args[3];
dataDigest = args[4];
generated = args[5];
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(data));
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(gravity));
out.mkdirs();
if (args.length > 6 && args[6].equals("time")) {
time();
return;
}
twoBody(out);
gravity(out, 8);
gravity(out, 20);
thirdBody(out);
radiation(out);
drag(out);
msisDrag(out);
}
/** EGM96 at a given truncation, which is the first fixture with a model to match. */
static void gravity(File out, int degree) throws Exception {
GravityFieldFactory.clearPotentialCoefficientsReaders();
GravityFieldFactory.addPotentialCoefficientsReader(
new EGMFormatReader("egm96_to360\\.ascii", false));
NormalizedSphericalHarmonicsProvider provider =
GravityFieldFactory.getNormalizedProvider(degree, degree);
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame inertial = FramesFactory.getGCRF();
// ITRF93 specifically, which is the realization this crate reads from its SPICE
// binary PCK. Orekit defaults to the latest, and the difference is real.
Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_1993, IERSConventions.IERS_2010, true);
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.978e6, 0.0, 0.0), new Vector3D(0.0, 7.56e3, 0.0)),
inertial, epoch, provider.getMu());
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
// Holmes-Featherstone carries the central term, so there is no separate point
// mass beside it; adding one would count the monopole twice.
propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
propagator.setInitialState(new SpacecraftState(orbit));
write(new File(out, String.format(Locale.ROOT, "gravity-%d-low-earth.json", degree)),
String.format(Locale.ROOT, "gravity-%d-low-earth", degree),
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] { String.format(Locale.ROOT,
"HolmesFeatherstoneAttractionModel(EGM96 %dx%d, ae = %.1f, mu = %.9e)",
degree, degree, provider.getAe(), provider.getMu()) },
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
"Coefficients read from the same EGM96 file this crate generates its own from.",
"Inertial frame is GCRF and body-fixed is ITRF93 under IERS 2010 with EOP"
+ " from orekit-data, chosen to match what this crate reads from"
+ " ANISE. The two rotations still differ by about 4 mas, which is"
+ " the floor this comparison sits on.",
"No third body, no atmosphere, no radiation pressure.",
},
states(propagator, epoch));
}
/** The Sun and the Moon, over a point mass, from the DE-440 ephemeris.
*
* No field: the point is to isolate the third-body term, which is a difference of
* two nearly equal attractions and so the one most easily got wrong.
*/
static void thirdBody(File out) throws Exception {
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame inertial = FramesFactory.getGCRF();
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.978e6, 0.0, 0.0), new Vector3D(0.0, 7.56e3, 0.0)),
inertial, epoch, MU);
CelestialBody sun = CelestialBodyFactory.getSun();
CelestialBody moon = CelestialBodyFactory.getMoon();
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
propagator.addForceModel(new NewtonianAttraction(MU));
propagator.addForceModel(new ThirdBodyAttraction(moon));
propagator.addForceModel(new ThirdBodyAttraction(sun));
propagator.setInitialState(new SpacecraftState(orbit));
write(new File(out, "third-body-low-earth.json"), "third-body-low-earth",
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] {
String.format(Locale.ROOT, "NewtonianAttraction(mu = %.9e)", MU),
String.format(Locale.ROOT, "ThirdBodyAttraction(Moon, mu = %.12e)",
moon.getGM()),
String.format(Locale.ROOT, "ThirdBodyAttraction(Sun, mu = %.12e)",
sun.getGM()),
},
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
"Ephemeris is DE-440, from orekit-data's lnxp1990.440; this crate reads"
+ " the same model from de440s.bsp through ANISE.",
"The gravitational parameters above are the ephemeris header's, which are"
+ " not necessarily the ones this crate carries for those bodies.",
"No field, no atmosphere, no radiation pressure.",
},
states(propagator, epoch));
}
/** Radiation pressure on a cannonball, over a point mass.
*
* Cannonball rather than a faceted body: a faceted one would need the two tools to
* agree on the attitude as well as the geometry, which is a larger matching problem
* than the force itself.
*/
static void radiation(File out) throws Exception {
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame inertial = FramesFactory.getGCRF();
Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_1993, IERSConventions.IERS_2010, true);
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.978e6, 0.0, 0.0), new Vector3D(0.0, 7.56e3, 0.0)),
inertial, epoch, MU);
// Orekit's default is 4.56e-6 N/m^2, which is 1367 W/m^2; this crate uses 1361.
// Set explicitly rather than left to either tool's idea of the solar constant.
double pressure = IRRADIANCE_AT_1AU / LIGHT_SPEED;
// Spherical, because this crate's eclipse geometry treats the occulting body as a
// sphere of its equatorial radius. Orekit's Earth is oblate unless told otherwise.
OneAxisEllipsoid occulter = new OneAxisEllipsoid(EARTH_RADIUS, 0.0, itrf);
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
propagator.addForceModel(new NewtonianAttraction(MU));
// Distance first, then pressure. Passing them the other way round is accepted
// silently and makes the force about 1e-16 of what it should be.
propagator.addForceModel(new SolarRadiationPressure(ASTRONOMICAL_UNIT, pressure,
(ExtendedPositionProvider) CelestialBodyFactory.getSun(), occulter,
new IsotropicRadiationSingleCoefficient(AREA, CR)));
propagator.setInitialState(new SpacecraftState(orbit).withMass(MASS));
write(new File(out, "radiation-low-earth.json"), "radiation-low-earth",
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] {
String.format(Locale.ROOT, "NewtonianAttraction(mu = %.9e)", MU),
String.format(Locale.ROOT,
"SolarRadiationPressure(pRef = %.9e N/m^2 at %.9e m, "
+ "IsotropicRadiationSingleCoefficient(area = %.1f, cr = %.1f))",
pressure, ASTRONOMICAL_UNIT, AREA, CR),
},
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
String.format(Locale.ROOT, "Spacecraft mass %.1f kg.", MASS),
"Reference pressure set from 1361 W/m^2, which is what this crate uses;"
+ " Orekit's own default is 1367.",
"The occulting Earth is a sphere of 6378137 m, matching this crate's"
+ " eclipse geometry; Orekit would otherwise use an oblate one.",
"Orekit's solar radius is 6.95508e8 m against this crate's 6.957e8, which"
+ " differs only across the penumbra.",
"No field, no third body, no atmosphere.",
},
states(propagator, epoch));
}
/** Drag on a cannonball through Harris-Priester, over a point mass.
*
* Lower than the other scenarios, at 400 km rather than 600, because drag at 600 km
* is too weak over a day to measure a disagreement against.
*/
static void drag(File out) throws Exception {
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame inertial = FramesFactory.getGCRF();
Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_1993, IERSConventions.IERS_2010, true);
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.778e6, 0.0, 0.0), new Vector3D(0.0, 7.668e3, 0.0)),
inertial, epoch, MU);
// The real ellipsoid here, not the sphere the eclipse wanted: altitude into a
// density model is geodetic, and this crate computes it against WGS84.
OneAxisEllipsoid earth = new OneAxisEllipsoid(EARTH_RADIUS, EARTH_FLATTENING, itrf);
// Exponent 2, this crate's low-inclination value. Orekit's own default is 4.
HarrisPriester atmosphere = new HarrisPriester(
(ExtendedPositionProvider) CelestialBodyFactory.getSun(), earth, BULGE_EXPONENT);
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
propagator.addForceModel(new NewtonianAttraction(MU));
propagator.addForceModel(new DragForce(atmosphere, new IsotropicDrag(AREA, CD)));
propagator.setInitialState(new SpacecraftState(orbit).withMass(MASS));
write(new File(out, "drag-low-earth.json"), "drag-low-earth",
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] {
String.format(Locale.ROOT, "NewtonianAttraction(mu = %.9e)", MU),
String.format(Locale.ROOT,
"DragForce(HarrisPriester(n = %.1f), IsotropicDrag(area = %.1f, cd = %.1f))",
BULGE_EXPONENT, AREA, CD),
},
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
String.format(Locale.ROOT, "Spacecraft mass %.1f kg, starting at 400 km.", MASS),
"Harris-Priester with the same tabulated profiles and the same 30 degree"
+ " bulge lag, but Orekit interpolates the table exponentially while"
+ " this crate fits a monotone cubic to the logarithm of the density."
+ " That is a deliberate difference, worth about 0.2% of the density"
+ " in this band, and the disagreement here is what it costs.",
"Altitude is geodetic against WGS84 on both sides.",
"No field, no third body, no radiation pressure.",
},
states(propagator, epoch));
}
/** Drag through NRLMSISE-00, which both tools implement from the same reference.
*
* Unlike Harris-Priester, where this crate deliberately interpolates the table
* differently, there is no intended difference here — so what the comparison measures
* is the two implementations and the space weather they are given.
*/
static void msisDrag(File out) throws Exception {
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame inertial = FramesFactory.getGCRF();
Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_1993, IERSConventions.IERS_2010, true);
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.778e6, 0.0, 0.0), new Vector3D(0.0, 7.668e3, 0.0)),
inertial, epoch, MU);
OneAxisEllipsoid earth = new OneAxisEllipsoid(EARTH_RADIUS, EARTH_FLATTENING, itrf);
CssiSpaceWeatherData weather =
new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES);
// Two knobs, both of which have to be moved off Orekit's defaults to match.
//
// MEAN local solar time: Orekit defaults to APPARENT, which carries the equation
// of time, and this crate derives the time from UT and longitude as the reference
// does. Worth up to 0.6% of the density.
//
// Switch 9 at -1: the reference's default is to take a single daily Ap, and this
// crate passes the seven-element history the model can also accept. Worth 1.1%,
// and the larger of the two.
NRLMSISE00 atmosphere = new NRLMSISE00(new UnroundedWeather(weather),
(ExtendedPositionProvider) CelestialBodyFactory.getSun(), earth)
.withLocalSolarTimeMode(NRLMSISE00.LocalSolarTimeMode.MEAN)
.withSwitch(9, -1);
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
propagator.addForceModel(new NewtonianAttraction(MU));
propagator.addForceModel(new DragForce(atmosphere, new IsotropicDrag(AREA, CD)));
propagator.setInitialState(new SpacecraftState(orbit).withMass(MASS));
write(new File(out, "msis-drag-low-earth.json"), "msis-drag-low-earth",
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] {
String.format(Locale.ROOT, "NewtonianAttraction(mu = %.9e)", MU),
String.format(Locale.ROOT,
"DragForce(NRLMSISE00(CSSI space weather), "
+ "IsotropicDrag(area = %.1f, cd = %.1f))", AREA, CD),
},
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
String.format(Locale.ROOT, "Spacecraft mass %.1f kg, starting at 400 km.", MASS),
"Mean local solar time and switch 9 at -1, neither of which is Orekit's"
+ " default: it uses apparent solar time, worth 0.6% here, and the"
+ " reference's single daily Ap rather than the seven-element"
+ " history this crate passes, worth 1.1%.",
"Space weather is CSSI's, out of orekit-data; this crate reads GFZ's. For"
+ " this epoch the daily flux and every three-hourly ap agree"
+ " exactly, and the 81-day averages differ only because CSSI"
+ " publishes it rounded to a tenth. That rounding is recovered"
+ " here by averaging CSSI's own daily fluxes over the same"
+ " centered 81 days, so the two models see the same number.",
"Altitude is geodetic against WGS84 on both sides.",
"No field, no third body, no radiation pressure.",
},
states(propagator, epoch));
}
/** Times each scenario, for comparison against this crate's own propagator.
*
* The JVM is the awkward part: a cold run measures the bytecode interpreter rather
* than the compiled code, and the difference is often an order of magnitude. Each
* scenario is therefore run WARMUP times and discarded before anything is recorded,
* which is enough for the hot paths to reach C2, and the figure reported is the
* minimum rather than the mean — a minimum is the run that was least disturbed by
* garbage collection and rescheduling, which is what a like-for-like comparison
* against a Rust binary wants.
*/
static void time() throws Exception {
String[] names = { "two-body", "gravity-8", "gravity-20", "third-body", "radiation",
"drag-harris-priester", "drag-msis", "drag-msis-plain" };
for (String name : names) {
// Build it fresh each run: a propagator carries its state, and reusing one
// would time a different thing the second time round.
for (int i = 0; i < WARMUP; i++) {
scenario(name).propagate(EPOCH().shiftedBy(ARC_SECONDS));
}
long best = Long.MAX_VALUE;
long total = 0;
for (int i = 0; i < RUNS; i++) {
NumericalPropagator p = scenario(name);
long started = System.nanoTime();
p.propagate(EPOCH().shiftedBy(ARC_SECONDS));
long took = System.nanoTime() - started;
best = Math.min(best, took);
total += took;
}
// How much work it did, so a ratio of times can be told apart from a ratio
// of tolerances. Two tools asked for different accuracy take different
// numbers of steps and their times are then not comparable.
NumericalPropagator counted = scenario(name);
final int[] steps = { 0 };
counted.setStepHandler(interpolator -> steps[0]++);
counted.propagate(EPOCH().shiftedBy(ARC_SECONDS));
System.out.printf(Locale.ROOT,
"TIMING %-22s min %8.2f ms mean %8.2f ms steps %5d%n",
name, best / 1.0e6, total / 1.0e6 / RUNS, steps[0]);
}
}
static final int WARMUP = 10;
static final int RUNS = 30;
static AbsoluteDate EPOCH() {
return new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, TimeScalesFactory.getUTC());
}
/** One scenario's propagator, configured exactly as its fixture was. */
static NumericalPropagator scenario(String name) throws Exception {
AbsoluteDate epoch = EPOCH();
Frame inertial = FramesFactory.getGCRF();
Frame itrf = FramesFactory.getITRF(ITRFVersion.ITRF_1993, IERSConventions.IERS_2010, true);
boolean low = name.startsWith("drag");
Vector3D r = low ? new Vector3D(6.778e6, 0.0, 0.0) : new Vector3D(6.978e6, 0.0, 0.0);
Vector3D v = low ? new Vector3D(0.0, 7.668e3, 0.0) : new Vector3D(0.0, 7.56e3, 0.0);
Orbit orbit = new CartesianOrbit(new PVCoordinates(r, v), inertial, epoch, MU);
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
NumericalPropagator p = new NumericalPropagator(
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]));
p.setOrbitType(OrbitType.CARTESIAN);
if (name.startsWith("gravity")) {
int degree = Integer.parseInt(name.substring(name.indexOf('-') + 1));
GravityFieldFactory.clearPotentialCoefficientsReaders();
GravityFieldFactory.addPotentialCoefficientsReader(
new EGMFormatReader("egm96_to360\\.ascii", false));
p.addForceModel(new HolmesFeatherstoneAttractionModel(itrf,
GravityFieldFactory.getNormalizedProvider(degree, degree)));
} else {
p.addForceModel(new NewtonianAttraction(MU));
}
if (name.equals("third-body")) {
p.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
p.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
}
if (name.equals("radiation")) {
p.addForceModel(new SolarRadiationPressure(ASTRONOMICAL_UNIT,
IRRADIANCE_AT_1AU / LIGHT_SPEED,
(ExtendedPositionProvider) CelestialBodyFactory.getSun(),
new OneAxisEllipsoid(EARTH_RADIUS, 0.0, itrf),
new IsotropicRadiationSingleCoefficient(AREA, CR)));
}
if (name.equals("drag-harris-priester")) {
p.addForceModel(new DragForce(new HarrisPriester(
(ExtendedPositionProvider) CelestialBodyFactory.getSun(),
new OneAxisEllipsoid(EARTH_RADIUS, EARTH_FLATTENING, itrf), BULGE_EXPONENT),
new IsotropicDrag(AREA, CD)));
}
if (name.equals("drag-msis-plain")) {
// The same thing with CSSI's published average, to show what the recomputing
// wrapper costs rather than assuming it is free.
p.addForceModel(new DragForce(new NRLMSISE00(
new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES),
(ExtendedPositionProvider) CelestialBodyFactory.getSun(),
new OneAxisEllipsoid(EARTH_RADIUS, EARTH_FLATTENING, itrf))
.withLocalSolarTimeMode(NRLMSISE00.LocalSolarTimeMode.MEAN)
.withSwitch(9, -1),
new IsotropicDrag(AREA, CD)));
}
if (name.equals("drag-msis")) {
p.addForceModel(new DragForce(new NRLMSISE00(
new UnroundedWeather(new CssiSpaceWeatherData(
CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES)),
(ExtendedPositionProvider) CelestialBodyFactory.getSun(),
new OneAxisEllipsoid(EARTH_RADIUS, EARTH_FLATTENING, itrf))
.withLocalSolarTimeMode(NRLMSISE00.LocalSolarTimeMode.MEAN)
.withSwitch(9, -1),
new IsotropicDrag(AREA, CD)));
}
p.setInitialState(new SpacecraftState(orbit).withMass(MASS));
return p;
}
/** Runs an already-configured propagator across the arc and formats the states.
*
* Sampled with a fixed step handler, which interpolates inside the integrator's own
* steps, rather than by propagating to each output time in turn. The latter forces
* the controller to land exactly on every sample and changes the answer: at degree
* 70 over a day it moved the endpoint by 7 cm, which is the size of the disagreement
* these fixtures exist to measure.
*/
static String states(NumericalPropagator propagator, AbsoluteDate epoch) {
StringBuilder states = new StringBuilder();
propagator.setStepHandler(STEP_SECONDS, state -> {
Vector3D p = state.getPosition();
Vector3D v = state.getPVCoordinates().getVelocity();
if (states.length() > 0) {
states.append(",\n");
}
states.append(String.format(Locale.ROOT,
" { \"elapsed\": %.1f, \"position\": [%.9e, %.9e, %.9e], "
+ "\"velocity\": [%.9e, %.9e, %.9e] }",
state.getDate().durationFrom(epoch),
p.getX(), p.getY(), p.getZ(), v.getX(), v.getY(), v.getZ()));
});
propagator.propagate(epoch, epoch.shiftedBy(ARC_SECONDS));
return states.toString();
}
/** A point mass and nothing else: the case that validates the comparison itself. */
static void twoBody(File out) throws Exception {
TimeScale utc = TimeScalesFactory.getUTC();
AbsoluteDate epoch = new AbsoluteDate(2026, 1, 1, 0, 0, 0.0, utc);
Frame frame = FramesFactory.getEME2000();
Orbit orbit = new CartesianOrbit(
new PVCoordinates(new Vector3D(6.978e6, 0.0, 0.0), new Vector3D(0.0, 7.56e3, 0.0)),
frame, epoch, MU);
double[][] tol = NumericalPropagator.tolerances(TARGET_METERS, orbit, OrbitType.CARTESIAN);
DormandPrince853Integrator integrator =
new DormandPrince853Integrator(1.0e-6, 600.0, tol[0], tol[1]);
NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(OrbitType.CARTESIAN);
propagator.addForceModel(new NewtonianAttraction(MU));
propagator.setInitialState(new SpacecraftState(orbit));
write(new File(out, "two-body-low-earth.json"), "two-body-low-earth",
"2026-01-01T00:00:00 UTC", "GCRF",
new String[] { "NewtonianAttraction(mu = 3.986004415e14)" },
String.format(Locale.ROOT,
"DormandPrince853, tolerances from NumericalPropagator.tolerances(%.0e m)",
TARGET_METERS),
new String[] {
"GCRF, which is ICRF-aligned and so matches the inertial frame this"
+ " crate gets from ANISE. Orekit's EME2000 is the FK5 mean equator"
+ " and equinox instead, 23 mas away, which is the largest frame"
+ " difference there is between the two tools.",
"Inertial throughout, so no Earth orientation is involved.",
"No third body, no atmosphere, no radiation pressure.",
"mu is EGM96's, matching this crate's Earth.",
},
states(propagator, epoch));
}
static void write(File file, String name, String epoch, String frame, String[] forces,
String integrator, String[] notes, String states) throws Exception {
try (PrintWriter w = new PrintWriter(file, "UTF-8")) {
w.println("{");
w.printf(" \"name\": \"%s\",%n", name);
w.println(" \"setup\": {");
w.println(" \"tool\": \"orekit\",");
w.printf(" \"version\": \"%s\",%n", version);
w.printf(" \"generated\": \"%s\",%n", generated);
w.printf(" \"forces\": [%s],%n", quoted(forces));
w.printf(" \"integrator\": \"%s\",%n", integrator);
String[] all = new String[notes.length + 1];
System.arraycopy(notes, 0, all, 0, notes.length);
all[notes.length] = "orekit-data sha256 " + dataDigest;
w.printf(" \"notes\": [%s]%n", quoted(all));
w.println(" },");
w.printf(" \"epoch\": \"%s\",%n", epoch);
w.printf(" \"frame\": \"%s\",%n", frame);
w.println(" \"states\": [");
w.println(states);
w.println(" ]");
w.println("}");
}
System.out.println("wrote " + file.getName());
}
static String quoted(String[] items) {
StringBuilder b = new StringBuilder();
for (String item : items) {
if (b.length() > 0) {
b.append(", ");
}
b.append('"').append(item.replace("\"", "\\\"")).append('"');
}
return b.toString();
}
}
JAVA
javac -nowarn -cp "$work/lib/*" -d "$work/classes" "$work/Generate.java"
mkdir -p "$out"
java -cp "$work/lib/*:$work/classes" Generate \
"$work/gravity" "$data" "$out" "$OREKIT" "$data_digest" "$(date +%Y-%m-%d)" "$mode"
if [ "$mode" != "time" ]; then
echo "fixtures in ${out#"$root"/}"
fi
Reference in new issue View Git Blame Copy Permalink