Need Help Understanding ApsideDetector

Hello,

I have been unsuccessful in implementing the ApsideDetector to detect perigee and apogee for a given ephemeris. I have an externally generated ephemeris that is decimated every 180 seconds. I am using OrekitEphemeris to “read” the external ephemeris into an Orekit object and using the Apside Detector to detect apogee and perigee events.

The problem I am facing is that the logged apogee and perigee events don’t seem to be correct. For example, when plotting the dot product of the position and velocity vectors (using numpy, not the switching function value), I see that a switch from negative to positive occurs between 00:15 and 00:21. I would expect that my apsis detector would detect a perigee crossing somewhere within that time interval. However, I get a time value of 00:22, where the detector completely “misses” the switching interval. The errors in time are seemingly random in magnitude, sometimes seconds and other times minutes.

I have tried changing the MaxCheck and Threshold values on the detector but I can’t seem to find a sweet spot where the detector correctly identifies an apogee or perigee. Could someone help me understand why this behavior could be occurring? Or if I am missing something obvious?

I hope that I’ve described my problem adequately, but I’m happy to provide any code snippets to give additional context.

Thank you for any help you could provide!

Hello @oSharifi and welcome to the Orekit forum !

You already seem familiar with how the EventDetector works. The default settings should be enough as apogee and perigee events do not occur that often.

Hence i think that the problem probably lies in the implementation ? Could you provide a code snippet to show how you are trying to do this ?

Cheers,
Vincent

Hello Vincent! Thank you for your message!

Here is the code I’ve written up. I’m using the Python wrapper and followed this forum post somewhat.

I’m setting up my ephemeris like so:

UTC = TimeScalesFactory.getUTC()
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, False)  # get reference frame
GCRF = FramesFactory.getGCRF()
mu = Constants.WGS84_EARTH_MU

# Make OrekitEphemeris by reading values from provided ephemeris
orekit_sc_states = ArrayList()
init_t = eph.data[0].t
init_state = eph.data[0]
init_pv = PVCoordinates(Vector3D(init_state.x, init_state.y, init_state.z), Vector3D(init_state.xdot,init_state.ydot,init_state.zdot))
initial_time = AbsoluteDate(init_t.year, init_t.month, init_t.day, init_t.hour, init_t.minute, init_t.second+init_t.microsecond/1e6, UTC)
initial_orbit = CartesianOrbit(init_pv, GCRF, initial_time , mu)
r_dot_v = []
t = []
for state in eph.data:
    t_gcrf = state.t
    date_gcrf = AbsoluteDate(t_gcrf.year, t_gcrf.month, t_gcrf.day, t_gcrf.hour, t_gcrf.minute, t_gcrf.second+t_gcrf.microsecond/1e6, UTC)
    pos_gcrf = Vector3D(state.x, state.y, state.z)
    vel_gcrf = Vector3D(state.xdot, state.ydot, state.zdot)
    coords_GCRF = PVCoordinates(pos_gcrf, vel_gcrf)
    abs_coords_GCRF = AbsolutePVCoordinates(GCRF, date_gcrf, coords_GCRF)
    orekit_sc_states.add(SpacecraftState(abs_coords_GCRF)) 

    # Independently check r dot v for zero-crossings
    r = np.array([state.x, state.y, state.z])
    v = np.array([state.xdot, state.ydot, state.zdot])
    r_dot_v.append(np.dot(r,v))
    t.append(t_gcrf)

The events detector and logger setup are below. Here I was manually overriding the maxCheck and Threshold values from the initial CartesianOrbit to try and get the detector to give correct values

apsides_logger = EventsLogger()
orekit_ephemeris = OrekitEphemeris(orekit_sc_states,3)  
apsides_detector = ApsideDetector(initial_orbit).withHandler(RecordAndContinue()).withMaxCheck(180.0 * 2).withThreshold(1.0e-12)
orekit_ephemeris.addEventDetector(apsides_detector)
orekit_ephemeris.addEventDetector(apsides_logger.monitorDetector(apsides_detector))

Finally viewing and storing the logged events:

result = orekit_ephemeris.propagate(orekit_ephemeris.getMaxDate())
apsides_events = apsides_logger.getLoggedEvents()
apogee_epochs = []
perigee_epochs = []
for apsis_event in apsides_events:
    date = datetime.strptime(apsis_event.getState().getDate().toString()[:-9], '%Y-%m-%dT%H:%M:%S.%f')
    event_type = "Periapsis" if apsis_event.isIncreasing() else "Apoapsis"
    if event_type == "Periapsis":
        perigee_epochs.append(date.replace(tzinfo=timezone.utc))
    else:
        apogee_epochs.append(date.replace(tzinfo=timezone.utc))

Here is a plot that I’ve used in my analysis:
image
This is a plot of \vec{r} \cdot \vec{v} that I compute using numpy. Here, there is a zero crossing between 17:42 and 17:45, implying apoapsis sometime during this interval. However, the events logger gives me a time of 17:46 for this particular apoapsis event. This isn’t limited to just apoapsis events; the periapsis events have a similar error.

You are welcome :slight_smile: !

I have made a script based on your code with a fake Ephemeris and the ApsideDetector and r_dot_v values matches :thinking:.

Hence a few questions :slight_smile: !

  • Are you certain that the ephemeris data that you have is in GCRF ?
  • Similarly, are you sure that the dates are in UTC ?
  • Is there a constant time shift between recorded events in both methods ?

Hi Vincent,

  • I have confirmed that it is GCRF ephemeris
  • Also confirmed that dates are UTC
  • There is no pattern to the time offset between a detected perigee/apogee and the r_dot_v crossing.

With the default settings, the apogee/perigee events are so off from their expected values. The only way I’ve gotten somewhat close values is by changing the MaxCheck value to double the ephemeris time step (180s). Is your implementation of this detector completely different from mine?

Hello @oSharifi,

Implementation should be quite similar as i have started from your code snippets and added an ephemeris generation beforehand. You can check my code below.

Would it be possible for you to provide a sample of your ephemeris data to be closer to your case ?

Code sample

Please change the orekit_data = Path(dir_path, “PATH/TO/OREKIT/DATA”) accordingly.

from datetime import datetime, timezone
from pathlib import Path

import numpy as np
import orekit
from matplotlib import pyplot as plt
from orekit.pyhelpers import setup_orekit_curdir
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.frames import FramesFactory
from org.orekit.orbits import CartesianOrbit
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.analytical import Ephemeris, KeplerianPropagator
from org.orekit.propagation.events import EventsLogger, ApsideDetector
from org.orekit.propagation.events.handlers import RecordAndContinue
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import AbsolutePVCoordinates, PVCoordinates, Constants
from java.util import ArrayList
import plotly.express as px

vm = orekit.initVM()
dir_path = Path(__file__).parent
orekit_data = Path(dir_path, "PATH/TO/OREKIT/DATA")
setup_orekit_curdir(str(orekit_data))

UTC = TimeScalesFactory.getUTC()
GCRF = FramesFactory.getGCRF()
mu = Constants.WGS84_EARTH_MU

orekit_sc_states = ArrayList()
initial_time = AbsoluteDate()
initial_orbit = CartesianOrbit(PVCoordinates(Vector3D(6378e3 + 500e3, 0., 0.),
                                             Vector3D(0., 7613., 0.)),
                               GCRF, initial_time, mu)

# Generate fake ephemeris
kep_propagator = KeplerianPropagator(initial_orbit)
ephemeris = [kep_propagator.propagate(initial_time.shiftedBy(i * 60.)) for i in range(1440)]

# Make OrekitEphemeris by reading values from provided ephemeris
r_dot_v = []
t = []
for state in ephemeris:
    date_gcrf = state.getDate()
    coords_GCRF = state.getPVCoordinates()
    pos = coords_GCRF.getPosition()
    vel = coords_GCRF.getVelocity()
    abs_coords_GCRF = AbsolutePVCoordinates(GCRF, date_gcrf, coords_GCRF)
    orekit_sc_states.add(SpacecraftState(abs_coords_GCRF))

    # Independently check r dot v for zero-crossings
    r = np.array([pos.x, pos.y, pos.z])
    v = np.array([vel.x, vel.y, vel.z])
    r_dot_v.append(np.dot(r, v))
    t.append(date_gcrf.durationFrom(initial_time))

apsides_logger = EventsLogger()
orekit_ephemeris = Ephemeris(orekit_sc_states, 3)
apsides_detector = (ApsideDetector(initial_orbit).withHandler(RecordAndContinue())
                    .withMaxCheck(180.0 * 2)
                    .withThreshold(1.0e-12))

orekit_ephemeris.addEventDetector(apsides_detector)
orekit_ephemeris.addEventDetector(apsides_logger.monitorDetector(apsides_detector))

result = orekit_ephemeris.propagate(orekit_ephemeris.getMaxDate())
apsides_events = apsides_logger.getLoggedEvents()
apogee_epochs = []
perigee_epochs = []

for apsis_event in apsides_events:
    date = datetime.strptime(apsis_event.getState().getDate().toString()[:-9], '%Y-%m-%dT%H:%M:%S.%f')
    event_type = "Periapsis" if apsis_event.isIncreasing() else "Apoapsis"
    if event_type == "Periapsis":
        perigee_epochs.append(apsis_event.getDate().durationFrom(initial_time))
    else:
        apogee_epochs.append(apsis_event.getDate().durationFrom(initial_time))

print("Perigee epochs")
print(perigee_epochs)
print()
print("Apogee epochs")
print(apogee_epochs)

fig = px.scatter(x=t, y=r_dot_v)
fig.show()

Cheers,
Vincent