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