Angular parameters of orbit (eg anomaly, longitude argument) get normalized after maneuvers

Hello Orekit Community,

Several orbital parameters in angular units propagate in a non-cyclic and continuously increasing fashion, by design in Orekit. For example anomaly in most orbit and propagation types does not cycle to zero after crossing 360 degrees or 2Pi rads.

I’ve noticed however through tracking such values that doing a maneuver normalizes these quantities back to their 2*Pi modulus, eg. if the spacecrafts pre-maneuver anomaly was 2Pi + 1, in the post maneuver state it will be just 1.

Is this a bug or intended behavrior? It’s requiring special treatment to deal with and to avoid huge steps in the data.

Below is a sample maneuver program in python to reproduce the effect, and a sample of the output printing Epoch and True Longitude Argument from an equinoctial orbit.


Output sample:

Epoch                   True Longitude Argument
2021-01-21T13:31:04.000 110.00577364096642
2021-01-21T13:31:23.695 maneuver trigged by mean anomaly
2021-01-21T13:31:29.000 111.4276796208567  *** OK
2021-01-21T13:31:54.000 112.84885714550067


2021-01-21T15:16:04.000 469.6007670760501
2021-01-21T15:16:29.000 471.0228353012965
2021-01-21T15:16:30.815 maneuver trigged by mean anomaly
2021-01-21T15:16:54.000 112.44424542187713 *** 'Normalized' 112 ~= 471 % 360
2021-01-21T15:17:19.000 113.86493004033836


# Bare-bones maneuver trigger on mean anomaly with out put of epoch and argument of longitude
# to demonstrate the normalization of this angle after a maneuver.

import orekit
from orekit.pyhelpers import setup_orekit_curdir
from org.orekit.attitudes import LofOffset
from org.orekit.utils import Constants, IERSConventions
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.maneuvers import ImpulseManeuver
from org.orekit.frames import FramesFactory, LOFType
from org.orekit.propagation import PropagationType, SpacecraftState
from org.orekit.propagation.conversion import DSSTPropagatorBuilder, DormandPrince853IntegratorBuilder, NumericalPropagatorBuilder, EcksteinHechlerPropagatorBuilder
from import PositionAngleDetector, EventDetector, EventEnablingPredicateFilter
from import EventHandler, PythonEventHandler
from org.orekit.time import TimeScalesFactory, AbsoluteDate, UT1Scale, TimeScales
from org.orekit.orbits import KeplerianOrbit, PositionAngle, OrbitType, EquinoctialOrbit
from org.orekit.propagation.sampling import OrekitFixedStepHandler, PythonOrekitFixedStepHandler, PythonOrekitStepHandler
from org.hipparchus.geometry.euclidean.threed import Vector3D
from import Action
import math

class FsHandler(PythonOrekitFixedStepHandler):
    # Handler to print for Fixed Step propagation mode
    def init(self, s0, t, step):
        # needs to be stated to fulfill the interface specification

    def handleStep(self, currentState, isLast):

class Handler(PythonOrekitStepHandler):
    # Handler to print for variable step propagation mode
    def init(self, s0, t):      # this is different from fixed step handler!
        # needs to be stated to fulfill the interface specification

    def handleStep(self, interpolator, isLast):
        currentState = interpolator.getCurrentState()

class ManHandle(PythonEventHandler):
    # Optional: Do something additional when maneuver trigger condition is found
    # For this example it's just to print in-line when it's triggered
    def init(self, initialstate, target):

    def eventOccurred(self, s, T, increasing):
        kep_orb = KeplerianOrbit(s.getOrbit())
        mean_anom_deg = round(math.degrees(kep_orb.getMeanAnomaly()), 3)
        #print(f'** mean_anom trigger ({trigger_anom_deg}) occurred  at {mean_anom_deg} - lets get a maneuver ***')
        print(f'**\n{s.getDate()} maneuver trigged by mean anomaly\n**')
        return Action.STOP  # this does NOT stop the propagation. Needs to be this value to trigger

def handler_printer(sc_state):
    # This will be called by the different types of handlers, common reporting tasks
    # - Only prints near the trigger condition ie. mean anom (condenses output)
    kep_orb = KeplerianOrbit(sc_state.getOrbit())
    mean_anom_deg = round(math.degrees(kep_orb.getMeanAnomaly()), 1)

    if abs(mean_anom_deg - trigger_anom_deg) <= 5.0:    # Only print data around maneuver trigger point
        sma = round(kep_orb.getA() / 1000, 3)  # km's
        inc_new = round(math.degrees(kep_orb.getI()), 5)
        #print(str(mean_anom_deg) + '\t' + str(sma) + '\t' + str(inc_new))

vm = orekit.initVM()

eme2000 = FramesFactory.getEME2000()
utc = TimeScalesFactory.getUTC()
mu = Constants.EIGEN5C_EARTH_MU
sma_start = 1000.0*1000 + ae        # meters
inc = math.radians(50.0)

startDate = AbsoluteDate(2021, 1, 21, 12, 58, 59.00, utc)

initialOrbit2 = KeplerianOrbit(sma_start,
                               0.01,  # e
                               inc,  # inc
                               math.radians(0.0),   # pa
                               math.radians(10.0),  # raan
                               math.radians(-11.0), # anomaly

Set sat attitude frame:
- VNC:  Velocity - Normal - Co-normal frame (X axis aligned with velocity, Y axis aligned with orbital momentum)
attitudeOverride = LofOffset(eme2000, LOFType.VNC)
sc_state2 = SpacecraftState(initialOrbit2)

# Propagator:
step_sec = 500.0        # propagator, seconds
span_range = 2.0 * initialOrbit2.getKeplerianPeriod()  # propagation, seconds
# - integrator:
min_step = 0.1
max_step = step_sec / 2.0        # integrator, seconds
pos_err = 10.0

# Sat 2:
#prop = 'DSST'
prop = 'NUMERICAL'
#prop = 'ECK'
if prop == 'DSST':
    prop_builder = DSSTPropagatorBuilder(initialOrbit2,
                                         DormandPrince853IntegratorBuilder(min_step, max_step, pos_err),
                                         1.0,  # position 'scaling factor'. Typically standard dev of pos err
                                         PropagationType.MEAN,  # type OF THE ORBIT used for prop
                                         PropagationType.MEAN  # type of elements used to define orbital state
elif prop == 'NUMERICAL':
    prop_builder = NumericalPropagatorBuilder(initialOrbit2,
                                              DormandPrince853IntegratorBuilder(min_step, max_step, pos_err),
                                              PositionAngle.MEAN,  # type OF THE ORBIT used for prop
                                              1.0  # position 'scaling factor'. Typically standard dev of pos err
elif prop == 'ECK':
    prop_builder = EcksteinHechlerPropagatorBuilder(initialOrbit2,
                                                    GravityFieldFactory.getUnnormalizedProvider(6, 6),  # min so far is 6,6

propagator = prop_builder.buildPropagator(prop_builder.getSelectedNormalizedParameters())

# Create maneuver trigger:
trigger_anom_deg = 100.0
print(f'- Trigger MA angle, deg: {trigger_anom_deg} \n')
angle = math.radians(trigger_anom_deg)       # mean anomaly value
trigger = PositionAngleDetector(OrbitType.KEPLERIAN, PositionAngle.MEAN, angle)

# Optionally add custom handler to the trigger (event detector1):
trigger = trigger.withHandler(ManHandle())

# Create maneuver with trigger
dv_vec = Vector3D.PLUS_J.scalarMultiply(0.01)
maneuver = ImpulseManeuver(trigger,
                           dv_vec,  # delta-v dir, along x-axis (i,j,k)
                           100.0  # ISP - required but NOT used for impuslive.

# add maneuver to propagator

use_fixed_handler = False
if use_fixed_handler:
    propagator.setMasterMode(step_sec, FsHandler())

# Set up and run propagation:
time_steps = 500
step_size = round(span_range / time_steps)

for t in range(1, time_steps):
    finalState = propagator.propagate(startDate.shiftedBy(float(t*step_size + step_size)))
    orbit = finalState.getOrbit()
    Lv = EquinoctialOrbit(orbit).getLv()
    epoch = orbit.getDate()

    # Print elements of interest at fixed time steps:
    print(epoch, math.degrees(Lv))

Hi @ShikaDing

There is no bug and everything is expected. You can find a wonderful explanation about the increasing of True/Mean/Eccentric anomalies during orbit propagation in the following post:

For the maneuver, I think this is due to different uses atan2 function which normalize the angle between [-pi, +pi]. To be consistent when printing the anomaly values, maybe you can normalize them before printing.

Best regards,

Hi @bcazabonne ,

I actually had that thread in mind when writing the post - and I agree fully with the design intent of angles evolving over time without limitation/cycling.

However I do find that the resetting of the angles after a maneuver to be counter to this design concept. Yes I could normalize them myself, but then we get away from the benefit of having them evolve continuously as they do without maneuvers. It just seems like a strange exception to have to work around.

Understood that it’s just a feature of the way things work together, so I’ll keep it in mind and account for it.


Thinking about it, I agree with you.

One solution we could propose is to add somewhere (perhaps in OrbitType) a normalizing method in the spirit of MathUtils.normalizeAngle that would take two orbits and bring the angles of the second one back within ±π of the angles of the first one that would be considered a reference. This utility method could then be called at the end of the maneuver computation to avoid too large discontinuties. Of course impulsive maneuvers do introduce discontinuties anyways.

Could you open an issue asking for this new feature?

Done, link to issue.
Thank you.