Comparison to SP3 precision ephemeris

Orekit usage
1

Hey folks, I’ve been trying to build an Orekit propagator that closely matches the SLR precision ephemeris pulled from CDDIS for the Larets spacecraft (< 20m error per day of propagation). My propagator as it currently stands is only able to get to < 10km of error per day of propagation (see below).

image
image1598×732 64.2 KB

I loaded in the SLR precision ephemeris using the SP3 parsing / propagation framework like this:

// Parsing the SP3 file with the intent to keep in GCRF
        Path sp3File = Paths.get(this.getClass().getClassLoader().getResource(sp3DataFile).toURI());
        SP3Parser parser = new SP3Parser();
        SP3 file = parser.parse(
                new DataSource(sp3File.toString(), () -> getClass().getResourceAsStream(sp3DataFile)));

        // Pulling out the ephemeris for the Larets spacecraft
        SP3Ephemeris ephem = file.getSatellites().get(sp3SatName);

        // Converting the ephemeris into a propagator which will interpolate
        Propagator sp3Prop = ephem.getPropagator();

And then essentially built my Orekit propagator and compared the resulting propagation over the course of 24 hours. Since the documentation doesn’t cover more complex cases, it’s difficult to know if I’m setting up the software correctly or if I’m missing a small detail that may be causing this. The propagation setup I’ve built looks like this:

public Propagator getNumPropagator(TimeStampedPVCoordinates pvCoords, Frame frame, int gravDegree, int gravOrder,
            double coefAbsorb, double coefReflect, double crossSection, double dragCoef, double mass) {

        // Create orbit from input PV coords
        Orbit orbit = new CartesianOrbit(pvCoords, frame, OrekitUtils.MU);
        SpacecraftState initialState = new SpacecraftState(orbit, mass);

        // Create adaptive step integrator
        double minStep = 0.001;
        double maxstep = 1000.0;
        double positionTolerance = 1.0e-4;
        OrbitType propagationType = OrbitType.CARTESIAN;
        double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, propagationType);
        AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0],
                tolerances[1]);
        NumericalPropagator prop = new NumericalPropagator(integrator);

        // Create Gravity Model
        NormalizedSphericalHarmonicsProvider gravProvider = GravityFieldFactory.getNormalizedProvider(gravDegree,
                gravOrder);
        ForceModel gravModel = new HolmesFeatherstoneAttractionModel(frame, gravProvider);
        prop.addForceModel(gravModel);

        // Creating Lunar and Solar gravitational effects
        ThirdBodyAttraction lunarGrav = new ThirdBodyAttraction(CelestialBodyFactory.getMoon());
        ThirdBodyAttraction solarGrav = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
        prop.addForceModel(lunarGrav);
        prop.addForceModel(solarGrav);

        // Solid Earth Tides
        SolidTides tides = new SolidTides(frame,
                Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU,
                TideSystem.ZERO_TIDE,
                IERSConventions.IERS_2010,
                TimeScalesFactory.getUT1(IERSConventions.IERS_2010, false),
                CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
        prop.addForceModel(tides);

        // Relativity
        Relativity rel = new Relativity(Constants.EIGEN5C_EARTH_MU);
        prop.addForceModel(rel);

        // Solar Radiation Pressure Model
        if (!Double.isNaN(coefAbsorb) && !Double.isNaN(coefReflect)) {
            RadiationSensitive satSRPConfig = new IsotropicRadiationClassicalConvention(crossSection, coefAbsorb,
                    coefReflect);
            SolarRadiationPressure srp = new SolarRadiationPressure(CelestialBodyFactory.getSun(),
                    Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, satSRPConfig);
            srp.addOccultingBody(CelestialBodyFactory.getMoon(),
                    Constants.MOON_EQUATORIAL_RADIUS);
            prop.addForceModel(srp);
        }

        // Drag Model
        if (!Double.isNaN(dragCoef)) {
            BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                    Constants.WGS84_EARTH_FLATTENING, FramesFactory.getGCRF());
            CssiSpaceWeatherData cssi = new CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt");
            MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation(
                    MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
                    MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
            DTM2000 dragModel = new DTM2000(cssi, CelestialBodyFactory.getSun(), earth);
            // NRLMSISE00 dragModel = new NRLMSISE00(msafe, CelestialBodyFactory.getSun(), earth);
            DragSensitive satDragConfig = new IsotropicDrag(crossSection, dragCoef);
            DragForce drag = new DragForce(dragModel, satDragConfig);
            prop.addForceModel(drag);
        }

        prop.setInitialState(initialState);
        return prop;
    }

Is there something I’m configuring incorrectly, or are there better models / practices for a more accurate propagator? Any help you can provide would be appreciated!

Nick

2

One more thing… I verified the SP3 was loaded and being used properly by comparing to an alternate propagator which was able to keep the errors within 15m over the course of a day (and could probably be made better with some force model tweaks).

3

Could you try adding propagator.setOrbitType(propagationType) near the beginning of the propagator configuration?
By default the numerical propagator uses equinoxial parameters, so in your case I guess there is an inconsistency between the tolerances arrays and the state vector.

4

Another point (but it will be far below kilometers differences), rather than using a hard-coded TideSystem.ZERO_TIDE when you build the solid tide, you should use the value returned by gravModel.getTideSystem(), so you are sure both force models are consistent with each other.

5

I added this, and it didn’t end up making a large difference, but it is a good suggestion that I’ll keep in the code:

((HolmesFeatherstoneAttractionModel)gravModel).getTideSystem()
6

I added this and it didn’t make a noticeable difference unfortunately:

prop.setOrbitType(OrbitType.CARTESIAN);

Just so you can see what inputs I gave to the propagator method, here’s the call:

Propagator orekitProp = utils.getNumPropagator(initCoords, frame, 21, 21,
            Double.NaN, Double.NaN, 0.01, 2.2, 21);
7

Could you also try replacing OrekitUtils.MU with the mu from the gravity field?
What are the degree and order of your gravity field?

8

When you build earth model for the drag, use ITRF rather than GCRF for the body frame.

9

Ok modified to look like this, and there’s still a ~9.5km error after a day, although it did make about .5km difference, so were moving in the right direction:

        // Create Gravity Model
        NormalizedSphericalHarmonicsProvider gravProvider = GravityFieldFactory.getNormalizedProvider(gravDegree,
                gravOrder);
        ForceModel gravModel = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, false), gravProvider);

        // Create orbit from input PV coords
        Orbit orbit = new CartesianOrbit(pvCoords, frame, gravProvider.getMu());
        SpacecraftState initialState = new SpacecraftState(orbit, mass);
10

Oh, I think I get it. Here too you should use ITRF rather than GCRF. Your Earth (and its gravity field) is not rotating at all!

11

i just realized after I pasted it that I had put the wrong frame in there and then fixed it. So now I need to find what you meant originally

12

Ok I fixed the Earth in the drag model as well, and there wasn’t a significant change (< 100m):

BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                    Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, false));
13

Good news! That was the issue, but I didn’t save before the last run. That said, here’s the updated results! Seems pretty stable at around 18 meters max without much of an increasing trend. Thank you so much for your help!

image
image1598×732 81.8 KB

It turns out a rotating gravitational field is important!

14

This is a relief!

:rofl:

1 Like
15

Hi @nsabey

Does this simulation include solar radiation pressure? (I ask the question because I see in your code that there is a “if” condition)
Did you try adding OceanTides? It can maybe improve the accuracy by about 0.5-1 meter.

Bryan

16

Also, if you want to play with the limit of orbit propagation you can add KnockeRediffusedForceModel for Earth’s albedo and IR, and also LenseThirringRelativity and DeSitterRelativity for 2nd order relativity terms.

I’m not sure it will improve a lot the accuracy, but I think those models are used to generate sp3 files.

1 Like
17

great suggestions! Yes it does include solar radiation pressure in this case. I’ll try implementing your suggestions and see where it leaves the propagation.

18

To help you in your quest, please find below two documents presenting the propagator configuration used by GFZ and NSGF to generate their sp3 files from SLR orbit determination (see “ORBIT MODEL” section)

nsgf.dsc (13.0 KB)

gfz.dsc (12.2 KB)

19

Thank you, I’ll use these as a comparison to what I’ve gotten together!

20

After tweaking the force modeling a bit, I was able to get 2 days of propagation within 35 m of error, which is well within our needs. Thanks for all of your suggestions and help getting here!

image
image1862×1176 144 KB