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

Makss · March 9, 2021, 4:19pm

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)}

luc · March 9, 2021, 5:00pm

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).

Makss · March 9, 2021, 5:18pm

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 !

luc · March 9, 2021, 5:33pm

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.

gavineadie · March 10, 2021, 3:06pm

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?

luc · March 10, 2021, 3:18pm

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.

bcazabonne · March 10, 2021, 3:51pm

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

emasca · July 21, 2022, 3:23pm

Hello, I have a similar issue but unrelated to TLEs propagation.

I am using a NumericalPropagator with DormandPrince853 on Orekit 10.3 to propagate a 530x540km SSO.

I observed that the HolmesFeatherstoneAttractionModel is causing my propagation to mismatch from orbit navigation measures I have. In facts, I noticed that the PerigeeArgumentDot and the RightAscensionOfAscendingNodeDot taken from the KeplerianOrbit at every loop of propagation are very different with respect to the theoretical values calculated with the simplified formulas using J2 coefficient.

In particular I obtain values of this order of magnitude within 2 hours of propagation:

RAAN Drift
from propagator: 1.995089212078394 deg/day
from theory: 0.9935407050657382 deg/day
AOP Drift
from propagator: 8163.319249530916 deg/day
from theory: -3.469592309549168 deg/day

My implementation of the model is the following:

provider = GravityFieldFactory.getNormalizedProvider(120, 120)
zonal = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), provider)
propagator.addForceModel(zonal)

The Integrator tolerances, min and max steps seem not to have effect on this behavior, while the harmoncs order and degree only change the digits of the numbers I get.

Has anybody ever experienced this behavior?

Thanks in advance

luc · July 21, 2022, 3:34pm

Such huge errors are probably related to errors in the initial orbit.
Did you check the units? The Orekit API is in SI units only, i.e. meters and radians.

emasca · July 21, 2022, 3:53pm

Thanks for the super quick reply.

I create the SpacecraftState starting from ITRF PVCoordinates, converting to EME2000 and finally creating a KeplerianOrbit.

The units of ITRF coordinates seem correct to me

ITRF: {2022-11-07T15:54:24.654, P(711429.298, 6872315.899, 52751.354), V(1480.555, -211.137, 7532.18), A(0.0, 0.0, 0.0)}

EME2000: {2022-11-07T15:54:24.654, P(6819613.04813664, 1108644.8353221426, 37744.79895023305), V(118.68071673208637, -986.7960798768436, 7531.96617534098), A(0.18020885558178157, 0.02079268206280105, -3.961225944778964E-4)}

Keplerian parameters: {a: 6912534.274232921; e: 4.7616602313774544E-4; i: 97.5106244557543; pa: 1.1350980740766814; raan: 9.274900184088436; v: -0.8193843093970072;}

luc · July 21, 2022, 4:12pm

This orbit also seems correct to me (when using toString() to print an orbit for display to human beings, the units are converted to degrees, so inclination is correct here).

Do you use other force models? Do you include gravity field only once (I noticed you called it zonal when it fact it already contains both zonal and tesseral terms, so perhaps you have an additional tesseral force model elsewhere?). Could you also check which gravity field is loaded (looking at the output of DataContext.getDefault().getDataProvidersManager().getLoadedDataNames() at the end of the run, when you are sure all needed data as been loaded?

emasca · July 21, 2022, 4:34pm

Yes. we use other force models:

DragForce with NRLMSISE00 and IsotropicDrag
ThirdBodyAttraction for Sun and Moon
SolarRadiationPressure with IsotropicRadiationSingleCoefficient

I attach here what I get from the getLoadedDataNames()
data_providers.txt (27.7 KB)

luc · July 21, 2022, 6:00pm

I tried your input data, using the numerical propagation tutorial as the basis, and
adding the same force models as you.
I get a proper drift for RAAN, as shown with the following plot.
Beware that numerical propagator will get you osculating values, not mean values, and
that simplified formulas with J₂ are really simplified.

emasca · July 22, 2022, 9:19am

Ok, you are right. In fact I now calculated the Perigee drift as:

PerigeerArgumentDot + TrueAnomalyDot - KeplerianMeanMotion

and got a much more familiar result:

aop_dot

However, I am still concerned about the behavior of the Perigee and Apogee. In the following figure I plotted the altitude of my orbit (left Y axis, red) and the ApsideDetector g function (right Y axis, black):
apside_altitude
Altitude is calculated by taking the norm of the PVCoordinates Position in EME2000 and subtracting the WGS84 Earth equatorial radius.

As you may notice, the apogees and the perigees are about 45-50 minutes apart, while the KeplerianPeriod of this orbit is around 95 minutes.

I have orbit navigation measurements of real satellites in similar orbits and this behavior is not present.

Edit: let me specify that altitude trend in this short-term propagation is as expected without the Holmes-Featherstone model i.e. 95 minutes between two consecutive apogees or perigees.

luc · July 22, 2022, 10:05am

This plot seems correct to me. The g function of the apside detector does not represent an altitude but a dot product between position and velocity. It therefore crosses 0 at apsides when position and velocity are orthogonal to each other. The crossing is from negative values to positive values (i.e. an increasing event) at perigee, and the crossing is from positive values to negative values (i.e. a decreasing event) at apogee. This is indeed what we see in the plot: for example when the red curve reaches its min and max values, the black curve crosses its zero line.

emasca · July 22, 2022, 10:21am

What I don’t think it is correct, is the actual trend of the altitude. If you see it shows an apogee around 16:40 UTC, and another one around 17:30 UTC, so roughly 50 minutes after.

This is not coherent with the 95 minutes period the orbit should have.

luc · July 22, 2022, 11:22am

Oh, I see.

Don’t you think it is just the classical figure-8 behavior of eccentricity on almost frozen orbits?
Look at this post. The
osculating parameters in the special case of SSO and frozen eccentricity are really different from mean parameters. The eccentricity short periodic terms are larger than the mean eccentricity. This implies that the perigee revolves around orbit once per period. In this post, I explain that for such orbits “it is not the satellite that crosses its perigee once per orbit, it is the perigee that crosses the satellite once per orbit”. So getting two perigees per orbit is mind-blowing but true. It in fact happens as soon as you use J₂ (the global figure-8 shape), or J₂ + J₃ (the global figure-8 shape plus its shift along ey). What is counter-intuitive is that short periods are so large in this case, but it is because by design these orbits have collapsed all other terms down to zero, hence the short periodic terms show up.

Another consequence of this is that true anomaly does not increase steadily anymore. In fact, is oscillates and goes backward at some points! This is the result of two different effects: the large short period terms of eccentricity and the singularity of true anomaly for near singular orbits.

emasca · July 22, 2022, 1:46pm

That was super-interesting. However I was still not convinced because I ddon’t see this behaviour in real-life measurements of other SSO orbits.

I tested the same propagator on one of these real-life dataset and the behavior is not there anymore. The propagation follows the measurements quite steadly and altitude behaves intuitively.

Here a plot:
apside_altitude

So I see two possible explanations:

the SSO orbit I was testing before is a peculiar orbit (not very convinced by this)
there is some hidden issue in the model implementation (currently outside my field of expertise)

I can try to test the same propagator on other orbit cases. Do you have some suggestions?

luc · July 22, 2022, 2:13pm

Indeed, your orbit is special. It seems to be very far from a frozen eccentricity, as its ey component is close to 0 instead of the traditional 1.2e-3 (which is related to J₃).

When plotting the circular eccentricity vector, the figure-8 shape is almost symmetrical along ey.

Here the mean eccentricity (the center of the oval part of the curve) is far from frozen eccentricity. This implies that on the long term (order of 4 months), the full figure will drift a lot as it will rotate around the frozen eccentricity point.

Could you try shifting your orbit so it gets closer to frozen eccentricity? You can use the Phasing tutorial for that.

luc · July 22, 2022, 2:37pm

Here is the eccentricity plot on a 120 days period, showing how the curve moves.
If eccentricity was frozen, we would have only a slight fattening of the curve, but here
we have a big cloud of points. The frozen eccentricity at positive ey appears clearly,
and obviously your orbit is not there.

The color code of the curve is time from propagation start, in seconds.