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