Harmonic Acc Question

Orekit Python Wrapper usage
1

Hello!

I am trying to use a HarmonicAccelerationModel() to account for some forces in LEO. I have a GPS dataset that I trust and I’m trying to fit an orbit to it using BLS. In theory I think this should be fairly simple, fitting to a number of GPS PV measurements. However, I think there are some periodic forces that are not accounted for in the Orekit modeling/propagator.

I took the acceleration of the dataset and compared it with the accelerations in Orekit, subtracted them and plotted to see what forces I’m missing. Here is an example of the accelerations in one body frame axis (m/s^2)

image
image986×347 20.8 KB

Clearly, the force(s) I’m “missing” are periodic. The frequency of that signal is nearly exact to one half orbit.

I’m trying to account for this discrepancy by using a HarmonicAccelerationModel() in each satellite body axis (LVLH). The part I don’t understand is fundamentalPeriod vs harmonicMultiplier. In this example where the force has a period of 1/2 orbit, what should the fundamentalPeriod be? I cannot set the harmonicMultiplier to 0.5, it needs to be an integer (and why is that?). I may be misunderstanding the use case of HarmonicAccelerationModel().

2

Clarification:
When I took the accelerations to compare, I grabbed the PV of each GPS measurement and extracted propagated it for one timestep in Orekit, getting the Orekit acceleration at that point. I did this for each GPS point.
So the Orekit and the GPS orbits did not diverge and I was getting the accelerations at each point along the GPS dataset orbit.

3

Hi @dbyrne,

It seems there’s a bias (0.12 m/s²) and a harmonic signal (0.03 m/s²) in the difference plot. And it’s a lot!
What force models do you use?

The angular frequency used to compute the amplitude of the signal is \omega = \frac{2k\pi}{T}
With k the multiplier and T the fundamental. If you want to have an \omega consistent with half an orbit, then you need to use k=2 and T = keplerian period, if I’m not mistaken.

Hope this helps,
Maxime

1 Like
4

@MaximeJ
That plot is the Z-axis in the LVLH_CCSDS body frame, so it should be mostly gravity related terms. I’ll put my code for creating the Earth and gravity below, in case there are any glaring issues. Otherwise I am including NRLMSISE00 drag, moon and sun ThirdBodyAttraction, SolidTides and OceanTides, and SolarRadiationPressure. The drag and SRP act on an isotropic body.

I agree, it is quite a big force difference. In the X and Y axes the differences are a few orders of magnitude smaller. I’ve also tried adding a 0th order (constant) PolynomialAcceleration to account for that bias in the Z axis.

Maybe the issue is the way I am updating EOPs? Here’s the code for that too:

    eopfile = 'finals2000A.txt'
    urllib.request.urlretrieve(url, eopfile)
    current_dir = os.getcwd()
    directory   = File(current_dir)
    eopdata     = DirectoryCrawler(directory)
    DataContext.getDefault().getDataProvidersManager().addProvider(eopdata)

And my code for setting up the Earth object:

    earthConventions = IERSConventions.IERS_2010
    earthRadius      = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
    earthFlattening  = Constants.WGS84_EARTH_FLATTENING
    earthBodyFrame   = FramesFactory.getITRF(earthConventions, False)
    earthShape = OneAxisEllipsoid(earthRadius, earthFlattening, earthBodyFrame)
    earthShape = ReferenceEllipsoid.getWgs84(earthBodyFrame)
    earthMu    = Constants.WGS84_EARTH_MU
    gravity_provider = GravityFieldFactory.getNormalizedProvider(200, 200)
    gravity_model = HolmesFeatherstoneAttractionModel(earthBodyFrame, gravity_provider)```
5

Hi @dbyrne,

I don’t see any striking mistake in your code. Not sure that bad EOPs would give such a difference in acceleration. You can try without them and see what the difference is.

Is the acceleration from your GNSS “inertial”? Or is it also in an Earth-rotating frame?
I’m thinking of frame transforms and relative/absolute acceleration issues but it’s just a wild guess.

1 Like
6

@MaximeJ,
Thanks for looking that over. I’ve run into another question while looking through my code. When I initialize a NumericalPropagator with a CartesianOrbit it will automatically apply a NewtonionAttraction force. Is the HolmesFeatherstoneAttractionModel supposed to be applied on top of this, or replacing it?

Maybe I am initializing the propagator incorrectly, because then I cannot remove the NewtonionAttraction with removeForceModels()

        integrator = DormandPrince853Integrator(min_step, max_step, pos_error, vel_error)
        integrator.setInitialStepSize(step)
        propagator = NumericalPropagator(integrator)
        propagator.setOrbitType(OrbitType.CARTESIAN)
        state = SpacecraftState(orbit, mass)
        propagator.setMu(earth_params['earthMu'])
        propagator.setInitialState(state)
7

Hi there,

If your propagator is orbit based the Newtonian part is automatically added indeed.
If you want geopotential terms on top of it, you need to add the Holmes model yes.

Cheers,
Romain

1 Like
8

@Serrof
Does this also apply for a NumericalPropagatorBuilder that I give to an estimator? i.e. I need to add both the NewtonionAttraction and the HolmesFeatherstone to the builder?
Thanks for the input

9

Hi @dbyrne,

NewtonianAttraction is added on the fly when building the propagator (see here) if you didn’t add it manually.
Does that work for you?

Maxime

1 Like
10

@MaximeJ , @Serrof
Thanks for the combined help.

For those interested, my problem came from two places:
My misunderstanding of which gravity force models to use, as mentioned above.
A frame conversion issue, the SP3 file I used was already converted into an intertial/GCRF frame. I though SP3 convention was to use an ECEF frame and was so I was converting it again and fudging with the frame unnecessarily.

A new plot of the differenced accelerations in the body Z-axis: (still an interesting harmonic nature)

image
image1920×325 81 KB

1 Like