Thanks @petrus.hyvonen,
I’ve tried implementing Luc’s solution in my Python class:
from org.hipparchus.util import FastMath
from org.hipparchus.ode.events import Action
from org.orekit.frames import LOFType
from org.orekit.orbits import KeplerianOrbit
from org.orekit.attitudes import LofOffset
from org.orekit.propagation.events import DateDetector, NodeDetector, PythonEnablingPredicate, EventEnablingPredicateFilter, PythonAbstractDetector
from org.orekit.propagation.events.handlers import PythonEventHandler, StopOnEvent
from org.orekit.forces.maneuvers import ImpulseManeuver
from org.orekit.utils import Constants
from java.util import ArrayList as ArrayList
# constants
tformat = '%Y-%m-%dT%H:%M:%S'
dtsl = 19 # [int] date time string length
TAtol = 10.0 # [deg] True Anomaly tolerance for Apoapsis/Periapsis detection
# global variables
StopDetecting = False
enabledEvent = False
def add_impulsive_maneuver(self,inertialFrame,firingDate,mandv,mantype,isp):
'''
This method inserts an impulsive maneuver in the propagator.
Rules:
- in plane maneuvers trigger: firingDate
- out of plane maneuvers trigger: first Acending Node after firingDate
inertialFrame : [Frame] Frame for the attitudeOverride
firingDate : [AbsoluteDate] Maneuver minimum time (or exact time when mantype == 'T')
mandv : [3D vector] Maneuver delta v vector
mantype : [str] Maneuver type ('AN' Ascending Node; 'DN' Descending Node; 'A' @ Apogee; 'P' @ Perigee; 'T' Time triggered)
isp : [s] Specific impulse
'''
class MyPredicate(PythonEnablingPredicate): # MyPredicate class definition
def eventIsEnabled(self, state, detector, g): # eventIsEnabled method definition
global enabledEvent
statedate = state.getDate() # get state date
mindate = firingDate # set min date
maxdate = firingDate.shiftedBy(state.getKeplerianPeriod()) # set max date (min date + 1 orbit period)
cnd1 = statedate.isAfterOrEqualTo(mindate) # set condition 1 (state date >= min date)
cnd2 = statedate.isBefore(maxdate) # set condition 2 (state date < max date)
if mantype == 'AN': # If maneuver shall be executed at Ascending Node:
Vz = state.getPVCoordinates(inertialFrame).getVelocity().getZ() # get state Vz
cnd3 = Vz > 0.0 # set condition 3 (state positive Vz)
elif mantype == 'DN': # else (descending node):
Vz = state.getPVCoordinates(inertialFrame).getVelocity().getZ() # get state Vz
cnd3 = Vz < 0.0 # set condition 3 (state negative Vz)
elif mantype == 'A': # If maneuver shall be executed at Apogee:
SMA = state.getA() # get SMA
r = state.getPVCoordinates(inertialFrame).getPosition().getNorm() # get position radius r
cnd3 = r > SMA # true if radius is greater than SMA
elif mantype == 'P': # If maneuver shall be executed at Perigee:
SMA = state.getA() # get SMA
r = state.getPVCoordinates(inertialFrame).getPosition().getNorm() # get position radius r
cnd3 = r < SMA # true if radius is smaller than SMA
else:
raise Exception('{:s} is an unknown maneuver type. Known options are:\n''AN'' Ascending Node\n''DN'' Descending Node\n''A'' @ Apogee\n''P'' @ Perigee\n''T'' Time triggered'.format(mantype)) # raise exception
enabledEvent = cnd1 and cnd2 and cnd3
kep = KeplerianOrbit(state.getPVCoordinates(), state.getFrame(), state.getDate(), Constants.EIGEN5C_EARTH_MU) # Keplerian orbit from state
kTA = FastMath.toDegrees(kep.getTrueAnomaly()) # [deg] True Anomaly
kMAdot = FastMath.toDegrees(kep.getMeanAnomalyDot())*86400.0/360.0 # [rev/d] Mean Anomaly dot
print('[dbg] {:s} Searching: {:s} TA: {:6.1f} MAdot[rev/d]: {:6.1f} dotProd[km^2/s]: {:5.1f} cnd: {:d} {:d} {:d}'.format(state.getDate().toString()[0:dtsl],mantype,kTA,kMAdot,g*1e-6,cnd1,cnd2,cnd3))
return enabledEvent # return boolean: true if 3 conditions are met
class ApsideStopper(PythonEventHandler):
''' Eventhandler that stops at apside '''
def init(self, initialstate, target, detector):
pass
def eventOccurred(self, s, detector, increasing):
global StopDetecting
kep = KeplerianOrbit(s.getPVCoordinates(), s.getFrame(), s.getDate(), Constants.EIGEN5C_EARTH_MU)
kTA = FastMath.toDegrees(kep.getTrueAnomaly()) # [deg] True Anomaly
kMAdot = FastMath.toDegrees(kep.getMeanAnomalyDot())*86400.0/360.0 # [rev/day] Mean Anomaly dot
evDate = s.getDate().toString()[0:dtsl] # [str] event date
if (mantype == 'P' and abs(kTA) < TAtol) or \
(mantype == 'A' and abs(abs(kTA)-180.0) < TAtol):
StopDetecting = True
ation2return = Action.STOP
else:
ation2return = Action.CONTINUE
print('Time: {:s} type: {:s} TA: {:6.1f} MAdot[rev/d]: {:6.1f} increasing: {:d} action: {}'.format(evDate,mantype,kTA,kMAdot,increasing,ation2return))
return ation2return
def resetState(self, detector, oldState):
return oldState
class FirstEvenDetector(PythonAbstractDetector):
''' Detects only the first event '''
def __init__(self, orbit, newMaxCheck=None, newThreshHold=None, newMaxIter=None, handler=None):
self.orbit = orbit
self.firstIncreasingSeen = False
global StopDetecting, enabledEvent
StopDetecting = False
if newMaxCheck is None:
# The orbit is used only to set an upper bound for the max check interval to period/3
dmax = orbit.getKeplerianPeriod() / 3
else:
dmax = newMaxCheck
if newThreshHold is None:
# The orbit is used only to set the convergence threshold according to orbit size
dthresh = 1.0e-13 * orbit.getKeplerianPeriod()
else:
dthresh = newThreshHold
if newMaxIter is None:
dmaxiter = PythonAbstractDetector.DEFAULT_MAX_ITER
else:
dmaxiter = newMaxIter
if handler is None:
handler = StopOnEvent().of_(FirstEvenDetector)
super(FirstEvenDetector, self).__init__(dmax, dthresh, dmaxiter, handler) #super(maxCheck, threshold, maxIter, handler);
def init(self, *args, **kw):
pass
def setFirstIncreasingSeen(self, value):
self.firstIncreasingSeen = value
def g(self, s):
if StopDetecting:
print('[dbg] {:s} detection is stopped'.format(s.getDate().toString()[0:dtsl]))
return 1.0
else:
dotProd = s.getPVCoordinates().getPosition().dotProduct(s.getPVCoordinates().getVelocity()) # [m^2/s] positon velocity dot product
return dotProd
def create(self, newMaxCheck, newThreshHold, newMaxIter, newHandler):
return FirstEvenDetector(self.orbit, handler=newHandler)
if mantype == 'T': # If this is a Time triggered maneuver:
triggers = DateDetector(firingDate) # define trigger as date detector at firing date
elif mantype == 'A' or mantype == 'P': # If this is an Apogee or Perigee triggered maneuver:
detector = FirstEvenDetector(self.pv.getOrbit())
mystopper = ApsideStopper().of_(FirstEvenDetector)
detector = detector.withHandler(mystopper)
predicate = MyPredicate() # define predicate object
triggers = EventEnablingPredicateFilter(detector, predicate) # define trigger as first apside detector stopper with predicate
elif mantype == 'AN' or mantype == 'DN': # If this is an out-of-plane maneuver to be executed @ AN or DN:
ascendingNodeStopper = NodeDetector(inertialFrame) # stop at ascending nodes and continue at descending nodes
predicate = MyPredicate() # define predicate object
triggers = EventEnablingPredicateFilter(ascendingNodeStopper, predicate) # define trigger as ascending node stopper with predicate
else: # This is an unknown maneuver:
raise Exception('{:s} is an unknown maneuver type. Known options are:\n''AN'' Ascending Node\n''DN'' Descending Node\n''A'' @ Apogee\n''P'' @ Perigee\n''T'' Time triggered'.format(mantype)) # raise exception
ManLoggedDectector = self.ManLogger.monitorDetector(triggers) # init. monitor detector from triggers
self.prop.addEventDetector(ManLoggedDectector) # add monitor detector to propagator
attitudeOverride = LofOffset(inertialFrame,LOFType.VNC) # Align attitude with VNC frame
impMan = ImpulseManeuver(triggers, attitudeOverride, mandv, isp) # define impulse maneuver
self.prop.addEventDetector(impMan) # add maneuver to the propagator with event detector method
Do you see any issue in it?
Since I didn’t know how to set the attribute self.firstIncreasingSeen in class ApsideStopper, I had to use the global variable StopDetecting in order to communicate to class FirstEvenDetector that the event was found. Is there a more elegant way to set that flag?
Thank you very much.
Kind regards,
Alfredo