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:)