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

Hello @Vincent and @Serrof,
Sorry for the late answer I was out of the office.

Thank you for both of your messages! I tried reducing the number of points to 8 and it worked without any issues! I didn’t need to use the HermiteInterpolator—just the ephemeris, as shown below:

itrf     = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
eme2000  = FramesFactory.getEME2000()

for point in trajectory.keys():
    point_date = initial_date.shiftedBy(trajectory[point]["t"])

    # 1) PV in ITRF
    pv = PVCoordinates(
        Vector3D(trajectory[point]["X"],
                 trajectory[point]["Y"],
                 trajectory[point]["Z"]),
        Vector3D(float(trajectory[point]["Vx"]),
                 float(trajectory[point]["Vy"]),
                 float(trajectory[point]["Vz"])),
        Vector3D(float(trajectory[point]["Ax"]),
                 float(trajectory[point]["Ay"]),
                 float(trajectory[point]["Az"]))
    )

    # 2) Transform → EME2000
    transform = itrf.getTransformTo(eme2000, point_date)
    pv_eme    = transform.transformPVCoordinates(pv)

    # 3) Create the state list in order to create the ephemeris
    abs_pv = AbsolutePVCoordinates(eme2000, point_date, pv_eme)
    state  = SpacecraftState(abs_pv)
    aircraft_states.add(state)

ephemeris = Ephemeris(aircraft_states, len(traj))

I’m getting back the position with the following code :


def Handler_For_OEM():
    class HandlerForDate(PythonOrekitFixedStepHandler):
        def init(self, so, t, step):
            super().__init__()
            self.data = {
                "dates": [],
                "absdates": [],
                "dt": [],
                "Pos_X": [],
                "Pos_Y": [],
                "Pos_Z": [],
            }
            self.print_counter = 0

        def handleStep(self, state):
            date = state.getDate()
            pv_coordinates = state.getPVCoordinates(FramesFactory.getEME2000())
            x = pv_coordinates.getPosition().getX()
            y = pv_coordinates.getPosition().getY()
            z = pv_coordinates.getPosition().getZ()
            mjd = date.getMJD()
            mjd_day = int(mjd)  
            mjd_sec = (mjd - mjd_day) * 86400 
            self.data["absdates"].append(absolute_to_datetime(date))
            self.data["dates"].append(mjd_day)
            self.data["dt"].append(round(mjd_sec))
            self.data["Pos_X"].append(x / 1000.0)
            self.data["Pos_Y"].append(y / 1000.0)
            self.data["Pos_Z"].append(z / 1000.0)


            self.print_counter += 1

        def finish(self, finalState):
            pass

        def to_dataframe(self):
            return pd.DataFrame(self.data)

    return HandlerForDate()

This code allows propagation without any issues. However, I was wondering whether my frame transformation from ITRF to EME2000 is correct. This conversion is necessary to propagate the ephemeris.

I was also wondering whether a rotation-related correction might be missing here—specifically to account for Earth’s rotation—particularly when transforming velocity and acceleration in the PV list. This concern arises from the unusual ground track I’m seeing between the created points (only between, meaning there could be an issue with the velocity/acceleration transform) . I’m referring to the kind of correction discussed in Frames conversion .

Thanks for the help !

EM

Upping thread (no time right now :grimacing:)

Hi everyone,

I’m still working through this issue and wondering whether my Ephemeris creation is accurate, especially with the second parameter.

Let’s say I have 100 points defined by position-velocity (PV) coordinates in the EME2000 frame. The second point is 1000 seconds after the first, and the remaining points are spaced 1 second apart.

I’m using the Ephemeris class, which is no longer working with nb_points > 2.

ephemeris = Ephemeris(aircraft_states, nb_points)
  • For nb_points > 2, the trajectory becomes erratic, unless the spacecraft is secretly taking a detour past Mars between extrapolation steps !
  • For nb_points = 1, the trajectory appears to follow a straight line between interpolated points, except near the interpolated points themselves, where it’s curved. It’s better than the multi-point case, but still not ideal.
  • For nb_points = 2, the trajectory looks accurate overall, except for the first two points, which seem inconsistent or incorrect. The behaviour of the trajectory is like the above one in the following picture (Maybe it is due to the propagator trying to create a continous trajectory by using the velocities and the accelerations of the first two points):

I was expecting a trajectory like the second one.
My questions are the following :

Is nb_points related to the polynomial order in the Hermite interpolator associated with the ephemeris?
How should nb_points be used correctly?
Is it possible to set up a 2nd-order polynomial function between the interpolation points? Is it wrong only because of wrong velocities and accelerations ? (which is possible in my case).

Thanks

EM