Propagation with Holmes-Featherstone Attraction Model is too different than reality

Hi everyone !

I’m trying to estimate the errors caused by my numerical propagator and came up with the following method :
I recover the orbital elements from a TLE at a certain date, say the 1st of March, and propagate to another date, for example the 9th of March. I then compare the results of my propagation with the updated TLE of this same satellite at the exact date of the results.
I also compare the results of my propagator with and without the added ForceModels.
However, my results seem way off. I compute the relative errors on the position with \frac{abs(x_{th} - x_{exp})}{ x_{th}} and for some reason I have more error when propagating with the “realistic” ForceModels than without.

I brute forced an investigation by trying all the combinations of forces I had, and discovered that the HolmesFeatherstoneAttractionModel is the main cause of the problem. When i remove it from the propagator, my results are way closer to the updated TLE (the reference), and i don’t understand why or how…

some pseudo code to illustrate my thoughts :

# define all the characteristics of the spacecraft
satellite_mass = ...

# I know the following part is not a good thing to do, but i don't know how to do it properly
# get the TLE at the initial date, and extract the orbital parameters from it
sat_initial_tle = TLE(line1,line2)
sat_initial_orbit = user_defined_function_that_creates_a_KeplerianOrbit(sat_initial_TLE)
# get the TLE at the final date, the reference
sat_reference_tle = TLE(line1,line2)
sat_reference_orbit = user_defined_function_that_creates_a_KeplerianOrbit(sat_reference_TLE)

# Create all the ForceModels that are added to the numerical propagator
isotropic_drag = IsotropicDrag(cross_section, drag_coef)
list_of_forces.append(DragForce(SimpleExponentialAtmosphere(earth, 0.0004, 42000.0, 7500.0), isotropic_drag))
gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
list_of_forces.append(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))
list_of_forces.append(ThirdBodyAttraction(moon))
list_of_forces.append(ThirdBodyAttraction(sun))
list_of_forces.append(SolarRadiationPressure(sun, Constants.WGS84_EARTH_EQUATORIAL_RADIUS, IsotropicRadiationClassicalConvention(cross_section, absorption_coef, reflection_coef)))

# create a NumericalPropagator like in the examples
integrator = DormandPrince853Integrator(...)
propagator_num = NumericalPropagator(integrator)
for force in list_of_forces:
    propagator_num.addForceModel(force)

# propagate
final_state = propagator_num.propagate(start_date.shiftedBy(duration_between_both_TLE))

# investigate the results
pos_z_theory = sat_reference_orbit.getPVCoordinates().getPosition().z
pos_z_experimental = final_state.orbit.getPVCoordinates().getPosition().z
relative_error = abs(pos_z_theory - pos_e_experimental) / pos_z_theory * 100

# sample output
# relative error on x : 0.12 %
#  relative error on z : 1276%

# when doing the exact same calculations but without any ForceModels : 
# relative error on x : 0.081 %
#  relative error on z : 496%

# Again the same calculation with all ForceModels except HolmesFeatherStoneAttractionModel
# relative error on x : 0.088 %
#  relative error on z : 401%

So as you can see, those results are quite unsettling… One explanation could be that the satellite performed a maneuver between the two TLE I used, but I tried on many different satellites and despite the amplitude of the values being lower, the weird trend was still present.

Please let me know if i did something wrong, i’m still a beginner. Thank you all for your feedback.

PS : this is the output PVCoordinates with all forces included, the acceleration on the Z axis seems too high

Theo :  {2020-12-09T09:11:13.098, P(-2.724526618536018E7, 1.1576141046260186E7, 9908.796140895838), V(-822.2336477064472, -1936.9220117427658, 3006.0504020394133), A(0.41864020909753735, -0.17787450029587157, -1.5225472418245217E-4)}
Expe :  {2020-12-09T09:11:13.098, P(-2.727990503738411E7, 1.1493375183418749E7, 136348.45491551282), V(-804.3678414387626, -1944.106155548559, 3006.263174630185), A(0.4192015938477335, -0.17661915775838216, -0.002094830244918838)}

I would say that you are mainly observing the errors in the TLE themselves.

TLE are not an accurate model, by far. They are used because for many satellites this is the only available data, but it is very often way off reality.

I remember a presentation by T.S. Kelso where he showed an extreme case with about 60000km error (he confirmed it was really kilometers, so the satellite was on the wrong side of the Earth, and far from it, and there was a logical explanation for that).

1 Like

Thanks for your quick answer !

You’re probably right, i know TLEs are not reliable and the way i convert them into keplerian elements probably does not help, but i did not think it was that bad …
So can you confirm that the method is not wrong ? and if I find more accurate data it should work ?
Also what bothers me is the fact that we see this problem almost only on the Z axis and i can’t explain it.

I’ll go look for more data, thanks for the help !

As propagation errors mainly increase along track, I would suggest to compare position not in the inertial frame in which it is computed, but rather using one satellite as the reference and using its local orbital frame for the axes. So say you have the PVCoordinates pvRef and pvComputed for your two
evaluations at the same time. Then you could compute:

  Rotation toVNC         = LOFType.VNC.rotationToInertial(pvRef);
  Vector3D deltaInertial = pvComputed.getPosition().subtract(pvRef.getPosition());
  Vector3D deltaVNC      = toVNC.applyTo(deltaInertial);

Then you could plot the errors as the components of the deltaVNC vector as the X component would be along velocity, the Y component would be along orbital momentum and the Z component would be roughly radial.

It also seems to me your simple atmosphere model uses a reference altitude that is very low (42000m), so take care of its validity at orbital altitudes.

1 Like

Is it the case that TLEs, are an ‘accurate’ model only for SGP4 propagation and, even then, are not valid for longer than a week, or two?

TLE and SGP4/SDP4 are bound together. TLE can be propagated only by SGP4/SDP4, and SDGP4/SDP4 can propagate only TLE. One is the data structure, the other one is the algorithm that handles it. It is never accurate.

Sometimes, even one week or two is too much, particularly in LEO (i.e. SGP4) when drag is involved.

1 Like

In my opinion, beyond 24H or 48H, the accuracy of the TLE propagation is very limited.

1 Like