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