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