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.propagation.semianalytical.dsst;
18
19 import java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.Collection;
22 import java.util.Collections;
23 import java.util.HashSet;
24 import java.util.List;
25 import java.util.Set;
26
27 import org.hipparchus.CalculusFieldElement;
28 import org.hipparchus.Field;
29 import org.hipparchus.ode.FieldODEIntegrator;
30 import org.hipparchus.ode.FieldODEStateAndDerivative;
31 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
32 import org.hipparchus.ode.sampling.FieldODEStepHandler;
33 import org.hipparchus.util.MathArrays;
34 import org.orekit.annotation.DefaultDataContext;
35 import org.orekit.attitudes.AttitudeProvider;
36 import org.orekit.attitudes.FieldAttitude;
37 import org.orekit.data.DataContext;
38 import org.orekit.errors.OrekitException;
39 import org.orekit.errors.OrekitInternalError;
40 import org.orekit.errors.OrekitMessages;
41 import org.orekit.frames.Frame;
42 import org.orekit.orbits.FieldEquinoctialOrbit;
43 import org.orekit.orbits.FieldOrbit;
44 import org.orekit.orbits.OrbitType;
45 import org.orekit.orbits.PositionAngleType;
46 import org.orekit.propagation.CartesianToleranceProvider;
47 import org.orekit.propagation.FieldSpacecraftState;
48 import org.orekit.propagation.PropagationType;
49 import org.orekit.propagation.Propagator;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.propagation.ToleranceProvider;
52 import org.orekit.propagation.conversion.osc2mean.DSSTTheory;
53 import org.orekit.propagation.conversion.osc2mean.FixedPointConverter;
54 import org.orekit.propagation.conversion.osc2mean.MeanTheory;
55 import org.orekit.propagation.conversion.osc2mean.OsculatingToMeanConverter;
56 import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator;
57 import org.orekit.propagation.integration.FieldStateMapper;
58 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
59 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
60 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
61 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
62 import org.orekit.propagation.semianalytical.dsst.utilities.FieldFixedNumberInterpolationGrid;
63 import org.orekit.propagation.semianalytical.dsst.utilities.FieldInterpolationGrid;
64 import org.orekit.propagation.semianalytical.dsst.utilities.FieldMaxGapInterpolationGrid;
65 import org.orekit.time.AbsoluteDate;
66 import org.orekit.time.FieldAbsoluteDate;
67 import org.orekit.utils.FieldDataDictionary;
68 import org.orekit.utils.FieldPVCoordinates;
69 import org.orekit.utils.ParameterDriver;
70 import org.orekit.utils.ParameterObserver;
71 import org.orekit.utils.TimeSpanMap;
72
73 /**
74 * This class propagates {@link org.orekit.orbits.FieldOrbit orbits} using the DSST theory.
75 * <p>
76 * Whereas analytical propagators are configured only thanks to their various
77 * constructors and can be used immediately after construction, such a semianalytical
78 * propagator configuration involves setting several parameters between construction
79 * time and propagation time, just as numerical propagators.
80 * </p>
81 * <p>
82 * The configuration parameters that can be set are:
83 * </p>
84 * <ul>
85 * <li>the initial spacecraft state ({@link #setInitialState(FieldSpacecraftState)})</li>
86 * <li>the various force models ({@link #addForceModel(DSSTForceModel)},
87 * {@link #removeForceModels()})</li>
88 * <li>the discrete events that should be triggered during propagation (
89 * {@link #addEventDetector(org.orekit.propagation.events.FieldEventDetector)},
90 * {@link #clearEventsDetectors()})</li>
91 * <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
92 * </ul>
93 * <p>
94 * From these configuration parameters, only the initial state is mandatory.
95 * The default propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial}
96 * parameters with {@link PositionAngleType#TRUE true} longitude argument.
97 * The central attraction coefficient used to define the initial orbit will be used.
98 * However, specifying only the initial state would mean the propagator would use
99 * only Keplerian forces. In this case, the simpler
100 * {@link org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator}
101 * class would be more effective.
102 * </p>
103 * <p>
104 * The underlying numerical integrator set up in the constructor may also have
105 * its own configuration parameters. Typical configuration parameters for adaptive
106 * stepsize integrators are the min, max and perhaps start step size as well as
107 * the absolute and/or relative errors thresholds.
108 * </p>
109 * <p>
110 * The state that is seen by the integrator is a simple six elements double array.
111 * These six elements are:
112 * <ul>
113 * <li>the {@link org.orekit.orbits.FieldEquinoctialOrbit equinoctial orbit parameters}
114 * (a, e<sub>x</sub>, e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>m</sub>)
115 * in meters and radians,</li>
116 * </ul>
117 *
118 * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
119 * thus allowing a new propagation to be started from there without recomputing the part already performed.
120 * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
121 * </p>
122 * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
123 * thread-safe.</p>
124 *
125 * @see FieldSpacecraftState
126 * @see DSSTForceModel
127 * @author Romain Di Costanzo
128 * @author Pascal Parraud
129 * @since 10.0
130 * @param <T> type of the field elements
131 */
132 public class FieldDSSTPropagator<T extends CalculusFieldElement<T>> extends FieldAbstractIntegratedPropagator<T> {
133
134 /** Retrograde factor I.
135 * <p>
136 * DSST model needs equinoctial orbit as internal representation.
137 * Classical equinoctial elements have discontinuities when inclination
138 * is close to zero. In this representation, I = +1. <br>
139 * To avoid this discontinuity, another representation exists and equinoctial
140 * elements can be expressed in a different way, called "retrograde" orbit.
141 * This implies I = -1. <br>
142 * As Orekit doesn't implement the retrograde orbit, I is always set to +1.
143 * But for the sake of consistency with the theory, the retrograde factor
144 * has been kept in the formulas.
145 * </p>
146 */
147 private static final int I = 1;
148
149 /** Number of grid points per integration step to be used in interpolation of short periodics coefficients.*/
150 private static final int INTERPOLATION_POINTS_PER_STEP = 3;
151
152 /** Default value for epsilon. */
153 private static final double EPSILON_DEFAULT = 1.0e-13;
154
155 /** Default value for maxIterations. */
156 private static final int MAX_ITERATIONS_DEFAULT = 200;
157
158 /** Flag specifying whether the initial orbital state is given with osculating elements. */
159 private boolean initialIsOsculating;
160
161 /** Field used by this class.*/
162 private final Field<T> field;
163
164 /** Force models used to compute short periodic terms. */
165 private final List<DSSTForceModel> forceModels;
166
167 /** State mapper holding the force models. */
168 private FieldMeanPlusShortPeriodicMapper mapper;
169
170 /** Generator for the interpolation grid. */
171 private FieldInterpolationGrid<T> interpolationgrid;
172
173 /** Create a new instance of DSSTPropagator.
174 * <p>
175 * After creation, there are no perturbing forces at all.
176 * This means that if {@link #addForceModel addForceModel}
177 * is not called after creation, the integrated orbit will
178 * follow a Keplerian evolution only.
179 * </p>
180 *
181 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
182 *
183 * @param field field used by default
184 * @param integrator numerical integrator to use for propagation.
185 * @param propagationType type of orbit to output (mean or osculating).
186 * @see #FieldDSSTPropagator(Field, FieldODEIntegrator, PropagationType,
187 * AttitudeProvider)
188 */
189 @DefaultDataContext
190 public FieldDSSTPropagator(final Field<T> field, final FieldODEIntegrator<T> integrator, final PropagationType propagationType) {
191 this(field, integrator, propagationType,
192 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
193 }
194
195 /** Create a new instance of DSSTPropagator.
196 * <p>
197 * After creation, there are no perturbing forces at all.
198 * This means that if {@link #addForceModel addForceModel}
199 * is not called after creation, the integrated orbit will
200 * follow a Keplerian evolution only.
201 * </p>
202 * @param field field used by default
203 * @param integrator numerical integrator to use for propagation.
204 * @param propagationType type of orbit to output (mean or osculating).
205 * @param attitudeProvider attitude law to use.
206 * @since 10.1
207 */
208 public FieldDSSTPropagator(final Field<T> field,
209 final FieldODEIntegrator<T> integrator,
210 final PropagationType propagationType,
211 final AttitudeProvider attitudeProvider) {
212 super(field, integrator, propagationType);
213 this.field = field;
214 forceModels = new ArrayList<>();
215 initMapper(field);
216 // DSST uses only equinoctial orbits and mean longitude argument
217 setOrbitType(OrbitType.EQUINOCTIAL);
218 setPositionAngleType(PositionAngleType.MEAN);
219 setAttitudeProvider(attitudeProvider);
220 setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
221 }
222
223 /** Create a new instance of DSSTPropagator.
224 * <p>
225 * After creation, there are no perturbing forces at all.
226 * This means that if {@link #addForceModel addForceModel}
227 * is not called after creation, the integrated orbit will
228 * follow a Keplerian evolution only. Only the mean orbits
229 * will be generated.
230 * </p>
231 *
232 * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
233 *
234 * @param field fied used by default
235 * @param integrator numerical integrator to use for propagation.
236 * @see #FieldDSSTPropagator(Field, FieldODEIntegrator, AttitudeProvider)
237 */
238 @DefaultDataContext
239 public FieldDSSTPropagator(final Field<T> field, final FieldODEIntegrator<T> integrator) {
240 this(field, integrator,
241 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
242 }
243
244 /** Create a new instance of DSSTPropagator.
245 * <p>
246 * After creation, there are no perturbing forces at all.
247 * This means that if {@link #addForceModel addForceModel}
248 * is not called after creation, the integrated orbit will
249 * follow a Keplerian evolution only. Only the mean orbits
250 * will be generated.
251 * </p>
252 * @param field fied used by default
253 * @param integrator numerical integrator to use for propagation.
254 * @param attitudeProvider attitude law to use.
255 * @since 10.1
256 */
257 public FieldDSSTPropagator(final Field<T> field,
258 final FieldODEIntegrator<T> integrator,
259 final AttitudeProvider attitudeProvider) {
260 super(field, integrator, PropagationType.MEAN);
261 this.field = field;
262 forceModels = new ArrayList<>();
263 initMapper(field);
264 // DSST uses only equinoctial orbits and mean longitude argument
265 setOrbitType(OrbitType.EQUINOCTIAL);
266 setPositionAngleType(PositionAngleType.MEAN);
267 setAttitudeProvider(attitudeProvider);
268 setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
269 }
270
271 /** Set the central attraction coefficient μ.
272 * <p>
273 * Setting the central attraction coefficient is
274 * equivalent to {@link #addForceModel(DSSTForceModel) add}
275 * a {@link DSSTNewtonianAttraction} force model.
276 * </p>
277 * @param mu central attraction coefficient (m³/s²)
278 * @see #addForceModel(DSSTForceModel)
279 * @see #getAllForceModels()
280 */
281 @Override
282 public void setMu(final T mu) {
283 addForceModel(new DSSTNewtonianAttraction(mu.getReal()));
284 }
285
286 /** Set the central attraction coefficient μ only in upper class.
287 * @param mu central attraction coefficient (m³/s²)
288 */
289 private void superSetMu(final T mu) {
290 super.setMu(mu);
291 }
292
293 /** Check if Newtonian attraction force model is available.
294 * <p>
295 * Newtonian attraction is always the last force model in the list.
296 * </p>
297 * @return true if Newtonian attraction force model is available
298 */
299 private boolean hasNewtonianAttraction() {
300 final int last = forceModels.size() - 1;
301 return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
302 }
303
304 /** Set the initial state with osculating orbital elements.
305 * @param initialState initial state (defined with osculating elements)
306 */
307 public void setInitialState(final FieldSpacecraftState<T> initialState) {
308 setInitialState(initialState, PropagationType.OSCULATING);
309 }
310
311 /** Set the initial state.
312 * @param initialState initial state
313 * @param stateType defined if the orbital state is defined with osculating or mean elements
314 */
315 public void setInitialState(final FieldSpacecraftState<T> initialState,
316 final PropagationType stateType) {
317 resetInitialState(initialState, stateType);
318 }
319
320 /** Reset the initial state.
321 *
322 * @param state new initial state
323 */
324 @Override
325 public void resetInitialState(final FieldSpacecraftState<T> state) {
326 super.resetInitialState(state);
327 if (!hasNewtonianAttraction()) {
328 // use the state to define central attraction
329 setMu(state.getOrbit().getMu());
330 }
331 super.setStartDate(state.getDate());
332 }
333
334 /** {@inheritDoc}.
335 *
336 * <p>Change parameter {@link #initialIsOsculating()} accordingly
337 * @since 12.1.3
338 */
339 @Override
340 public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
341 // Reset initial state
342 resetInitialState(state);
343
344 // Change state of initial osculating, if needed
345 initialIsOsculating = stateType == PropagationType.OSCULATING;
346 }
347
348 /** Set the selected short periodic coefficients that must be stored as additional states.
349 * @param selectedCoefficients short periodic coefficients that must be stored as additional states
350 * (null means no coefficients are selected, empty set means all coefficients are selected)
351 */
352 public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
353 mapper.setSelectedCoefficients(selectedCoefficients == null ?
354 null : new HashSet<>(selectedCoefficients));
355 }
356
357 /** Get the selected short periodic coefficients that must be stored as additional states.
358 * @return short periodic coefficients that must be stored as additional states
359 * (null means no coefficients are selected, empty set means all coefficients are selected)
360 */
361 public Set<String> getSelectedCoefficients() {
362 final Set<String> set = mapper.getSelectedCoefficients();
363 return set == null ? null : Collections.unmodifiableSet(set);
364 }
365
366 /** Check if the initial state is provided in osculating elements.
367 * @return true if initial state is provided in osculating elements
368 */
369 public boolean initialIsOsculating() {
370 return initialIsOsculating;
371 }
372
373 /** Set the interpolation grid generator.
374 * <p>
375 * The generator will create an interpolation grid with a fixed
376 * number of points for each mean element integration step.
377 * </p>
378 * <p>
379 * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
380 * nor {@link #setInterpolationGridToMaxTimeGap(CalculusFieldElement)} has been called,
381 * by default the propagator is set as to 3 interpolations points per step.
382 * </p>
383 * @param interpolationPoints number of interpolation points at
384 * each integration step
385 * @see #setInterpolationGridToMaxTimeGap(CalculusFieldElement)
386 * @since 7.1
387 */
388 public void setInterpolationGridToFixedNumberOfPoints(final int interpolationPoints) {
389 interpolationgrid = new FieldFixedNumberInterpolationGrid<>(field, interpolationPoints);
390 }
391
392 /** Set the interpolation grid generator.
393 * <p>
394 * The generator will create an interpolation grid with a maximum
395 * time gap between interpolation points.
396 * </p>
397 * <p>
398 * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
399 * nor {@link #setInterpolationGridToMaxTimeGap(CalculusFieldElement)} has been called,
400 * by default the propagator is set as to 3 interpolations points per step.
401 * </p>
402 * @param maxGap maximum time gap between interpolation points (seconds)
403 * @see #setInterpolationGridToFixedNumberOfPoints(int)
404 * @since 7.1
405 */
406 public void setInterpolationGridToMaxTimeGap(final T maxGap) {
407 interpolationgrid = new FieldMaxGapInterpolationGrid<>(field, maxGap);
408 }
409
410 /** Add a force model to the global perturbation model.
411 * <p>
412 * If this method is not called at all,
413 * the integrated orbit will follow a Keplerian evolution only.
414 * </p>
415 * @param force perturbing {@link DSSTForceModel force} to add
416 * @see #removeForceModels()
417 * @see #setMu(CalculusFieldElement)
418 */
419 public void addForceModel(final DSSTForceModel force) {
420
421 if (force instanceof DSSTNewtonianAttraction) {
422 // we want to add the central attraction force model
423
424 try {
425 // ensure we are notified of any mu change
426 force.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
427 /** {@inheritDoc} */
428 @Override
429 public void valueChanged(final double previousValue, final ParameterDriver driver, final AbsoluteDate date) {
430 // mu PDriver should have only 1 span
431 superSetMu(field.getZero().newInstance(driver.getValue()));
432 }
433 /** {@inheritDoc} */
434 @Override
435 public void valueSpanMapChanged(final TimeSpanMap<Double> previousValue, final ParameterDriver driver) {
436 // mu PDriver should have only 1 span
437 superSetMu(field.getZero().newInstance(driver.getValue()));
438 }
439 });
440 } catch (OrekitException oe) {
441 // this should never happen
442 throw new OrekitInternalError(oe);
443 }
444
445 if (hasNewtonianAttraction()) {
446 // there is already a central attraction model, replace it
447 forceModels.set(forceModels.size() - 1, force);
448 } else {
449 // there are no central attraction model yet, add it at the end of the list
450 forceModels.add(force);
451 }
452 } else {
453 // we want to add a perturbing force model
454 if (hasNewtonianAttraction()) {
455 // insert the new force model before Newtonian attraction,
456 // which should always be the last one in the list
457 forceModels.add(forceModels.size() - 1, force);
458 } else {
459 // we only have perturbing force models up to now, just append at the end of the list
460 forceModels.add(force);
461 }
462 }
463
464 force.registerAttitudeProvider(getAttitudeProvider());
465
466 }
467
468 /** Remove all perturbing force models from the global perturbation model
469 * (except central attraction).
470 * <p>
471 * Once all perturbing forces have been removed (and as long as no new force model is added),
472 * the integrated orbit will follow a Keplerian evolution only.
473 * </p>
474 * @see #addForceModel(DSSTForceModel)
475 */
476 public void removeForceModels() {
477 final int last = forceModels.size() - 1;
478 if (hasNewtonianAttraction()) {
479 // preserve the Newtonian attraction model at the end
480 final DSSTForceModel newton = forceModels.get(last);
481 forceModels.clear();
482 forceModels.add(newton);
483 } else {
484 forceModels.clear();
485 }
486 }
487
488 /** Get all the force models, perturbing forces and Newtonian attraction included.
489 * @return list of perturbing force models, with Newtonian attraction being the
490 * last one
491 * @see #addForceModel(DSSTForceModel)
492 * @see #setMu(CalculusFieldElement)
493 */
494 public List<DSSTForceModel> getAllForceModels() {
495 return Collections.unmodifiableList(forceModels);
496 }
497
498 /** Get propagation parameter type.
499 * @return orbit type used for propagation
500 */
501 @Override
502 public OrbitType getOrbitType() {
503 return super.getOrbitType();
504 }
505
506 /** Get propagation parameter type.
507 * @return angle type to use for propagation
508 */
509 @Override
510 public PositionAngleType getPositionAngleType() {
511 return super.getPositionAngleType();
512 }
513
514 /** Conversion from mean to osculating orbit.
515 * <p>
516 * Compute osculating state <b>in a DSST sense</b>, corresponding to the
517 * mean SpacecraftState in input, and according to the Force models taken
518 * into account.
519 * </p><p>
520 * Since the osculating state is obtained by adding short-periodic variation
521 * of each force model, the resulting output will depend on the
522 * force models parameterized in input.
523 * </p>
524 * @param mean Mean state to convert
525 * @param forces Forces to take into account
526 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
527 * like atmospheric drag, radiation pressure or specific user-defined models)
528 * @param <T> type of the elements
529 * @return osculating state in a DSST sense
530 */
531 @SuppressWarnings("unchecked")
532 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeOsculatingState(final FieldSpacecraftState<T> mean,
533 final AttitudeProvider attitudeProvider,
534 final Collection<DSSTForceModel> forces) {
535
536 //Create the auxiliary object
537 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(mean.getOrbit(), I);
538
539 // Set the force models
540 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<>();
541 for (final DSSTForceModel force : forces) {
542 force.registerAttitudeProvider(attitudeProvider);
543 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters(mean.getDate().getField(), mean.getDate())));
544 force.updateShortPeriodTerms(force.getParametersAllValues(mean.getDate().getField()), mean);
545 }
546
547 final FieldEquinoctialOrbit<T> osculatingOrbit = computeOsculatingOrbit(mean, shortPeriodTerms);
548
549 return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
550 mean.getAdditionalDataValues(), null);
551
552 }
553
554 /** Conversion from osculating to mean orbit.
555 * <p>
556 * Compute mean state <b>in a DSST sense</b>, corresponding to the
557 * osculating SpacecraftState in input, and according to the Force models
558 * taken into account.
559 * </p><p>
560 * Since the osculating state is obtained with the computation of
561 * short-periodic variation of each force model, the resulting output will
562 * depend on the force models parameterized in input.
563 * </p><p>
564 * The computation is done through a fixed-point iteration process.
565 * </p>
566 * @param osculating Osculating state to convert
567 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
568 * like atmospheric drag, radiation pressure or specific user-defined models)
569 * @param forceModel Forces to take into account
570 * @param <T> type of the elements
571 * @return mean state in a DSST sense
572 */
573 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
574 final AttitudeProvider attitudeProvider,
575 final Collection<DSSTForceModel> forceModel) {
576 final OsculatingToMeanConverter converter = new FixedPointConverter(EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT,
577 FixedPointConverter.DEFAULT_DAMPING);
578 return computeMeanState(osculating, attitudeProvider, forceModel, converter);
579 }
580
581 /** Conversion from osculating to mean orbit.
582 * <p>
583 * Compute mean state <b>in a DSST sense</b>, corresponding to the
584 * osculating SpacecraftState in input, and according to the Force models
585 * taken into account.
586 * </p><p>
587 * Since the osculating state is obtained with the computation of
588 * short-periodic variation of each force model, the resulting output will
589 * depend on the force models parameterized in input.
590 * </p><p>
591 * The computation is done through a fixed-point iteration process.
592 * </p>
593 * @param osculating Osculating state to convert
594 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
595 * like atmospheric drag, radiation pressure or specific user-defined models)
596 * @param forceModel Forces to take into account
597 * @param epsilon convergence threshold for mean parameters conversion
598 * @param maxIterations maximum iterations for mean parameters conversion
599 * @return mean state in a DSST sense
600 * @param <T> type of the elements
601 * @since 10.1
602 */
603 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
604 final AttitudeProvider attitudeProvider,
605 final Collection<DSSTForceModel> forceModel,
606 final double epsilon,
607 final int maxIterations) {
608 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
609 FixedPointConverter.DEFAULT_DAMPING);
610 return computeMeanState(osculating, attitudeProvider, forceModel, converter);
611 }
612
613 /** Conversion from osculating to mean orbit.
614 * <p>
615 * Compute mean state <b>in a DSST sense</b>, corresponding to the
616 * osculating SpacecraftState in input, and according to the Force models
617 * taken into account.
618 * </p><p>
619 * Since the osculating state is obtained with the computation of
620 * short-periodic variation of each force model, the resulting output will
621 * depend on the force models parameterized in input.
622 * </p><p>
623 * The computation is done using the given osculating to mean orbit converter.
624 * </p>
625 * @param <T> type of the elements
626 * @param osculating osculating state to convert
627 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
628 * like atmospheric drag, radiation pressure or specific user-defined models)
629 * @param forceModel forces to take into account
630 * @param converter osculating to mean orbit converter
631 * @return mean state in a DSST sense
632 * @since 13.0
633 */
634 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
635 final AttitudeProvider attitudeProvider,
636 final Collection<DSSTForceModel> forceModel,
637 final OsculatingToMeanConverter converter) {
638
639 final MeanTheory theory = new DSSTTheory(forceModel, attitudeProvider, osculating.getMass().getReal());
640 converter.setMeanTheory(theory);
641 final FieldOrbit<T> meanOrbit = converter.convertToMean(osculating.getOrbit());
642 return new FieldSpacecraftState<>(meanOrbit, osculating.getAttitude(), osculating.getMass(),
643 osculating.getAdditionalDataValues(), null);
644 }
645
646 /** Override the default value of the parameter.
647 * <p>
648 * By default, if the initial orbit is defined as osculating,
649 * it will be averaged over 2 satellite revolutions.
650 * This can be changed by using this method.
651 * </p>
652 * @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
653 * elements
654 */
655 public void setSatelliteRevolution(final int satelliteRevolution) {
656 mapper.setSatelliteRevolution(satelliteRevolution);
657 }
658
659 /** Get the number of satellite revolutions to use for converting osculating to mean elements.
660 * @return number of satellite revolutions to use for converting osculating to mean elements
661 */
662 public int getSatelliteRevolution() {
663 return mapper.getSatelliteRevolution();
664 }
665
666 /** {@inheritDoc} */
667 @Override
668 public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
669 super.setAttitudeProvider(attitudeProvider);
670
671 //Register the attitude provider for each force model
672 for (final DSSTForceModel force : forceModels) {
673 force.registerAttitudeProvider(attitudeProvider);
674 }
675 }
676
677 /** Method called just before integration.
678 * <p>
679 * The default implementation does nothing, it may be specialized in subclasses.
680 * </p>
681 * @param initialState initial state
682 * @param tEnd target date at which state should be propagated
683 */
684 @SuppressWarnings("unchecked")
685 @Override
686 protected void beforeIntegration(final FieldSpacecraftState<T> initialState,
687 final FieldAbsoluteDate<T> tEnd) {
688
689 // check if only mean elements must be used
690 final PropagationType type = isMeanOrbit();
691
692 // compute common auxiliary elements
693 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(initialState.getOrbit(), I);
694
695 // initialize all perturbing forces
696 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<>();
697 for (final DSSTForceModel force : forceModels) {
698 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters(field, initialState.getDate())));
699 }
700 mapper.setShortPeriodTerms(shortPeriodTerms);
701
702 // if required, insert the special short periodics step handler
703 if (type == PropagationType.OSCULATING) {
704 final FieldShortPeriodicsHandler spHandler = new FieldShortPeriodicsHandler(forceModels);
705 // Compute short periodic coefficients for this point
706 for (DSSTForceModel forceModel : forceModels) {
707 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(field), initialState);
708
709 }
710 final Collection<FieldODEStepHandler<T>> stepHandlers = new ArrayList<>();
711 stepHandlers.add(spHandler);
712 final FieldODEIntegrator<T> integrator = getIntegrator();
713 final Collection<FieldODEStepHandler<T>> existing = integrator.getStepHandlers();
714 stepHandlers.addAll(existing);
715
716 integrator.clearStepHandlers();
717
718 // add back the existing handlers after the short periodics one
719 for (final FieldODEStepHandler<T> sp : stepHandlers) {
720 integrator.addStepHandler(sp);
721 }
722 }
723 }
724
725 /** {@inheritDoc} */
726 @Override
727 protected void afterIntegration() {
728 // remove the special short periodics step handler if added before
729 if (isMeanOrbit() == PropagationType.OSCULATING) {
730 final List<FieldODEStepHandler<T>> preserved = new ArrayList<>();
731 final FieldODEIntegrator<T> integrator = getIntegrator();
732
733 // clear the list
734 integrator.clearStepHandlers();
735
736 // add back the step handlers that were important for the user
737 for (final FieldODEStepHandler<T> sp : preserved) {
738 integrator.addStepHandler(sp);
739 }
740 }
741 }
742
743 /** Compute osculating state from mean state.
744 * <p>
745 * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
746 * </p>
747 * @param meanState initial mean state
748 * @param shortPeriodTerms short period terms
749 * @param <T> type of the elements
750 * @return osculating state
751 */
752 private static <T extends CalculusFieldElement<T>> FieldEquinoctialOrbit<T> computeOsculatingOrbit(final FieldSpacecraftState<T> meanState,
753 final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
754
755 final T[] mean = MathArrays.buildArray(meanState.getDate().getField(), 6);
756 final T[] meanDot = MathArrays.buildArray(meanState.getDate().getField(), 6);
757 OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngleType.MEAN, mean, meanDot);
758 final T[] y = mean.clone();
759 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
760 final T[] shortPeriodic = spt.value(meanState.getOrbit());
761 for (int i = 0; i < shortPeriodic.length; i++) {
762 y[i] = y[i].add(shortPeriodic[i]);
763 }
764 }
765 return (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
766 PositionAngleType.MEAN, meanState.getDate(),
767 meanState.getOrbit().getMu(), meanState.getFrame());
768 }
769
770 /** {@inheritDoc} */
771 @Override
772 protected FieldSpacecraftState<T> getInitialIntegrationState() {
773 if (initialIsOsculating) {
774 // the initial state is an osculating state,
775 // it must be converted to mean state
776 return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
777 } else {
778 // the initial state is already a mean state
779 return getInitialState();
780 }
781 }
782
783 /** {@inheritDoc}
784 * <p>
785 * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
786 * and position angle type is hardcoded to {@link PositionAngleType#MEAN}, so
787 * the corresponding parameters are ignored.
788 * </p>
789 */
790 @Override
791 protected FieldStateMapper<T> createMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
792 final OrbitType ignoredOrbitType, final PositionAngleType ignoredPositionAngleType,
793 final AttitudeProvider attitudeProvider, final Frame frame) {
794
795 // create a mapper with the common settings provided as arguments
796 final FieldMeanPlusShortPeriodicMapper newMapper =
797 new FieldMeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);
798
799 // copy the specific settings from the existing mapper
800 if (mapper != null) {
801 newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
802 newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
803 newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
804 }
805
806 mapper = newMapper;
807 return mapper;
808
809 }
810
811 /** Internal mapper using mean parameters plus short periodic terms. */
812 private class FieldMeanPlusShortPeriodicMapper extends FieldStateMapper<T> {
813
814 /** Short periodic coefficients that must be stored as additional states. */
815 private Set<String> selectedCoefficients;
816
817 /** Number of satellite revolutions in the averaging interval. */
818 private int satelliteRevolution;
819
820 /** Short period terms. */
821 private List<FieldShortPeriodTerms<T>> shortPeriodTerms;
822
823 /** Simple constructor.
824 * @param referenceDate reference date
825 * @param mu central attraction coefficient (m³/s²)
826 * @param attitudeProvider attitude provider
827 * @param frame inertial frame
828 */
829 FieldMeanPlusShortPeriodicMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
830 final AttitudeProvider attitudeProvider, final Frame frame) {
831
832 super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, attitudeProvider, frame);
833
834 this.selectedCoefficients = null;
835
836 // Default averaging period for conversion from osculating to mean elements
837 this.satelliteRevolution = 2;
838
839 this.shortPeriodTerms = Collections.emptyList();
840
841 }
842
843 /** {@inheritDoc} */
844 @Override
845 public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
846 final T[] y, final T[] yDot,
847 final PropagationType type) {
848
849 // add short periodic variations to mean elements to get osculating elements
850 // (the loop may not be performed if there are no force models and in the
851 // case we want to remain in mean parameters only)
852 final T[] elements = y.clone();
853 final FieldDataDictionary<T> coefficients;
854 switch (type) {
855 case MEAN:
856 coefficients = null;
857 break;
858 case OSCULATING:
859 final FieldOrbit<T> meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngleType.MEAN, date, getMu(), getFrame());
860 coefficients = selectedCoefficients == null ? null : new FieldDataDictionary<>(date.getField());
861 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
862 final T[] shortPeriodic = spt.value(meanOrbit);
863 for (int i = 0; i < shortPeriodic.length; i++) {
864 elements[i] = elements[i].add(shortPeriodic[i]);
865 }
866 if (selectedCoefficients != null) {
867 coefficients.putAllFields(spt.getCoefficients(date, selectedCoefficients));
868 }
869 }
870 break;
871 default:
872 throw new OrekitInternalError(null);
873 }
874
875 final T mass = elements[6];
876 if (mass.getReal() <= 0.0) {
877 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
878 }
879
880 final FieldOrbit<T> orbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngleType.MEAN, date, getMu(), getFrame());
881 final FieldAttitude<T> attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());
882
883 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients, null);
884
885 }
886
887 /** {@inheritDoc} */
888 @Override
889 public void mapStateToArray(final FieldSpacecraftState<T> state, final T[] y, final T[] yDot) {
890
891 OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, y, yDot);
892 y[6] = state.getMass();
893
894 }
895
896 /** Set the number of satellite revolutions to use for converting osculating to mean elements.
897 * <p>
898 * By default, if the initial orbit is defined as osculating,
899 * it will be averaged over 2 satellite revolutions.
900 * This can be changed by using this method.
901 * </p>
902 * @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
903 * elements
904 */
905 public void setSatelliteRevolution(final int satelliteRevolution) {
906 this.satelliteRevolution = satelliteRevolution;
907 }
908
909 /** Get the number of satellite revolutions to use for converting osculating to mean elements.
910 * @return number of satellite revolutions to use for converting osculating to mean elements
911 */
912 public int getSatelliteRevolution() {
913 return satelliteRevolution;
914 }
915
916 /** Set the selected short periodic coefficients that must be stored as additional states.
917 * @param selectedCoefficients short periodic coefficients that must be stored as additional states
918 * (null means no coefficients are selected, empty set means all coefficients are selected)
919 */
920 public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
921 this.selectedCoefficients = selectedCoefficients;
922 }
923
924 /** Get the selected short periodic coefficients that must be stored as additional states.
925 * @return short periodic coefficients that must be stored as additional states
926 * (null means no coefficients are selected, empty set means all coefficients are selected)
927 */
928 public Set<String> getSelectedCoefficients() {
929 return selectedCoefficients;
930 }
931
932 /** Set the short period terms.
933 * @param shortPeriodTerms short period terms
934 * @since 7.1
935 */
936 public void setShortPeriodTerms(final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
937 this.shortPeriodTerms = shortPeriodTerms;
938 }
939
940 /** Get the short period terms.
941 * @return shortPeriodTerms short period terms
942 * @since 7.1
943 */
944 public List<FieldShortPeriodTerms<T>> getShortPeriodTerms() {
945 return shortPeriodTerms;
946 }
947
948 }
949
950 /** {@inheritDoc} */
951 @Override
952 protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
953 return new Main(integrator);
954 }
955
956 /** Internal class for mean parameters integration. */
957 private class Main implements MainStateEquations<T> {
958
959 /** Derivatives array. */
960 private final T[] yDot;
961
962 /** Simple constructor.
963 * @param integrator numerical integrator to use for propagation.
964 */
965 Main(final FieldODEIntegrator<T> integrator) {
966 yDot = MathArrays.buildArray(field, 7);
967
968 // Setup event detectors for each force model
969 forceModels.forEach(dsstForceModel -> dsstForceModel.getFieldEventDetectors(field).
970 forEach(eventDetector -> setUpEventDetector(integrator, eventDetector)));
971 }
972
973 /** {@inheritDoc} */
974 @Override
975 public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
976 forceModels.forEach(fm -> fm.init(initialState, target));
977 }
978
979 /** {@inheritDoc} */
980 @Override
981 public T[] computeDerivatives(final FieldSpacecraftState<T> state) {
982
983 final T zero = state.getDate().getField().getZero();
984 Arrays.fill(yDot, zero);
985
986 // compute common auxiliary elements
987 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);
988
989 // compute the contributions of all perturbing forces
990 for (final DSSTForceModel forceModel : forceModels) {
991 final T[] daidt = elementRates(forceModel, state, auxiliaryElements, forceModel.getParametersAllValues(field));
992 for (int i = 0; i < daidt.length; i++) {
993 yDot[i] = yDot[i].add(daidt[i]);
994 }
995 }
996
997 return yDot.clone();
998 }
999
1000 /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
1001 * for a specific force model.
1002 * @param forceModel force to take into account
1003 * @param state current state
1004 * @param auxiliaryElements auxiliary elements related to the current orbit
1005 * @param parameters force model parameters
1006 * @return the mean equinoctial elements rates da<sub>i</sub> / dt
1007 */
1008 private T[] elementRates(final DSSTForceModel forceModel,
1009 final FieldSpacecraftState<T> state,
1010 final FieldAuxiliaryElements<T> auxiliaryElements,
1011 final T[] parameters) {
1012 return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
1013 }
1014
1015 }
1016
1017 /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
1018 * <p>
1019 * The errors are estimated from partial derivatives properties of orbits,
1020 * starting from a scalar position error specified by the user.
1021 * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
1022 * we get at constant energy (i.e. on a Keplerian trajectory):
1023 *
1024 * <pre>
1025 * V r² |dV| = mu |dr|
1026 * </pre>
1027 *
1028 * <p> So we deduce a scalar velocity error consistent with the position error. From here, we apply
1029 * orbits Jacobians matrices to get consistent errors on orbital parameters.
1030 *
1031 * <p>
1032 * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
1033 * local estimates, not global ones. So some care must be taken when using these tolerances.
1034 * Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
1035 * position after several orbits integration.
1036 * </p>
1037 * @param <T> elements type
1038 * @param dP user specified position error (m)
1039 * @param orbit reference orbit
1040 * @return a two rows array, row 0 being the absolute tolerance error
1041 * and row 1 being the relative tolerance error
1042 * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
1043 */
1044 @Deprecated
1045 public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final FieldOrbit<T> orbit) {
1046 // estimate the scalar velocity error
1047 final FieldPVCoordinates<T> pv = orbit.getPVCoordinates();
1048 final T r2 = pv.getPosition().getNormSq();
1049 final T v = pv.getVelocity().getNorm();
1050 final T dV = orbit.getMu().multiply(dP).divide(v.multiply(r2));
1051
1052 return tolerances(dP, dV, orbit);
1053 }
1054
1055 /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
1056 * <p>
1057 * The errors are estimated from partial derivatives properties of orbits,
1058 * starting from scalar position and velocity errors specified by the user.
1059 * <p>
1060 * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
1061 * local estimates, not global ones. So some care must be taken when using these tolerances.
1062 * Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
1063 * position after several orbits integration.
1064 * </p>
1065 *
1066 * @param <T> elements type
1067 * @param dP user specified position error (m)
1068 * @param dV user specified velocity error (m/s)
1069 * @param orbit reference orbit
1070 * @return a two rows array, row 0 being the absolute tolerance error
1071 * and row 1 being the relative tolerance error
1072 * @since 10.3
1073 * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
1074 */
1075 @Deprecated
1076 public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final T dV,
1077 final FieldOrbit<T> orbit) {
1078 return ToleranceProvider.of(CartesianToleranceProvider.of(dP.getReal(), dV.getReal(),
1079 CartesianToleranceProvider.DEFAULT_ABSOLUTE_MASS_TOLERANCE))
1080 .getTolerances(orbit, OrbitType.EQUINOCTIAL, PositionAngleType.MEAN);
1081 }
1082
1083 /** Step handler used to compute the parameters for the short periodic contributions.
1084 * @author Lucian Barbulescu
1085 */
1086 private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {
1087
1088 /** Force models used to compute short periodic terms. */
1089 private final List<DSSTForceModel> forceModels;
1090
1091 /** Constructor.
1092 * @param forceModels force models
1093 */
1094 FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
1095 this.forceModels = forceModels;
1096 }
1097
1098 /** {@inheritDoc} */
1099 @SuppressWarnings("unchecked")
1100 @Override
1101 public void handleStep(final FieldODEStateInterpolator<T> interpolator) {
1102
1103 // Get the grid points to compute
1104 final T[] interpolationPoints =
1105 interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
1106 interpolator.getCurrentState().getTime());
1107
1108 final FieldSpacecraftState<T>[] meanStates = new FieldSpacecraftState[interpolationPoints.length];
1109 for (int i = 0; i < interpolationPoints.length; ++i) {
1110
1111 // Build the mean state interpolated at grid point
1112 final T time = interpolationPoints[i];
1113 final FieldODEStateAndDerivative<T> sd = interpolator.getInterpolatedState(time);
1114 meanStates[i] = mapper.mapArrayToState(time,
1115 sd.getPrimaryState(),
1116 sd.getPrimaryDerivative(),
1117 PropagationType.MEAN);
1118
1119 }
1120
1121 // Compute short periodic coefficients for this step
1122 for (DSSTForceModel forceModel : forceModels) {
1123 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(field), meanStates);
1124 }
1125
1126 }
1127 }
1128
1129 }