Brouwer-Lyddane propagation precision issue

Orekit usage
1

Hi, when using the Brouwer-Lyddane propagator in Orekit to extrapolate the orbit of a low earth circular orbit satellite, the error reached 50 kilometers during a one-hour extrapolation. What could be the possible reasons? the reference value is close to the numerical propagated results.

final Frame inertialFrame = FramesFactory.getEME2000();
        AbsoluteDate initDate = new AbsoluteDate(2024, 12, 1, 0, 0, 0, TimeScalesFactory.getUTC());
        double timeshift = 3600 ;
        NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 0);

        // Initial orbit

        final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(new Vector3D(-5053374.356, -3419077.892, -4241149.605),
                new Vector3D(-4221.103996, -1055.944563, 5887.522421)),
                inertialFrame, initDate, provider.getMu());
        // Initial state definition
        final SpacecraftState initialState = new SpacecraftState(initialOrbit);

        //_______________________________________________________________________________________________
        // SET UP A REFERENCE NUMERICAL PROPAGATION
        //_______________________________________________________________________________________________

        // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
        final double minStep = 0.001;
        final double maxstep = 1000.0;
        final double positionTolerance = 10.0;
        final OrbitType propagationType = OrbitType.CIRCULAR;
        final double[][] tolerances =
                ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit, propagationType);
        final AdaptiveStepsizeIntegrator integrator =
                new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);

        // Numerical Propagator
        final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
        NumPropagator.setOrbitType(propagationType);

        final ForceModel holmesFeatherstone =
                new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
        NumPropagator.addForceModel(holmesFeatherstone);

        // Set up initial state in the propagator
        NumPropagator.setInitialState(initialState);

        // Extrapolate from the initial to the final date
        final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
        final TimeStampedPVCoordinates numPV = NumFinalState.getOrbit().getPVCoordinates();

        //_______________________________________________________________________________________________
        // SET UP A BROUWER LYDDANE PROPAGATION
        //_______________________________________________________________________________________________

        BrouwerLyddanePropagator BLextrapolator =
                new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);

        SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
        TimeStampedPVCoordinates blPv = BLFinalState.getPVCoordinates();
2

Hello @afei and welcome to the Orekit forum !

After further investigations, i believe your issue comes from the Orekit data you are using.

When i used your example, i had a difference of 17km. I then tried several settings and other propagators but once i updated my orekit data, it went down to 2.6 km.

Could you launch the update.sh script located in your orekit data and tell us how it goes ?

Cheers,
Vincent

3

Thank you very much for your reply. @Vincent I’m curious to know which files in the orekit-data might affect the accuracy of the BL propagator. Since the orekit-date was updated in 2025, the propagated data is 2024

4

I updated the orekit data, but the result accuracy hasn’t shown significant improvement…

5

Ok, i wanted to make sure since i did not have the same error to begin with.

I see that you based your code from an Orekit test which is a good start. I’ll keep you updated :+1: .

6

I apologize for not expressing the issue clearly. Using the parameter settings in the code for extrapolation, I obtained a deviation of 17,340 meters. During testing, I found that the maximum hourly extrapolation deviation can reach up to 50,000 meters. Even after updating Orekit-data, the accuracy still cannot be improved.

First, I would like to know the approximate accuracy range of the BLpropagator for low-Earth orbit (LEO) satellites. Second, since my M2 parameter is set to 0, which files in Orekit-data need to be updated when using the BLpropagator?

7

Brouwer-Lyddane is a very crude model, it only takes very few perturbations in account.

We used it a lot in the 90’s when computers where not as fast as they are now, but not so much anymore.

Is your satellite very low? Gravity field effects are more important when satellite is closer to Earth and as Brouwer-Lyddane only considers very few zonal parameters, its accuracy is very limited.

8

The orbital altitude of my satellite is 1050km, with an eccentricity <= 0.03 and an orbital inclination of 99 degrees

9

Hello @luc,

What i find strange is that the code is comparing a brouwer-lyddane propagation with a numerical one with only a HolmesFeatherstoneAttractionModel of degree=5 and order=0. Hence there are no other perturbations and i would expect to find almost identical results in such case. I must be missing something :confused: .

Cheers,
Vincent

10

mean versus osculating parameters, perhaps ?

11

I made a plot to better understand the error evolution:

image
image1657Ă—906 85.6 KB

I still don’t get why we have discrepencies with the numerical propagation given that it only has potential of degree = 5, order = 0. At least we do not have a secular term as with the Keplerian propagation.

Cheers,
Vincent

12

Vincent have you tried passing GCRF as the rotating frame for the geopotential (which is physically wrong, let’s be clear)?
It looks like the BL does not need a non-inertial frame anyway.

13

Did not think about this, unfortunately i don’t see much difference (~14m).

You will find my code below (quick script so please don’t mind the quality :grimacing:):

It expects orekit-data folder at home directory and requires an environment with orekit + plotly.

import os
from pathlib import Path

import orekit
import plotly.graph_objs as go
from orekit import JArray_double
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory, UnnormalizedSphericalHarmonicsProvider
from org.orekit.frames import FramesFactory
from org.orekit.orbits import KeplerianOrbit, OrbitType
from org.orekit.propagation import SpacecraftState, ToleranceProvider, Propagator, PropagationType
from org.orekit.propagation.analytical import BrouwerLyddanePropagator, KeplerianPropagator
from org.orekit.propagation.conversion.osc2mean import FixedPointConverter, BrouwerLyddaneTheory
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import PVCoordinates, IERSConventions

vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir

setup_orekit_curdir(os.path.join(Path.home(), "orekit-data"))

class RelativeBLPositionHandler(PythonOrekitFixedStepHandler):

    def __init__(self, t0: AbsoluteDate, other: Propagator, provider: UnnormalizedSphericalHarmonicsProvider):
        super().__init__()
        self.other = other
        self.t0 = t0
        self.positions = []
        self.dates = []

        theory = BrouwerLyddaneTheory(provider, 0.)
        self.converter = FixedPointConverter(1e-12, 500, FixedPointConverter.DEFAULT_DAMPING)
        self.converter.setMeanTheory(theory)

    def finalize(self):
        pass

    def finish(self, spacecraftState):
        pass

    def handleStep(self, spacecraftState):
        other_orbit = self.other.propagate(spacecraftState.getDate()).getOrbit()
        other_pv = self.converter.convertToMean(other_orbit).getPVCoordinates()

        current_pv = spacecraftState.getPVCoordinates()
        self.positions.append(PVCoordinates(current_pv, other_pv).getPosition().getNorm())
        self.dates.append(spacecraftState.getDate().durationFrom(self.t0))

    def init(self, spacecraftState, absoluteDate, double):
        pass


class RelativePositionHandler(PythonOrekitFixedStepHandler):

    def __init__(self, t0: AbsoluteDate, other: Propagator):
        super().__init__()
        self.other = other
        self.t0 = t0
        self.positions = []
        self.dates = []

    def finalize(self):
        pass

    def finish(self, spacecraftState):
        pass

    def handleStep(self, spacecraftState):
        other_pv = self.other.propagate(spacecraftState.getDate()).getPVCoordinates()
        current_pv = spacecraftState.getPVCoordinates()
        self.positions.append(PVCoordinates(current_pv, other_pv).getPosition().getNorm())
        self.dates.append(spacecraftState.getDate().durationFrom(self.t0))

    def init(self, spacecraftState, absoluteDate, double):
        pass


# Dates and frame
utc = TimeScalesFactory.getUTC()
init_date = AbsoluteDate(2024, 12, 1, 0, 0, 0.0, utc)
time_shift = 3600.0 * 6
final_date = init_date.shiftedBy(time_shift)
inertial_frame = FramesFactory.getEME2000()
rotating_frame = FramesFactory.getITRF(IERSConventions.IERS_2010, True)

# Gravity provider
normalized_provider = GravityFieldFactory.getNormalizedProvider(5, 0)
unnormalized_provider = GravityFieldFactory.getUnnormalizedProvider(normalized_provider)

# Initial Orbit
position = Vector3D(-5053374.356, -3419077.892, -4241149.605)
velocity = Vector3D(-4221.103996, -1055.944563, 5887.522421)
initial_orbit = KeplerianOrbit(PVCoordinates(position, velocity),
                               inertial_frame,
                               init_date,
                               normalized_provider.getMu())

# Initial state
initial_state = SpacecraftState(initial_orbit)

###################### BROUWER-LYDDANE ######################
# Brouwer-Lyddane propagator
# BL_extrapolator = BrouwerLyddanePropagator(
#     initial_orbit,
#     FrameAlignedProvider(inertial_frame),
#     1000.,
#     Constants.IERS2010_EARTH_EQUATORIAL_RADIUS,
#     Constants.IERS2010_EARTH_MU,
#     Constants.EIGEN5C_EARTH_C20,
#     Constants.EIGEN5C_EARTH_C30,
#     Constants.EIGEN5C_EARTH_C40,
#     Constants.EIGEN5C_EARTH_C50,
#     PropagationType.OSCULATING,
#     0.
# )

bl_extrapolator = (
    BrouwerLyddanePropagator(
        initial_orbit,
        unnormalized_provider,
        PropagationType.OSCULATING,
        0.))

kp_extrapolator = KeplerianPropagator(initial_orbit)

###################### NUMERICAL ######################
min_step = 0.001
max_step = 100.0
position_tolerance = 1e-3
propagation_type = OrbitType.EQUINOCTIAL
tolerances = ToleranceProvider.getDefaultToleranceProvider(position_tolerance).getTolerances(initial_orbit,
                                                                                             propagation_type)
integrator = DormandPrince853Integrator(min_step, max_step,
                                        JArray_double.cast_(tolerances[0]),
                                        JArray_double.cast_(tolerances[1]))

num_propagator = NumericalPropagator(integrator)

num_propagator.setInitialState(initial_state)

num_propagator.setOrbitType(propagation_type)

holmes_featherstone = HolmesFeatherstoneAttractionModel(rotating_frame, normalized_provider)
num_propagator.addForceModel(holmes_featherstone)

# Add handlers
bl_position_handler = RelativePositionHandler(init_date, bl_extrapolator)
kp_position_handler = RelativePositionHandler(init_date, kp_extrapolator)
bl_mean_position_handler = RelativeBLPositionHandler(init_date, bl_extrapolator, unnormalized_provider)

timestep = 60.0
num_propagator.getMultiplexer().add(timestep, bl_position_handler)
num_propagator.getMultiplexer().add(timestep, kp_position_handler)
num_propagator.getMultiplexer().add(timestep, bl_mean_position_handler)

# Propagate numerically
num_final_state = num_propagator.propagate(final_date)

# Plot differences
fig = go.Figure()

fig.add_trace(go.Scatter(
    x=bl_position_handler.dates,
    y=bl_position_handler.positions,
    mode='markers',
    name='Brouwer-Lyddane'
))

fig.add_trace(go.Scatter(
    x=bl_mean_position_handler.dates,
    y=bl_mean_position_handler.positions,
    mode='markers',
    name='Mean Brouwer-Lyddane'
))

fig.add_trace(go.Scatter(
    x=kp_position_handler.dates,
    y=kp_position_handler.positions,
    mode='lines',
    name='Keplerian'
))

fig.update_layout(
    title="Comparison w.r.t numeric propagation [5,0]",
    xaxis_title="Time from start [s]",
    yaxis_title="Position error L2 norm [m]"
)
fig.show()