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.analytical;
18
19 import java.io.Serializable;
20
21 import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.linear.RealMatrix;
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.MathUtils;
27 import org.hipparchus.util.SinCos;
28 import org.orekit.attitudes.AttitudeProvider;
29 import org.orekit.attitudes.InertialProvider;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitMessages;
32 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
33 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
34 import org.orekit.orbits.CartesianOrbit;
35 import org.orekit.orbits.CircularOrbit;
36 import org.orekit.orbits.Orbit;
37 import org.orekit.orbits.OrbitType;
38 import org.orekit.orbits.PositionAngle;
39 import org.orekit.propagation.AbstractMatricesHarvester;
40 import org.orekit.propagation.PropagationType;
41 import org.orekit.propagation.SpacecraftState;
42 import org.orekit.time.AbsoluteDate;
43 import org.orekit.utils.DoubleArrayDictionary;
44 import org.orekit.utils.TimeSpanMap;
45 import org.orekit.utils.TimeStampedPVCoordinates;
46
47 /** This class propagates a {@link org.orekit.propagation.SpacecraftState}
48 * using the analytical Eckstein-Hechler model.
49 * <p>The Eckstein-Hechler model is suited for near circular orbits
50 * (e < 0.1, with poor accuracy between 0.005 and 0.1) and inclination
51 * neither equatorial (direct or retrograde) nor critical (direct or
52 * retrograde).</p>
53 * <p>
54 * Note that before version 7.0, there was a large inconsistency in the generated
55 * orbits, and it was fixed as of version 7.0 of Orekit, with a visible side effect.
56 * The problems is that if the circular parameters produced by the Eckstein-Hechler
57 * model are used to build an orbit considered to be osculating, the velocity deduced
58 * from this orbit was <em>inconsistent with the position evolution</em>! The reason is
59 * that the model includes non-Keplerian effects but it does not include a corresponding
60 * circular/Cartesian conversion. As a consequence, all subsequent computation involving
61 * velocity were wrong. This includes attitude modes like yaw compensation and Doppler
62 * effect. As this effect was considered serious enough and as accurate velocities were
63 * considered important, the propagator now generates {@link CartesianOrbit Cartesian
64 * orbits} which are built in a special way to ensure consistency throughout propagation.
65 * A side effect is that if circular parameters are rebuilt by user from these propagated
66 * Cartesian orbit, the circular parameters will generally <em>not</em> match the initial
67 * orbit (differences in semi-major axis can exceed 120 m). The position however <em>will</em>
68 * match to sub-micrometer level, and this position will be identical to the positions
69 * that were generated by previous versions (in other words, the internals of the models
70 * have not been changed, only the output parameters have been changed). The correctness
71 * of the initialization has been assessed and is good, as it allows the subsequent orbit
72 * to remain close to a numerical reference orbit.
73 * </p>
74 * <p>
75 * If users need a more definitive initialization of an Eckstein-Hechler propagator, they
76 * should consider using a {@link org.orekit.propagation.conversion.PropagatorConverter
77 * propagator converter} to initialize their Eckstein-Hechler propagator using a complete
78 * sample instead of just a single initial orbit.
79 * </p>
80 * @see Orbit
81 * @author Guylaine Prat
82 */
83 public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator {
84
85 /** Initial Eckstein-Hechler model. */
86 private EHModel initialModel;
87
88 /** All models. */
89 private transient TimeSpanMap<EHModel> models;
90
91 /** Reference radius of the central body attraction model (m). */
92 private double referenceRadius;
93
94 /** Central attraction coefficient (m³/s²). */
95 private double mu;
96
97 /** Un-normalized zonal coefficients. */
98 private double[] ck0;
99
100 /** Build a propagator from orbit and potential provider.
101 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
102 *
103 * <p>Using this constructor, an initial osculating orbit is considered.</p>
104 *
105 * @param initialOrbit initial orbit
106 * @param provider for un-normalized zonal coefficients
107 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider,
108 * UnnormalizedSphericalHarmonicsProvider)
109 * @see #EcksteinHechlerPropagator(Orbit, UnnormalizedSphericalHarmonicsProvider,
110 * PropagationType)
111 */
112 public EcksteinHechlerPropagator(final Orbit initialOrbit,
113 final UnnormalizedSphericalHarmonicsProvider provider) {
114 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
115 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
116 }
117
118 /**
119 * Private helper constructor.
120 * <p>Using this constructor, an initial osculating orbit is considered.</p>
121 * @param initialOrbit initial orbit
122 * @param attitude attitude provider
123 * @param mass spacecraft mass
124 * @param provider for un-normalized zonal coefficients
125 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
126 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
127 * UnnormalizedSphericalHarmonicsProvider,
128 * UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
129 * PropagationType)
130 */
131 public EcksteinHechlerPropagator(final Orbit initialOrbit,
132 final AttitudeProvider attitude,
133 final double mass,
134 final UnnormalizedSphericalHarmonicsProvider provider,
135 final UnnormalizedSphericalHarmonics harmonics) {
136 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
137 harmonics.getUnnormalizedCnm(2, 0),
138 harmonics.getUnnormalizedCnm(3, 0),
139 harmonics.getUnnormalizedCnm(4, 0),
140 harmonics.getUnnormalizedCnm(5, 0),
141 harmonics.getUnnormalizedCnm(6, 0));
142 }
143
144 /** Build a propagator from orbit and potential.
145 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
146 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
147 * are related to both the normalized coefficients
148 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
149 * and the J<sub>n</sub> one as follows:</p>
150 *
151 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
152 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
153 *
154 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
155 *
156 * <p>Using this constructor, an initial osculating orbit is considered.</p>
157 *
158 * @param initialOrbit initial orbit
159 * @param referenceRadius reference radius of the Earth for the potential model (m)
160 * @param mu central attraction coefficient (m³/s²)
161 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
162 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
163 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
164 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
165 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
166 * @see org.orekit.utils.Constants
167 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
168 * double, double, double, double, double)
169 */
170 public EcksteinHechlerPropagator(final Orbit initialOrbit,
171 final double referenceRadius, final double mu,
172 final double c20, final double c30, final double c40,
173 final double c50, final double c60) {
174 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
175 DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
176 }
177
178 /** Build a propagator from orbit, mass and potential provider.
179 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
180 *
181 * <p>Using this constructor, an initial osculating orbit is considered.</p>
182 *
183 * @param initialOrbit initial orbit
184 * @param mass spacecraft mass
185 * @param provider for un-normalized zonal coefficients
186 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
187 * UnnormalizedSphericalHarmonicsProvider)
188 */
189 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
190 final UnnormalizedSphericalHarmonicsProvider provider) {
191 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
192 mass, provider, provider.onDate(initialOrbit.getDate()));
193 }
194
195 /** Build a propagator from orbit, mass and potential.
196 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
197 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
198 * are related to both the normalized coefficients
199 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
200 * and the J<sub>n</sub> one as follows:</p>
201 *
202 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
203 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
204 *
205 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
206 *
207 * <p>Using this constructor, an initial osculating orbit is considered.</p>
208 *
209 * @param initialOrbit initial orbit
210 * @param mass spacecraft mass
211 * @param referenceRadius reference radius of the Earth for the potential model (m)
212 * @param mu central attraction coefficient (m³/s²)
213 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
214 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
215 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
216 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
217 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
218 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
219 * double, double, double, double, double)
220 */
221 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
222 final double referenceRadius, final double mu,
223 final double c20, final double c30, final double c40,
224 final double c50, final double c60) {
225 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
226 mass, referenceRadius, mu, c20, c30, c40, c50, c60);
227 }
228
229 /** Build a propagator from orbit, attitude provider and potential provider.
230 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
231 * <p>Using this constructor, an initial osculating orbit is considered.</p>
232 * @param initialOrbit initial orbit
233 * @param attitudeProv attitude provider
234 * @param provider for un-normalized zonal coefficients
235 */
236 public EcksteinHechlerPropagator(final Orbit initialOrbit,
237 final AttitudeProvider attitudeProv,
238 final UnnormalizedSphericalHarmonicsProvider provider) {
239 this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
240 }
241
242 /** Build a propagator from orbit, attitude provider and potential.
243 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
244 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
245 * are related to both the normalized coefficients
246 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
247 * and the J<sub>n</sub> one as follows:</p>
248 *
249 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
250 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
251 *
252 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
253 *
254 * <p>Using this constructor, an initial osculating orbit is considered.</p>
255 *
256 * @param initialOrbit initial orbit
257 * @param attitudeProv attitude provider
258 * @param referenceRadius reference radius of the Earth for the potential model (m)
259 * @param mu central attraction coefficient (m³/s²)
260 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
261 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
262 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
263 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
264 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
265 */
266 public EcksteinHechlerPropagator(final Orbit initialOrbit,
267 final AttitudeProvider attitudeProv,
268 final double referenceRadius, final double mu,
269 final double c20, final double c30, final double c40,
270 final double c50, final double c60) {
271 this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
272 }
273
274 /** Build a propagator from orbit, attitude provider, mass and potential provider.
275 * <p>Using this constructor, an initial osculating orbit is considered.</p>
276 * @param initialOrbit initial orbit
277 * @param attitudeProv attitude provider
278 * @param mass spacecraft mass
279 * @param provider for un-normalized zonal coefficients
280 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
281 * UnnormalizedSphericalHarmonicsProvider, PropagationType)
282 */
283 public EcksteinHechlerPropagator(final Orbit initialOrbit,
284 final AttitudeProvider attitudeProv,
285 final double mass,
286 final UnnormalizedSphericalHarmonicsProvider provider) {
287 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()));
288 }
289
290 /** Build a propagator from orbit, attitude provider, mass and potential.
291 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
292 * are related to both the normalized coefficients
293 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
294 * and the J<sub>n</sub> one as follows:</p>
295 *
296 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
297 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
298 *
299 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
300 *
301 * <p>Using this constructor, an initial osculating orbit is considered.</p>
302 *
303 * @param initialOrbit initial orbit
304 * @param attitudeProv attitude provider
305 * @param mass spacecraft mass
306 * @param referenceRadius reference radius of the Earth for the potential model (m)
307 * @param mu central attraction coefficient (m³/s²)
308 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
309 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
310 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
311 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
312 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
313 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
314 * double, double, double, double, double, PropagationType)
315 */
316 public EcksteinHechlerPropagator(final Orbit initialOrbit,
317 final AttitudeProvider attitudeProv,
318 final double mass,
319 final double referenceRadius, final double mu,
320 final double c20, final double c30, final double c40,
321 final double c50, final double c60) {
322 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, c60,
323 PropagationType.OSCULATING);
324 }
325
326
327 /** Build a propagator from orbit and potential provider.
328 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
329 *
330 * <p>Using this constructor, it is possible to define the initial orbit as
331 * a mean Eckstein-Hechler orbit or an osculating one.</p>
332 *
333 * @param initialOrbit initial orbit
334 * @param provider for un-normalized zonal coefficients
335 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
336 * @since 10.2
337 */
338 public EcksteinHechlerPropagator(final Orbit initialOrbit,
339 final UnnormalizedSphericalHarmonicsProvider provider,
340 final PropagationType initialType) {
341 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
342 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType);
343 }
344
345 /** Build a propagator from orbit, attitude provider, mass and potential provider.
346 * <p>Using this constructor, it is possible to define the initial orbit as
347 * a mean Eckstein-Hechler orbit or an osculating one.</p>
348 * @param initialOrbit initial orbit
349 * @param attitudeProv attitude provider
350 * @param mass spacecraft mass
351 * @param provider for un-normalized zonal coefficients
352 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
353 * @since 10.2
354 */
355 public EcksteinHechlerPropagator(final Orbit initialOrbit,
356 final AttitudeProvider attitudeProv,
357 final double mass,
358 final UnnormalizedSphericalHarmonicsProvider provider,
359 final PropagationType initialType) {
360 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType);
361 }
362
363 /**
364 * Private helper constructor.
365 * <p>Using this constructor, it is possible to define the initial orbit as
366 * a mean Eckstein-Hechler orbit or an osculating one.</p>
367 * @param initialOrbit initial orbit
368 * @param attitude attitude provider
369 * @param mass spacecraft mass
370 * @param provider for un-normalized zonal coefficients
371 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
372 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
373 * @since 10.2
374 */
375 public EcksteinHechlerPropagator(final Orbit initialOrbit,
376 final AttitudeProvider attitude,
377 final double mass,
378 final UnnormalizedSphericalHarmonicsProvider provider,
379 final UnnormalizedSphericalHarmonics harmonics,
380 final PropagationType initialType) {
381 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
382 harmonics.getUnnormalizedCnm(2, 0),
383 harmonics.getUnnormalizedCnm(3, 0),
384 harmonics.getUnnormalizedCnm(4, 0),
385 harmonics.getUnnormalizedCnm(5, 0),
386 harmonics.getUnnormalizedCnm(6, 0),
387 initialType);
388 }
389
390 /** Build a propagator from orbit, attitude provider, mass and potential.
391 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
392 * are related to both the normalized coefficients
393 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
394 * and the J<sub>n</sub> one as follows:</p>
395 *
396 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
397 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
398 *
399 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
400 *
401 * <p>Using this constructor, it is possible to define the initial orbit as
402 * a mean Eckstein-Hechler orbit or an osculating one.</p>
403 *
404 * @param initialOrbit initial orbit
405 * @param attitudeProv attitude provider
406 * @param mass spacecraft mass
407 * @param referenceRadius reference radius of the Earth for the potential model (m)
408 * @param mu central attraction coefficient (m³/s²)
409 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
410 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
411 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
412 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
413 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
414 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
415 * @since 10.2
416 */
417 public EcksteinHechlerPropagator(final Orbit initialOrbit,
418 final AttitudeProvider attitudeProv,
419 final double mass,
420 final double referenceRadius, final double mu,
421 final double c20, final double c30, final double c40,
422 final double c50, final double c60,
423 final PropagationType initialType) {
424
425 super(attitudeProv);
426
427 // store model coefficients
428 this.referenceRadius = referenceRadius;
429 this.mu = mu;
430 this.ck0 = new double[] {
431 0.0, 0.0, c20, c30, c40, c50, c60
432 };
433
434 // compute mean parameters if needed
435 // transform into circular adapted parameters used by the Eckstein-Hechler model
436 resetInitialState(new SpacecraftState(initialOrbit,
437 attitudeProv.getAttitude(initialOrbit,
438 initialOrbit.getDate(),
439 initialOrbit.getFrame()),
440 mass),
441 initialType);
442
443 }
444
445 /** {@inheritDoc}
446 * <p>The new initial state to consider
447 * must be defined with an osculating orbit.</p>
448 * @see #resetInitialState(SpacecraftState, PropagationType)
449 */
450 public void resetInitialState(final SpacecraftState state) {
451 resetInitialState(state, PropagationType.OSCULATING);
452 }
453
454 /** Reset the propagator initial state.
455 * @param state new initial state to consider
456 * @param stateType mean Eckstein-Hechler orbit or osculating orbit
457 * @since 10.2
458 */
459 public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
460 super.resetInitialState(state);
461 final CircularOrbit circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
462 this.initialModel = (stateType == PropagationType.MEAN) ?
463 new EHModel(circular, state.getMass(), referenceRadius, mu, ck0) :
464 computeMeanParameters(circular, state.getMass());
465 this.models = new TimeSpanMap<EHModel>(initialModel);
466 }
467
468 /** {@inheritDoc} */
469 protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
470 final EHModel newModel = computeMeanParameters((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit()),
471 state.getMass());
472 if (forward) {
473 models.addValidAfter(newModel, state.getDate(), false);
474 } else {
475 models.addValidBefore(newModel, state.getDate(), false);
476 }
477 stateChanged(state);
478 }
479
480 /** Compute mean parameters according to the Eckstein-Hechler analytical model.
481 * @param osculating osculating orbit
482 * @param mass constant mass
483 * @return Eckstein-Hechler mean model
484 */
485 private EHModel computeMeanParameters(final CircularOrbit osculating, final double mass) {
486
487 // sanity check
488 if (osculating.getA() < referenceRadius) {
489 throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
490 osculating.getA());
491 }
492
493 // rough initialization of the mean parameters
494 EHModel current = new EHModel(osculating, mass, referenceRadius, mu, ck0);
495
496 // threshold for each parameter
497 final double epsilon = 1.0e-13;
498 final double thresholdA = epsilon * (1 + FastMath.abs(current.mean.getA()));
499 final double thresholdE = epsilon * (1 + current.mean.getE());
500 final double thresholdAngles = epsilon * FastMath.PI;
501
502 int i = 0;
503 while (i++ < 100) {
504
505 // recompute the osculating parameters from the current mean parameters
506 final UnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate());
507
508 // adapted parameters residuals
509 final double deltaA = osculating.getA() - parameters[0].getValue();
510 final double deltaEx = osculating.getCircularEx() - parameters[1].getValue();
511 final double deltaEy = osculating.getCircularEy() - parameters[2].getValue();
512 final double deltaI = osculating.getI() - parameters[3].getValue();
513 final double deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
514 parameters[4].getValue(),
515 0.0);
516 final double deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM() - parameters[5].getValue(), 0.0);
517
518 // update mean parameters
519 current = new EHModel(new CircularOrbit(current.mean.getA() + deltaA,
520 current.mean.getCircularEx() + deltaEx,
521 current.mean.getCircularEy() + deltaEy,
522 current.mean.getI() + deltaI,
523 current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
524 current.mean.getAlphaM() + deltaAlphaM,
525 PositionAngle.MEAN,
526 current.mean.getFrame(),
527 current.mean.getDate(), mu),
528 mass, referenceRadius, mu, ck0);
529
530 // check convergence
531 if (FastMath.abs(deltaA) < thresholdA &&
532 FastMath.abs(deltaEx) < thresholdE &&
533 FastMath.abs(deltaEy) < thresholdE &&
534 FastMath.abs(deltaI) < thresholdAngles &&
535 FastMath.abs(deltaRAAN) < thresholdAngles &&
536 FastMath.abs(deltaAlphaM) < thresholdAngles) {
537 return current;
538 }
539
540 }
541
542 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECKSTEIN_HECHLER_MEAN_PARAMETERS, i);
543
544 }
545
546 /** {@inheritDoc} */
547 public CartesianOrbit propagateOrbit(final AbsoluteDate date) {
548 // compute Cartesian parameters, taking derivatives into account
549 // to make sure velocity and acceleration are consistent
550 final EHModel current = models.get(date);
551 return new CartesianOrbit(toCartesian(date, current.propagateParameters(date)),
552 current.mean.getFrame(), mu);
553 }
554
555 /**
556 * Get the central attraction coefficient μ.
557 * @return mu central attraction coefficient (m³/s²)
558 * @since 11.1
559 */
560 public double getMu() {
561 return mu;
562 }
563
564 /**
565 * Get the un-normalized zonal coefficients.
566 * @return the un-normalized zonal coefficients
567 * @since 11.1
568 */
569 public double[] getCk0() {
570 return ck0.clone();
571 }
572
573 /**
574 * Get the reference radius of the central body attraction model.
575 * @return the reference radius in meters
576 * @since 11.1
577 */
578 public double getReferenceRadius() {
579 return referenceRadius;
580 }
581
582 /** {@inheritDoc} */
583 @Override
584 protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
585 final DoubleArrayDictionary initialJacobianColumns) {
586 return new EcksteinHechlerHarvester(this, stmName, initialStm, initialJacobianColumns);
587 }
588
589 /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
590 private static class EHModel implements Serializable {
591
592 /** Serializable UID. */
593 private static final long serialVersionUID = 20160115L;
594
595 /** Mean orbit. */
596 private final CircularOrbit mean;
597
598 /** Constant mass. */
599 private final double mass;
600
601 // CHECKSTYLE: stop JavadocVariable check
602
603 // preprocessed values
604 private final double xnotDot;
605 private final double rdpom;
606 private final double rdpomp;
607 private final double eps1;
608 private final double eps2;
609 private final double xim;
610 private final double ommD;
611 private final double rdl;
612 private final double aMD;
613
614 private final double kh;
615 private final double kl;
616
617 private final double ax1;
618 private final double ay1;
619 private final double as1;
620 private final double ac2;
621 private final double axy3;
622 private final double as3;
623 private final double ac4;
624 private final double as5;
625 private final double ac6;
626
627 private final double ex1;
628 private final double exx2;
629 private final double exy2;
630 private final double ex3;
631 private final double ex4;
632
633 private final double ey1;
634 private final double eyx2;
635 private final double eyy2;
636 private final double ey3;
637 private final double ey4;
638
639 private final double rx1;
640 private final double ry1;
641 private final double r2;
642 private final double r3;
643 private final double rl;
644
645 private final double iy1;
646 private final double ix1;
647 private final double i2;
648 private final double i3;
649 private final double ih;
650
651 private final double lx1;
652 private final double ly1;
653 private final double l2;
654 private final double l3;
655 private final double ll;
656
657 // CHECKSTYLE: resume JavadocVariable check
658
659 /** Create a model for specified mean orbit.
660 * @param mean mean orbit
661 * @param mass constant mass
662 * @param referenceRadius reference radius of the central body attraction model (m)
663 * @param mu central attraction coefficient (m³/s²)
664 * @param ck0 un-normalized zonal coefficients
665 */
666 EHModel(final CircularOrbit mean, final double mass,
667 final double referenceRadius, final double mu, final double[] ck0) {
668
669 this.mean = mean;
670 this.mass = mass;
671
672 // preliminary processing
673 double q = referenceRadius / mean.getA();
674 double ql = q * q;
675 final double g2 = ck0[2] * ql;
676 ql *= q;
677 final double g3 = ck0[3] * ql;
678 ql *= q;
679 final double g4 = ck0[4] * ql;
680 ql *= q;
681 final double g5 = ck0[5] * ql;
682 ql *= q;
683 final double g6 = ck0[6] * ql;
684
685 final SinCos sc = FastMath.sinCos(mean.getI());
686 final double cosI1 = sc.cos();
687 final double sinI1 = sc.sin();
688 final double sinI2 = sinI1 * sinI1;
689 final double sinI4 = sinI2 * sinI2;
690 final double sinI6 = sinI2 * sinI4;
691
692 if (sinI2 < 1.0e-10) {
693 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
694 FastMath.toDegrees(mean.getI()));
695 }
696
697 if (FastMath.abs(sinI2 - 4.0 / 5.0) < 1.0e-3) {
698 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
699 FastMath.toDegrees(mean.getI()));
700 }
701
702 if (mean.getE() > 0.1) {
703 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
704 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
705 mean.getE());
706 }
707
708 xnotDot = FastMath.sqrt(mu / mean.getA()) / mean.getA();
709
710 rdpom = -0.75 * g2 * (4.0 - 5.0 * sinI2);
711 rdpomp = 7.5 * g4 * (1.0 - 31.0 / 8.0 * sinI2 + 49.0 / 16.0 * sinI4) -
712 13.125 * g6 * (1.0 - 8.0 * sinI2 + 129.0 / 8.0 * sinI4 - 297.0 / 32.0 * sinI6);
713
714 q = 3.0 / (32.0 * rdpom);
715 eps1 = q * g4 * sinI2 * (30.0 - 35.0 * sinI2) -
716 175.0 * q * g6 * sinI2 * (1.0 - 3.0 * sinI2 + 2.0625 * sinI4);
717 q = 3.0 * sinI1 / (8.0 * rdpom);
718 eps2 = q * g3 * (4.0 - 5.0 * sinI2) - q * g5 * (10.0 - 35.0 * sinI2 + 26.25 * sinI4);
719
720 xim = mean.getI();
721 ommD = cosI1 * (1.50 * g2 - 2.25 * g2 * g2 * (2.5 - 19.0 / 6.0 * sinI2) +
722 0.9375 * g4 * (7.0 * sinI2 - 4.0) +
723 3.28125 * g6 * (2.0 - 9.0 * sinI2 + 8.25 * sinI4));
724
725 rdl = 1.0 - 1.50 * g2 * (3.0 - 4.0 * sinI2);
726 aMD = rdl +
727 2.25 * g2 * g2 * (9.0 - 263.0 / 12.0 * sinI2 + 341.0 / 24.0 * sinI4) +
728 15.0 / 16.0 * g4 * (8.0 - 31.0 * sinI2 + 24.5 * sinI4) +
729 105.0 / 32.0 * g6 * (-10.0 / 3.0 + 25.0 * sinI2 - 48.75 * sinI4 + 27.5 * sinI6);
730
731 final double qq = -1.5 * g2 / rdl;
732 final double qA = 0.75 * g2 * g2 * sinI2;
733 final double qB = 0.25 * g4 * sinI2;
734 final double qC = 105.0 / 16.0 * g6 * sinI2;
735 final double qD = -0.75 * g3 * sinI1;
736 final double qE = 3.75 * g5 * sinI1;
737 kh = 0.375 / rdpom;
738 kl = kh / sinI1;
739
740 ax1 = qq * (2.0 - 3.5 * sinI2);
741 ay1 = qq * (2.0 - 2.5 * sinI2);
742 as1 = qD * (4.0 - 5.0 * sinI2) +
743 qE * (2.625 * sinI4 - 3.5 * sinI2 + 1.0);
744 ac2 = qq * sinI2 +
745 qA * 7.0 * (2.0 - 3.0 * sinI2) +
746 qB * (15.0 - 17.5 * sinI2) +
747 qC * (3.0 * sinI2 - 1.0 - 33.0 / 16.0 * sinI4);
748 axy3 = qq * 3.5 * sinI2;
749 as3 = qD * 5.0 / 3.0 * sinI2 +
750 qE * 7.0 / 6.0 * sinI2 * (1.0 - 1.125 * sinI2);
751 ac4 = qA * sinI2 +
752 qB * 4.375 * sinI2 +
753 qC * 0.75 * (1.1 * sinI4 - sinI2);
754
755 as5 = qE * 21.0 / 80.0 * sinI4;
756
757 ac6 = qC * -11.0 / 80.0 * sinI4;
758
759 ex1 = qq * (1.0 - 1.25 * sinI2);
760 exx2 = qq * 0.5 * (3.0 - 5.0 * sinI2);
761 exy2 = qq * (2.0 - 1.5 * sinI2);
762 ex3 = qq * 7.0 / 12.0 * sinI2;
763 ex4 = qq * 17.0 / 8.0 * sinI2;
764
765 ey1 = qq * (1.0 - 1.75 * sinI2);
766 eyx2 = qq * (1.0 - 3.0 * sinI2);
767 eyy2 = qq * (2.0 * sinI2 - 1.5);
768 ey3 = qq * 7.0 / 12.0 * sinI2;
769 ey4 = qq * 17.0 / 8.0 * sinI2;
770
771 q = -qq * cosI1;
772 rx1 = 3.5 * q;
773 ry1 = -2.5 * q;
774 r2 = -0.5 * q;
775 r3 = 7.0 / 6.0 * q;
776 rl = g3 * cosI1 * (4.0 - 15.0 * sinI2) -
777 2.5 * g5 * cosI1 * (4.0 - 42.0 * sinI2 + 52.5 * sinI4);
778
779 q = 0.5 * qq * sinI1 * cosI1;
780 iy1 = q;
781 ix1 = -q;
782 i2 = q;
783 i3 = q * 7.0 / 3.0;
784 ih = -g3 * cosI1 * (4.0 - 5.0 * sinI2) +
785 2.5 * g5 * cosI1 * (4.0 - 14.0 * sinI2 + 10.5 * sinI4);
786
787 lx1 = qq * (7.0 - 77.0 / 8.0 * sinI2);
788 ly1 = qq * (55.0 / 8.0 * sinI2 - 7.50);
789 l2 = qq * (1.25 * sinI2 - 0.5);
790 l3 = qq * (77.0 / 24.0 * sinI2 - 7.0 / 6.0);
791 ll = g3 * (53.0 * sinI2 - 4.0 - 57.5 * sinI4) +
792 2.5 * g5 * (4.0 - 96.0 * sinI2 + 269.5 * sinI4 - 183.75 * sinI6);
793
794 }
795
796 /** Extrapolate an orbit up to a specific target date.
797 * @param date target date for the orbit
798 * @return propagated parameters
799 */
800 public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) {
801
802 // Keplerian evolution
803 final UnivariateDerivative2 dt = new UnivariateDerivative2(date.durationFrom(mean.getDate()), 1.0, 0.0);
804 final UnivariateDerivative2 xnot = dt.multiply(xnotDot);
805
806 // secular effects
807
808 // eccentricity
809 final UnivariateDerivative2 x = xnot.multiply(rdpom + rdpomp);
810 final UnivariateDerivative2 cx = x.cos();
811 final UnivariateDerivative2 sx = x.sin();
812 final UnivariateDerivative2 exm = cx.multiply(mean.getCircularEx()).
813 add(sx.multiply(eps2 - (1.0 - eps1) * mean.getCircularEy()));
814 final UnivariateDerivative2 eym = sx.multiply((1.0 + eps1) * mean.getCircularEx()).
815 add(cx.multiply(mean.getCircularEy() - eps2)).
816 add(eps2);
817
818 // no secular effect on inclination
819
820 // right ascension of ascending node
821 final UnivariateDerivative2 omm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ommD * xnot.getValue(),
822 FastMath.PI),
823 ommD * xnotDot,
824 0.0);
825
826 // latitude argument
827 final UnivariateDerivative2 xlm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getAlphaM() + aMD * xnot.getValue(),
828 FastMath.PI),
829 aMD * xnotDot,
830 0.0);
831
832 // periodical terms
833 final UnivariateDerivative2 cl1 = xlm.cos();
834 final UnivariateDerivative2 sl1 = xlm.sin();
835 final UnivariateDerivative2 cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
836 final UnivariateDerivative2 sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
837 final UnivariateDerivative2 cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
838 final UnivariateDerivative2 sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
839 final UnivariateDerivative2 cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
840 final UnivariateDerivative2 sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
841 final UnivariateDerivative2 cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
842 final UnivariateDerivative2 sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
843 final UnivariateDerivative2 cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));
844
845 final UnivariateDerivative2 qh = eym.subtract(eps2).multiply(kh);
846 final UnivariateDerivative2 ql = exm.multiply(kl);
847
848 final UnivariateDerivative2 exmCl1 = exm.multiply(cl1);
849 final UnivariateDerivative2 exmSl1 = exm.multiply(sl1);
850 final UnivariateDerivative2 eymCl1 = eym.multiply(cl1);
851 final UnivariateDerivative2 eymSl1 = eym.multiply(sl1);
852 final UnivariateDerivative2 exmCl2 = exm.multiply(cl2);
853 final UnivariateDerivative2 exmSl2 = exm.multiply(sl2);
854 final UnivariateDerivative2 eymCl2 = eym.multiply(cl2);
855 final UnivariateDerivative2 eymSl2 = eym.multiply(sl2);
856 final UnivariateDerivative2 exmCl3 = exm.multiply(cl3);
857 final UnivariateDerivative2 exmSl3 = exm.multiply(sl3);
858 final UnivariateDerivative2 eymCl3 = eym.multiply(cl3);
859 final UnivariateDerivative2 eymSl3 = eym.multiply(sl3);
860 final UnivariateDerivative2 exmCl4 = exm.multiply(cl4);
861 final UnivariateDerivative2 exmSl4 = exm.multiply(sl4);
862 final UnivariateDerivative2 eymCl4 = eym.multiply(cl4);
863 final UnivariateDerivative2 eymSl4 = eym.multiply(sl4);
864
865 // semi major axis
866 final UnivariateDerivative2 rda = exmCl1.multiply(ax1).
867 add(eymSl1.multiply(ay1)).
868 add(sl1.multiply(as1)).
869 add(cl2.multiply(ac2)).
870 add(exmCl3.add(eymSl3).multiply(axy3)).
871 add(sl3.multiply(as3)).
872 add(cl4.multiply(ac4)).
873 add(sl5.multiply(as5)).
874 add(cl6.multiply(ac6));
875
876 // eccentricity
877 final UnivariateDerivative2 rdex = cl1.multiply(ex1).
878 add(exmCl2.multiply(exx2)).
879 add(eymSl2.multiply(exy2)).
880 add(cl3.multiply(ex3)).
881 add(exmCl4.add(eymSl4).multiply(ex4));
882 final UnivariateDerivative2 rdey = sl1.multiply(ey1).
883 add(exmSl2.multiply(eyx2)).
884 add(eymCl2.multiply(eyy2)).
885 add(sl3.multiply(ey3)).
886 add(exmSl4.subtract(eymCl4).multiply(ey4));
887
888 // ascending node
889 final UnivariateDerivative2 rdom = exmSl1.multiply(rx1).
890 add(eymCl1.multiply(ry1)).
891 add(sl2.multiply(r2)).
892 add(eymCl3.subtract(exmSl3).multiply(r3)).
893 add(ql.multiply(rl));
894
895 // inclination
896 final UnivariateDerivative2 rdxi = eymSl1.multiply(iy1).
897 add(exmCl1.multiply(ix1)).
898 add(cl2.multiply(i2)).
899 add(exmCl3.add(eymSl3).multiply(i3)).
900 add(qh.multiply(ih));
901
902 // latitude argument
903 final UnivariateDerivative2 rdxl = exmSl1.multiply(lx1).
904 add(eymCl1.multiply(ly1)).
905 add(sl2.multiply(l2)).
906 add(exmSl3.subtract(eymCl3).multiply(l3)).
907 add(ql.multiply(ll));
908
909 // osculating parameters
910 return new UnivariateDerivative2[] {
911 rda.add(1.0).multiply(mean.getA()),
912 rdex.add(exm),
913 rdey.add(eym),
914 rdxi.add(xim),
915 rdom.add(omm),
916 rdxl.add(xlm)
917 };
918
919 }
920
921 }
922
923 /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
924 * @param date date of the orbital parameters
925 * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
926 * @return Cartesian coordinates consistent with values and derivatives
927 */
928 private TimeStampedPVCoordinates toCartesian(final AbsoluteDate date, final UnivariateDerivative2[] parameters) {
929
930 // evaluate coordinates in the orbit canonical reference frame
931 final UnivariateDerivative2 cosOmega = parameters[4].cos();
932 final UnivariateDerivative2 sinOmega = parameters[4].sin();
933 final UnivariateDerivative2 cosI = parameters[3].cos();
934 final UnivariateDerivative2 sinI = parameters[3].sin();
935 final UnivariateDerivative2 alphaE = meanToEccentric(parameters[5], parameters[1], parameters[2]);
936 final UnivariateDerivative2 cosAE = alphaE.cos();
937 final UnivariateDerivative2 sinAE = alphaE.sin();
938 final UnivariateDerivative2 ex2 = parameters[1].multiply(parameters[1]);
939 final UnivariateDerivative2 ey2 = parameters[2].multiply(parameters[2]);
940 final UnivariateDerivative2 exy = parameters[1].multiply(parameters[2]);
941 final UnivariateDerivative2 q = ex2.add(ey2).subtract(1).negate().sqrt();
942 final UnivariateDerivative2 beta = q.add(1).reciprocal();
943 final UnivariateDerivative2 bx2 = beta.multiply(ex2);
944 final UnivariateDerivative2 by2 = beta.multiply(ey2);
945 final UnivariateDerivative2 bxy = beta.multiply(exy);
946 final UnivariateDerivative2 u = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
947 final UnivariateDerivative2 v = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
948 final UnivariateDerivative2 x = parameters[0].multiply(u);
949 final UnivariateDerivative2 y = parameters[0].multiply(v);
950
951 // canonical orbit reference frame
952 final FieldVector3D<UnivariateDerivative2> p =
953 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
954 x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
955 y.multiply(sinI));
956
957 // dispatch derivatives
958 final Vector3D p0 = new Vector3D(p.getX().getValue(),
959 p.getY().getValue(),
960 p.getZ().getValue());
961 final Vector3D p1 = new Vector3D(p.getX().getFirstDerivative(),
962 p.getY().getFirstDerivative(),
963 p.getZ().getFirstDerivative());
964 final Vector3D p2 = new Vector3D(p.getX().getSecondDerivative(),
965 p.getY().getSecondDerivative(),
966 p.getZ().getSecondDerivative());
967 return new TimeStampedPVCoordinates(date, p0, p1, p2);
968
969 }
970
971 /** Computes the eccentric latitude argument from the mean latitude argument.
972 * @param alphaM = M + Ω mean latitude argument (rad)
973 * @param ex e cos(Ω), first component of circular eccentricity vector
974 * @param ey e sin(Ω), second component of circular eccentricity vector
975 * @return the eccentric latitude argument.
976 */
977 private UnivariateDerivative2 meanToEccentric(final UnivariateDerivative2 alphaM,
978 final UnivariateDerivative2 ex,
979 final UnivariateDerivative2 ey) {
980 // Generalization of Kepler equation to circular parameters
981 // with alphaE = PA + E and
982 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
983 UnivariateDerivative2 alphaE = alphaM;
984 UnivariateDerivative2 shift = alphaM.getField().getZero();
985 UnivariateDerivative2 alphaEMalphaM = alphaM.getField().getZero();
986 UnivariateDerivative2 cosAlphaE = alphaE.cos();
987 UnivariateDerivative2 sinAlphaE = alphaE.sin();
988 int iter = 0;
989 do {
990 final UnivariateDerivative2 f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
991 final UnivariateDerivative2 f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
992 final UnivariateDerivative2 f0 = alphaEMalphaM.subtract(f2);
993
994 final UnivariateDerivative2 f12 = f1.multiply(2);
995 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
996
997 alphaEMalphaM = alphaEMalphaM.subtract(shift);
998 alphaE = alphaM.add(alphaEMalphaM);
999 cosAlphaE = alphaE.cos();
1000 sinAlphaE = alphaE.sin();
1001
1002 } while (++iter < 50 && FastMath.abs(shift.getValue()) > 1.0e-12);
1003
1004 return alphaE;
1005
1006 }
1007
1008 /** {@inheritDoc} */
1009 protected double getMass(final AbsoluteDate date) {
1010 return models.get(date).mass;
1011 }
1012
1013 }