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:
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.
!
.