Station keeping using low thrust

[HugoBocc]

Hi, I have been trying to implement a station keeping script using low thrust maneuvers. Two models could be used , the ConfigurableLowThrust , and constantThrust . The issue is that, with configurableLowThrustManeuver I can’t make PythonStartStopEventTrigger work using altitude detector in order to start the maneuver at a given altitude until i reach the second altitude then stop , then keep propagating till i reach again the lower threshold altitude and again.
The error is the following one :

    propagator.propagate(final_date)
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
java.lang.RuntimeException: AttributeError
        at org.orekit.forces.maneuvers.trigger.PythonStartStopEventsTrigger.getParametersDrivers(Native Method)
        at org.orekit.forces.maneuvers.Maneuver.getParametersDrivers(Maneuver.java:275)
        at org.orekit.utils.ParameterDriversProvider.getParameters(ParameterDriversProvider.java:82)
        at org.orekit.forces.maneuvers.Maneuver.addContribution(Maneuver.java:159)
        at org.orekit.propagation.numerical.NumericalPropagator$Main.computeDerivatives(NumericalPropagator.java:979)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:794)
        at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
        at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:259)
        at org.hipparchus.ode.AbstractIntegrator.initIntegration(AbstractIntegrator.java:204)
        at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:195)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:493)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:440)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:400)

This is the current code :

simDur = 365 * 3600 * 24*1.0



initial_date = AbsoluteDate(2020, 8, 23, 12, 0, 0.0, TimeScalesFactory.getUTC())

############################ Orbit Definition ######################################
####################################################################################

ra = 500.0 * 1000  # Apogee in meters
rp = 500.0 * 1000  
i = 87.0  # Inclination (degres)
omega = 20.0  # Perigee argument in radians
raan = 10.0    # Right ascension of ascending node in radians
lv = 0.0     # True anomaly in radians
a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0    
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a

mu = Constants.EGM96_EARTH_MU 
# Setting the initial orbit : 
initialOrbit = KeplerianOrbit(a, e, math.radians(i), math.radians(omega), math.radians(raan), math.radians(lv),
                              PositionAngleType.MEAN, FramesFactory.getEME2000(), initial_date, mu)



############################ Satellite Definition ######################################
########################################################################################


satellite_mass = 6.0
area = 0.15


############################ Environment Definition ######################################
########################################################################################


gravityProvider = GravityFieldFactory.getNormalizedProvider(60, 60)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                         Constants.WGS84_EARTH_FLATTENING,
                         FramesFactory.getITRF(IERSConventions.IERS_2010, True))

inertialFrame = FramesFactory.getEME2000()
sun = CelestialBodyFactory.getSun()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)
dragCoef = 2.2
drag = IsotropicDrag(area, dragCoef)
strengthLevel = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
supportedNames = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
inputParameters_Drag = MarshallSolarActivityFutureEstimation(supportedNames, strengthLevel)
atmosphere = NRLMSISE00(inputParameters_Drag, sun, earth)
orbitType = OrbitType.CARTESIAN


# Attitude
attitude_provider = LofOffset(inertialFrame, LOFType.NTW)
attitude = attitude_provider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame())

# Initial state with attitude
initialState = SpacecraftState(initialOrbit, attitude, satellite_mass)


############################ Reentry Definition ######################################
########################################################################################################

Altitude_reentry = 200.0 * 1000.0
altitude_low=650*1000.0
altitude_high=700*1000.0
altitudeDetector =  AltitudeDetector(Altitude_reentry, earth).withHandler(StopOnEvent())


############################ Integrator and propagator Definition ######################################
########################################################################################################


# integrator
minStep = 0.00001
maxStep = 1000.0
initStep = 60.0
positionTolerance = 1.0
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit,orbitType)
integrator = DormandPrince853Integrator(minStep, maxStep,
                                        JArray_double.cast_(tolerances[0]), JArray_double.cast_(tolerances[1]))

integrator.setInitialStepSize(initStep)

# Propagator
propagator = NumericalPropagator(integrator)
propagator.addEventDetector(altitudeDetector)
# AddforceModel
propagator.addForceModel(DragForce(atmosphere, drag))
#propagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))


############################ Handler Definition ########################################################
########################################################################################################
sma_limite=68800000
# Paramètres de la manœuvre à poussée constante
thrust = 0.01  # Poussée en Newtons
isp = 300.0  # Impulsion spécifique (s)
duration = 6000.0  # Durée de la manœuvre en secondes
direction = Vector3D(1.0, 0.0, 0.0)  # Direction de la poussée (dans le référentiel LVLH)
prop_model=BasicConstantThrustPropulsionModel(thrust,isp,direction,"low thrust")

class HandlerForDate(PythonOrekitFixedStepHandler):


    def init(self, so, t, step):
        super().__init__()
        self.smas = []
        self.dates = []
        self.incs = []
        self.eccs = []


    def handleStep(self, state):
        orbit = OrbitType.KEPLERIAN.convertType(state.getOrbit())
        sma = orbit.getA()
        inc = orbit.getI()
        date = state.getDate()
        ecc = state.getE()
        #print(sma)
        print(date)
        self.smas.append(sma / 1000 - Constants.WGS84_EARTH_EQUATORIAL_RADIUS / 1000)  # Convert to km
        self.dates.append(date.durationFrom(initial_date) / (3600 * 24 * 365))  # Convert to years
        self.incs.append(math.degrees(inc))
        self.eccs.append(ecc)
        maneuver = ConstantThrustManeuver(date, duration,attitude_provider, prop_model)
        print(sma,sma_limite)
      
           
    def finish(self, finalState):
        orbit = KeplerianOrbit(finalState.getOrbit())


time_step = 1000.0
#Handler_Date = HandlerForDate()
#Propagator.cast_(propagator).setStepHandler(time_step,Handler_Date)

direction = Vector3D(float(1), float(0), float(0))
thrust = 0.06
thruster_axis = Vector3D(float(1), float(0), float(0))

thrust_direction_provider = ConstantThrustDirectionProvider(direction)
thrust_direction_and_attitude_provider = ThrustDirectionAndAttitudeProvider.buildFromDirectionInLOF(
    LOFType.QSW,
    thrust_direction_provider,
    thruster_axis)
trigger_start = altitudeDetector =  AltitudeDetector(altitude_low, earth).withHandler(StopOnEvent())
#trigger_start.init(propagator.getInitialState(), date_end)
trigger_end = altitudeDetector =  AltitudeDetector(altitude_high, earth).withHandler(StopOnEvent())
#trigger_end.init(propagator.getInitialState(), date_end)

trigg=PythonStartStopEventsTrigger(trigger_start,trigger_end)
configurable_low_thrust_maneuver = ConfigurableLowThrustManeuver(thrust_direction_and_attitude_provider,
                                                                 trigg,
                                                                 thrust,
                                                                 isp)
propagator.addForceModel(configurable_low_thrust_maneuver)


############################ Propagate #################################################################
########################################################################################################


propagator.setOrbitType(orbitType)
propagator.setInitialState(initialState)
#final_state = propagator.propagate(initial_date, initial_date.shiftedBy(simDur))
final_date=initial_date.shiftedBy(simDur)
propagator.propagate(final_date)

Thanks in advance for you answers

[Emilien_mrlt]

I have the same issue except that i’m using impulsive maneuvers, don’t know how to handle it …

[Emilien_mrlt]

I believe this issue comes from the subclass process of StartStopEventTrigger.

In the java doc, convertStartDetector and convertStopDetector are sometimes used, meaning your trigger aren’t considered by configurableLowThrustManeuver.

Would someone agree ?

Cheers,

[Vincent]

Hello @HugoBocc and welcome to the orekit forum @Emilien_mrlt,

This message to let you know that i’m starting to make an example to see what’s wrong. I’ll keep you updated.

Cheers,
Vincent

[Vincent]

Ok so after stumbling on the same error, i realized that this is because the StartStopEventsTrigger is an abstract class and so, must be extended to be used.

Coming back when my example is ready

[HugoBocc]

Hi vincent, really appreciate you help. Thanks for your time

[Vincent]

Ok it has been a long time since i last used the Python wrapper but here is the example below.

I used a ridiculously high drag coefficient so that we can see the maneuvers in action in a single day .

Description

This example shows the use of a custom ConfigurableLowThrustManeuver used with a custom semi major apsis event detector (called SMADetector within the script).

The maneuver is activated if the semi major apsis is below 490 km and stopped when above 510 km.

Beware : You need to change ‘PATH_TO_OREKIT_DATA’ to your orekit data directory. Some may use “orekit-data.zip” as well.

Code

import math

import matplotlib.pyplot as plt
import numpy as np
import orekit
from java.util import Collections
from orekit import JArray_double
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.attitudes import LofOffset
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.forces.drag import DragForce
from org.orekit.forces.drag import IsotropicDrag
from org.orekit.forces.maneuvers import ConfigurableLowThrustManeuver
from org.orekit.forces.maneuvers.propulsion import ThrustDirectionAndAttitudeProvider
from org.orekit.forces.maneuvers.trigger import PythonStartStopEventsTrigger
from org.orekit.frames import FramesFactory
from org.orekit.frames import LOFType
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.orbits import KeplerianOrbit
from org.orekit.orbits import OrbitType
from org.orekit.orbits import PositionAngleType
from org.orekit.propagation import Propagator
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.events import
from org.orekit.propagation.events import AbstractDetector
from org.orekit.propagation.events import BooleanDetector
from org.orekit.propagation.events import PythonAbstractDetector
from org.orekit.propagation.events.handlers import ContinueOnEvent
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate
from org.orekit.time import TimeScalesFactory
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions

vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir

setup_orekit_curdir('PATH_TO_OREKIT_DATA')

### CONSTANTS
UTC = TimeScalesFactory.getUTC()
INERTIAL = FramesFactory.getEME2000()
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
EARTH_BODYSHAPE = ReferenceEllipsoid.getIers2010(ITRF)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
SUN = CelestialBodyFactory.getSun()
INITIAL_DATE = AbsoluteDate(2024, 8, 30, 15, 0, 00.000, UTC)


### CLASSES
class SMAHandler(PythonOrekitFixedStepHandler):
    def __init__(self):
        super().__init__()
        self.values = []
        self.dates = []

    def init(self, spacecraftState, absoluteDate, double):
        pass

    def handleStep(self, spacecraftState):
        self.values.append(spacecraftState.getA())
        self.dates.append(spacecraftState.getDate().durationFrom(INITIAL_DATE))

    def finish(self, spacecraftState):
        pass


class ActualStartStopTriggerEvents(PythonStartStopEventsTrigger):
    def __init__(self, start_detector, stop_detector):
        super().__init__(start_detector, stop_detector)

    def getParametersDrivers(self):
        # Needs to return an empty java list so that the code does not fail !
        return Collections.emptyList()


class SMADetector(PythonAbstractDetector):

    def __init__(self, sma_threshold,
                 max_check=AbstractDetector.DEFAULT_MAXCHECK,
                 threshold=AbstractDetector.DEFAULT_THRESHOLD,
                 max_iter=AbstractDetector.DEFAULT_MAX_ITER,
                 event_handler=ContinueOnEvent()):
        super().__init__(max_check, threshold, max_iter, event_handler)
        self.sma_threshold = sma_threshold

    def g(self, spacecraftState):
        return spacecraftState.getA() - self.sma_threshold

    def create(self, max_check, threshold, max_iter, event_handler):
        return SMADetector(self.sma_threshold, max_check.currentInterval(None), threshold, max_iter, event_handler)


### METHODS
def getNumericalPropagator(initialOrbit):
    minStep = 0.001
    maxstep = 1000.0
    initStep = 60.0
    positionTolerance = 1.0

    tolerances = NumericalPropagator.tolerances(positionTolerance,
                                                initialOrbit,
                                                initialOrbit.getType())

    integrator = DormandPrince853Integrator(minStep, maxstep,

                                            JArray_double.cast_(tolerances[0]),
                                            # Double array of doubles needs to be casted in Python
                                            JArray_double.cast_(tolerances[1]))
    integrator.setInitialStepSize(initStep)

    return NumericalPropagator(integrator)


def getOrbit():
    # Define orbit
    ra = 500 * 1000  # Apogee
    rp = 498 * 1000  # Perigee
    i = math.radians(87.0)  # inclination
    omega = math.radians(20.0)  # perigee argument
    raan = math.radians(10.0)  # right ascension of ascending node
    lv = math.radians(0.0)  # True anomaly

    a = (rp + ra + 2 * EARTH_RADIUS) / 2.0
    e = 1.0 - (rp + EARTH_RADIUS) / a

    return KeplerianOrbit(a, e, i, omega, raan, lv,
                          PositionAngleType.TRUE,
                          INERTIAL, INITIAL_DATE, Constants.WGS84_EARTH_MU)


def getSMABasedManeuverTriggers():
    """
    Create a Semi Major Apsis based maneuver triggers which START maneuver when sma is below 490 km + earth radius
      and STOP maneuver when sma is above 510 km + earth radius.
    :return: Altitude based maneuver trigger
    """
    # Define start detector
    start_detector = BooleanDetector.notCombine(SMADetector(EARTH_RADIUS + 490 * 1000).withMaxCheck(float(10)))

    # Define stop detector
    stop_detector = (SMADetector(EARTH_RADIUS + 510 * 1000).withMaxCheck(float(10)))

    return ActualStartStopTriggerEvents(AbstractDetector.cast_(start_detector), AbstractDetector.cast_(stop_detector))


def getSMABasedManeuvers():
    # Define thrust direction provider (thrust along velocity axis to increase sma)
    thrust_direction_provider = (
        ThrustDirectionAndAttitudeProvider.buildFromCustomAttitude(LofOffset(INERTIAL, LOFType.TNW), Vector3D.PLUS_I))

    # Define maneuver triggers
    maneuver_triggers = getSMABasedManeuverTriggers()

    # Define thrust
    thrust = float(1)  # N, High thrust on purpose

    # Define isp
    isp = float(2000)  # s, High isp on purpose

    return ConfigurableLowThrustManeuver(thrust_direction_provider, maneuver_triggers, thrust, isp)


def getDragForce():
    # Get atmosphere
    atmosphere = NRLMSISE00(CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES), SUN, EARTH_BODYSHAPE)

    # Get drag sensitive (ridiculously high drag coeff on purpose)
    cross_section = float(2)
    drag_coeff = float(2000)
    drag_sensitive = IsotropicDrag(cross_section, drag_coeff)

    return DragForce(atmosphere, drag_sensitive)


### MAIN
if __name__ == '__main__':
    ## Get initial orbit
    initialOrbit = getOrbit()

    # Get Numerical propagator
    propagator_num = getNumericalPropagator(initialOrbit)

    # Create initial state
    satellite_mass = 1000.0  # kg.
    initialState = SpacecraftState(initialOrbit, satellite_mass)

    # Set orbit type & initial state
    propagator_num.setOrbitType(OrbitType.CARTESIAN)
    propagator_num.setInitialState(initialState)

    # Add drag so that the orbit loses energy
    drag_force = getDragForce()
    propagator_num.addForceModel(drag_force)

    # Add configurable low thrust
    sma_maneuvers = getSMABasedManeuvers()
    propagator_num.addForceModel(sma_maneuvers)

    # Add custom handler to register semi major axis evolution throughout the propagation
    sma_handler = SMAHandler()
    Propagator.cast_(propagator_num).setStepHandler(float(10), sma_handler)

    # Propagate for one day
    propagator_num.propagate(initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY))

    # Plot semi major axis minus earth radius to see maneuver effects
    fig, ax = plt.subplots()
    ax.scatter(np.array(sma_handler.dates) / 3600, (np.array(sma_handler.values) - EARTH_RADIUS) / 1000)
    ax.set(xlabel='time since start (hour)', ylabel='SMA - Earth radius (km)')
    ax.grid()
    plt.show()

Plot

Below is the plot for those who do not have matplotlib :

[HugoBocc]

Thanks a lot for your help it’s working perfectly fine. It’s nice to see how involved is the orekit community. Have a nice day

[Emilien_mrlt]

Hello everyone, it’s a pleasure to see how reactive the Orekit community is !

I combined the SMA based trigger with a MissionDuration event which desable any maneuvers after the propagator detect the end of mission date, which is working great. (Use of the Boolean.andCombine([SMADetector, MissionDurationDetector])

To do so, I create a subclass exactly like SMADetector.

Now I’m trying to create the maneuver only when the spacecraft isn’t in Eclipse in order to Combine the start maneuver which will be triggered by (SMA detector and MissionDurationDetector and notInEclipse.) In the same way, my stoptrigger is now
Trigger_start = BooleanDetector.andCombine([sma_start_detector,eclipse_detect_penumbra, MissionDuration_detector])
Trigger_stop = BooleanDetector.orCombine([sma_stop_detector, eclipse_detect_umbra])

However, the eclipse aren’t detected. I think this is due to the fact that EclipseDetector isn’t an abstractDetector (but MissionDurationDetector and SMADetector are AbstractDetector). So it seems that casting detector that aren’t the same type is not possible right ?
Would you have any idea on how to handle it ?

I have a code if needed to let you know how I combined this EventsDetector.

Thanks for your help !

[Vincent]

Hello @Emilien_mrlt,

First i think you should use this :

Trigger_start = BooleanDetector.andCombine([sma_start_detector,eclipse_detect_penumbra, MissionDuration_detector])
Trigger_stop = BooleanDetector.andCombine([sma_stop_detector, eclipse_detect_penumbra])

@Emilien_mrlt

However, the eclipse aren’t detected

Do you mean that the maneuvers do not even start ?

@Emilien_mrlt

I think this is due to the fact that EclipseDetector isn’t an abstractDetector (but MissionDurationDetector and SMADetector are AbstractDetector)

EclipseDetector is indeed an AbstractDetector so casting should work.

Keep us updated !

Cheers,
Vincent

[Emilien_mrlt]

Hi Vincent,
To explain a bit, the following code aims to create station keeping maneuver between 480km and 500km. The mission duration is here setted to be 3 years and the propagation duration 4 years (meaning that the maneuver should not be triggered after 3 years.)
Without any Eclipse duration, here is the result, which I think is working nice.

Now I’m trying to allow these differents raising SMA maneuvers only in when the satellite is not in eclipse, by adding EclipseDetector to the trigger maneuver and to stop it too. Only a natural decay is allowed in umbra.
the code is the following :

def getSMABasedManeuverTriggers():
        # Define detectors
        sma_start_detector = BooleanDetector.notCombine(SMA_f_Detector(SK_SMA_LINF).withMaxCheck(float(600)))    # Detect when the lower altitude is reached (which will begin the maneuver)
        sma_stop_detector = BooleanDetector.notCombine(SMA_f_Detector(SK_SMA_SUP).withMaxCheck(float(600)))                                 # Detect when the upper altitude is reached (which will end the maneuver)

        # Defin Eclipse detector
        sun = CelestialBodyFactory.getSun()
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
        earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)

        Eclipse_Detect_Penumbra = BooleanDetector.notCombine(EclipseDetector(sun, Constants.SUN_RADIUS, earth).withPenumbra())
        Eclipse_Detect_Umbra = BooleanDetector.notCombine(EclipseDetector(sun, Constants.SUN_RADIUS, earth).withUmbra())
        # Define mission duration detector
        mission_duration_detector = BooleanDetector.notCombine(MissionDurationDetector(initial_date.shiftedBy(MD)).withMaxCheck(float(10))) # Detect when the end of mission date is passed by the propagator.
        
        # Cast Detectors
        sma_start_detector = AbstractDetector.cast_(sma_start_detector)
        mission_duration_detector = AbstractDetector.cast_(mission_duration_detector)
        sma_stop_detector = AbstractDetector.cast_(sma_stop_detector)   
        Eclipse_Detect_Penumbra = AbstractDetector.cast_(Eclipse_Detect_Penumbra)
        Eclipse_Detect_Umbra = AbstractDetector.cast_(Eclipse_Detect_Umbra)

        # Define detectors for logging events based on SMA limits
        First_Possible_trigger = BooleanDetector.andCombine([sma_start_detector, mission_duration_detector, Eclipse_Detect_Penumbra])
        Second_Possible_trigger = BooleanDetector.andCombine([sma_start_detector, mission_duration_detector,Eclipse_Detect_Penumbra])
        
        # Combine detectors: Only execute SMA maneuvers if within mission duration

        Detector_events_logg_low  =  BooleanDetector.notCombine(SMA_f_Detector(SK_SMA_LINF))
        Mission_duration_event_logg_low = AbstractDetector.cast_(BooleanDetector.notCombine(MissionDurationDetector(initial_date.shiftedBy(MD)).withMaxCheck(float(10))))
        combined_start_detector_loggs = BooleanDetector.andCombine([Detector_events_logg_low, Mission_duration_event_logg_low])
        eventsLogger_lower = EventsLogger()     # Logger to count and record maneuver start events
        eventsLogger_upper = EventsLogger()     # Logger to count and record maneuver end events
        ev_logs_lower = eventsLogger_lower.monitorDetector(combined_start_detector_loggs)  # Monitor events and store them in the logs
        ev_logs_upper = eventsLogger_upper.monitorDetector(SMA_f_Detector(SK_SMA_SUP).withMaxCheck(float(600)))
        
        

        # Stop Trigger: Since my station keeping is defined by severals finite burns to raise the SMA, the finite maneuvers end when there is an eclipse or when the Upper bound of the SMA is reached. 
        Trigger_stop = BooleanDetector.orCombine([sma_stop_detector, Eclipse_Detect_Umbra])
        start_trigger=BooleanDetector.andCombine([First_Possible_trigger, Second_Possible_trigger])
        return ev_logs_lower, eventsLogger_lower, ev_logs_upper, eventsLogger_upper, ActualStartStopTriggerEvents(AbstractDetector.cast_(start_trigger), AbstractDetector.cast_(Trigger_stop))

    def getSMABasedManeuvers():
        thrust_direction_provider = (
            ThrustDirectionAndAttitudeProvider.buildFromCustomAttitude(LofOffset(inertialFrame, LOFType.TNW), Vector3D.PLUS_I)) # Define thrust direction provider for the maneuver

        ev_logs, eventsLogger, ev_logs_upper, eventslogger_upper, maneuver_triggers = getSMABasedManeuverTriggers() # Retrieve event logs, loggers, and maneuver triggers

        return ev_logs, eventsLogger,ev_logs_upper, eventslogger_upper, ConfigurableLowThrustManeuver(thrust_direction_provider, maneuver_triggers, Thrust, ISP)
        # Creation of the low thrust maneuver
    return getSMABasedManeuvers()

The logs are used here to get several informations on the maneuvers (date, duration, …)
I obtained the following output

By zooming, the maneuver isn’t triggered at 480 km (a bit under this threshold) and as you can see, the eclipseDetector isn’t working neither

I divided the start trigger in two separates one because if the second one should be triggered only if the first is (The other maneuvers should only take in consideration mission duration and EclipseDetector, and not the sma Threshold which is the case in here …)
I know this maneuver is a bit tricky to understand as it is described here and I really apologize for it.
the maneuver logic is the following

MANEUVER triggered by ((SMA detector AND time < mission duration AND Penumbra) OR (first maneuver triggered AND time < MD AND Penumbra) UNTIL (SMA_Upper_bound OR Umbra)

Thanks foryour help and thanks for your reply @Vincent !

[Emilien_mrlt]

Trigger_start = BooleanDetector.andCombine([sma_start_detector,eclipse_detect_penumbra, MissionDuration_detector])
Trigger_stop = BooleanDetector.andCombine([sma_stop_detector, eclipse_detect_penumbra])

By doing this, the maneuver will be triggered when all the three detector are setted to True but I want it to be true for the first one and then the maneuver should be triggered when not in eclipse until the sma upper threshold is detected (Without going back to the lower SMA threshold.)

That is why I thought to use BooleanDetector.orCombine([sma_stop_detector, eclipse_detect_penumbra]) to stop the maneuver. (not using a andCombine), but it’s indeed not working as shown in the previous post.

Thanks !

[Vincent]

@Emilien_mrlt

By zooming, the maneuver isn’t triggered at 480 km (a bit under this threshold)

This is (probably) normal as your plot shows the altitude and not the SMA assuming that you are not using a spherical Earth or an equatorial orbit.

@Emilien_mrlt

and as you can see, the eclipseDetector isn’t working neither

It’s hard to see whether it is working or not given the scale on the x axis. I would advise you to purposely increase the drag to see this kind of plots over several days only. It will be :

• Faster to simulate
• Easier to see eclipses

@Emilien_mrlt

By doing this, the maneuver will be triggered when all the three detector are setted to True but I want it to be true for the first one and then the maneuver should be triggered when not in eclipse until the sma upper threshold is detected (Without going back to the lower SMA threshold.)

Ok thank you for adding this part, it makes things much clearer. You are going to need to do something a bit more elaborate to have this kind of behaviour.

Basically the start trigger is only a trigger so you need to store somewhere the fact that it has been triggered at one point for your second maneuvers to continue their thrust until target upper SMA is reached.

You could create a copy of the start event trigger of your first maneuvers that would be wrapped inside a custom detector. The goal of this custom detector would be to switch to a positive value when the wrapped detector is triggered. This way, the second maneuvers would thrusts when outside eclipse until target upper SMA is reached.

[Emilien_mrlt]

Hi Vincent, thanks for your reply, actually it helped me solve the problem.
I switched the area to 2 m² in order to increase the drag and reduce the propagation time.

I created a custom detector that switch a flag to true when the first maneuver has been triggered.
I also made a BIG mistake (if you want to perform a maneuver when the satellite isn’t in eclipse, the start trigger is (NegateDetector(EclipseDetector(sun, Constants.SUN_RADIUS, earth).withPenumbra()).

Here is the new code that solve the problem and perform maneuvers when in sunlight, between two sma threshold and before mission duration date has occured.

def Sk_fin_burn_eclipse_Maneuvers_Definition(inertialFrame, SK_SMA_LINF, SK_SMA_SUP, MD, initial_date, Thrust, ISP):


    class Detector_first_trigger(PythonAbstractDetector):
        def __init__(self, sma_lower_threshold, sma_upper_threshold,
                    max_check=AbstractDetector.DEFAULT_MAXCHECK,
                    threshold=AbstractDetector.DEFAULT_THRESHOLD,
                    max_iter=AbstractDetector.DEFAULT_MAX_ITER,
                    event_handler=ContinueOnEvent()):
            super().__init__(max_check, threshold, max_iter, event_handler)
            self.sma_lower_threshold = sma_lower_threshold
            self.sma_upper_threshold = sma_upper_threshold

        def g(self, spacecraftState):
            global maneuverLOW_executed

            # Current semi-major axis
            current_sma = spacecraftState.getA()
            
            # If below the lower threshold and maneuver has not been executed, set the flag to True
            if current_sma < self.sma_lower_threshold and not maneuverLOW_executed:
                maneuverLOW_executed = True

            # If above the upper threshold, reset the flag
            if current_sma > self.sma_upper_threshold:
                maneuverLOW_executed = False

            # Detector's function g returns positive when maneuverLOW_executed is True
            return 1.0 if maneuverLOW_executed else -1.0

        def create(self, max_check, threshold, max_iter, event_handler):
            return Detector_first_trigger(self.sma_lower_threshold, self.sma_upper_threshold, max_check, threshold, max_iter, event_handler)



    class ActualStartStopTriggerEvents(PythonStartStopEventsTrigger):
        def __init__(self, start_detector, stop_detector):
            super().__init__(start_detector, stop_detector)

        def getParametersDrivers(self):
            return Collections.emptyList()


    class SMA_f_Detector(PythonAbstractDetector): # Class to detect when the spacecraft reaches a specific semi-major axis (SMA)
        def __init__(self, sma_threshold,
                     max_check=AbstractDetector.DEFAULT_MAXCHECK,
                     threshold=AbstractDetector.DEFAULT_THRESHOLD,
                     max_iter=AbstractDetector.DEFAULT_MAX_ITER,
                     event_handler=ContinueOnEvent()):
            super().__init__(max_check, threshold, max_iter, event_handler)
            self.sma_threshold = sma_threshold

        def g(self, spacecraftState): # Define the g function to detect events: returns the difference between the current SMA and the threshold
            return spacecraftState.getA() - self.sma_threshold

        def create(self, max_check, threshold, max_iter, event_handler): # Create a new instance of the detector with the same threshold and modified parameters
            return SMA_f_Detector(self.sma_threshold, max_check.currentInterval(None), threshold, max_iter, event_handler)


    class MissionDurationDetector(PythonAbstractDetector):          #This class is used to detect the date when the elapsed time is the mission duration
        def __init__(self, max_date, 
                     max_check=AbstractDetector.DEFAULT_MAXCHECK,
                     threshold=AbstractDetector.DEFAULT_THRESHOLD,
                     max_iter=AbstractDetector.DEFAULT_MAX_ITER,
                     event_handler=StopOnEvent()):
            super().__init__(max_check, threshold, max_iter, event_handler)
            self.max_date = max_date

        def g(self, spacecraftState):
            # Calculate the time difference between current spacecraft state date and the max_date
            current_date = spacecraftState.getDate()
            return current_date.durationFrom(self.max_date)

        def create(self, max_check, threshold, max_iter, event_handler):
            return MissionDurationDetector(self.max_date, max_check, threshold, max_iter, event_handler)
        
        
    def getSMABasedManeuverTriggers():

        # Define detectors
        sma_start_detector = BooleanDetector.notCombine(SMA_f_Detector(SK_SMA_LINF)).withMaxCheck(100.0)                            # Detect when the lower altitude is reached (which will begin the maneuver)
        sma_stop_detector = SMA_f_Detector(SK_SMA_SUP)                              # Detect when the upper altitude is reached (which will end the maneuver)

        detectorX = Detector_first_trigger(SK_SMA_LINF, SK_SMA_SUP)
        # Define mission duration detector
        mission_duration_detector = BooleanDetector.notCombine(MissionDurationDetector(initial_date.shiftedBy(MD))) # Detect when the end of mission date is passed by the propagator.


        sun = CelestialBodyFactory.getSun()
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
        earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)
        EclipseDetector_penumbra = BooleanDetector.notCombine(EclipseDetector(sun, Constants.SUN_RADIUS, earth).withPenumbra())


        First_trigger = BooleanDetector.andCombine([NegateDetector(EclipseDetector_penumbra), mission_duration_detector, sma_start_detector])
        Second_trigger = BooleanDetector.andCombine([NegateDetector(EclipseDetector_penumbra),mission_duration_detector, detectorX])
        Maneuvers_trigger = BooleanDetector.orCombine([First_trigger,Second_trigger])

        Stop_trigger = BooleanDetector.orCombine([EclipseDetector_penumbra, sma_stop_detector])


        eventsLogger_Exiting_Eclipse = EventsLogger()     # Logger to count and record maneuver start events
        eventsLogger_Entering_Eclipse = EventsLogger()     # Logger to count and record maneuver start events
        eventsLogger_Starting_maneuver = EventsLogger()
        eventsLogger_Ending_maneuver = EventsLogger()
        ev_logs_exiting_eclipse = eventsLogger_Exiting_Eclipse.monitorDetector(BooleanDetector.andCombine([EclipseDetector_penumbra,detectorX]))  
        ev_logs_entering_eclipse = eventsLogger_Entering_Eclipse.monitorDetector(BooleanDetector.andCombine([NegateDetector(EclipseDetector_penumbra),detectorX]))
        ev_logs_starting_maneuver = eventsLogger_Starting_maneuver.monitorDetector(BooleanDetector.notCombine(First_trigger.withMaxCheck(10.0)))
        #ev_logs_ending_maneuver = eventsLogger_Ending_maneuver.monitorDetector(BooleanDetector.notCombine(sma_stop_detector))
        ev_logs_ending_maneuver = eventsLogger_Ending_maneuver.monitorDetector(BooleanDetector.andCombine([sma_stop_detector,NegateDetector(EclipseDetector_penumbra)]))

        return (ev_logs_exiting_eclipse, eventsLogger_Exiting_Eclipse, ev_logs_entering_eclipse, eventsLogger_Entering_Eclipse,
                ev_logs_starting_maneuver,eventsLogger_Starting_maneuver, ev_logs_ending_maneuver, eventsLogger_Ending_maneuver,
                ActualStartStopTriggerEvents(AbstractDetector.cast_(Maneuvers_trigger), AbstractDetector.cast_(Stop_trigger)))

    def getSMABasedManeuvers():
        thrust_direction_provider = (
            ThrustDirectionAndAttitudeProvider.buildFromCustomAttitude(LofOffset(inertialFrame, LOFType.TNW), Vector3D.PLUS_I)) # Define thrust direction provider for the maneuver
        


        (ev_logs_exiting_eclipse, eventsLogger_Exiting_Eclipse, ev_logs_entering_eclipse, eventsLogger_Entering_Eclipse,ev_logs_starting_maneuver,eventsLogger_Starting_maneuver,
        ev_logs_ending_maneuver, eventsLogger_Ending_maneuver, maneuver_triggers)= getSMABasedManeuverTriggers() # Retrieve event logs, loggers, and maneuver triggers
        

        return (ev_logs_exiting_eclipse, eventsLogger_Exiting_Eclipse, ev_logs_entering_eclipse, eventsLogger_Entering_Eclipse,ev_logs_starting_maneuver,
        eventsLogger_Starting_maneuver, ev_logs_ending_maneuver, eventsLogger_Ending_maneuver, 
        ConfigurableLowThrustManeuver(thrust_direction_provider, maneuver_triggers, Thrust, ISP))
    
 
    return getSMABasedManeuvers()

Maybe this code could help other people.
Here is the result :

The eclipse and sunlight times are quite accurate with analytical calculations.

I would like to thank the Orekit community and you @Vincent for your help !

Emilien

[Vincent]

Hello @Emilien_mrlt,

Good job ! Thank you for providing the code which will be of help to other users !

Cheers,
Vincent

[dominic-mt]

Hello @Vincent ,

The above solution is perfect for people using the JCC-based Orekit package. However, for those using the orekit-jpype package, how can we create a StartStopTrigger? In the jpype version, it is not possible to inherit Java classes in Python, right? I believe creating the Java class and compiling it into a JAR file is necessary in such cases. Could you provide a similar solution using Java?

[Vincent]

Hello @dominic-mt and welcome to the Orekit forum !

@dominic-mt

In the jpype version, it is not possible to inherit Java classes in Python, right?

From my understanding of Jpype, you are right indeed.

@dominic-mt

Could you provide a similar solution using Java?

I’ll try to do so in my spare time but i cannot promise anything regarding when i will be able to do it .

Cheers,
Vincent

[dominic-mt]

Greetings @Vincent ,

Hope you’re doing great! Below is a Java class I crafted with my bare-minimum Java wizardry. In Orekit, to implement StartStopEventsTrigger, we must primarily handle two abstract methods from its base class: convertStartDetector and convertStopDetector.

The current implementation, as you’ll notice, is pretty much hardcoded—and yes, I know it’s incorrect . Could you help me with a more general and proper implementation?

Thanks in advance for being awesome!

import java.util.List;
import java.util.ArrayList;

import org.orekit.forces.maneuvers.trigger.StartStopEventsTrigger;
import org.orekit.propagation.events.AbstractDetector;
import org.orekit.propagation.events.ParameterDrivenDateIntervalDetector;
import org.orekit.propagation.events.FieldParameterDrivenDateIntervalDetector;
import org.orekit.propagation.events.FieldEventDetectionSettings;
import org.orekit.propagation.events.FieldAbstractDetector;
import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
import org.orekit.propagation.events.FieldAdaptableInterval;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.ParameterDriver;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;


public class OrekitStartStopEventsTrigger<A extends AbstractDetector<A>, O extends AbstractDetector<O>>
    extends StartStopEventsTrigger<A, O> {

    public OrekitStartStopEventsTrigger(A prototypeStartDetector, O prototypeStopDetector) {
        super(prototypeStartDetector, prototypeStopDetector);
    }


    protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>>
         FieldAbstractDetector<D, S> convertStopDetector(final Field<S> field, final O detector)
    {
        return convertDetector(field, detector);
    }

    // Need to implement this because the function is abstract in StartStopEventsTrigger
    protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>>
         FieldAbstractDetector<D, S> convertStartDetector(final Field<S> field, final A detector)
    {
        return convertDetector(field, detector);
    }

    protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>>
    FieldAbstractDetector<D, S> convertDetector(final Field<S> field, final AbstractDetector detector) {
        AbsoluteDate refstart = new AbsoluteDate(2024, 11, 19, 0, 0, 0, TimeScalesFactory.getUTC());
        AbsoluteDate refend = new AbsoluteDate(2024, 11, 19, 0, 0, 0, TimeScalesFactory.getUTC());
        return new FieldParameterDrivenDateIntervalDetector(field, "", refstart, refend);
    }


    public List<ParameterDriver> getParametersDrivers() {
        List<ParameterDriver> drivers = new ArrayList<>();
        return drivers;
    }
}

[Vincent]

Hello @dominic-mt,

I figured i would give a full Java example instead. It is slightly different from what @Emilien_mrlt did but the result is equivalent.

@dominic-mt

Hope you’re doing great! Below is a Java class I crafted with my bare-minimum Java wizardry. In Orekit, to implement StartStopEventsTrigger, we must primarily handle two abstract methods from its base class: convertStartDetector and convertStopDetector.

To answer your question, you do not really need to implement these methods as they are only used when performing operation with CalculusFieldElement (advanced usage of Orekit i would say). Please see the code below for a concrete example in Java without implementing these functions.

Otherwise you would have to implement these methods to convert input start/stop detectors into their equivalent Field class.

Logic

In my version, i decided to implement some kind of TriggeredEvent which is a wrapping event taking one raw EventDetector to be triggered and two other EventDetector to start and stop checking the g function of the raw detector. Perhaps the same things could have been done with a EventEnablingPredicateFilter but i wanted to write something simple to understand.

By doing so, the wrapped eclipse detector that is used to trigger the maneuvers is only activated when the lower SMA threshold is triggered. Once the target upper SMA threshold is reached, the wrapped eclipse detector is deactivated until lower SMA threshold is reached again.

Java code sample

BEWARE 1: Check the loadOrekitData method as your Orekit data could be located elsewhere on your system

BEWARE 2: I used a ridiculously high drag coefficient again to avoid having to propagate for weeks/months

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.events.Action;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.LofOffset;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DataProvider;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.drag.DragForce;
import org.orekit.forces.drag.DragSensitive;
import org.orekit.forces.drag.IsotropicDrag;
import org.orekit.forces.maneuvers.Maneuver;
import org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel;
import org.orekit.forces.maneuvers.propulsion.PropulsionModel;
import org.orekit.forces.maneuvers.trigger.ManeuverTriggers;
import org.orekit.forces.maneuvers.trigger.StartStopEventsTrigger;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.LOFType;
import org.orekit.models.earth.ReferenceEllipsoid;
import org.orekit.models.earth.atmosphere.Atmosphere;
import org.orekit.models.earth.atmosphere.NRLMSISE00;
import org.orekit.models.earth.atmosphere.data.CssiSpaceWeatherData;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.PositionAngleType;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.AbstractDetector;
import org.orekit.propagation.events.AdaptableInterval;
import org.orekit.propagation.events.BooleanDetector;
import org.orekit.propagation.events.EclipseDetector;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldAbstractDetector;
import org.orekit.propagation.events.handlers.ContinueOnEvent;
import org.orekit.propagation.events.handlers.EventHandler;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.ParameterDriver;

import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Collections;
import java.util.List;

public class StationKeepingNotInEclipse {

    public static void main(String[] args) {
        // Load orekit data
        loadOrekitData();

        // Initialize constants
        final Frame            GCRF  = FramesFactory.getGCRF();
        final Frame            ITRF  = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
        final OneAxisEllipsoid EARTH = ReferenceEllipsoid.getIers2010(ITRF);

        // Create some initial orbit
        final AbsoluteDate initialDate = new AbsoluteDate();

        final double a     = 6950e3;
        final double e     = 0.001;
        final double i     = FastMath.toRadians(30.0);
        final double omega = FastMath.toRadians(87.0);
        final double raan  = FastMath.toRadians(45.0);
        final double nu    = FastMath.toRadians(0);

        final KeplerianOrbit initialOrbit = new KeplerianOrbit(a,
                                                               e,
                                                               i,
                                                               omega,
                                                               raan,
                                                               nu,
                                                               PositionAngleType.TRUE,
                                                               FramesFactory.getGCRF(),
                                                               initialDate,
                                                               Constants.EIGEN5C_EARTH_MU);

        // Create propagator
        final NumericalPropagator numericalPropagator = getNumericalPropagator(initialOrbit);

        // Create initial state
        final SpacecraftState initialState = new SpacecraftState(initialOrbit, 150);
        numericalPropagator.setInitialState(initialState);

        // Add drag force
        final double    crossSection = 1;
        final double    dragCoef     = 20000; // Ridiculously high drag coefficient
        final DragForce dragForce    = createDragForce(crossSection, dragCoef, EARTH);
        numericalPropagator.addForceModel(dragForce);

        // Add station keeping maneuver
        // Create attitude provider override
        final AttitudeProvider attitudeProvider = new LofOffset(GCRF, LOFType.TNW);

        // Create maneuver trigger
        final double           lowSMAThreshold  = 6900e3;
        final double           highSMAThreshold = 7000e3;
        final ManeuverTriggers maneuverTriggers = getManeuverTriggers(lowSMAThreshold, highSMAThreshold, EARTH);

        // Create propulsion model
        final double thrust = 1;
        final double isp    = 300;
        final PropulsionModel propulsionModel =
                new BasicConstantThrustPropulsionModel(thrust, isp, Vector3D.PLUS_I, "station_keeping");

        final Maneuver stationKeepingManeuver = new Maneuver(attitudeProvider, maneuverTriggers, propulsionModel);
        numericalPropagator.addForceModel(stationKeepingManeuver);

        // Add eclipse with handler for better understanding of displayed values
        final EclipseDetector eclipseDetector = getEclipseDetector(EARTH).withHandler((s, detector, increasing) -> {
            if (increasing) {
                System.out.println("Exiting eclipse");
            } else {
                System.out.println("Entering eclipse");
            }
            return Action.CONTINUE;
        });
        numericalPropagator.addEventDetector(eclipseDetector);

        // Add step handler for semi major axis display & output during propagation
        final StringBuilder dataBuilder = new StringBuilder();
        dataBuilder.append("timeFromStart,a\n");
        numericalPropagator.setStepHandler(60, (s) -> {

            // Display date and current SMA
            System.out.format("Date: %s | Semi-major axis: %3.3f km\n", s.getDate(), s.getA() / 1e3);

            // Add data for output
            dataBuilder.append(s.getDate().durationFrom(initialDate));
            dataBuilder.append(",");
            dataBuilder.append(s.getA());
            dataBuilder.append("\n");
        });

        // Propagate
        numericalPropagator.propagate(initialDate.shiftedBy(2 * Constants.JULIAN_DAY));

        // Store results in csv file
        writeResults(dataBuilder);
    }

    private static void loadOrekitData() {

        //Loading Data
        final File rootFile   = new File(System.getProperty("user.home"));
        final File orekitData = new File(rootFile, "orekit-data");
        if (!orekitData.exists()) {
            System.out.format("File %s not found.%n", orekitData.getAbsolutePath());
        }
        final DataProvider dirCrawler = new DirectoryCrawler(orekitData);
        DataContext.getDefault().getDataProvidersManager().addProvider(dirCrawler);
    }

    private static NumericalPropagator getNumericalPropagator(Orbit initialOrbit) throws OrekitException {

        // Create integrator
        double minStep           = 0.001;
        double maxStep           = 60.;
        double positionTolerance = 1e-3;
        double[][] tolerances =
                NumericalPropagator.tolerances(positionTolerance, initialOrbit, initialOrbit.getType());
        DormandPrince853Integrator integrator =
                new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);

        // Create propagator
        return new NumericalPropagator(integrator);
    }

    private static DragForce createDragForce(final double crossSection,
                                             final double dragCoef,
                                             final BodyShape earth) {

        // Create atmosphere
        final Atmosphere atmosphere =
                new NRLMSISE00(new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES),
                               CelestialBodyFactory.getSun(), earth);

        // Create drag sensitive
        final DragSensitive dragSensitive = new IsotropicDrag(crossSection, dragCoef);


        return new DragForce(atmosphere, dragSensitive);
    }

    private static ManeuverTriggers getManeuverTriggers(final double lowSMAThreshold,
                                                        final double highSMAThreshold,
                                                        final OneAxisEllipsoid earth) {
        // Create low and high sma thresholds
        final EventDetector smaLowCrossingDetector =
                BooleanDetector.notCombine(new SemiMajorAxisCrossingDetector(lowSMAThreshold));
        final EventDetector smaHighCrossingDetector = new SemiMajorAxisCrossingDetector(highSMAThreshold);

        // Create eclipse detector
        final EventDetector rawEclipseDetector = getEclipseDetector(earth);

        // Create triggered eclipse detector
        final TriggeredEvent triggeredEclipseDetector =
                new TriggeredEvent(rawEclipseDetector, smaLowCrossingDetector, smaHighCrossingDetector);

        // Create stop event trigger
        final EventDetector reversedEclipseDetector =
                BooleanDetector.notCombine(getEclipseDetector(earth));

        final BooleanDetector stopEventTrigger =
                BooleanDetector.orCombine(smaHighCrossingDetector, reversedEclipseDetector);

        return new StationKeepingTrigger(triggeredEclipseDetector, stopEventTrigger);
    }

    private static EclipseDetector getEclipseDetector(OneAxisEllipsoid earth) {
        return new EclipseDetector(CelestialBodyFactory.getSun(),
                                   Constants.SUN_RADIUS,
                                   earth);
    }

    private static void writeResults(StringBuilder dataBuilder) {
        final File root   = new File(System.getProperty("user.home"));
        final File output = new File(root, "output.csv");
        try {
            final FileWriter fileWriter = new FileWriter(output);
            fileWriter.write(dataBuilder.toString());
            fileWriter.close();
        } catch (IOException e) {
            throw new RuntimeException(e);
        }
    }

    private static class SemiMajorAxisCrossingDetector extends AbstractDetector<SemiMajorAxisCrossingDetector> {

        private final double threshold;

        public SemiMajorAxisCrossingDetector(double threshold) {
            super(AbstractDetector.DEFAULT_MAXCHECK,
                  AbstractDetector.DEFAULT_THRESHOLD,
                  AbstractDetector.DEFAULT_MAX_ITER,
                  new ContinueOnEvent());
            this.threshold = threshold;
        }

        public SemiMajorAxisCrossingDetector(double maxCheck,
                                             double threshold,
                                             int maxIter,
                                             EventHandler handler,
                                             double smaThreshold) {
            super(maxCheck, threshold, maxIter, handler);
            this.threshold = smaThreshold;
        }

        @Override
        protected SemiMajorAxisCrossingDetector create(AdaptableInterval newMaxCheck,
                                                       double newThreshold,
                                                       int newMaxIter,
                                                       EventHandler newHandler) {
            return new SemiMajorAxisCrossingDetector(newMaxCheck.currentInterval(null),
                                                     newThreshold, newMaxIter, newHandler, threshold);
        }

        @Override
        public double g(SpacecraftState s) {
            return s.getA() - threshold;
        }
    }

    private static class TriggeredEvent extends AbstractDetector<TriggeredEvent> {

        private final EventDetector triggeredDetector;
        private final EventDetector startTriggerDetector;
        private final EventDetector stopTriggerDetector;

        private boolean isTriggered = false;

        public TriggeredEvent(EventDetector triggeredDetector,
                              EventDetector startTriggerDetector,
                              EventDetector stopTriggerDetector) {
            super(AbstractDetector.DEFAULT_MAXCHECK,
                  AbstractDetector.DEFAULT_THRESHOLD,
                  AbstractDetector.DEFAULT_MAX_ITER,
                  new ContinueOnEvent());
            this.triggeredDetector    = triggeredDetector;
            this.startTriggerDetector = startTriggerDetector;
            this.stopTriggerDetector  = stopTriggerDetector;
        }

        public TriggeredEvent(double maxCheck,
                              double threshold,
                              int maxIter,
                              EventHandler handler,
                              EventDetector triggeredDetector,
                              EventDetector startTriggerDetector,
                              EventDetector stopTriggerDetector) {
            super(maxCheck, threshold, maxIter, handler);
            this.triggeredDetector    = triggeredDetector;
            this.startTriggerDetector = startTriggerDetector;
            this.stopTriggerDetector  = stopTriggerDetector;
        }

        @Override
        protected TriggeredEvent create(AdaptableInterval newMaxCheck,
                                        double newThreshold,
                                        int newMaxIter,
                                        EventHandler newHandler) {
            return new TriggeredEvent(newMaxCheck.currentInterval(null),
                                      newThreshold,
                                      newMaxIter,
                                      newHandler,
                                      triggeredDetector,
                                      startTriggerDetector,
                                      stopTriggerDetector);
        }

        @Override
        public double g(SpacecraftState s) {

            if (startTriggerDetector.g(s) > 0) {
                isTriggered = true;
            } else if (stopTriggerDetector.g(s) > 0) {
                isTriggered = false;
            }

            if (isTriggered) {
                return triggeredDetector.g(s);
            }

            return -1;
        }
    }

    private static class StationKeepingTrigger extends StartStopEventsTrigger<TriggeredEvent, BooleanDetector> {

        public StationKeepingTrigger(TriggeredEvent prototypeStartDetector,
                                     BooleanDetector prototypeStopDetector) {
            super(prototypeStartDetector, prototypeStopDetector);
        }

        @Override
        public List<ParameterDriver> getParametersDrivers() {
            return Collections.emptyList();
        }

        @Override
        protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>> FieldAbstractDetector<D, S> convertStartDetector(
                Field<S> field,
                TriggeredEvent detector) {
            return null;
        }

        @Override
        protected <D extends FieldAbstractDetector<D, S>, S extends CalculusFieldElement<S>> FieldAbstractDetector<D, S> convertStopDetector(
                Field<S> field,
                BooleanDetector detector) {
            return null;
        }
    }
}

Example plot