Hello @bcazabonne !
Thanks for looking into my issue
Yes, please find attached the code (file upload didnt work for some reason)
it is standalone (we will need though to update path to orekit-data.zip) and only requires numpy. I’m using orekit 11.3
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 20 11:13:38 2023
@author: languille
"""
from pathlib import Path, WindowsPath
import orekit
from orekit import pyhelpers
pathToOreKitData = WindowsPath('c:/users/languille/team/team/resources')
orekit_data_name = "orekit-data.zip"
# Check if a java engine has been started
vm = orekit.getVMEnv()
if vm is None:
# Init orekit VM
vm = orekit.initVM()
else:
vm.attachCurrentThread()
# loads the orekit-data.zip
pyhelpers.setup_orekit_curdir(filename=str(pathToOreKitData / orekit_data_name))
import time
from org.orekit.attitudes import Attitude, InertialProvider, BodyCenterPointing, LofOffset, AttitudesSequence, AttitudeProvider, TargetPointing
from org.orekit.attitudes import CelestialBodyPointed, YawSteering, PythonAttitudeProvider, FieldAttitude, PythonAttitudeProviderModifier
from org.orekit.frames import FramesFactory, LOFType, ITRFVersion
from org.orekit.propagation.events import DateDetector, EventSlopeFilter, FilterType
from org.orekit.propagation.events.handlers import ContinueOnEvent
from org.hipparchus.geometry.euclidean.threed import Vector3D, FieldVector3D, FieldRotation, Rotation, RotationOrder, RotationConvention
from org.orekit.frames import FramesFactory, Frame
from org.orekit.bodies import CelestialBodyFactory, GeodeticPoint
from orekit.pyhelpers import absolutedate_to_datetime, datetime_to_absolutedate
from org.orekit.utils import AngularDerivativesFilter, AngularCoordinates, IERSConventions
from org.orekit.utils import PVCoordinates, TimeStampedPVCoordinates, FieldPVCoordinates, TimeStampedFieldAngularCoordinates
from org.orekit.utils import TimeStampedAngularCoordinates, PVCoordinatesProvider, FieldPVCoordinatesProvider
from org.orekit.time import AbsoluteDate, FieldAbsoluteDate
from org.orekit.models.earth import ReferenceEllipsoid
from org.hipparchus import CalculusFieldElement
import inspect
from orekit.pyhelpers import absolutedate_to_datetime, datetime_to_absolutedate
from org.orekit.orbits import OrbitType, PositionAngle, EquinoctialOrbit
import datetime
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator, ClassicalRungeKuttaIntegrator
from orekit import JArray_double
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.events import PythonAbstractDetector, PythonEventDetector
from org.hipparchus.ode.events import Action
from org.orekit.utils import PVCoordinates
from org.hipparchus.geometry.euclidean.threed import Vector3D
import numpy as np
class BlindingPositionDetector(PythonEventDetector):
''' this detector is used to detect when satellite arrives in a position where blinding
occurs on selected sensor'''
def __init__(self, blindingPvProvider, prop, exclAng,
pointingLaw):
self.blindingPvProvider = blindingPvProvider
self.prop = prop
self.exclAng = exclAng
self.pointingLaw = pointingLaw
super(BlindingPositionDetector, self).__init__()
def init(self, s, T):
pass
def getThreshold(self):
return float(PythonAbstractDetector.DEFAULT_THRESHOLD)
def getMaxCheckInterval(self):
return float(3600)
def getMaxIterationCount(self):
return PythonAbstractDetector.DEFAULT_MAX_ITER
def g(self, s):
''' this function computes the separation angle ie. margin
before being blinded
'''
absoluteDate = s.getDate()
frame = s.getFrame()
sunDirection_J2000 = PVCoordinates(self.prop.getPVCoordinates(absoluteDate, frame),
self.blindingPvProvider.getPVCoordinates(absoluteDate, frame))
LVLH_att = self.pointingLaw.getAttitude(self.prop, absoluteDate, frame).getRotation()
losDir_B = np.array([1,0,0])
sensorLos_J2000 = LVLH_att.applyInverseTo(
Vector3D(float(losDir_B[0]), float(losDir_B[1]), float(losDir_B[2])))
sepAng = Vector3D.angle(sensorLos_J2000, sunDirection_J2000.getPosition())
return float(sepAng - self.exclAng)
def eventOccurred(self, s, increasing):
return Action.CONTINUE
#%% first attitude, Earth Pointed
earthEllipse = ReferenceEllipsoid(6378e6,float(0.),
FramesFactory.getEME2000(),
398600441800000.0, 7.29211585530666e-05)
earthPointed = BodyCenterPointing(FramesFactory.getEME2000(),
earthEllipse)
#%% second attitude,
sat_YS = YawSteering(FramesFactory.getEME2000(),earthPointed,
PVCoordinatesProvider.cast_(CelestialBodyFactory.getSun()),
Vector3D.PLUS_I)
#%% create SC
startDate = datetime_to_absolutedate(datetime.datetime(2018, 1, 1, 0, 0, 0))
orbit = EquinoctialOrbit(
float(42e6),
float(0),#ex
float(0),#ey
float(0),#hx
float(0),#hy
float(0),#lv
PositionAngle.TRUE,
FramesFactory.getEME2000(),
startDate,
float(398600441800000.0)
)
minStep = float(0.001)
maxStep = float(5000)
initStep = float(60)
positionTolerance = float(1.)
tol = NumericalPropagator.tolerances(positionTolerance, orbit, OrbitType.EQUINOCTIAL)
# select integration scheme
integrator = DormandPrince853Integrator(minStep, maxStep,
JArray_double.cast_(tol[0]),
JArray_double.cast_(tol[1]))
integrator.setInitialStepSize(initStep)
prop = NumericalPropagator(integrator, earthPointed)
#%% create detectors
sun = PVCoordinatesProvider.cast_(CelestialBodyFactory.getSun())
exclAng = np.deg2rad(65)
det1 = BlindingPositionDetector(sun, prop, exclAng, earthPointed)
det2 = BlindingPositionDetector(sun, prop, exclAng, earthPointed)
blindingStart = EventSlopeFilter(det1, FilterType.TRIGGER_ONLY_DECREASING_EVENTS)
blindingEnd = EventSlopeFilter(det2, FilterType.TRIGGER_ONLY_INCREASING_EVENTS)
#%% register
attProviderList = [earthPointed]
switchTriggers = [None]
attProviderList.append(sat_YS)
switchTriggers.append(blindingStart)
attProviderList.append(earthPointed)
switchTriggers.append(blindingEnd)
seq = AttitudesSequence()
for idx, attProv_i in enumerate(attProviderList):
# skip first row as it is initial attitude
if idx == 0:
continue
# create trigger
trigger = switchTriggers[idx]
# determine previous attitude law
prevAttProvider = attProviderList[idx-1]
# create switch
seq.addSwitchingCondition(prevAttProvider, attProviderList[idx], trigger, True,
True, float(1000), AngularDerivativesFilter.USE_RR,
None)
# register all events
prop.setAttitudeProvider(seq)
seq.registerSwitchEvents(prop)
#%% run
initialState = SpacecraftState(orbit, float(5000))
prop.setInitialState(initialState)
prop.setOrbitType(orbit.getType())
curDate = startDate.shiftedBy(float(100000))
curState = prop.propagate(curDate)