Dear Orekit team,
This is my first post on this forum, thanks already for your time reading this.
I am trying to implement a custom attitude provider for a number of different applications. I am currently using Orekit in Python via your Python wrapper.
While testing my implementation, I noticed an unexpected behaviour of a custom (Python)AttitudeProviderModifier. So I have two questions for you.
The first one is - can you please tell me if this is the correct implementation of a custom attitude provider modifier? It simply wraps a YawCompensation provider, returning the original attitude (I am using this just to test that everything is in order).
### ATTITUDE PROVIDER MODIFIER - CUSTOM IMPLEMENTATION
class TestAttProviderModifier(PythonAttitudeProviderModifier):
EME2000 = FramesFactory.getEME2000()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)
nadir_pointing_att_provider = NadirPointing(EME2000, earth)
yaw_comp_att_provider = YawCompensation(EME2000, nadir_pointing_att_provider)
def getAttitude(self, pv_coordinates_provider, date, frame):
match (pv_coordinates_provider):
case PVCoordinatesProvider():
original_attitude = TestAttProviderModifier.yaw_comp_att_provider.getAttitude(pv_coordinates_provider, date, frame)
return original_attitude
case FieldPVCoordinatesProvider():
original_attitude = TestAttProviderModifier.yaw_comp_att_provider.getAttitude(pv_coordinates_provider, date, frame)
return FieldAttitude(date, frame, original_attitude.getRotation(), original_attitude.getRotationAcceleration())
case _:
raise RuntimeError(f'Not supported type of PVCoordinatesProvider: {type(pv_coordinates_provider).__name__}')
def getUnderlyingAttitudeProvider(self):
return TestAttProviderModifier.yaw_comp_att_provider
Then, assuming the implementation above is at least correct (but feel free to let me know if there is a nicer way to do it), I have found the following unexpected behaviour, seemingly only appearing in one specific condition, that is, if I propagate a non-equatorial orbit with a low-thrust manoeuvre taking place at a certain point during the propagation. In these conditions, it looks like the attitude using my custom AttitudeProviderModifier diverges from the underlying attitude even though I am just returning the same underlying attitude. And since the thrust direction is, of course, attitude dependent, the orbit also tends to diverge given time, especially if further low thrust manoeuvres take place at a later stage.
I have managed to reproduce this behaviour with a short code, which I am attaching down here. Let me know if you can have a look and spot any mistake, or if this is something you expect, for any reason.
Here you can see the “rotation distance” between the two attitudes. This was generated using the code below, where you can find all the details, and it’s for a LEO with a 1h-long low thrust manoeuvre happening after 5 days of propagation. 20 days are propagated in total.
Thank you very much for your help!
Here is the code (remember to change the path to the orekit data)
import orekit
from orekit import JArray_double
from pathlib import Path
from org.orekit.attitudes import NadirPointing, YawCompensation, PythonAttitudeProviderModifier, FieldAttitude
from org.orekit.utils import PVCoordinatesProvider, FieldPVCoordinatesProvider
from org.orekit.bodies import OneAxisEllipsoid
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.util import FastMath
from org.orekit.frames import FramesFactory
from org.orekit.orbits import EquinoctialOrbit, PositionAngleType, OrbitType
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions
from org.orekit.forces import maneuvers
from org.orekit.forces.maneuvers import trigger
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.forces.maneuvers.propulsion import ThrustDirectionAndAttitudeProvider as thrust_dir_prov
from org.hipparchus.geometry.euclidean.threed import Rotation
import matplotlib.pyplot as plt
vm = orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir
dir_path = Path(__file__).parent
orekit_data = Path(dir_path, <PATH_TO_OREKIT_DATA>)
setup_orekit_curdir(str(orekit_data))
### CONSTANTS
UTC = TimeScalesFactory.getUTC()
INERTIAL = FramesFactory.getEME2000()
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
EARTH_BODYSHAPE = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ITRF)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
INITIAL_DATE = AbsoluteDate(2024, 8, 30, 15, 0, 00.000, UTC)
### ATTITUDE PROVIDER MODIFIER - CUSTOM IMPLEMENTATION
class TestAttProviderModifier(PythonAttitudeProviderModifier):
EME2000 = FramesFactory.getEME2000()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf)
nadir_pointing_att_provider = NadirPointing(EME2000, earth)
yaw_comp_att_provider = YawCompensation(EME2000, nadir_pointing_att_provider)
def getAttitude(self, pv_coordinates_provider, date, frame):
match (pv_coordinates_provider):
case PVCoordinatesProvider():
original_attitude = TestAttProviderModifier.yaw_comp_att_provider.getAttitude(pv_coordinates_provider, date, frame)
return original_attitude
case FieldPVCoordinatesProvider():
original_attitude = TestAttProviderModifier.yaw_comp_att_provider.getAttitude(pv_coordinates_provider, date, frame)
return FieldAttitude(date, frame, original_attitude.getRotation(), original_attitude.getRotationAcceleration())
case _:
raise RuntimeError(f'Not supported type of PVCoordinatesProvider: {type(pv_coordinates_provider).__name__}')
def getUnderlyingAttitudeProvider(self):
return TestAttProviderModifier.yaw_comp_att_provider
### HELPER METHODS
def get_numerical_propagator(initial_orbit, attitude_provider):
# Standard numerical propagator, nothing fancy here
min_step = 0.001
maxstep = 1000.0
init_step = 60.0
position_tolerance = 0.000001
tolerances = NumericalPropagator.tolerances(
position_tolerance,
initial_orbit,
OrbitType.EQUINOCTIAL
)
integrator = DormandPrince853Integrator(
min_step, maxstep,
JArray_double.cast_(tolerances[0]),
JArray_double.cast_(tolerances[1])
)
integrator.setInitialStepSize(init_step)
return NumericalPropagator(integrator, attitude_provider)
def get_orbit():
# LEO with non-null inclination
a = 500. * 1000. + EARTH_RADIUS
return EquinoctialOrbit(
a, 0.001, 0.001, 0.7, 0.7,
FastMath.toRadians(30.), PositionAngleType.TRUE, INERTIAL,
INITIAL_DATE, Constants.WGS84_EARTH_MU
)
def get_manoeuvre(man_start_epoch : AbsoluteDate):
# Get a 1h long low thrust manoeuvre without attitude override, aligned with body x direction, starting at a given date.
maneuver_trigger = trigger.DateBasedManeuverTriggers(man_start_epoch, 3600.)
thrust_direction = Vector3D(1., 0., 0.)
thrust_direction_provider = thrust_dir_prov.buildFromFixedDirectionInSatelliteFrame(thrust_direction)
thrust = 0.001
isp = 3000.
low_thrust_man = maneuvers.ConfigurableLowThrustManeuver(
thrust_direction_provider,
maneuver_trigger,
thrust,
isp
)
return low_thrust_man
### MAIN
if __name__ == '__main__':
# Get initial orbit
initial_orbit = get_orbit()
# Get two identical attitudes, but one is generated using TestAttProviderModifier(), see definition above.
# Expected behaviour: the two attitude should behave exactly in the same manner, since TestAttProviderModifier
# simply returns the attitude generated using YawCompensation. (right?)
# Yaw compensation
nadir_pointing_att_provider = NadirPointing(INERTIAL, EARTH_BODYSHAPE)
yaw_comp_att_provider = YawCompensation(INERTIAL, nadir_pointing_att_provider)
# Wrapped yaw compensation
CUSTOM_attitude_provider = TestAttProviderModifier()
# Initialize two numerical propagators using the two attitude providers just initialized (and the same orbit)
propagator_yaw_compensation = get_numerical_propagator(initial_orbit, yaw_comp_att_provider)
propagator_CUSTOM_attitude = get_numerical_propagator(initial_orbit, CUSTOM_attitude_provider)
# Create initial state
initial_state = SpacecraftState(initial_orbit)
# Get a low thrust manoeuvre (1h long, starting 5 days after propagation start)
man = get_manoeuvre(initial_orbit.getDate().shiftedBy(5* 24*3600.))
# Compute date list for propagation and plots
propagation_target_date = initial_orbit.getDate().shiftedBy(20* 24*3600.)
date_list = [initial_orbit.getDate().shiftedBy(h*100.) for h in range(36*24*20)] # for plotting
# Setup both propagators
propagator_yaw_compensation.setOrbitType(OrbitType.EQUINOCTIAL)
propagator_CUSTOM_attitude.setOrbitType( OrbitType.EQUINOCTIAL)
propagator_yaw_compensation.setInitialState(initial_state)
propagator_CUSTOM_attitude.setInitialState( initial_state)
propagator_yaw_compensation.addForceModel(man)
propagator_CUSTOM_attitude.addForceModel( man)
eph_generator_yaw_compensation = propagator_yaw_compensation.getEphemerisGenerator()
eph_generator_CUSTOM_attitude = propagator_CUSTOM_attitude.getEphemerisGenerator()
# Propagate and extract attitude rotations at all times
propagator_yaw_compensation.propagate(propagation_target_date)
bound_prop_yaw_compensation = eph_generator_yaw_compensation.getGeneratedEphemeris()
propagator_CUSTOM_attitude.propagate(propagation_target_date)
bound_prop_CUSTOM_attitude = eph_generator_CUSTOM_attitude.getGeneratedEphemeris()
rotations_yaw_compensation = [bound_prop_yaw_compensation.propagate(t).getAttitude().getRotation() for t in date_list]
rotations_CUSTOM_attitude = [bound_prop_CUSTOM_attitude.propagate(t).getAttitude().getRotation() for t in date_list]
rotation_distance = [FastMath.toDegrees(Rotation.distance(rotations_yaw_compensation[i], rotations_CUSTOM_attitude[i])) for i in range(len(date_list))]
# Plot rotation distances
time_vec = [i*100/(3600*24) for i in range(len(date_list))]
plt.plot(time_vec, rotation_distance)
.
.