1 /* Copyright 2002-2024 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.analytical;
18
19 import java.util.Collections;
20 import java.util.List;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.analysis.differentiation.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.FieldKeplerianOrbit;
36 import org.orekit.orbits.FieldOrbit;
37 import org.orekit.orbits.FieldKeplerianAnomalyUtility;
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.time.AbsoluteDate;
44 import org.orekit.time.FieldAbsoluteDate;
45 import org.orekit.utils.FieldTimeSpanMap;
46 import org.orekit.utils.ParameterDriver;
47
48 /** This class propagates a {@link org.orekit.propagation.FieldSpacecraftState}
49 * using the analytical Brouwer-Lyddane model (from J2 to J5 zonal harmonics).
50 * <p>
51 * At the opposite of the {@link FieldEcksteinHechlerPropagator}, the Brouwer-Lyddane model is
52 * suited for elliptical orbits, there is no problem having a rather small eccentricity or inclination
53 * (Lyddane helped to solve this issue with the Brouwer model). Singularity for the critical
54 * inclination i = 63.4° is avoided using the method developed in Warren Phipps' 1992 thesis.
55 * <p>
56 * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
57 * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
58 * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
59 * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient
60 * {@link #M2Driver}.
61 *
62 * Usually, M2 is adjusted during an orbit determination process and it represents the
63 * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
64 * The behavior of M2 is close to the {@link FieldTLE#getBStar()} parameter for the TLE.
65 *
66 * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track secular
67 * effects are not considered in the dynamical model. Typical values for M2 are not known.
68 * 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).
69 * The unit of M2 is rad/s².
70 *
71 * The along-track effects, represented by the secular rates of the mean semi-major axis
72 * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis.
73 *
74 * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag.
75 * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959."
76 *
77 * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the
78 * artificial satellite. The Astronomical Journal 68 (1963): 555."
79 *
80 * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center
81 * (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992."
82 *
83 * @author Melina Vanel
84 * @author Bryan Cazabonne
85 * @since 11.1
86 * @param <T> type of the field elements
87 */
88 public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> extends FieldAbstractAnalyticalPropagator<T> {
89
90 /** Default convergence threshold for mean parameters conversion. */
91 private static final double EPSILON_DEFAULT = 1.0e-13;
92
93 /** Default value for maxIterations. */
94 private static final int MAX_ITERATIONS_DEFAULT = 200;
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, -20);
103
104 /** Beta constant used by T2 function. */
105 private static final double BETA = FastMath.scalb(100, -11);
106
107 /** Initial Brouwer-Lyddane model. */
108 private FieldBLModel<T> initialModel;
109
110 /** All models. */
111 private transient FieldTimeSpanMap<FieldBLModel<T>, T> models;
112
113 /** Reference radius of the central body attraction model (m). */
114 private double referenceRadius;
115
116 /** Central attraction coefficient (m³/s²). */
117 private T mu;
118
119 /** Un-normalized zonal coefficients. */
120 private double[] ck0;
121
122 /** Empirical coefficient used in the drag modeling. */
123 private final ParameterDriver M2Driver;
124
125 /** Build a propagator from orbit and potential provider.
126 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
127 *
128 * <p>Using this constructor, an initial osculating orbit is considered.</p>
129 *
130 * @param initialOrbit initial orbit
131 * @param provider for un-normalized zonal coefficients
132 * @param M2 value of empirical drag coefficient in rad/s².
133 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
134 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
135 */
136 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
137 final UnnormalizedSphericalHarmonicsProvider provider,
138 final double M2) {
139 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
140 initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
141 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
142 }
143
144 /**
145 * Private helper constructor.
146 * <p>Using this constructor, an initial osculating orbit is considered.</p>
147 * @param initialOrbit initial orbit
148 * @param attitude attitude provider
149 * @param mass spacecraft mass
150 * @param provider for un-normalized zonal coefficients
151 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
152 * @param M2 value of empirical drag coefficient in rad/s².
153 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
154 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement,
155 * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType, double)
156 */
157 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
158 final AttitudeProvider attitude,
159 final T mass,
160 final UnnormalizedSphericalHarmonicsProvider provider,
161 final UnnormalizedSphericalHarmonics harmonics,
162 final double M2) {
163 this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
164 harmonics.getUnnormalizedCnm(2, 0),
165 harmonics.getUnnormalizedCnm(3, 0),
166 harmonics.getUnnormalizedCnm(4, 0),
167 harmonics.getUnnormalizedCnm(5, 0),
168 M2);
169 }
170
171 /** Build a propagator from orbit and potential.
172 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
173 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
174 * are related to both the normalized coefficients
175 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
176 * and the J<sub>n</sub> one as follows:</p>
177 *
178 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
179 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
180 *
181 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
182 *
183 * <p>Using this constructor, an initial osculating orbit is considered.</p>
184 *
185 * @param initialOrbit initial orbit
186 * @param referenceRadius reference radius of the Earth for the potential model (m)
187 * @param mu central attraction coefficient (m³/s²)
188 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
189 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
190 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
191 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
192 * @param M2 value of empirical drag coefficient in rad/s².
193 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
194 * @see org.orekit.utils.Constants
195 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, double, CalculusFieldElement, double, double, double, double, double)
196 */
197 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
198 final double referenceRadius, final T mu,
199 final double c20, final double c30, final double c40,
200 final double c50, final double M2) {
201 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
202 initialOrbit.getMu().newInstance(DEFAULT_MASS),
203 referenceRadius, mu, c20, c30, c40, c50, M2);
204 }
205
206 /** Build a propagator from orbit, mass and potential provider.
207 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
208 *
209 * <p>Using this constructor, an initial osculating orbit is considered.</p>
210 *
211 * @param initialOrbit initial orbit
212 * @param mass spacecraft mass
213 * @param provider for un-normalized zonal coefficients
214 * @param M2 value of empirical drag coefficient in rad/s².
215 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
216 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, double)
217 */
218 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit, final T mass,
219 final UnnormalizedSphericalHarmonicsProvider provider,
220 final double M2) {
221 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
222 mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
223 }
224
225 /** Build a propagator from orbit, mass and potential.
226 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
227 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
228 * are related to both the normalized coefficients
229 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
230 * and the J<sub>n</sub> one as follows:</p>
231 *
232 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
233 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
234 *
235 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
236 *
237 * <p>Using this constructor, an initial osculating orbit is considered.</p>
238 *
239 * @param initialOrbit initial orbit
240 * @param mass spacecraft mass
241 * @param referenceRadius reference radius of the Earth for the potential model (m)
242 * @param mu central attraction coefficient (m³/s²)
243 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
244 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
245 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
246 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
247 * @param M2 value of empirical drag coefficient in rad/s².
248 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
249 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, double)
250 */
251 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit, final T mass,
252 final double referenceRadius, final T mu,
253 final double c20, final double c30, final double c40,
254 final double c50, final double M2) {
255 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
256 mass, referenceRadius, mu, c20, c30, c40, c50, M2);
257 }
258
259 /** Build a propagator from orbit, attitude provider and potential provider.
260 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
261 * <p>Using this constructor, an initial osculating orbit is considered.</p>
262 * @param initialOrbit initial orbit
263 * @param attitudeProv attitude provider
264 * @param provider for un-normalized zonal coefficients
265 * @param M2 value of empirical drag coefficient in rad/s².
266 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
267 */
268 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
269 final AttitudeProvider attitudeProv,
270 final UnnormalizedSphericalHarmonicsProvider provider,
271 final double M2) {
272 this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
273 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
274 }
275
276 /** Build a propagator from orbit, attitude provider and potential.
277 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
278 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
279 * are related to both the normalized coefficients
280 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
281 * and the J<sub>n</sub> one as follows:</p>
282 *
283 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
284 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
285 *
286 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
287 *
288 * <p>Using this constructor, an initial osculating orbit is considered.</p>
289 *
290 * @param initialOrbit initial orbit
291 * @param attitudeProv attitude provider
292 * @param referenceRadius reference radius of the Earth for the potential model (m)
293 * @param mu central attraction coefficient (m³/s²)
294 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
295 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
296 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
297 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
298 * @param M2 value of empirical drag coefficient in rad/s².
299 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
300 */
301 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
302 final AttitudeProvider attitudeProv,
303 final double referenceRadius, final T mu,
304 final double c20, final double c30, final double c40,
305 final double c50, final double M2) {
306 this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
307 referenceRadius, mu, c20, c30, c40, c50, M2);
308 }
309
310 /** Build a propagator from orbit, attitude provider, mass and potential provider.
311 * <p>Using this constructor, an initial osculating orbit is considered.</p>
312 * @param initialOrbit initial orbit
313 * @param attitudeProv attitude provider
314 * @param mass spacecraft mass
315 * @param provider for un-normalized zonal coefficients
316 * @param M2 value of empirical drag coefficient in rad/s².
317 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
318 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType, double)
319 */
320 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
321 final AttitudeProvider attitudeProv,
322 final T mass,
323 final UnnormalizedSphericalHarmonicsProvider provider,
324 final double M2) {
325 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
326 }
327
328 /** Build a propagator from orbit, attitude provider, mass and potential.
329 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
330 * are related to both the normalized coefficients
331 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
332 * and the J<sub>n</sub> one as follows:</p>
333 *
334 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
335 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
336 *
337 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
338 *
339 * <p>Using this constructor, an initial osculating orbit is considered.</p>
340 *
341 * @param initialOrbit initial orbit
342 * @param attitudeProv attitude provider
343 * @param mass spacecraft mass
344 * @param referenceRadius reference radius of the Earth for the potential model (m)
345 * @param mu central attraction coefficient (m³/s²)
346 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
347 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
348 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
349 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
350 * @param M2 value of empirical drag coefficient in rad/s².
351 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
352 * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, PropagationType, double)
353 */
354 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
355 final AttitudeProvider attitudeProv,
356 final T mass,
357 final double referenceRadius, final T mu,
358 final double c20, final double c30, final double c40,
359 final double c50, final double M2) {
360 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, PropagationType.OSCULATING, M2);
361 }
362
363
364 /** Build a propagator from orbit and potential provider.
365 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
366 *
367 * <p>Using this constructor, it is possible to define the initial orbit as
368 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
369 *
370 * @param initialOrbit initial orbit
371 * @param provider for un-normalized zonal coefficients
372 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
373 * @param M2 value of empirical drag coefficient in rad/s².
374 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
375 */
376 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
377 final UnnormalizedSphericalHarmonicsProvider provider,
378 final PropagationType initialType,
379 final double M2) {
380 this(initialOrbit, FrameAlignedProvider.of(initialOrbit.getFrame()),
381 initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
382 provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2);
383 }
384
385 /** Build a propagator from orbit, attitude provider, mass and potential provider.
386 * <p>Using this constructor, it is possible to define the initial orbit as
387 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
388 * @param initialOrbit initial orbit
389 * @param attitudeProv attitude provider
390 * @param mass spacecraft mass
391 * @param provider for un-normalized zonal coefficients
392 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
393 * @param M2 value of empirical drag coefficient in rad/s².
394 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
395 */
396 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
397 final AttitudeProvider attitudeProv,
398 final T mass,
399 final UnnormalizedSphericalHarmonicsProvider provider,
400 final PropagationType initialType,
401 final double M2) {
402 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2);
403 }
404
405 /**
406 * Private helper constructor.
407 * <p>Using this constructor, it is possible to define the initial orbit as
408 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
409 * @param initialOrbit initial orbit
410 * @param attitude attitude provider
411 * @param mass spacecraft mass
412 * @param provider for un-normalized zonal coefficients
413 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
414 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
415 * @param M2 value of empirical drag coefficient in rad/s².
416 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
417 */
418 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
419 final AttitudeProvider attitude,
420 final T mass,
421 final UnnormalizedSphericalHarmonicsProvider provider,
422 final UnnormalizedSphericalHarmonics harmonics,
423 final PropagationType initialType,
424 final double M2) {
425 this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
426 harmonics.getUnnormalizedCnm(2, 0),
427 harmonics.getUnnormalizedCnm(3, 0),
428 harmonics.getUnnormalizedCnm(4, 0),
429 harmonics.getUnnormalizedCnm(5, 0),
430 initialType, M2);
431 }
432
433 /** Build a propagator from orbit, attitude provider, mass and potential.
434 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
435 * are related to both the normalized coefficients
436 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
437 * and the J<sub>n</sub> one as follows:</p>
438 *
439 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
440 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
441 *
442 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
443 *
444 * <p>Using this constructor, it is possible to define the initial orbit as
445 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
446 *
447 * @param initialOrbit initial orbit
448 * @param attitudeProv attitude provider
449 * @param mass spacecraft mass
450 * @param referenceRadius reference radius of the Earth for the potential model (m)
451 * @param mu central attraction coefficient (m³/s²)
452 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
453 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
454 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
455 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
456 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
457 * @param M2 value of empirical drag coefficient in rad/s².
458 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
459 */
460 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
461 final AttitudeProvider attitudeProv,
462 final T mass,
463 final double referenceRadius, final T mu,
464 final double c20, final double c30, final double c40,
465 final double c50,
466 final PropagationType initialType,
467 final double M2) {
468 this(initialOrbit, attitudeProv, mass, referenceRadius, mu,
469 c20, c30, c40, c50, initialType, M2, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
470 }
471
472 /** Build a propagator from orbit, attitude provider, mass and potential.
473 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
474 * are related to both the normalized coefficients
475 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
476 * and the J<sub>n</sub> one as follows:</p>
477 *
478 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
479 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
480 *
481 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
482 *
483 * <p>Using this constructor, it is possible to define the initial orbit as
484 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
485 *
486 * @param initialOrbit initial orbit
487 * @param attitudeProv attitude provider
488 * @param mass spacecraft mass
489 * @param referenceRadius reference radius of the Earth for the potential model (m)
490 * @param mu central attraction coefficient (m³/s²)
491 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
492 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
493 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
494 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
495 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
496 * @param M2 value of empirical drag coefficient in rad/s².
497 * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed
498 * @param epsilon convergence threshold for mean parameters conversion
499 * @param maxIterations maximum iterations for mean parameters conversion
500 * @since 11.2
501 */
502 public FieldBrouwerLyddanePropagator(final FieldOrbit<T> initialOrbit,
503 final AttitudeProvider attitudeProv,
504 final T mass,
505 final double referenceRadius, final T mu,
506 final double c20, final double c30, final double c40,
507 final double c50,
508 final PropagationType initialType,
509 final double M2, final double epsilon, final int maxIterations) {
510
511 super(mass.getField(), attitudeProv);
512
513 // store model coefficients
514 this.referenceRadius = referenceRadius;
515 this.mu = mu;
516 this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50};
517
518 // initialize M2 driver
519 this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE,
520 Double.NEGATIVE_INFINITY,
521 Double.POSITIVE_INFINITY);
522
523 // compute mean parameters if needed
524 resetInitialState(new FieldSpacecraftState<>(initialOrbit,
525 attitudeProv.getAttitude(initialOrbit,
526 initialOrbit.getDate(),
527 initialOrbit.getFrame()),
528 mass),
529 initialType, epsilon, maxIterations);
530
531 }
532
533 /** Conversion from osculating to mean orbit.
534 * <p>
535 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
536 * osculating SpacecraftState in input.
537 * </p>
538 * <p>
539 * Since the osculating orbit is obtained with the computation of
540 * short-periodic variation, the resulting output will depend on
541 * both the gravity field parameterized in input and the
542 * atmospheric drag represented by the {@code m2} parameter.
543 * </p>
544 * <p>
545 * The computation is done through a fixed-point iteration process.
546 * </p>
547 * @param <T> type of the filed elements
548 * @param osculating osculating orbit to convert
549 * @param provider for un-normalized zonal coefficients
550 * @param harmonics {@code provider.onDate(osculating.getDate())}
551 * @param M2Value value of empirical drag coefficient in rad/s².
552 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
553 * @return mean orbit in a Brouwer-Lyddane sense
554 * @since 11.2
555 */
556 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
557 final UnnormalizedSphericalHarmonicsProvider provider,
558 final UnnormalizedSphericalHarmonics harmonics,
559 final double M2Value) {
560 return computeMeanOrbit(osculating, provider, harmonics, M2Value, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
561 }
562
563 /** Conversion from osculating to mean orbit.
564 * <p>
565 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
566 * osculating SpacecraftState in input.
567 * </p>
568 * <p>
569 * Since the osculating orbit is obtained with the computation of
570 * short-periodic variation, the resulting output will depend on
571 * both the gravity field parameterized in input and the
572 * atmospheric drag represented by the {@code m2} parameter.
573 * </p>
574 * <p>
575 * The computation is done through a fixed-point iteration process.
576 * </p>
577 * @param <T> type of the filed elements
578 * @param osculating osculating orbit to convert
579 * @param provider for un-normalized zonal coefficients
580 * @param harmonics {@code provider.onDate(osculating.getDate())}
581 * @param M2Value value of empirical drag coefficient in rad/s².
582 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
583 * @param epsilon convergence threshold for mean parameters conversion
584 * @param maxIterations maximum iterations for mean parameters conversion
585 * @return mean orbit in a Brouwer-Lyddane sense
586 * @since 11.2
587 */
588 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
589 final UnnormalizedSphericalHarmonicsProvider provider,
590 final UnnormalizedSphericalHarmonics harmonics,
591 final double M2Value,
592 final double epsilon, final int maxIterations) {
593 return computeMeanOrbit(osculating,
594 provider.getAe(), provider.getMu(),
595 harmonics.getUnnormalizedCnm(2, 0),
596 harmonics.getUnnormalizedCnm(3, 0),
597 harmonics.getUnnormalizedCnm(4, 0),
598 harmonics.getUnnormalizedCnm(5, 0),
599 M2Value, epsilon, maxIterations);
600 }
601
602 /** Conversion from osculating to mean orbit.
603 * <p>
604 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
605 * osculating SpacecraftState in input.
606 * </p>
607 * <p>
608 * Since the osculating orbit is obtained with the computation of
609 * short-periodic variation, the resulting output will depend on
610 * both the gravity field parameterized in input and the
611 * atmospheric drag represented by the {@code m2} parameter.
612 * </p>
613 * <p>
614 * The computation is done through a fixed-point iteration process.
615 * </p>
616 * @param <T> type of the filed elements
617 * @param osculating osculating orbit to convert
618 * @param referenceRadius reference radius of the Earth for the potential model (m)
619 * @param mu central attraction coefficient (m³/s²)
620 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
621 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
622 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
623 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
624 * @param M2Value value of empirical drag coefficient in rad/s².
625 * If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
626 * @param epsilon convergence threshold for mean parameters conversion
627 * @param maxIterations maximum iterations for mean parameters conversion
628 * @return mean orbit in a Brouwer-Lyddane sense
629 * @since 11.2
630 */
631 public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
632 final double referenceRadius, final double mu,
633 final double c20, final double c30, final double c40,
634 final double c50, final double M2Value,
635 final double epsilon, final int maxIterations) {
636 final FieldBrouwerLyddanePropagator<T> propagator =
637 new FieldBrouwerLyddanePropagator<>(osculating,
638 FrameAlignedProvider.of(osculating.getFrame()),
639 osculating.getMu().newInstance(DEFAULT_MASS),
640 referenceRadius, osculating.getMu().newInstance(mu),
641 c20, c30, c40, c50,
642 PropagationType.OSCULATING, M2Value,
643 epsilon, maxIterations);
644 return propagator.initialModel.mean;
645 }
646
647 /** {@inheritDoc}
648 * <p>The new initial state to consider
649 * must be defined with an osculating orbit.</p>
650 * @see #resetInitialState(FieldSpacecraftState, PropagationType)
651 */
652 @Override
653 public void resetInitialState(final FieldSpacecraftState<T> state) {
654 resetInitialState(state, PropagationType.OSCULATING);
655 }
656
657 /** Reset the propagator initial state.
658 * @param state new initial state to consider
659 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
660 */
661 public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
662 resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
663 }
664
665 /** Reset the propagator initial state.
666 * @param state new initial state to consider
667 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
668 * @param epsilon convergence threshold for mean parameters conversion
669 * @param maxIterations maximum iterations for mean parameters conversion
670 * @since 11.2
671 */
672 public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType,
673 final double epsilon, final int maxIterations) {
674 super.resetInitialState(state);
675 final FieldKeplerianOrbit<T> keplerian = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit());
676 this.initialModel = (stateType == PropagationType.MEAN) ?
677 new FieldBLModel<>(keplerian, state.getMass(), referenceRadius, mu, ck0) :
678 computeMeanParameters(keplerian, state.getMass(), epsilon, maxIterations);
679 this.models = new FieldTimeSpanMap<>(initialModel, state.getA().getField());
680 }
681
682 /** {@inheritDoc} */
683 @Override
684 protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward) {
685 resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
686 }
687
688 /** Reset an intermediate state.
689 * @param state new intermediate state to consider
690 * @param forward if true, the intermediate state is valid for
691 * propagations after itself
692 * @param epsilon convergence threshold for mean parameters conversion
693 * @param maxIterations maximum iterations for mean parameters conversion
694 * @since 11.2
695 */
696 protected void resetIntermediateState(final FieldSpacecraftState<T> state, final boolean forward,
697 final double epsilon, final int maxIterations) {
698 final FieldBLModel<T> newModel = computeMeanParameters((FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(state.getOrbit()),
699 state.getMass(), epsilon, maxIterations);
700 if (forward) {
701 models.addValidAfter(newModel, state.getDate());
702 } else {
703 models.addValidBefore(newModel, state.getDate());
704 }
705 stateChanged(state);
706 }
707
708 /** Compute mean parameters according to the Brouwer-Lyddane analytical model computation
709 * in order to do the propagation.
710 * @param osculating osculating orbit
711 * @param mass constant mass
712 * @param epsilon convergence threshold for mean parameters conversion
713 * @param maxIterations maximum iterations for mean parameters conversion
714 * @return Brouwer-Lyddane mean model
715 */
716 private FieldBLModel<T> computeMeanParameters(final FieldKeplerianOrbit<T> osculating, final T mass,
717 final double epsilon, final int maxIterations) {
718
719 // sanity check
720 if (osculating.getA().getReal() < referenceRadius) {
721 throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
722 osculating.getA());
723 }
724
725 final Field<T> field = mass.getField();
726 final T one = field.getOne();
727 final T zero = field.getZero();
728 // rough initialization of the mean parameters
729 FieldBLModel<T> current = new FieldBLModel<>(osculating, mass, referenceRadius, mu, ck0);
730
731 // threshold for each parameter
732 final T thresholdA = current.mean.getA().abs().add(1.0).multiply(epsilon);
733 final T thresholdE = current.mean.getE().add(1.0).multiply(epsilon);
734 final T thresholdAngles = one.getPi().multiply(epsilon);
735
736 int i = 0;
737 while (i++ < maxIterations) {
738
739 // recompute the osculating parameters from the current mean parameters
740 final FieldKeplerianOrbit<T> parameters = current.propagateParameters(current.mean.getDate(), getParameters(mass.getField(), current.mean.getDate()));
741
742 // adapted parameters residuals
743 final T deltaA = osculating.getA() .subtract(parameters.getA());
744 final T deltaE = osculating.getE() .subtract(parameters.getE());
745 final T deltaI = osculating.getI() .subtract(parameters.getI());
746 final T deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument().subtract(parameters.getPerigeeArgument()), zero);
747 final T deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(parameters.getRightAscensionOfAscendingNode()), zero);
748 final T deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly().subtract(parameters.getMeanAnomaly()), zero);
749
750
751 // update mean parameters
752 current = new FieldBLModel<>(new FieldKeplerianOrbit<>(current.mean.getA() .add(deltaA),
753 FastMath.max(current.mean.getE().add(deltaE), zero),
754 current.mean.getI() .add(deltaI),
755 current.mean.getPerigeeArgument() .add(deltaOmega),
756 current.mean.getRightAscensionOfAscendingNode().add(deltaRAAN),
757 current.mean.getMeanAnomaly() .add(deltaAnom),
758 PositionAngleType.MEAN, PositionAngleType.MEAN,
759 current.mean.getFrame(),
760 current.mean.getDate(), mu),
761 mass, referenceRadius, mu, ck0);
762 // check convergence
763 if (FastMath.abs(deltaA.getReal()) < thresholdA.getReal() &&
764 FastMath.abs(deltaE.getReal()) < thresholdE.getReal() &&
765 FastMath.abs(deltaI.getReal()) < thresholdAngles.getReal() &&
766 FastMath.abs(deltaOmega.getReal()) < thresholdAngles.getReal() &&
767 FastMath.abs(deltaRAAN.getReal()) < thresholdAngles.getReal() &&
768 FastMath.abs(deltaAnom.getReal()) < thresholdAngles.getReal()) {
769 return current;
770 }
771 }
772 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_BROUWER_LYDDANE_MEAN_PARAMETERS, i);
773 }
774
775
776 /** {@inheritDoc} */
777 public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
778 // compute Cartesian parameters, taking derivatives into account
779 final FieldBLModel<T> current = models.get(date);
780 return current.propagateParameters(date, parameters);
781 }
782
783 /**
784 * Get the value of the M2 drag parameter.
785 * @return the value of the M2 drag parameter
786 */
787 public double getM2() {
788 return M2Driver.getValue();
789 }
790
791 /**
792 * Get the value of the M2 drag parameter.
793 * @param date date at which the model parameters want to be known
794 * @return the value of the M2 drag parameter
795 */
796 public double getM2(final AbsoluteDate date) {
797 return M2Driver.getValue(date);
798 }
799
800 /** Local class for Brouwer-Lyddane model. */
801 private static class FieldBLModel<T extends CalculusFieldElement<T>> {
802
803 /** Mean orbit. */
804 private final FieldKeplerianOrbit<T> mean;
805
806 /** Constant mass. */
807 private final T mass;
808
809 /** Central attraction coefficient. */
810 private final T mu;
811
812 // CHECKSTYLE: stop JavadocVariable check
813
814 // preprocessed values
815
816 // Constant for secular terms l'', g'', h''
817 // l standing for mean anomaly, g for perigee argument and h for raan
818 private final T xnotDot;
819 private final T n;
820 private final T lt;
821 private final T gt;
822 private final T ht;
823
824
825 // Long periodT
826 private final T dei3sg;
827 private final T de2sg;
828 private final T deisg;
829 private final T de;
830
831
832 private final T dlgs2g;
833 private final T dlgc3g;
834 private final T dlgcg;
835
836
837 private final T dh2sgcg;
838 private final T dhsgcg;
839 private final T dhcg;
840
841
842 // Short perioTs
843 private final T aC;
844 private final T aCbis;
845 private final T ac2g2f;
846
847
848 private final T eC;
849 private final T ecf;
850 private final T e2cf;
851 private final T e3cf;
852 private final T ec2f2g;
853 private final T ecfc2f2g;
854 private final T e2cfc2f2g;
855 private final T e3cfc2f2g;
856 private final T ec2gf;
857 private final T ec2g3f;
858
859
860 private final T ide;
861 private final T isfs2f2g;
862 private final T icfc2f2g;
863 private final T ic2f2g;
864
865
866 private final T glf;
867 private final T gll;
868 private final T glsf;
869 private final T glosf;
870 private final T gls2f2g;
871 private final T gls2gf;
872 private final T glos2gf;
873 private final T gls2g3f;
874 private final T glos2g3f;
875
876
877 private final T hf;
878 private final T hl;
879 private final T hsf;
880 private final T hcfs2g2f;
881 private final T hs2g2f;
882 private final T hsfc2g2f;
883
884
885 private final T edls2g;
886 private final T edlcg;
887 private final T edlc3g;
888 private final T edlsf;
889 private final T edls2gf;
890 private final T edls2g3f;
891
892 // Drag terms
893 private final T aRate;
894 private final T eRate;
895
896 // CHECKSTYLE: resume JavadocVariable check
897
898 /** Create a model for specified mean orbit.
899 * @param mean mean Fieldorbit
900 * @param mass constant mass
901 * @param referenceRadius reference radius of the central body attraction model (m)
902 * @param mu central attraction coefficient (m³/s²)
903 * @param ck0 un-normalized zonal coefficients
904 */
905 FieldBLModel(final FieldKeplerianOrbit<T> mean, final T mass,
906 final double referenceRadius, final T mu, final double[] ck0) {
907
908 this.mean = mean;
909 this.mass = mass;
910 this.mu = mu;
911 final T one = mass.getField().getOne();
912
913 final T app = mean.getA();
914 xnotDot = mu.divide(app).sqrt().divide(app);
915
916 // preliminary processing
917 final T q = app.divide(referenceRadius).reciprocal();
918 T ql = q.square();
919 final T y2 = ql.multiply(-0.5 * ck0[2]);
920
921 n = ((mean.getE().square().negate()).add(1.0)).sqrt();
922 final T n2 = n.square();
923 final T n3 = n2.multiply(n);
924 final T n4 = n2.square();
925 final T n6 = n4.multiply(n2);
926 final T n8 = n4.square();
927 final T n10 = n8.multiply(n2);
928
929 final T yp2 = y2.divide(n4);
930 ql = ql.multiply(q);
931 final T yp3 = ql.multiply(ck0[3]).divide(n6);
932 ql = ql.multiply(q);
933 final T yp4 = ql.multiply(0.375 * ck0[4]).divide(n8);
934 ql = ql.multiply(q);
935 final T yp5 = ql.multiply(ck0[5]).divide(n10);
936
937
938 final FieldSinCos<T> sc = FastMath.sinCos(mean.getI());
939 final T sinI1 = sc.sin();
940 final T sinI2 = sinI1.square();
941 final T cosI1 = sc.cos();
942 final T cosI2 = cosI1.square();
943 final T cosI3 = cosI2.multiply(cosI1);
944 final T cosI4 = cosI2.square();
945 final T cosI6 = cosI4.multiply(cosI2);
946 final T C5c2 = T2(cosI1).reciprocal();
947 final T C3c2 = cosI2.multiply(3.0).subtract(1.0);
948
949 final T epp = mean.getE();
950 final T epp2 = epp.square();
951 final T epp3 = epp2.multiply(epp);
952 final T epp4 = epp2.square();
953
954 if (epp.getReal() >= 1) {
955 // Only for elliptical (e < 1) orbits
956 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
957 mean.getE().getReal());
958 }
959
960 // secular multiplicative
961 lt = one.add(yp2.multiply(n).multiply(C3c2).multiply(1.5)).
962 add(yp2.multiply(0.09375).multiply(yp2).multiply(n).multiply(n2.multiply(25.0).add(n.multiply(16.0)).add(-15.0).
963 add(cosI2.multiply(n2.multiply(-90.0).add(n.multiply(-96.0)).add(30.0))).
964 add(cosI4.multiply(n2.multiply(25.0).add(n.multiply(144.0)).add(105.0))))).
965 add(yp4.multiply(0.9375).multiply(n).multiply(epp2).multiply(cosI4.multiply(35.0).add(cosI2.multiply(-30.0)).add(3.0)));
966
967 gt = yp2.multiply(-1.5).multiply(C5c2).
968 add(yp2.multiply(0.09375).multiply(yp2).multiply(n2.multiply(25.0).add(n.multiply(24.0)).add(-35.0).
969 add(cosI2.multiply(n2.multiply(-126.0).add(n.multiply(-192.0)).add(90.0))).
970 add(cosI4.multiply(n2.multiply(45.0).add(n.multiply(360.0)).add(385.0))))).
971 add(yp4.multiply(0.3125).multiply(n2.multiply(-9.0).add(21.0).
972 add(cosI2.multiply(n2.multiply(126.0).add(-270.0))).
973 add(cosI4.multiply(n2.multiply(-189.0).add(385.0)))));
974
975 ht = yp2.multiply(-3.0).multiply(cosI1).
976 add(yp2.multiply(0.375).multiply(yp2).multiply(cosI1.multiply(n2.multiply(9.0).add(n.multiply(12.0)).add(-5.0)).
977 add(cosI3.multiply(n2.multiply(-5.0).add(n.multiply(-36.0)).add(-35.0))))).
978 add(yp4.multiply(1.25).multiply(cosI1).multiply(n2.multiply(-3.0).add(5.0)).multiply(cosI2.multiply(-7.0).add(3.0)));
979
980 final T cA = one.subtract(cosI2.multiply(11.0)).subtract(cosI4.multiply(40.0).divide(C5c2));
981 final T cB = one.subtract(cosI2.multiply(3.0)).subtract(cosI4.multiply(8.0).divide(C5c2));
982 final T cC = one.subtract(cosI2.multiply(9.0)).subtract(cosI4.multiply(24.0).divide(C5c2));
983 final T cD = one.subtract(cosI2.multiply(5.0)).subtract(cosI4.multiply(16.0).divide(C5c2));
984
985 final T qyp2_4 = yp2.multiply(3.0).multiply(yp2).multiply(cA).
986 subtract(yp4.multiply(10.0).multiply(cB));
987 final T qyp52 = cosI1.multiply(epp3).multiply(cD.multiply(0.5).divide(sinI1).
988 add(sinI1.multiply(cosI4.divide(C5c2).divide(C5c2).multiply(80.0).
989 add(cosI2.divide(C5c2).multiply(32.0).
990 add(5.0)))));
991 final T qyp22 = epp2.add(2.0).subtract(epp2.multiply(3.0).add(2.0).multiply(11.0).multiply(cosI2)).
992 subtract(epp2.multiply(5.0).add(2.0).multiply(40.0).multiply(cosI4.divide(C5c2))).
993 subtract(epp2.multiply(400.0).multiply(cosI6).divide(C5c2).divide(C5c2));
994 final T qyp42 = one.divide(5.0).multiply(qyp22.
995 add(one.newInstance(4.0).multiply(epp2.
996 add(2.0).
997 subtract(cosI2.multiply(epp2.multiply(3.0).add(2.0))))));
998 final T qyp52bis = cosI1.multiply(sinI1).multiply(epp).multiply(epp2.multiply(3.0).add(4.0)).
999 multiply(cosI4.divide(C5c2).divide(C5c2).multiply(40.0).
1000 add(cosI2.divide(C5c2).multiply(16.0)).
1001 add(3.0));
1002 // long periodic multiplicative
1003 dei3sg = yp5.divide(yp2).multiply(35.0 / 96.0).multiply(epp2).multiply(n2).multiply(cD).multiply(sinI1);
1004 de2sg = qyp2_4.divide(yp2).multiply(epp).multiply(n2).multiply(-1.0 / 12.0);
1005 deisg = sinI1.multiply(yp5.multiply(-35.0 / 128.0).divide(yp2).multiply(epp2).multiply(n2).multiply(cD).
1006 add(n2.multiply(0.25).divide(yp2).multiply(yp3.add(yp5.multiply(5.0 / 16.0).multiply(epp2.multiply(3.0).add(4.0)).multiply(cC)))));
1007 de = epp2.multiply(n2).multiply(qyp2_4).divide(24.0).divide(yp2);
1008
1009 final T qyp52quotient = epp.multiply(epp4.multiply(81.0).add(-32.0)).divide(n.multiply(epp2.multiply(9.0).add(4.0)).add(epp2.multiply(3.0)).add(4.0));
1010 dlgs2g = yp2.multiply(48.0).reciprocal().multiply(yp2.multiply(-3.0).multiply(yp2).multiply(qyp22).
1011 add(yp4.multiply(10.0).multiply(qyp42))).
1012 add(n3.divide(yp2).multiply(qyp2_4).divide(24.0));
1013 dlgc3g = yp5.multiply(35.0 / 384.0).divide(yp2).multiply(n3).multiply(epp).multiply(cD).multiply(sinI1).
1014 add(yp5.multiply(35.0 / 1152.0).divide(yp2).multiply(qyp52.multiply(2.0).multiply(cosI1).subtract(epp.multiply(cD).multiply(sinI1).multiply(epp2.multiply(2.0).add(3.0)))));
1015 dlgcg = yp3.negate().multiply(cosI2).multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).
1016 add(yp5.divide(yp2).multiply(0.078125).multiply(cC).multiply(cosI2.divide(sinI1).multiply(epp.negate()).multiply(epp2.multiply(3.0).add(4.0)).
1017 add(sinI1.multiply(epp2).multiply(epp2.multiply(9.0).add(26.0)))).
1018 subtract(yp5.divide(yp2).multiply(0.46875).multiply(qyp52bis).multiply(cosI1)).
1019 add(yp3.divide(yp2).multiply(0.25).multiply(sinI1).multiply(epp).divide(n3.add(1.0)).multiply(epp2.multiply(epp2.subtract(3.0)).add(3.0))).
1020 add(yp5.divide(yp2).multiply(0.078125).multiply(n2).multiply(cC).multiply(qyp52quotient).multiply(sinI1)));
1021 final T qyp24 = yp2.multiply(yp2).multiply(3.0).multiply(cosI4.divide(sinI2).multiply(200.0).add(cosI2.divide(sinI1).multiply(80.0)).add(11.0)).
1022 subtract(yp4.multiply(10.0).multiply(cosI4.divide(sinI2).multiply(40.0).add(cosI2.divide(sinI1).multiply(16.0)).add(3.0)));
1023 dh2sgcg = yp5.divide(yp2).multiply(qyp52).multiply(35.0 / 144.0);
1024 dhsgcg = qyp24.multiply(cosI1).multiply(epp2.negate()).divide(yp2.multiply(12.0));
1025 dhcg = yp5.divide(yp2).multiply(qyp52).multiply(-35.0 / 576.0).
1026 add(cosI1.multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).multiply(yp3.add(yp5.multiply(0.3125).multiply(cC).multiply(epp2.multiply(3.0).add(4.0))))).
1027 add(yp5.multiply(qyp52bis).multiply(1.875).divide(yp2.multiply(4.0)));
1028
1029 // short periodic multiplicative
1030 aC = yp2.negate().multiply(C3c2).multiply(app).divide(n3);
1031 aCbis = y2.multiply(app).multiply(C3c2);
1032 ac2g2f = y2.multiply(app).multiply(sinI2).multiply(3.0);
1033
1034 T qe = y2.multiply(C3c2).multiply(0.5).multiply(n2).divide(n6);
1035 eC = qe.multiply(epp).divide(n3.add(1.0)).multiply(epp2.multiply(epp2.subtract(3.0)).add(3.0));
1036 ecf = qe.multiply(3.0);
1037 e2cf = qe.multiply(3.0).multiply(epp);
1038 e3cf = qe.multiply(epp2);
1039 qe = y2.multiply(0.5).multiply(n2).multiply(3.0).multiply(cosI2.negate().add(1.0)).divide(n6);
1040 ec2f2g = qe.multiply(epp);
1041 ecfc2f2g = qe.multiply(3.0);
1042 e2cfc2f2g = qe.multiply(3.0).multiply(epp);
1043 e3cfc2f2g = qe.multiply(epp2);
1044 qe = yp2.multiply(-0.5).multiply(n2).multiply(cosI2.negate().add(1.0));
1045 ec2gf = qe.multiply(3.0);
1046 ec2g3f = qe;
1047
1048 T qi = yp2.multiply(epp).multiply(cosI1).multiply(sinI1);
1049 ide = cosI1.multiply(epp.negate()).divide(sinI1.multiply(n2));
1050 isfs2f2g = qi;
1051 icfc2f2g = qi.multiply(2.0);
1052 qi = yp2.multiply(cosI1).multiply(sinI1);
1053 ic2f2g = qi.multiply(1.5);
1054
1055 T qgl1 = yp2.multiply(0.25);
1056 T qgl2 = yp2.multiply(epp).multiply(n2).multiply(0.25).divide(n.add(1.0));
1057 glf = qgl1.multiply(C5c2).multiply(-6.0);
1058 gll = qgl1.multiply(C5c2).multiply(6.0);
1059 glsf = qgl1.multiply(C5c2).multiply(-6.0).multiply(epp).
1060 add(qgl2.multiply(C3c2).multiply(2.0));
1061 glosf = qgl2.multiply(C3c2).multiply(2.0);
1062 qgl1 = qgl1.multiply(cosI2.multiply(-5.0).add(3.0));
1063 qgl2 = qgl2.multiply(3.0).multiply(cosI2.negate().add(1.0));
1064 gls2f2g = qgl1.multiply(3.0);
1065 gls2gf = qgl1.multiply(3.0).multiply(epp).
1066 add(qgl2);
1067 glos2gf = qgl2.negate();
1068 gls2g3f = qgl1.multiply(epp).
1069 add(qgl2.multiply(1.0 / 3.0));
1070 glos2g3f = qgl2;
1071
1072 final T qh = yp2.multiply(cosI1).multiply(3.0);
1073 hf = qh.negate();
1074 hl = qh;
1075 hsf = qh.multiply(epp).negate();
1076 hcfs2g2f = yp2.multiply(cosI1).multiply(epp).multiply(2.0);
1077 hs2g2f = yp2.multiply(cosI1).multiply(1.5);
1078 hsfc2g2f = yp2.multiply(cosI1).multiply(epp).negate();
1079
1080 final T qedl = yp2.multiply(n3).multiply(-0.25);
1081 edls2g = qyp2_4.multiply(1.0 / 24.0).multiply(epp).multiply(n3).divide(yp2);
1082 edlcg = yp3.divide(yp2).multiply(-0.25).multiply(n3).multiply(sinI1).
1083 subtract(yp5.divide(yp2).multiply(0.078125).multiply(n3).multiply(sinI1).multiply(cC).multiply(epp2.multiply(9.0).add(4.0)));
1084 edlc3g = yp5.divide(yp2).multiply(n3).multiply(epp2).multiply(cD).multiply(sinI1).multiply(35.0 / 384.0);
1085 edlsf = qedl.multiply(C3c2).multiply(2.0);
1086 edls2gf = qedl.multiply(3.0).multiply(cosI2.negate().add(1.0));
1087 edls2g3f = qedl.multiply(1.0 / 3.0);
1088
1089 // secular rates of the mean semi-major axis and eccentricity
1090 // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
1091 aRate = app.multiply(-4.0).divide(xnotDot.multiply(3.0));
1092 eRate = epp.multiply(n).multiply(n).multiply(-4.0).divide(xnotDot.multiply(3.0));
1093
1094 }
1095
1096 /**
1097 * Gets eccentric anomaly from mean anomaly.
1098 * @param mk the mean anomaly (rad)
1099 * @return the eccentric anomaly (rad)
1100 */
1101 private FieldUnivariateDerivative1<T> getEccentricAnomaly(final FieldUnivariateDerivative1<T> mk) {
1102
1103 final T zero = mean.getE().getField().getZero();
1104
1105 // reduce M to [-PI PI] interval
1106 final FieldUnivariateDerivative1<T> reducedM = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mk.getValue(), zero ),
1107 mk.getFirstDerivative());
1108
1109 final FieldUnivariateDerivative1<T> meanE = new FieldUnivariateDerivative1<>(mean.getE(), zero);
1110 FieldUnivariateDerivative1<T> ek = FieldKeplerianAnomalyUtility.ellipticMeanToEccentric(meanE, mk);
1111
1112 // expand the result back to original range
1113 ek = ek.add(mk.getValue().subtract(reducedM.getValue()));
1114
1115 // Returns the eccentric anomaly
1116 return ek;
1117 }
1118
1119 /**
1120 * This method is used in Brouwer-Lyddane model to avoid singularity at the
1121 * critical inclination (i = 63.4°).
1122 * <p>
1123 * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
1124 * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity)
1125 * by a function, named T2 in the thesis.
1126 * </p>
1127 * @param cosInc cosine of the mean inclination
1128 * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term
1129 */
1130 private T T2(final T cosInc) {
1131
1132 // X = (1.0 - 5.0 * cos²(inc))
1133 final T x = cosInc.square().multiply(-5.0).add(1.0);
1134 final T x2 = x.square();
1135
1136 // Eq. 2.48
1137 T sum = x.getField().getZero();
1138 for (int i = 0; i <= 12; i++) {
1139 final double sign = i % 2 == 0 ? +1.0 : -1.0;
1140 sum = sum.add(FastMath.pow(x2, i).multiply(FastMath.pow(BETA, i)).multiply(sign).divide(CombinatoricsUtils.factorialDouble(i + 1)));
1141 }
1142
1143 // Right term of equation 2.47
1144 T product = x.getField().getOne();
1145 for (int i = 0; i <= 10; i++) {
1146 product = product.multiply(FastMath.exp(x2.multiply(BETA).multiply(FastMath.scalb(-1.0, i))).add(1.0));
1147 }
1148
1149 // Return (Eq. 2.47)
1150 return x.multiply(BETA).multiply(sum).multiply(product);
1151
1152 }
1153
1154 /** Extrapolate an orbit up to a specific target date.
1155 * @param date target date for the orbit
1156 * @param parameters model parameters
1157 * @return propagated parameters
1158 */
1159 public FieldKeplerianOrbit<T> propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {
1160
1161 // Field
1162 final Field<T> field = date.getField();
1163 final T one = field.getOne();
1164 final T zero = field.getZero();
1165
1166 // Empirical drag coefficient M2
1167 final T m2 = parameters[0];
1168
1169 // Keplerian evolution
1170 final FieldUnivariateDerivative1<T> dt = new FieldUnivariateDerivative1<>(date.durationFrom(mean.getDate()), one);
1171 final FieldUnivariateDerivative1<T> xnot = dt.multiply(xnotDot);
1172
1173 //____________________________________
1174 // secular effects
1175
1176 // mean mean anomaly
1177 final FieldUnivariateDerivative1<T> dtM2 = dt.multiply(m2);
1178 final FieldUnivariateDerivative1<T> dt2M2 = dt.multiply(dtM2);
1179 final FieldUnivariateDerivative1<T> lpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getMeanAnomaly().add(lt.multiply(xnot.getValue())).add(dt2M2.getValue()), zero),
1180 lt.multiply(xnotDot).add(dtM2.multiply(2.0).getValue()));
1181 // mean argument of perigee
1182 final FieldUnivariateDerivative1<T> gpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getPerigeeArgument().add(gt.multiply(xnot.getValue())), zero),
1183 gt.multiply(xnotDot));
1184 // mean longitude of ascending node
1185 final FieldUnivariateDerivative1<T> hpp = new FieldUnivariateDerivative1<>(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode().add(ht.multiply(xnot.getValue())), zero),
1186 ht.multiply(xnotDot));
1187
1188 // ________________________________________________
1189 // secular rates of the mean semi-major axis and eccentricity
1190
1191 // semi-major axis
1192 final FieldUnivariateDerivative1<T> appDrag = dt.multiply(aRate.multiply(m2));
1193
1194 // eccentricity
1195 final FieldUnivariateDerivative1<T> eppDrag = dt.multiply(eRate.multiply(m2));
1196
1197 //____________________________________
1198 // Long periodical terms
1199 final FieldSinCos<FieldUnivariateDerivative1<T>> sinCosGpp = gpp.sinCos();
1200 final FieldUnivariateDerivative1<T> cg1 = sinCosGpp.cos();
1201 final FieldUnivariateDerivative1<T> sg1 = sinCosGpp.sin();
1202 final FieldUnivariateDerivative1<T> c2g = cg1.multiply(cg1).subtract(sg1.multiply(sg1));
1203 final FieldUnivariateDerivative1<T> s2g = cg1.multiply(sg1).add(sg1.multiply(cg1));
1204 final FieldUnivariateDerivative1<T> c3g = c2g.multiply(cg1).subtract(s2g.multiply(sg1));
1205 final FieldUnivariateDerivative1<T> sg2 = sg1.square();
1206 final FieldUnivariateDerivative1<T> sg3 = sg1.multiply(sg2);
1207
1208
1209 // de eccentricity
1210 final FieldUnivariateDerivative1<T> d1e = sg3.multiply(dei3sg).
1211 add(sg1.multiply(deisg)).
1212 add(sg2.multiply(de2sg)).
1213 add(de);
1214
1215 // l' + g'
1216 final FieldUnivariateDerivative1<T> lpPGp = s2g.multiply(dlgs2g).
1217 add(c3g.multiply(dlgc3g)).
1218 add(cg1.multiply(dlgcg)).
1219 add(lpp).
1220 add(gpp);
1221
1222 // h'
1223 final FieldUnivariateDerivative1<T> hp = sg2.multiply(cg1).multiply(dh2sgcg).
1224 add(sg1.multiply(cg1).multiply(dhsgcg)).
1225 add(cg1.multiply(dhcg)).
1226 add(hpp);
1227
1228 // Short periodical terms
1229 // eccentric anomaly
1230 final FieldUnivariateDerivative1<T> Ep = getEccentricAnomaly(lpp);
1231 final FieldSinCos<FieldUnivariateDerivative1<T>> sinCosEp = Ep.sinCos();
1232 final FieldUnivariateDerivative1<T> cf1 = (sinCosEp.cos().subtract(mean.getE())).divide(sinCosEp.cos().multiply(mean.getE().negate()).add(1.0));
1233 final FieldUnivariateDerivative1<T> sf1 = (sinCosEp.sin().multiply(n)).divide(sinCosEp.cos().multiply(mean.getE().negate()).add(1.0));
1234 final FieldUnivariateDerivative1<T> f = FastMath.atan2(sf1, cf1);
1235
1236 final FieldUnivariateDerivative1<T> cf2 = cf1.square();
1237 final FieldUnivariateDerivative1<T> c2f = cf2.subtract(sf1.multiply(sf1));
1238 final FieldUnivariateDerivative1<T> s2f = cf1.multiply(sf1).add(sf1.multiply(cf1));
1239 final FieldUnivariateDerivative1<T> c3f = c2f.multiply(cf1).subtract(s2f.multiply(sf1));
1240 final FieldUnivariateDerivative1<T> s3f = c2f.multiply(sf1).add(s2f.multiply(cf1));
1241 final FieldUnivariateDerivative1<T> cf3 = cf1.multiply(cf2);
1242
1243 final FieldUnivariateDerivative1<T> c2g1f = cf1.multiply(c2g).subtract(sf1.multiply(s2g));
1244 final FieldUnivariateDerivative1<T> c2g2f = c2f.multiply(c2g).subtract(s2f.multiply(s2g));
1245 final FieldUnivariateDerivative1<T> c2g3f = c3f.multiply(c2g).subtract(s3f.multiply(s2g));
1246 final FieldUnivariateDerivative1<T> s2g1f = cf1.multiply(s2g).add(c2g.multiply(sf1));
1247 final FieldUnivariateDerivative1<T> s2g2f = c2f.multiply(s2g).add(c2g.multiply(s2f));
1248 final FieldUnivariateDerivative1<T> s2g3f = c3f.multiply(s2g).add(c2g.multiply(s3f));
1249
1250 final FieldUnivariateDerivative1<T> eE = (sinCosEp.cos().multiply(mean.getE().negate()).add(1.0)).reciprocal();
1251 final FieldUnivariateDerivative1<T> eE3 = eE.square().multiply(eE);
1252 final FieldUnivariateDerivative1<T> sigma = eE.multiply(n.square()).multiply(eE).add(eE);
1253
1254 // Semi-major axis
1255 final FieldUnivariateDerivative1<T> a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())).
1256 add(aC).
1257 add(eE3.multiply(c2g2f).multiply(ac2g2f));
1258
1259 // Eccentricity
1260 final FieldUnivariateDerivative1<T> e = d1e.add(eppDrag.add(mean.getE())).
1261 add(eC).
1262 add(cf1.multiply(ecf)).
1263 add(cf2.multiply(e2cf)).
1264 add(cf3.multiply(e3cf)).
1265 add(c2g2f.multiply(ec2f2g)).
1266 add(c2g2f.multiply(cf1).multiply(ecfc2f2g)).
1267 add(c2g2f.multiply(cf2).multiply(e2cfc2f2g)).
1268 add(c2g2f.multiply(cf3).multiply(e3cfc2f2g)).
1269 add(c2g1f.multiply(ec2gf)).
1270 add(c2g3f.multiply(ec2g3f));
1271
1272 // Inclination
1273 final FieldUnivariateDerivative1<T> i = d1e.multiply(ide).
1274 add(mean.getI()).
1275 add(sf1.multiply(s2g2f.multiply(isfs2f2g))).
1276 add(cf1.multiply(c2g2f.multiply(icfc2f2g))).
1277 add(c2g2f.multiply(ic2f2g));
1278
1279 // Argument of perigee + mean anomaly
1280 final FieldUnivariateDerivative1<T> gPL = lpPGp.add(f.multiply(glf)).
1281 add(lpp.multiply(gll)).
1282 add(sf1.multiply(glsf)).
1283 add(sigma.multiply(sf1).multiply(glosf)).
1284 add(s2g2f.multiply(gls2f2g)).
1285 add(s2g1f.multiply(gls2gf)).
1286 add(sigma.multiply(s2g1f).multiply(glos2gf)).
1287 add(s2g3f.multiply(gls2g3f)).
1288 add(sigma.multiply(s2g3f).multiply(glos2g3f));
1289
1290 // Longitude of ascending node
1291 final FieldUnivariateDerivative1<T> h = hp.add(f.multiply(hf)).
1292 add(lpp.multiply(hl)).
1293 add(sf1.multiply(hsf)).
1294 add(cf1.multiply(s2g2f).multiply(hcfs2g2f)).
1295 add(s2g2f.multiply(hs2g2f)).
1296 add(c2g2f.multiply(sf1).multiply(hsfc2g2f));
1297
1298 final FieldUnivariateDerivative1<T> edl = s2g.multiply(edls2g).
1299 add(cg1.multiply(edlcg)).
1300 add(c3g.multiply(edlc3g)).
1301 add(sf1.multiply(edlsf)).
1302 add(s2g1f.multiply(edls2gf)).
1303 add(s2g3f.multiply(edls2g3f)).
1304 add(sf1.multiply(sigma).multiply(edlsf)).
1305 add(s2g1f.multiply(sigma).multiply(edls2gf.negate())).
1306 add(s2g3f.multiply(sigma).multiply(edls2g3f.multiply(3.0)));
1307
1308 final FieldUnivariateDerivative1<T> A = e.multiply(lpp.cos()).subtract(edl.multiply(lpp.sin()));
1309 final FieldUnivariateDerivative1<T> B = e.multiply(lpp.sin()).add(edl.multiply(lpp.cos()));
1310
1311 // Mean anomaly
1312 final FieldUnivariateDerivative1<T> l = FastMath.atan2(B, A);
1313
1314 // Argument of perigee
1315 final FieldUnivariateDerivative1<T> g = gPL.subtract(l);
1316
1317 // Return a Keplerian orbit
1318 return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
1319 g.getValue(), h.getValue(), l.getValue(),
1320 a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
1321 g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
1322 PositionAngleType.MEAN, mean.getFrame(), date, this.mu);
1323
1324 }
1325 }
1326
1327 /** {@inheritDoc} */
1328 @Override
1329 protected T getMass(final FieldAbsoluteDate<T> date) {
1330 return models.get(date).mass;
1331 }
1332
1333 /** {@inheritDoc} */
1334 @Override
1335 public List<ParameterDriver> getParametersDrivers() {
1336 return Collections.singletonList(M2Driver);
1337 }
1338
1339 }