While developing features for our service using Orekit propagation via the Python wrapper, I identified a discrepancy with a known orbit that I am unable to explain. Specifically, the example is a simplified version of the ISS’s, reduced to a perfectly equatorial, circular orbit in the J2K frame. The initial PV coords are:
Modelling this orbit in STK produces the expected output, a nearly perfectly periodic orbit with a period of ~93 min at a near constant altitude of 420 km. However, running it through Orekit (DormandPrince) produces the following:
As you can see, the orbit period is correct, but the altitude dips by 20 km rather than remaining perfectly circular (STK reported 419.999 km). To dispel some of the potential issues I have tracked down:
This altitude is not measured in a smart way (it’s just the norm of the ECI pos minus the earth’s radius), so this isn’t a result of surface curvature. The orbit should be perfectly circular.
It also isn’t a result of atmospheric drag or SRP effects, as the feature persists with both of those forces disabled.
If I change the altitude of the orbit or make it more eccentric, the low 20 km error persists at perigee.
The state values are consistent all the way back to where the propagator multiplexor returns them, so this is happening within Orekit somewhere.
The absolute tolerance is fixed to a very low value, so this isn’t just internal integrator wiggle.
I’d appreciate any suggestions of what might be causing this effect! At the very least, if someone could send my orbit state through their own install and confirm if they also observe the error, I would greatly appreciate it!
Are the inclination & RAAN staying put or are they moving around too? (That’s what I use to work out if there’s a force hanging around that I didn’t want )
I’ll try to be brief, but I am not sure if J_2 effects are included in the force models we use. I’ll add all of them in case any of the others draw attention:
Gravity: HolmesFeatherstoneAttractionModel() with NormalizedSphericalHarmonicsProvider() and a OneAxisEllipsoid() for the Earth.
Tides: SolidTides() and OceanTides()
Drag: NRLMSISE00()
SRP: SolarRadiationPressure() with IsotropicRadiationSingleCoefficient()
3rd Body: Sun and Moon
There is no explicit use of any of the J_2 providers I see in the API, but perhaps they are added automatically through the harmonics provider (order = degree = 70)? An additional gravity term would be the conservative force effect we are observing, so it is a good suggestion!
I’ll add those terms explicitly in STK and check if that creates the effect.
Thank you Luc, it was indeed J_2 effects that were squishing the orbit. I didn’t think these were implemented, but reducing the order/degree of the gravity harmonics to 0 produced the expected circular orbit with minor drift from other forces:
I guess the giveaway for this was that the mystery effect was acting in a conservative way, squishing the orbit slightly but ultimately returning to the same starting altitude/velocity.
The J₂ effect is in fact part of the Holmes-Featherstone attraction model. This is a representation of the irregularities of the gravity field in spherical harmonics. J₂ is the largest of these spherical harmonics and for low Earth orbits like the ISS, it is the most important force (together with atmospheric drag) after central attraction. When you reduce the degree to 0, you end up with central attraction only (note that the degree1 harmonics always have zero amplitude when the origin of the frame is taken at the center of mass of the Earth, so you end up with only the 0 degree term which is the central attraction). When you do that, you surely should also ignore tides, SRP and 3rd body, because they are smaller and it would be strange to ignore the big J₂ term and consider other forces that are smaller than that, but I would in fact rather recommend keeping these forces and putting J₂ (and many other spherical harmonic terms) back into the simulation, in both Orekit and STK. For realistic simulations, especially in low Earth orbit, these terms are really needed.
ISS case is however a little special because it is so low and so big, its drag effect is difficult to model properly.
Looking further at your initial post, rather than the J₂ term, it is probably the C₂₋₂ term that is the culprit, but reducing degree to 0 also cancelled it. The reason I think it is this term is because you said the orbit was equatorial, and on an equatorial orbit you are always at the same latitude, and the J₂ term effect depends only on latitude, so it is constant on equatorial orbits. The C₂₋₂ term, on the other hand, depends on longitude, so I guess it created the effect you noticed.
You are absolutely correct. Apologies for my imprecise language, I meant to say that higher order perturbations (J_2 and up) were the cause of the issue. At the time of writing, I hadn’t yet verified which specific term was the largest contributor. As you suggested, it was in fact the C_{2-2} term.
)