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.
Thanks
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
Code:
# 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 org.orekit.propagation.events import PositionAngleDetector, EventDetector, EventEnablingPredicateFilter
from org.orekit.propagation.events.handlers 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 org.hipparchus.ode.events 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
pass
def handleStep(self, currentState, isLast):
handler_printer(currentState)
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
pass
def handleStep(self, interpolator, isLast):
currentState = interpolator.getCurrentState()
handler_printer(currentState)
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):
pass
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()
setup_orekit_curdir()
eme2000 = FramesFactory.getEME2000()
utc = TimeScalesFactory.getUTC()
mu = Constants.EIGEN5C_EARTH_MU
ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
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
PositionAngle.MEAN,
eme2000,
startDate,
mu
)
'''
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
PositionAngle.MEAN,
1.0
)
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,
attitudeOverride,
dv_vec, # delta-v dir, along x-axis (i,j,k)
100.0 # ISP - required but NOT used for impuslive.
)
# add maneuver to propagator
propagator.addEventDetector(maneuver)
use_fixed_handler = False
if use_fixed_handler:
propagator.setMasterMode(step_sec, FsHandler())
else:
propagator.setMasterMode(Handler())
# 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))