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