Hello,
Noticing unexpected behavior in python with a simple scenario where triggering a maneuver based on event detector (mean anomaly from PositionAngleDetector) does not trigger on the first occurrence of the condition (on the first orbit), only on the n+1 occurrences. Is this by design?
Interestingly, if I add a custom handler to the event detector using .withHandler(), to simply print a message upon event detection and explicitly return an Action.STOP, it works on every orbit.
The basic flow in this example uses a fixed or variable step handler (configurable for testing) on the DSST propagator, in master mode, to print orbit values as we approach the trigger value, which also crudely confirms whether or not the maneuver actually happens (changes in orbit elements visible, etc).
In propagating for exactly 2 orbits, we get one maneuver on the 2nd orbit without this simple additional handler:
results_without_handler.txt (420 Bytes)
And two maneuvers with the additional handler (via trigger.withHandler(ManHandle), which I expected to be the default behavior:
results_with_handler.txt (609 Bytes)
And below is the code:
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.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
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
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 ***')
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
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.1, # 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 = 50.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_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
)
prop_builder.setMass(200.0)
propagator = prop_builder.buildPropagator(prop_builder.getSelectedNormalizedParameters())
# Create maneuver trigger:
trigger_anom_deg = 1.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 detector):
trigger = trigger.withHandler(ManHandle())
# Create maneuver with trigger
dv_dir = Vector3D.PLUS_J.scalarMultiply(15.0)
maneuver = ImpulseManeuver(trigger,
attitudeOverride,
dv_dir, # 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())
print('MeanAnom\tSMA\tInc.')
finalState = propagator.propagate(startDate.shiftedBy(span_range))
finalOrbit = finalState.getOrbit()
finalOrb_Kep = KeplerianOrbit(finalOrbit)
final_inc = round(finalOrb_Kep.getI(), 5)
final_SMA = round(finalOrb_Kep.getA() / 1000, 5)
final_RAAN = round(math.degrees(finalOrb_Kep.getRightAscensionOfAscendingNode()), 5)
inc_delta = round(math.degrees(final_inc - inc), 3)
sma_delta = round(final_SMA - sma_start/1000, 3)
print(f"\nSMA delta, km: {sma_delta}")
print(f"Inc. delta, deg: {inc_delta}")
print(f"RAAN: {final_RAAN}")