1 /* Copyright 2002-2026 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.FieldUnivariateDerivative1;
25 import org.hipparchus.util.CombinatoricsUtils;
26 import org.hipparchus.util.FastMath;
27 import org.hipparchus.util.FieldSinCos;
28 import org.hipparchus.util.MathUtils;
29 import org.orekit.attitudes.AttitudeProvider;
30 import org.orekit.attitudes.FrameAlignedProvider;
31 import org.orekit.errors.OrekitException;
32 import org.orekit.errors.OrekitMessages;
33 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
34 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
35 import org.orekit.orbits.FieldKeplerianAnomalyUtility;
36 import org.orekit.orbits.FieldKeplerianOrbit;
37 import org.orekit.orbits.FieldOrbit;
38 import org.orekit.orbits.OrbitParamsType;
39 import org.orekit.orbits.PositionAngleType;
40 import org.orekit.propagation.FieldSpacecraftState;
41 import org.orekit.propagation.PropagationType;
42 import org.orekit.propagation.analytical.tle.FieldTLE;
43 import org.orekit.propagation.conversion.osc2mean.BrouwerLyddaneTheory;
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.time.TimeInterval;
49 import org.orekit.utils.FieldTimeSpanMap;
50 import org.orekit.utils.drivers.ParameterDriver;
51
52 /** This class propagates a {@link org.orekit.propagation.FieldSpacecraftState}
53 * using the analytical Brouwer-Lyddane model (from J2 to J5 zonal harmonics).
54 * <p>
55 * At the opposite of the {@link FieldEcksteinHechlerPropagator}, the Brouwer-Lyddane model is
56 * suited for elliptical orbits, there is no problem having a rather small eccentricity or inclination
57 * (Lyddane helped to solve this issue with the Brouwer model). Singularity for the critical
58 * inclination i = 63.4° is avoided using the method developed in Warren Phipps' 1992 thesis.
59 * <p>
60 * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
61 * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
62 * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
63 * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient
64 * {@link #M2Driver}.
65 * <p>
66 * Usually, M2 is adjusted during an orbit determination process and it represents the
67 * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
68 * The behavior of M2 is close to the {@link FieldTLE#getBStar()} parameter for the TLE.
69 * </p>
70 * <p>
71 * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track secular
72 * effects are not considered in the dynamical model. Typical values for M2 are not known.
73 * It depends on the orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15).
74 * The unit of M2 is rad/s².
75 * </p>
76 * <p>
77 * The along-track effects, represented by the secular rates of the mean semi-major axis
78 * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis.
79 * </p>
80 * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag.
81 * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959."
82 * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the
83 * artificial satellite. The Astronomical Journal 68 (1963): 555."
84 * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center
85 * (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992."
86 * @see "Solomon, Daniel, THE NAVSPASUR Satellite Motion Model,
87 * Naval Research Laboratory, August 8, 1991."
88 * @author Melina Vanel
89 * @author Bryan Cazabonne
90 * @author Pascal Parraud
91 * @since 11.1
92 * @param <T> type of the field elements
93 */
94 public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> extends FieldAbstractAnalyticalPropagator<T> {
95
96 /** Parameters scaling factor.
97 * <p>
98 * We use a power of 2 to avoid numeric noise introduction
99 * in the multiplications/divisions sequences.
100 * </p>
101 */
102 private static final double SCALE = FastMath.scalb(1.0, -32);
103
104 /** Beta constant used by T2 function. */
105 private static final double BETA = FastMath.scalb(100, -11);
106
107 /** Max value for the eccentricity. */
108 private static final double MAX_ECC = 0.999999;
109
110 /** Initial Brouwer-Lyddane model. */
111 private FieldBLModel<T> initialModel;
112
113 /** All models. */
114 private FieldTimeSpanMap<FieldBLModel<T>, T> models;
115
116 /** Reference radius of the central body attraction model (m). */
117 private final double referenceRadius;
118
119 /** Central attraction coefficient (m³/s²). */
120 private final T mu;
121
122 /** Un-normalized zonal coefficients. */
123 private final double[] ck0;
124
125 /** Empirical coefficient used in the drag modeling. */
126 private final ParameterDriver M2Driver;
127
128 /** Build a propagator from orbit and potential provider.
129 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
130 *
131 * <p>Using this constructor, an initial osculating orbit is considered.</p>
132 *
133 * @param initialOrbit initial orbit
134 * @param provider for un-normalized zonal coefficients
135 * @param m2Value value of empirical drag coefficient in rad/s².
136 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
137 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
138 */
139 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
140 final UnnormalizedSphericalHarmonicsProvider provider,
141 final double m2Value) {
142 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
143 initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
144 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), m2Value);
145 }
146
147 /**
148 * Private helper constructor.
149 * <p>Using this constructor, an initial osculating orbit is considered.</p>
150 * @param initialOrbit initial orbit
151 * @param attitude attitude provider
152 * @param mass spacecraft mass
153 * @param provider for un-normalized zonal coefficients
154 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
155 * @param m2Value value of empirical drag coefficient in rad/s².
156 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
157 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement,
158 * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType, double)
159 */
160 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
161 final AttitudeProvider attitude,
162 final T mass,
163 final UnnormalizedSphericalHarmonicsProvider provider,
164 final UnnormalizedSphericalHarmonics harmonics,
165 final double m2Value) {
166 this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
167 harmonics.getUnnormalizedCnm(2, 0),
168 harmonics.getUnnormalizedCnm(3, 0),
169 harmonics.getUnnormalizedCnm(4, 0),
170 harmonics.getUnnormalizedCnm(5, 0),
171 m2Value);
172 }
173
174 /** Build a propagator from orbit and potential.
175 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
176 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
177 * are related to both the normalized coefficients
178 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
179 * and the J<sub>n</sub> one as follows:</p>
180 *
181 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
182 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
183 *
184 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
185 *
186 * <p>Using this constructor, an initial osculating orbit is considered.</p>
187 *
188 * @param initialOrbit initial orbit
189 * @param referenceRadius reference radius of the Earth for the potential model (m)
190 * @param mu central attraction coefficient (m³/s²)
191 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
192 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
193 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
194 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
195 * @param m2Value value of empirical drag coefficient in rad/s².
196 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
197 * @see org.orekit.utils.Constants
198 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, double, CalculusFieldElement, double, double, double, double, double)
199 */
200 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
201 final double referenceRadius,
202 final T mu,
203 final double c20,
204 final double c30,
205 final double c40,
206 final double c50,
207 final double m2Value) {
208 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
209 initialOrbit.getMu().newInstance(DEFAULT_MASS),
210 referenceRadius, mu, c20, c30, c40, c50, m2Value);
211 }
212
213 /** Build a propagator from orbit, mass and potential provider.
214 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
215 *
216 * <p>Using this constructor, an initial osculating orbit is considered.</p>
217 *
218 * @param initialOrbit initial orbit
219 * @param mass spacecraft mass
220 * @param provider for un-normalized zonal coefficients
221 * @param m2Value value of empirical drag coefficient in rad/s².
222 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
223 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, double)
224 */
225 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
226 final T mass,
227 final UnnormalizedSphericalHarmonicsProvider provider,
228 final double m2Value) {
229 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
230 mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), m2Value);
231 }
232
233 /** Build a propagator from orbit, mass and potential.
234 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
235 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
236 * are related to both the normalized coefficients
237 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
238 * and the J<sub>n</sub> one as follows:</p>
239 *
240 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
241 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
242 *
243 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
244 *
245 * <p>Using this constructor, an initial osculating orbit is considered.</p>
246 *
247 * @param initialOrbit initial orbit
248 * @param mass spacecraft mass
249 * @param referenceRadius reference radius of the Earth for the potential model (m)
250 * @param mu central attraction coefficient (m³/s²)
251 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
252 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
253 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
254 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
255 * @param m2Value value of empirical drag coefficient in rad/s².
256 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
257 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, double)
258 */
259 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit, final T mass,
260 final double referenceRadius, final T mu,
261 final double c20, final double c30, final double c40,
262 final double c50, final double m2Value) {
263 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
264 mass, referenceRadius, mu, c20, c30, c40, c50, m2Value);
265 }
266
267 /** Build a propagator from orbit, attitude provider and potential provider.
268 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
269 * <p>Using this constructor, an initial osculating orbit is considered.</p>
270 * @param initialOrbit initial orbit
271 * @param attitudeProv attitude provider
272 * @param provider for un-normalized zonal coefficients
273 * @param m2Value value of empirical drag coefficient in rad/s².
274 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
275 */
276 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
277 final AttitudeProvider attitudeProv,
278 final UnnormalizedSphericalHarmonicsProvider provider,
279 final double m2Value) {
280 this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
281 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), m2Value);
282 }
283
284 /** Build a propagator from orbit, attitude provider and potential.
285 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
286 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
287 * are related to both the normalized coefficients
288 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
289 * and the J<sub>n</sub> one as follows:</p>
290 *
291 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
292 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
293 *
294 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
295 *
296 * <p>Using this constructor, an initial osculating orbit is considered.</p>
297 *
298 * @param initialOrbit initial orbit
299 * @param attitudeProv attitude provider
300 * @param referenceRadius reference radius of the Earth for the potential model (m)
301 * @param mu central attraction coefficient (m³/s²)
302 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
303 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
304 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
305 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
306 * @param m2Value value of empirical drag coefficient in rad/s².
307 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
308 */
309 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
310 final AttitudeProvider attitudeProv,
311 final double referenceRadius, final T mu,
312 final double c20, final double c30, final double c40,
313 final double c50, final double m2Value) {
314 this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
315 referenceRadius, mu, c20, c30, c40, c50, m2Value);
316 }
317
318 /** Build a propagator from orbit, attitude provider, mass and potential provider.
319 * <p>Using this constructor, an initial osculating orbit is considered.</p>
320 * @param initialOrbit initial orbit
321 * @param attitudeProv attitude provider
322 * @param mass spacecraft mass
323 * @param provider for un-normalized zonal coefficients
324 * @param m2Value value of empirical drag coefficient in rad/s².
325 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
326 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
327 */
328 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
329 final AttitudeProvider attitudeProv,
330 final T mass,
331 final UnnormalizedSphericalHarmonicsProvider provider,
332 final double m2Value) {
333 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), m2Value);
334 }
335
336 /** Build a propagator from orbit, attitude provider, mass and potential.
337 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
338 * are related to both the normalized coefficients
339 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
340 * and the J<sub>n</sub> one as follows:</p>
341 *
342 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
343 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
344 *
345 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
346 *
347 * <p>Using this constructor, an initial osculating orbit is considered.</p>
348 *
349 * @param initialOrbit initial orbit
350 * @param attitudeProv attitude provider
351 * @param mass spacecraft mass
352 * @param referenceRadius reference radius of the Earth for the potential model (m)
353 * @param mu central attraction coefficient (m³/s²)
354 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
355 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
356 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
357 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
358 * @param m2Value value of empirical drag coefficient in rad/s².
359 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
360 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, PropagationType, double)
361 */
362 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
363 final AttitudeProvider attitudeProv,
364 final T mass,
365 final double referenceRadius, final T mu,
366 final double c20, final double c30, final double c40,
367 final double c50, final double m2Value) {
368 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, PropagationType.OSCULATING, m2Value);
369 }
370
371
372 /** Build a propagator from orbit and potential provider.
373 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
374 *
375 * <p>Using this constructor, it is possible to define the initial orbit as
376 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
377 *
378 * @param initialOrbit initial orbit
379 * @param provider for un-normalized zonal coefficients
380 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
381 * @param m2Value value of empirical drag coefficient in rad/s².
382 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
383 */
384 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
385 final UnnormalizedSphericalHarmonicsProvider provider,
386 final PropagationType initialType,
387 final double m2Value) {
388 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
389 initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
390 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, m2Value);
391 }
392
393 /** Build a propagator from orbit, attitude provider, mass and potential provider.
394 * <p>Using this constructor, it is possible to define the initial orbit as
395 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
396 * @param initialOrbit initial orbit
397 * @param attitudeProv attitude provider
398 * @param mass spacecraft mass
399 * @param provider for un-normalized zonal coefficients
400 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
401 * @param m2Value value of empirical drag coefficient in rad/s².
402 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
403 */
404 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
405 final AttitudeProvider attitudeProv,
406 final T mass,
407 final UnnormalizedSphericalHarmonicsProvider provider,
408 final PropagationType initialType,
409 final double m2Value) {
410 this(initialOrbit, attitudeProv, mass, provider,
411 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, m2Value);
412 }
413
414 /**
415 * Private helper constructor.
416 * <p>Using this constructor, it is possible to define the initial orbit as
417 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
418 * @param initialOrbit initial orbit
419 * @param attitude attitude provider
420 * @param mass spacecraft mass
421 * @param provider for un-normalized zonal coefficients
422 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
423 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
424 * @param m2Value value of empirical drag coefficient in rad/s².
425 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
426 */
427 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
428 final AttitudeProvider attitude,
429 final T mass,
430 final UnnormalizedSphericalHarmonicsProvider provider,
431 final UnnormalizedSphericalHarmonics harmonics,
432 final PropagationType initialType,
433 final double m2Value) {
434 this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
435 harmonics.getUnnormalizedCnm(2, 0),
436 harmonics.getUnnormalizedCnm(3, 0),
437 harmonics.getUnnormalizedCnm(4, 0),
438 harmonics.getUnnormalizedCnm(5, 0),
439 initialType, m2Value);
440 }
441
442 /** Build a propagator from orbit, attitude provider, mass and potential.
443 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
444 * are related to both the normalized coefficients
445 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
446 * and the J<sub>n</sub> one as follows:</p>
447 *
448 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
449 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
450 *
451 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
452 *
453 * <p>Using this constructor, it is possible to define the initial orbit as
454 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
455 *
456 * @param initialOrbit initial orbit
457 * @param attitudeProv attitude provider
458 * @param mass spacecraft mass
459 * @param referenceRadius reference radius of the Earth for the potential model (m)
460 * @param mu central attraction coefficient (m³/s²)
461 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
462 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
463 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
464 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
465 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
466 * @param m2Value value of empirical drag coefficient in rad/s².
467 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
468 */
469 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
470 final AttitudeProvider attitudeProv,
471 final T mass,
472 final double referenceRadius, final T mu,
473 final double c20, final double c30, final double c40,
474 final double c50,
475 final PropagationType initialType,
476 final double m2Value) {
477 this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
478 c20, c30, c40, c50, initialType, m2Value,
479 new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
480 BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
481 FixedPointConverter.DEFAULT_DAMPING));
482 }
483
484 /** Build a propagator from orbit, attitude provider, mass and potential.
485 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
486 * are related to both the normalized coefficients
487 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
488 * and the J<sub>n</sub> one as follows:</p>
489 *
490 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
491 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
492 *
493 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
494 *
495 * <p>Using this constructor, it is possible to define the initial orbit as
496 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
497 *
498 * @param initialOrbit initial orbit
499 * @param attitudeProv attitude provider
500 * @param mass spacecraft mass
501 * @param referenceRadius reference radius of the Earth for the potential model (m)
502 * @param mu central attraction coefficient (m³/s²)
503 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
504 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
505 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
506 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
507 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
508 * @param m2Value value of empirical drag coefficient in rad/s².
509 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
510 * @param epsilon convergence threshold for mean parameters conversion
511 * @param maxIterations maximum iterations for mean parameters conversion
512 * @since 11.2
513 */
514 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
515 final AttitudeProvider attitudeProv,
516 final T mass,
517 final double referenceRadius,
518 final T mu,
519 final double c20,
520 final double c30,
521 final double c40,
522 final double c50,
523 final PropagationType initialType,
524 final double m2Value,
525 final double epsilon,
526 final int maxIterations) {
527 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50,
528 initialType, m2Value, new FixedPointConverter(epsilon, maxIterations,
529 FixedPointConverter.DEFAULT_DAMPING));
530 }
531
532 /** Build a propagator from orbit, attitude provider, mass and potential.
533 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
534 * are related to both the normalized coefficients
535 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
536 * and the J<sub>n</sub> one as follows:</p>
537 *
538 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
539 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
540 *
541 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
542 *
543 * <p>Using this constructor, it is possible to define the initial orbit as
544 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
545 *
546 * @param initialOrbit initial orbit
547 * @param attitudeProv attitude provider
548 * @param mass spacecraft mass
549 * @param referenceRadius reference radius of the Earth for the potential model (m)
550 * @param mu central attraction coefficient (m³/s²)
551 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
552 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
553 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
554 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
555 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
556 * @param m2Value value of empirical drag coefficient in rad/s².
557 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
558 * @param converter osculating to mean orbit converter
559 * @since 13.0
560 */
561 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
562 final AttitudeProvider attitudeProv,
563 final T mass,
564 final double referenceRadius,
565 final T mu,
566 final double c20,
567 final double c30,
568 final double c40,
569 final double c50,
570 final PropagationType initialType,
571 final double m2Value,
572 final OsculatingToMeanConverter converter) {
573
574 super(mass.getField(), attitudeProv);
575
576 // store model coefficients
577 this.referenceRadius = referenceRadius;
578 this.mu = mu;
579 this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};
580
581 // initialize M2 driver
582 this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, m2Value, SCALE,
583 Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY,
584 TimeInterval.UNLIMITED);
585
586 // compute mean parameters if needed
587 resetInitialState(new FieldSpacecraftState<>(initialOrbit,
588 attitudeProv.getAttitude(initialOrbit,
589 initialOrbit.getDate(),
590 initialOrbit.getFrame())).withMass(mass),
591 initialType, converter);
592
593 }
594
595 /** Conversion from osculating to mean orbit.
596 * <p>
597 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
598 * osculating SpacecraftState in input.
599 * </p>
600 * <p>
601 * Since the osculating orbit is obtained with the computation of
602 * short-periodic variation, the resulting output will depend on
603 * both the gravity field parameterized in input and the
604 * atmospheric drag represented by the {@code m2} parameter.
605 * </p>
606 * <p>
607 * The computation is done through a fixed-point iteration process.
608 * </p>
609 * @param <T> type of the filed elements
610 * @param osculating osculating orbit to convert
611 * @param provider for un-normalized zonal coefficients
612 * @param harmonics {@code provider.onDate(osculating.getDate())}
613 * @param m2Value value of empirical drag coefficient in rad/s².
614 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
615 * @return mean orbit in a Brouwer-Lyddane sense
616 * @since 11.2
617 */
618 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
619 final UnnormalizedSphericalHarmonicsProvider provider,
620 final UnnormalizedSphericalHarmonics harmonics,
621 final double m2Value) {
622 return computeMeanOrbit(osculating, provider, harmonics, m2Value,
623 BrouwerLyddanePropagator.EPSILON_DEFAULT,
624 BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT);
625 }
626
627 /** Conversion from osculating to mean orbit.
628 * <p>
629 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
630 * osculating SpacecraftState in input.
631 * </p>
632 * <p>
633 * Since the osculating orbit is obtained with the computation of
634 * short-periodic variation, the resulting output will depend on
635 * both the gravity field parameterized in input and the
636 * atmospheric drag represented by the {@code m2} parameter.
637 * </p>
638 * <p>
639 * The computation is done through a fixed-point iteration process.
640 * </p>
641 * @param <T> type of the filed elements
642 * @param osculating osculating orbit to convert
643 * @param provider for un-normalized zonal coefficients
644 * @param harmonics {@code provider.onDate(osculating.getDate())}
645 * @param m2Value value of empirical drag coefficient in rad/s².
646 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
647 * @param epsilon convergence threshold for mean parameters conversion
648 * @param maxIterations maximum iterations for mean parameters conversion
649 * @return mean orbit in a Brouwer-Lyddane sense
650 * @since 11.2
651 */
652 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
653 final UnnormalizedSphericalHarmonicsProvider provider,
654 final UnnormalizedSphericalHarmonics harmonics,
655 final double m2Value,
656 final double epsilon,
657 final int maxIterations) {
658 return computeMeanOrbit(osculating,
659 provider.getAe(), provider.getMu(),
660 harmonics.getUnnormalizedCnm(2, 0),
661 harmonics.getUnnormalizedCnm(3, 0),
662 harmonics.getUnnormalizedCnm(4, 0),
663 harmonics.getUnnormalizedCnm(5, 0),
664 m2Value, epsilon, maxIterations);
665 }
666
667 /** Conversion from osculating to mean orbit.
668 * <p>
669 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
670 * osculating SpacecraftState in input.
671 * </p>
672 * <p>
673 * Since the osculating orbit is obtained with the computation of
674 * short-periodic variation, the resulting output will depend on
675 * both the gravity field parameterized in input and the
676 * atmospheric drag represented by the {@code m2} parameter.
677 * </p>
678 * <p>
679 * The computation is done through a fixed-point iteration process.
680 * </p>
681 * @param <T> type of the filed elements
682 * @param osculating osculating orbit to convert
683 * @param referenceRadius reference radius of the Earth for the potential model (m)
684 * @param mu central attraction coefficient (m³/s²)
685 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
686 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
687 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
688 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
689 * @param m2Value value of empirical drag coefficient in rad/s².
690 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
691 * @param epsilon convergence threshold for mean parameters conversion
692 * @param maxIterations maximum iterations for mean parameters conversion
693 * @return mean orbit in a Brouwer-Lyddane sense
694 * @since 11.2
695 */
696 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
697 final double referenceRadius,
698 final double mu,
699 final double c20,
700 final double c30,
701 final double c40,
702 final double c50,
703 final double m2Value,
704 final double epsilon,
705 final int maxIterations) {
706 // Build a fixed-point converter
707 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
708 FixedPointConverter.DEFAULT_DAMPING);
709 return computeMeanOrbit(osculating, referenceRadius, mu, c20, c30, c40, c50, m2Value, converter);
710 }
711
712 /** Conversion from osculating to mean orbit.
713 * <p>
714 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
715 * osculating SpacecraftState in input.
716 * </p>
717 * <p>
718 * Since the osculating orbit is obtained with the computation of
719 * short-periodic variation, the resulting output will depend on
720 * both the gravity field parameterized in input and the
721 * atmospheric drag represented by the {@code m2} parameter.
722 * </p>
723 * <p>
724 * The computation is done through the given osculating to mean orbit converter.
725 * </p>
726 * @param <T> type of the filed elements
727 * @param osculating osculating orbit to convert
728 * @param referenceRadius reference radius of the Earth for the potential model (m)
729 * @param mu central attraction coefficient (m³/s²)
730 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
731 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
732 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
733 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
734 * @param m2Value value of empirical drag coefficient in rad/s².
735 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
736 * @param converter osculating to mean orbit converter
737 * @return mean orbit in a Brouwer-Lyddane sense
738 * @since 13.0
739 */
740 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
741 final double referenceRadius,
742 final double mu,
743 final double c20,
744 final double c30,
745 final double c40,
746 final double c50,
747 final double m2Value,
748 final OsculatingToMeanConverter converter) {
749 // Set BL as the mean theory for converting
750 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu, c20, c30, c40, c50, m2Value);
751 converter.setMeanTheory(theory);
752 return (FieldKeplerianOrbit<T>) OrbitParamsType.KEPLERIAN.convertType(converter.convertToMean(osculating));
753 }
754
755 /** {@inheritDoc}
756 * <p>The new initial state to consider
757 * must be defined with an osculating orbit.</p>
758 * @see #resetInitialState(FieldSpacecraftState, PropagationType)
759 */
760 @Override
761 public void resetInitialState(final FieldSpacecraftState<T> state) {
762 resetInitialState(state, PropagationType.OSCULATING);
763 }
764
765 /** Reset the propagator initial state.
766 * @param state new initial state to consider
767 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
768 */
769 public void resetInitialState(final FieldSpacecraftState<T> state,
770 final PropagationType stateType) {
771 final OsculatingToMeanConverter converter = new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
772 BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
773 FixedPointConverter.DEFAULT_DAMPING);
774 resetInitialState(state, stateType, converter);
775 }
776
777 /** Reset the propagator initial state.
778 * @param state new initial state to consider
779 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
780 * @param epsilon convergence threshold for mean parameters conversion
781 * @param maxIterations maximum iterations for mean parameters conversion
782 * @since 11.2
783 */
784 public void resetInitialState(final FieldSpacecraftState<T> state,
785 final PropagationType stateType,
786 final double epsilon,
787 final int maxIterations) {
788 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
789 FixedPointConverter.DEFAULT_DAMPING);
790 resetInitialState(state, stateType, converter);
791 }
792
793 /** Reset the propagator initial state.
794 * @param state new initial state to consider
795 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
796 * @param converter osculating to mean orbit converter
797 * @since 13.0
798 */
799 public void resetInitialState(final FieldSpacecraftState<T> state,
800 final PropagationType stateType,
801 final OsculatingToMeanConverter converter) {
802 super.resetInitialState(state);
803 FieldKeplerianOrbit<T> keplerian = (FieldKeplerianOrbit<T>) OrbitParamsType.KEPLERIAN.convertType(state.getOrbit());
804 if (stateType == PropagationType.OSCULATING) {
805 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu.getReal(),
806 ck0[2], ck0[3], ck0[4], ck0[5],
807 getM2());
808 converter.setMeanTheory(theory);
809 keplerian = (FieldKeplerianOrbit<T>) OrbitParamsType.KEPLERIAN.convertType(converter.convertToMean(keplerian));
810 }
811 this.initialModel = new FieldBLModel<>(keplerian, state.getMass(), referenceRadius, mu, ck0);
812 this.models = new FieldTimeSpanMap<>(initialModel, state.getMass().getField());
813 }
814
815 /** {@inheritDoc} */
816 @Override
817 protected void resetIntermediateState(final FieldSpacecraftState<T> state,
818 final boolean forward) {
819 final OsculatingToMeanConverter converter = new FixedPointConverter(BrouwerLyddanePropagator.EPSILON_DEFAULT,
820 BrouwerLyddanePropagator.MAX_ITERATIONS_DEFAULT,
821 FixedPointConverter.DEFAULT_DAMPING);
822 resetIntermediateState(state, forward, converter);
823 }
824
825 /** Reset an intermediate state.
826 * @param state new intermediate state to consider
827 * @param forward if true, the intermediate state is valid for
828 * propagations after itself
829 * @param epsilon convergence threshold for mean parameters conversion
830 * @param maxIterations maximum iterations for mean parameters conversion
831 * @since 11.2
832 */
833 protected void resetIntermediateState(final FieldSpacecraftState<T> state,
834 final boolean forward,
835 final double epsilon,
836 final int maxIterations) {
837 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
838 FixedPointConverter.DEFAULT_DAMPING);
839 resetIntermediateState(state, forward, converter);
840 }
841
842 /** Reset an intermediate state.
843 * @param state new intermediate state to consider
844 * @param forward if true, the intermediate state is valid for
845 * propagations after itself
846 * @param converter osculating to mean orbit converter
847 * @since 13.0
848 */
849 protected void resetIntermediateState(final FieldSpacecraftState<T> state,
850 final boolean forward,
851 final OsculatingToMeanConverter converter) {
852 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu.getReal(),
853 ck0[2], ck0[3], ck0[4], ck0[5],
854 getM2());
855 converter.setMeanTheory(theory);
856 final FieldKeplerianOrbit<T> mean = (FieldKeplerianOrbit<T>) OrbitParamsType.KEPLERIAN.convertType(converter.convertToMean(state.getOrbit()));
857 final FieldBLModel<T> newModel = new FieldBLModel<>(mean, state.getMass(), referenceRadius, mu, ck0);
858 if (forward) {
859 models.addValidAfter(newModel, state.getDate(), false);
860 } else {
861 models.addValidBefore(newModel, state.getDate(), false);
862 }
863 stateChanged(state);
864 }
865
866 /** {@inheritDoc} */
867 public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
868 // compute Cartesian parameters, taking derivatives into account
869 final FieldBLModel<T> current = models.get(date);
870 return current.propagateParameters(date, parameters);
871 }
872
873 /**
874 * Get the value of the M2 drag parameter.
875 * @return the value of the M2 drag parameter
876 */
877 public double getM2() {
878 return M2Driver.getValue();
879 }
880
881 /** Local class for Brouwer-Lyddane model. */
882 private static class FieldBLModel<T extends CalculusFieldElement<T>> {
883
884 /** Constant mass. */
885 private final T mass;
886
887 /** Central attraction coefficient. */
888 private final T mu;
889
890 /** Brouwer-Lyddane mean orbit. */
891 private final FieldKeplerianOrbit<T> mean;
892
893 // Preprocessed values
894
895 /** Mean mean motion: n0 = √(μ/a")/a". */
896 private final T n0;
897
898 /** η = √(1 - e"²). */
899 private final T n;
900 /** η². */
901 private final T n2;
902 /** η³. */
903 private final T n3;
904 /** η + 1 / (1 + η). */
905 private final T t8;
906
907 /** Secular correction for mean anomaly l: δ<sub>s</sub>l. */
908 private final T dsl;
909 /** Secular correction for periapsis argument g: δ<sub>s</sub>g. */
910 private final T dsg;
911 /** Secular correction for raan h: δ<sub>s</sub>h. */
912 private final T dsh;
913
914 /** Secular rate of change of semi-major axis due to drag. */
915 private final T aRate;
916 /** Secular rate of change of eccentricity due to drag. */
917 private final T eRate;
918
919 // CHECKSTYLE: stop JavadocVariable check
920
921 // Storage for speed-up
922 private final T yp2;
923 private final T ci;
924 private final T si;
925 private final T oneMci2;
926 private final T ci2X3M1;
927
928 // Long periodic corrections factors
929 private final T vle1;
930 private final T vle2;
931 private final T vle3;
932 private final T vli1;
933 private final T vli2;
934 private final T vli3;
935 private final T vll2;
936 private final T vlh1I;
937 private final T vlh2I;
938 private final T vlh3I;
939 private final T vls1;
940 private final T vls2;
941 private final T vls3;
942
943 // CHECKSTYLE: resume JavadocVariable check
944
945 /** Create a model for specified mean orbit.
946 * @param mean mean Fieldorbit
947 * @param mass constant mass
948 * @param referenceRadius reference radius of the central body attraction model (m)
949 * @param mu central attraction coefficient (m³/s²)
950 * @param ck0 un-normalized zonal coefficients
951 */
952 FieldBLModel(final FieldKeplerianOrbit<T> mean, final T mass,
953 final double referenceRadius, final T mu, final double[] ck0) {
954
955 this.mass = mass;
956 this.mu = mu;
957
958 // mean orbit
959 this.mean = mean;
960
961 final T one = mass.getField().getOne();
962
963 // mean eccentricity e"
964 final T epp = mean.getE();
965 if (epp.getReal() >= 1) {
966 // Only for elliptical (e < 1) orbits
967 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
968 epp.getReal());
969 }
970 final T epp2 = epp.square();
971
972 // η
973 n2 = one.subtract(epp2);
974 n = n2.sqrt();
975 n3 = n2.multiply(n);
976 t8 = n.add(one.add(n).reciprocal());
977
978 // mean semi-major axis a"
979 final T app = mean.getA();
980
981 // mean mean motion
982 n0 = mu.divide(app).sqrt().divide(app);
983
984 // ae/a"
985 final T q = app.divide(referenceRadius).reciprocal();
986
987 // γ2'
988 T ql = q.square();
989 T nl = n2.square();
990 yp2 = ql.multiply(-0.5 * ck0[2]).divide(nl);
991 final T yp22 = yp2.square();
992
993 // γ3'
994 ql = ql.multiply(q);
995 nl = nl.multiply(n2);
996 final T yp3 = ql.multiply(ck0[3]).divide(nl);
997
998 // γ4'
999 ql = ql.multiply(q);
1000 nl = nl.multiply(n2);
1001 final T yp4 = ql.multiply(0.375 * ck0[4]).divide(nl);
1002
1003 // γ5'
1004 ql = ql.multiply(q);
1005 nl = nl.multiply(n2);
1006 final T yp5 = ql.multiply(ck0[5]).divide(nl);
1007
1008 // mean inclination I" sin & cos
1009 final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
1010 si = sc.sin();
1011 ci = sc.cos();
1012 final T ci2 = ci.square();
1013 oneMci2 = one.subtract(ci2);
1014 ci2X3M1 = ci2.multiply(3.).subtract(one);
1015 final T ci2X5M1 = ci2.multiply(5.).subtract(one);
1016
1017 // secular corrections
1018 // true anomaly
1019 final T dsl1 = yp2.multiply(n).multiply(1.5);
1020 final T dsl2a = n.multiply(n.multiply(25.).add(16.)).subtract(15.);
1021 final T dsl2b = n.multiply(n.multiply(90.).add(96.)).negate().add(30.);
1022 final T dsl2c = n.multiply(n.multiply(25.).add(144.)).add(105.);
1023 final T dsl21 = dsl2a.add(ci2.multiply(dsl2b.add(ci2.multiply(dsl2c))));
1024 final T dsl2 = ci2X3M1.add(yp2.multiply(0.0625).multiply(dsl21));
1025 final T dsl3 = yp4.multiply(n).multiply(epp2).multiply(0.9375).
1026 multiply(ci2.multiply(35.0).subtract(30.0).multiply(ci2).add(3.));
1027 dsl = dsl1.multiply(dsl2).add(dsl3);
1028
1029 // periapsis argument
1030 final T dsg1 = yp2.multiply(1.5).multiply(ci2X5M1);
1031 final T dsg2a = n.multiply(25.).add(24.).multiply(n).add(-35.);
1032 final T dsg2b = n.multiply(126.).add(192.).multiply(n).negate().add(90.);
1033 final T dsg2c = n.multiply(45.).add(360.).multiply(n).add(385.);
1034 final T dsg21 = dsg2a.add(ci2.multiply(dsg2b.add(ci2.multiply(dsg2c))));
1035 final T dsg2 = yp22.multiply(0.09375).multiply(dsg21);
1036 final T dsg3a = n2.multiply(-9.).add(21.);
1037 final T dsg3b = n2.multiply(126.).add(-270.);
1038 final T dsg3c = n2.multiply(-189.).add(385.);
1039 final T dsg31 = dsg3a.add(ci2.multiply(dsg3b.add(ci2.multiply(dsg3c))));
1040 final T dsg3 = yp4.multiply(0.3125).multiply(dsg31);
1041 dsg = dsg1.add(dsg2).add(dsg3);
1042
1043 // right ascension of ascending node
1044 final T dsh1 = yp2.multiply(-3.);
1045 final T dsh2a = n.multiply(9.).add(12.).multiply(n).add(-5.);
1046 final T dsh2b = n.multiply(5.).add(36.).multiply(n).add(35.);
1047 final T dsh21 = dsh2a.subtract(ci2.multiply(dsh2b));
1048 final T dsh2 = yp22.multiply(0.375).multiply(dsh21);
1049 final T dsh31 = n2.multiply(3.).subtract(5.);
1050 final T dsh32 = ci2.multiply(7.).subtract(3.);
1051 final T dsh3 = yp4.multiply(1.25).multiply(dsh31).multiply(dsh32);
1052 dsh = ci.multiply(dsh1.add(dsh2).add(dsh3));
1053
1054 // secular rates of change due to drag
1055 // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
1056 final T coef = n0.multiply(one.add(dsl)).multiply(3.).reciprocal().multiply(-4);
1057 aRate = coef.multiply(app);
1058 eRate = coef.multiply(epp).multiply(n2);
1059
1060 // singular term 1/(1 - 5 * cos²(I")) replaced by T2 function
1061 final T t2 = T2(ci);
1062
1063 // factors for long periodic corrections
1064 final T fs12 = yp3.divide(yp2);
1065 final T fs13 = yp4.multiply(10).divide(yp2.multiply(3));
1066 final T fs14 = yp5.divide(yp2);
1067
1068 final T ci2Xt2 = ci2.multiply(t2);
1069 final T cA = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(11.)));
1070 final T cB = one.subtract(ci2.multiply(ci2Xt2.multiply(8.) .add(3.)));
1071 final T cC = one.subtract(ci2.multiply(ci2Xt2.multiply(24.) .add(9.)));
1072 final T cD = one.subtract(ci2.multiply(ci2Xt2.multiply(16.) .add(5.)));
1073 final T cE = one.subtract(ci2.multiply(ci2Xt2.multiply(200.).add(33.)));
1074 final T cF = one.subtract(ci2.multiply(ci2Xt2.multiply(40.) .add(9.)));
1075
1076 final T p5p = one.add(ci2Xt2.multiply(ci2Xt2.multiply(20.).add(8.)));
1077 final T p5p2 = one.add(p5p.multiply(2.));
1078 final T p5p4 = one.add(p5p.multiply(4.));
1079 final T p5p10 = one.add(p5p.multiply(10.));
1080
1081 final T e2X3P4 = epp2.multiply(3.).add(4.);
1082 final T ciO1Pci = ci.divide(one.add(ci));
1083 final T oneMci = one.subtract(ci);
1084
1085 final T q1 = (yp2.multiply(cA).subtract(fs13.multiply(cB))).
1086 multiply(0.125);
1087 final T q2 = (yp2.multiply(p5p10).subtract(fs13.multiply(p5p2))).
1088 multiply(epp2).multiply(ci).multiply(0.125);
1089 final T q5 = (fs12.add(e2X3P4.multiply(fs14).multiply(cC).multiply(0.3125))).
1090 multiply(0.25);
1091 final T p2 = p5p2.multiply(epp).multiply(ci).multiply(si).multiply(e2X3P4).multiply(fs14).
1092 multiply(0.46875);
1093 final T p3 = epp.multiply(si).multiply(fs14).multiply(cC).
1094 multiply(0.15625);
1095 final double kf = 35. / 1152.;
1096 final T p4 = epp.multiply(fs14).multiply(cD).
1097 multiply(kf);
1098 final T p5 = epp.multiply(epp2).multiply(ci).multiply(si).multiply(fs14).multiply(p5p4).
1099 multiply(2. * kf);
1100
1101 vle1 = epp.multiply(n2).multiply(q1);
1102 vle2 = n2.multiply(si).multiply(q5);
1103 vle3 = epp.multiply(n2).multiply(si).multiply(p4).multiply(-3.0);
1104
1105 vli1 = epp.multiply(q1).divide(si).negate();
1106 vli2 = epp.multiply(ci).multiply(q5).negate();
1107 vli3 = epp2.multiply(ci).multiply(p4).multiply(-3.0);
1108
1109 vll2 = vle2.add(epp.multiply(n2).multiply(p3).multiply(3.0));
1110
1111 vlh1I = si.multiply(q2).negate();
1112 vlh2I = epp.multiply(ci).multiply(q5).add(si.multiply(p2));
1113 vlh3I = (epp2.multiply(ci).multiply(p4).add(si.multiply(p5))).negate();
1114
1115 vls1 = q1.multiply(n3.subtract(one)).
1116 subtract(q2).
1117 add(epp2.multiply(ci2).multiply(ci2Xt2).multiply(ci2Xt2).
1118 multiply(yp2.subtract(fs13.multiply(0.2))).multiply(25.0)).
1119 subtract(epp2.multiply(yp2.multiply(cE).subtract(fs13.multiply(cF))).multiply(0.0625));
1120
1121 vls2 = epp.multiply(si).multiply(t8.add(ciO1Pci)).multiply(q5).
1122 add((epp2.subtract(n3).multiply(3.).add(11.)).multiply(p3)).
1123 add(oneMci.multiply(p2));
1124
1125 vls3 = si.multiply(p4).multiply(n3.subtract(one).multiply(3.).
1126 subtract(epp2.multiply(ciO1Pci.add(2.)))).
1127 subtract(oneMci.multiply(p5));
1128 }
1129
1130 /**
1131 * Get true anomaly from mean anomaly.
1132 * @param lM the mean anomaly (rad)
1133 * @param ecc the eccentricity
1134 * @return the true anomaly (rad)
1135 */
1136 private FieldUnivariateDerivative1<T> getTrueAnomaly(final FieldUnivariateDerivative1<T> lM,
1137 final FieldUnivariateDerivative1<T> ecc) {
1138
1139 final T zero = mean.getE().getField().getZero();
1140
1141 // reduce M to [-PI PI] interval
1142 final FieldUnivariateDerivative1<T> reducedM = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(lM.getValue(), zero),
1143 lM.getFirstDerivative());
1144
1145 // compute the true anomaly
1146 FieldUnivariateDerivative1<T> lV = FieldKeplerianAnomalyUtility.ellipticMeanToTrue(ecc, lM);
1147
1148 // expand the result back to original range
1149 lV = lV.add(lM.getValue().subtract(reducedM.getValue()));
1150
1151 // Returns the true anomaly
1152 return lV;
1153 }
1154
1155 /**
1156 * This method is used in Brouwer-Lyddane model to avoid singularity at the
1157 * critical inclination (i = 63.4°).
1158 * <p>
1159 * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
1160 * approximate the factor (1.0 - 5.0 * cos²(i))<sup>-1</sup> (causing the singularity)
1161 * by a function, named T2 in the thesis.
1162 * </p>
1163 * @param cosI cosine of the mean inclination
1164 * @return an approximation of (1.0 - 5.0 * cos²(i))<sup>-1</sup> term
1165 */
1166 private T T2(final T cosI) {
1167
1168 // X = (1.0 - 5.0 * cos²(i))
1169 final T x = cosI.square().multiply(-5.0).add(1.0);
1170 final T x2 = x.square();
1171 final T xb = x2.multiply(BETA);
1172
1173 // Eq. 2.48
1174 T sum = x.getField().getZero();
1175 for (int i = 0; i <= 12; i++) {
1176 final double sign = i % 2 == 0 ? +1.0 : -1.0;
1177 sum = sum.add(FastMath.pow(x2, i).
1178 multiply(FastMath.pow(BETA, i)).
1179 multiply(sign).
1180 divide(CombinatoricsUtils.factorialDouble(i + 1)));
1181 }
1182
1183 // Right term of equation 2.47
1184 final T one = x.getField().getOne();
1185 T product = one;
1186 for (int i = 0; i <= 10; i++) {
1187 product = product.multiply(one.add(FastMath.exp(xb.multiply(FastMath.scalb(-1.0, i)))));
1188 }
1189
1190 // Return (Eq. 2.47)
1191 return x.multiply(BETA).multiply(sum).multiply(product);
1192 }
1193
1194 /** Extrapolate an orbit up to a specific target date.
1195 * @param date target date for the orbit
1196 * @param parameters model parameters
1197 * @return propagated parameters
1198 */
1199 public FieldKeplerianOrbit<T> propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {
1200
1201 // Field
1202 final Field<T> field = date.getField();
1203 final T one = field.getOne();
1204 final T zero = field.getZero();
1205
1206 // Empirical drag coefficient M2
1207 final T m2 = parameters[0];
1208
1209 // Keplerian evolution
1210 final FieldUnivariateDerivative1<T> dt = new FieldUnivariateDerivative1<>(date.durationFrom(mean.getDate()), one);
1211 final FieldUnivariateDerivative1<T> not = dt.multiply(n0);
1212
1213 final FieldUnivariateDerivative1<T> dtM2 = dt.multiply(m2);
1214 final FieldUnivariateDerivative1<T> dt2M2 = dt.multiply(dtM2);
1215
1216 // Secular corrections
1217 // -------------------
1218
1219 // semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis)
1220 final FieldUnivariateDerivative1<T> app = dtM2.multiply(aRate).add(mean.getA());
1221
1222 // eccentricity (with drag Eq. 2.45 of Phipps' 1992 thesis) reduced to [0, 1[
1223 final FieldUnivariateDerivative1<T> tmp = dtM2.multiply(eRate).add(mean.getE());
1224 final FieldUnivariateDerivative1<T> epp = FastMath.max(FastMath.min(tmp, MAX_ECC), 0.);
1225
1226 // mean argument of periapsis
1227 final T gp0 = MathUtils.normalizeAngle(mean.getPeriapsisArgument().add(dsg.multiply(not.getValue())), zero);
1228 final T gp1 = dsg.multiply(n0);
1229 final FieldUnivariateDerivative1<T> gpp = new FieldUnivariateDerivative1<>(gp0, gp1);
1230
1231 // mean longitude of ascending node
1232 final T hp0 = MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(dsh.multiply(not.getValue())), zero);
1233 final T hp1 = dsh.multiply(n0);
1234 final FieldUnivariateDerivative1<T> hpp = new FieldUnivariateDerivative1<>(hp0, hp1);
1235
1236 // mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis)
1237 final T lp0 = MathUtils.normalizeAngle(mean.getMeanAnomaly().add(dsl.add(one).multiply(not.getValue())).add(dt2M2.getValue()), zero);
1238 final T lp1 = dsl.add(one).multiply(n0).add(dtM2.multiply(2.0).getValue());
1239 final FieldUnivariateDerivative1<T> lpp = new FieldUnivariateDerivative1<>(lp0, lp1);
1240
1241 // Long period corrections
1242 //------------------------
1243 final FieldSinCos<FieldUnivariateDerivative1<T>> scgpp = gpp.sinCos();
1244 final FieldUnivariateDerivative1<T> cgpp = scgpp.cos();
1245 final FieldUnivariateDerivative1<T> sgpp = scgpp.sin();
1246 final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gpp = gpp.multiply(2).sinCos();
1247 final FieldUnivariateDerivative1<T> c2gpp = sc2gpp.cos();
1248 final FieldUnivariateDerivative1<T> s2gpp = sc2gpp.sin();
1249 final FieldSinCos<FieldUnivariateDerivative1<T>> sc3gpp = gpp.multiply(3).sinCos();
1250 final FieldUnivariateDerivative1<T> c3gpp = sc3gpp.cos();
1251 final FieldUnivariateDerivative1<T> s3gpp = sc3gpp.sin();
1252
1253 // δ1e
1254 final FieldUnivariateDerivative1<T> d1e = c2gpp.multiply(vle1).
1255 add(sgpp.multiply(vle2)).
1256 add(s3gpp.multiply(vle3));
1257
1258 // δ1I
1259 FieldUnivariateDerivative1<T> d1I = sgpp.multiply(vli2).
1260 add(s3gpp.multiply(vli3));
1261 // Pseudo singular term, not to add if I" is zero
1262 if (Double.isFinite(vli1.getReal())) {
1263 d1I = d1I.add(c2gpp.multiply(vli1));
1264 }
1265
1266 // e"δ1l
1267 final FieldUnivariateDerivative1<T> eppd1l = s2gpp.multiply(vle1).
1268 subtract(cgpp.multiply(vll2)).
1269 subtract(c3gpp.multiply(vle3)).
1270 multiply(n);
1271
1272 // sinI"δ1h
1273 final FieldUnivariateDerivative1<T> sIppd1h = s2gpp.multiply(vlh1I).
1274 add(cgpp.multiply(vlh2I)).
1275 add(c3gpp.multiply(vlh3I));
1276
1277 // δ1z = δ1l + δ1g + δ1h
1278 final FieldUnivariateDerivative1<T> d1z = s2gpp.multiply(vls1).
1279 add(cgpp.multiply(vls2)).
1280 add(c3gpp.multiply(vls3));
1281
1282 // Short period corrections
1283 // ------------------------
1284
1285 // true anomaly
1286 final FieldUnivariateDerivative1<T> fpp = getTrueAnomaly(lpp, epp);
1287 final FieldSinCos<FieldUnivariateDerivative1<T>> scfpp = fpp.sinCos();
1288 final FieldUnivariateDerivative1<T> cfpp = scfpp.cos();
1289 final FieldUnivariateDerivative1<T> sfpp = scfpp.sin();
1290
1291 // e"sin(f')
1292 final FieldUnivariateDerivative1<T> eppsfpp = epp.multiply(sfpp);
1293 // e"cos(f')
1294 final FieldUnivariateDerivative1<T> eppcfpp = epp.multiply(cfpp);
1295 // 1 + e"cos(f')
1296 final FieldUnivariateDerivative1<T> eppcfppP1 = eppcfpp.add(1);
1297 // 2 + e"cos(f')
1298 final FieldUnivariateDerivative1<T> eppcfppP2 = eppcfpp.add(2);
1299 // 3 + e"cos(f')
1300 final FieldUnivariateDerivative1<T> eppcfppP3 = eppcfpp.add(3);
1301 // (1 + e"cos(f'))³
1302 final FieldUnivariateDerivative1<T> eppcfppP1_3 = eppcfppP1.square().multiply(eppcfppP1);
1303
1304 // 2g"
1305 final FieldUnivariateDerivative1<T> g2 = gpp.multiply(2);
1306
1307 // 2g" + f"
1308 final FieldUnivariateDerivative1<T> g2f = g2.add(fpp);
1309 final FieldSinCos<FieldUnivariateDerivative1<T>> sc2gf = g2f.sinCos();
1310 final FieldUnivariateDerivative1<T> c2gf = sc2gf.cos();
1311 final FieldUnivariateDerivative1<T> s2gf = sc2gf.sin();
1312 final FieldUnivariateDerivative1<T> eppc2gf = epp.multiply(c2gf);
1313 final FieldUnivariateDerivative1<T> epps2gf = epp.multiply(s2gf);
1314
1315 // 2g" + 2f"
1316 final FieldUnivariateDerivative1<T> g2f2 = g2.add(fpp.multiply(2));
1317 final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g2f = g2f2.sinCos();
1318 final FieldUnivariateDerivative1<T> c2g2f = sc2g2f.cos();
1319 final FieldUnivariateDerivative1<T> s2g2f = sc2g2f.sin();
1320
1321 // 2g" + 3f"
1322 final FieldUnivariateDerivative1<T> g2f3 = g2.add(fpp.multiply(3));
1323 final FieldSinCos<FieldUnivariateDerivative1<T>> sc2g3f = g2f3.sinCos();
1324 final FieldUnivariateDerivative1<T> c2g3f = sc2g3f.cos();
1325 final FieldUnivariateDerivative1<T> s2g3f = sc2g3f.sin();
1326
1327 // e"cos(2g" + 3f")
1328 final FieldUnivariateDerivative1<T> eppc2g3f = epp.multiply(c2g3f);
1329 // e"sin(2g" + 3f")
1330 final FieldUnivariateDerivative1<T> epps2g3f = epp.multiply(s2g3f);
1331
1332 // f" + e"sin(f") - l"
1333 final FieldUnivariateDerivative1<T> w17 = fpp.add(eppsfpp).subtract(lpp);
1334
1335 // ((e"cos(f") + 3)e"cos(f") + 3)cos(f")
1336 final FieldUnivariateDerivative1<T> w20 = cfpp.multiply(eppcfppP3.multiply(eppcfpp).add(3.));
1337
1338 // 3sin(2g" + 2f") + 3e"sin(2g" + f") + e"sin(2g" + f")
1339 final FieldUnivariateDerivative1<T> w21 = s2g2f.add(epps2gf).multiply(3).add(epps2g3f);
1340
1341 // (1 + e"cos(f"))(2 + e"cos(f"))/η²
1342 final FieldUnivariateDerivative1<T> w22 = eppcfppP1.multiply(eppcfppP2).divide(n2);
1343
1344 // sinCos(I"/2)
1345 final FieldSinCos<T> sci = FastMath.sinCos(mean.getI().divide(2.));
1346 final T siO2 = sci.sin();
1347 final T ciO2 = sci.cos();
1348
1349 // δ2a
1350 final FieldUnivariateDerivative1<T> d2a = app.multiply(yp2).divide(n2).
1351 multiply(eppcfppP1_3.subtract(n3).multiply(ci2X3M1).
1352 add(c2g2f.multiply(eppcfppP1_3).multiply(oneMci2).multiply(3.)));
1353
1354 // δ2e
1355 final FieldUnivariateDerivative1<T> d2e = (w20.add(epp.multiply(t8))).multiply(ci2X3M1).
1356 add((w20.add(epp.multiply(c2g2f))).multiply(oneMci2.multiply(3))).
1357 subtract((eppc2gf.multiply(3).add(eppc2g3f)).multiply(oneMci2.multiply(n2))).
1358 multiply(yp2.multiply(0.5));
1359
1360 // δ2I
1361 final FieldUnivariateDerivative1<T> d2I = ((c2g2f.add(eppc2gf)).multiply(3).add(eppc2g3f)).
1362 multiply(yp2.divide(2.).multiply(ci).multiply(si));
1363
1364 // e"δ2l
1365 final FieldUnivariateDerivative1<T> eppd2l = (w22.add(1).multiply(sfpp).multiply(oneMci2).multiply(2.).
1366 add((w22.subtract(1).negate().multiply(s2gf)).
1367 add(w22.add(1. / 3.).multiply(s2g3f)).
1368 multiply(oneMci2.multiply(3.)))).
1369 multiply(yp2.divide(4.).multiply(n3)).negate();
1370
1371 // sinI"δ2h
1372 final FieldUnivariateDerivative1<T> sIppd2h = (w21.subtract(w17.multiply(6))).
1373 multiply(yp2).multiply(ci).multiply(si).divide(2.);
1374
1375 // δ2z = δ2l + δ2g + δ2h
1376 final T ttt = one.add(ci.multiply(ci.multiply(-5).add(2.)));
1377 final FieldUnivariateDerivative1<T> d2z = (epp.multiply(eppd2l).multiply(t8.subtract(one)).divide(n3).
1378 add(w17.multiply(ttt).multiply(6).subtract(w21.multiply(ttt.add(2.))).
1379 multiply(yp2.divide(4.)))).
1380 negate();
1381
1382 // Assembling elements
1383 // -------------------
1384
1385 // e" + δe
1386 final FieldUnivariateDerivative1<T> de = epp.add(d1e).add(d2e);
1387
1388 // e"δl
1389 final FieldUnivariateDerivative1<T> dl = eppd1l.add(eppd2l);
1390
1391 // sin(I"/2)δh = sin(I")δh / cos(I"/2) (singular for I" = π, very unlikely)
1392 final FieldUnivariateDerivative1<T> dh = sIppd1h.add(sIppd2h).divide(ciO2.multiply(2.));
1393
1394 // δI
1395 final FieldUnivariateDerivative1<T> di = d1I.add(d2I).multiply(ciO2).divide(2.).add(siO2);
1396
1397 // z = l" + g" + h" + δ1z + δ2z
1398 final FieldUnivariateDerivative1<T> z = lpp.add(gpp).add(hpp).add(d1z).add(d2z);
1399
1400 // Osculating elements
1401 // -------------------
1402
1403 // Semi-major axis
1404 final FieldUnivariateDerivative1<T> a = app.add(d2a);
1405
1406 // Eccentricity
1407 final FieldUnivariateDerivative1<T> e = FastMath.sqrt(de.square().add(dl.square()));
1408
1409 // Mean anomaly
1410 final FieldSinCos<FieldUnivariateDerivative1<T>> sclpp = lpp.sinCos();
1411 final FieldUnivariateDerivative1<T> clpp = sclpp.cos();
1412 final FieldUnivariateDerivative1<T> slpp = sclpp.sin();
1413 final FieldUnivariateDerivative1<T> l = FastMath.atan2(de.multiply(slpp).add(dl.multiply(clpp)),
1414 de.multiply(clpp).subtract(dl.multiply(slpp)));
1415
1416 // Inclination
1417 final FieldUnivariateDerivative1<T> i = FastMath.acos(di.square().add(dh.square()).multiply(2).negate().add(1.));
1418
1419 // Longitude of ascending node
1420 final FieldSinCos<FieldUnivariateDerivative1<T>> schpp = hpp.sinCos();
1421 final FieldUnivariateDerivative1<T> chpp = schpp.cos();
1422 final FieldUnivariateDerivative1<T> shpp = schpp.sin();
1423 final FieldUnivariateDerivative1<T> h = FastMath.atan2(di.multiply(shpp).add(dh.multiply(chpp)),
1424 di.multiply(chpp).subtract(dh.multiply(shpp)));
1425
1426 // Argument of periapsis
1427 final FieldUnivariateDerivative1<T> g = z.subtract(l).subtract(h);
1428
1429 // Return a Keplerian orbit
1430 return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
1431 g.getValue(), h.getValue(), l.getValue(),
1432 a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
1433 g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
1434 PositionAngleType.MEAN, mean.getFrame(), date, this.mu);
1435 }
1436 }
1437
1438 /** {@inheritDoc} */
1439 @Override
1440 protected T getMass(final FieldAbsoluteDate<T> date) {
1441 return models.get(date).mass;
1442 }
1443
1444 /** {@inheritDoc} */
1445 @Override
1446 public List<ParameterDriver> getParametersDrivers() {
1447 return Collections.singletonList(M2Driver);
1448 }
1449
1450 }