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