1 /* Copyright 2002-2024 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.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.FrameAlignedProvider;
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.PositionAngleType;
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 /** Default convergence threshold for mean parameters conversion. */
86 private static final double EPSILON_DEFAULT = 1.0e-11;
87
88 /** Default value for maxIterations. */
89 private static final int MAX_ITERATIONS_DEFAULT = 100;
90
91 /** Initial Eckstein-Hechler model. */
92 private EHModel initialModel;
93
94 /** All models. */
95 private transient TimeSpanMap<EHModel> models;
96
97 /** Reference radius of the central body attraction model (m). */
98 private double referenceRadius;
99
100 /** Central attraction coefficient (m³/s²). */
101 private double mu;
102
103 /** Un-normalized zonal coefficients. */
104 private double[] ck0;
105
106 /** Build a propagator from orbit and potential provider.
107 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
108 *
109 * <p>Using this constructor, an initial osculating orbit is considered.</p>
110 *
111 * @param initialOrbit initial orbit
112 * @param provider for un-normalized zonal coefficients
113 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider,
114 * UnnormalizedSphericalHarmonicsProvider)
115 * @see #EcksteinHechlerPropagator(Orbit, UnnormalizedSphericalHarmonicsProvider,
116 * PropagationType)
117 */
118 public EcksteinHechlerPropagator(final Orbit initialOrbit,
119 final UnnormalizedSphericalHarmonicsProvider provider) {
120 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
121 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
122 }
123
124 /**
125 * Private helper constructor.
126 * <p>Using this constructor, an initial osculating orbit is considered.</p>
127 * @param initialOrbit initial orbit
128 * @param attitude attitude provider
129 * @param mass spacecraft mass
130 * @param provider for un-normalized zonal coefficients
131 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
132 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
133 * UnnormalizedSphericalHarmonicsProvider,
134 * UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
135 * PropagationType)
136 */
137 public EcksteinHechlerPropagator(final Orbit initialOrbit,
138 final AttitudeProvider attitude,
139 final double mass,
140 final UnnormalizedSphericalHarmonicsProvider provider,
141 final UnnormalizedSphericalHarmonics harmonics) {
142 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
143 harmonics.getUnnormalizedCnm(2, 0),
144 harmonics.getUnnormalizedCnm(3, 0),
145 harmonics.getUnnormalizedCnm(4, 0),
146 harmonics.getUnnormalizedCnm(5, 0),
147 harmonics.getUnnormalizedCnm(6, 0));
148 }
149
150 /** Build a propagator from orbit and potential.
151 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
152 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
153 * are related to both the normalized coefficients
154 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
155 * and the J<sub>n</sub> one as follows:</p>
156 *
157 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
158 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
159 *
160 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
161 *
162 * <p>Using this constructor, an initial osculating orbit is considered.</p>
163 *
164 * @param initialOrbit initial orbit
165 * @param referenceRadius reference radius of the Earth for the potential model (m)
166 * @param mu central attraction coefficient (m³/s²)
167 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
168 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
169 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
170 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
171 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
172 * @see org.orekit.utils.Constants
173 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
174 * double, double, double, double, double)
175 */
176 public EcksteinHechlerPropagator(final Orbit initialOrbit,
177 final double referenceRadius, final double mu,
178 final double c20, final double c30, final double c40,
179 final double c50, final double c60) {
180 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
181 DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
182 }
183
184 /** Build a propagator from orbit, mass and potential provider.
185 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
186 *
187 * <p>Using this constructor, an initial osculating orbit is considered.</p>
188 *
189 * @param initialOrbit initial orbit
190 * @param mass spacecraft mass
191 * @param provider for un-normalized zonal coefficients
192 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
193 * UnnormalizedSphericalHarmonicsProvider)
194 */
195 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
196 final UnnormalizedSphericalHarmonicsProvider provider) {
197 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
198 mass, provider, provider.onDate(initialOrbit.getDate()));
199 }
200
201 /** Build a propagator from orbit, mass and potential.
202 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
203 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
204 * are related to both the normalized coefficients
205 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
206 * and the J<sub>n</sub> one as follows:</p>
207 *
208 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
209 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
210 *
211 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
212 *
213 * <p>Using this constructor, an initial osculating orbit is considered.</p>
214 *
215 * @param initialOrbit initial orbit
216 * @param mass spacecraft mass
217 * @param referenceRadius reference radius of the Earth for the potential model (m)
218 * @param mu central attraction coefficient (m³/s²)
219 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
220 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
221 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
222 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
223 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
224 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
225 * double, double, double, double, double)
226 */
227 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
228 final double referenceRadius, final double mu,
229 final double c20, final double c30, final double c40,
230 final double c50, final double c60) {
231 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
232 mass, referenceRadius, mu, c20, c30, c40, c50, c60);
233 }
234
235 /** Build a propagator from orbit, attitude provider and potential provider.
236 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
237 * <p>Using this constructor, an initial osculating orbit is considered.</p>
238 * @param initialOrbit initial orbit
239 * @param attitudeProv attitude provider
240 * @param provider for un-normalized zonal coefficients
241 */
242 public EcksteinHechlerPropagator(final Orbit initialOrbit,
243 final AttitudeProvider attitudeProv,
244 final UnnormalizedSphericalHarmonicsProvider provider) {
245 this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
246 }
247
248 /** Build a propagator from orbit, attitude provider and potential.
249 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
250 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
251 * are related to both the normalized coefficients
252 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
253 * and the J<sub>n</sub> one as follows:</p>
254 *
255 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
256 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
257 *
258 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
259 *
260 * <p>Using this constructor, an initial osculating orbit is considered.</p>
261 *
262 * @param initialOrbit initial orbit
263 * @param attitudeProv attitude provider
264 * @param referenceRadius reference radius of the Earth for the potential model (m)
265 * @param mu central attraction coefficient (m³/s²)
266 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
267 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
268 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
269 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
270 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
271 */
272 public EcksteinHechlerPropagator(final Orbit initialOrbit,
273 final AttitudeProvider attitudeProv,
274 final double referenceRadius, final double mu,
275 final double c20, final double c30, final double c40,
276 final double c50, final double c60) {
277 this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
278 }
279
280 /** Build a propagator from orbit, attitude provider, mass and potential provider.
281 * <p>Using this constructor, an initial osculating orbit is considered.</p>
282 * @param initialOrbit initial orbit
283 * @param attitudeProv attitude provider
284 * @param mass spacecraft mass
285 * @param provider for un-normalized zonal coefficients
286 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
287 * UnnormalizedSphericalHarmonicsProvider, PropagationType)
288 */
289 public EcksteinHechlerPropagator(final Orbit initialOrbit,
290 final AttitudeProvider attitudeProv,
291 final double mass,
292 final UnnormalizedSphericalHarmonicsProvider provider) {
293 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()));
294 }
295
296 /** Build a propagator from orbit, attitude provider, mass and potential.
297 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
298 * are related to both the normalized coefficients
299 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
300 * and the J<sub>n</sub> one as follows:</p>
301 *
302 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
303 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
304 *
305 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
306 *
307 * <p>Using this constructor, an initial osculating orbit is considered.</p>
308 *
309 * @param initialOrbit initial orbit
310 * @param attitudeProv attitude provider
311 * @param mass spacecraft mass
312 * @param referenceRadius reference radius of the Earth for the potential model (m)
313 * @param mu central attraction coefficient (m³/s²)
314 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
315 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
316 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
317 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
318 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
319 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
320 * double, double, double, double, double, PropagationType)
321 */
322 public EcksteinHechlerPropagator(final Orbit initialOrbit,
323 final AttitudeProvider attitudeProv,
324 final double mass,
325 final double referenceRadius, final double mu,
326 final double c20, final double c30, final double c40,
327 final double c50, final double c60) {
328 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, c60,
329 PropagationType.OSCULATING);
330 }
331
332
333 /** Build a propagator from orbit and potential provider.
334 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
335 *
336 * <p>Using this constructor, it is possible to define the initial orbit as
337 * a mean Eckstein-Hechler orbit or an osculating one.</p>
338 *
339 * @param initialOrbit initial orbit
340 * @param provider for un-normalized zonal coefficients
341 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
342 * @since 10.2
343 */
344 public EcksteinHechlerPropagator(final Orbit initialOrbit,
345 final UnnormalizedSphericalHarmonicsProvider provider,
346 final PropagationType initialType) {
347 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
348 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType);
349 }
350
351 /** Build a propagator from orbit, attitude provider, mass and potential provider.
352 * <p>Using this constructor, it is possible to define the initial orbit as
353 * a mean Eckstein-Hechler orbit or an osculating one.</p>
354 * @param initialOrbit initial orbit
355 * @param attitudeProv attitude provider
356 * @param mass spacecraft mass
357 * @param provider for un-normalized zonal coefficients
358 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
359 * @since 10.2
360 */
361 public EcksteinHechlerPropagator(final Orbit initialOrbit,
362 final AttitudeProvider attitudeProv,
363 final double mass,
364 final UnnormalizedSphericalHarmonicsProvider provider,
365 final PropagationType initialType) {
366 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType);
367 }
368
369 /**
370 * Private helper constructor.
371 * <p>Using this constructor, it is possible to define the initial orbit as
372 * a mean Eckstein-Hechler orbit or an osculating one.</p>
373 * @param initialOrbit initial orbit
374 * @param attitude attitude provider
375 * @param mass spacecraft mass
376 * @param provider for un-normalized zonal coefficients
377 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
378 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
379 * @since 10.2
380 */
381 public EcksteinHechlerPropagator(final Orbit initialOrbit,
382 final AttitudeProvider attitude,
383 final double mass,
384 final UnnormalizedSphericalHarmonicsProvider provider,
385 final UnnormalizedSphericalHarmonics harmonics,
386 final PropagationType initialType) {
387 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
388 harmonics.getUnnormalizedCnm(2, 0),
389 harmonics.getUnnormalizedCnm(3, 0),
390 harmonics.getUnnormalizedCnm(4, 0),
391 harmonics.getUnnormalizedCnm(5, 0),
392 harmonics.getUnnormalizedCnm(6, 0),
393 initialType);
394 }
395
396 /** Build a propagator from orbit, attitude provider, mass and potential.
397 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
398 * are related to both the normalized coefficients
399 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
400 * and the J<sub>n</sub> one as follows:</p>
401 *
402 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
403 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
404 *
405 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
406 *
407 * <p>Using this constructor, it is possible to define the initial orbit as
408 * a mean Eckstein-Hechler orbit or an osculating one.</p>
409 *
410 * @param initialOrbit initial orbit
411 * @param attitudeProv attitude provider
412 * @param mass spacecraft mass
413 * @param referenceRadius reference radius of the Earth for the potential model (m)
414 * @param mu central attraction coefficient (m³/s²)
415 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
416 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
417 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
418 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
419 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
420 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
421 * @since 10.2
422 */
423 public EcksteinHechlerPropagator(final Orbit initialOrbit,
424 final AttitudeProvider attitudeProv,
425 final double mass,
426 final double referenceRadius, final double mu,
427 final double c20, final double c30, final double c40,
428 final double c50, final double c60,
429 final PropagationType initialType) {
430 this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
431 c20, c30, c40, c50, c60, initialType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
432 }
433
434 /** Build a propagator from orbit, attitude provider, mass and potential.
435 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
436 * are related to both the normalized coefficients
437 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
438 * and the J<sub>n</sub> one as follows:</p>
439 *
440 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
441 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
442 *
443 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
444 *
445 * <p>Using this constructor, it is possible to define the initial orbit as
446 * a mean Eckstein-Hechler orbit or an osculating one.</p>
447 *
448 * @param initialOrbit initial orbit
449 * @param attitudeProv attitude provider
450 * @param mass spacecraft mass
451 * @param referenceRadius reference radius of the Earth for the potential model (m)
452 * @param mu central attraction coefficient (m³/s²)
453 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
454 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
455 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
456 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
457 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
458 * @param epsilon convergence threshold for mean parameters conversion
459 * @param maxIterations maximum iterations for mean parameters conversion
460 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
461 * @since 11.2
462 */
463 public EcksteinHechlerPropagator(final Orbit initialOrbit,
464 final AttitudeProvider attitudeProv,
465 final double mass,
466 final double referenceRadius, final double mu,
467 final double c20, final double c30, final double c40,
468 final double c50, final double c60,
469 final PropagationType initialType,
470 final double epsilon, final int maxIterations) {
471 super(attitudeProv);
472
473 // store model coefficients
474 this.referenceRadius = referenceRadius;
475 this.mu = mu;
476 this.ck0 = new double[] {
477 0.0, 0.0, c20, c30, c40, c50, c60
478 };
479
480 // compute mean parameters if needed
481 // transform into circular adapted parameters used by the Eckstein-Hechler model
482 resetInitialState(new SpacecraftState(initialOrbit,
483 attitudeProv.getAttitude(initialOrbit,
484 initialOrbit.getDate(),
485 initialOrbit.getFrame()),
486 mass),
487 initialType, epsilon, maxIterations);
488
489 }
490
491 /** Conversion from osculating to mean orbit.
492 * <p>
493 * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
494 * osculating SpacecraftState in input.
495 * </p>
496 * <p>
497 * Since the osculating orbit is obtained with the computation of
498 * short-periodic variation, the resulting output will depend on
499 * the gravity field parameterized in input.
500 * </p>
501 * <p>
502 * The computation is done through a fixed-point iteration process.
503 * </p>
504 * @param osculating osculating orbit to convert
505 * @param provider for un-normalized zonal coefficients
506 * @param harmonics {@code provider.onDate(osculating.getDate())}
507 * @return mean orbit in a Eckstein-Hechler sense
508 * @since 11.2
509 */
510 public static CircularOrbit computeMeanOrbit(final Orbit osculating,
511 final UnnormalizedSphericalHarmonicsProvider provider,
512 final UnnormalizedSphericalHarmonics harmonics) {
513 return computeMeanOrbit(osculating, provider, harmonics, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
514 }
515
516 /** Conversion from osculating to mean orbit.
517 * <p>
518 * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
519 * osculating SpacecraftState in input.
520 * </p>
521 * <p>
522 * Since the osculating orbit is obtained with the computation of
523 * short-periodic variation, the resulting output will depend on
524 * the gravity field parameterized in input.
525 * </p>
526 * <p>
527 * The computation is done through a fixed-point iteration process.
528 * </p>
529 * @param osculating osculating orbit to convert
530 * @param provider for un-normalized zonal coefficients
531 * @param harmonics {@code provider.onDate(osculating.getDate())}
532 * @param epsilon convergence threshold for mean parameters conversion
533 * @param maxIterations maximum iterations for mean parameters conversion
534 * @return mean orbit in a Eckstein-Hechler sense
535 * @since 11.2
536 */
537 public static CircularOrbit computeMeanOrbit(final Orbit osculating,
538 final UnnormalizedSphericalHarmonicsProvider provider,
539 final UnnormalizedSphericalHarmonics harmonics,
540 final double epsilon, final int maxIterations) {
541 return computeMeanOrbit(osculating,
542 provider.getAe(), provider.getMu(),
543 harmonics.getUnnormalizedCnm(2, 0),
544 harmonics.getUnnormalizedCnm(3, 0),
545 harmonics.getUnnormalizedCnm(4, 0),
546 harmonics.getUnnormalizedCnm(5, 0),
547 harmonics.getUnnormalizedCnm(6, 0),
548 epsilon, maxIterations);
549 }
550
551 /** Conversion from osculating to mean orbit.
552 * <p>
553 * Compute mean orbit <b>in a Eckstein-Hechler sense</b>, corresponding to the
554 * osculating SpacecraftState in input.
555 * </p>
556 * <p>
557 * Since the osculating orbit is obtained with the computation of
558 * short-periodic variation, the resulting output will depend on
559 * the gravity field parameterized in input.
560 * </p>
561 * <p>
562 * The computation is done through a fixed-point iteration process.
563 * </p>
564 * @param osculating osculating orbit to convert
565 * @param referenceRadius reference radius of the Earth for the potential model (m)
566 * @param mu central attraction coefficient (m³/s²)
567 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
568 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
569 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
570 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
571 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
572 * @param epsilon convergence threshold for mean parameters conversion
573 * @param maxIterations maximum iterations for mean parameters conversion
574 * @return mean orbit in a Eckstein-Hechler sense
575 * @since 11.2
576 */
577 public static CircularOrbit computeMeanOrbit(final Orbit osculating,
578 final double referenceRadius, final double mu,
579 final double c20, final double c30, final double c40,
580 final double c50, final double c60,
581 final double epsilon, final int maxIterations) {
582 final EcksteinHechlerPropagator propagator =
583 new EcksteinHechlerPropagator(osculating,
584 FrameAlignedProvider.of(osculating.getFrame()),
585 DEFAULT_MASS,
586 referenceRadius, mu, c20, c30, c40, c50, c60,
587 PropagationType.OSCULATING, epsilon, maxIterations);
588 return propagator.initialModel.mean;
589 }
590
591 /** {@inheritDoc}
592 * <p>The new initial state to consider
593 * must be defined with an osculating orbit.</p>
594 * @see #resetInitialState(SpacecraftState, PropagationType)
595 */
596 public void resetInitialState(final SpacecraftState state) {
597 resetInitialState(state, PropagationType.OSCULATING);
598 }
599
600 /** Reset the propagator initial state.
601 * @param state new initial state to consider
602 * @param stateType mean Eckstein-Hechler orbit or osculating orbit
603 * @since 10.2
604 */
605 public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
606 resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
607 }
608
609 /** Reset the propagator initial state.
610 * @param state new initial state to consider
611 * @param stateType mean Eckstein-Hechler orbit or osculating orbit
612 * @param epsilon convergence threshold for mean parameters conversion
613 * @param maxIterations maximum iterations for mean parameters conversion
614 * @since 11.2
615 */
616 public void resetInitialState(final SpacecraftState state, final PropagationType stateType,
617 final double epsilon, final int maxIterations) {
618 super.resetInitialState(state);
619 final CircularOrbit circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
620 this.initialModel = (stateType == PropagationType.MEAN) ?
621 new EHModel(circular, state.getMass(), referenceRadius, mu, ck0) :
622 computeMeanParameters(circular, state.getMass(), epsilon, maxIterations);
623 this.models = new TimeSpanMap<EHModel>(initialModel);
624 }
625
626 /** {@inheritDoc} */
627 protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
628 resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
629 }
630
631 /** Reset an intermediate state.
632 * @param state new intermediate state to consider
633 * @param forward if true, the intermediate state is valid for
634 * propagations after itself
635 * @param epsilon convergence threshold for mean parameters conversion
636 * @param maxIterations maximum iterations for mean parameters conversion
637 * @since 11.2
638 */
639 protected void resetIntermediateState(final SpacecraftState state, final boolean forward,
640 final double epsilon, final int maxIterations) {
641 final EHModel newModel = computeMeanParameters((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit()),
642 state.getMass(), epsilon, maxIterations);
643 if (forward) {
644 models.addValidAfter(newModel, state.getDate(), false);
645 } else {
646 models.addValidBefore(newModel, state.getDate(), false);
647 }
648 stateChanged(state);
649 }
650
651
652 /** Compute mean parameters according to the Eckstein-Hechler analytical model.
653 * @param osculating osculating orbit
654 * @param mass constant mass
655 * @param epsilon convergence threshold for mean parameters conversion
656 * @param maxIterations maximum iterations for mean parameters conversion
657 * @return Eckstein-Hechler mean model
658 */
659 private EHModel computeMeanParameters(final CircularOrbit osculating, final double mass,
660 final double epsilon, final int maxIterations) {
661
662 // sanity check
663 if (osculating.getA() < referenceRadius) {
664 throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
665 osculating.getA());
666 }
667
668 // rough initialization of the mean parameters
669 EHModel current = new EHModel(osculating, mass, referenceRadius, mu, ck0);
670
671 // threshold for each parameter
672 final double thresholdA = epsilon * (1 + FastMath.abs(current.mean.getA()));
673 final double thresholdE = epsilon * (1 + current.mean.getE());
674 final double thresholdAngles = epsilon * MathUtils.TWO_PI;
675
676 int i = 0;
677 while (i++ < maxIterations) {
678
679 // recompute the osculating parameters from the current mean parameters
680 final UnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate());
681
682 // adapted parameters residuals
683 final double deltaA = osculating.getA() - parameters[0].getValue();
684 final double deltaEx = osculating.getCircularEx() - parameters[1].getValue();
685 final double deltaEy = osculating.getCircularEy() - parameters[2].getValue();
686 final double deltaI = osculating.getI() - parameters[3].getValue();
687 final double deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
688 parameters[4].getValue(),
689 0.0);
690 final double deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM() - parameters[5].getValue(), 0.0);
691
692 // update mean parameters
693 current = new EHModel(new CircularOrbit(current.mean.getA() + deltaA,
694 current.mean.getCircularEx() + deltaEx,
695 current.mean.getCircularEy() + deltaEy,
696 current.mean.getI() + deltaI,
697 current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
698 current.mean.getAlphaM() + deltaAlphaM,
699 PositionAngleType.MEAN,
700 current.mean.getFrame(),
701 current.mean.getDate(), mu),
702 mass, referenceRadius, mu, ck0);
703
704 // check convergence
705 if (FastMath.abs(deltaA) < thresholdA &&
706 FastMath.abs(deltaEx) < thresholdE &&
707 FastMath.abs(deltaEy) < thresholdE &&
708 FastMath.abs(deltaI) < thresholdAngles &&
709 FastMath.abs(deltaRAAN) < thresholdAngles &&
710 FastMath.abs(deltaAlphaM) < thresholdAngles) {
711 return current;
712 }
713
714 }
715
716 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECKSTEIN_HECHLER_MEAN_PARAMETERS, i);
717
718 }
719
720 /** {@inheritDoc} */
721 public CartesianOrbit propagateOrbit(final AbsoluteDate date) {
722 // compute Cartesian parameters, taking derivatives into account
723 // to make sure velocity and acceleration are consistent
724 final EHModel current = models.get(date);
725 return new CartesianOrbit(toCartesian(date, current.propagateParameters(date)),
726 current.mean.getFrame(), mu);
727 }
728
729 /**
730 * Get the central attraction coefficient μ.
731 * @return mu central attraction coefficient (m³/s²)
732 * @since 11.1
733 */
734 public double getMu() {
735 return mu;
736 }
737
738 /**
739 * Get the un-normalized zonal coefficients.
740 * @return the un-normalized zonal coefficients
741 * @since 11.1
742 */
743 public double[] getCk0() {
744 return ck0.clone();
745 }
746
747 /**
748 * Get the reference radius of the central body attraction model.
749 * @return the reference radius in meters
750 * @since 11.1
751 */
752 public double getReferenceRadius() {
753 return referenceRadius;
754 }
755
756 /** {@inheritDoc} */
757 @Override
758 protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
759 final DoubleArrayDictionary initialJacobianColumns) {
760 // Create the harvester
761 final EcksteinHechlerHarvester harvester = new EcksteinHechlerHarvester(this, stmName, initialStm, initialJacobianColumns);
762 // Update the list of additional state provider
763 addAdditionalStateProvider(harvester);
764 // Return the configured harvester
765 return harvester;
766 }
767
768 /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
769 private static class EHModel implements Serializable {
770
771 /** Serializable UID. */
772 private static final long serialVersionUID = 20160115L;
773
774 /** Mean orbit. */
775 private final CircularOrbit mean;
776
777 /** Constant mass. */
778 private final double mass;
779
780 // CHECKSTYLE: stop JavadocVariable check
781
782 // preprocessed values
783 private final double xnotDot;
784 private final double rdpom;
785 private final double rdpomp;
786 private final double eps1;
787 private final double eps2;
788 private final double xim;
789 private final double ommD;
790 private final double rdl;
791 private final double aMD;
792
793 private final double kh;
794 private final double kl;
795
796 private final double ax1;
797 private final double ay1;
798 private final double as1;
799 private final double ac2;
800 private final double axy3;
801 private final double as3;
802 private final double ac4;
803 private final double as5;
804 private final double ac6;
805
806 private final double ex1;
807 private final double exx2;
808 private final double exy2;
809 private final double ex3;
810 private final double ex4;
811
812 private final double ey1;
813 private final double eyx2;
814 private final double eyy2;
815 private final double ey3;
816 private final double ey4;
817
818 private final double rx1;
819 private final double ry1;
820 private final double r2;
821 private final double r3;
822 private final double rl;
823
824 private final double iy1;
825 private final double ix1;
826 private final double i2;
827 private final double i3;
828 private final double ih;
829
830 private final double lx1;
831 private final double ly1;
832 private final double l2;
833 private final double l3;
834 private final double ll;
835
836 // CHECKSTYLE: resume JavadocVariable check
837
838 /** Create a model for specified mean orbit.
839 * @param mean mean orbit
840 * @param mass constant mass
841 * @param referenceRadius reference radius of the central body attraction model (m)
842 * @param mu central attraction coefficient (m³/s²)
843 * @param ck0 un-normalized zonal coefficients
844 */
845 EHModel(final CircularOrbit mean, final double mass,
846 final double referenceRadius, final double mu, final double[] ck0) {
847
848 this.mean = mean;
849 this.mass = mass;
850
851 // preliminary processing
852 double q = referenceRadius / mean.getA();
853 double ql = q * q;
854 final double g2 = ck0[2] * ql;
855 ql *= q;
856 final double g3 = ck0[3] * ql;
857 ql *= q;
858 final double g4 = ck0[4] * ql;
859 ql *= q;
860 final double g5 = ck0[5] * ql;
861 ql *= q;
862 final double g6 = ck0[6] * ql;
863
864 final SinCos sc = FastMath.sinCos(mean.getI());
865 final double cosI1 = sc.cos();
866 final double sinI1 = sc.sin();
867 final double sinI2 = sinI1 * sinI1;
868 final double sinI4 = sinI2 * sinI2;
869 final double sinI6 = sinI2 * sinI4;
870
871 if (sinI2 < 1.0e-10) {
872 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
873 FastMath.toDegrees(mean.getI()));
874 }
875
876 if (FastMath.abs(sinI2 - 4.0 / 5.0) < 1.0e-3) {
877 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
878 FastMath.toDegrees(mean.getI()));
879 }
880
881 if (mean.getE() > 0.1) {
882 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
883 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
884 mean.getE());
885 }
886
887 xnotDot = FastMath.sqrt(mu / mean.getA()) / mean.getA();
888
889 rdpom = -0.75 * g2 * (4.0 - 5.0 * sinI2);
890 rdpomp = 7.5 * g4 * (1.0 - 31.0 / 8.0 * sinI2 + 49.0 / 16.0 * sinI4) -
891 13.125 * g6 * (1.0 - 8.0 * sinI2 + 129.0 / 8.0 * sinI4 - 297.0 / 32.0 * sinI6);
892
893 q = 3.0 / (32.0 * rdpom);
894 eps1 = q * g4 * sinI2 * (30.0 - 35.0 * sinI2) -
895 175.0 * q * g6 * sinI2 * (1.0 - 3.0 * sinI2 + 2.0625 * sinI4);
896 q = 3.0 * sinI1 / (8.0 * rdpom);
897 eps2 = q * g3 * (4.0 - 5.0 * sinI2) - q * g5 * (10.0 - 35.0 * sinI2 + 26.25 * sinI4);
898
899 xim = mean.getI();
900 ommD = cosI1 * (1.50 * g2 - 2.25 * g2 * g2 * (2.5 - 19.0 / 6.0 * sinI2) +
901 0.9375 * g4 * (7.0 * sinI2 - 4.0) +
902 3.28125 * g6 * (2.0 - 9.0 * sinI2 + 8.25 * sinI4));
903
904 rdl = 1.0 - 1.50 * g2 * (3.0 - 4.0 * sinI2);
905 aMD = rdl +
906 2.25 * g2 * g2 * (9.0 - 263.0 / 12.0 * sinI2 + 341.0 / 24.0 * sinI4) +
907 15.0 / 16.0 * g4 * (8.0 - 31.0 * sinI2 + 24.5 * sinI4) +
908 105.0 / 32.0 * g6 * (-10.0 / 3.0 + 25.0 * sinI2 - 48.75 * sinI4 + 27.5 * sinI6);
909
910 final double qq = -1.5 * g2 / rdl;
911 final double qA = 0.75 * g2 * g2 * sinI2;
912 final double qB = 0.25 * g4 * sinI2;
913 final double qC = 105.0 / 16.0 * g6 * sinI2;
914 final double qD = -0.75 * g3 * sinI1;
915 final double qE = 3.75 * g5 * sinI1;
916 kh = 0.375 / rdpom;
917 kl = kh / sinI1;
918
919 ax1 = qq * (2.0 - 3.5 * sinI2);
920 ay1 = qq * (2.0 - 2.5 * sinI2);
921 as1 = qD * (4.0 - 5.0 * sinI2) +
922 qE * (2.625 * sinI4 - 3.5 * sinI2 + 1.0);
923 ac2 = qq * sinI2 +
924 qA * 7.0 * (2.0 - 3.0 * sinI2) +
925 qB * (15.0 - 17.5 * sinI2) +
926 qC * (3.0 * sinI2 - 1.0 - 33.0 / 16.0 * sinI4);
927 axy3 = qq * 3.5 * sinI2;
928 as3 = qD * 5.0 / 3.0 * sinI2 +
929 qE * 7.0 / 6.0 * sinI2 * (1.0 - 1.125 * sinI2);
930 ac4 = qA * sinI2 +
931 qB * 4.375 * sinI2 +
932 qC * 0.75 * (1.1 * sinI4 - sinI2);
933
934 as5 = qE * 21.0 / 80.0 * sinI4;
935
936 ac6 = qC * -11.0 / 80.0 * sinI4;
937
938 ex1 = qq * (1.0 - 1.25 * sinI2);
939 exx2 = qq * 0.5 * (3.0 - 5.0 * sinI2);
940 exy2 = qq * (2.0 - 1.5 * sinI2);
941 ex3 = qq * 7.0 / 12.0 * sinI2;
942 ex4 = qq * 17.0 / 8.0 * sinI2;
943
944 ey1 = qq * (1.0 - 1.75 * sinI2);
945 eyx2 = qq * (1.0 - 3.0 * sinI2);
946 eyy2 = qq * (2.0 * sinI2 - 1.5);
947 ey3 = qq * 7.0 / 12.0 * sinI2;
948 ey4 = qq * 17.0 / 8.0 * sinI2;
949
950 q = -qq * cosI1;
951 rx1 = 3.5 * q;
952 ry1 = -2.5 * q;
953 r2 = -0.5 * q;
954 r3 = 7.0 / 6.0 * q;
955 rl = g3 * cosI1 * (4.0 - 15.0 * sinI2) -
956 2.5 * g5 * cosI1 * (4.0 - 42.0 * sinI2 + 52.5 * sinI4);
957
958 q = 0.5 * qq * sinI1 * cosI1;
959 iy1 = q;
960 ix1 = -q;
961 i2 = q;
962 i3 = q * 7.0 / 3.0;
963 ih = -g3 * cosI1 * (4.0 - 5.0 * sinI2) +
964 2.5 * g5 * cosI1 * (4.0 - 14.0 * sinI2 + 10.5 * sinI4);
965
966 lx1 = qq * (7.0 - 77.0 / 8.0 * sinI2);
967 ly1 = qq * (55.0 / 8.0 * sinI2 - 7.50);
968 l2 = qq * (1.25 * sinI2 - 0.5);
969 l3 = qq * (77.0 / 24.0 * sinI2 - 7.0 / 6.0);
970 ll = g3 * (53.0 * sinI2 - 4.0 - 57.5 * sinI4) +
971 2.5 * g5 * (4.0 - 96.0 * sinI2 + 269.5 * sinI4 - 183.75 * sinI6);
972
973 }
974
975 /** Extrapolate an orbit up to a specific target date.
976 * @param date target date for the orbit
977 * @return propagated parameters
978 */
979 public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) {
980
981 // Keplerian evolution
982 final UnivariateDerivative2 dt = new UnivariateDerivative2(date.durationFrom(mean.getDate()), 1.0, 0.0);
983 final UnivariateDerivative2 xnot = dt.multiply(xnotDot);
984
985 // secular effects
986
987 // eccentricity
988 final UnivariateDerivative2 x = xnot.multiply(rdpom + rdpomp);
989 final UnivariateDerivative2 cx = x.cos();
990 final UnivariateDerivative2 sx = x.sin();
991 final UnivariateDerivative2 exm = cx.multiply(mean.getCircularEx()).
992 add(sx.multiply(eps2 - (1.0 - eps1) * mean.getCircularEy()));
993 final UnivariateDerivative2 eym = sx.multiply((1.0 + eps1) * mean.getCircularEx()).
994 add(cx.multiply(mean.getCircularEy() - eps2)).
995 add(eps2);
996
997 // no secular effect on inclination
998
999 // right ascension of ascending node
1000 final UnivariateDerivative2 omm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ommD * xnot.getValue(),
1001 FastMath.PI),
1002 ommD * xnotDot,
1003 0.0);
1004
1005 // latitude argument
1006 final UnivariateDerivative2 xlm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getAlphaM() + aMD * xnot.getValue(),
1007 FastMath.PI),
1008 aMD * xnotDot,
1009 0.0);
1010
1011 // periodical terms
1012 final UnivariateDerivative2 cl1 = xlm.cos();
1013 final UnivariateDerivative2 sl1 = xlm.sin();
1014 final UnivariateDerivative2 cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
1015 final UnivariateDerivative2 sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
1016 final UnivariateDerivative2 cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
1017 final UnivariateDerivative2 sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
1018 final UnivariateDerivative2 cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
1019 final UnivariateDerivative2 sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
1020 final UnivariateDerivative2 cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
1021 final UnivariateDerivative2 sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
1022 final UnivariateDerivative2 cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));
1023
1024 final UnivariateDerivative2 qh = eym.subtract(eps2).multiply(kh);
1025 final UnivariateDerivative2 ql = exm.multiply(kl);
1026
1027 final UnivariateDerivative2 exmCl1 = exm.multiply(cl1);
1028 final UnivariateDerivative2 exmSl1 = exm.multiply(sl1);
1029 final UnivariateDerivative2 eymCl1 = eym.multiply(cl1);
1030 final UnivariateDerivative2 eymSl1 = eym.multiply(sl1);
1031 final UnivariateDerivative2 exmCl2 = exm.multiply(cl2);
1032 final UnivariateDerivative2 exmSl2 = exm.multiply(sl2);
1033 final UnivariateDerivative2 eymCl2 = eym.multiply(cl2);
1034 final UnivariateDerivative2 eymSl2 = eym.multiply(sl2);
1035 final UnivariateDerivative2 exmCl3 = exm.multiply(cl3);
1036 final UnivariateDerivative2 exmSl3 = exm.multiply(sl3);
1037 final UnivariateDerivative2 eymCl3 = eym.multiply(cl3);
1038 final UnivariateDerivative2 eymSl3 = eym.multiply(sl3);
1039 final UnivariateDerivative2 exmCl4 = exm.multiply(cl4);
1040 final UnivariateDerivative2 exmSl4 = exm.multiply(sl4);
1041 final UnivariateDerivative2 eymCl4 = eym.multiply(cl4);
1042 final UnivariateDerivative2 eymSl4 = eym.multiply(sl4);
1043
1044 // semi major axis
1045 final UnivariateDerivative2 rda = exmCl1.multiply(ax1).
1046 add(eymSl1.multiply(ay1)).
1047 add(sl1.multiply(as1)).
1048 add(cl2.multiply(ac2)).
1049 add(exmCl3.add(eymSl3).multiply(axy3)).
1050 add(sl3.multiply(as3)).
1051 add(cl4.multiply(ac4)).
1052 add(sl5.multiply(as5)).
1053 add(cl6.multiply(ac6));
1054
1055 // eccentricity
1056 final UnivariateDerivative2 rdex = cl1.multiply(ex1).
1057 add(exmCl2.multiply(exx2)).
1058 add(eymSl2.multiply(exy2)).
1059 add(cl3.multiply(ex3)).
1060 add(exmCl4.add(eymSl4).multiply(ex4));
1061 final UnivariateDerivative2 rdey = sl1.multiply(ey1).
1062 add(exmSl2.multiply(eyx2)).
1063 add(eymCl2.multiply(eyy2)).
1064 add(sl3.multiply(ey3)).
1065 add(exmSl4.subtract(eymCl4).multiply(ey4));
1066
1067 // ascending node
1068 final UnivariateDerivative2 rdom = exmSl1.multiply(rx1).
1069 add(eymCl1.multiply(ry1)).
1070 add(sl2.multiply(r2)).
1071 add(eymCl3.subtract(exmSl3).multiply(r3)).
1072 add(ql.multiply(rl));
1073
1074 // inclination
1075 final UnivariateDerivative2 rdxi = eymSl1.multiply(iy1).
1076 add(exmCl1.multiply(ix1)).
1077 add(cl2.multiply(i2)).
1078 add(exmCl3.add(eymSl3).multiply(i3)).
1079 add(qh.multiply(ih));
1080
1081 // latitude argument
1082 final UnivariateDerivative2 rdxl = exmSl1.multiply(lx1).
1083 add(eymCl1.multiply(ly1)).
1084 add(sl2.multiply(l2)).
1085 add(exmSl3.subtract(eymCl3).multiply(l3)).
1086 add(ql.multiply(ll));
1087
1088 // osculating parameters
1089 return new UnivariateDerivative2[] {
1090 rda.add(1.0).multiply(mean.getA()),
1091 rdex.add(exm),
1092 rdey.add(eym),
1093 rdxi.add(xim),
1094 rdom.add(omm),
1095 rdxl.add(xlm)
1096 };
1097
1098 }
1099
1100 }
1101
1102 /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
1103 * @param date date of the orbital parameters
1104 * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
1105 * @return Cartesian coordinates consistent with values and derivatives
1106 */
1107 private TimeStampedPVCoordinates toCartesian(final AbsoluteDate date, final UnivariateDerivative2[] parameters) {
1108
1109 // evaluate coordinates in the orbit canonical reference frame
1110 final UnivariateDerivative2 cosOmega = parameters[4].cos();
1111 final UnivariateDerivative2 sinOmega = parameters[4].sin();
1112 final UnivariateDerivative2 cosI = parameters[3].cos();
1113 final UnivariateDerivative2 sinI = parameters[3].sin();
1114 final UnivariateDerivative2 alphaE = meanToEccentric(parameters[5], parameters[1], parameters[2]);
1115 final UnivariateDerivative2 cosAE = alphaE.cos();
1116 final UnivariateDerivative2 sinAE = alphaE.sin();
1117 final UnivariateDerivative2 ex2 = parameters[1].square();
1118 final UnivariateDerivative2 ey2 = parameters[2].square();
1119 final UnivariateDerivative2 exy = parameters[1].multiply(parameters[2]);
1120 final UnivariateDerivative2 q = ex2.add(ey2).subtract(1).negate().sqrt();
1121 final UnivariateDerivative2 beta = q.add(1).reciprocal();
1122 final UnivariateDerivative2 bx2 = beta.multiply(ex2);
1123 final UnivariateDerivative2 by2 = beta.multiply(ey2);
1124 final UnivariateDerivative2 bxy = beta.multiply(exy);
1125 final UnivariateDerivative2 u = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
1126 final UnivariateDerivative2 v = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
1127 final UnivariateDerivative2 x = parameters[0].multiply(u);
1128 final UnivariateDerivative2 y = parameters[0].multiply(v);
1129
1130 // canonical orbit reference frame
1131 final FieldVector3D<UnivariateDerivative2> p =
1132 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
1133 x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
1134 y.multiply(sinI));
1135
1136 // dispatch derivatives
1137 final Vector3D p0 = new Vector3D(p.getX().getValue(),
1138 p.getY().getValue(),
1139 p.getZ().getValue());
1140 final Vector3D p1 = new Vector3D(p.getX().getFirstDerivative(),
1141 p.getY().getFirstDerivative(),
1142 p.getZ().getFirstDerivative());
1143 final Vector3D p2 = new Vector3D(p.getX().getSecondDerivative(),
1144 p.getY().getSecondDerivative(),
1145 p.getZ().getSecondDerivative());
1146 return new TimeStampedPVCoordinates(date, p0, p1, p2);
1147
1148 }
1149
1150 /** Computes the eccentric latitude argument from the mean latitude argument.
1151 * @param alphaM = M + Ω mean latitude argument (rad)
1152 * @param ex e cos(Ω), first component of circular eccentricity vector
1153 * @param ey e sin(Ω), second component of circular eccentricity vector
1154 * @return the eccentric latitude argument.
1155 */
1156 private UnivariateDerivative2 meanToEccentric(final UnivariateDerivative2 alphaM,
1157 final UnivariateDerivative2 ex,
1158 final UnivariateDerivative2 ey) {
1159 // Generalization of Kepler equation to circular parameters
1160 // with alphaE = PA + E and
1161 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
1162 UnivariateDerivative2 alphaE = alphaM;
1163 UnivariateDerivative2 shift = alphaM.getField().getZero();
1164 UnivariateDerivative2 alphaEMalphaM = alphaM.getField().getZero();
1165 UnivariateDerivative2 cosAlphaE = alphaE.cos();
1166 UnivariateDerivative2 sinAlphaE = alphaE.sin();
1167 int iter = 0;
1168 do {
1169 final UnivariateDerivative2 f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
1170 final UnivariateDerivative2 f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
1171 final UnivariateDerivative2 f0 = alphaEMalphaM.subtract(f2);
1172
1173 final UnivariateDerivative2 f12 = f1.multiply(2);
1174 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
1175
1176 alphaEMalphaM = alphaEMalphaM.subtract(shift);
1177 alphaE = alphaM.add(alphaEMalphaM);
1178 cosAlphaE = alphaE.cos();
1179 sinAlphaE = alphaE.sin();
1180
1181 } while (++iter < 50 && FastMath.abs(shift.getValue()) > 1.0e-12);
1182
1183 return alphaE;
1184
1185 }
1186
1187 /** {@inheritDoc} */
1188 protected double getMass(final AbsoluteDate date) {
1189 return models.get(date).mass;
1190 }
1191
1192 }