1 /* Copyright 2002-2021 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.semianalytical.dsst;
18
19 import java.util.ArrayList;
20 import java.util.Arrays;
21 import java.util.Collection;
22 import java.util.Collections;
23 import java.util.HashMap;
24 import java.util.HashSet;
25 import java.util.List;
26 import java.util.Map;
27 import java.util.Set;
28
29 import org.hipparchus.Field;
30 import org.hipparchus.CalculusFieldElement;
31 import org.hipparchus.ode.FieldODEIntegrator;
32 import org.hipparchus.ode.FieldODEStateAndDerivative;
33 import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
34 import org.hipparchus.ode.sampling.FieldODEStepHandler;
35 import org.hipparchus.util.FastMath;
36 import org.hipparchus.util.MathArrays;
37 import org.hipparchus.util.MathUtils;
38 import org.orekit.annotation.DefaultDataContext;
39 import org.orekit.attitudes.AttitudeProvider;
40 import org.orekit.attitudes.FieldAttitude;
41 import org.orekit.data.DataContext;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.errors.OrekitInternalError;
44 import org.orekit.errors.OrekitMessages;
45 import org.orekit.frames.Frame;
46 import org.orekit.orbits.FieldEquinoctialOrbit;
47 import org.orekit.orbits.FieldOrbit;
48 import org.orekit.orbits.OrbitType;
49 import org.orekit.orbits.PositionAngle;
50 import org.orekit.propagation.FieldSpacecraftState;
51 import org.orekit.propagation.PropagationType;
52 import org.orekit.propagation.Propagator;
53 import org.orekit.propagation.SpacecraftState;
54 import org.orekit.propagation.events.FieldEventDetector;
55 import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator;
56 import org.orekit.propagation.integration.FieldStateMapper;
57 import org.orekit.propagation.numerical.FieldNumericalPropagator;
58 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
59 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
60 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
61 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
62 import org.orekit.propagation.semianalytical.dsst.utilities.FieldFixedNumberInterpolationGrid;
63 import org.orekit.propagation.semianalytical.dsst.utilities.FieldInterpolationGrid;
64 import org.orekit.propagation.semianalytical.dsst.utilities.FieldMaxGapInterpolationGrid;
65 import org.orekit.time.AbsoluteDate;
66 import org.orekit.time.FieldAbsoluteDate;
67 import org.orekit.utils.ParameterDriver;
68 import org.orekit.utils.ParameterObserver;
69
70 /**
71 * This class propagates {@link org.orekit.orbits.FieldOrbit 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(FieldSpacecraftState)})</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.FieldEventDetector)},
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 PositionAngle#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.FieldEquinoctialOrbit 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 FieldSpacecraftState
123 * @see DSSTForceModel
124 * @author Romain Di Costanzo
125 * @author Pascal Parraud
126 * @since 10.0
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<DSSTForceModel>();
211 initMapper(field);
212 // DSST uses only equinoctial orbits and mean longitude argument
213 setOrbitType(OrbitType.EQUINOCTIAL);
214 setPositionAngleType(PositionAngle.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<DSSTForceModel>();
259 initMapper(field);
260 // DSST uses only equinoctial orbits and mean longitude argument
261 setOrbitType(OrbitType.EQUINOCTIAL);
262 setPositionAngleType(PositionAngle.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) {
421 superSetMu(field.getZero().add(driver.getValue()));
422 }
423 });
424 } catch (OrekitException oe) {
425 // this should never happen
426 throw new OrekitInternalError(oe);
427 }
428
429 if (hasNewtonianAttraction()) {
430 // there is already a central attraction model, replace it
431 forceModels.set(forceModels.size() - 1, force);
432 } else {
433 // there are no central attraction model yet, add it at the end of the list
434 forceModels.add(force);
435 }
436 } else {
437 // we want to add a perturbing force model
438 if (hasNewtonianAttraction()) {
439 // insert the new force model before Newtonian attraction,
440 // which should always be the last one in the list
441 forceModels.add(forceModels.size() - 1, force);
442 } else {
443 // we only have perturbing force models up to now, just append at the end of the list
444 forceModels.add(force);
445 }
446 }
447
448 force.registerAttitudeProvider(getAttitudeProvider());
449
450 }
451
452 /** Remove all perturbing force models from the global perturbation model
453 * (except central attraction).
454 * <p>
455 * Once all perturbing forces have been removed (and as long as no new force model is added),
456 * the integrated orbit will follow a Keplerian evolution only.
457 * </p>
458 * @see #addForceModel(DSSTForceModel)
459 */
460 public void removeForceModels() {
461 final int last = forceModels.size() - 1;
462 if (hasNewtonianAttraction()) {
463 // preserve the Newtonian attraction model at the end
464 final DSSTForceModel newton = forceModels.get(last);
465 forceModels.clear();
466 forceModels.add(newton);
467 } else {
468 forceModels.clear();
469 }
470 }
471
472 /** Get all the force models, perturbing forces and Newtonian attraction included.
473 * @return list of perturbing force models, with Newtonian attraction being the
474 * last one
475 * @see #addForceModel(DSSTForceModel)
476 * @see #setMu(CalculusFieldElement)
477 */
478 public List<DSSTForceModel> getAllForceModels() {
479 return Collections.unmodifiableList(forceModels);
480 }
481
482 /** Get propagation parameter type.
483 * @return orbit type used for propagation
484 */
485 public OrbitType getOrbitType() {
486 return super.getOrbitType();
487 }
488
489 /** Get propagation parameter type.
490 * @return angle type to use for propagation
491 */
492 public PositionAngle getPositionAngleType() {
493 return super.getPositionAngleType();
494 }
495
496 /** Conversion from mean to osculating orbit.
497 * <p>
498 * Compute osculating state <b>in a DSST sense</b>, corresponding to the
499 * mean SpacecraftState in input, and according to the Force models taken
500 * into account.
501 * </p><p>
502 * Since the osculating state is obtained by adding short-periodic variation
503 * of each force model, the resulting output will depend on the
504 * force models parameterized in input.
505 * </p>
506 * @param mean Mean state to convert
507 * @param forces Forces to take into account
508 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
509 * like atmospheric drag, radiation pressure or specific user-defined models)
510 * @param <T> type of the elements
511 * @return osculating state in a DSST sense
512 */
513 @SuppressWarnings("unchecked")
514 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeOsculatingState(final FieldSpacecraftState<T> mean,
515 final AttitudeProvider attitudeProvider,
516 final Collection<DSSTForceModel> forces) {
517
518 //Create the auxiliary object
519 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(mean.getOrbit(), I);
520
521 // Set the force models
522 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
523 for (final DSSTForceModel force : forces) {
524 final T[] parameters = force.getParameters(mean.getDate().getField());
525 force.registerAttitudeProvider(attitudeProvider);
526 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, parameters));
527 force.updateShortPeriodTerms(parameters, mean);
528 }
529
530 final FieldEquinoctialOrbit<T> osculatingOrbit = computeOsculatingOrbit(mean, shortPeriodTerms);
531
532 return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
533 mean.getAdditionalStates());
534
535 }
536
537 /** Conversion from osculating to mean orbit.
538 * <p>
539 * Compute mean state <b>in a DSST sense</b>, corresponding to the
540 * osculating SpacecraftState in input, and according to the Force models
541 * taken into account.
542 * </p><p>
543 * Since the osculating state is obtained with the computation of
544 * short-periodic variation of each force model, the resulting output will
545 * depend on the force models parameterized in input.
546 * </p><p>
547 * The computation is done through a fixed-point iteration process.
548 * </p>
549 * @param osculating Osculating state to convert
550 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
551 * like atmospheric drag, radiation pressure or specific user-defined models)
552 * @param forceModel Forces to take into account
553 * @param <T> type of the elements
554 * @return mean state in a DSST sense
555 */
556 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
557 final AttitudeProvider attitudeProvider,
558 final Collection<DSSTForceModel> forceModel) {
559 return computeMeanState(osculating, attitudeProvider, forceModel, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
560 }
561
562 /** Conversion from osculating to mean orbit.
563 * <p>
564 * Compute mean state <b>in a DSST sense</b>, corresponding to the
565 * osculating SpacecraftState in input, and according to the Force models
566 * taken into account.
567 * </p><p>
568 * Since the osculating state is obtained with the computation of
569 * short-periodic variation of each force model, the resulting output will
570 * depend on the force models parameterized in input.
571 * </p><p>
572 * The computation is done through a fixed-point iteration process.
573 * </p>
574 * @param osculating Osculating state to convert
575 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
576 * like atmospheric drag, radiation pressure or specific user-defined models)
577 * @param forceModel Forces to take into account
578 * @param epsilon convergence threshold for mean parameters conversion
579 * @param maxIterations maximum iterations for mean parameters conversion
580 * @return mean state in a DSST sense
581 * @param <T> type of the elements
582 * @since 10.1
583 */
584 public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
585 final AttitudeProvider attitudeProvider,
586 final Collection<DSSTForceModel> forceModel,
587 final double epsilon,
588 final int maxIterations) {
589 final FieldOrbit<T> meanOrbit = computeMeanOrbit(osculating, attitudeProvider, forceModel, epsilon, maxIterations);
590 return new FieldSpacecraftState<>(meanOrbit, osculating.getAttitude(), osculating.getMass(), osculating.getAdditionalStates());
591 }
592
593 /** Override the default value of the parameter.
594 * <p>
595 * By default, if the initial orbit is defined as osculating,
596 * it will be averaged over 2 satellite revolutions.
597 * This can be changed by using this method.
598 * </p>
599 * @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
600 * elements
601 */
602 public void setSatelliteRevolution(final int satelliteRevolution) {
603 mapper.setSatelliteRevolution(satelliteRevolution);
604 }
605
606 /** Get the number of satellite revolutions to use for converting osculating to mean elements.
607 * @return number of satellite revolutions to use for converting osculating to mean elements
608 */
609 public int getSatelliteRevolution() {
610 return mapper.getSatelliteRevolution();
611 }
612
613 /** {@inheritDoc} */
614 @Override
615 public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
616 super.setAttitudeProvider(attitudeProvider);
617
618 //Register the attitude provider for each force model
619 for (final DSSTForceModel force : forceModels) {
620 force.registerAttitudeProvider(attitudeProvider);
621 }
622 }
623
624 /** Method called just before integration.
625 * <p>
626 * The default implementation does nothing, it may be specialized in subclasses.
627 * </p>
628 * @param initialState initial state
629 * @param tEnd target date at which state should be propagated
630 */
631 @SuppressWarnings("unchecked")
632 @Override
633 protected void beforeIntegration(final FieldSpacecraftState<T> initialState,
634 final FieldAbsoluteDate<T> tEnd) {
635
636 // check if only mean elements must be used
637 final PropagationType type = isMeanOrbit();
638
639 // compute common auxiliary elements
640 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(initialState.getOrbit(), I);
641
642 // initialize all perturbing forces
643 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
644 for (final DSSTForceModel force : forceModels) {
645 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters(field)));
646 }
647 mapper.setShortPeriodTerms(shortPeriodTerms);
648
649 // if required, insert the special short periodics step handler
650 if (type == PropagationType.OSCULATING) {
651 final FieldShortPeriodicsHandler spHandler = new FieldShortPeriodicsHandler(forceModels);
652 // Compute short periodic coefficients for this point
653 for (DSSTForceModel forceModel : forceModels) {
654 forceModel.updateShortPeriodTerms(forceModel.getParameters(field), initialState);
655
656 }
657 final Collection<FieldODEStepHandler<T>> stepHandlers = new ArrayList<FieldODEStepHandler<T>>();
658 stepHandlers.add(spHandler);
659 final FieldODEIntegrator<T> integrator = getIntegrator();
660 final Collection<FieldODEStepHandler<T>> existing = integrator.getStepHandlers();
661 stepHandlers.addAll(existing);
662
663 integrator.clearStepHandlers();
664
665 // add back the existing handlers after the short periodics one
666 for (final FieldODEStepHandler<T> sp : stepHandlers) {
667 integrator.addStepHandler(sp);
668 }
669 }
670 }
671
672 /** {@inheritDoc} */
673 @Override
674 protected void afterIntegration() {
675 // remove the special short periodics step handler if added before
676 if (isMeanOrbit() == PropagationType.OSCULATING) {
677 final List<FieldODEStepHandler<T>> preserved = new ArrayList<FieldODEStepHandler<T>>();
678 final FieldODEIntegrator<T> integrator = getIntegrator();
679
680 // clear the list
681 integrator.clearStepHandlers();
682
683 // add back the step handlers that were important for the user
684 for (final FieldODEStepHandler<T> sp : preserved) {
685 integrator.addStepHandler(sp);
686 }
687 }
688 }
689
690 /** Compute mean state from osculating state.
691 * <p>
692 * Compute in a DSST sense the mean state corresponding to the input osculating state.
693 * </p><p>
694 * The computing is done through a fixed-point iteration process.
695 * </p>
696 * @param osculating initial osculating state
697 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
698 * like atmospheric drag, radiation pressure or specific user-defined models)
699 * @param forceModel force models
700 * @param epsilon convergence threshold for mean parameters conversion
701 * @param maxIterations maximum iterations for mean parameters conversion
702 * @param <T> type of the elements
703 * @return mean state
704 * @since 10.1
705 */
706 @SuppressWarnings("unchecked")
707 private static <T extends CalculusFieldElement<T>> FieldOrbit<T> computeMeanOrbit(final FieldSpacecraftState<T> osculating, final AttitudeProvider attitudeProvider, final Collection<DSSTForceModel> forceModel,
708 final double epsilon, final int maxIterations) {
709
710 // zero
711 final T zero = osculating.getDate().getField().getZero();
712
713 // rough initialization of the mean parameters
714 FieldEquinoctialOrbit<T> meanOrbit = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.convertType(osculating.getOrbit());
715
716 // threshold for each parameter
717 final T epsilonT = zero.add(epsilon);
718 final T thresholdA = epsilonT.multiply(FastMath.abs(meanOrbit.getA()).add(1.));
719 final T thresholdE = epsilonT.multiply(meanOrbit.getE().add(1.));
720 final T thresholdI = epsilonT.multiply(meanOrbit.getI().add(1.));
721 final T thresholdL = epsilonT.multiply(zero.getPi());
722
723 // ensure all Gaussian force models can rely on attitude
724 for (final DSSTForceModel force : forceModel) {
725 force.registerAttitudeProvider(attitudeProvider);
726 }
727
728 int i = 0;
729 while (i++ < maxIterations) {
730
731 final FieldSpacecraftState<T> meanState = new FieldSpacecraftState<>(meanOrbit, osculating.getAttitude(), osculating.getMass());
732
733 //Create the auxiliary object
734 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(meanOrbit, I);
735
736 // Set the force models
737 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
738 for (final DSSTForceModel force : forceModel) {
739 final T[] parameters = force.getParameters(osculating.getDate().getField());
740 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, parameters));
741 force.updateShortPeriodTerms(parameters, meanState);
742 }
743
744 // recompute the osculating parameters from the current mean parameters
745 final FieldEquinoctialOrbit<T> rebuilt = computeOsculatingOrbit(meanState, shortPeriodTerms);
746
747 // adapted parameters residuals
748 final T deltaA = osculating.getA().subtract(rebuilt.getA());
749 final T deltaEx = osculating.getEquinoctialEx().subtract(rebuilt.getEquinoctialEx());
750 final T deltaEy = osculating.getEquinoctialEy().subtract(rebuilt.getEquinoctialEy());
751 final T deltaHx = osculating.getHx().subtract(rebuilt.getHx());
752 final T deltaHy = osculating.getHy().subtract(rebuilt.getHy());
753 final T deltaLv = MathUtils.normalizeAngle(osculating.getLv().subtract(rebuilt.getLv()), zero);
754
755 // check convergence
756 if (FastMath.abs(deltaA).getReal() < thresholdA.getReal() &&
757 FastMath.abs(deltaEx).getReal() < thresholdE.getReal() &&
758 FastMath.abs(deltaEy).getReal() < thresholdE.getReal() &&
759 FastMath.abs(deltaHx).getReal() < thresholdI.getReal() &&
760 FastMath.abs(deltaHy).getReal() < thresholdI.getReal() &&
761 FastMath.abs(deltaLv).getReal() < thresholdL.getReal()) {
762 return meanOrbit;
763 }
764
765 // update mean parameters
766 meanOrbit = new FieldEquinoctialOrbit<>(meanOrbit.getA().add(deltaA),
767 meanOrbit.getEquinoctialEx().add(deltaEx),
768 meanOrbit.getEquinoctialEy().add(deltaEy),
769 meanOrbit.getHx().add(deltaHx),
770 meanOrbit.getHy().add(deltaHy),
771 meanOrbit.getLv().add(deltaLv),
772 PositionAngle.TRUE, meanOrbit.getFrame(),
773 meanOrbit.getDate(), meanOrbit.getMu());
774 }
775
776 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
777 }
778
779 /** Compute osculating state from mean state.
780 * <p>
781 * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
782 * </p>
783 * @param meanState initial mean state
784 * @param shortPeriodTerms short period terms
785 * @param <T> type of the elements
786 * @return osculating state
787 */
788 private static <T extends CalculusFieldElement<T>> FieldEquinoctialOrbit<T> computeOsculatingOrbit(final FieldSpacecraftState<T> meanState,
789 final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
790
791 final T[] mean = MathArrays.buildArray(meanState.getDate().getField(), 6);
792 final T[] meanDot = MathArrays.buildArray(meanState.getDate().getField(), 6);
793 OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngle.MEAN, mean, meanDot);
794 final T[] y = mean.clone();
795 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
796 final T[] shortPeriodic = spt.value(meanState.getOrbit());
797 for (int i = 0; i < shortPeriodic.length; i++) {
798 y[i] = y[i].add(shortPeriodic[i]);
799 }
800 }
801 return (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
802 PositionAngle.MEAN, meanState.getDate(),
803 meanState.getMu(), meanState.getFrame());
804 }
805
806 /** {@inheritDoc} */
807 @Override
808 protected FieldSpacecraftState<T> getInitialIntegrationState() {
809 if (initialIsOsculating) {
810 // the initial state is an osculating state,
811 // it must be converted to mean state
812 return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
813 } else {
814 // the initial state is already a mean state
815 return getInitialState();
816 }
817 }
818
819 /** {@inheritDoc}
820 * <p>
821 * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
822 * and position angle type is hardcoded to {@link PositionAngle#MEAN}, so
823 * the corresponding parameters are ignored.
824 * </p>
825 */
826 @Override
827 protected FieldStateMapper<T> createMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
828 final OrbitType ignoredOrbitType, final PositionAngle ignoredPositionAngleType,
829 final AttitudeProvider attitudeProvider, final Frame frame) {
830
831 // create a mapper with the common settings provided as arguments
832 final FieldMeanPlusShortPeriodicMapper newMapper =
833 new FieldMeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);
834
835 // copy the specific settings from the existing mapper
836 if (mapper != null) {
837 newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
838 newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
839 newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
840 }
841
842 mapper = newMapper;
843 return mapper;
844
845 }
846
847 /** Internal mapper using mean parameters plus short periodic terms. */
848 private class FieldMeanPlusShortPeriodicMapper extends FieldStateMapper<T> {
849
850 /** Short periodic coefficients that must be stored as additional states. */
851 private Set<String> selectedCoefficients;
852
853 /** Number of satellite revolutions in the averaging interval. */
854 private int satelliteRevolution;
855
856 /** Short period terms. */
857 private List<FieldShortPeriodTerms<T>> shortPeriodTerms;
858
859 /** Simple constructor.
860 * @param referenceDate reference date
861 * @param mu central attraction coefficient (m³/s²)
862 * @param attitudeProvider attitude provider
863 * @param frame inertial frame
864 */
865 FieldMeanPlusShortPeriodicMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
866 final AttitudeProvider attitudeProvider, final Frame frame) {
867
868 super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngle.MEAN, attitudeProvider, frame);
869
870 this.selectedCoefficients = null;
871
872 // Default averaging period for conversion from osculating to mean elements
873 this.satelliteRevolution = 2;
874
875 this.shortPeriodTerms = Collections.emptyList();
876
877 }
878
879 /** {@inheritDoc} */
880 @Override
881 public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
882 final T[] y, final T[] yDot,
883 final PropagationType type) {
884
885 // add short periodic variations to mean elements to get osculating elements
886 // (the loop may not be performed if there are no force models and in the
887 // case we want to remain in mean parameters only)
888 final T[] elements = y.clone();
889 final Map<String, T[]> coefficients;
890 switch (type) {
891 case MEAN:
892 coefficients = null;
893 break;
894 case OSCULATING:
895 final FieldOrbit<T> meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngle.MEAN, date, getMu(), getFrame());
896 coefficients = selectedCoefficients == null ? null : new HashMap<String, T[]>();
897 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
898 final T[] shortPeriodic = spt.value(meanOrbit);
899 for (int i = 0; i < shortPeriodic.length; i++) {
900 elements[i] = elements[i].add(shortPeriodic[i]);
901 }
902 if (selectedCoefficients != null) {
903 coefficients.putAll(spt.getCoefficients(date, selectedCoefficients));
904 }
905 }
906 break;
907 default:
908 throw new OrekitInternalError(null);
909 }
910
911 final T mass = elements[6];
912 if (mass.getReal() <= 0.0) {
913 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
914 }
915
916 final FieldOrbit<T> orbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngle.MEAN, date, getMu(), getFrame());
917 final FieldAttitude<T> attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());
918
919 if (coefficients == null) {
920 return new FieldSpacecraftState<>(orbit, attitude, mass);
921 } else {
922 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients);
923 }
924
925 }
926
927 /** {@inheritDoc} */
928 @Override
929 public void mapStateToArray(final FieldSpacecraftState<T> state, final T[] y, final T[] yDot) {
930
931 OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, y, yDot);
932 y[6] = state.getMass();
933
934 }
935
936 /** Set the number of satellite revolutions to use for converting osculating to mean elements.
937 * <p>
938 * By default, if the initial orbit is defined as osculating,
939 * it will be averaged over 2 satellite revolutions.
940 * This can be changed by using this method.
941 * </p>
942 * @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
943 * elements
944 */
945 public void setSatelliteRevolution(final int satelliteRevolution) {
946 this.satelliteRevolution = satelliteRevolution;
947 }
948
949 /** Get the number of satellite revolutions to use for converting osculating to mean elements.
950 * @return number of satellite revolutions to use for converting osculating to mean elements
951 */
952 public int getSatelliteRevolution() {
953 return satelliteRevolution;
954 }
955
956 /** Set the selected short periodic coefficients that must be stored as additional states.
957 * @param selectedCoefficients short periodic coefficients that must be stored as additional states
958 * (null means no coefficients are selected, empty set means all coefficients are selected)
959 */
960 public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
961 this.selectedCoefficients = selectedCoefficients;
962 }
963
964 /** Get the selected short periodic coefficients that must be stored as additional states.
965 * @return short periodic coefficients that must be stored as additional states
966 * (null means no coefficients are selected, empty set means all coefficients are selected)
967 */
968 public Set<String> getSelectedCoefficients() {
969 return selectedCoefficients;
970 }
971
972 /** Set the short period terms.
973 * @param shortPeriodTerms short period terms
974 * @since 7.1
975 */
976 public void setShortPeriodTerms(final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
977 this.shortPeriodTerms = shortPeriodTerms;
978 }
979
980 /** Get the short period terms.
981 * @return shortPeriodTerms short period terms
982 * @since 7.1
983 */
984 public List<FieldShortPeriodTerms<T>> getShortPeriodTerms() {
985 return shortPeriodTerms;
986 }
987
988 }
989
990 /** {@inheritDoc} */
991 @Override
992 protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
993 return new Main(integrator);
994 }
995
996 /** Internal class for mean parameters integration. */
997 private class Main implements MainStateEquations<T> {
998
999 /** Derivatives array. */
1000 private final T[] yDot;
1001
1002 /** Simple constructor.
1003 * @param integrator numerical integrator to use for propagation.
1004 */
1005 Main(final FieldODEIntegrator<T> integrator) {
1006 yDot = MathArrays.buildArray(field, 7);
1007
1008 for (final DSSTForceModel forceModel : forceModels) {
1009 final FieldEventDetector<T>[] modelDetectors = forceModel.getFieldEventsDetectors(field);
1010 if (modelDetectors != null) {
1011 for (final FieldEventDetector<T> detector : modelDetectors) {
1012 setUpEventDetector(integrator, detector);
1013 }
1014 }
1015 }
1016
1017 }
1018
1019 /** {@inheritDoc} */
1020 @Override
1021 public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
1022 final SpacecraftState stateD = initialState.toSpacecraftState();
1023 final AbsoluteDate targetD = target.toAbsoluteDate();
1024 for (final DSSTForceModel forceModel : forceModels) {
1025 forceModel.init(stateD, targetD);
1026 }
1027 }
1028
1029 /** {@inheritDoc} */
1030 @Override
1031 public T[] computeDerivatives(final FieldSpacecraftState<T> state) {
1032
1033 final T zero = state.getDate().getField().getZero();
1034 Arrays.fill(yDot, zero);
1035
1036 // compute common auxiliary elements
1037 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);
1038
1039 // compute the contributions of all perturbing forces
1040 for (final DSSTForceModel forceModel : forceModels) {
1041 final T[] daidt = elementRates(forceModel, state, auxiliaryElements, forceModel.getParameters(field));
1042 for (int i = 0; i < daidt.length; i++) {
1043 yDot[i] = yDot[i].add(daidt[i]);
1044 }
1045 }
1046
1047 return yDot.clone();
1048 }
1049
1050 /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
1051 * for a specific force model.
1052 * @param forceModel force to take into account
1053 * @param state current state
1054 * @param auxiliaryElements auxiliary elements related to the current orbit
1055 * @param parameters force model parameters
1056 * @return the mean equinoctial elements rates da<sub>i</sub> / dt
1057 */
1058 private T[] elementRates(final DSSTForceModel forceModel,
1059 final FieldSpacecraftState<T> state,
1060 final FieldAuxiliaryElements<T> auxiliaryElements,
1061 final T[] parameters) {
1062 return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
1063 }
1064
1065 }
1066
1067 /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
1068 * <p>
1069 * The errors are estimated from partial derivatives properties of orbits,
1070 * starting from a scalar position error specified by the user.
1071 * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
1072 * we get at constant energy (i.e. on a Keplerian trajectory):
1073 *
1074 * <pre>
1075 * V² r |dV| = mu |dr|
1076 * </pre>
1077 *
1078 * <p> So we deduce a scalar velocity error consistent with the position error. From here, we apply
1079 * orbits Jacobians matrices to get consistent errors on orbital parameters.
1080 *
1081 * <p>
1082 * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
1083 * local estimates, not global ones. So some care must be taken when using these tolerances.
1084 * Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
1085 * position after several orbits integration.
1086 * </p>
1087 * @param <T> elements type
1088 * @param dP user specified position error (m)
1089 * @param orbit reference orbit
1090 * @return a two rows array, row 0 being the absolute tolerance error
1091 * and row 1 being the relative tolerance error
1092 */
1093 public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final FieldOrbit<T> orbit) {
1094 return FieldNumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
1095 }
1096
1097 /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
1098 * <p>
1099 * The errors are estimated from partial derivatives properties of orbits,
1100 * starting from scalar position and velocity errors specified by the user.
1101 * <p>
1102 * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
1103 * local estimates, not global ones. So some care must be taken when using these tolerances.
1104 * Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
1105 * position after several orbits integration.
1106 * </p>
1107 *
1108 * @param <T> elements type
1109 * @param dP user specified position error (m)
1110 * @param dV user specified velocity error (m/s)
1111 * @param orbit reference orbit
1112 * @return a two rows array, row 0 being the absolute tolerance error
1113 * and row 1 being the relative tolerance error
1114 * @since 10.3
1115 */
1116 public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final T dV,
1117 final FieldOrbit<T> orbit) {
1118 return FieldNumericalPropagator.tolerances(dP, dV, orbit, OrbitType.EQUINOCTIAL);
1119 }
1120
1121 /** Step handler used to compute the parameters for the short periodic contributions.
1122 * @author Lucian Barbulescu
1123 */
1124 private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {
1125
1126 /** Force models used to compute short periodic terms. */
1127 private final List<DSSTForceModel> forceModels;
1128
1129 /** Constructor.
1130 * @param forceModels force models
1131 */
1132 FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
1133 this.forceModels = forceModels;
1134 }
1135
1136 /** {@inheritDoc} */
1137 @SuppressWarnings("unchecked")
1138 @Override
1139 public void handleStep(final FieldODEStateInterpolator<T> interpolator) {
1140
1141 // Get the grid points to compute
1142 final T[] interpolationPoints =
1143 interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
1144 interpolator.getCurrentState().getTime());
1145
1146 final FieldSpacecraftState<T>[] meanStates = new FieldSpacecraftState[interpolationPoints.length];
1147 for (int i = 0; i < interpolationPoints.length; ++i) {
1148
1149 // Build the mean state interpolated at grid point
1150 final T time = interpolationPoints[i];
1151 final FieldODEStateAndDerivative<T> sd = interpolator.getInterpolatedState(time);
1152 meanStates[i] = mapper.mapArrayToState(time,
1153 sd.getPrimaryState(),
1154 sd.getPrimaryDerivative(),
1155 PropagationType.MEAN);
1156
1157 }
1158
1159 // Compute short periodic coefficients for this step
1160 for (DSSTForceModel forceModel : forceModels) {
1161 forceModel.updateShortPeriodTerms(forceModel.getParameters(field), meanStates);
1162 }
1163
1164 }
1165 }
1166
1167 }