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