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