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