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 }