Propagation with certain earth gravity field

Hi all.
I’m new in using orekit. I’ve a general question:
How can I propagate mean and osculating elements, with mathematical model that takes into account the gravity potential up to a certain spherical harmonics order and degree?

Thanks in advice.

Hi @Bryan welcome,

You can take a look at the tutorials page.

Hi @luc, and thanks for your answer. I’m trying to use the numerical propagator. I want to retrieve position/velocity in time (between start and final date), and in different frames.
I’m trying to use EphemerisGenerator. Is what I’ve done correct? (I attached the python code lines below).

Thanks in advance,


#Integrator options
minStep = 0.001 # [s]
maxstep = 100.0 # [s]
initStep = 1.0 # [s]

positionTolerance = 0.001 # [m]

propagationType = OrbitType.KEPLERIAN
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType)

#The actual integrator, in this case DormandPrince853, is part of the Hipparchos library.
#Note that the tolerances needs casting in Python to an array of doubles (floats).

integrator = DormandPrince853Integrator(minStep, maxstep, JArray_double.cast_(tolerances[0]), JArray_double.cast_(tolerances[1]))
# Double array of doubles needs to be casted in Python


satellite_mass = 100.0 # The models need a spacecraft mass, unit kg.
initialState = SpacecraftState(initialOrbit, satellite_mass)

propagator = NumericalPropagator(integrator)

#For the propagator to make sense it needs some forces acting on the satellite.
#Here we are adding a gravity field model.For a more detailed propagation, other force models can be added.

gravityProvider = GravityFieldFactory.getNormalizedProvider(4, 4)
propagator.addForceModel(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, False), gravityProvider))

#Set the propagator to ephemeris mode
generator = propagator.getEphemerisGenerator()

end_state = propagator.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))

ephemeris = generator.getGeneratedEphemeris()

intermediateState = ephemeris.propagate(initialDate.shiftedBy(3600.0 * 4))

pv_ECEF = intermediateState.getPVCoordinates(ECI_J2000)

pv_ECI = intermediateState.getPVCoordinates(ITRF)

I add another point, what does the NumericalPropagator expect to receive as input? Osculating elements?
So, if I want to propagate mean elements I can’t give them directly to the NumericalPropagator, but I’ve to convert them to osculating?


From a quick look, this code seems correct to me.
Note however that if you want to retrieve positions at several different times, you can do this on the fly rather than building an ephemeris and then picking up the points from it. You have to set up a step handler (most certainly a fixed step handler as it is more straightforward to use) that will be called during the propagation. In this step handler, you can do whatever you want with the current state, including extracting the embedded position in several frames.

You are right about the numerical propagator, it must be initialized with osculating parameters. Beware that in order to convert mean parameters into osculating parameters, you have to know the theory against which mean parameters have been built. Mean parameters with respect to Eckstein-Hechler are not the same as mean parameters against Brouwer-Lyddane, which are not the same as mean parameters with respect to DSST, which are not the same as mean parameters with respect to SGP4/SDP4, which are not the same as mean parameters obtained by simple numerical smoothing, which are not the same as, … well any other mean parameters you can dream of.

Thanks a lot @luc for your clarification. So if I want to propagate directly mean elements is not possible to do?

Moreover, I’ve a problem in the code I’ve shown above.

I’ve propagated osculating parameters ( of Sentinel-1A but when I plot the ground-track the starting point and the point after 12 days (repeating period) are far away from each other. I’ve tried to play with tolerance and steps, but nothing has changed.

Where could the problem be? In the precision of the elements? Is there another better site to take them?

Thanks in advance

You can propagator in mean elements using for example DSST propagator. This propagator allows you to specify if your initial elements are osculating or mean, and to specify if you want osculating or mean as output. Beware that if you propagate in mean elements, you cannot compare positions with another run performed in osculating elements, you will see several kilometers differences due to short periodic terms (but these differences should remain centered).

Your force model is insufficient. A 4x4 gravity field is much too low, I would use at least 60x60. You also don’t take into account drag and 3rd body attraction from Moon and Sun, which are not negligible.

I agree that the model I’ve used is not so accurate. But the displacement I see is about 80° for Latitude and 50° for Longitude, so I don’t think that is an accuracy problem (infact if I set 60,60 harmonics the result is still so bad)

What does adding drag change?
Is your wrong position still on the right orbit (i.e. somewhere in the same plane)? You can check that using the orbital momenum (with getPVCoordinates().getMomentum()). If you are still in the same plane, the difference is along the track in inertial frame (not in Earth frame), and this would point to a semi-major axis error. A semi-major axis error comes either from a mean/osculating conversion not performed properly (because the short period terms on semi-major axis are huge) or a drag modelling problem.

These elements do not look as osculating elements to me. They look as extracted from TLE (mainly because of the lack of semi major axis replaced by a mean revolution per day).

So, maybe the wrong result is due to wrong osculating elements. Could you suggest me a site where to get mean/osculating elements of artificial satellites? I’ll change initial conditions and see what happens

I’m greatefull for your availability

You can get the TLE from spacetrack, then you would need the TLEPropagator to use them directly.
For sentinel 3 and sentinel 6, you can get precise orbit in SP3 format from CDDIS.

Ok, with TLE I had no problem. I need mean and osculating to propagate them. If I will not find them, If I’m not wrong, there is a converter in orekit (from TLE), right?

Converting from/to TLE is not straightforward. You can look at TLE.stateToTLE which uses a TleGenerationAlgorithm, with two implementations provided by Orekit. For the other way round (from TLE to state), you have to use TLEPropagator.

When I say use a TLEPropagator, I do not mean just use it to get one osculating point, I mean use it in a PropagatorConverter that will use several points spread over at least one orbit. Converting on a single point will fail.

So, to get my initial osculating state to give as input to numerical propagator I have to propagate with TLEPropagator for at least one orbit, and then with PropagatorConverter I can get osculating elements of the initial state?

In fact the PropagatorConvertor will give you a completely configured numerical propagator that you can use immediately.

Maybe I’m wrong but, if from TLEPropagator I extract position/velocity at a certain time, from these I can compute osculating elements at that time. I can’t do in this way?

No, it would be correct at this date, but as soon as you propagate it will diverge, and it will diverge fast.
Converting TLE must do done over a sufficient time range because TLE propagation model (i.e. SGP4/SDP4) has its own set of perturbations, modeled by its own methods, so some effects are present in the models and some are not: when you compare it to any other model, you will see differences in both short periods and harmonics, so using only one point and matching two models on that single point never works.

Hi @luc, can I define initial conditions to give as input to NumericalPropagator as (x,y,z,vx,vy,vz) in ECEF, using SpacecraftState?