Differences in numerical propagation between Python and Java

Orekit Python Wrapper usage
1

Hi folks,

during the GTOC I encountered a strange behaviour that I’m reporting there.
In a nutshell, if I define the same exact integration settings in Python (jpype wrapper with additional .jar) and Java, I don’t get the same propagated Orbit at the end, even in Keplerian motion.
The most worrying thing is that it’s the Python version that is the closest to the “truth” from analytical propagation (which seems to be the same in both langages). I’ve been able to understand that in Java, the iintegration steps performed are much larger than in Python, as though the tolerances were less tight.

Maybe there is a subtle difference in the code that I missed, or something else is afoot.

Python code:

DAY = 86400

def build_propagator(orbit) -> NumericalPropagator:
    orbit_type = OrbitType.CARTESIAN
    min_step = 1.e-3
    max_step = DAY * 100.
    tol_pos, tol_vel = 1e1, 1e-5  # m, m/s
    rel_tol = 1e-6
    cartesian_tolerance = GtocCartesianToleranceProvider.of(rel_tol, tol_pos, tol_vel)
    integrator_builder = DormandPrince853IntegratorBuilder(min_step, max_step, ToleranceProvider.of(cartesian_tolerance))
    integrator = integrator_builder.buildIntegrator(orbit, orbit_type)
    integrator.setInitialStepSize(max_step / 10.)
    propagator = NumericalPropagator(integrator, CheapAttitudeProvider())
    propagator.setOrbitType(orbit_type)
    propagator.setInitialState(SpacecraftState(orbit))
    propagator.setResetAtEnd(True)
    return propagator

Java:

    private static final double MIN_STEP = 1e-3;
    private static final double MAX_STEP = 86400. * 100.;
    private static final double REL_TOL = 1e-6;

    public static NumericalPropagator getPropagator(final Orbit orbit) {
        return getPropagator(orbit, OrbitType.CARTESIAN);
    }

    public static NumericalPropagator getPropagator(final Orbit orbit, final OrbitType orbitType) {
        final DormandPrince853IntegratorBuilder integratorBuilder = new DormandPrince853IntegratorBuilder(MIN_STEP, MAX_STEP,
                getToleranceProvider());
        final DormandPrince853Integrator integrator = integratorBuilder.buildIntegrator(orbit, orbitType);
        integrator.setInitialStepSize(MAX_STEP / 10.);
        final NumericalPropagator propagator = new NumericalPropagator(integrator, new CheapAttitudeProvider());
        propagator.setOrbitType(orbitType);
        propagator.setInitialState(new SpacecraftState(orbit));
        propagator.setResetAtEnd(true);
        return propagator;
    }

    static ToleranceProvider getToleranceProvider() {
        return ToleranceProvider.of(new GtocCartesianToleranceProvider(REL_TOL, 1e1, 1e-5));
    }

public class GtocCartesianToleranceProvider implements CartesianToleranceProvider {

    private final double relativeTolerance;
    private final double absoluteTolerancePosition;
    private final double absoluteToleranceVelocity;

    GtocCartesianToleranceProvider(final double relativeTolerance, final double absoluteTolerancePosition,
                                   final double absoluteToleranceVelocity) {
        this.relativeTolerance = relativeTolerance;
        this.absoluteTolerancePosition = absoluteTolerancePosition;
        this.absoluteToleranceVelocity = absoluteToleranceVelocity;
    }

    @Override
    public double[][] getTolerances(Vector3D position, Vector3D velocity) {
        final double[] absTol = new double[7];
        final double[] relTol = absTol.clone();
        Arrays.fill(absTol, 0, 3, absoluteTolerancePosition);
        Arrays.fill(absTol, 3, 6, absoluteToleranceVelocity);
        absTol[6] = 1.;
        Arrays.fill(relTol, 0, 7, relativeTolerance);
        return new double[][] { absTol, relTol };
    }
}

The MU used in Python is 139348062043.343 * 1e9 and in Java 1.39348062043343E20.
Initial conditions:

        final double[] initialKm = new double[] {-1009643439.3646595,-1875136175.8096392,	13802415.843832716,
                8.897037079480848,	6.806874624952856,	-0.07805306755722974};
        final RealVector vector = MatrixUtils.createRealVector(initialKm).mapMultiply(1e3);
        final double tof =222604888.3859024;
        final CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(new Vector3D(vector.getSubVector(0, 3).toArray()),
                new Vector3D(vector.getSubVector(3, 3).toArray())), FramesFactory.getGCRF(), AbsoluteDate.ARBITRARY_EPOCH,
                GM_SI);

Cheers,
Romain.

2

That is indeed surprising. Are you getting “significantly” different results? There is always a risk in float representation, but what I have understood the double in java is very close to a python float, but I would not assume they are bit-identical. But this shouldn’t make a big difference in the propagation.

3

Yes the difference is significant, even after one integration step. I’ll try to make a more thorough comparison.

Cheers,
Romain

4

OK so I found the problem…

In Python, instead of calling the constructor of GtocCartesianToleranceProvider, I was calling the static method of inherited from the interface. So basically my implementation was totally ignored, and worse the object had random inputs. I thought my integration settings were fine when in actual fact, the tolerance were pretty sketchy and probably not appropriate for some orbit regimes. Oh well…

5 Likes