Hi everyone,
I encountered a bug related to the osculating/mean conversion: starting from the same osculating orbit (defined in EME2000 frame in the test below), the output mean orbit (expressed in EME2000) shows a difference, mainly in semi-major axis (~24 m), if converted directly in EME2000 vs. a conversion in TOD frame.
Such discrepancy appears even when considering only J2 as perturbation.
Issue #1104 seems to confirm this bug. In this particular case, though, no propagation is performed and only the simple conversion shows the discrepancy.
This forces to convert the osculating orbit to TOD, convert it to mean, and convert the output mean orbit to the initial frame again to obtain a result closer to numerical propagation (according to #1104).
I leave here below a test highlighting the difference.
public static final double J2 = 0.0010826262155533902;
public static final int DEGREE = 2;
public static final double EARTH_MU = 3.986004415E14;
public static final double EARTH_RADIUS = 6378136.46;
@Test
void testDSSTConversion() {
// Set data context
DataContext dataContext = createDataContext();
Frame eme2000 = dataContext.getFrames().getEME2000();
Frame tod = dataContext.getFrames().getTOD(false);
UnnormalizedSphericalHarmonicsProvider provider = dataContext.getGravityFields().getUnnormalizedProvider(DEGREE, 0);
double earthMu = provider.getMu();
// Initialize osculating orbit in EME2000
AbsoluteDate date = new AbsoluteDate(2024, 9, 9, 0, 0, 0.0, dataContext.getTimeScales().getTAI());
KeplerianOrbit osculatingKepOrbit = new KeplerianOrbit(7000e3, 0.1, Math.toRadians(100), Math.toRadians(4), Math.toRadians(10), Math.toRadians(8), PositionAngleType.TRUE,
eme2000, date, earthMu);
// Initialize DSST converter
List<DSSTForceModel> forceModels = List.of(new DSSTZonal(provider));
OsculatingToMeanConverter converter = new FixedPointConverter(1.0E-11, 200, 1.0);
AttitudeProvider dummyUnusedProvider = FrameAlignedProvider.of(eme2000);
DSSTTheory meanTheory = new DSSTTheory(forceModels, dummyUnusedProvider, 0);
converter.setMeanTheory(meanTheory);
// Convert orbit to mean (output always in EME2000)
KeplerianOrbit meanOrbitConvertedInEme2000OutputInEme2000 = (KeplerianOrbit) converter.convertToMean(osculatingKepOrbit);
KeplerianOrbit meanOrbitConvertedInTodOutputInEme2000 = (KeplerianOrbit) converter.convertToMean(osculatingKepOrbit.inFrame(tod)).inFrame(eme2000);
// Check differences in orbital parameters
Assertions.assertEquals(23.45318115130067, meanOrbitConvertedInTodOutputInEme2000.getA() - meanOrbitConvertedInEme2000OutputInEme2000.getA());
Assertions.assertEquals(3.5301352648287043E-6, meanOrbitConvertedInTodOutputInEme2000.getEquinoctialEx() - meanOrbitConvertedInEme2000OutputInEme2000.getEquinoctialEx());
Assertions.assertEquals(-1.3023342389055503E-6, meanOrbitConvertedInTodOutputInEme2000.getEquinoctialEy() - meanOrbitConvertedInEme2000OutputInEme2000.getEquinoctialEy());
Assertions.assertEquals(-5.262301536745895E-7, meanOrbitConvertedInTodOutputInEme2000.getHx() - meanOrbitConvertedInEme2000OutputInEme2000.getHx());
Assertions.assertEquals(1.5438937456258017E-7, meanOrbitConvertedInTodOutputInEme2000.getHy() - meanOrbitConvertedInEme2000OutputInEme2000.getHy());
Assertions.assertEquals(-2.8807629653959665E-7, meanOrbitConvertedInTodOutputInEme2000.getLv() - meanOrbitConvertedInEme2000OutputInEme2000.getLv());
}
private DataContext createDataContext() {
return new DataContext() {
@Override
public TimeScales getTimeScales() {
return new LazyLoadedDataContext().getTimeScales();
}
@Override
public org.orekit.frames.Frames getFrames() {
return new LazyLoadedDataContext().getFrames();
}
@Override
public CelestialBodies getCelestialBodies() {
return null;
}
@Override
public GravityFields getGravityFields() {
return new GravityFields() {
@Override
public NormalizedSphericalHarmonicsProvider getConstantNormalizedProvider(int degree, int order, AbsoluteDate freezingDate) {
return null;
}
@Override
public NormalizedSphericalHarmonicsProvider getNormalizedProvider(int degree, int order) {
return null;
}
@Override
public UnnormalizedSphericalHarmonicsProvider getConstantUnnormalizedProvider(int degree, int order, AbsoluteDate freezingDate) {
return null;
}
@Override
public UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(int degree, int order) {
return new UnnormalizedSphericalHarmonicsProvider() {
@Override
public UnnormalizedSphericalHarmonics onDate(AbsoluteDate date) {
return new UnnormalizedSphericalHarmonics() {
@Override
public double getUnnormalizedCnm(int n, int m) {
return n == DEGREE ? -J2 : 0.0;
}
@Override
public double getUnnormalizedSnm(int n, int m) {
return 0;
}
@Override
public AbsoluteDate getDate() {
return null;
}
};
}
@Override
public int getMaxDegree() {
return DEGREE;
}
@Override
public int getMaxOrder() {
return 0;
}
@Override
public double getMu() {
return EARTH_MU;
}
@Override
public double getAe() {
return EARTH_RADIUS;
}
@Override
public AbsoluteDate getReferenceDate() {
return null;
}
@Override
public TideSystem getTideSystem() {
return null;
}
};
}
@Override
public List<OceanTidesWave> getOceanTidesWaves(int degree, int order) {
return null;
}
};
}
@Override
public GeoMagneticFields getGeoMagneticFields() {
return null;
}
};
}