Hi @emery,
After several tries, i have finally found a way to make it work :
import orekit
from java.util import Collections
from java.util.stream import Stream
from orekit import JArray_double
from orekit.pyhelpers import setup_orekit_curdir
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.attitudes import FrameAlignedProvider
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.maneuvers import Maneuver
from org.orekit.forces.maneuvers.propulsion import ConstantThrustDirectionProvider
from org.orekit.forces.maneuvers.propulsion import PythonPropulsionModel
from org.orekit.forces.maneuvers.trigger import DateBasedManeuverTriggers
from org.orekit.frames import FramesFactory
from org.orekit.orbits import KeplerianOrbit
from org.orekit.orbits import OrbitType
from org.orekit.orbits import PositionAngleType
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.time import AbsoluteDate
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions
vm = orekit.initVM()
setup_orekit_curdir("orekit-data.zip")
class PropulsionModel(PythonPropulsionModel):
def __init__(self, thrustDirectionProvider, isp):
super().__init__()
self.thrustDirectionProvider = thrustDirectionProvider
self.isp = isp
def init(self, spacecraftState, absoluteDate):
pass
def getParametersDrivers(self):
return Collections.emptyList()
def getEventDetectors(self):
return Stream.empty()
def getAcceleration(self, spacecraftState, attitude, doubleArray):
return Vector3D(1., 0., 0.)
def getDirection(self, s):
direction = Vector3D.MINUS_I
return direction
def getFlowRate(self, s):
thrust_mag, _ = self.thrustDirectionProvider.computeThrust(s)
flow_rate = - thrust_mag / (self.isp * 9.80665)
return float(flow_rate)
def getIsp(self, s):
return float(self.isp)
def getThrust(self, s):
thrust_mag, _ = self.thrustDirectionProvider.computeThrust(s)
return thrust_mag
def getMassDerivatives(self, spacecraftState, doubleArray):
return 0.
def getThrustVector(self, s):
thrust_mag, _ = self.thrustDirectionProvider.computeThrust(s)
# propulsion vector in spacecraft frame
thrust_vec_u = [-1, 0, 0]
thrust_vec = thrust_mag * thrust_vec_u
if thrust_mag == 0:
thrust_vec = Vector3D(-1e-20, 0.0, 0.0)
else:
thrust_vec = Vector3D(float(thrust_vec[0]), float(thrust_vec[1]), float(thrust_vec[2]))
return thrust_vec
mu = Constants.IERS2010_EARTH_MU
initialDate = AbsoluteDate()
inertialFrame = FramesFactory.getEME2000()
initialOrbit = KeplerianOrbit(6878000., 0.001, 1., 0., 0., 0., PositionAngleType.MEAN, inertialFrame, initialDate, mu)
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)
satellite_mass = 100.0 # The models need a spacecraft mass, unit kg.
initialState = SpacecraftState(initialOrbit, satellite_mass)
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(OrbitType.CARTESIAN)
propagator_num.setInitialState(initialState)
gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
propagator_num.addForceModel(
HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))
# Simple attitude provider
attitudeProvider = FrameAlignedProvider(inertialFrame)
# Maneuver trigger
propDate = initialDate.shiftedBy(500.)
propDuration = 100.
maneuverTrigger = DateBasedManeuverTriggers(propDate, propDuration)
# Custom propulsion model
customPropulsionModel = PropulsionModel(ConstantThrustDirectionProvider(Vector3D(1., 0., 0.)), 300.)
# Create custom maneuver
customManeuver = Maneuver(attitudeProvider, maneuverTrigger, customPropulsionModel)
# Add custom maneuver
propagator_num.addForceModel(customManeuver)
# Propagate
end_state = propagator_num.propagate(initialDate, initialDate.shiftedBy(3600.0 * 1))
It turns out that using the PythonThrustPropulsionModel was leading to a null pointer exception being thrown no matter what i did so i switched to PythonPropulsionModel
I hope it suits your needs !
Cheers,
Vincent