1   /* Copyright 2002-2020 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS 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.analytical;
18  
19  import java.util.ArrayList;
20  import java.util.Collection;
21  import java.util.Collections;
22  import java.util.Comparator;
23  import java.util.List;
24  import java.util.PriorityQueue;
25  import java.util.Queue;
26  
27  import org.hipparchus.exception.MathRuntimeException;
28  import org.hipparchus.ode.events.Action;
29  import org.hipparchus.util.FastMath;
30  import org.orekit.attitudes.Attitude;
31  import org.orekit.attitudes.AttitudeProvider;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.errors.OrekitInternalError;
34  import org.orekit.frames.Frame;
35  import org.orekit.orbits.Orbit;
36  import org.orekit.propagation.AbstractPropagator;
37  import org.orekit.propagation.AdditionalStateProvider;
38  import org.orekit.propagation.BoundedPropagator;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.events.EventDetector;
41  import org.orekit.propagation.events.EventState;
42  import org.orekit.propagation.events.EventState.EventOccurrence;
43  import org.orekit.propagation.sampling.OrekitStepInterpolator;
44  import org.orekit.time.AbsoluteDate;
45  import org.orekit.utils.PVCoordinatesProvider;
46  import org.orekit.utils.TimeStampedPVCoordinates;
47  
48  /** Common handling of {@link org.orekit.propagation.Propagator} methods for analytical propagators.
49   * <p>
50   * This abstract class allows to provide easily the full set of {@link
51   * org.orekit.propagation.Propagator Propagator} methods, including all propagation
52   * modes support and discrete events support for any simple propagation method. Only
53   * two methods must be implemented by derived classes: {@link #propagateOrbit(AbsoluteDate)}
54   * and {@link #getMass(AbsoluteDate)}. The first method should perform straightforward
55   * propagation starting from some internally stored initial state up to the specified target date.
56   * </p>
57   * @author Luc Maisonobe
58   */
59  public abstract class AbstractAnalyticalPropagator extends AbstractPropagator {
60  
61      /** Provider for attitude computation. */
62      private PVCoordinatesProvider pvProvider;
63  
64      /** Start date of last propagation. */
65      private AbsoluteDate lastPropagationStart;
66  
67      /** End date of last propagation. */
68      private AbsoluteDate lastPropagationEnd;
69  
70      /** Initialization indicator of events states. */
71      private boolean statesInitialized;
72  
73      /** Indicator for last step. */
74      private boolean isLastStep;
75  
76      /** Event steps. */
77      private final Collection<EventState<?>> eventsStates;
78  
79      /** Build a new instance.
80       * @param attitudeProvider provider for attitude computation
81       */
82      protected AbstractAnalyticalPropagator(final AttitudeProvider attitudeProvider) {
83          setAttitudeProvider(attitudeProvider);
84          pvProvider               = new LocalPVProvider();
85          lastPropagationStart     = AbsoluteDate.PAST_INFINITY;
86          lastPropagationEnd       = AbsoluteDate.FUTURE_INFINITY;
87          statesInitialized        = false;
88          eventsStates             = new ArrayList<EventState<?>>();
89      }
90  
91      /** {@inheritDoc} */
92      public BoundedPropagator getGeneratedEphemeris() {
93          return new BoundedPropagatorView(lastPropagationStart, lastPropagationEnd);
94      }
95  
96      /** {@inheritDoc} */
97      public <T extends EventDetector> void addEventDetector(final T detector) {
98          eventsStates.add(new EventState<T>(detector));
99      }
100 
101     /** {@inheritDoc} */
102     public Collection<EventDetector> getEventsDetectors() {
103         final List<EventDetector> list = new ArrayList<EventDetector>();
104         for (final EventState<?> state : eventsStates) {
105             list.add(state.getEventDetector());
106         }
107         return Collections.unmodifiableCollection(list);
108     }
109 
110     /** {@inheritDoc} */
111     public void clearEventsDetectors() {
112         eventsStates.clear();
113     }
114 
115     /** {@inheritDoc} */
116     public SpacecraftState propagate(final AbsoluteDateAbsoluteDate">AbsoluteDate start, final AbsoluteDate target) {
117         try {
118 
119             initializePropagation();
120 
121             lastPropagationStart = start;
122 
123             final double dt       = target.durationFrom(start);
124             final double epsilon  = FastMath.ulp(dt);
125             SpacecraftState state = updateAdditionalStates(basicPropagate(start));
126 
127             // evaluate step size
128             final double stepSize;
129             if (getMode() == MASTER_MODE) {
130                 if (Double.isNaN(getFixedStepSize())) {
131                     stepSize = FastMath.copySign(state.getKeplerianPeriod() / 100, dt);
132                 } else {
133                     stepSize = FastMath.copySign(getFixedStepSize(), dt);
134                 }
135             } else {
136                 stepSize = dt;
137             }
138 
139             // initialize event detectors
140             for (final EventState<?> es : eventsStates) {
141                 es.init(state, target);
142             }
143 
144             // initialize step handler
145             if (getStepHandler() != null) {
146                 getStepHandler().init(state, target);
147             }
148 
149             // iterate over the propagation range
150             statesInitialized = false;
151             isLastStep = false;
152             do {
153 
154                 // go ahead one step size
155                 final SpacecraftState previous = state;
156                 AbsoluteDate t = previous.getDate().shiftedBy(stepSize);
157                 if ((dt == 0) || ((dt > 0) ^ (t.compareTo(target) <= 0)) ||
158                         (FastMath.abs(target.durationFrom(t)) <= epsilon)) {
159                     // current step exceeds target
160                     // or is target to within double precision
161                     t = target;
162                 }
163                 final SpacecraftState current = updateAdditionalStates(basicPropagate(t));
164                 final OrekitStepInterpolator interpolator = new BasicStepInterpolator(dt >= 0, previous, current);
165 
166 
167                 // accept the step, trigger events and step handlers
168                 state = acceptStep(interpolator, target, epsilon);
169 
170                 // Update the potential changes in the spacecraft state due to the events
171                 // especially the potential attitude transition
172                 state = updateAdditionalStates(basicPropagate(state.getDate()));
173 
174             } while (!isLastStep);
175 
176             // return the last computed state
177             lastPropagationEnd = state.getDate();
178             setStartDate(state.getDate());
179             return state;
180 
181         } catch (MathRuntimeException mrte) {
182             throw OrekitException.unwrap(mrte);
183         }
184     }
185 
186     /** Accept a step, triggering events and step handlers.
187      * @param interpolator interpolator for the current step
188      * @param target final propagation time
189      * @param epsilon threshold for end date detection
190      * @return state at the end of the step
191           * @exception MathRuntimeException if an event cannot be located
192      */
193     protected SpacecraftState acceptStep(final OrekitStepInterpolator interpolator,
194                                          final AbsoluteDate target, final double epsilon)
195         throws MathRuntimeException {
196 
197         SpacecraftState       previous = interpolator.getPreviousState();
198         final SpacecraftState current  = interpolator.getCurrentState();
199         OrekitStepInterpolator restricted = interpolator;
200 
201 
202         // initialize the events states if needed
203         if (!statesInitialized) {
204 
205             if (!eventsStates.isEmpty()) {
206                 // initialize the events states
207                 for (final EventState<?> state : eventsStates) {
208                     state.reinitializeBegin(interpolator);
209                 }
210             }
211 
212             statesInitialized = true;
213 
214         }
215 
216         // search for next events that may occur during the step
217         final int orderingSign = interpolator.isForward() ? +1 : -1;
218         final Queue<EventState<?>> occurringEvents = new PriorityQueue<>(new Comparator<EventState<?>>() {
219             /** {@inheritDoc} */
220             @Override
221             public int compare(final EventState<?> es0, final EventState<?> es1) {
222                 return orderingSign * es0.getEventDate().compareTo(es1.getEventDate());
223             }
224         });
225 
226         boolean doneWithStep = false;
227         resetEvents:
228         do {
229 
230             // Evaluate all event detectors for events
231             occurringEvents.clear();
232             for (final EventState<?> state : eventsStates) {
233                 if (state.evaluateStep(interpolator)) {
234                     // the event occurs during the current step
235                     occurringEvents.add(state);
236                 }
237             }
238 
239             do {
240 
241                 eventLoop:
242                 while (!occurringEvents.isEmpty()) {
243 
244                     // handle the chronologically first event
245                     final EventState<?> currentEvent = occurringEvents.poll();
246 
247                     // get state at event time
248                     SpacecraftState eventState = restricted.getInterpolatedState(currentEvent.getEventDate());
249 
250                     // try to advance all event states to current time
251                     for (final EventState<?> state : eventsStates) {
252                         if (state != currentEvent && state.tryAdvance(eventState, interpolator)) {
253                             // we need to handle another event first
254                             // remove event we just updated to prevent heap corruption
255                             occurringEvents.remove(state);
256                             // add it back to update its position in the heap
257                             occurringEvents.add(state);
258                             // re-queue the event we were processing
259                             occurringEvents.add(currentEvent);
260                             continue eventLoop;
261                         }
262                     }
263                     // all event detectors agree we can advance to the current event time
264 
265                     final EventOccurrence occurrence = currentEvent.doEvent(eventState);
266                     final Action action = occurrence.getAction();
267                     isLastStep = action == Action.STOP;
268 
269                     if (isLastStep) {
270                         // ensure the event is after the root if it is returned STOP
271                         // this lets the user integrate to a STOP event and then restart
272                         // integration from the same time.
273                         eventState = interpolator.getInterpolatedState(occurrence.getStopDate());
274                         restricted = restricted.restrictStep(previous, eventState);
275                     }
276 
277                     // handle the first part of the step, up to the event
278                     if (getStepHandler() != null) {
279                         getStepHandler().handleStep(restricted, isLastStep);
280                     }
281 
282                     if (isLastStep) {
283                         // the event asked to stop integration
284                         return eventState;
285                     }
286 
287                     if (action == Action.RESET_DERIVATIVES || action == Action.RESET_STATE) {
288                         // some event handler has triggered changes that
289                         // invalidate the derivatives, we need to recompute them
290                         final SpacecraftState resetState = occurrence.getNewState();
291                         if (resetState != null) {
292                             resetIntermediateState(resetState, interpolator.isForward());
293                             return resetState;
294                         }
295                     }
296                     // at this point action == Action.CONTINUE or Action.RESET_EVENTS
297 
298                     // prepare handling of the remaining part of the step
299                     previous = eventState;
300                     restricted = new BasicStepInterpolator(restricted.isForward(), eventState, current);
301 
302                     if (action == Action.RESET_EVENTS) {
303                         continue resetEvents;
304                     }
305 
306                     // at this point action == Action.CONTINUE
307                     // check if the same event occurs again in the remaining part of the step
308                     if (currentEvent.evaluateStep(restricted)) {
309                         // the event occurs during the current step
310                         occurringEvents.add(currentEvent);
311                     }
312 
313                 }
314 
315                 // last part of the step, after the last event. Advance all detectors to
316                 // the end of the step. Should only detect a new event here if an event
317                 // modified the g function of another detector. Detecting such events here
318                 // is unreliable and RESET_EVENTS should be used instead. Might as well
319                 // re-check here because we have to loop through all the detectors anyway
320                 // and the alternative is to throw an exception.
321                 for (final EventState<?> state : eventsStates) {
322                     if (state.tryAdvance(current, interpolator)) {
323                         occurringEvents.add(state);
324                     }
325                 }
326 
327             } while (!occurringEvents.isEmpty());
328 
329             doneWithStep = true;
330         } while (!doneWithStep);
331 
332         isLastStep = target.equals(current.getDate());
333 
334         // handle the remaining part of the step, after all events if any
335         if (getStepHandler() != null) {
336             getStepHandler().handleStep(interpolator, isLastStep);
337         }
338 
339         return current;
340 
341     }
342 
343     /** Get the mass.
344      * @param date target date for the orbit
345      * @return mass mass
346      */
347     protected abstract double getMass(AbsoluteDate date);
348 
349     /** Get PV coordinates provider.
350      * @return PV coordinates provider
351      */
352     public PVCoordinatesProvider getPvProvider() {
353         return pvProvider;
354     }
355 
356     /** Reset an intermediate state.
357      * @param state new intermediate state to consider
358      * @param forward if true, the intermediate state is valid for
359      * propagations after itself
360      */
361     protected abstract void resetIntermediateState(SpacecraftState state, boolean forward);
362 
363     /** Extrapolate an orbit up to a specific target date.
364      * @param date target date for the orbit
365      * @return extrapolated parameters
366      */
367     protected abstract Orbit propagateOrbit(AbsoluteDate date);
368 
369     /** Propagate an orbit without any fancy features.
370      * <p>This method is similar in spirit to the {@link #propagate} method,
371      * except that it does <strong>not</strong> call any handler during
372      * propagation, nor any discrete events, not additional states. It always
373      * stop exactly at the specified date.</p>
374      * @param date target date for propagation
375      * @return state at specified date
376      */
377     protected SpacecraftState basicPropagate(final AbsoluteDate date) {
378         try {
379 
380             // evaluate orbit
381             final Orbit orbit = propagateOrbit(date);
382 
383             // evaluate attitude
384             final Attitude attitude =
385                 getAttitudeProvider().getAttitude(pvProvider, date, orbit.getFrame());
386 
387             // build raw state
388             return new SpacecraftState(orbit, attitude, getMass(date));
389 
390         } catch (OrekitException oe) {
391             throw new OrekitException(oe);
392         }
393     }
394 
395     /** Internal PVCoordinatesProvider for attitude computation. */
396     private class LocalPVProvider implements PVCoordinatesProvider {
397 
398         /** {@inheritDoc} */
399         public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
400             return propagateOrbit(date).getPVCoordinates(frame);
401         }
402 
403     }
404 
405     /** {@link BoundedPropagator} view of the instance. */
406     private class BoundedPropagatorView extends AbstractAnalyticalPropagator implements BoundedPropagator {
407 
408         /** Min date. */
409         private final AbsoluteDate minDate;
410 
411         /** Max date. */
412         private final AbsoluteDate maxDate;
413 
414         /** Simple constructor.
415          * @param startDate start date of the propagation
416          * @param endDate end date of the propagation
417          */
418         BoundedPropagatorView(final AbsoluteDateluteDate">AbsoluteDate startDate, final AbsoluteDate endDate) {
419             super(AbstractAnalyticalPropagator.this.getAttitudeProvider());
420             if (startDate.compareTo(endDate) <= 0) {
421                 minDate = startDate;
422                 maxDate = endDate;
423             } else {
424                 minDate = endDate;
425                 maxDate = startDate;
426             }
427 
428             try {
429                 // copy the same additional state providers as the original propagator
430                 for (AdditionalStateProvider provider : AbstractAnalyticalPropagator.this.getAdditionalStateProviders()) {
431                     addAdditionalStateProvider(provider);
432                 }
433             } catch (OrekitException oe) {
434                 // as the providers are already compatible with each other,
435                 // this should never happen
436                 throw new OrekitInternalError(null);
437             }
438 
439         }
440 
441         /** {@inheritDoc} */
442         public AbsoluteDate getMinDate() {
443             return minDate;
444         }
445 
446         /** {@inheritDoc} */
447         public AbsoluteDate getMaxDate() {
448             return maxDate;
449         }
450 
451         /** {@inheritDoc} */
452         protected Orbit propagateOrbit(final AbsoluteDate target) {
453             return AbstractAnalyticalPropagator.this.propagateOrbit(target);
454         }
455 
456         /** {@inheritDoc} */
457         public double getMass(final AbsoluteDate date) {
458             return AbstractAnalyticalPropagator.this.getMass(date);
459         }
460 
461         /** {@inheritDoc} */
462         public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
463             return propagate(date).getPVCoordinates(frame);
464         }
465 
466         /** {@inheritDoc} */
467         public void resetInitialState(final SpacecraftState state) {
468             AbstractAnalyticalPropagator.this.resetInitialState(state);
469         }
470 
471         /** {@inheritDoc} */
472         protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
473             AbstractAnalyticalPropagator.this.resetIntermediateState(state, forward);
474         }
475 
476         /** {@inheritDoc} */
477         public SpacecraftState getInitialState() {
478             return AbstractAnalyticalPropagator.this.getInitialState();
479         }
480 
481         /** {@inheritDoc} */
482         public Frame getFrame() {
483             return AbstractAnalyticalPropagator.this.getFrame();
484         }
485 
486     }
487 
488     /** Internal class for local propagation. */
489     private class BasicStepInterpolator implements OrekitStepInterpolator {
490 
491         /** Previous state. */
492         private final SpacecraftState previousState;
493 
494         /** Current state. */
495         private final SpacecraftState currentState;
496 
497         /** Forward propagation indicator. */
498         private final boolean forward;
499 
500         /** Simple constructor.
501          * @param isForward integration direction indicator
502          * @param previousState start of the step
503          * @param currentState end of the step
504          */
505         BasicStepInterpolator(final boolean isForward,
506                               final SpacecraftState previousState,
507                               final SpacecraftState currentState) {
508             this.forward         = isForward;
509             this.previousState   = previousState;
510             this.currentState    = currentState;
511         }
512 
513         /** {@inheritDoc} */
514         @Override
515         public SpacecraftState getPreviousState() {
516             return previousState;
517         }
518 
519         /** {@inheritDoc} */
520         @Override
521         public boolean isPreviousStateInterpolated() {
522             // no difference in analytical propagators
523             return false;
524         }
525 
526         /** {@inheritDoc} */
527         @Override
528         public SpacecraftState getCurrentState() {
529             return currentState;
530         }
531 
532         /** {@inheritDoc} */
533         @Override
534         public boolean isCurrentStateInterpolated() {
535             // no difference in analytical propagators
536             return false;
537         }
538 
539         /** {@inheritDoc} */
540         @Override
541         public SpacecraftState getInterpolatedState(final AbsoluteDate date) {
542 
543             // compute the basic spacecraft state
544             final SpacecraftState basicState = basicPropagate(date);
545 
546             // add the additional states
547             return updateAdditionalStates(basicState);
548 
549         }
550 
551         /** {@inheritDoc} */
552         @Override
553         public boolean isForward() {
554             return forward;
555         }
556 
557         /** {@inheritDoc} */
558         @Override
559         public BasicStepInterpolator restrictStep(final SpacecraftState newPreviousState,
560                                                   final SpacecraftState newCurrentState) {
561             return new BasicStepInterpolator(forward, newPreviousState, newCurrentState);
562         }
563 
564     }
565 
566 }