Hi all,
I’ve been trying to count the amount of time that the satellite has spent in each eclipse in order to generate an event when the cumulative time in the current eclipse reaches a certain threshold. I’ve used this post as a starting guide and ultimately have been able to come up with the below minimum working example, which seems to be counting the eclipse durations accurately. However, I don’t know how to use it in order to trigger events after the ElapsedTimeCounter.elapsedSec
reaches the ElapsedTimeCounter.limit
. I suppose this is related to the ElapsedTimeCounter.apply
method introduced in the above post but I simply don’t know how to leverage this. Any suggestions?
Also, a bunch of additional questions that arose while I was working on this MWE and I haven’t been able to find answers to (marked with #TODO
in the MWE):
- The order of adding event detectors and handlers seems to sometimes affect the results that one gets. Is there any “golden rule” on how to do this safely to guarantee robust results? Does the order as in the
#%% Configure the propagator events
cell below make sense? - Why if I move
eclipseDetectorWithPassCounter
up or down in the#%% Configure the propagator events
cell the number of events found byloggedDetector
changes? - Why are the
PassCounter.init
andElapsedTimeCounter.init
methods never called? When exactly would they be called?
import numpy
MAX_ECLIPSES = 5 # Max no. eclipses after which propagation will be stopped.
#%% Start Orekit and configure frames & bodies
from orekit.pyhelpers import download_orekit_data_curdir, setup_orekit_curdir
if not os.path.isfile('orekit-data.zip'):
download_orekit_data_curdir()
setup_orekit_curdir()
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions
from org.orekit.frames import FramesFactory
from org.orekit.bodies import OneAxisEllipsoid, CelestialBodyFactory
from orekit import JArray_double
from org.orekit.orbits import OrbitType
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.python import PythonEventHandler
from org.hipparchus.ode.events import Action
from org.orekit.propagation.events import EclipseDetector, EventsLogger
from org.orekit.propagation.events.handlers import ContinueOnEvent
utc = TimeScalesFactory.getUTC()
j2000 = FramesFactory.getEME2000()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING, itrf)
sun = CelestialBodyFactory.getSun()
#%% Custom event handlers
class PassCounter(PythonEventHandler):
""" Counts events. Will stop propagation after first n events."""
def __init__(self, limit):
self.limit = limit
self.passesInc = 0
self.passesDec = 0
super(PassCounter, self).__init__()
def init(self, initialState, target):
# Not used, just to try and port functionality from ElapsedTimeCounter.
self.startEpoch = initialState.getDate() #TODO 'PassCounter' object has no attribute 'startEpoch'?
self.elapsedSec = 0.0 #TODO: 'PassCounter' object has no attribute 'elapsedSec'?
super(PassCounter, self).init() # Doesn't affect the odd attribute errors.
#TODO for whatever reason, init is never called. What does this method do exactly?
def eventOccurred(self, s, T, increasing):
if increasing:
self.passesInc = self.passesInc + 1
print(f'Increasing pass at {s.getDate()}, passes = {self.passesInc}')
else:
self.passesDec = self.passesDec + 1
print(f'Decreasing pass at {s.getDate()}, passes = {self.passesDec}')
# Just as an example, stop propagation after some no. events.
if self.passesInc >= self.limit and self.passesDec >= self.limit:
return Action.STOP
else:
return Action.CONTINUE
def resetState(self, detector, oldState):
return oldState
class ElapsedTimeCounter(PythonEventHandler):
""" Counts the elapsed time since a decreasing event, and triggers and event
when the elapsed time reaches a predefined limit. Inspired by this:
https://forum.orekit.org/t/eventdetection-based-on-time-in/293
"""
def __init__(self, limit, initialState):
self.limit = limit
self.startEpoch = initialState.getDate() # Start counting from this epoch.
self.elapsedSec = 0.0
super(ElapsedTimeCounter, self).__init__()
def init(self, initialState, target):
pass
def eventOccurred(self, s, T, increasing):
if not increasing: # Eclipse start.
self.startEpoch = s.getDate()
self.elapsedSec = 0.0
else: # increasing => eclipse stop.
self.elapsedSec += s.getDate().durationFrom(self.startEpoch)
print(f'Elapsed {self.elapsedSec} sec since {self.startEpoch}')
self.startEpoch = None
return Action.CONTINUE
def apply(self, s):
""" Defines g function of a new Functional detector. """
print('apply') #TODO This is never called, how to trigger events from here?
if self.startEpoch is None:
return self.elapsedSec - self.limit
else:
return (self.elapsedSec +
s.getDate().durationFrom(self.startEpoch)) - self.limit
def resetState(self, detector, oldState):
return oldState
#%% Initial state of the satellite
ra = 500 * 1000. # Apogee
rp = 400 * 1000. # Perigee
i = float(numpy.deg2rad(87.0)) # inclination
omega = float(numpy.deg2rad(20.0)) # argument of perigee
raan = float(numpy.deg2rad(10.0)) # right ascension of ascending node
ta = float(numpy.deg2rad(0.0)) # True anomaly
epoch = AbsoluteDate(2020, 1, 1, 0, 0, 00.000, utc)
a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a
initialMass = 1.0 # kg
## Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, ta, PositionAngle.TRUE,
j2000, epoch, Constants.WGS84_EARTH_MU)
#%% Configure the propagator
minIntegStep = 0.0001
maxIntegStep = 1000.0
initIntegStep = 1.0
posTolerance = 0.1
orbitType = OrbitType.CARTESIAN
tol = NumericalPropagator.tolerances(posTolerance, initialOrbit, orbitType)
# Setup the integrator.
integrator = DormandPrince853Integrator(minIntegStep,
maxIntegStep,
JArray_double.cast_(tol[0]),
JArray_double.cast_(tol[1]))
integrator.setInitialStepSize(initIntegStep)
# Setup the propagator.
propagator = NumericalPropagator(integrator)
propagator.setOrbitType(orbitType)
# Set the initial state.
initialState = SpacecraftState(initialOrbit, initialMass)
propagator.setInitialState(initialState)
#%% Configure the propagator events
passCounter = PassCounter(limit=MAX_ECLIPSES)
timeElapsedCounter = ElapsedTimeCounter(100.0, initialState)
logger = EventsLogger()
eclipseDetectorWithElapsedCounter = EclipseDetector(sun, Constants.SUN_RADIUS, earth).withUmbra().withHandler(timeElapsedCounter)
propagator.addEventDetector(eclipseDetectorWithElapsedCounter)
#TODO eclipseDetectorWithPassCounter has to be here for eclipseDetectorWithLogger to work properly. Otherwise, we get 9 logged events. Why?
eclipseDetectorWithPassCounter = EclipseDetector(sun, Constants.SUN_RADIUS, earth).withUmbra().withHandler(passCounter)
propagator.addEventDetector(eclipseDetectorWithPassCounter)
eclipseDetectorWithLogger = EclipseDetector(sun, Constants.SUN_RADIUS, earth).withUmbra().withHandler(ContinueOnEvent())
loggedDetector = logger.monitorDetector(eclipseDetectorWithLogger)
propagator.addEventDetector(loggedDetector)
#%% Propagate the orbit
epochList = [epoch, epoch.shiftedBy(36000.0)]
orekitStates = []
for currentEpoch in epochList:
currentState = propagator.propagate(currentEpoch)
orekitStates.append(currentState)
propEvents = logger.getLoggedEvents()
print(f'Found {propEvents.size()/2} eclipses out of maximum no. of {MAX_ECLIPSES}.')