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