Why the Eph generated by my orekit numeric propagator with drag model is not good?

Orekit usage
1

Hello everyone!
I use orekit numeric propagator to generate sat eph which about 66 hours away from the epoch of the element.
drag and solar pressure force model, third body attraction are included.
the DragCD I used is 2.1，DragSF is 9.7 m2

today I got some telemetry data of the target sat, there are gps information in it.

I compared the gps data and the eph，the distance between them is about 40KM, this is out of my estimate。

I tried some adjustment to improve the result. I found that when I removed the drag and solar pressure force model from propagator, the distance between the telemetry gps data and eph is
narrowed to 1.5KM

I tried STK HPOP to make the eph，the result is very similar to my orekit result, with drag model ,the result is not good ,without drag model,the result is very close to my telemetry gps data.

I’m very confused about that. How did this happen,Is there something I did wrong.

Hope some guy can help me to make clear of that.

Thanks!

p.s the target sat mass is about 1000kg，the DragSF is about 10m2. which are official data.

2

Hi @youngcle

Is it possible to have a sample of code? Especially to see how you configure your orbit propagator and how you initialize the force models.

Thank you,
Bryan

3

Comparing directly propagation against measurements works only if your dynamical model is really finely tuned. Drag is modelled using an atmosphere model whose evolution is difficult to predict.
How did you get the 2.1 drag coefficient? Does it come from an orbit determination from previous period data? How are the solar activity coefficients updated in the atmosphere model? Do you use an accurate modelling for the spacecraft shape, including solar arrays with an associated attitude model or a simpler sphere model?

4
  1. The drag coefficient “2.1” I used is told by the satellite provider. It’s a fixed value，so I do not believe it come from an orbit determination.
  2. The solar activity coefficients were fixed value. I made static fixed input params for nls
    both daily and average of F10.7 solar flux is set to 140,which I think is similar as STK force model configuration “manual input”
    Although the atm model setting is not perfect, I think the result is out of my expect.
  3. I have no accurate model about the satellite, so I used IsotropicDrag model and set the SF to 9.6 m2 which was told by satellite provider. In Stk Scene, I used sphere model too.

The both result made by my code and stk use the similar drag model are not good compare to telemetry gps data(66 hours away from the elements epoch,the difference is 80km, which cannot be used in my app), when I removed the drag model.Both of the result are very good(the difference is about 1km ). This confused me very much.

The following is the code I used to define the propagator.

        final double minStep_num = 0.001;
        final double maxstep_num = 20.0;
        final double positionTolerance = 1.0;
        final OrbitType propagationType = OrbitType.KEPLERIAN;
        final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit_loaded,
                propagationType);
        AdaptiveStepsizeIntegrator integrator_adpative = new DormandPrince853Integrator(minStep_num, maxstep_num,
                tolerances[0], tolerances[1]);

        // Propagator
        NumericalPropagator propagator_num = new NumericalPropagator(integrator_adpative);

        propagator_num.setOrbitType(propagationType);

        // Force Model (reduced to perturbing gravity field)
        // 3rd body (SUN)
        propagator_num.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
        // 3rd body (MOON)
        propagator_num.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));

        // ## Drag coefficient
        double DRAG_CD = drag_cd
        // Drag area (m^2)
        double DRAG_SF = drag_sf
        // Drag
        final double ae = normalized.getAe();

        final OneAxisEllipsoid earth_atm = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame);

       

        //here,I make a static fixed input params for nls
        //both daily and average of F10.7 solar flux  is set to 140,which I think is similar as STK forcemodel "manual input"
        StaticInputParams_NLS staticInputParams_NLS = new StaticInputParams_NLS();
        final Atmosphere atm =new NRLMSISE00(staticInputParams_NLS, CelestialBodyFactory.getSun(),
         earth_atm); //
        
        DragSensitive dragShap = new IsotropicDrag(DRAG_SF, DRAG_CD );

        propagator_num.addForceModel(new DragForce(atm, dragShap));

        
        // Solar Radiation Pressure
        propagator_num.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae,
                new IsotropicRadiationSingleCoefficient(SOLAR_RADIATION_PRESSURE_SF, SOLAR_RADIATION_PRESSURE_CR)));
        // ocean tide
        UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);

        OceanTides oceanTides = new OceanTides(earthFrame, normalized.getAe(), normalized.getMu(), true,
                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS, 50, 50, IERSConventions.IERS_2010, ut1);//
        propagator_num.addForceModel(oceanTides);

        // solid tide
        SolidTides solidTides = new SolidTides(earthFrame, normalized.getAe(), normalized.getMu(),
                normalized.getTideSystem(), IERSConventions.IERS_2010, ut1, CelestialBodyFactory.getSun(),
                CelestialBodyFactory.getMoon());
        propagator_num.addForceModel(solidTides);

        // relativity
        final Relativity forceModel = new Relativity(Constants.EIGEN5C_EARTH_MU);
        propagator_num.addForceModel(forceModel);

        ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(earthFrame, normalized);

        // Add force model to the propagator
        propagator_num.addForceModel(holmesFeatherstone);
        propagator_num.setInitialState(new SpacecraftState(orbit_loaded, mass));
1 Like
5

Hi @youngcle,

What are the dates you are propagating for ?
Maybe a flux of 140 is too high for these dates?

You can check the flux values in the CSSI Space Weather data file for example and set it accordingly to your dates.
A version of this file is available on the orekit-data repository here, and the latest version is here.

For example the latest prediction for the F10.7 flux for july 2020 is 68.
Note that with Orekit version 10.2 you will be able to use this file to set up your atmospheric model.

Maxime

6

Is there a big difference between your results and STK’s HPOP results

7

From your conclusion, I think the main reason is the Atmosphere model wrong used, maybe you should check the input parameter carefully, maybe you could put all of your code on the forum to CHECK.

8

Very big drag coefficient

9

Four years later, I got another problem, I’m here for some help.

I have a sat, and the propagation precision about 24 hours is very bed( and at same program the other’s is good, precision about 1km).

here’s my sample orbit. I use these to verify the propagator’s precision. these orbit is collect from the satellite onboard orbit determination unit, they are computed by onboard gps data, which could be trusted in precision about 20 meters.

epoch a e i raan mean anormaly pa
2024-10-10T08:09:27 6876309.074281642 0.0016968127588000805 97.40846087585038 358.27004563297686 71.18672546431036 81.52521973643408
2024-10-11T08:33:28 6873657.369565798 0.0012595511732739545 97.4116928917632 359.2760187502814 222.80674750963234 102.29252510119082
2024-10-12T08:59:21 6873506.9612276675 0.0016978517562921804 97.41501173669005 0.28179640709381515 76.67027998179596 68.19769199758626
2024-10-13T08:43:23 6878366.920891901 0.0016874493888859622 97.41458745784337 1.2569278168086473 73.29568242381583 92.18120178660028

I use one of the orbits as the initial state and propagate state to next orbit’s epoch,about 24 hours and caculate the euclidean distance of the two state.

the following is the result:

2024-10-11T08:33:28Z propagateTimeDis:24 pos distance:15285.566503

sample:{5,641,849.683862683; 437,784.2678928641; -3,913,115.329677807}
work:{5,633,256.782776439; 439,542.7488851628; -3,925,634.0492141857}

2024-10-12T08:59:21Z propagateTimeDis:24 pos distance:873.502079

sample:{-5,629,636.735888331; -535,567.581005029; 3,902,400.9011060437}
work:{-5,630,120.063441826; -535,479.0984061974; 3,901,678.7022498716}

2024-10-13T08:43:23Z propagateTimeDis:23 pos distance:9432.690144

sample:{-6,654,486.504124772; -365,766.8072691014; 1,688,288.6726203945}
work:{-6,656,614.755638698; -364,637.5940812194; 1,679,168.8553841766}

the orekit numeric propagator running with new cssi space weather data at 20241101

and the drag model is solar array and box

drag SF is 8 m2

the satellite is an earth observe satellite at 500KM solar synchronous orbit, it’s attitude aim to sun when go above latitude ±60, and
aim to earth between ±lat 60.

the result is too bad for my application.

Any ideas why these happens and how to fix it?

p.s. with STK HPOP, same data, the propagate precision is about 4KM

which also not good.

the expected precision should within 1km.

at the same time,some other satellite’s result is good enough.

10

Instead of using a single osculating points from the on-board GNSS device, I would suggest to get a range covering about one orbit or more and fit it with Orekit orbit determination. I will help ensure orbits are closer even despite probably different dynamical model.

11

thanks for reply

These samples are already the orbit determine results which computed onboard the satellite, use kalman filter algorithm

I found some TLE orbits, datetime close to the samples to verify the reliablity. The difference between the sample and tle results is very close,about 200meters

So I believe the sample orbit is reliable.

I doubt that the high solar activity and the attitude policy cause the result

But, not all the satellites in my system got the bad precision.

I’m very confused and have no idea to fix these.

12

I understand, but what I was telling is that the dynamical models used on-board and in your own orbit determination may be different. Using one reference point only is often a cause of problems. This is in fact the same reason why when converting between TLE and anything else we recommend to fit over an arc.
When you fit, you may “smooth out” some differences between dynamical models, this will mainly appear in linear and quadratic trends along the orbit by slightly adjusting semi major axis and drag coefficient, and in short periodic terms in just about all orbital parameters.

It seems error is not steadily augmenting, but has some harmonics (the point for 2024-10-12 is much closer than the two other points). One thing you could try is plot the error with a much finer time step and look if it exhibits some periodic effect. You could also display this error not has a 3D distance but has a signed offset in some local orbital frame (radial, tangential, normal). Such plots may help identify the culprit. For example if you have something wrong in the orbital plane orientation, you will get an harmonic error at orbital period in the normal component. If you have something wrong in 3rd body modelling you would get a drift in ascending node plus an harmonic at half orbital period. If you have something wrong with drag, it would show up as a quadratic buildup in the tangential direction.

13

thanks.

I believe it is very worth to do.

One thing you could try is plot the error with a much finer time step and look if it exhibits some periodic effect. You could also display this error not has a 3D distance but has a signed offset in some local orbital frame (radial, tangential, normal). Such plots may help identify the culprit. For example if you have something wrong in the orbital plane orientation, you will get an harmonic error at orbital period in the normal component. If you have something wrong in 3rd body modelling you would get a drift in ascending node plus an harmonic at half orbital period. If you have something wrong with drag, it would show up as a quadratic buildup in the tangential direction.

If I got some progress I will post here.