Problem with Numerical Propagator Precision

#1

Hi, all,
I am trying to use NumericalPropagator to propagate Envisat’s trajectory for 24 hours, but compared with its CPF ephemeris, the errors of the result are really large. There must be something wrong in my code, but I don’t know where.
The initial Date is 2018-10-25 00:00:00.000. I added OceanTides, solid tides, third ocean attraction, solar radiation pressure, and Attitude modes.
Since Envisat is about 780km high, 8211kg, with dimensions 26m10m5m, I figured several parameters used in the code: cd = 2.2, area = 32.0 ( in fact the area = 2610, but as the default mass is 1000kg, I used area = 2610/8.211), cr =1.88
Here is my code:

public class MasterMode {

/** Program entry point.
 * @param args program arguments (unused here)
 */
public static void main(String[] args) {
    try {

        // configure Orekit
        File home       = new File(System.getProperty("user.home"));
        File orekitData = new File(home, "orekit-data");
        if (!orekitData.exists()) {
            System.err.format(Locale.US, "Failed to find %s folder%n",
                              orekitData.getAbsolutePath());
            System.err.format(Locale.US, "You need to download %s from %s, unzip it in %s and rename it 'orekit-data' for this tutorial to work%n",
                              "orekit-data-master.zip", "https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip",
                              home.getAbsolutePath());
            System.exit(1);
        }
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.addProvider(new DirectoryCrawler(orekitData));
        
        // gravitation coefficient
        double mu =  3.986004415e+14;

        // inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();

        // Initial date
        AbsoluteDate initialDate = new AbsoluteDate(2016, 10, 25, 0, 0, 0.000,
        											TimeScalesFactory.getUTC());


        
        Vector3D p= new Vector3D(-3963665.96200000,	3057637.69600000,	5089564.99900000);
        Vector3D v= new Vector3D(-2974.79114922258,	4691.50122498370,	-5122.50267471032);
        final double r2 = p.getNormSq();
        final double r3 = r2 * FastMath.sqrt(r2);
        
        PVCoordinates pv = new PVCoordinates(p,v);
        System.out.println("current dir:" + System.getProperty("user.dir"));
        Frame ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        
        Frame j2000 = FramesFactory.getEME2000();
        PVCoordinates pv1 = ITRF.getTransformTo(j2000, initialDate).transformPVCoordinates(pv);
        
        TimeStampedPVCoordinates pv2 = new TimeStampedPVCoordinates(initialDate, pv1.getPosition(), pv1.getVelocity(), new Vector3D(-mu / r3, pv1.getPosition()));
        System.out.println(" orbit0 j2000: " + pv1);
        
        Orbit initialOrbit = new CartesianOrbit(pv1, j2000, initialDate, mu);

// // Initial state definition
SpacecraftState initialState = new SpacecraftState(initialOrbit);

        // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
        final double minStep = 0.001;
        final double maxstep = 100.0;
        final double positionTolerance = 5.0;
        
        final OrbitType propagationType = OrbitType.CARTESIAN;
        final double[][] tolerances =
                NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
        AdaptiveStepsizeIntegrator integrator =
                new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);

        

        // Propagator
        NumericalPropagator propagator = new NumericalPropagator(integrator);
        propagator.setOrbitType(propagationType);

        // Force Model (reduced to perturbing gravity field)
        final NormalizedSphericalHarmonicsProvider provider = 
                 GravityFieldFactory.getNormalizedProvider(200, 200);
        		         		
        ForceModel holmesFeatherstone =
                new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,
                                                                            true),
                                                      provider);

        // Add force model to the propagator
		propagator.addForceModel(holmesFeatherstone);
        
     // ocean tides force model
        final Frame bodyFrame;
        bodyFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true); 

        final double equatorialRadius;

        equatorialRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;

        final double flattening;
        flattening = Constants.WGS84_EARTH_FLATTENING;


        propagator.addForceModel(new OceanTides((new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getBodyFrame(),
        		provider.getAe(), provider.getMu(),
        		1, 1, IERSConventions.IERS_2010,
        		TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true)));

        
        
        
     // solid tides force model
        List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
        
        solidTidesBodies.add(CelestialBodyFactory.getSun());


        solidTidesBodies.add(CelestialBodyFactory.getMoon());
        

        solidTidesBodies.add( CelestialBodyFactory.getJupiter());

//
// solidTidesBodies.add( CelestialBodyFactory.getMars());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getMercury());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getNeptune());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getPluto());
//
//
// solidTidesBodies.add( CelestialBodyFactory.getSaturn());
solidTidesBodies.add( CelestialBodyFactory.getVenus());

        propagator.addForceModel(new SolidTides((new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getBodyFrame(),
        		provider.getAe(), provider.getMu(),
        		provider.getTideSystem(), IERSConventions.IERS_2010,
        		TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true),
        		solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
        

        // third body attraction
        
        propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
        propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));

        
     // drag
            final double  cd          = 2.2;
            final double  area        = 32.0;
            final boolean cdEstimated = true;

            MarshallSolarActivityFutureEstimation msafe =
                            new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)",
            MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);

// DataProvidersManager manager = DataProvidersManager.getInstance();
manager.feed(msafe.getSupportedNames(), msafe);
Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame));
propagator.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));

        // solar radiation pressure
              
            final double  cr          = 1.88;
            propagator.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(),
                                                                       (new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame)).getEquatorialRadius(),
                                                                       new IsotropicRadiationSingleCoefficient(area, cr)));
           
        
        
        
        
         // post-Newtonian correction force due to general relativity
            
          propagator.addForceModel(new Relativity(provider.getMu()));
            

            
           
          
          final AttitudeMode mode;
              mode = AttitudeMode.NADIR_POINTING_WITH_YAW_COMPENSATION;
            
        // Set up initial state in the propagator
        propagator.setInitialState(initialState);

        // Set up operating mode for the propagator as master mode
        // with fixed step and specialized step handler
        propagator.setMasterMode(60., new TutorialStepHandler());
        
        // Extrapolate from the initial to the final date
        SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(86400.));
        KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
        System.out.format(Locale.US, "Final state:%n%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
                          finalState.getDate(),
                          o.getA(), o.getE(),
                          FastMath.toDegrees(o.getI()),
                          FastMath.toDegrees(o.getPerigeeArgument()),
                          FastMath.toDegrees(o.getRightAscensionOfAscendingNode()),
                          FastMath.toDegrees(o.getTrueAnomaly()));

    } catch (OrekitException oe) {
        System.err.println(oe.getMessage());
    }
}

/** Specialized step handler.
 * <p>This class extends the step handler in order to print on the output stream at the given step.<p>
 * @author Pascal Parraud
 */
private static class TutorialStepHandler implements OrekitFixedStepHandler {

	private PrintStream outputStream;
	private final File outputFile = new File(new File(System.getProperty("user.home")), "result047.txt");
	
    private TutorialStepHandler() {
        //private constructor
    }

    public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
    	try{
    		outputStream        = new PrintStream(outputFile, "UTF-8");
    	}catch (IOException ioe) {
            throw new OrekitException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE, ioe.getLocalizedMessage());
        }
        System.out.println("          date                a           e" +
                           "           i         \u03c9          \u03a9" +
                           "          \u03bd");
    }

    public void handleStep(SpacecraftState currentState, boolean isLast) {
        Frame ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        AbsoluteDate initialDate = currentState.getDate();
        
        
    		outputStream.format(Locale.US, "%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n",
                currentState.getDate(),
                o.getPVCoordinates(initialDate, ITRF).getPosition().getX(), 
                o.getPVCoordinates(initialDate, ITRF).getPosition().getY(),
                o.getPVCoordinates(initialDate, ITRF).getPosition().getZ(),
                o.getPVCoordinates(initialDate, ITRF).getVelocity().getX(), 
                o.getPVCoordinates(initialDate, ITRF).getVelocity().getY(),
                o.getPVCoordinates(initialDate, ITRF).getVelocity().getZ());
        if (isLast) {
            System.out.println("this was the last step ");
            System.out.println();
        }
    }

}

/** Attitude modes. */
private static enum AttitudeMode {
NADIR_POINTING_WITH_YAW_COMPENSATION() {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new YawCompensation(inertialFrame, new NadirPointing(inertialFrame, body));
}
},
CENTER_POINTING_WITH_YAW_STEERING {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new YawSteering(inertialFrame,
new BodyCenterPointing(inertialFrame, body),
CelestialBodyFactory.getSun(),
Vector3D.PLUS_I);
}
},
LOF_ALIGNED_LVLH {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.LVLH);
}
},
LOF_ALIGNED_QSW {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.QSW);
}
},
LOF_ALIGNED_TNW {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.TNW);
}
},
LOF_ALIGNED_VNC {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.VNC);
}
},
LOF_ALIGNED_VVLH {
public AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
{
return new LofOffset(inertialFrame, LOFType.VVLH);
}
};

    public abstract AttitudeProvider getProvider(final Frame inertialFrame, final OneAxisEllipsoid body)
       ;

}

}
The results written in the “result_047.txt” are in the ITRF frame. As Envisat is a SLR satellite, I compared the results with its CPF Ephemeris.The errors of x,y,z(ITRF) can reach up to almost 160m after 24 hrs. I am now trying to decrease the errors, but it seems there is nothing more I can do. Is there any perturbation I missed? Shall I change my propagator? Can anyone give me some advices?

By the way, I tried another satellite named Ajisai(altitude 1488km, weight 685 kg ), the errors of x,y,z(ITRF) can reach up to almost 20km after 24 hrs. I used parameters cd = 2.2, area = 3.2*0.685 ( it is a sphere with a diameter of 2.15m ), cr =1.88. The errors are really large……… I am clearly getting something wrong. I would be very grateful for any hints provided!

#2

Hi,

I have a few remarks about your model:

  • 8211kg was the launch mass, Envisat’s current mass is slightly lower because some of the hydrazine was burned (though I don’t know exactly how much)
  • You should set properly the spacecraft mass in the propagator
  • The real cross section is lower than 260m², because the spacecraft is not 10m wide everywhere. And because it is tumbling, the cross section is variable. There are several papers out there which analyzed Envisat’s spin properties using laser ranging: https://ieeexplore.ieee.org/abstract/document/6809997
  • Instead of the DTM2000 atmosphere model, you could give it a try with NRLMSISE00
  • Is your EOP data up-to-date? For 2018-10-25, I guess it is, but I am just asking
  • You used the getITRF method with the simpleEOP=true argument, i.e. tidal effects are ignored when interpolating EOP: https://www.orekit.org/static/apidocs/org/orekit/frames/FramesFactory.html#getITRF(org.orekit.utils.IERSConventions,boolean)

About the results themselves, a difference of a hundred meters after 24 hours is not shocking if the propagation model used for the orbit determination behind the CPF is different than yours.
I have been doing a little bit of orbit determination using SLR data (https://nbviewer.jupyter.org/github/GorgiAstro/laser-orbit-determination/blob/6cafcef83dbc03a61d64417d0aeb0977caf0e064/02-orbit-determination-example.ipynb), and noticed that I get a difference of several tens of meters after 24 hours (see second-to-last plot in my notebook), it depends on the satellite.
I would like to know which model is used by the different CPF providers to perform their OD. Apart from the mass and cross section, the cd and cr coefficients are hard to properly estimate.

As for Ajisai, 20km difference is really a lot.

Cheers
yzokras

#3

Hi @miaomiao ,

I agree with @yzokras, 160m of error for a 24hours propagation is not that big.
It translates to an error of only ~20 milliseconds on the orbit track.

20km on the other hand is indeed a lot. You probably have something wrong in your setup.

For Envisat you could maybe get closer results by using a better modelling of the solar activity.
MSAFE (Marshall Solar Activity Future Estimation) files will give you monthly averaged data.
You can find daily data on the internet that will help you improve your atmosphere model.
Here is an example from AGI ftp. The problem is that you will have to implement your own parser and your own versions of the interfaces NRLMSISE00InputParameters or DTM2000InputParameters since they are not available in Orekit yet.

Cheers,
Maxime

#6

Hi, @MaximeJ @yzokras,
Would you please give me some more hints? I tried to reconfigure my parameters, but the error still reached up to 12km. Do I need to change my GravityField model? or the Atomsphere drag model? Thanks a lot!

#7

Hi @miaomiao,

You are talking about your “Ajisai” test right ?

If you downloaded the latest orekit-data on Orekit website then I think your gravity field should be good enough.

At this altitude (1488km) the drag should have little effect. I’m not even sure you need to plug the model.

Did you follow @yzokras recommendations to improve your accuracy ?

I may be mistaken but it seems you did add not the Earth solid tides in your force models. You should probably do it.
Also the ocean tides order/degree is only (1,1). You could increase that, say to (7,7), but I’m not enough of an expert to tell you the increase in precision you’d get from that.
To speed up your propagation I advice you not to use the yaw steering attitude provider. It is a complex and computation intensive model. Asijai is a “ball-like” satellite so there is no need to give an accurate attitude provider here. The default provider (aligned with GCRF frame) should be enough.

Did you get these parameters from the CPF ephemeris website ?

Now if you are still stuck with your problem you could answer this post and join the latest version of your program and the CPF ephemeris you’re trying to follow so we can try to reproduce your results and investigate a little more.

Cheers,
Maxime

#8

Hi, @MaximeJ
Thank you very much for your reply! I followed your instructions, but it doesn’t work.
The error of position is still so high as almost 12km after 24 hrs.
I attached my code and the corresponding cpf ephemeris here. The initial time is 2018-12-31,00:00:00 with the time Interval is 60s. And the frame is ECEF.
In the model, I added gravity, ocean tides, solid tides, third body attraction, solar radiation pressure and post-Newtonian correction force due to general relativity. Is there anything I missed?
Looking forward to your reply! I really appreciate it!

MasterMode.java (17.8 KB)
cpf_ajisai.txt (157.6 KB)

#9

Hi all,

With your error of 12km maybe the first thing to do is simplify a little bit the problem.

  1. As @MaximeJ said, you can first remove the drag force. At ~ 1500 km of altitude, this effect is negligible.
  2. Such as the drag force, you can remove the relativity force model. Its impact is low compared to the other force models.
  3. Reduce the degree and the order of the gravity field: 20 * 20 is sufficient.

Moreover, when you use getITRF method to build the initial pv coordinates you put simpleEOP argument to false whereas when you call this method inside the StepHandler to get the current pv coordinates this parameter is equal to true. Is it normal ?

For the solar radiation pressure, I agree with the value of 1.0 for the reflection coefficient. However, for the cross section the value is about 3.6 m^2 (~ pi * (2.15 / 2)^2) not 6.6.

In order to reduce the error in position maybe you can investigate the initialization of the integrator. For instance, play with the values of minStep, maxStep and positionTolerance. Or try with another integrator.
With minStep = 0.001; maxStep = 1000.0 and positionTolerance = 10.0 I have a position error of 3km.

Best regards,
Bryan

1 Like
#10

Yeah,I get your point!It seems the most important thing is the steps and tolerance of the propagator. Thank you very much!
Best regards,
miaomiao

1 Like