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