Hi Serrof
The sample code below generates an error on line 186 because Orekit 13 removes the registerSwitchEvents method.
To resolve this, the line should be commented out.
The code works correctly in Orekit 12.2.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Aug 12 23:43:29 2025
"""
from orekit import initVM, JArray_double
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
from org.orekit.attitudes import AlignedAndConstrained, AttitudesSequence, LofOffset, PredefinedTarget
from org.orekit.bodies import CelestialBodyFactory, OneAxisEllipsoid
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.frames import FramesFactory, LOFType
from org.orekit.orbits import KeplerianOrbit, PositionAngleType
from org.orekit.propagation import SpacecraftState, PropagationType, Propagator
from org.orekit.propagation.events import EclipseDetector
from org.orekit.propagation.events.handlers import PythonEventHandler
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.propagation.semianalytical.dsst import DSSTPropagator
from org.orekit.propagation.semianalytical.dsst.forces import DSSTTesseral, DSSTThirdBody, DSSTZonal
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import AngularDerivativesFilter, Constants, IERSConventions, PVCoordinatesProvider
from org.hipparchus.geometry import Vector
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.events import Action
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
import matplotlib.pyplot as plt
import os
import sys
from math import radians, degrees
import numpy as np
here = os.getcwd()
up = os.path.dirname(here)
cond1 = os.path.basename(here) == 'python'
cond2 = os.path.basename(up) == 'python'
if cond1:
data = os.path.join(here, 'orekit-data-master.zip')
if not here in sys.path:
sys.path.append(here)
elif cond2:
data = os.path.join(up, 'orekit-data-master.zip')
if not up in sys.path:
sys.path.append(up)
initVM()
setup_orekit_curdir(data)
class EclipseEventHandler(PythonEventHandler):
def __init__(self):
self.eclipse_in = []
self.eclipse_out = []
super().__init__()
def init(self, s, target, detector):
pass
def eventOccurred(self, s, detector, increasing):
if increasing: # sunlit
self.eclipse_out.append(s.getDate())
else: # eclipse
self.eclipse_in.append(s.getDate())
return Action.CONTINUE
def finish(self, final_state, detector):
pass
class DefaultOrbitHandler(PythonOrekitFixedStepHandler):
def __init__(self):
self.state = []
self.time = []
self.sun_vec = []
super().__init__()
def init(self, s0, t, step):
pass
def handleStep(self, s):
eci2body = s.getAttitude().getRotation()
r_sun_eci = PVCoordinatesProvider.cast_(SUN).getPVCoordinates(s.getDate(), ECI).getPosition()
e_r_rel_eci = Vector.cast_(r_sun_eci.subtract(s.getPosition())).normalize()
e_r_rel_body = eci2body.applyTo(e_r_rel_eci)
self.time.append(s.getDate())
self.state.append(s)
self.sun_vec.append(e_r_rel_body)
def finish(self, s):
pass
#%% Earth frames
ECI = FramesFactory.getEME2000()
ECEF = FramesFactory.getITRF(IERSConventions.IERS_2003, True)
ELLIPSOID = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ECEF)
#%% Perturbations for mean elements
forces = []
# Gravity
gravity = GravityFieldFactory.getUnnormalizedProvider(30, 0)
zonal_force = DSSTZonal(gravity)
tesseral_force = DSSTTesseral(ECEF, Constants.WGS84_EARTH_ANGULAR_VELOCITY, gravity)
[forces.append(gf) for gf in [zonal_force, tesseral_force]]
# Third bodies
SUN = CelestialBodyFactory.getSun()
MOON = CelestialBodyFactory.getMoon()
[forces.append(DSSTThirdBody(b, b.getGM())) for b in [SUN, MOON]]
#%% Initial states
# Target orbit
INITIAL_DATE = AbsoluteDate(2025, 1, 1, 10, 30, 0., TimeScalesFactory.getUTC())
SMA0 = 6879376.415393126
ECC0 = 0.0012533557095235504
INC0 = 1.6997385137705967
AOP0 = 1.5707963267948966
RAAN0 = 4.51186235457749
TA0 = -1.5707963267948966
INITIAL_ORBIT = KeplerianOrbit(SMA0, ECC0, INC0, AOP0, RAAN0, TA0, PositionAngleType.TRUE,
ECI, INITIAL_DATE, Constants.EIGEN5C_EARTH_MU)
# Attitude
ac_day_law = AlignedAndConstrained(Vector3D.MINUS_K, PredefinedTarget.SUN,
Vector3D.MINUS_J, PredefinedTarget.MOMENTUM,
SUN, ELLIPSOID)
ac_night_law = AlignedAndConstrained(Vector3D.PLUS_K, PredefinedTarget.EARTH,
Vector3D.MINUS_J, PredefinedTarget.MOMENTUM,
SUN, ELLIPSOID)
INITIAL_STATE = SpacecraftState(INITIAL_ORBIT)
#%% Propagation
# Time
OUT_STEP = 60.
DURATION = 86400.
# Event detectors and step handlers
eclipse_handler = EclipseEventHandler()
eclipse_detector = EclipseDetector(SUN, Constants.IAU_2015_NOMINAL_SOLAR_RADIUS, ELLIPSOID).withHandler(eclipse_handler)
orbit_handler = DefaultOrbitHandler()
# Attitudes sequence
attitudes_seq = AttitudesSequence()
attitudes_seq.addSwitchingCondition(ac_day_law, ac_night_law, eclipse_detector, False, True, 60., AngularDerivativesFilter.USE_RRA, None)
attitudes_seq.addSwitchingCondition(ac_night_law, ac_day_law, eclipse_detector, True, False, 60., AngularDerivativesFilter.USE_RRA, None)
# attitudes_seq = AttitudesSwitcher()
# attitudes_seq.addSwitchingCondition(ac_day_law, ac_night_law, eclipse_detector, False, True, None)
# attitudes_seq.addSwitchingCondition(ac_night_law, ac_day_law, eclipse_detector, True, False, None)
if eclipse_detector.g(INITIAL_STATE) >= 0:
# initial position is in daytime
attitudes_seq.resetActiveProvider(ac_day_law)
else:
# initial position is in nighttime
attitudes_seq.resetActiveProvider(ac_night_law)
# Propagator
tolerances = DSSTPropagator.tolerances(1e-10, INITIAL_ORBIT)
integrator = DormandPrince853Integrator(1e-6, 86400.,
JArray_double.cast_(tolerances[0]), JArray_double.cast_(tolerances[1]))
dsst_prop = DSSTPropagator(integrator)
[dsst_prop.addForceModel(f) for f in forces]
dsst_prop.addEventDetector(eclipse_detector)
Propagator.cast_(dsst_prop).setStepHandler(OUT_STEP, orbit_handler)
dsst_prop.setInitialState(INITIAL_STATE, PropagationType.MEAN)
dsst_prop.setAttitudeProvider(attitudes_seq)
attitudes_seq.registerSwitchEvents(dsst_prop)
dsst_prop.propagate(INITIAL_DATE, INITIAL_DATE.shiftedBy(DURATION))
#%% Results
sun_angle = [degrees(Vector3D.angle(Vector3D.MINUS_K, sv)) for sv in orbit_handler.sun_vec]
time = [absolutedate_to_datetime(t) for t in orbit_handler.time]
eclipse_in = [absolutedate_to_datetime(t) for t in eclipse_handler.eclipse_in]
eclipse_out = [absolutedate_to_datetime(t) for t in eclipse_handler.eclipse_out]
if eclipse_detector.g(INITIAL_STATE) < 0:
eclipse_in.insert(0, time[0])
if len(eclipse_in) > len(eclipse_out):
eclipse_out.append(time[-1])
plt.close('all')
_, axs = plt.subplots(1, sharex=True)
axs.plot(time, sun_angle)
axs.set_xlabel('Time')
axs.set_ylabel('Sun Angle (deg)')
[axs.axvspan(s, e, alpha=0.3, color='grey') for s, e in zip(eclipse_in, eclipse_out)]
axs.grid()
Please review the result plots below:
As shown, without registerSwitchEvents in Orekit 13.1, the attitude does not switch, resulting in a constant sun angle of zero.