Propagating Backward in Time Seems Like Giving More Positional Error

Orekit usage
1

Hi everyone,

I am trying to calculate the close approach date and miss distance value for a given conjunction. Therefore, I am propagating the given close approach state to backward in time lets say 2 days and then re-calculating the close approaches with ExtremumApproachEventDetector. Yet, I am ending up a slightly different close approach which means that I cannot reach the initial state that I have provided. It seems like the propagator adding some error while propagating backward therefore changing the orbital state, then propagating from this perturbed state to the known close approach ending up a different time and miss distance value.
Expected Time of Closest Approach and Miss Distance is:
TCA = 2014-07-16T22:32:26.716, MISS DISTANCE = 91.455902 m
Calculated Time of Closest Approach and Miss Distance is:
TCA = 2014-07-16T22:32:26.71574653075333Z MISS DISTANCE: 90.9158464184037 m
In the below, one can find the example of this scenario:

        String tcaSTR = "2014-07-16T22:32:26.716";
        AbsoluteDate tca =  new AbsoluteDate(tcaSTR, TimeScalesFactory.getUTC());

        Vector3D pPrimary = new Vector3D(6703.392053e3, -1969.223265e3,3276.274276e3);
        Vector3D vPrimary = new Vector3D( 3.52773277e3, 2.322656571e3, -5.815904333e3);

        Vector3D pSecondary = new Vector3D(6703.301055e3,-1969.231944e3,3276.271432e3);
        Vector3D vSecondary = new Vector3D(2.538967996e3,6.496491661e3,5.104787983e3);

        TimeStampedPVCoordinates pvPrimary = new TimeStampedPVCoordinates(tca,pPrimary,vPrimary);
        TimeStampedPVCoordinates pvSecondary = new TimeStampedPVCoordinates(tca,pSecondary,vSecondary);

        CartesianOrbit oPrimary = new CartesianOrbit(pvPrimary, FramesFactory.getEME2000(), Constants.WGS84_EARTH_MU);
        CartesianOrbit oSecondary = new CartesianOrbit(pvSecondary, FramesFactory.getEME2000(),Constants.WGS84_EARTH_MU);

        SpacecraftState statePrimary = new SpacecraftState(oPrimary);
        SpacecraftState stateSecondary = new SpacecraftState(oSecondary);

        final NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(70,70);
        HolmesFeatherstoneAttractionModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,false),provider);

        double minStep = 1e-3;
        double maxStep = 60;
        double positionTolerance = 1.0;
        OrbitType propagationType = OrbitType.CARTESIAN;
        double[][] tolerances1 = NumericalPropagator.tolerances(positionTolerance, statePrimary.getOrbit(), propagationType);
        double[][] tolerances2 = NumericalPropagator.tolerances(positionTolerance, stateSecondary.getOrbit(), propagationType);

        DormandPrince853Integrator integratorPrimary = new DormandPrince853Integrator(minStep, maxStep, tolerances1[0], tolerances1[1]);
        DormandPrince853Integrator integratorSecondary = new DormandPrince853Integrator(minStep, maxStep, tolerances2[0], tolerances2[1]);

        NumericalPropagator propagatorPrimary = new NumericalPropagator(integratorPrimary);
        propagatorPrimary.addForceModel(holmesFeatherstone);
        propagatorPrimary.setOrbitType(propagationType);
        propagatorPrimary.setInitialState(statePrimary);

        NumericalPropagator propagatorSecondary = new NumericalPropagator(integratorSecondary);
        propagatorSecondary.addForceModel(holmesFeatherstone);
        propagatorSecondary.setOrbitType(propagationType);
        propagatorSecondary.setInitialState(stateSecondary);

        EventHandler handler = new EventHandler() {
            @Override
            public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
                if(increasing){
                    AbsoluteDate date = s.getDate();
                    Vector3D pos1 = s.getPVCoordinates().getPosition();
                    Vector3D pos2 = propagatorSecondary.getPVCoordinates(date,s.getFrame()).getPosition();

                    System.out.println("CLOSE APPROACH DATE: " + date + " MISS DISTANCE: " + pos1.subtract(pos2).getNorm());
                }
                return Action.CONTINUE;
            }
        };

        ExtremumApproachDetector detector = new ExtremumApproachDetector(propagatorSecondary).withHandler(handler);

        double nbDaysFromTCA = -2*Constants.JULIAN_DAY;
        double nbDaysFromStart = -nbDaysFromTCA + oPrimary.getKeplerianPeriod();
        AbsoluteDate startDate = tca.shiftedBy(nbDaysFromTCA);
        AbsoluteDate endDate = startDate.shiftedBy(nbDaysFromStart);
        propagatorPrimary.addEventDetector(detector);
        propagatorPrimary.propagate(startDate,endDate);

If we propagate more in backward in time, the error is increasing. What is the cause of this phenomenon? Am I setting the numerical propagator tolerances wrongly? Is it possible to avoid this kind of event, because I am calculating the same operation with well-known applications such as STK, FreeFlyer and GMTA. The error is insignificant for those applications.

2

Hi there,

As written in the documentation, you shouldn’t use this detector with a non-analytical propagator for the secondary. Instead, perform propagation for it beforehand on a large enough time window and then pass the IntegratedEphemeris/BoundedPropagator (using the EphemerisGenerator) to the detector.

Cheers,
Romain.

3

Hi Romain,

Thank you for your answer,

Unfortunately, I have tried the ephemeris approach but it yields to even more error. I suppose the ephemeris is suggested in order to avoid infinite recursions. Eventhough, it helps to reduce the computation time, again, one must propagate backward in time and this is doubling the integration error since we have propagating backward and forward for the same point.

4

Indeed, setting the tolerances to 1.0m, propagating for several days, and expecting the close approach to be below 1m is almost certain to fail.
Tolerances are a difficult subject. They are used as a limit the integrator should comply with when it estimates the error before deciding to accept or reject one step. There are two major points in the previous sentence: estimates and one step.
The integrator cannot compute exactly the error, it can only estimate it. If it could compute the error, it would do it and then compensate it and the compensated computation would be perfect. This is not possible. It only tries to estimate the order of magnitude of the error based on the evolution of the state vector with various assumptions.
The integrator estimates the error only for one step, not throughout the integration, so it will accept a step when its error estimation is that the error on the current step is below the tolerance.
Here, you tell it that each step is allowed to have a 1m error, this may build up a lot to a large error after only a few minutes, so after two days I guess the orbit could almost be random. In fact, you were lucky to have results still at meter level with these settings.
There are no ways to estimate a global error from the local step tolerance settings, it depends on the differential equations. Some equations would exhibit very stable behavior with local and global errors remaining of same order of magnitude, whereas other equations would exhibit divergence.

As a rule of thumb, our current recommendation is to use 0.001m for numerical integration tolerances.