Usage of IntersatDirectViewDetector with Hermite Interpolator

Dear Orekit community,

after jumping on the following topic : Computation between spacecraft and moving objects, I tried to compute such calculations in order to detect when a random moving object, like an aircraft, can be seen from a spacecraft, circulating in a keplerian orbit.

The idea I tried is the following :

  1. Create the primary object, orbit and propagator (I think this is good, defining a Keplerian Orbit and a Keplerian propagator)
  2. Create the aircraft trajectory, using 3 waypoints.
  3. Create an ephemeris from this trajectory in order to use it in the InterSatDirectViewDetector
  4. add this detector to the primary propagator, propagate and retrieve the detector logged events

The thing is, I’m not sure how to use the created Hermite interpolator in order to define the aircraft trajectory. The aircraft trajectory is here defined using 3 waypoints defined by PVcoordinates.

I have the following code that is running without erros (using the python wrapper), but which is not really well implemented regarding the Ephemeris and Interpolator. What I thought right was to use a set of interpolated point to create the ephemeris (not yet coded here). Is it the right way to create the Ephemeris with the Hermite Interpolator ? Please not that the aicraft waypoints defining the trajectory is quite random for the moment …


# Start date of the propagation
initial_date = AbsoluteDate(2025, 1, 1, 0, 0, 0.0, TimeScalesFactory.getUTC())

ra = 800 * 1000.0  # Apogee in meters
rp = 800 * 1000.0  # Perigee in meters
i = 60.0  # Inclination in deg
omega = float(0)  # Perigee argument in deg
raan = float(0)  # Right ascension of ascending node in deg
lv = float(0)  # True anomaly in deg

a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0  # SMA
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a  # Ecc


# Creation of the primary orbit
initial_orbit = KeplerianOrbit(a,e,math.radians(i),math.radians(omega),math.radians(raan),math.radians(lv),PositionAngleType.TRUE,FramesFactory.getEME2000(),
        initial_date,Constants.WGS84_EARTH_MU)

propagator_1 = KeplerianPropagator(initial_orbit)
Propagation_duration = 500.0       # propagation duration in seconds


# Definition of the aircraft trajectory
aircraft_states = ArrayList()
aircraft_pv_coordinates = ArrayList()
for k in range(3):
    dt = k * 250.0  
    # Define PV coordinates every 250 s
    date = initial_date.shiftedBy(dt) 
    position = Vector3D(7078000.0 , 0.0 + k * 800*100.0, 0.0 + k * 600*100.0)
    velocity = Vector3D(0.0, 3725.0, 6450.0)
    acceleration = Vector3D(-7.0, 0.0, 0.0)
    # Create the absolutePVcoordinates
    abs_pv = AbsolutePVCoordinates(FramesFactory.getEME2000(), date,position , velocity,acceleration)
    aircraft_pv_coordinates.add(TimeStampedPVCoordinates(date, position, velocity, acceleration))
    # Create the spacecraft (so here aircraft) state in order to use it in the ephemeris
    state = SpacecraftState(abs_pv)
    aircraft_states.add(state)


interpolator = TimeStampedPVCoordinatesHermiteInterpolator(3)

target_date = initial_date.shiftedBy(30.0)
interpolated_point = interpolator.interpolate(target_date, aircraft_pv_coordinates) # Where to use ? in the following ephemeris ?

ephemeris = Ephemeris(aircraft_states,3)
detector  = InterSatDirectViewDetector(earth, ephemeris)
# Adding the IntersatDirectView detector to the primary propagator
propagator_1.addEventDetector(detector)

# Log all events
logger = EventsLogger()
logged_detector = logger.monitorDetector(detector)
propagator_1.addEventDetector(logged_detector)


# Propagate the Primary orbit
final_state = propagator_1.propagate(initial_date, initial_date.shiftedBy(Propagation_duration))

get_visu_event = logger.getLoggedEvents()
print(get_visu_event)

Thank you for the help you could provide me with !

EM

Hi everyone,

I changed a bit my code to implement the Ephemeris with the following :


interpolator = TimeStampedPVCoordinatesHermiteInterpolator(3)

interpolated_aircraft_states = ArrayList()
for k in range(0,500):
    target_date = initial_date.shiftedBy(float(k))
    state_k = interpolator.interpolate(target_date, aircraft_pv_coordinates)
    abs_pv_k = AbsolutePVCoordinates(FramesFactory.getEME2000(),state_k.getDate() , state_k.getPosition() , state_k.getVelocity(),state_k.getAcceleration())
    state_k_new = SpacecraftState(abs_pv_k)
    interpolated_aircraft_states.add(state_k_new)


ephemeris_aircraft = Ephemeris(interpolated_aircraft_states,500)

detector  = InterSatDirectViewDetector(earth, ephemeris_aircraft)
logger = EventsLogger()
logged_detector = logger.monitorDetector(detector)
propagator_1.addEventDetector(logged_detector)

# Propagate the Primary orbit
final_state = propagator_1.propagate(initial_date, initial_date.shiftedBy(499.0))


get_visu_event = logger.getLoggedEvents()
print(get_visu_event)

basically, there are more points defining the Ephemeris. However, by doing, this, I got the following issue :

final_state = propagator_1.propagate(initial_date, initial_date.shiftedBy(499.0))
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
org.orekit.errors.OrekitException: l'intervalle n'encadre pas une racine : f(0E0) = �, f(1E0) = �
        at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
        at org.orekit.propagation.analytical.AbstractAnalyticalPropagator.propagate(AbstractAnalyticalPropagator.java:172)
Caused by: org.hipparchus.exception.MathIllegalArgumentException: l'intervalle n'encadre pas une racine : f(0E0) = �, f(1E0) = �
        at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.doSolveInterval(BracketingNthOrderBrentSolver.java:209)
        at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.doSolve(BracketingNthOrderBrentSolver.java:150)
        at org.hipparchus.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:201)
        at org.hipparchus.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:207)
        at org.orekit.bodies.OneAxisEllipsoid.lowestAltitudeIntermediate(OneAxisEllipsoid.java:878)
        at org.orekit.propagation.events.InterSatDirectViewDetector.g(InterSatDirectViewDetector.java:181)
        at org.orekit.propagation.events.EventsLogger$LoggingWrapper.g(EventsLogger.java:215)
        at org.orekit.propagation.events.EventState.g(EventState.java:161)
        at org.orekit.propagation.events.EventState.evaluateStep(EventState.java:224)
        at org.orekit.propagation.analytical.AbstractAnalyticalPropagator.acceptStep(AbstractAnalyticalPropagator.java:232)
        at org.orekit.propagation.analytical.AbstractAnalyticalPropagator.propagate(AbstractAnalyticalPropagator.java:153)

I believe this means that the g() function in my InterSatDirectViewDetector is likely the source of the error. I tried debugging by propagating the ephemeris with a step handler to retrieve the intermediate states and I obtained some NaN during the propagation which may be related to the previous error …"

{2025-01-01T00:00:00.000, P(1036623.1099034841, 3984882.414086635, -5851533.291402861), V(278.91972960035093, -83.56215743722947, 10.023579250634649), A(-0.006687025011755266, -0.019560547445742173, -0.029999935572681435)}
{2025-01-01T00:00:25.000, P(1.0464638857470095E7, 2.0621753741947186E8, 4.950896883781035E8), V(-1.6539840736021843E7, 2.4707326858521926E8, 8.76555405985677E8), A(5.395575524759574E8, 3.2726064121191144E9, 5.676559208873535E9)}
{2025-01-01T00:00:50.000, P(2.217231944377835E43, -7.228534801295975E43, 6.043112211554776E43), V(5.445376620710503E42, -1.7515578746906846E44, 3.489764483946072E44), A(3.583993447965185E44, 6.697203908108107E44, -1.3555236148075946E44)}
{2025-01-01T00:01:15.000, P(NaN, NaN, NaN), V(NaN, NaN, NaN), A(NaN, NaN, NaN)}
{2025-01-01T00:02:05.000, P(NaN, NaN, NaN), V(NaN, NaN, NaN), A(NaN, NaN, NaN)}
{2025-01-01T00:02:30.000, P(NaN, NaN, NaN), V(NaN, NaN, NaN), A(NaN, NaN, NaN)}
...

If you have any ideas about the source of this issue or suggestions on how to modify my ephemeris definition, I’d appreciate it !

Thanks,

EM

Hi there,

as you say, you’re giving pretty random waypoints. In particular, the velocity is not the derivative of the position vector. I’m not sure how the Hermite interpolator can cope with this. Could you try giving something more physical?

Cheers,
Romain.

Hello @Emilien_mrlt,

I noticed the following issue:

ephemeris_aircraft = Ephemeris(interpolated_aircraft_states,500)

You are telling the ephemeris to perform interpolation using 500 interpolation points which is huge. You are going to encounter what is called the Runge’s phenomenon because you are building a fitting polynomial of degree ~1000. So i would highly recommend to use a smaller number (say 10~20):

ephemeris_aircraft = Ephemeris(interpolated_aircraft_states,10)

Also, not giving it “realistic” points should not be an issue.

Cheers,
Vincent