Hello,
I am trying to propagate the covariance matrix using a DSST propagator. I noticed that introducing some force models (in my case the third body attraction and the zonal effect) the program crashes with an error saying that the “hansen” is null. This happens only if introducing the computation of the MatricesHarvester (as needed for the covariance propagation). The code below reproduces this behaviour:
// Settings
double a = 6860176.23765572;
double e = 0.0010160613326433334;
double i = 1.7044014447593043;
double om = 1.439233350035424;
double OM = 1.3884502116745805;
double th = 1.1625681691480643;
AbsoluteDate startingEpoch = new AbsoluteDate(2022, 01, 01, 12, 00, 0.0, TimeScalesFactory.getTAI());
AbsoluteDate finalEpoch = new AbsoluteDate(2022, 01, 02, 12, 00, 0.0, TimeScalesFactory.getTAI());
Orbit keplerianOrbit = new KeplerianOrbit(a, e, i, om, OM, th, PositionAngle.TRUE, FramesFactory.getGCRF(), startingEpoch, Constants.IERS2010_EARTH_MU);
// Propagator
DSSTPropagator dsstPropagator = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(60 * 20), PropagationType.MEAN);
dsstPropagator.setInitialState(new SpacecraftState(keplerianOrbit));
dsstPropagator.addForceModel(new DSSTNewtonianAttraction(Constants.IERS2010_EARTH_MU));
dsstPropagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getJupiter(), Constants.JPL_SSD_JUPITER_SYSTEM_GM));
MatricesHarvester matricesHarvester = dsstPropagator.setupMatricesComputation("stm", null, null);
// Propagate
dsstPropagator.propagate(finalEpoch);
Is there a reason why I get this error? Does the DSST propagation have some thoretical limits (like, for example, constraints in the orbital parameters, force models to be introduced, etc.) beyond which it is not safe to use it, at least when dealing with the MatricesHarvester?
You have the issue because the short periodic terms are not initialized for the “field” equations of the DSST. Because you perform covariance propagation, you need to compute partial derivatives to compute the state transition matrix. In Orekit, partial derivatives are not computed using finite differences or analytical formula. They are computed using automatic differentiation. To perform automatic differentiation, we had to implement “field” version of the equations. A “field” is an “extended” number. In other word, it contains the value of the parameter with additional information, like its partial derivatives with respect to user-defined parameters.
For more information about automatic differentiation in Orekit, you can look at these references [1] [2] [3].
The solution to you problem is an ugly, but working, work around. You have to do the following procedure:
Please not that we implemented for the 11.3 version a class allowing to perform covariance propagation with different type of orbit propagators in Orekit (numerical, DSST, brouwer-lyddane, eckstein-hechler, and keplerian). It avoid the previous workaround. The class is StateCovarianceMatrixProvider and it can be added to the orbit propagator using .addAdditionalState(...) method of the propagator. There is a tutorial showing how to use it: CovariancePropagation.java. Please not that this functionality is not officially released. It is currently available in the development version of Orekit.
DSST has the same limitations as the numerical propagator. It supports degree and order of the geopotential over 100 and 100. It supports all the third body attractions of all the bodies available in Orekit. It can model drag and solar radiation pressure. It is a really powerful tool.
I have few comments regarding the code:
The central attraction coefficient used to initialize the DSSTThirdBody is the one of the Earth, not the one of the planet. Therefore, you have to do the following fix:
You could add more force models like DSSTZonal and DSSTTesseral to represent the Earth potential. You could also add the attraction of the Moon and the Sun that are much more important than the one of Jupiter for this orbit. I recommend you to look at the following graph to help you having a consistent perturbation model. Below an example:
You could improve the computation time by using a variable step integrator. I recommend you to use a Dormand-Prince integrator with a min step equal to the orbital period and a max step equal to 10 orbital period. The integrator will compute the best step taking into account the accuracy and the computation time. Also, DSST is a semi-analytical propagator. Therefore, you can use step sizes equal to several orbit periods.
final double[][] tol = DSSTPropagator.tolerances(10.0, keplerianOrbit);
final ODEIntegrator integrator = new DormandPrince853Integrator(keplerianOrbit.getKeplerianPeriod(), 10 * keplerianOrbit.getKeplerianPeriod(), tol[0], tol[1]);
DSSTPropagator dsstPropagator = new DSSTPropagator(integrator, PropagationType.MEAN);
That’s something we have to discuss with the development team. The objective is to have a release by the end of the month (finger crossed).
We will provide updates about the release planning very soon
