First of all, I wanted to congratulate everyone that has contributed to the development of Orekit, it is an amazing tool! I am new to both Java and Orekit so I am advancing quite slowly, but I love what I have seen so far.
I have a question regarding the event detector. I am using the code below to compute accesses to London from a satellite flying at some orbit. The code is working, but I do not really understand how it really works. I have gone through Orekit’s documentation for the elevation detector but I still do not really understand it. Could someone please explain to me the function of s, detector, increasing, maxcheck, the actions STOP and CONTINUE…
Apologies beforehand if this is a silly question, as I said I am new to all of this and many times my lack of understanding simply derives from my lack of experience in coding.
Thank you very much!
Iliass
indent preformatted text by 4 spaces
EventDetector londonVisibility =
new ElevationDetector(access.getMaxCheck(), access.getThreshold(), londonFrame).
withConstantElevation(access.getMinimumElevation()).
withHandler(
(s, detector, increasing) -> {
System.out.println(
" Visibility on " +
detector.getTopocentricFrame().getName() +
(increasing ? " begins at " : " ends at ") +
s.getDate()
);
return increasing ? Action.CONTINUE : Action.STOP;
}
);
Orekit detectors are designed in the same way as Hipparchus event detectors for Ordinary Differential Equations, which in turn extended a principle established decades ago for ODE and known as G-stop functions. There purpose is to embed the notion of discrete events inside a continuous process.
The trick to achieve this feature is to define a g continuous function dedicated to the event we are interested in so that the sign of the function value changes when the event occurs. The g function is a function of time and state. In Hipparchus world, these are simply a double value (the independent variable in the ODE) and an array of doubles. In Orekit world, everything is packed in a SpacecraftState which contains current date, orbit, attitude and even derivatives and additional state parameters if the user added some custom equations and state parameters. The value of the g function by itself is not really meaningful, but it should be continuous at least around zero crossings so the root finding algorithms used internally work properly (they do work with discontinuous functions at zero crossings, but this degrades a lot performances as root finding falls back to dichotomy in this case).
During propagation/ODE integration, the propagator/integrator looks ahead to what will happen next, based on its local modelling of the state. At this stage, the next step is not yet accepted, but the propagator already knows how to compute any state in this possible next step. It then evaluates all the g functions that have been registered and checks if one or several detectors change sign during this possible next step. It locates the roots, sort them according to propagation direction (forward or backward) and when it finally identifies the first g function that will change sign at date tk, then it knows that it is allowed to accept this step up to time tk (i.e. consider the computation is valid and its model for the state from current time to tk is valid). It then calls the handler associated with the detector to notify it that the discrete event that was modeled using the g function is triggered at that time and with this state, and asking it what to do next.
This notification corresponds to the handler eventOccurred method to be called. It is called with several parameters:
s is the current SpacecraftState (date, orbit, attitude, derivatives, additional state parameters),
detector is the detector that was triggered,
increasing is a flag that shows is the event corresponds to values going from negative to positive values or the reverse if propagation was forward (i.e. if an event g function is negative for t<10 and positive for t>10, then the flag would be set to increasing=true or the event triggered at t=10 regardless of the propagation being forward or backward).
After having processed the discrete event, the handler eventOccurred method returns, and it should return an indicator telling the propagator what to do next. If the handler returns Action.CONTINUE, then propagation should continue, at least until next event or propagation target time. If the handler returns Action.STOP, then propagation should stop immediately, and everything already computed for dates after the event is silently dropped, as if it was never computed. There are some other special values that can be returned, if propagator should continue but reset some internal parts (reset derivatives or reset state).
As explained in my previous post, the propagator/ODE integrator checks the registered g functions (by the way they are sometime called switching functions) to see if they change sign during the possible next step.The check is performed by evaluating the g function value and compare it to the value we got from previous check and see if sign changed between the two checks. If the sign changed, then we know that there is at least one zero crossing in between, and then we run a root finding algorithm to locate this crossing accurately (i.e. down to the tolerance setting). We will not run the root finding if the sign is the same as the sign of the previous check.
Checking the g function only once during each step would miss some zero crossings if the next step is too long (in fact if an even number of roots occur during the step, we would miss them). This is important in numerical propagation if the numerical propagator uses an adaptive step size integrator with a large maxStep, and it is even more critical in analytical propagators which depending on the configuration may use only one step from start to end.
So in order to work around this problem, the maxCheck parameter is used. It means that the propagator/ODE integrator must check the g function value at least once every maxCheck seconds (and also at least once every step, even if steps are shorter than maxCheck). The g function of ElevationDetector is the elevation angle minus the current min elevation, so it is positive when spacecraft is above horizon, negative when spacecraft is below horizon, and zero at raising/setting time). If this detector uses a 60s maxCheck and the propagator has a current step size of 300s, it implies the propagator will check the g function value (i.e. the spacecraft elevation above current horizon) at least 5 times during the step. So as long as raising/setting events are separated by at least 60 seconds, we will not miss them as we will notice the sign change of the associated g function, start the root finding and locate the horizon crossing to millisecond accuracy if tolerance was set to 10⁻³ s). If the spacecraft trajectory is such that it raises barely above horizon and immediately sets back below it in less than one minute, we will probably miss the event pair (we may find it if the check date luckily occurs during the short above horizon period, but it is not likely).
For ElevationDetector, this default 60 seconds maxCheck corresponds to classical operational constraints: you generally don’t care about too short visibility intervals. They are too short to do anything useful and too close to the limits so are expected to have a noisy radio link. Not bothering about missing too close events pairs is therefore fine there. For other detectors, we may be concerned about missing a pair of events, so we have to set up a shorter maxCheck, but the price is high: we will evaluate the g function many times, even when it is in a stable domain far from roots.
There is no way we can guarantee for any continuous function that we will always find all crossings and never miss a pair. In fact, it is simple to design a g function that has arbitrarily close pairs of roots just to demonstrate that we are doomed. Setting maxCheck properly is a compromise.
It is important to understand the difference between maxCheck and tolerance. maxCheck is an often quite large value (order of magnitude of minutes or hours in space flight dynamics) that is used to bracket roots and start a root finding algorithm. tolerance is an often quite small value (order of magnitude of milliseconds or microseconds) that is used as a convergence criterion for the root finding algorithm.
Thank you so much Luc! I am so sorry I did not answer this before, I actually stayed away from orekit for the past few months to dedicate my time to other projects but I am back now and very motivated to learn more. Thank you very much for this very informative answer!