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 :