Station keeping using low thrust

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 :slight_smile: .

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 :

6 Likes