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