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