1   /* Copyright 2002-2025 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.forces.drag;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  import java.util.stream.Stream;
22  
23  import org.hipparchus.CalculusFieldElement;
24  import org.hipparchus.Field;
25  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.util.MathArrays;
28  import org.orekit.annotation.DefaultDataContext;
29  import org.orekit.frames.Frame;
30  import org.orekit.models.earth.atmosphere.Atmosphere;
31  import org.orekit.propagation.FieldSpacecraftState;
32  import org.orekit.propagation.SpacecraftState;
33  import org.orekit.propagation.events.DateDetector;
34  import org.orekit.propagation.events.EventDetector;
35  import org.orekit.propagation.events.FieldDateDetector;
36  import org.orekit.propagation.events.FieldEventDetector;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.FieldAbsoluteDate;
39  import org.orekit.time.TimeScale;
40  import org.orekit.time.TimeScalesFactory;
41  import org.orekit.utils.ParameterDriver;
42  import org.orekit.utils.ParameterDriversProvider;
43  import org.orekit.utils.TimeSpanMap;
44  import org.orekit.utils.TimeSpanMap.Span;
45  import org.orekit.utils.TimeSpanMap.Transition;
46  
47  
48  /** Time span atmospheric drag force model.
49   *  <p>
50   *  This class is closely related to {@link org.orekit.forces.drag.DragForce DragForce} class.<br>
51   *  The difference is that it has a {@link TimeSpanMap} of {@link DragSensitive} objects as attribute
52   *  instead of a single {@link DragSensitive} object. <br>
53   *  The idea behind this model is to allow the user to design a drag force model that can see its physical parameters
54   *  (drag coefficient and lift ratio) change with time, at dates chosen by the user. <br>
55   *  </p>
56   *  <p>
57   *  This is a behavior that can be sought in operational orbit determination.<br>
58   *  Indeed the solar activity has a strong influence on the local atmospheric density, and thus on the drag force effect.<br>
59   *  Solar activity is a physical phenomenon that is difficult to model and predict. <br>
60   *  The errors induced by this incomplete modeling can be estimated through the drag coefficients.<br>
61   *  Being able to define and estimate drag coefficients depending on user-chosen dates in a piecewise fashion allows for
62   *  a better  modeling of solar activity uncertainties.
63   *  </p>
64   *  <p>
65   *  A typical operational use case is to have a daily solar activity with three-hourly magnetic indexes provided by an
66   *  international organization (NOAA for example).<br>
67   *  Given this input, a user can define a piecewise drag force model with daily or three-hourly drag coefficients.<br>
68   *  Each timed coefficient will absorb a part of the uncertainties in the solar activity and will allow for a more accurate
69   *  orbit determination
70   *  </p>
71   *  <b>Usage</b>:<ul>
72   *  <li><u>Construction</u>: constructor takes an atmospheric model and a DragSensitive model.<br>
73   *  This last model will be your initial DragSensitive model and it will be initially valid for the whole time line.<br>
74   *  The real validity of this first entry will be truncated as other DragSensitive models are added.
75   *  <li><u>Time spans</u>: DragSensitive models are added using methods {@link #addDragSensitiveValidAfter(DragSensitive, AbsoluteDate)}
76   *   or {@link #addDragSensitiveValidBefore(DragSensitive, AbsoluteDate)}.<br>
77   *   Recommendations are the same than the ones in {@link TimeSpanMap}, meaning: <ul>
78   *   <li>As an entry is added, it truncates the validity of the neighboring entries already present in the map;
79   *   <li><b>The transition dates should be entered only once</b>. Repeating a transition date will lead to unexpected result and is not supported;
80   *   <li>It is advised to order your DragSensitive models chronologically when adding them to avoid any confusion.
81   *   </ul>
82   *   <li><u>Naming the parameter drivers</u>: It is strongly advised to give a custom name to the {@link ParameterDriver}(s)
83   *   of each DragSensitive model that is added to the object. This will allow you keeping track of the evolution of your models.<br>
84   *   Different names are mandatory to differentiate the different drivers.<br>
85   *   If you do not specify a name, a default name will be chosen. Example for the drag coefficient:<ul>
86   *   <li>Initial DragSensitive model: the driver's default name is "{@link DragSensitive#DRAG_COEFFICIENT}";
87   *   <li>Using {@link #addDragSensitiveValidAfter(DragSensitive, AbsoluteDate)}: the driver's default name is
88   *   "{@link DragSensitive#DRAG_COEFFICIENT} + {@link #DATE_AFTER} + date.toString()"
89   *   <li>Using {@link #addDragSensitiveValidBefore(DragSensitive, AbsoluteDate)}: the driver's default name is
90   *   "{@link DragSensitive#DRAG_COEFFICIENT} + {@link #DATE_BEFORE} + date.toString()"
91   *   </ul>
92   *   </ul>
93   *  <b>Example following previous recommendations</b>:<ul>
94   *  <li>Given:
95   *  <ul>
96   *  <li><code>atmosphere</code>: an {@link Atmosphere atmospheric model};
97   *  <li><code>isotropicDrag0, 1 and 2</code>: three {@link org.orekit.forces.drag.IsotropicDrag IsotropicDrag} models;
98   *  <li><code>date</code>: an {@link AbsoluteDate}.
99   *  </ul>
100  *  <li>Name the drivers:<br>
101  *  <code>isotropicDrag0.getDragParametersDrivers()[0].setName = "Cd0";</code><br>
102  *  <code>isotropicDrag1.getDragParametersDrivers()[0].setName = "Cd1";</code><br>
103  *  <code>isotropicDrag2.getDragParametersDrivers()[0].setName = "Cd2";</code><br>
104  *  <li>Initialize the model: <br>
105  *  <code>TimeSpanDragForce force = new TimeSpanDragForce(atmosphere, isotropicDrag0);</code>
106  *  <li>Set the second and third model one Julian day apart each:<br>
107  *  <code>force.addDragSensitiveValidAfter(isotropicDrag1, date.shiftedBy(Constants.JULIAN_DAY));</code><br>
108  *  <code>force.addDragSensitiveValidAfter(isotropicDrag2, date.shiftedBy(2 * Constants.JULIAN_DAY));</code><br>
109  *  <li>With this, your model will have the following properties:
110  *  <ul>
111  *  <li>t in ]-∞, date + 1 day [ / Cd = Cd0
112  *  <li>t in [date + 1 day, date + 2days [ / Cd = Cd1
113  *  <li>t in [date + 2 days, +∞ [ / Cd = Cd2
114  *  </ul>
115  *  </ul>
116  *  <p>
117  *  <b>Warning</b>:<br> The TimeSpanDragForce model is versatile and you could end up with non-physical modeling.<br>
118  *  For example you could add 2 {@link org.orekit.forces.drag.IsotropicDrag IsotropicDrag} models with different areas,
119  *  or one {@link org.orekit.forces.drag.IsotropicDrag IsotropicDrag} model and then one
120  *  {@link org.orekit.forces.BoxAndSolarArraySpacecraft BoxAndSolarArraySpacecraft} model.<br>
121  *  It is up to you to ensure that your models are consistent with each other, Orekit will not perform any check for that.
122  *  </p>
123  * @author Maxime Journot
124  * @since 10.2
125  */
126 public class TimeSpanDragForce extends AbstractDragForceModel {
127 
128     /** Prefix for dates before in the parameter drivers' name. */
129     public static final String DATE_BEFORE = " - Before ";
130 
131     /** Prefix for dates after in the parameter drivers' name. */
132     public static final String DATE_AFTER = " - After ";
133 
134     /** TimeSpanMap of DragSensitive objects. */
135     private final TimeSpanMap<DragSensitive> dragSensitiveTimeSpanMap;
136 
137     /** Time scale used for the default names of the drag parameter drivers. */
138     private final TimeScale timeScale;
139 
140     /** Constructor with default UTC time scale for the default names of the drag parameter drivers.
141      * @param atmosphere atmospheric model
142      * @param spacecraft Time scale used for the default names of the drag parameter drivers
143      */
144     @DefaultDataContext
145     public TimeSpanDragForce(final Atmosphere atmosphere,
146                              final DragSensitive spacecraft) {
147         super(atmosphere);
148         this.dragSensitiveTimeSpanMap = new TimeSpanMap<>(spacecraft);
149         this.timeScale = TimeScalesFactory.getUTC();
150     }
151 
152     /** Constructor.
153      * @param atmosphere atmospheric model
154      * @param spacecraft the initial object physical and geometric information
155      * @param timeScale Time scale used for the default names of the drag parameter drivers
156      */
157     public TimeSpanDragForce(final Atmosphere atmosphere,
158                              final DragSensitive spacecraft,
159                              final TimeScale timeScale) {
160         super(atmosphere);
161         this.dragSensitiveTimeSpanMap = new TimeSpanMap<>(spacecraft);
162         this.timeScale = timeScale;
163     }
164 
165     /** Add a DragSensitive entry valid before a limit date.<br>
166      * Using <code>addDragSensitiveValidBefore(entry, t)</code> will make <code>entry</code>
167      * valid in ]-∞, t[ (note the open bracket).
168      * @param dragSensitive DragSensitive entry
169      * @param latestValidityDate date before which the entry is valid
170      * (must be different from <b>all</b> dates already used for transitions)
171      */
172     public void addDragSensitiveValidBefore(final DragSensitive dragSensitive, final AbsoluteDate latestValidityDate) {
173         dragSensitiveTimeSpanMap.addValidBefore(changeDragParameterDriversNames(dragSensitive,
174                                                                                 latestValidityDate,
175                                                                                 DATE_BEFORE),
176                                                 latestValidityDate, false);
177     }
178 
179     /** Add a DragSensitive entry valid after a limit date.<br>
180      * Using <code>addDragSensitiveValidAfter(entry, t)</code> will make <code>entry</code>
181      * valid in [t, +∞[ (note the closed bracket).
182      * @param dragSensitive DragSensitive entry
183      * @param earliestValidityDate date after which the entry is valid
184      * (must be different from <b>all</b> dates already used for transitions)
185      */
186     public void addDragSensitiveValidAfter(final DragSensitive dragSensitive, final AbsoluteDate earliestValidityDate) {
187         dragSensitiveTimeSpanMap.addValidAfter(changeDragParameterDriversNames(dragSensitive,
188                                                                                earliestValidityDate,
189                                                                                DATE_AFTER),
190                                                earliestValidityDate, false);
191     }
192 
193     /** Get the {@link DragSensitive} model valid at a date.
194      * @param date the date of validity
195      * @return the DragSensitive model valid at date
196      */
197     public DragSensitive getDragSensitive(final AbsoluteDate date) {
198         return dragSensitiveTimeSpanMap.get(date);
199     }
200 
201     /** Get the {@link DragSensitive} {@link Span} containing a specified date.
202      * @param date date belonging to the desired time span
203      * @return the DragSensitive time span containing the specified date
204      */
205     public Span<DragSensitive> getDragSensitiveSpan(final AbsoluteDate date) {
206         return dragSensitiveTimeSpanMap.getSpan(date);
207     }
208 
209     /** Extract a range of the {@link DragSensitive} map.
210      * <p>
211      * The object returned will be a new independent instance that will contain
212      * only the transitions that lie in the specified range.
213      * </p>
214      * See the {@link TimeSpanMap#extractRange TimeSpanMap.extractRange method} for more.
215      * @param start earliest date at which a transition is included in the range
216      * (may be set to {@link AbsoluteDate#PAST_INFINITY} to keep all early transitions)
217      * @param end latest date at which a transition is included in the r
218      * (may be set to {@link AbsoluteDate#FUTURE_INFINITY} to keep all late transitions)
219      * @return a new TimeSpanMap instance of DragSensitive with all transitions restricted to the specified range
220      */
221     public TimeSpanMap<DragSensitive> extractDragSensitiveRange(final AbsoluteDate start, final AbsoluteDate end) {
222         return dragSensitiveTimeSpanMap.extractRange(start, end);
223     }
224 
225     /** Get the first {@link Span time span} of the drag sensitive time span map.
226      * @return the first {@link Span time span} of the drag sensitive time span map
227      * @since 11.1
228      */
229     public Span<DragSensitive> getFirstSpan() {
230         return dragSensitiveTimeSpanMap.getFirstSpan();
231     }
232 
233     /** {@inheritDoc} */
234     @Override
235     public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {
236 
237         // Local atmospheric density
238         final AbsoluteDate date     = s.getDate();
239         final Frame        frame    = s.getFrame();
240         final Vector3D     position = s.getPosition();
241         final double rho    = getAtmosphere().getDensity(date, position, frame);
242 
243         // Spacecraft relative velocity with respect to the atmosphere
244         final Vector3D vAtm = getAtmosphere().getVelocity(date, position, frame);
245         final Vector3D relativeVelocity = vAtm.subtract(s.getPVCoordinates().getVelocity());
246 
247         // Extract the proper parameters valid at date from the input array
248         final double[] extractedParameters = extractParameters(parameters, date);
249 
250         // Compute and return drag acceleration
251         return getDragSensitive(date).dragAcceleration(s, rho, relativeVelocity, extractedParameters);
252 
253     }
254 
255     /** {@inheritDoc} */
256     @Override
257     public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
258                                                                              final T[] parameters) {
259         // Density and its derivatives
260         final T rho = getFieldDensity(s);
261 
262         // Spacecraft relative velocity with respect to the atmosphere
263         final FieldAbsoluteDate<T> date     = s.getDate();
264         final Frame                frame    = s.getFrame();
265         final FieldVector3D<T>     position = s.getPosition();
266         final FieldVector3D<T> vAtm = getAtmosphere().getVelocity(date, position, frame);
267         final FieldVector3D<T> relativeVelocity = vAtm.subtract(s.getPVCoordinates().getVelocity());
268 
269         // Extract the proper parameters valid at date from the input array
270         final T[] extractedParameters = extractParameters(parameters, date);
271 
272         // Compute and return drag acceleration
273         return getDragSensitive(date.toAbsoluteDate()).dragAcceleration(s, rho, relativeVelocity, extractedParameters);
274     }
275 
276     /**{@inheritDoc}
277      * <p>
278      * A date detector is used to cleanly stop the propagator and reset
279      * the state derivatives at transition dates.
280      * </p>
281      */
282     @Override
283     public Stream<EventDetector> getEventDetectors() {
284 
285         // Get the transitions' dates from the TimeSpanMap
286         final AbsoluteDate[] transitionDates = getTransitionDates();
287 
288         // create detector and return it in Stream
289         final DateDetector detector = getDateDetector(transitionDates);
290         return Stream.of(detector);
291     }
292 
293     /** {@inheritDoc}
294      * <p>
295      * A date detector is used to cleanly stop the propagator and reset
296      * the state derivatives at transition dates.
297      * </p>
298      */
299     @Override
300     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
301 
302         // Get the transitions' dates from the TimeSpanMap
303         final AbsoluteDate[] transitionDates = getTransitionDates();
304 
305         // create detector and return it in Stream
306         final FieldDateDetector<T> detector = getFieldDateDetector(field, transitionDates);
307         return Stream.of(detector);
308     }
309 
310     /** {@inheritDoc}
311      * <p>
312      * All the parameter drivers of all DragSensitive models are returned in an array.
313      * Models are ordered chronologically.
314      * </p>
315      */
316     @Override
317     public List<ParameterDriver> getParametersDrivers() {
318 
319         // Get all transitions from the TimeSpanMap
320         final List<ParameterDriver> listParameterDrivers = new ArrayList<>();
321 
322         // Loop on the spans
323         for (Span<DragSensitive> span = getFirstSpan(); span != null; span = span.next()) {
324             // Add all the parameter drivers of the span
325             for (ParameterDriver driver : span.getData().getDragParametersDrivers()) {
326                 // Add the driver only if the name does not exist already
327                 if (!findByName(listParameterDrivers, driver.getName())) {
328                     listParameterDrivers.add(driver);
329                 }
330             }
331         }
332 
333         // Return an array of parameter drivers with no duplicated name
334         return listParameterDrivers;
335 
336     }
337 
338     /** Extract the proper parameter drivers' values from the array in input of the
339      * {@link #acceleration(SpacecraftState, double[]) acceleration} method.
340      *  Parameters are filtered given an input date.
341      * @param parameters the input parameters array
342      * @param date the date
343      * @return the parameters given the date
344      */
345     public double[] extractParameters(final double[] parameters, final AbsoluteDate date) {
346 
347         // Get the drag parameter drivers of the date
348         final List<ParameterDriver> dragParameterDriver = getDragSensitive(date).getDragParametersDrivers();
349 
350         // Find out the indexes of the parameters in the whole array of parameters
351         final List<ParameterDriver> allParameters = getParametersDrivers();
352         final double[] outParameters = new double[dragParameterDriver.size()];
353         int index = 0;
354         for (int i = 0; i < allParameters.size(); i++) {
355             final String driverName = allParameters.get(i).getName();
356             for (ParameterDriver dragDriver : dragParameterDriver) {
357                 if (dragDriver.getName().equals(driverName)) {
358                     outParameters[index++] = parameters[i];
359                 }
360             }
361         }
362         return outParameters;
363     }
364 
365     /** Extract the proper parameter drivers' values from the array in input of the
366      * {@link #acceleration(FieldSpacecraftState, CalculusFieldElement[]) acceleration} method.
367      *  Parameters are filtered given an input date.
368      * @param parameters the input parameters array
369      * @param date the date
370      * @param <T> extends CalculusFieldElement
371      * @return the parameters given the date
372      */
373     public <T extends CalculusFieldElement<T>> T[] extractParameters(final T[] parameters,
374                                                                  final FieldAbsoluteDate<T> date) {
375 
376         // Get the drag parameter drivers of the date
377         final List<ParameterDriver> dragPD = getDragSensitive(date.toAbsoluteDate()).getDragParametersDrivers();
378 
379         // Find out the indexes of the parameters in the whole array of parameters
380         final List<ParameterDriver> allParameters = getParametersDrivers();
381         final T[] outParameters = MathArrays.buildArray(date.getField(), dragPD.size());
382         int index = 0;
383         for (int i = 0; i < allParameters.size(); i++) {
384             final String driverName = allParameters.get(i).getName();
385             for (ParameterDriver dragDriver : dragPD) {
386                 if (dragDriver.getName().equals(driverName)) {
387                     outParameters[index++] = parameters[i];
388                 }
389             }
390         }
391         return outParameters;
392     }
393 
394     /** Find if a parameter driver with a given name already exists in a list of parameter drivers.
395      * @param driversList the list of parameter drivers
396      * @param name the parameter driver's name to filter with
397      * @return true if the name was found, false otherwise
398      */
399     private boolean findByName(final List<ParameterDriver> driversList, final String name) {
400         return ParameterDriversProvider.findByName(driversList, name);
401     }
402 
403     /** Get the dates of the transitions for the drag sensitive models {@link TimeSpanMap}.
404      * @return dates of the transitions for the drag sensitive models {@link TimeSpanMap}
405      */
406     private AbsoluteDate[] getTransitionDates() {
407 
408         // Get all transitions
409         final List<AbsoluteDate> listDates = new ArrayList<>();
410 
411         // Extract all the transitions' dates
412         for (Transition<DragSensitive> transition = getFirstSpan().getEndTransition(); transition != null; transition = transition.next()) {
413             listDates.add(transition.getDate());
414         }
415         // Return the array of transition dates
416         return listDates.toArray(new AbsoluteDate[0]);
417     }
418 
419     /** Change the parameter drivers names of a {@link DragSensitive} model, if needed.
420      * <p>
421      * This is done to avoid that several parameter drivers have the same name.<br>
422      * It is done only if the user hasn't modify the DragSensitive parameter drivers default names.
423      * </p>
424      * @param dragSensitive the DragSensitive model
425      * @param date the date used in the parameter driver's name
426      * @param datePrefix the date prefix used in the parameter driver's name
427      * @return the DragSensitive with its drivers' names changed
428      */
429     private DragSensitive changeDragParameterDriversNames(final DragSensitive dragSensitive,
430                                                           final AbsoluteDate date,
431                                                           final String datePrefix) {
432         // Loop on the parameter drivers of the DragSensitive model
433         for (ParameterDriver driver: dragSensitive.getDragParametersDrivers()) {
434             final String driverName = driver.getName();
435 
436             // If the name is the default name for DragSensitive parameter drivers
437             // Modify the name to add the prefix and the date
438             if (driverName.equals(DragSensitive.GLOBAL_DRAG_FACTOR) ||
439                 driverName.equals(DragSensitive.DRAG_COEFFICIENT) ||
440                 driverName.equals(DragSensitive.LIFT_RATIO)) {
441                 driver.setName(driverName + datePrefix + date.toString(timeScale));
442             }
443         }
444         return dragSensitive;
445     }
446 
447 }