1 /* Copyright 2002-2025 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.analytical;
18
19 import java.util.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.OrbitType;
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.utils.DoubleArrayDictionary;
52 import org.orekit.utils.ParameterDriver;
53 import org.orekit.utils.ParameterDriversProvider;
54 import org.orekit.utils.TimeSpanMap;
55 import org.orekit.utils.TimeSpanMap.Span;
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 * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics.
67 * However, for low Earth orbits, the magnitude of the perturbative acceleration due to
68 * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric
69 * drag by time derivatives of the <i>mean</i> mean anomaly using the catch-all coefficient
70 * {@link #M2Driver}. Beware that M2Driver must have only 1 span on its TimeSpanMap value.
71 *
72 * Usually, M2 is adjusted during an orbit determination process and it represents the
73 * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag).
74 * The behavior of M2 is close to the {@link TLE#getBStar()} parameter for the TLE.
75 *
76 * If the value of M2 is equal to {@link #M2 0.0}, the along-track secular effects are not
77 * considered in the dynamical model. Typical values for M2 are not known. It depends on the
78 * orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15).
79 * The unit of M2 is rad/s².
80 *
81 * The along-track effects, represented by the secular rates of the mean semi-major axis
82 * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis.
83 *
84 * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag.
85 * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959."
86 *
87 * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the
88 * artificial satellite. The Astronomical Journal 68 (1963): 555."
89 *
90 * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center
91 * (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992."
92 *
93 * @see "Solomon, Daniel, THE NAVSPASUR Satellite Motion Model,
94 * Naval Research Laboratory, August 8, 1991."
95 *
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 transient 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
565 // compute mean parameters if needed
566 resetInitialState(new SpacecraftState(initialOrbit,
567 attitudeProv.getAttitude(initialOrbit,
568 initialOrbit.getDate(),
569 initialOrbit.getFrame())).withMass(mass),
570 initialType, converter);
571
572 }
573
574 /**
575 * Private helper constructor.
576 * <p>Using this constructor, an initial osculating orbit is considered.</p>
577 * @param initialOrbit initial orbit
578 * @param attitude attitude provider
579 * @param mass spacecraft mass
580 * @param provider for un-normalized zonal coefficients
581 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
582 * @param m2Value value of empirical drag coefficient in rad/s².
583 * If equal to {@link #M2} drag is not computed
584 * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double,
585 * UnnormalizedSphericalHarmonicsProvider,
586 * UnnormalizedSphericalHarmonics,
587 * PropagationType, double)
588 */
589 private BrouwerLyddanePropagator(final Orbit initialOrbit,
590 final AttitudeProvider attitude,
591 final double mass,
592 final UnnormalizedSphericalHarmonicsProvider provider,
593 final UnnormalizedSphericalHarmonics harmonics,
594 final double m2Value) {
595 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
596 harmonics.getUnnormalizedCnm(2, 0),
597 harmonics.getUnnormalizedCnm(3, 0),
598 harmonics.getUnnormalizedCnm(4, 0),
599 harmonics.getUnnormalizedCnm(5, 0),
600 m2Value);
601 }
602
603 /**
604 * Private helper constructor.
605 * <p>Using this constructor, it is possible to define the initial orbit as
606 * a mean Brouwer-Lyddane orbit or an osculating one.</p>
607 * @param initialOrbit initial orbit
608 * @param attitude attitude provider
609 * @param mass spacecraft mass
610 * @param provider for un-normalized zonal coefficients
611 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
612 * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit)
613 * @param m2Value value of empirical drag coefficient in rad/s².
614 * If equal to {@link #M2} drag is not computed
615 */
616 private BrouwerLyddanePropagator(final Orbit initialOrbit,
617 final AttitudeProvider attitude,
618 final double mass,
619 final UnnormalizedSphericalHarmonicsProvider provider,
620 final UnnormalizedSphericalHarmonics harmonics,
621 final PropagationType initialType,
622 final double m2Value) {
623 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
624 harmonics.getUnnormalizedCnm(2, 0),
625 harmonics.getUnnormalizedCnm(3, 0),
626 harmonics.getUnnormalizedCnm(4, 0),
627 harmonics.getUnnormalizedCnm(5, 0),
628 initialType, m2Value);
629 }
630
631 /** Conversion from osculating to mean orbit.
632 * <p>
633 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
634 * osculating SpacecraftState in input.
635 * </p>
636 * <p>
637 * Since the osculating orbit is obtained with the computation of
638 * short-periodic variation, the resulting output will depend on
639 * both the gravity field parameterized in input and the
640 * atmospheric drag represented by the {@code m2Value} parameter.
641 * </p>
642 * <p>
643 * The computation is done through a fixed-point iteration process.
644 * </p>
645 * @param osculating osculating orbit to convert
646 * @param provider for un-normalized zonal coefficients
647 * @param harmonics {@code provider.onDate(osculating.getDate())}
648 * @param m2Value value of empirical drag coefficient in rad/s².
649 * If equal to {@link #M2} drag is not considered
650 * @return mean orbit in a Brouwer-Lyddane sense
651 * @since 11.2
652 */
653 public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
654 final UnnormalizedSphericalHarmonicsProvider provider,
655 final UnnormalizedSphericalHarmonics harmonics,
656 final double m2Value) {
657 return computeMeanOrbit(osculating, provider, harmonics, m2Value,
658 EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
659 }
660
661 /** Conversion from osculating to mean orbit.
662 * <p>
663 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
664 * osculating SpacecraftState in input.
665 * </p>
666 * <p>
667 * Since the osculating orbit is obtained with the computation of
668 * short-periodic variation, the resulting output will depend on
669 * both the gravity field parameterized in input and the
670 * atmospheric drag represented by the {@code m2Value} parameter.
671 * </p>
672 * <p>
673 * The computation is done through a fixed-point iteration process.
674 * </p>
675 * @param osculating osculating orbit to convert
676 * @param provider for un-normalized zonal coefficients
677 * @param harmonics {@code provider.onDate(osculating.getDate())}
678 * @param m2Value value of empirical drag coefficient in rad/s².
679 * If equal to {@link #M2} drag is not considered
680 * @param epsilon convergence threshold for mean parameters conversion
681 * @param maxIterations maximum iterations for mean parameters conversion
682 * @return mean orbit in a Brouwer-Lyddane sense
683 * @since 11.2
684 */
685 public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
686 final UnnormalizedSphericalHarmonicsProvider provider,
687 final UnnormalizedSphericalHarmonics harmonics,
688 final double m2Value,
689 final double epsilon,
690 final int maxIterations) {
691 return computeMeanOrbit(osculating,
692 provider.getAe(), provider.getMu(),
693 harmonics.getUnnormalizedCnm(2, 0),
694 harmonics.getUnnormalizedCnm(3, 0),
695 harmonics.getUnnormalizedCnm(4, 0),
696 harmonics.getUnnormalizedCnm(5, 0),
697 m2Value, epsilon, maxIterations);
698 }
699
700 /** Conversion from osculating to mean orbit.
701 * <p>
702 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
703 * osculating SpacecraftState in input.
704 * </p>
705 * <p>
706 * Since the osculating orbit is obtained with the computation of
707 * short-periodic variation, the resulting output will depend on
708 * both the gravity field parameterized in input and the
709 * atmospheric drag represented by the {@code m2Value} parameter.
710 * </p>
711 * <p>
712 * The computation is done through a fixed-point iteration process.
713 * </p>
714 * @param osculating osculating orbit to convert
715 * @param referenceRadius reference radius of the Earth for the potential model (m)
716 * @param mu central attraction coefficient (m³/s²)
717 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
718 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
719 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
720 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
721 * @param m2Value value of empirical drag coefficient in rad/s².
722 * If equal to {@link #M2} drag is not considered
723 * @param epsilon convergence threshold for mean parameters conversion
724 * @param maxIterations maximum iterations for mean parameters conversion
725 * @return mean orbit in a Brouwer-Lyddane sense
726 * @since 11.2
727 */
728 public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
729 final double referenceRadius,
730 final double mu,
731 final double c20,
732 final double c30,
733 final double c40,
734 final double c50,
735 final double m2Value,
736 final double epsilon,
737 final int maxIterations) {
738 // Build a fixed-point converter
739 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
740 FixedPointConverter.DEFAULT_DAMPING);
741 return computeMeanOrbit(osculating, referenceRadius, mu, c20, c30, c40, c50, m2Value, converter);
742 }
743
744 /** Conversion from osculating to mean orbit.
745 * <p>
746 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
747 * osculating SpacecraftState in input.
748 * </p>
749 * <p>
750 * Since the osculating orbit is obtained with the computation of
751 * short-periodic variation, the resulting output will depend on
752 * both the gravity field parameterized in input and the
753 * atmospheric drag represented by the {@code m2Value} parameter.
754 * </p>
755 * <p>
756 * The computation is done through the given osculating to mean orbit converter.
757 * </p>
758 * @param osculating osculating orbit to convert
759 * @param referenceRadius reference radius of the Earth for the potential model (m)
760 * @param mu central attraction coefficient (m³/s²)
761 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
762 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
763 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
764 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
765 * @param m2Value value of empirical drag coefficient in rad/s².
766 * If equal to {@link #M2} drag is not considered
767 * @param converter osculating to mean orbit converter
768 * @return mean orbit in a Brouwer-Lyddane sense
769 * @since 13.0
770 */
771 public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
772 final double referenceRadius,
773 final double mu,
774 final double c20,
775 final double c30,
776 final double c40,
777 final double c50,
778 final double m2Value,
779 final OsculatingToMeanConverter converter) {
780 // Set BL as the mean theory for converting
781 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu, c20, c30, c40, c50, m2Value);
782 converter.setMeanTheory(theory);
783 return (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converter.convertToMean(osculating));
784 }
785
786 /** Conversion from osculating to mean orbit.
787 * <p>
788 * Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
789 * osculating SpacecraftState in input.
790 * </p>
791 * <p>
792 * Since the osculating orbit is obtained with the computation of
793 * short-periodic variation, the resulting output will depend on
794 * both the gravity field parameterized in input and the
795 * atmospheric drag represented by the {@code m2Value} parameter.
796 * </p>
797 * <p>
798 * The computation is done through the given osculating to mean orbit converter.
799 * </p>
800 * @param osculating osculating orbit to convert
801 * @param provider for un-normalized zonal coefficients
802 * @param m2Value value of empirical drag coefficient in rad/s².
803 * If equal to {@link #M2} drag is not considered
804 * @param converter osculating to mean orbit converter
805 * @return mean orbit in a Brouwer-Lyddane sense
806 * @since 13.0
807 */
808 public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
809 final UnnormalizedSphericalHarmonicsProvider provider,
810 final double m2Value,
811 final OsculatingToMeanConverter converter) {
812 // Set BL as the mean theory for converting
813 final MeanTheory theory = new BrouwerLyddaneTheory(provider, m2Value);
814 converter.setMeanTheory(theory);
815 return (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converter.convertToMean(osculating));
816 }
817
818 /** {@inheritDoc}
819 * <p>The new initial state to consider
820 * must be defined with an osculating orbit.</p>
821 * @see #resetInitialState(SpacecraftState, PropagationType)
822 */
823 @Override
824 public void resetInitialState(final SpacecraftState state) {
825 resetInitialState(state, PropagationType.OSCULATING);
826 }
827
828 /** Reset the propagator initial state.
829 * @param state new initial state to consider
830 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
831 */
832 public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
833 resetInitialState(state, stateType, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
834 }
835
836 /** Reset the propagator initial state.
837 * @param state new initial state to consider
838 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
839 * @param epsilon convergence threshold for mean parameters conversion
840 * @param maxIterations maximum iterations for mean parameters conversion
841 * @since 11.2
842 */
843 public void resetInitialState(final SpacecraftState state,
844 final PropagationType stateType,
845 final double epsilon,
846 final int maxIterations) {
847 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
848 FixedPointConverter.DEFAULT_DAMPING);
849 resetInitialState(state, stateType, converter);
850 }
851
852 /** Reset the propagator initial state.
853 * @param state new initial state to consider
854 * @param stateType mean Brouwer-Lyddane orbit or osculating orbit
855 * @param converter osculating to mean orbit converter
856 * @since 13.0
857 */
858 public void resetInitialState(final SpacecraftState state,
859 final PropagationType stateType,
860 final OsculatingToMeanConverter converter) {
861 super.resetInitialState(state);
862 KeplerianOrbit keplerian = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit());
863 if (stateType == PropagationType.OSCULATING) {
864 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu,
865 ck0[2], ck0[3], ck0[4], ck0[5],
866 getM2());
867 converter.setMeanTheory(theory);
868 keplerian = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converter.convertToMean(keplerian));
869 }
870 this.initialModel = new BLModel(keplerian, state.getMass(), referenceRadius, mu, ck0);
871 this.models = new TimeSpanMap<>(initialModel);
872 }
873
874 /** {@inheritDoc} */
875 protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
876 resetIntermediateState(state, forward, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
877 }
878
879 /** Reset an intermediate state.
880 * @param state new intermediate state to consider
881 * @param forward if true, the intermediate state is valid for
882 * propagations after itself
883 * @param epsilon convergence threshold for mean parameters conversion
884 * @param maxIterations maximum iterations for mean parameters conversion
885 * @since 11.2
886 */
887 protected void resetIntermediateState(final SpacecraftState state,
888 final boolean forward,
889 final double epsilon,
890 final int maxIterations) {
891 final OsculatingToMeanConverter converter = new FixedPointConverter(epsilon, maxIterations,
892 FixedPointConverter.DEFAULT_DAMPING);
893 resetIntermediateState(state, forward, converter);
894 }
895
896 /** Reset an intermediate state.
897 * @param state new intermediate state to consider
898 * @param forward if true, the intermediate state is valid for
899 * propagations after itself
900 * @param converter osculating to mean orbit converter
901 * @since 13.0
902 */
903 protected void resetIntermediateState(final SpacecraftState state,
904 final boolean forward,
905 final OsculatingToMeanConverter converter) {
906 final MeanTheory theory = new BrouwerLyddaneTheory(referenceRadius, mu,
907 ck0[2], ck0[3], ck0[4], ck0[5],
908 getM2());
909 converter.setMeanTheory(theory);
910 final KeplerianOrbit mean = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(converter.convertToMean(state.getOrbit()));
911 final BLModel newModel = new BLModel(mean, state.getMass(), referenceRadius, mu, ck0);
912 if (forward) {
913 models.addValidAfter(newModel, state.getDate(), false);
914 } else {
915 models.addValidBefore(newModel, state.getDate(), false);
916 }
917 stateChanged(state);
918 }
919
920 /** {@inheritDoc} */
921 public KeplerianOrbit propagateOrbit(final AbsoluteDate date) {
922 // compute Keplerian parameters, taking derivatives into account
923 final BLModel current = models.get(date);
924 return current.propagateParameters(date);
925 }
926
927 /**
928 * Get the value of the M2 drag parameter. Beware that M2Driver
929 * must have only 1 span on its TimeSpanMap value (that is
930 * to say setPeriod method should not be called)
931 * @return the value of the M2 drag parameter
932 */
933 public double getM2() {
934 // As Brouwer Lyddane is an analytical propagator, for now it is not possible for
935 // M2Driver to have several values estimated
936 return M2Driver.getValue();
937 }
938
939 /**
940 * Get the central attraction coefficient μ.
941 * @return mu central attraction coefficient (m³/s²)
942 */
943 public double getMu() {
944 return mu;
945 }
946
947 /**
948 * Get the un-normalized zonal coefficients.
949 * @return the un-normalized zonal coefficients
950 */
951 public double[] getCk0() {
952 return ck0.clone();
953 }
954
955 /**
956 * Get the reference radius of the central body attraction model.
957 * @return the reference radius in meters
958 */
959 public double getReferenceRadius() {
960 return referenceRadius;
961 }
962
963 /**
964 * Get the parameters driver for propagation model.
965 * @return drivers for propagation model
966 */
967 public List<ParameterDriver> getParametersDrivers() {
968 return Collections.singletonList(M2Driver);
969 }
970
971 /** {@inheritDoc} */
972 @Override
973 protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
974 final DoubleArrayDictionary initialJacobianColumns) {
975 // Create the harvester
976 final BrouwerLyddaneHarvester harvester = new BrouwerLyddaneHarvester(this, stmName, initialStm, initialJacobianColumns);
977 // Update the list of additional state provider
978 addAdditionalDataProvider(harvester);
979 // Return the configured harvester
980 return harvester;
981 }
982
983 /**
984 * Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
985 * @return names of the parameters (i.e. columns) of the Jacobian matrix
986 */
987 @Override
988 protected List<String> getJacobiansColumnsNames() {
989 final List<String> columnsNames = new ArrayList<>();
990 for (final ParameterDriver driver : getParametersDrivers()) {
991 if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
992 for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
993 columnsNames.add(span.getData());
994 }
995 }
996 }
997 Collections.sort(columnsNames);
998 return columnsNames;
999 }
1000
1001 /** Local class for Brouwer-Lyddane model. */
1002 private class BLModel {
1003
1004 /** Constant mass. */
1005 private final double mass;
1006
1007 /** Brouwer-Lyddane mean orbit. */
1008 private final KeplerianOrbit mean;
1009
1010 // Preprocessed values
1011
1012 /** Mean mean motion: n0 = √(μ/a")/a". */
1013 private final double n0;
1014
1015 /** η = √(1 - e"²). */
1016 private final double n;
1017 /** η². */
1018 private final double n2;
1019 /** η³. */
1020 private final double n3;
1021 /** η + 1 / (1 + η). */
1022 private final double t8;
1023
1024 /** Secular correction for mean anomaly l: δ<sub>s</sub>l. */
1025 private final double dsl;
1026 /** Secular correction for perigee argument g: δ<sub>s</sub>g. */
1027 private final double dsg;
1028 /** Secular correction for raan h: δ<sub>s</sub>h. */
1029 private final double dsh;
1030
1031 /** Secular rate of change of semi-major axis due to drag. */
1032 private final double aRate;
1033 /** Secular rate of change of eccentricity due to drag. */
1034 private final double eRate;
1035
1036 // CHECKSTYLE: stop JavadocVariable check
1037
1038 // Storage for speed-up
1039 private final double yp2;
1040 private final double ci;
1041 private final double si;
1042 private final double oneMci2;
1043 private final double ci2X3M1;
1044
1045 // Long periodic corrections factors
1046 private final double vle1;
1047 private final double vle2;
1048 private final double vle3;
1049 private final double vli1;
1050 private final double vli2;
1051 private final double vli3;
1052 private final double vll2;
1053 private final double vlh1I;
1054 private final double vlh2I;
1055 private final double vlh3I;
1056 private final double vls1;
1057 private final double vls2;
1058 private final double vls3;
1059
1060 // CHECKSTYLE: resume JavadocVariable check
1061
1062 /** Create a model for specified mean orbit.
1063 * @param mean mean orbit
1064 * @param mass constant mass
1065 * @param referenceRadius reference radius of the central body attraction model (m)
1066 * @param mu central attraction coefficient (m³/s²)
1067 * @param ck0 un-normalized zonal coefficients
1068 */
1069 BLModel(final KeplerianOrbit mean, final double mass,
1070 final double referenceRadius, final double mu, final double[] ck0) {
1071
1072 this.mass = mass;
1073
1074 // mean orbit
1075 this.mean = mean;
1076
1077 // mean eccentricity e"
1078 final double epp = mean.getE();
1079 if (epp >= 1) {
1080 // Only for elliptical (e < 1) orbits
1081 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
1082 epp);
1083 }
1084 final double epp2 = epp * epp;
1085
1086 // η
1087 n2 = 1. - epp2;
1088 n = FastMath.sqrt(n2);
1089 n3 = n2 * n;
1090 t8 = n + 1. / (1. + n);
1091
1092 // mean semi-major axis a"
1093 final double app = mean.getA();
1094
1095 // mean mean motion
1096 n0 = FastMath.sqrt(mu / app) / app;
1097
1098 // ae/a"
1099 final double q = referenceRadius / app;
1100
1101 // γ2'
1102 double ql = q * q;
1103 double nl = n2 * n2;
1104 yp2 = -0.5 * ck0[2] * ql / nl;
1105 final double yp22 = yp2 * yp2;
1106
1107 // γ3'
1108 ql *= q;
1109 nl *= n2;
1110 final double yp3 = ck0[3] * ql / nl;
1111
1112 // γ4'
1113 ql *= q;
1114 nl *= n2;
1115 final double yp4 = 0.375 * ck0[4] * ql / nl;
1116
1117 // γ5'
1118 ql *= q;
1119 nl *= n2;
1120 final double yp5 = ck0[5] * ql / nl;
1121
1122 // mean inclination I" sin & cos
1123 final SinCos sci = FastMath.sinCos(mean.getI());
1124 si = sci.sin();
1125 ci = sci.cos();
1126 final double ci2 = ci * ci;
1127 oneMci2 = 1.0 - ci2;
1128 ci2X3M1 = 3.0 * ci2 - 1.0;
1129 final double ci2X5M1 = 5.0 * ci2 - 1.0;
1130
1131 // secular corrections
1132 // true anomaly
1133 dsl = 1.5 * yp2 * n * (ci2X3M1 +
1134 0.0625 * yp2 * (-15.0 + n * (16.0 + 25.0 * n) +
1135 ci2 * (30.0 - n * (96.0 + 90.0 * n) +
1136 ci2 * (105.0 + n * (144.0 + 25.0 * n))))) +
1137 0.9375 * yp4 * n * epp2 * (3.0 - ci2 * (30.0 - 35.0 * ci2));
1138 // perigee argument
1139 dsg = 1.5 * yp2 * ci2X5M1 +
1140 0.09375 * yp22 * (-35.0 + n * (24.0 + 25.0 * n) +
1141 ci2 * (90.0 - n * (192.0 + 126.0 * n) +
1142 ci2 * (385.0 + n * (360.0 + 45.0 * n)))) +
1143 0.3125 * yp4 * (21.0 - 9.0 * n2 + ci2 * (-270.0 + 126.0 * n2 +
1144 ci2 * (385.0 - 189.0 * n2)));
1145 // right ascension of ascending node
1146 dsh = (-3.0 * yp2 +
1147 0.375 * yp22 * (-5.0 + n * (12.0 + 9.0 * n) -
1148 ci2 * (35.0 + n * (36.0 + 5.0 * n))) +
1149 1.25 * yp4 * (5.0 - 3.0 * n2) * (3.0 - 7.0 * ci2)) * ci;
1150
1151 // secular rates of change due to drag
1152 // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis
1153 final double coef = -4.0 / (3.0 * n0 * (1 + dsl));
1154 aRate = coef * app;
1155 eRate = coef * epp * n2;
1156
1157 // singular term 1/(1 - 5 * cos²(I")) replaced by T2 function
1158 final double t2 = T2(ci);
1159
1160 // factors for long periodic corrections
1161 final double fs12 = yp3 / yp2;
1162 final double fs13 = 10. * yp4 / (3. * yp2);
1163 final double fs14 = yp5 / yp2;
1164
1165 final double ci2Xt2 = ci2 * t2;
1166 final double cA = 1. - ci2 * (11. + 40. * ci2Xt2);
1167 final double cB = 1. - ci2 * ( 3. + 8. * ci2Xt2);
1168 final double cC = 1. - ci2 * ( 9. + 24. * ci2Xt2);
1169 final double cD = 1. - ci2 * ( 5. + 16. * ci2Xt2);
1170 final double cE = 1. - ci2 * (33. + 200. * ci2Xt2);
1171 final double cF = 1. - ci2 * ( 9. + 40. * ci2Xt2);
1172
1173 final double p5p = 1. + ci2Xt2 * (8. + 20 * ci2Xt2);
1174 final double p5p2 = 1. + 2. * p5p;
1175 final double p5p4 = 1. + 4. * p5p;
1176 final double p5p10 = 1. + 10. * p5p;
1177
1178 final double e2X3P4 = 4. + 3. * epp2;
1179 final double ciO1Pci = ci / (1. + ci);
1180
1181 final double q1 = 0.125 * (yp2 * cA - fs13 * cB);
1182 final double q2 = 0.125 * epp2 * ci * (yp2 * p5p10 - fs13 * p5p2);
1183 final double q5 = 0.25 * (fs12 + 0.3125 * e2X3P4 * fs14 * cC);
1184 final double p2 = 0.46875 * p5p2 * epp * ci * si * e2X3P4 * fs14;
1185 final double p3 = 0.15625 * epp * si * fs14 * cC;
1186 final double kf = 35. / 1152.;
1187 final double p4 = kf * epp * fs14 * cD;
1188 final double p5 = 2. * kf * epp * epp2 * ci * si * fs14 * p5p4;
1189
1190 vle1 = epp * n2 * q1;
1191 vle2 = n2 * si * q5;
1192 vle3 = -3.0 * epp * n2 * si * p4;
1193
1194 vli1 = -epp * q1 / si;
1195 vli2 = -epp * ci * q5;
1196 vli3 = -3.0 * epp2 * ci * p4;
1197
1198 vll2 = vle2 + 3.0 * epp * n2 * p3;
1199
1200 vlh1I = -si * q2;
1201 vlh2I = epp * ci * q5 + si * p2;
1202 vlh3I = -epp2 * ci * p4 - si * p5;
1203
1204 vls1 = (n3 - 1.0) * q1 -
1205 q2 +
1206 25.0 * epp2 * ci2 * ci2Xt2 * ci2Xt2 * (yp2 - 0.2 * fs13) -
1207 0.0625 * epp2 * (yp2 * cE - fs13 * cF);
1208
1209 vls2 = epp * si * (t8 + ciO1Pci) * q5 +
1210 (11.0 + 3.0 * (epp2 - n3)) * p3 +
1211 (1.0 - ci) * p2;
1212
1213 vls3 = si * p4 * (3.0 * (n3 - 1.0) - epp2 * (2.0 + ciO1Pci)) -
1214 (1.0 - ci) * p5;
1215 }
1216
1217 /**
1218 * Get true anomaly from mean anomaly.
1219 * @param lM the mean anomaly (rad)
1220 * @param ecc the eccentricity
1221 * @return the true anomaly (rad)
1222 */
1223 private UnivariateDerivative1 getTrueAnomaly(final UnivariateDerivative1 lM,
1224 final UnivariateDerivative1 ecc) {
1225 // reduce M to [-PI PI] interval
1226 final double reducedM = MathUtils.normalizeAngle(lM.getValue(), 0.);
1227
1228 // compute the true anomaly
1229 UnivariateDerivative1 lV = FieldKeplerianAnomalyUtility.ellipticMeanToTrue(ecc, lM);
1230
1231 // expand the result back to original range
1232 lV = lV.add(lM.getValue() - reducedM);
1233
1234 // Returns the true anomaly
1235 return lV;
1236 }
1237
1238 /**
1239 * This method is used in Brouwer-Lyddane model to avoid singularity at the
1240 * critical inclination (i = 63.4°).
1241 * <p>
1242 * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48),
1243 * approximate the factor (1.0 - 5.0 * cos²(i))<sup>-1</sup> (causing the singularity)
1244 * by a function, named T2 in the thesis.
1245 * </p>
1246 * @param cosI cosine of the mean inclination
1247 * @return an approximation of (1.0 - 5.0 * cos²(i))<sup>-1</sup> term
1248 */
1249 private double T2(final double cosI) {
1250
1251 // X = 1.0 - 5.0 * cos²(i)
1252 final double x = 1.0 - 5.0 * cosI * cosI;
1253 final double x2 = x * x;
1254
1255 // Eq. 2.48
1256 double sum = 0.0;
1257 for (int i = 0; i <= 12; i++) {
1258 final double sign = i % 2 == 0 ? +1.0 : -1.0;
1259 sum += sign * FastMath.pow(BETA, i) * FastMath.pow(x2, i) / CombinatoricsUtils.factorialDouble(i + 1);
1260 }
1261
1262 // Right term of equation 2.47
1263 double product = 1.0;
1264 for (int i = 0; i <= 10; i++) {
1265 product *= 1 + FastMath.exp(FastMath.scalb(-1.0, i) * BETA * x2);
1266 }
1267
1268 // Return (Eq. 2.47)
1269 return BETA * x * sum * product;
1270 }
1271
1272 /** Extrapolate an orbit up to a specific target date.
1273 * @param date target date for the orbit
1274 * @return propagated parameters
1275 */
1276 public KeplerianOrbit propagateParameters(final AbsoluteDate date) {
1277
1278 // Empirical drag coefficient M2
1279 final double m2 = getM2();
1280
1281 // Keplerian evolution
1282 final UnivariateDerivative1 dt = new UnivariateDerivative1(date.durationFrom(mean.getDate()), 1.0);
1283 final UnivariateDerivative1 not = dt.multiply(n0);
1284
1285 final UnivariateDerivative1 dtM2 = dt.multiply(m2);
1286 final UnivariateDerivative1 dt2M2 = dt.multiply(dtM2);
1287
1288 // Secular corrections
1289 // -------------------
1290
1291 // semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis)
1292 final UnivariateDerivative1 app = dtM2.multiply(aRate).add(mean.getA());
1293
1294 // eccentricity (with drag Eq. 2.45 of Phipps' 1992 thesis) reduced to [0, 1[
1295 final UnivariateDerivative1 tmp = dtM2.multiply(eRate).add(mean.getE());
1296 final UnivariateDerivative1 epp = tmp.withValue(FastMath.max(0., FastMath.min(tmp.getValue(), MAX_ECC)));
1297
1298 // argument of perigee
1299 final double gppVal = mean.getPerigeeArgument() + dsg * not.getValue();
1300 final UnivariateDerivative1 gpp = new UnivariateDerivative1(MathUtils.normalizeAngle(gppVal, 0.),
1301 dsg * n0);
1302
1303 // longitude of ascending node
1304 final double hppVal = mean.getRightAscensionOfAscendingNode() + dsh * not.getValue();
1305 final UnivariateDerivative1 hpp = new UnivariateDerivative1(MathUtils.normalizeAngle(hppVal, 0.),
1306 dsh * n0);
1307
1308 // mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis)
1309 final double lppVal = mean.getMeanAnomaly() + (1. + dsl) * not.getValue() + dt2M2.getValue();
1310 final double dlppdt = (1. + dsl) * n0 + 2.0 * dtM2.getValue();
1311 final UnivariateDerivative1 lpp = new UnivariateDerivative1(MathUtils.normalizeAngle(lppVal, 0.),
1312 dlppdt);
1313
1314 // Long period corrections
1315 //------------------------
1316 final FieldSinCos<UnivariateDerivative1> scgpp = gpp.sinCos();
1317 final UnivariateDerivative1 cgpp = scgpp.cos();
1318 final UnivariateDerivative1 sgpp = scgpp.sin();
1319 final FieldSinCos<UnivariateDerivative1> sc2gpp = gpp.multiply(2).sinCos();
1320 final UnivariateDerivative1 c2gpp = sc2gpp.cos();
1321 final UnivariateDerivative1 s2gpp = sc2gpp.sin();
1322 final FieldSinCos<UnivariateDerivative1> sc3gpp = gpp.multiply(3).sinCos();
1323 final UnivariateDerivative1 c3gpp = sc3gpp.cos();
1324 final UnivariateDerivative1 s3gpp = sc3gpp.sin();
1325
1326 // δ1e
1327 final UnivariateDerivative1 d1e = c2gpp.multiply(vle1).
1328 add(sgpp.multiply(vle2)).
1329 add(s3gpp.multiply(vle3));
1330
1331 // δ1I
1332 UnivariateDerivative1 d1I = sgpp.multiply(vli2).
1333 add(s3gpp.multiply(vli3));
1334 // Pseudo singular term, not to add if Ipp is zero
1335 if (Double.isFinite(vli1)) {
1336 d1I = d1I.add(c2gpp.multiply(vli1));
1337 }
1338
1339 // e"δ1l
1340 final UnivariateDerivative1 eppd1l = s2gpp.multiply(vle1).
1341 subtract(cgpp.multiply(vll2)).
1342 subtract(c3gpp.multiply(vle3)).
1343 multiply(n);
1344
1345 // δ1h
1346 final UnivariateDerivative1 sIppd1h = s2gpp.multiply(vlh1I).
1347 add(cgpp.multiply(vlh2I)).
1348 add(c3gpp.multiply(vlh3I));
1349
1350 // δ1z = δ1l + δ1g + δ1h
1351 final UnivariateDerivative1 d1z = s2gpp.multiply(vls1).
1352 add(cgpp.multiply(vls2)).
1353 add(c3gpp.multiply(vls3));
1354
1355 // Short period corrections
1356 // ------------------------
1357
1358 // true anomaly
1359 final UnivariateDerivative1 fpp = getTrueAnomaly(lpp, epp);
1360 final FieldSinCos<UnivariateDerivative1> scfpp = fpp.sinCos();
1361 final UnivariateDerivative1 cfpp = scfpp.cos();
1362 final UnivariateDerivative1 sfpp = scfpp.sin();
1363
1364 // e"sin(f')
1365 final UnivariateDerivative1 eppsfpp = epp.multiply(sfpp);
1366 // e"cos(f')
1367 final UnivariateDerivative1 eppcfpp = epp.multiply(cfpp);
1368 // 1 + e"cos(f')
1369 final UnivariateDerivative1 eppcfppP1 = eppcfpp.add(1.);
1370 // 2 + e"cos(f')
1371 final UnivariateDerivative1 eppcfppP2 = eppcfpp.add(2.);
1372 // 3 + e"cos(f')
1373 final UnivariateDerivative1 eppcfppP3 = eppcfpp.add(3.);
1374 // (1 + e"cos(f'))³
1375 final UnivariateDerivative1 eppcfppP1_3 = eppcfppP1.square().multiply(eppcfppP1);
1376
1377 // 2g"
1378 final UnivariateDerivative1 g2 = gpp.multiply(2);
1379
1380 // 2g" + f"
1381 final UnivariateDerivative1 g2f = g2.add(fpp);
1382 final FieldSinCos<UnivariateDerivative1> sc2gf = g2f.sinCos();
1383 final UnivariateDerivative1 c2gf = sc2gf.cos();
1384 final UnivariateDerivative1 s2gf = sc2gf.sin();
1385 final UnivariateDerivative1 eppc2gf = epp.multiply(c2gf);
1386 final UnivariateDerivative1 epps2gf = epp.multiply(s2gf);
1387
1388 // 2g" + 2f"
1389 final UnivariateDerivative1 g2f2 = g2.add(fpp.multiply(2));
1390 final FieldSinCos<UnivariateDerivative1> sc2g2f = g2f2.sinCos();
1391 final UnivariateDerivative1 c2g2f = sc2g2f.cos();
1392 final UnivariateDerivative1 s2g2f = sc2g2f.sin();
1393
1394 // 2g" + 3f"
1395 final UnivariateDerivative1 g2f3 = g2.add(fpp.multiply(3));
1396 final FieldSinCos<UnivariateDerivative1> sc2g3f = g2f3.sinCos();
1397 final UnivariateDerivative1 c2g3f = sc2g3f.cos();
1398 final UnivariateDerivative1 s2g3f = sc2g3f.sin();
1399
1400 // e"cos(2g" + 3f")
1401 final UnivariateDerivative1 eppc2g3f = epp.multiply(c2g3f);
1402 // e"sin(2g" + 3f")
1403 final UnivariateDerivative1 epps2g3f = epp.multiply(s2g3f);
1404
1405 // f" + e"sin(f") - l"
1406 final UnivariateDerivative1 w17 = fpp.add(eppsfpp).subtract(lpp);
1407
1408 // ((e"cos(f") + 3)e"cos(f") + 3)cos(f")
1409 final UnivariateDerivative1 w20 = cfpp.multiply(eppcfppP3.multiply(eppcfpp).add(3.));
1410
1411 // 3sin(2g" + 2f") + 3e"sin(2g" + f") + e"sin(2g" + f")
1412 final UnivariateDerivative1 w21 = s2g2f.add(epps2gf).multiply(3).add(epps2g3f);
1413
1414 // (1 + e"cos(f"))(2 + e"cos(f"))/η²
1415 final UnivariateDerivative1 w22 = eppcfppP1.multiply(eppcfppP2).divide(n2);
1416
1417 // sinCos(I"/2)
1418 final SinCos sci = FastMath.sinCos(0.5 * mean.getI());
1419 final double siO2 = sci.sin();
1420 final double ciO2 = sci.cos();
1421
1422 // δ2a
1423 final UnivariateDerivative1 d2a = app.multiply(yp2 / n2).
1424 multiply(eppcfppP1_3.subtract(n3).multiply(ci2X3M1).
1425 add(eppcfppP1_3.multiply(c2g2f).multiply(3 * oneMci2)));
1426
1427 // δ2e
1428 final UnivariateDerivative1 d2e = (w20.add(epp.multiply(t8))).multiply(ci2X3M1).
1429 add((w20.add(epp.multiply(c2g2f))).multiply(3 * oneMci2)).
1430 subtract((eppc2gf.multiply(3).add(eppc2g3f)).multiply(n2 * oneMci2)).
1431 multiply(0.5 * yp2);
1432
1433 // δ2I
1434 final UnivariateDerivative1 d2I = ((c2g2f.add(eppc2gf)).multiply(3).add(eppc2g3f)).
1435 multiply(0.5 * yp2 * ci * si);
1436
1437 // e"δ2l
1438 final UnivariateDerivative1 eppd2l = (w22.add(1).multiply(sfpp).multiply(2 * oneMci2).
1439 add((w22.subtract(1).negate().multiply(s2gf)).
1440 add(w22.add(1. / 3.).multiply(s2g3f)).
1441 multiply(3 * oneMci2))).
1442 multiply(0.25 * yp2 * n3).negate();
1443
1444 // sinI"δ2h
1445 final UnivariateDerivative1 sIppd2h = (w21.subtract(w17.multiply(6))).
1446 multiply(0.5 * yp2 * ci * si);
1447
1448 // δ2z = δ2l + δ2g + δ2h
1449 final UnivariateDerivative1 d2z = (epp.multiply(eppd2l).multiply(t8 - 1.).divide(n3).
1450 add(w17.multiply(6. * (1. + ci * (2 - 5. * ci)))
1451 .subtract(w21.multiply(3. + ci * (2 - 5. * ci))).multiply(0.25 * yp2))).
1452 negate();
1453
1454 // Assembling elements
1455 // -------------------
1456
1457 // e" + δe
1458 final UnivariateDerivative1 de = epp.add(d1e).add(d2e);
1459
1460 // e"δl
1461 final UnivariateDerivative1 dl = eppd1l.add(eppd2l);
1462
1463 // sin(I"/2)δh = sin(I")δh / cos(I"/2) (singular for I" = π, very unlikely)
1464 final UnivariateDerivative1 dh = sIppd1h.add(sIppd2h).divide(2. * ciO2);
1465
1466 // δI
1467 final UnivariateDerivative1 di = d1I.add(d2I).multiply(0.5 * ciO2).add(siO2);
1468
1469 // z = l" + g" + h" + δ1z + δ2z
1470 final UnivariateDerivative1 z = lpp.add(gpp).add(hpp).add(d1z).add(d2z);
1471
1472 // Osculating elements
1473 // -------------------
1474
1475 // Semi-major axis
1476 final UnivariateDerivative1 a = app.add(d2a);
1477
1478 // Eccentricity
1479 final UnivariateDerivative1 e = FastMath.sqrt(de.square().add(dl.square()));
1480
1481 // Mean anomaly
1482 final FieldSinCos<UnivariateDerivative1> sclpp = lpp.sinCos();
1483 final UnivariateDerivative1 clpp = sclpp.cos();
1484 final UnivariateDerivative1 slpp = sclpp.sin();
1485 final UnivariateDerivative1 l = FastMath.atan2(de.multiply(slpp).add(dl.multiply(clpp)),
1486 de.multiply(clpp).subtract(dl.multiply(slpp)));
1487
1488 // Inclination
1489 final UnivariateDerivative1 i = FastMath.acos(di.square().add(dh.square()).multiply(2).negate().add(1.));
1490
1491 // Longitude of ascending node
1492 final FieldSinCos<UnivariateDerivative1> schpp = hpp.sinCos();
1493 final UnivariateDerivative1 chpp = schpp.cos();
1494 final UnivariateDerivative1 shpp = schpp.sin();
1495 final UnivariateDerivative1 h = FastMath.atan2(di.multiply(shpp).add(dh.multiply(chpp)),
1496 di.multiply(chpp).subtract(dh.multiply(shpp)));
1497
1498 // Argument of perigee
1499 final UnivariateDerivative1 g = z.subtract(l).subtract(h);
1500
1501 // Return a Keplerian orbit
1502 return new KeplerianOrbit(a.getValue(), e.getValue(), i.getValue(),
1503 g.getValue(), h.getValue(), l.getValue(),
1504 a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
1505 g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
1506 PositionAngleType.MEAN, mean.getFrame(), date, mu);
1507
1508 }
1509
1510 }
1511
1512 /** {@inheritDoc} */
1513 protected double getMass(final AbsoluteDate date) {
1514 return models.get(date).mass;
1515 }
1516
1517 }
1518