1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.events;
18  
19  import org.hipparchus.analysis.UnivariateFunction;
20  import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
21  import org.hipparchus.analysis.solvers.BracketedUnivariateSolver.Interval;
22  import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
23  import org.hipparchus.exception.MathRuntimeException;
24  import org.hipparchus.ode.events.Action;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.Precision;
27  import org.orekit.errors.OrekitException;
28  import org.orekit.errors.OrekitInternalError;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.propagation.SpacecraftState;
31  import org.orekit.propagation.sampling.OrekitStepInterpolator;
32  import org.orekit.time.AbsoluteDate;
33  
34  /** This class handles the state for one {@link EventDetector
35   * event detector} during integration steps.
36   *
37   * <p>This class is heavily based on the class with the same name from the
38   * Hipparchus library. The changes performed consist in replacing
39   * raw types (double and double arrays) with space dynamics types
40   * ({@link AbsoluteDate}, {@link SpacecraftState}).</p>
41   * <p>Each time the propagator proposes a step, the event detector
42   * should be checked. This class handles the state of one detector
43   * during one propagation step, with references to the state at the
44   * end of the preceding step. This information is used to determine if
45   * the detector should trigger an event or not during the proposed
46   * step (and hence the step should be reduced to ensure the event
47   * occurs at a bound rather than inside the step).</p>
48   * @author Luc Maisonobe
49   * @param <T> class type for the generic version
50   */
51  public class EventState<T extends EventDetector> {
52  
53      /** Event detector. */
54      private T detector;
55  
56      /** Time of the previous call to g. */
57      private AbsoluteDate lastT;
58  
59      /** Value from the previous call to g. */
60      private double lastG;
61  
62      /** Time at the beginning of the step. */
63      private AbsoluteDate t0;
64  
65      /** Value of the event detector at the beginning of the step. */
66      private double g0;
67  
68      /** Simulated sign of g0 (we cheat when crossing events). */
69      private boolean g0Positive;
70  
71      /** Indicator of event expected during the step. */
72      private boolean pendingEvent;
73  
74      /** Occurrence time of the pending event. */
75      private AbsoluteDate pendingEventTime;
76  
77      /**
78       * Time to stop propagation if the event is a stop event. Used to enable stopping at
79       * an event and then restarting after that event.
80       */
81      private AbsoluteDate stopTime;
82  
83      /** Time after the current event. */
84      private AbsoluteDate afterEvent;
85  
86      /** Value of the g function after the current event. */
87      private double afterG;
88  
89      /** The earliest time considered for events. */
90      private AbsoluteDate earliestTimeConsidered;
91  
92      /** Integration direction. */
93      private boolean forward;
94  
95      /** Variation direction around pending event.
96       *  (this is considered with respect to the integration direction)
97       */
98      private boolean increasing;
99  
100     /** Simple constructor.
101      * @param detector monitored event detector
102      */
103     public EventState(final T detector) {
104         this.detector     = detector;
105 
106         // some dummy values ...
107         lastT                  = AbsoluteDate.PAST_INFINITY;
108         lastG                  = Double.NaN;
109         t0                     = null;
110         g0                     = Double.NaN;
111         g0Positive             = true;
112         pendingEvent           = false;
113         pendingEventTime       = null;
114         stopTime               = null;
115         increasing             = true;
116         earliestTimeConsidered = null;
117         afterEvent             = null;
118         afterG                 = Double.NaN;
119 
120     }
121 
122     /** Get the underlying event detector.
123      * @return underlying event detector
124      */
125     public T getEventDetector() {
126         return detector;
127     }
128 
129     /** Initialize event handler at the start of a propagation.
130      * <p>
131      * This method is called once at the start of the propagation. It
132      * may be used by the event handler to initialize some internal data
133      * if needed.
134      * </p>
135      * @param s0 initial state
136      * @param t target time for the integration
137      *
138      */
139     public void init(final SpacecraftState s0,
140                      final AbsoluteDate t) {
141         detector.init(s0, t);
142         lastT = AbsoluteDate.PAST_INFINITY;
143         lastG = Double.NaN;
144     }
145 
146     /** Compute the value of the switching function.
147      * This function must be continuous (at least in its roots neighborhood),
148      * as the integrator will need to find its roots to locate the events.
149      * @param s the current state information: date, kinematics, attitude
150      * @return value of the switching function
151      */
152     private double g(final SpacecraftState s) {
153         if (!s.getDate().equals(lastT)) {
154             lastG = detector.g(s);
155             lastT = s.getDate();
156         }
157         return lastG;
158     }
159 
160     /** Reinitialize the beginning of the step.
161      * @param interpolator interpolator valid for the current step
162      */
163     public void reinitializeBegin(final OrekitStepInterpolator interpolator) {
164         forward = interpolator.isForward();
165         final SpacecraftState s0 = interpolator.getPreviousState();
166         this.t0 = s0.getDate();
167         g0 = g(s0);
168         while (g0 == 0) {
169             // extremely rare case: there is a zero EXACTLY at interval start
170             // we will use the sign slightly after step beginning to force ignoring this zero
171             // try moving forward by half a convergence interval
172             final double dt = (forward ? 0.5 : -0.5) * detector.getThreshold();
173             AbsoluteDate startDate = t0.shiftedBy(dt);
174             // if convergence is too small move an ulp
175             if (t0.equals(startDate)) {
176                 startDate = nextAfter(startDate);
177             }
178             t0 = startDate;
179             g0 = g(interpolator.getInterpolatedState(t0));
180         }
181         g0Positive = g0 > 0;
182         // "last" event was increasing
183         increasing = g0Positive;
184     }
185 
186     /** Evaluate the impact of the proposed step on the event detector.
187      * @param interpolator step interpolator for the proposed step
188      * @return true if the event detector triggers an event before
189      * the end of the proposed step (this implies the step should be
190      * rejected)
191      * @exception MathRuntimeException if an event cannot be located
192      */
193     public boolean evaluateStep(final OrekitStepInterpolator interpolator)
194         throws MathRuntimeException {
195 
196         forward = interpolator.isForward();
197         final SpacecraftState s1 = interpolator.getCurrentState();
198         final AbsoluteDate t1 = s1.getDate();
199         final double dt = t1.durationFrom(t0);
200         if (FastMath.abs(dt) < detector.getThreshold()) {
201             // we cannot do anything on such a small step, don't trigger any events
202             return false;
203         }
204         // number of points to check in the current step
205         final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / detector.getMaxCheckInterval()));
206         final double h = dt / n;
207 
208 
209         AbsoluteDate ta = t0;
210         double ga = g0;
211         for (int i = 0; i < n; ++i) {
212 
213             // evaluate handler value at the end of the substep
214             final AbsoluteDate tb = (i == n - 1) ? t1 : t0.shiftedBy((i + 1) * h);
215             final double gb = g(interpolator.getInterpolatedState(tb));
216 
217             // check events occurrence
218             if (gb == 0.0 || (g0Positive ^ (gb > 0))) {
219                 // there is a sign change: an event is expected during this step
220                 if (findRoot(interpolator, ta, ga, tb, gb)) {
221                     return true;
222                 }
223             } else {
224                 // no sign change: there is no event for now
225                 ta = tb;
226                 ga = gb;
227             }
228 
229         }
230 
231         // no event during the whole step
232         pendingEvent     = false;
233         pendingEventTime = null;
234         return false;
235 
236     }
237 
238     /**
239      * Find a root in a bracketing interval.
240      *
241      * <p> When calling this method one of the following must be true. Either ga == 0, gb
242      * == 0, (ga < 0  and gb > 0), or (ga > 0 and gb < 0).
243      *
244      * @param interpolator that covers the interval.
245      * @param ta           earliest possible time for root.
246      * @param ga           g(ta).
247      * @param tb           latest possible time for root.
248      * @param gb           g(tb).
249      * @return if a zero crossing was found.
250      */
251     private boolean findRoot(final OrekitStepInterpolator interpolator,
252                              final AbsoluteDate ta, final double ga,
253                              final AbsoluteDate tb, final double gb) {
254         // check there appears to be a root in [ta, tb]
255         check(ga == 0.0 || gb == 0.0 || ga > 0.0 && gb < 0.0 || ga < 0.0 && gb > 0.0);
256 
257         final double convergence = detector.getThreshold();
258         final int maxIterationCount = detector.getMaxIterationCount();
259         final BracketedUnivariateSolver<UnivariateFunction> solver =
260                 new BracketingNthOrderBrentSolver(0, convergence, 0, 5);
261 
262         // event time, just at or before the actual root.
263         AbsoluteDate beforeRootT = null;
264         double beforeRootG = Double.NaN;
265         // time on the other side of the root.
266         // Initialized the the loop below executes once.
267         AbsoluteDate afterRootT = ta;
268         double afterRootG = 0.0;
269 
270         // check for some conditions that the root finders don't like
271         // these conditions cannot not happen in the loop below
272         // the ga == 0.0 case is handled by the loop below
273         if (ta.equals(tb)) {
274             // both non-zero but times are the same. Probably due to reset state
275             beforeRootT = ta;
276             beforeRootG = ga;
277             afterRootT = shiftedBy(beforeRootT, convergence);
278             afterRootG = g(interpolator.getInterpolatedState(afterRootT));
279         } else if (ga != 0.0 && gb == 0.0) {
280             // hard: ga != 0.0 and gb == 0.0
281             // look past gb by up to convergence to find next sign
282             // throw an exception if g(t) = 0.0 in [tb, tb + convergence]
283             beforeRootT = tb;
284             beforeRootG = gb;
285             afterRootT = shiftedBy(beforeRootT, convergence);
286             afterRootG = g(interpolator.getInterpolatedState(afterRootT));
287         } else if (ga != 0.0) {
288             final double newGa = g(interpolator.getInterpolatedState(ta));
289             if (ga > 0 != newGa > 0) {
290                 // both non-zero, step sign change at ta, possibly due to reset state
291                 beforeRootT = ta;
292                 beforeRootG = newGa;
293                 afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
294                 afterRootG = g(interpolator.getInterpolatedState(afterRootT));
295             }
296         }
297 
298         // loop to skip through "fake" roots, i.e. where g(t) = g'(t) = 0.0
299         // executed once if we didn't hit a special case above
300         AbsoluteDate loopT = ta;
301         double loopG = ga;
302         while ((afterRootG == 0.0 || afterRootG > 0.0 == g0Positive) &&
303                 strictlyAfter(afterRootT, tb)) {
304             if (loopG == 0.0) {
305                 // ga == 0.0 and gb may or may not be 0.0
306                 // handle the root at ta first
307                 beforeRootT = loopT;
308                 beforeRootG = loopG;
309                 afterRootT = minTime(shiftedBy(beforeRootT, convergence), tb);
310                 afterRootG = g(interpolator.getInterpolatedState(afterRootT));
311             } else {
312                 // both non-zero, the usual case, use a root finder.
313                 // time zero for evaluating the function f. Needs to be final
314                 final AbsoluteDate fT0 = loopT;
315                 final UnivariateFunction f = dt -> {
316                     return g(interpolator.getInterpolatedState(fT0.shiftedBy(dt)));
317                 };
318                 // tb as a double for use in f
319                 final double tbDouble = tb.durationFrom(fT0);
320                 if (forward) {
321                     try {
322                         final Interval interval =
323                                 solver.solveInterval(maxIterationCount, f, 0, tbDouble);
324                         beforeRootT = fT0.shiftedBy(interval.getLeftAbscissa());
325                         beforeRootG = interval.getLeftValue();
326                         afterRootT = fT0.shiftedBy(interval.getRightAbscissa());
327                         afterRootG = interval.getRightValue();
328                         // CHECKSTYLE: stop IllegalCatch check
329                     } catch (RuntimeException e) {
330                         // CHECKSTYLE: resume IllegalCatch check
331                         throw new OrekitException(e, OrekitMessages.FIND_ROOT,
332                                 detector, loopT, loopG, tb, gb, lastT, lastG);
333                     }
334                 } else {
335                     try {
336                         final Interval interval =
337                                 solver.solveInterval(maxIterationCount, f, tbDouble, 0);
338                         beforeRootT = fT0.shiftedBy(interval.getRightAbscissa());
339                         beforeRootG = interval.getRightValue();
340                         afterRootT = fT0.shiftedBy(interval.getLeftAbscissa());
341                         afterRootG = interval.getLeftValue();
342                         // CHECKSTYLE: stop IllegalCatch check
343                     } catch (RuntimeException e) {
344                         // CHECKSTYLE: resume IllegalCatch check
345                         throw new OrekitException(e, OrekitMessages.FIND_ROOT,
346                                 detector, tb, gb, loopT, loopG, lastT, lastG);
347                     }
348                 }
349             }
350             // tolerance is set to less than 1 ulp
351             // assume tolerance is 1 ulp
352             if (beforeRootT.equals(afterRootT)) {
353                 afterRootT = nextAfter(afterRootT);
354                 afterRootG = g(interpolator.getInterpolatedState(afterRootT));
355             }
356             // check loop is making some progress
357             check(forward && afterRootT.compareTo(beforeRootT) > 0 ||
358                   !forward && afterRootT.compareTo(beforeRootT) < 0);
359             // setup next iteration
360             loopT = afterRootT;
361             loopG = afterRootG;
362         }
363 
364         // figure out the result of root finding, and return accordingly
365         if (afterRootG == 0.0 || afterRootG > 0.0 == g0Positive) {
366             // loop gave up and didn't find any crossing within this step
367             return false;
368         } else {
369             // real crossing
370             check(beforeRootT != null && !Double.isNaN(beforeRootG));
371             // variation direction, with respect to the integration direction
372             increasing = !g0Positive;
373             pendingEventTime = beforeRootT;
374             stopTime = beforeRootG == 0.0 ? beforeRootT : afterRootT;
375             pendingEvent = true;
376             afterEvent = afterRootT;
377             afterG = afterRootG;
378 
379             // check increasing set correctly
380             check(afterG > 0 == increasing);
381             check(increasing == gb >= ga);
382 
383             return true;
384         }
385 
386     }
387 
388     /**
389      * Get the next number after the given number in the current propagation direction.
390      *
391      * @param t input time
392      * @return t +/- 1 ulp depending on the direction.
393      */
394     private AbsoluteDate nextAfter(final AbsoluteDate t) {
395         return t.shiftedBy(forward ? +Precision.EPSILON : -Precision.EPSILON);
396     }
397 
398     /** Get the occurrence time of the event triggered in the current
399      * step.
400      * @return occurrence time of the event triggered in the current
401      * step.
402      */
403     public AbsoluteDate getEventDate() {
404         return pendingEventTime;
405     }
406 
407     /**
408      * Try to accept the current history up to the given time.
409      *
410      * <p> It is not necessary to call this method before calling {@link
411      * #doEvent(SpacecraftState)} with the same state. It is necessary to call this
412      * method before you call {@link #doEvent(SpacecraftState)} on some other event
413      * detector.
414      *
415      * @param state        to try to accept.
416      * @param interpolator to use to find the new root, if any.
417      * @return if the event detector has an event it has not detected before that is on or
418      * before the same time as {@code state}. In other words {@code false} means continue
419      * on while {@code true} means stop and handle my event first.
420      */
421     public boolean tryAdvance(final SpacecraftState state,
422                               final OrekitStepInterpolator interpolator) {
423         final AbsoluteDate t = state.getDate();
424         // check this is only called before a pending event.
425         check(!pendingEvent || !strictlyAfter(pendingEventTime, t));
426 
427         final boolean meFirst;
428 
429         if (strictlyAfter(t, earliestTimeConsidered)) {
430             // just found an event and we know the next time we want to search again
431             meFirst = false;
432         } else {
433             // check g function to see if there is a new event
434             final double g = g(state);
435             final boolean positive = g > 0;
436 
437             if (positive == g0Positive) {
438                 // g function has expected sign
439                 g0 = g; // g0Positive is the same
440                 meFirst = false;
441             } else {
442                 // found a root we didn't expect -> find precise location
443                 final AbsoluteDate oldPendingEventTime = pendingEventTime;
444                 final boolean foundRoot = findRoot(interpolator, t0, g0, t, g);
445                 // make sure the new root is not the same as the old root, if one exists
446                 meFirst = foundRoot && !pendingEventTime.equals(oldPendingEventTime);
447             }
448         }
449 
450         if (!meFirst) {
451             // advance t0 to the current time so we can't find events that occur before t
452             t0 = t;
453         }
454 
455         return meFirst;
456     }
457 
458     /**
459      * Notify the user's listener of the event. The event occurs wholly within this method
460      * call including a call to {@link EventDetector#resetState(SpacecraftState)}
461      * if necessary.
462      *
463      * @param state the state at the time of the event. This must be at the same time as
464      *              the current value of {@link #getEventDate()}.
465      * @return the user's requested action and the new state if the action is {@link
466      * Action#RESET_STATE Action.RESET_STATE}.
467      * Otherwise the new state is {@code state}. The stop time indicates what time propagation
468      * should stop if the action is {@link Action#STOP Action.STOP}.
469      * This guarantees the integration will stop on or after the root, so that integration
470      * may be restarted safely.
471      */
472     public EventOccurrence doEvent(final SpacecraftState state) {
473         // check event is pending and is at the same time
474         check(pendingEvent);
475         check(state.getDate().equals(this.pendingEventTime));
476 
477         final Action action = detector.eventOccurred(state, increasing == forward);
478         final SpacecraftState newState;
479         if (action == Action.RESET_STATE) {
480             newState = detector.resetState(state);
481         } else {
482             newState = state;
483         }
484         // clear pending event
485         pendingEvent     = false;
486         pendingEventTime = null;
487         // setup for next search
488         earliestTimeConsidered = afterEvent;
489         t0 = afterEvent;
490         g0 = afterG;
491         g0Positive = increasing;
492         // check g0Positive set correctly
493         check(g0 == 0.0 || g0Positive == g0 > 0);
494         return new EventOccurrence(action, newState, stopTime);
495     }
496 
497     /**
498      * Shift a time value along the current integration direction: {@link #forward}.
499      *
500      * @param t     the time to shift.
501      * @param delta the amount to shift.
502      * @return t + delta if forward, else t - delta. If the result has to be rounded it
503      * will be rounded to be before the true value of t + delta.
504      */
505     private AbsoluteDate shiftedBy(final AbsoluteDate t, final double delta) {
506         if (forward) {
507             final AbsoluteDate ret = t.shiftedBy(delta);
508             if (ret.durationFrom(t) > delta) {
509                 return ret.shiftedBy(-Precision.EPSILON);
510             } else {
511                 return ret;
512             }
513         } else {
514             final AbsoluteDate ret = t.shiftedBy(-delta);
515             if (t.durationFrom(ret) > delta) {
516                 return ret.shiftedBy(+Precision.EPSILON);
517             } else {
518                 return ret;
519             }
520         }
521     }
522 
523     /**
524      * Get the time that happens first along the current propagation direction: {@link
525      * #forward}.
526      *
527      * @param a first time
528      * @param b second time
529      * @return min(a, b) if forward, else max (a, b)
530      */
531     private AbsoluteDate minTime(final AbsoluteDate a, final AbsoluteDate b) {
532         return (forward ^ (a.compareTo(b) > 0)) ? a : b;
533     }
534 
535     /**
536      * Check the ordering of two times.
537      *
538      * @param t1 the first time.
539      * @param t2 the second time.
540      * @return true if {@code t2} is strictly after {@code t1} in the propagation
541      * direction.
542      */
543     private boolean strictlyAfter(final AbsoluteDate t1, final AbsoluteDate t2) {
544         if (t1 == null || t2 == null) {
545             return false;
546         } else {
547             return forward ? t1.compareTo(t2) < 0 : t2.compareTo(t1) < 0;
548         }
549     }
550 
551     /**
552      * Same as keyword assert, but throw a {@link MathRuntimeException}.
553      *
554      * @param condition to check
555      * @throws MathRuntimeException if {@code condition} is false.
556      */
557     private void check(final boolean condition) throws MathRuntimeException {
558         if (!condition) {
559             throw new OrekitInternalError(null);
560         }
561     }
562 
563     /**
564      * Class to hold the data related to an event occurrence that is needed to decide how
565      * to modify integration.
566      */
567     public static class EventOccurrence {
568 
569         /** User requested action. */
570         private final Action action;
571         /** New state for a reset action. */
572         private final SpacecraftState newState;
573         /** The time to stop propagation if the action is a stop event. */
574         private final AbsoluteDate stopDate;
575 
576         /**
577          * Create a new occurrence of an event.
578          *
579          * @param action   the user requested action.
580          * @param newState for a reset event. Should be the current state unless the
581          *                 action is {@link Action#RESET_STATE}.
582          * @param stopDate to stop propagation if the action is {@link Action#STOP}. Used
583          *                 to move the stop time to just after the root.
584          */
585         EventOccurrence(final Action action,
586                         final SpacecraftState newState,
587                         final AbsoluteDate stopDate) {
588             this.action = action;
589             this.newState = newState;
590             this.stopDate = stopDate;
591         }
592 
593         /**
594          * Get the user requested action.
595          *
596          * @return the action.
597          */
598         public Action getAction() {
599             return action;
600         }
601 
602         /**
603          * Get the new state for a reset action.
604          *
605          * @return the new state.
606          */
607         public SpacecraftState getNewState() {
608             return newState;
609         }
610 
611         /**
612          * Get the new time for a stop action.
613          *
614          * @return when to stop propagation.
615          */
616         public AbsoluteDate getStopDate() {
617             return stopDate;
618         }
619 
620     }
621 
622 }