Behavior of DSSTForceModel

icav · October 29, 2025, 1:36pm

I obtained a strange result in performing a propagation tests with the DSSTPropagator and I think the problem is related to the class DSSTForceModel. I initialized the propagator and added several DSSTForceModels (DSSTZonal, DSSTTesseral, DSSTThirdBody, DSSTAtmosphericDrag, DSSTSolarRadiationPressure) to it. After performing a first propagation up to a target date, I reset the same initial state as before in the propagator and repeated the same propagation again. The outcomes of the two propagations are different. I think this is related to the DSSTForceModel as consequence of a second test I performed. I repeated the same steps as before up to the second propagation. Before performing the second propagation, I removed the force models from the propagator, I recreated the DSSTForceModels and added them to the propagation. In this case, the outcomes of the two propagations are the same.

In both tests, I used the same inputs:

////////////////////////////////////////////////////////////////////////
/// INPUTS        ////////////////////////////////////////////////////////////////////////
        
// CONTEXT
IERSConventions conventions = IERSConventions.IERS_2010;
TimeScale timeScale = TimeScalesFactory.getUTC();
Frame pseudoInertialFrame = FramesFactory.getTOD(conventions,   false);
Frame rotatingEarthFrame = FramesFactory.getGTOD(conventions, false);
double earthAngularVelocity = Constants.IERS2010_EARTH_ANGULAR_VELOCITY;
double earthRadius = Constants.IERS2010_EARTH_EQUATORIAL_RADIUS;
double earthFlattening = Constants.IERS2010_EARTH_FLATTENING;
OneAxisEllipsoid earth = new OneAxisEllipsoid(earthRadius,                     earthFlattening,rotatingEarthFrame);

// FORCE MODEL SETTINGS
// tesseral and zonal harmonics - gravity field
int degree = 10; 
int order  = 10;
// third bodies
List<String> thirdBodiesNames = new ArrayList();
thirdBodiesNames.add("MOON");
thirdBodiesNames.add("SUN");
// drag
double cd = 2.0;
double dragArea = 3.6;
Atmosphere atm = new DTM2000(new MarshallSolarActivityFutureEstimation(MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,                                                                  MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE),CelestialBodyFactory.getSun(), earth);
// srp
double cR = 1.2; 
double srpArea = 3.6;
        
// GRAVITATIONAL PARAMETER
UnnormalizedSphericalHarmonicsProvider unnormalized =  GravityFieldFactory.getUnnormalizedProvider(degree, order);
double mu = unnormalized.getMu();
         
        
// INITIAL MEAN ORBIT        
AbsoluteDate epoch = new AbsoluteDate(2026,1,1,0,0,0.0,timeScale);
double sma =   6926293.6178500970000000;
double eX = 0.0000123429246566; 
double eY = 0.0007713595726267;
double incl = FastMath.toRadians(49.0022185134797250); 
double raan = FastMath.toRadians(0.0026519771157681); 
double AOL = FastMath.toRadians(-0.0029642631972070);
        
// INITIAL MASS 
double mass = 289.0;

I display here the first test java code:

////////////////////////////////////////////////////////////////////////
/// INITIALIZATION        ///////////////////////////////////////////////////////////////////////
// INITIALIZATION SPACECRAFT STATE        
CircularOrbit meanCO = new CircularOrbit(sma, eX, eY, incl,raan,AOL,
            PositionAngleType.TRUE, pseudoInertialFrame,epoch,mu); 
SpacecraftState initialMeanState = new SpacecraftState(meanCO,mass);

// INITIALIZATION FORCE MODELS
List<DSSTForceModel> dsstForces = new ArrayList(); 
      
dsstForces.add(new DSSTZonal(unnormalized));
dsstForces.add(new DSSTTesseral(rotatingEarthFrame, earthAngularVelocity, unnormalized));
for(int ii=0; ii<thirdBodiesNames.size(); ii++){
            dsstForces.add(new DSSTThirdBody(CelestialBodyFactory.getBody(thirdBodiesNames.get(ii)), mu));
}
dsstForces.add(new DSSTAtmosphericDrag(atm, cd,dragArea, mu));
dsstForces.add(new DSSTSolarRadiationPressure(cR,srpArea,CelestialBodyFactory.getSun(), earth, mu));
        
// INITIALIZATION DSST PROPAGATOR
double minStep = 1e-3; 
double maxStep = 864000.0; 
double positionError = 1e-3; 
double positionScale = 1.0;    
DormandPrince853IntegratorBuilder integratorBuilder = 
                new DormandPrince853IntegratorBuilder(minStep,
                maxStep, positionError);
DSSTPropagatorBuilder dsstPropagatorBuilder = new DSSTPropagatorBuilder(meanCO,integratorBuilder,positionScale,PropagationType.MEAN,PropagationType.MEAN);
DSSTPropagator prop = dsstPropagatorBuilder.buildPropagator(dsstPropagatorBuilder.getSelectedNormalizedParameters());     
for(int i=0; i<dsstForces.size(); i++){
  prop.addForceModel(dsstForces.get(i));
}   
           ////////////////////////////////////////////////////////////////////////
/// PROPAGATION        ////////////////////////////////////////////////////////////////////////
// PROPAGATION 1
prop.setInitialState(initialMeanState,PropagationType.MEAN);
SpacecraftState st1 = prop.propagate(initialMeanState.getDate(),     initialMeanState.getDate().shiftedBy(3.0*Constants.JULIAN_DAY));
        
System.out.println("OUTCOMES PROPAGATION 1");
System.out.println("sma "+ st1.getA());
System.out.println("equinoctial ex "+ st1.getEquinoctialEx());
System.out.println("equinoctial ey "+ st1.getEquinoctialEy());
System.out.println("hx "+ st1.getHx());
System.out.println("hy "+ st1.getHy());
System.out.println("lv "+ st1.getLv());
      
       
        
// PROPAGATION 2        
prop.setInitialState(initialMeanState,PropagationType.MEAN);
SpacecraftState st2 = prop.propagate(initialMeanState.getDate(),
                initialMeanState.getDate().shiftedBy(3.0*Constants.JULIAN_DAY));
System.out.println("OUTCOMES PROPAGATION 2");
System.out.println("sma "+ st2.getA());
System.out.println("equinoctial ex "+ st2.getEquinoctialEx());
System.out.println("equinoctial ey "+ st2.getEquinoctialEy());
System.out.println("hx "+ st2.getHx());
System.out.println("hy "+ st2.getHy());
System.out.println("lv "+ st2.getLv());

The outcomes I obtained are:

OUTCOMES PROPAGATION 1
sma 6926182.174564137
equinoctial ex 2.049349492819461E-4
equinoctial ey 7.502247052818768E-4
hx 0.4408591973930216
hy -0.11555193292657281
lv 283.92149739293023
OUTCOMES PROPAGATION 2
sma 6926182.112382858
equinoctial ex 2.0492907020660124E-4
equinoctial ey 7.502325907695503E-4
hx 0.4408591964314755
hy -0.11555193612854477
lv 283.92150095673

I display here the second test java code:

////////////////////////////////////////////////////////////////////////
/// INITIALIZATION
///////////////////////////////////////////////////////////////////////

// INITIALIZATION SPACECRAFT STATE        
CircularOrbit meanCO = new CircularOrbit(sma, eX, eY, incl,raan,AOL,PositionAngleType.TRUE, pseudoInertialFrame,epoch,mu); 
SpacecraftState initialMeanState = new SpacecraftState(meanCO,mass);

// INITIALIZATION FORCE MODELS
List<DSSTForceModel> dsstForces1 = new ArrayList(); 
dsstForces1.add(new DSSTZonal(unnormalized));
dsstForces1.add(new DSSTTesseral(rotatingEarthFrame, earthAngularVelocity, unnormalized));
for(int ii=0; ii<thirdBodiesNames.size(); ii++){
        dsstForces1.add(new DSSTThirdBody(CelestialBodyFactory.getBody(thirdBodiesNames.get(ii)), mu));
}
dsstForces1.add(new DSSTAtmosphericDrag(atm, cd, dragArea, mu));
dsstForces1.add(new DSSTSolarRadiationPressure(cR,srpArea, CelestialBodyFactory.getSun(), earth, mu));
    
    
List<DSSTForceModel> dsstForces2 = new ArrayList(); 
dsstForces2.add(new DSSTZonal(unnormalized));
dsstForces2.add(new DSSTTesseral(rotatingEarthFrame, earthAngularVelocity,unnormalized));
for(int ii=0; ii<thirdBodiesNames.size(); ii++){
       dsstForces2.add(new DSSTThirdBody(CelestialBodyFactory.getBody(thirdBodiesNames.get(ii)), mu));
}
dsstForces2.add(new DSSTAtmosphericDrag(atm, cd,dragArea, mu));
dsstForces2.add(new DSSTSolarRadiationPressure(cR,srpArea, CelestialBodyFactory.getSun(), earth, mu));
    
// INITIALIZATION DSST PROPAGATOR
double minStep = 1e-3; 
double maxStep = 864000.0; 
double positionError = 1e-3; 
double positionScale = 1.0; 
    
DormandPrince853IntegratorBuilder integratorBuilder = new DormandPrince853IntegratorBuilder(minStep,maxStep, positionError);
DSSTPropagatorBuilder dsstPropagatorBuilder = new DSSTPropagatorBuilder(meanCO,integratorBuilder,positionScale,PropagationType.MEAN,PropagationType.MEAN);
DSSTPropagator prop = dsstPropagatorBuilder.buildPropagator(dsstPropagatorBuilder.getSelectedNormalizedParameters());        
    
    
    
////////////////////////////////////////////////////////////////////////
/// PROPAGATION
////////////////////////////////////////////////////////////////////////
// PROPAGATION 1
prop.setInitialState(initialMeanState,PropagationType.MEAN);
for(int i=0; i<dsstForces1.size(); i++){
   prop.addForceModel(dsstForces1.get(i));
}
SpacecraftState st1 = prop.propagate(initialMeanState.getDate(),initialMeanState.getDate().shiftedBy(3.0*Constants.JULIAN_DAY));
System.out.println("OUTCOMES PROPAGATION 1");
System.out.println("sma "+ st1.getA());
System.out.println("equinoctial ex "+ st1.getEquinoctialEx());
System.out.println("equinoctial ey "+ st1.getEquinoctialEy());
System.out.println("hx "+ st1.getHx());
System.out.println("hy "+ st1.getHy());
System.out.println("lv "+ st1.getLv());
    
prop.removeForceModels();

// PROPAGATION 2        
prop.setInitialState(initialMeanState,PropagationType.MEAN);
for(int i=0; i<dsstForces2.size(); i++){
    prop.addForceModel(dsstForces2.get(i));
}
SpacecraftState st2 = prop.propagate(initialMeanState.getDate(),initialMeanState.getDate().shiftedBy(3.0*Constants.JULIAN_DAY));
System.out.println("OUTCOMES PROPAGATION 2");
System.out.println("sma "+ st2.getA());
System.out.println("equinoctial ex "+ st2.getEquinoctialEx());
System.out.println("equinoctial ey "+ st2.getEquinoctialEy());
System.out.println("hx "+ st2.getHx());
System.out.println("hy "+ st2.getHy());
System.out.println("lv "+ st2.getLv());

The outcomes I obtained are:

OUTCOMES PROPAGATION 1
sma 6926182.174564137
equinoctial ex 2.049349492819461E-4
equinoctial ey 7.502247052818768E-4
hx 0.4408591973930216
hy -0.11555193292657281
lv 283.92149739293023
OUTCOMES PROPAGATION 2
sma 6926182.174564137
equinoctial ex 2.049349492819461E-4
equinoctial ey 7.502247052818768E-4
hx 0.4408591973930216
hy -0.11555193292657281
lv 283.92149739293023

I would like to know whether I’m doing something wrong in using the DSSTPropagator: is there maybe a way to “re-initializate” the DSSTForceModels? Or is it necessary to recreate them at each propagation?

Thank you in advance

Serrof · November 26, 2025, 12:34pm

Hi,

thanks for your analysis.
It is possible that some computations are cached or something and create this slight difference.
One would need to dig into the DSST code to be sure…
Have you tried using a subset of forces to see if its an orbital perturbation in particular that introduces the descrepancy?

Cheers,
Romain.

icav · December 3, 2025, 1:19pm

Hi Romain,

Thank you for your answer. I did the test you suggested. Apparently, the issue is caused by the introduction of the DSSTAtmosphericDrag. Do you have any idea why this occurs?

Thank you

Irene