Hi 
I am writing some software that runs a numerical propagation using osculating elements, and can include propulsive manoeuvres. To aid understanding of how these manoeuvres affect the resultant orbit, I am working on converting the osculating elements to mean elements using the DSST theory, and exposing the mean elements as an output.
I am finding, however, that the results from DSSTPropagator.computeMeanState are not what I expect. The code below shows a minimal example of running a numerical propagation (with no manoeuvres, and a simple gravity field with degree and order = 10), and within a step handler capturing the osculating and mean elements (from both DSST and Eckstein-Hechler) and plotting the semi-major axis evolution. Using Eckstein-Hechler gives me the results I expected - I was surprised to find that DSST follows the osculating values almost exactly instead of giving a more “mean” result.
Code example:
{
// Load Orekit data etc. here
// gravitation coefficient
final double mu = 3.986004415e+14;
// rotation rate
final double earthOmega = 7.292115e-5;
// inertial frame
final Frame inertialFrame = FramesFactory.getEME2000();
// body frame
final Frame bodyFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
// Initial date
final AbsoluteDate initialDate =
new AbsoluteDate(2004, 01, 01, 23, 30, 0, TimeScalesFactory.getUTC());
// Initial orbit
final double a = 7000000; // semi major axis in meters
final double e = 0.01; // eccentricity
final double i = FastMath.toRadians(97); // inclination
final double omega = FastMath.toRadians(0); // perigee argument
final double raan = FastMath.toRadians(0); // right ascention of ascending node
final double lM = 0; // mean anomaly
final Orbit initialOrbit =
new KeplerianOrbit(
a, e, i, omega, raan, lM, PositionAngleType.MEAN, inertialFrame, initialDate, mu);
// Initial state definition
final SpacecraftState initialState = new SpacecraftState(initialOrbit);
// Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
final double minStep = 0.001;
final double maxstep = 1000.0;
final double positionTolerance = 10.0;
final OrbitType propagationType = OrbitType.KEPLERIAN;
final ToleranceProvider toleranceProvider =
ToleranceProvider.of(CartesianToleranceProvider.of(positionTolerance));
final double[][] tolerances = toleranceProvider.getTolerances(initialOrbit, propagationType);
final AdaptiveStepsizeIntegrator integrator =
new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
// Propagator
final NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);
// Force Model (reduced to perturbing gravity field)
final NormalizedSphericalHarmonicsProvider provider =
GravityFieldFactory.getNormalizedProvider(10, 10);
final ForceModel holmesFeatherstone =
new HolmesFeatherstoneAttractionModel(
FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
// Equivalent DSST gravity force
final DSSTTesseral dsstTesseral =
new DSSTTesseral(
bodyFrame, earthOmega, GravityFieldFactory.getUnnormalizedProvider(provider));
// Add force model to the propagator
propagator.addForceModel(holmesFeatherstone);
// Set up initial state in the propagator
propagator.setInitialState(initialState);
final List<Double> durations = new ArrayList<>();
final List<KeplerianOrbit> osculatingOrbits = new ArrayList<>();
final List<KeplerianOrbit> dsstOrbits = new ArrayList<>();
final List<KeplerianOrbit> ehOrbits = new ArrayList<>();
final OrekitFixedStepHandler stepHandler =
new OrekitFixedStepHandler() {
/** {@inheritDoc} */
@Override
public void handleStep(SpacecraftState s) {
durations.add(s.getDate().durationFrom(initialDate));
// Compute osculating orbits
Orbit osculatingOrbit = s.getOrbit();
osculatingOrbits.add((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(osculatingOrbit));
// Compute mean DSST orbits
SpacecraftState meanDSSTState =
DSSTPropagator.computeMeanState(
s,
null,
List.of(dsstTesseral),
FixedPointOsculatingToAveragedConverter.DEFAULT_EPSILON,
FixedPointOsculatingToAveragedConverter.DEFAULT_MAX_ITERATIONS);
Orbit meanDSSTOrbit = meanDSSTState.getOrbit();
dsstOrbits.add((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(meanDSSTOrbit));
// Compute mean Eckstein-Hechler orbits
CircularOrbit meanEHOrbit =
EcksteinHechlerPropagator.computeMeanOrbit(
s.getOrbit(),
GravityFieldFactory.getUnnormalizedProvider(provider),
GravityFieldFactory.getUnnormalizedProvider(provider).onDate(s.getDate()));
ehOrbits.add((KeplerianOrbit) OrbitType.KEPLERIAN.convertType(meanEHOrbit));
}
};
propagator.setStepHandler(120, stepHandler);
propagator.propagate(initialDate.shiftedBy(86400 / 2));
plotMeanVsOsculatingSemiMajorAxis(
durations,
osculatingOrbits,
Map.of("DSST", dsstOrbits, "Eckstein-Hechler", ehOrbits),
"mean_vs_osculating_semi_major_axis");
}
private static void plotMeanVsOsculatingSemiMajorAxis(
List<Double> durations,
List<KeplerianOrbit> osculatingOrbits,
Map<String, List<KeplerianOrbit>> meanOrbits,
String filename) {
List<Double> osculatingSMA = osculatingOrbits.stream().map(KeplerianOrbit::getA).toList();
XYChart chart =
new XYChartBuilder()
.width(1600)
.height(800)
.title("Mean vs Osculating Semi-Major Axis")
.xAxisTitle("Time (seconds)")
.yAxisTitle("Semi-Major Axis (m)")
.build();
XYSeries oscSeries = chart.addSeries("Osculating", durations, osculatingSMA);
oscSeries.setMarker(SeriesMarkers.CIRCLE);
for (String key : meanOrbits.keySet()) {
List<Double> meanSMA = meanOrbits.get(key).stream().map(KeplerianOrbit::getA).toList();
XYSeries series = chart.addSeries(key, durations, meanSMA);
series.setMarker(SeriesMarkers.CIRCLE);
}
try {
BitmapEncoder.saveBitmap(
chart,
filename,
BitmapEncoder.BitmapFormat.PNG);
} catch (IOException e) {
e.printStackTrace();
}
}
Plot:
Can anyone help me understand If I am calculating the mean DSST elements wrongly, or if there is some underlying issue here?
Thank you!
