Hello again, sorry for the late reply 
I did try what you suggested, interestingly enough, if I compare the equinoctial elements of the two states they do indeed match.
However, once again, if I try to convert the orbit to a Keplerian one, or if I use the getI() method of SpacecraftState, the inclination doesn’t match.
I obtain the same result if I use Hx and Hy to manually compute the inclination. I have literally no clue as to what’s wrong.
Best regards,
Emiliano
Hello Emiliano,
Are you using mean elements in your conversion?
If so, then I believe that is the problem. The “average” property is bound to the type of coordinates. So if you convert them in another system, you introduce an error. I think we can summarize this by “non linear maps and integrals are not interchangeable”. This is why comparing Cartesian parameters showed a discrepancy.
Best,
Romain.
Hello Romain,
Yes, what I meant before was that I was converting the mean equinoctial elements to Keplerian ones.
I’m not very familiar with the DSST theory, but, given your answer, I get that mean Keplerian elements don’t really exist and the nomenclature of “mean elements” is strictly bounded to the equinoctial elements. Is that correct?
Best regards,
Emiliano
Hi @Emiliano,
That’s very weird indeed…
Could you please zoom in on hx/hy plots so that we can better see the differences between mean and osculating, and also plot the variation of right ascension at ascending node?
Thanks,
Maxime
Hello Maxime,
Here’s the plots of hx and hy over 14 days:
This is the same plot but for a propagation of 1 year (not much to see ):
These instead are the plots of the difference between osculating and mean hx and hy, again over 14 days and 1 year:
Finally, these are the plots for the RAAN, again for 14 days and 1 year:
I also would like to mention that all these plots are generated considering only the geopotential expansion (4,4) as perturbing force. Indeed, by previous analyses I noticed that this is the term that causes this behaviour. Indeed, the DSSTPropagator was perfectly able to capture the dynamic evolution caused by the other perturbations. Moreover, this error becomes less and less visible with increasing altitude, again suggesting that the cause is the geopotential expansion.
Thank you very much for helping me.
Best regards,
Emiliano
Yes, that is my understanding. To be more precise, the only rigorous comparison are between osculating elements or between mean ones, but only when obtained straight by averaging the osculating ones according to the corresponding theory (so it can only be on the type of coordinates associated to it e.g. equinoctial for DSST). I think that anything else would in general introduce an “algebraic” error. That being said, it can be small, and maybe even negligible in many cases.
On another note, let’s keep in mind that the mean-to-osculating conversion (and the reciprocal by dependence) is only valid to a certain extent as there is always some approximation somewhere, in particular with complex perturbations. A good illustration is what Bryan said about Week Time Dependent terms in Orekit’s implementation of Third bodies. Also the second-order effects of J_2 are only available in the develop branch (have you tried it?). So you might just be seeing the limits of the DSST’s short-term geopotential, especially the tesseral harmonics given their simplified Earth rotation pointed out by Luc. It’s actually a shame (and perhaps not a coincidence) that we cannot simulate a 2x2 or 3x3 expansion.
Romain.
I see, thanks for the explaination.
Also, thanks for pointing me towards some known issues about the DSST implementation. Maybe you’re right, and what I’m seeing is just another example of those issues. Honestly, I have no idea
BTW, no I haven’t tried the develop branch because I think that’s only available for the Java version, and unfortunately I’m not very good in Java, so I tend to stay away from it as much as I can 
Best regards,
Emiliano
Is this still a problem in Orekit 11.3.2? I would say yes as I have been testing the DSST propagator to compare against a numerical one with similar perturbations, and adding the Harmonics (even just the J2) makes a big difference (at least when comparing the position):
I have checked the (Osculating) equinoctial elements to see where are the differences, and it appears to be mainly in Hx and Hy (this is after 5 days):
Some clarifications:
setInitialState. Also, I include PropagationType.OSCULATING in the creation of the DSSTPropagator to get the output in osculating elements and compare with the other propagators (I hope this does not create the type of errors in position that I am seeing…)I have concluded that the problem, if there is one, has to be with the harmonics part of the propagation, without it I get:
which now could be attributed to the tolerances being greater, or just a product of the aproximation of the DSST theory (which I don’t completely understand at the moment).
The purpose of these test are to find a compromise in the propagator that allows me to get faster approximations of the prediction for maneuver detection. And the most strange thing to me is that I am getting higher propagation times with the DSST:
I think the increase in computation comes only from the addition of the Atmosphere model (which is an empirical one and expensive to compute), but I don’t understand why it affects so much the DSST propagator. This are the computation times without atmosphere perturbation (the rest of the perturbations are Moon, Sun and SRP):
If it is the case that DSST has not changed in 11.3.2, and that 12.0 will inlcude the “fixes” for the difference with the numerical one, then I probably won’t be able to use it for my current project (I don’t know how to compile the develop branch for use in Matlab).
Thanks in advance for any acclaration.
Hi @astror,
I recently opened issue 1104 related to the discrepancy in osculating elements.
Could you please share your initial orbit(s), force models etc. so we can reproduce your test.
Or a Junit test, that would be the best.
Yes the atmospheric drag computation is not very performant with DSST. Still, I’m surprised you get such an overhead with the DSST.
Maxime
Hi there,
As you suspect, the J2 square thing will only be available in the next major release. This and other approximation/errors in the DSST model might explain the discrepancy.
As for performance, non conservative forces are treated by quadrature so should be much slower than the others. You should be able to tune the number of quadrature points if I’m not wrong.
Moreover, in my experience, if you require many intermediate points and you have many perturbations taken into account, propagation in osculating elements is computationally expensive. It can make you lose the time advantage gained by having large integration steps. If mean elements are enough for you, I would suggest propagating only them. Note also that I’m not sure what happens with DSST when you have events resetting the state or it’s derivatives e.g. impulsive maneuvers.
Cheers,
Romain
Hi again,
I would not recommend RESET_STATE with an osculating DSST until this behavior is fixed.
I have some leads on how to solve this but I’m not very happy with it yet.
Maxime
Hi, this is a part of the Matlab script that I used to get those results:
% TESTING PROPAGATORS
import org.orekit.attitudes.*
import org.orekit.bodies.*
import org.orekit.forces.gravity.potential.*
import org.orekit.frames.*
import org.orekit.orbits.*
import org.orekit.time.*
import org.orekit.utils.*
import org.hipparchus.geometry.euclidean.threed.*
import org.orekit.models.earth.atmosphere.data.*
import org.orekit.propagation.numerical.*
import org.orekit.forces.*
import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation.*
import org.hipparchus.ode.nonstiff.*
import org.orekit.propagation.*
import org.orekit.models.earth.*
import org.orekit.propagation.semianalytical.dsst.DSSTPropagator
import org.orekit.propagation.semianalytical.dsst.forces.*
crawler = org.orekit.data.ZipJarCrawler('orekit-data.zip');
DM = org.orekit.data.DataContext.getDefault().getDataProvidersManager();
DM.addProvider(crawler);
% Obtain constants
mu = Constants.WGS84_EARTH_MU;
Re = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
f = Constants.WGS84_EARTH_FLATTENING;
% Define constants for SRP
dRef = 149597870000.0;
pRef = 4.5605*1E-06;
% Define frames and time
ITRF = FramesFactory.getITRF(ITRFVersion.ITRF_2000,IERSConventions.IERS_2010, 1); % ITRF_2000 as in the OEMs
earth = org.orekit.bodies.OneAxisEllipsoid(Re, f, ITRF);
rotatingFrame = earth.getBodyFrame();
inertialFrame = FramesFactory.getEME2000();
% Initial data
t0 = datetime(2023,07,28,9,30,0);
tref = date_conversion(t0,'matlab2orekit');
a = 6900e3;
e = 0.01;
i = 60*pi/180;
w = 10*pi/180;
ra = 30*pi/180;
v = 0;
% % INITIAL STATE
orb0 = KeplerianOrbit(a,e,i,w,ra,v,PositionAngle.TRUE,inertialFrame,tref,mu);
orb0_eq = EquinoctialOrbit(orb0);
orb0_car = CartesianOrbit(orb0);
pv0 = orb0.getPVCoordinates();
% REFERENCE NUMERICAL ORBIT
dP = 0.00001; % POSITION TOLERANCE
minStep = 0.001;
maxStep = 500;
initStep = 60;
% tolerance_kep = NumericalPropagator.tolerances(dP, orb0, OrbitType.KEPLERIAN);
tolerance_equ = NumericalPropagator.tolerances(dP, orb0_eq, OrbitType.EQUINOCTIAL);
% tolerance_car = NumericalPropagator.tolerances(dP, orb0_car, OrbitType.CARTESIAN);
TOL = tolerance_equ;
integrator = DormandPrince853Integrator(minStep, maxStep, TOL(1,:), TOL(2,:));
integrator.setInitialStepSize(initStep);
propagator_ref1 = NumericalPropagator(integrator);
% REFERENCE NUMERICAL ORBIT 2
dP = 0.0001;
minStep = 0.001;
maxStep = 500;
initStep = 60;
% tolerance_kep = NumericalPropagator.tolerances(dP, orb0, OrbitType.KEPLERIAN);
tolerance_equ = NumericalPropagator.tolerances(dP, orb0_eq, OrbitType.EQUINOCTIAL);
% tolerance_car = NumericalPropagator.tolerances(dP, orb0_car, OrbitType.CARTESIAN);
TOL = tolerance_equ;
integrator = DormandPrince853Integrator(minStep, maxStep, TOL(1,:), TOL(2,:));
integrator.setInitialStepSize(initStep);
propagator_ref2 = NumericalPropagator(integrator);
% INITIAL SPACECRAFT
SS0 = SpacecraftState(orb0_eq,1000);
% DEFINE PERTURBATIONS FOR NUMERICAL PROPAGATOR
CD = 0.8; S = 10; SRPCoef = 1.4; Ssrp = 10;
% atmosphere
msafe = MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,javaMethod('valueOf','org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation$StrengthLevel', 'AVERAGE'));
earthAtmosphere = atmosphere.NRLMSISE00(msafe, CelestialBodyFactory.getSun(),earth);
dragF = drag.IsotropicDrag(S,CD);
dragForce = drag.DragForce(earthAtmosphere,dragF);
% earthAtmosphere_SIMPLE = atmosphere.HarrisPriester(CelestialBodyFactory.getSun(),earth);
% dragForce_SIMPLE = drag.DragForce(earthAtmosphere_SIMPLE,dragF);
sunPert = gravity.ThirdBodyAttraction(CelestialBodyFactory.getSun());
moonPert = gravity.ThirdBodyAttraction(CelestialBodyFactory.getMoon());
degree = 2; order = 0;
norm_provider = get_grav_provider(degree,order,true);
harmonicsPert1 = gravity.HolmesFeatherstoneAttractionModel(ITRF, norm_provider);
norm_provider_2 = get_grav_provider(degree,order,true);
harmonicsPert2 = gravity.HolmesFeatherstoneAttractionModel(ITRF, norm_provider_2);
rad = radiation.IsotropicRadiationSingleCoefficient(SRPCoef,Ssrp);
srpPert = radiation.SolarRadiationPressure(dRef,pRef,CelestialBodyFactory.getSun(),Re,rad);
% REF1
propagator_ref1.setInitialState(SS0);
propagator_ref1.setMu(mu);
propagator_ref1.setOrbitType(orb0_eq.getType());
propagator_ref1.addForceModel(harmonicsPert1);
propagator_ref1.addForceModel(srpPert);
propagator_ref1.addForceModel(dragForce);
propagator_ref1.addForceModel(sunPert);
propagator_ref1.addForceModel(moonPert);
% REF2
propagator_ref2.setInitialState(SS0);
propagator_ref2.setMu(mu);
propagator_ref2.setOrbitType(orb0_eq.getType());
propagator_ref2.addForceModel(harmonicsPert2);
propagator_ref2.addForceModel(srpPert);
propagator_ref2.addForceModel(dragForce);
propagator_ref2.addForceModel(sunPert);
propagator_ref2.addForceModel(moonPert);
% % % DSST PROPAGATOR
minStep = SS0.getKeplerianPeriod()/10;
% minStep = 10;
maxStep = 100. * minStep;
tol = DSSTPropagator.tolerances(dP*10, orb0_eq);
integrator = DormandPrince853Integrator(minStep, maxStep, tol(1,:), tol(2,:));
% dsstProp = DSSTPropagator(integrator);
dsstProp = DSSTPropagator(integrator, PropagationType.OSCULATING);
dsstProp.setInitialState(SS0, PropagationType.OSCULATING);
% dsstProp.setInitialState(SS0);
% Central Body geopotential
unnorm_provider = get_grav_provider(degree,order,false);
dsstProp.setMu(mu);
dsstProp.addForceModel(DSSTZonal(unnorm_provider));
dsstProp.addForceModel(DSSTTesseral(rotatingFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, unnorm_provider));
dsstProp.addForceModel(DSSTAtmosphericDrag(earthAtmosphere, CD, S, unnorm_provider.getMu()));
dsstProp.addForceModel(DSSTSolarRadiationPressure(dRef,pRef, SRPCoef, Ssrp, CelestialBodyFactory.getSun(), Re, unnorm_provider.getMu()));
dsstProp.addForceModel(DSSTThirdBody(CelestialBodyFactory.getSun(), unnorm_provider.getMu()));
dsstProp.addForceModel(DSSTThirdBody(CelestialBodyFactory.getMoon(), unnorm_provider.getMu()));
% Propagate REFERENCE PROPAGATORS
tf = t0 + days(5);
tf_ref = date_conversion(tf,'matlab2orekit');
tic
generator1 = propagator_ref1.getEphemerisGenerator();
propagator_ref1.propagate(tf_ref);
ephemerisref1 = generator1.getGeneratedEphemeris();
t_prop_ref1 = toc;
tic
generator2 = propagator_ref2.getEphemerisGenerator();
propagator_ref2.propagate(tf_ref);
ephemerisref2 = generator2.getGeneratedEphemeris();
t_prop_ref2 = toc;
tic
generator3 = dsstProp.getEphemerisGenerator();
p = dsstProp.propagate(tf_ref);
ephemerisref3 = generator3.getGeneratedEphemeris();
t_prop_ref3 = toc;
fprintf('Reference propagator 1 took: %.2f seconds \n',t_prop_ref1)
fprintf('Reference propagator 2 took: %.2f seconds \n',t_prop_ref2)
fprintf('Reference propagator 3 (DSST) took: %.2f seconds \n',t_prop_ref3)
N = 1500;
t_vector = linspace(t0,tf,N);
pos_error = zeros(size(t_vector));
pos_error_dsst = zeros(size(t_vector));
ref1_elems = zeros(N,6);
ref2_elems = zeros(N,6);
ref3_elems = zeros(N,6);
for i = 1:N
ti = t_vector(i);
tref_i = date_conversion(ti,'matlab2orekit');
sc1 = ephemerisref1.propagate(tref_i);
sc2 = ephemerisref2.propagate(tref_i);
sc3 = ephemerisref3.propagate(tref_i);
ref1_elems(i,:) = [sc1.getA(),sc1.getEquinoctialEx(),sc1.getEquinoctialEy(),sc1.getHx(),sc1.getHy(),sc1.getLv()];
ref2_elems(i,:) = [sc2.getA(),sc2.getEquinoctialEx(),sc2.getEquinoctialEy(),sc2.getHx(),sc2.getHy(),sc2.getLv()];
ref3_elems(i,:) = [sc3.getA(),sc3.getEquinoctialEx(),sc3.getEquinoctialEy(),sc3.getHx(),sc3.getHy(),sc3.getLv()];
pos1 = sc1.getPVCoordinates.getPosition.toArray();
pos2 = sc2.getPVCoordinates.getPosition.toArray();
pos3 = sc3.getPVCoordinates.getPosition.toArray();
pos_dif = pos2-pos1;
pos_dif_dsst = pos3-pos1;
pos_error(i) = sqrt(sum(pos_dif.^2));
pos_error_dsst(i) = sqrt(sum(pos_dif_dsst.^2));
end
I know there are a couple functions in there that are not specified anywhere, but one is for date conversion between Matlab and Orekit format, and the other is just a shortcut for calling the gravity field provider with the degree and order specified (either normalized or unnormalized). Aside from that is just a script for comparing the error in instantaneous position (and osculating elements) between DSST and a numerical propagator, all with the same perturbations (or I thought they were the same).
If something else is needed to reproduce the results just ask for it.
Thanks in advance.
Thank you @Serrof , actually this test was my first try at the DSST propagator so I really did not know what I would find. Actually, I don’t really need many intermediate points, but ploting many intermediate states between start and end of the propagation is always what I do when I want to get a feel at how the error grows. In reality I would only need the state (not in mean, but osculating) at the end of the propagation (and an estimation of the covariance, but that is another story), and wanted to see if this propagator would be a good alternative/compromise between precision and speed. For my application I need both the harmonics and a good atmosphere model (the orbits are too low to ignore those), so maybe this DSST propagator is not a good fit for my needs.