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