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