Apside Detection: Numerical Propagator (Osculating state) vs DSSTPropagator

Dear Orekit Community,

I’m using the Python Wrapper for orekit.
My orekit version is 11.3

My goal is to detect an apside (Perigee or Apogee).
So I have this code

    state = prop.getInitialState()
    state_orbit = state.getOrbit()
    apside_detector = ApsideDetector(0.001, state_orbit)
    if apside == 'Apogee':
        apside_detector = apside_detector.withHandler(StopOnDecreasing())
    prop.addEventDetector(apside_detector)
    orbit_period = state_orbit.getKeplerianPeriod()
    apside_state = prop.propagate(state.getDate().shiftedBy(1.1 * orbit_period))

the variable apside is an input (string) that can be ‘Apogee’ or ‘Perigee’, depending on the apside that I want to detect.
The variable prop is a propagator. This propagator can be a NumericalPropagator or a DSSTPropagator

The way I create a NumericalPropagator is here below

    integrator = DormandPrince853Integrator(0.0000001, 2700.0, 9.999999999999999e-12, 9.999999999999999e-12 )
    propagator = NumericalPropagator(integrator)

    cd = 2.2 
    cr = 1.7
    drag_area = 1.358
    srp_area = 1.358
    harmonics_degree = 64 
    harmonics_order = 64
    strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE

    # DRAG
    supported_names = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
    input_data = MarshallSolarActivityFutureEstimation(supported_names, strength_level)
    atmosphere = NRLMSISE00(input_data, constants.sun, constants.earth)
    isotropic_drag = IsotropicDrag(drag_area, cd)
    drag = DragForce(atmosphere, isotropic_drag)
    propagator.addForceModel(drag)

    # SRP
    radiation_sensitive = IsotropicRadiationSingleCoefficient(srp_area, cr)
    srp = SolarRadiationPressure(constants.sun, constants.R_E, radiation_sensitive)
    propagator.addForceModel(srp)

    # Third Body
    moon3b = ThirdBodyAttraction(constants.moon)
    sun3b = ThirdBodyAttraction(constants.sun)
    propagator.addForceModel(moon3b)
    propagator.addForceModel(sun3b)

    # Zonal Harmonics
    provider = GravityFieldFactory.getNormalizedProvider(harmonics_degree, harmonics_order)
    zonal = HolmesFeatherstoneAttractionModel(constants.fixedFrame, provider)
    propagator.addForceModel(zonal)

The way I create a DSSTPropagator is here below

    integrator = DormandPrince853Integrator(0.0000001, 2700.0, 9.999999999999999e-12, 9.999999999999999e-12 )
    dsst_propagator = DSSTPropagator(integrator, PropagationType.MEAN)

    cd = 2.2 
    cr = 1.7
    drag_area = 1.358
    srp_area = 1.358
    harmonics_degree = 64 
    harmonics_order = 64
    strength_level = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE

    # DRAG
    supported_names = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
    input_data = MarshallSolarActivityFutureEstimation(supported_names, strength_level)
    dsst_atmosphere = NRLMSISE00(input_data, constants.sun, constants.earth)
    dsst_drag = DSSTAtmosphericDrag(dsst_atmosphere, cd, drag_area, constants.mu)
    dsst_propagator.addForceModel(dsst_drag)

    # SRP
    dsst_srp = DSSTSolarRadiationPressure(cr, srp_area, constants.sun, constants.R_E, constants.mu)
    dsst_propagator.addForceModel(dsst_srp)

    # Third Body
    dsst_moon = DSSTThirdBody(constants.moon, constants.mu)
    dsst_sun = DSSTThirdBody(constants.sun, constants.mu)
    dsst_propagator.addForceModel(dsst_moon)
    dsst_propagator.addForceModel(dsst_sun)

    # Zonal and Tesseral Harmonics
    provider = GravityFieldFactory.getUnnormalizedProvider(harmonics_degree, harmonics_order)
    dsst_zonal = DSSTZonal(provider)
    dsst_propagator.addForceModel(dsst_zonal)

    dsst_tesseral = DSSTTesseral(constants.fixedFrame, constants.w_E, provider)
    dsst_propagator.addForceModel(dsst_tesseral)

The initial state of a NumericalPropagator is always an osculating state.
The initial state of a DSSTPropagator is always a mean state

Example:
When I detect the perigee using the NumericalPropagator and an osculating state,
Input:

State Type: OSCULATING
State: EME2000 {2023-12-09T08:00:00.000, P(4030126.8284191545, 4921843.180766312, -2799991.6722682836), V(2774.1232318262637, 1630.0034760832289, 6851.526650257949), A(-4.785816425329783, -5.844849373511939, 3.334164249424359)}
Sat Mass [kg]: 236.32
SMA [km]: 6943.723235401812
Eccentricity: 0.0010076720807084784
Inclination [°]: 97.73925894828648

Output (Perigee detected):

State Type: OSCULATING
State: EME2000 {2023-12-09T08:41:32.66744242998948, P(-2623975.7447196837, -3869357.272155609, 5122457.127232371), V(-4334.009454997956, -3690.7982749619277, -5008.017089338541), A(3.128110025173131, 4.6130243048460615, -6.1237018358396735)}
Sat Mass [kg]: 236.32
SMA [km]: 6936.477215477286
Eccentricity: 0.00018680782661333016
Inclination [°]: 97.74396675829782


When I detect the perigee using the DSSTPropagator and a mean state,
Input:

State Type: MEAN
State: EME2000 {2023-12-09T08:00:00.000, P(4030104.56652526, 4929036.597456111, -2773780.2724690065), V(2753.755555065786, 1612.8769826599641, 6866.525882432539), A(-4.790062485596439, -5.857554535534209, 3.2966029808461985)}
Sat Mass [kg]: 236.32
SMA [km]: 6937.43806449925
Eccentricity: 0.0010694630404250913
Inclination [°]: 97.73027302870626

Output (Perigee detected):

State Type: MEAN
State: EME2000 {2023-12-09T08:48:24.20415838127064, P(-4086382.6875770474, -4959621.763395213, 2593947.8119772924), V(-2636.5061570213365, -1468.9675435211052, -6962.089380103862), A(4.888210975166692, 5.931909264413749, -3.102750345982305)}
Sat Mass [kg]: 236.32
SMA [km]: 6937.437203379718
Eccentricity: 0.001071656748111149
Inclination [°]: 97.73028535976285

The question is: what can explain the difference in time on the Perigee detected?
On NumericalPropagator(Osculating State) the perigee occurs at 08:41:32.667
On DSSTPropagator(Mean State) the perigee occurs at 08:48:24.204
7 min difference is a considerable difference, higher than I was expecting.
Which method (NumericalPropagator or DSSTPropagator) should be used to detect the apside?

Thank you in advance,
Looking forward for your feedback

Hi @jpedro500

That’s a known issue. Because DSST integrates the mean elements, the event detection will be based on the mean state. So, I’m not surprised to see differences, even big depending the dynamical model used to compute the mean parameters.

Because the numerical propagator computes events based on the oscillating state, I can recommend you to use it.

Best regards,
Bryan