State transition matrix parametrization

Hello,

I am propagating the covariance matrix making use of the state transition matrix derived from the MatricesHarvester.
I just would like to make sure that the following is true:

• If KeplerianPropagator / EcksteinHechlerPropagator → the state transition matrix is in cartesian parametrization;
• if NumericalPropagator → the state transition matrix parametrization is function of the setOrbitType and setPositionAngleType methods;
• if DSSTPropagator → the state transition matrix is in equinoctial / mean anomaly parametrization.

Since I need to implement a general purpose code where the input covariance may be expressed in any parametrization, do I need therefore to transform the input covariance in the correct parametrization (depending on the propagator, through the getJacobianWrt___ methods), propagate it, and then transform it back? Is there any other simpler way to do this?

Thank you!

We recently implemented a class for that. It is currently available in the development version of Orekit and will be officially released in the 11.3 version.
This class is StateCovarianceMatrixProvider.java. It can be used with different orbit propagators (DSST, numerical, Eckstein-Hechler, Brouwer-Lyddane, SGP4, and Keplerian) and it can provide the covariance in different frames and orbit types.

You can find examples in the test class of this function: StateCovarianceMatrixProviderTest.java

Best regards,
Bryan

Hi @bcazabonne and thank you!

I’ll make use of that class then, as soon as it is available. Just for my understanding, is anyway true what I said?

Thank you.

True

By default it is equinoctial / true anomaly but yes you can update it using setOrbitType(OrbitType) and setPositionAngleType(PositionAngle) methods.

True

When computing the propagated covariance, the initial one must be indeed in the same orbit type than the state transition matrix. For the orbit type transformation, you can look the implementations of the changeCovarianceType method in StateCovarianceMatrixProvider.
So, if your state transition matrix is computed in cartesian elements but your input covariance uses keplerian elements, you need:

1. Transform the input covariance Cov_t0 in cartesian
2. Compute the propagated covariance using Cov_tk = STM * Cov_t0 * STM^T
3. Transform the propagated covariance Cov_tk in keplerian

Bryan

Thank you!

Hi,
I just have a last question…
I am trying to propagate using a DSST propagator. As soon as I invoke the method

MatricesHarvester matricesHarvester = dsstPropagator.setupMatricesComputation("stm", null, null);


the propagation does not work anymore.
Am I not supposed to extract the MatricesHarvester (needed for stm retrieval) in this way when propagating with a DSST propagator?

Thank you

Hi,

Could you try this:

    final DSSTHarvester harvester = (DSSTHarvester) propagator.setupMatricesComputation(stmName, null, null);
harvester.initializeFieldShortPeriodTerms(initialState);


Best regards,
Bryan

Thank you @bcazabonne,

this actually works, but it really takes too long if setting some dsst force models besides the newtonian attraction. Isn’t there any other alternative to this initializeFieldShortPeriodTerms?

Thank you!

Unfortunately no… I’m thinking about a solution.

Bryan

I am using a classical Runge Kutta with step 0.5 and I am propagating the orbit for more or less one hour.
If there is no alternative, I’ll just use a numerical propagator then. Thank you!

Marco

Ok, using DSST you can increase the integration step because in DSST theory, the integrated elements are the mean orbital elements. Typically, the integration step can be equal to several orbital periods.

Bryan

Hi,

I am confused by your answer here @bcazabonne. In my understanding, the variables propagated by an instance of KeplerianPropagator are the ones of the underlying initial state (via its Orbit and the corresponding shiftedBy method). So I assumed that the STM would be using the same coordinates. But you’re saying that it is always Cartesian?

Cheers,
Romain.

Hi Romain,

For analytical orbit propagators the type of the STM is defined in the AbstractAnalyticalMatricesHarvester class. If you look at the implementation of the setReferenceState method you will see that we use the cartesian derivatives.

We decided to use cartesian parameters for two reasons : (1) consistency between propagators (Keplerian propagator uses keplerian elements, Eckstein-Hechler uses the cartesian ones and Brouwer-Lyddane uses keplerian) and (2) to avoid possible singularities.

Now, based on your comment, what could be interesting is to add a method getOrbitType() and getPositionAngle() in the MatricesHarvester interface to help users knowing the type used.

Best regards,
Bryan

Hi Bryan,

thank you for the clarification. I would indeed appreciate some getters as I find this behaviour a bit counter-intuitive.
However, you confused me even more by stating that a Keplerian propagator uses Keplerian elements!!!
Again, I thought it used the shiftedBy method of the underlying Orbit, so if it is a CartesianOrbit, it uses shiftPVElliptic for non-hyperbolic cases.

Best,
Romain.

The main problem is that we don’t have  getOrbitType() and/or getPositionAngle() methods for analytical orbit propagators. So we had to make a choice and we opted for the Cartesian for the two previous reasons.

I’ll create an issue for adding the getter. It will only be applicable for the next major version because MatricesHarvester is an interface and having default implementations is not a good idea for these methods.

Right! I thought the return type of the propgateOrbit() method was a KeplerianOrbit. But it is indeed an Orbit, so there is no imposed type as for the Eckstein-Hechler or Brouwer-Lyddane propagators which return a CartesianOrbit or a KeplerianOrbit respectively.
My apologies.

Best regards,
Bryan

No worries. And to be clear, I appreciate that the architecture of Orekit is not a trivial problem and that you guys are tackling it as best as you can. This is a great library! Thanks for the amazing work.

Romain.

Hi @bcazabonne, I exploit this post to add a small point.
I tried to replicate the test “testWithAnalyticalPropagator()” in src/test/java/org/orekit/propagation/StateCovarianceMatrixProviderTest.java · develop · Orekit / Orekit · GitLab. The code I wrote is the following:

        // Initial state
AbsoluteDate initEpoch = new AbsoluteDate("2016-02-13T16:00:00.000", TimeScalesFactory.getUTC());
Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(new Vector3D(7526993.581890527, -9646310.10026971, 1464110.4928112086),
new Vector3D(3033.79456099698, 1715.265069098717, -4447.658745923895)),
FramesFactory.getEME2000(),
initEpoch,
Constants.WGS84_EARTH_MU);

double[][] initCov = new double[][] {{8.651816029e+01,   5.689987127e+01, -2.763870764e+01, -2.435617201e-02,  2.058274137e-02, -5.872883051e-03},
{5.689987127e+01,   7.070624321e+01,  1.367120909e+01, -6.112622013e-03,  7.623626008e-03, -1.239413190e-02},
{-2.763870764e+01,  1.367120909e+01,  1.811858898e+02,  3.143798992e-02, -4.963106559e-02, -7.420114385e-04},
{-2.435617201e-02, -6.112622013e-03,  3.143798992e-02,  4.657077389e-05,  1.469943634e-05,  3.328475593e-05},
{2.058274137e-02,   7.623626008e-03, -4.963106559e-02,  1.469943634e-05,  3.950715934e-05,  2.516044258e-05},
{-5.872883051e-03, -1.239413190e-02, -7.420114385e-04,  3.328475593e-05,  2.516044258e-05,  3.547466120e-05}};
RealMatrix initCovReal = MatrixUtils.createRealMatrix(initCov);

// Propagator
EcksteinHechlerPropagator prop = new EcksteinHechlerPropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(6, 0));
MatricesHarvester harvester = prop.setupMatricesComputation("stmName", null, null);

// Propagate and compute covariance
SpacecraftState propagatedState = prop.propagate(initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY));

RealMatrix dYdY0 = harvester.getStateTransitionMatrix(propagatedState);
RealMatrix propCov = dYdY0.multiply(initCovReal.multiplyTransposed(dYdY0));

System.out.println(propCov.getEntry(0, 0));


Running it, I get that the first element of the propagated covariance matrix is 577.3144349040394, which is not within the threshold of 5.0e-4 to your reference 5.770543135e+02. Is this just a numerical error or did I make a mistake in my code?
Thank you.

Hi @MarcoTosato,

I take the liberty to reply in place of @bcazabonne
If you take a close look at method StateCovarianceMatrixProviderTest#compareCovariance you’ll see that the threshold is relative.

Assert.assertEquals(reference.getEntry(row, column), computed.getEntry(row, column), FastMath.abs(threshold * reference.getEntry(row, column)));


So in the case of the first element of the matrix the absolute threshold is 5.10^{-4} \times 577.3144349040394 = 0.28866

Hope this helps,
Maxime

Thank you @MaximeJ