1 /* Copyright 2002-2021 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.io.Serializable;
20
21 import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.MathUtils;
26 import org.hipparchus.util.SinCos;
27 import org.orekit.attitudes.AttitudeProvider;
28 import org.orekit.attitudes.InertialProvider;
29 import org.orekit.errors.OrekitException;
30 import org.orekit.errors.OrekitMessages;
31 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
32 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
33 import org.orekit.orbits.CartesianOrbit;
34 import org.orekit.orbits.CircularOrbit;
35 import org.orekit.orbits.Orbit;
36 import org.orekit.orbits.OrbitType;
37 import org.orekit.orbits.PositionAngle;
38 import org.orekit.propagation.PropagationType;
39 import org.orekit.propagation.SpacecraftState;
40 import org.orekit.time.AbsoluteDate;
41 import org.orekit.utils.TimeSpanMap;
42 import org.orekit.utils.TimeStampedPVCoordinates;
43
44 /** This class propagates a {@link org.orekit.propagation.SpacecraftState}
45 * using the analytical Eckstein-Hechler model.
46 * <p>The Eckstein-Hechler model is suited for near circular orbits
47 * (e < 0.1, with poor accuracy between 0.005 and 0.1) and inclination
48 * neither equatorial (direct or retrograde) nor critical (direct or
49 * retrograde).</p>
50 * <p>
51 * Note that before version 7.0, there was a large inconsistency in the generated
52 * orbits, and it was fixed as of version 7.0 of Orekit, with a visible side effect.
53 * The problems is that if the circular parameters produced by the Eckstein-Hechler
54 * model are used to build an orbit considered to be osculating, the velocity deduced
55 * from this orbit was <em>inconsistent with the position evolution</em>! The reason is
56 * that the model includes non-Keplerian effects but it does not include a corresponding
57 * circular/Cartesian conversion. As a consequence, all subsequent computation involving
58 * velocity were wrong. This includes attitude modes like yaw compensation and Doppler
59 * effect. As this effect was considered serious enough and as accurate velocities were
60 * considered important, the propagator now generates {@link CartesianOrbit Cartesian
61 * orbits} which are built in a special way to ensure consistency throughout propagation.
62 * A side effect is that if circular parameters are rebuilt by user from these propagated
63 * Cartesian orbit, the circular parameters will generally <em>not</em> match the initial
64 * orbit (differences in semi-major axis can exceed 120 m). The position however <em>will</em>
65 * match to sub-micrometer level, and this position will be identical to the positions
66 * that were generated by previous versions (in other words, the internals of the models
67 * have not been changed, only the output parameters have been changed). The correctness
68 * of the initialization has been assessed and is good, as it allows the subsequent orbit
69 * to remain close to a numerical reference orbit.
70 * </p>
71 * <p>
72 * If users need a more definitive initialization of an Eckstein-Hechler propagator, they
73 * should consider using a {@link org.orekit.propagation.conversion.PropagatorConverter
74 * propagator converter} to initialize their Eckstein-Hechler propagator using a complete
75 * sample instead of just a single initial orbit.
76 * </p>
77 * @see Orbit
78 * @author Guylaine Prat
79 */
80 public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator {
81
82 /** Initial Eckstein-Hechler model. */
83 private EHModel initialModel;
84
85 /** All models. */
86 private transient TimeSpanMap<EHModel> models;
87
88 /** Reference radius of the central body attraction model (m). */
89 private double referenceRadius;
90
91 /** Central attraction coefficient (m³/s²). */
92 private double mu;
93
94 /** Un-normalized zonal coefficients. */
95 private double[] ck0;
96
97 /** Build a propagator from orbit and potential provider.
98 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
99 *
100 * <p>Using this constructor, an initial osculating orbit is considered.</p>
101 *
102 * @param initialOrbit initial orbit
103 * @param provider for un-normalized zonal coefficients
104 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider,
105 * UnnormalizedSphericalHarmonicsProvider)
106 * @see #EcksteinHechlerPropagator(Orbit, UnnormalizedSphericalHarmonicsProvider,
107 * PropagationType)
108 */
109 public EcksteinHechlerPropagator(final Orbit initialOrbit,
110 final UnnormalizedSphericalHarmonicsProvider provider) {
111 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
112 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
113 }
114
115 /**
116 * Private helper constructor.
117 * <p>Using this constructor, an initial osculating orbit is considered.</p>
118 * @param initialOrbit initial orbit
119 * @param attitude attitude provider
120 * @param mass spacecraft mass
121 * @param provider for un-normalized zonal coefficients
122 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
123 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
124 * UnnormalizedSphericalHarmonicsProvider,
125 * UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics,
126 * PropagationType)
127 */
128 public EcksteinHechlerPropagator(final Orbit initialOrbit,
129 final AttitudeProvider attitude,
130 final double mass,
131 final UnnormalizedSphericalHarmonicsProvider provider,
132 final UnnormalizedSphericalHarmonics harmonics) {
133 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
134 harmonics.getUnnormalizedCnm(2, 0),
135 harmonics.getUnnormalizedCnm(3, 0),
136 harmonics.getUnnormalizedCnm(4, 0),
137 harmonics.getUnnormalizedCnm(5, 0),
138 harmonics.getUnnormalizedCnm(6, 0));
139 }
140
141 /** Build a propagator from orbit and potential.
142 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
143 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
144 * are related to both the normalized coefficients
145 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
146 * and the J<sub>n</sub> one as follows:</p>
147 *
148 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
149 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
150 *
151 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
152 *
153 * <p>Using this constructor, an initial osculating orbit is considered.</p>
154 *
155 * @param initialOrbit initial orbit
156 * @param referenceRadius reference radius of the Earth for the potential model (m)
157 * @param mu central attraction coefficient (m³/s²)
158 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
159 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
160 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
161 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
162 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
163 * @see org.orekit.utils.Constants
164 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
165 * double, double, double, double, double)
166 */
167 public EcksteinHechlerPropagator(final Orbit initialOrbit,
168 final double referenceRadius, final double mu,
169 final double c20, final double c30, final double c40,
170 final double c50, final double c60) {
171 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
172 DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
173 }
174
175 /** Build a propagator from orbit, mass and potential provider.
176 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
177 *
178 * <p>Using this constructor, an initial osculating orbit is considered.</p>
179 *
180 * @param initialOrbit initial orbit
181 * @param mass spacecraft mass
182 * @param provider for un-normalized zonal coefficients
183 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
184 * UnnormalizedSphericalHarmonicsProvider)
185 */
186 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
187 final UnnormalizedSphericalHarmonicsProvider provider) {
188 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
189 mass, provider, provider.onDate(initialOrbit.getDate()));
190 }
191
192 /** Build a propagator from orbit, mass and potential.
193 * <p>Attitude law is set to an unspecified non-null arbitrary value.</p>
194 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
195 * are related to both the normalized coefficients
196 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
197 * and the J<sub>n</sub> one as follows:</p>
198 *
199 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
200 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
201 *
202 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
203 *
204 * <p>Using this constructor, an initial osculating orbit is considered.</p>
205 *
206 * @param initialOrbit initial orbit
207 * @param mass spacecraft mass
208 * @param referenceRadius reference radius of the Earth for the potential model (m)
209 * @param mu central attraction coefficient (m³/s²)
210 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
211 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
212 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
213 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
214 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
215 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
216 * double, double, double, double, double)
217 */
218 public EcksteinHechlerPropagator(final Orbit initialOrbit, final double mass,
219 final double referenceRadius, final double mu,
220 final double c20, final double c30, final double c40,
221 final double c50, final double c60) {
222 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
223 mass, referenceRadius, mu, c20, c30, c40, c50, c60);
224 }
225
226 /** Build a propagator from orbit, attitude provider and potential provider.
227 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
228 * <p>Using this constructor, an initial osculating orbit is considered.</p>
229 * @param initialOrbit initial orbit
230 * @param attitudeProv attitude provider
231 * @param provider for un-normalized zonal coefficients
232 */
233 public EcksteinHechlerPropagator(final Orbit initialOrbit,
234 final AttitudeProvider attitudeProv,
235 final UnnormalizedSphericalHarmonicsProvider provider) {
236 this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()));
237 }
238
239 /** Build a propagator from orbit, attitude provider and potential.
240 * <p>Mass is set to an unspecified non-null arbitrary value.</p>
241 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
242 * are related to both the normalized coefficients
243 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
244 * and the J<sub>n</sub> one as follows:</p>
245 *
246 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
247 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
248 *
249 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
250 *
251 * <p>Using this constructor, an initial osculating orbit is considered.</p>
252 *
253 * @param initialOrbit initial orbit
254 * @param attitudeProv attitude provider
255 * @param referenceRadius reference radius of the Earth for the potential model (m)
256 * @param mu central attraction coefficient (m³/s²)
257 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
258 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
259 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
260 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
261 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
262 */
263 public EcksteinHechlerPropagator(final Orbit initialOrbit,
264 final AttitudeProvider attitudeProv,
265 final double referenceRadius, final double mu,
266 final double c20, final double c30, final double c40,
267 final double c50, final double c60) {
268 this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, c60);
269 }
270
271 /** Build a propagator from orbit, attitude provider, mass and potential provider.
272 * <p>Using this constructor, an initial osculating orbit is considered.</p>
273 * @param initialOrbit initial orbit
274 * @param attitudeProv attitude provider
275 * @param mass spacecraft mass
276 * @param provider for un-normalized zonal coefficients
277 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double,
278 * UnnormalizedSphericalHarmonicsProvider, PropagationType)
279 */
280 public EcksteinHechlerPropagator(final Orbit initialOrbit,
281 final AttitudeProvider attitudeProv,
282 final double mass,
283 final UnnormalizedSphericalHarmonicsProvider provider) {
284 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()));
285 }
286
287 /** Build a propagator from orbit, attitude provider, mass and potential.
288 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
289 * are related to both the normalized coefficients
290 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
291 * and the J<sub>n</sub> one as follows:</p>
292 *
293 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
294 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
295 *
296 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
297 *
298 * <p>Using this constructor, an initial osculating orbit is considered.</p>
299 *
300 * @param initialOrbit initial orbit
301 * @param attitudeProv attitude provider
302 * @param mass spacecraft mass
303 * @param referenceRadius reference radius of the Earth for the potential model (m)
304 * @param mu central attraction coefficient (m³/s²)
305 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
306 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
307 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
308 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
309 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
310 * @see #EcksteinHechlerPropagator(Orbit, AttitudeProvider, double, double, double,
311 * double, double, double, double, double, PropagationType)
312 */
313 public EcksteinHechlerPropagator(final Orbit initialOrbit,
314 final AttitudeProvider attitudeProv,
315 final double mass,
316 final double referenceRadius, final double mu,
317 final double c20, final double c30, final double c40,
318 final double c50, final double c60) {
319 this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, c60,
320 PropagationType.OSCULATING);
321 }
322
323
324 /** Build a propagator from orbit and potential provider.
325 * <p>Mass and attitude provider are set to unspecified non-null arbitrary values.</p>
326 *
327 * <p>Using this constructor, it is possible to define the initial orbit as
328 * a mean Eckstein-Hechler orbit or an osculating one.</p>
329 *
330 * @param initialOrbit initial orbit
331 * @param provider for un-normalized zonal coefficients
332 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
333 * @since 10.2
334 */
335 public EcksteinHechlerPropagator(final Orbit initialOrbit,
336 final UnnormalizedSphericalHarmonicsProvider provider,
337 final PropagationType initialType) {
338 this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
339 DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType);
340 }
341
342 /** Build a propagator from orbit, attitude provider, mass and potential provider.
343 * <p>Using this constructor, it is possible to define the initial orbit as
344 * a mean Eckstein-Hechler orbit or an osculating one.</p>
345 * @param initialOrbit initial orbit
346 * @param attitudeProv attitude provider
347 * @param mass spacecraft mass
348 * @param provider for un-normalized zonal coefficients
349 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
350 * @since 10.2
351 */
352 public EcksteinHechlerPropagator(final Orbit initialOrbit,
353 final AttitudeProvider attitudeProv,
354 final double mass,
355 final UnnormalizedSphericalHarmonicsProvider provider,
356 final PropagationType initialType) {
357 this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType);
358 }
359
360 /**
361 * Private helper constructor.
362 * <p>Using this constructor, it is possible to define the initial orbit as
363 * a mean Eckstein-Hechler orbit or an osculating one.</p>
364 * @param initialOrbit initial orbit
365 * @param attitude attitude provider
366 * @param mass spacecraft mass
367 * @param provider for un-normalized zonal coefficients
368 * @param harmonics {@code provider.onDate(initialOrbit.getDate())}
369 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
370 * @since 10.2
371 */
372 public EcksteinHechlerPropagator(final Orbit initialOrbit,
373 final AttitudeProvider attitude,
374 final double mass,
375 final UnnormalizedSphericalHarmonicsProvider provider,
376 final UnnormalizedSphericalHarmonics harmonics,
377 final PropagationType initialType) {
378 this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(),
379 harmonics.getUnnormalizedCnm(2, 0),
380 harmonics.getUnnormalizedCnm(3, 0),
381 harmonics.getUnnormalizedCnm(4, 0),
382 harmonics.getUnnormalizedCnm(5, 0),
383 harmonics.getUnnormalizedCnm(6, 0),
384 initialType);
385 }
386
387 /** Build a propagator from orbit, attitude provider, mass and potential.
388 * <p>The C<sub>n,0</sub> coefficients are the denormalized zonal coefficients, they
389 * are related to both the normalized coefficients
390 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
391 * and the J<sub>n</sub> one as follows:</p>
392 *
393 * <p> C<sub>n,0</sub> = [(2-δ<sub>0,m</sub>)(2n+1)(n-m)!/(n+m)!]<sup>½</sup>
394 * <span style="text-decoration: overline">C</span><sub>n,0</sub>
395 *
396 * <p> C<sub>n,0</sub> = -J<sub>n</sub>
397 *
398 * <p>Using this constructor, it is possible to define the initial orbit as
399 * a mean Eckstein-Hechler orbit or an osculating one.</p>
400 *
401 * @param initialOrbit initial orbit
402 * @param attitudeProv attitude provider
403 * @param mass spacecraft mass
404 * @param referenceRadius reference radius of the Earth for the potential model (m)
405 * @param mu central attraction coefficient (m³/s²)
406 * @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
407 * @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
408 * @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
409 * @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
410 * @param c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
411 * @param initialType initial orbit type (mean Eckstein-Hechler orbit or osculating orbit)
412 * @since 10.2
413 */
414 public EcksteinHechlerPropagator(final Orbit initialOrbit,
415 final AttitudeProvider attitudeProv,
416 final double mass,
417 final double referenceRadius, final double mu,
418 final double c20, final double c30, final double c40,
419 final double c50, final double c60,
420 final PropagationType initialType) {
421
422 super(attitudeProv);
423
424 // store model coefficients
425 this.referenceRadius = referenceRadius;
426 this.mu = mu;
427 this.ck0 = new double[] {
428 0.0, 0.0, c20, c30, c40, c50, c60
429 };
430
431 // compute mean parameters if needed
432 // transform into circular adapted parameters used by the Eckstein-Hechler model
433 resetInitialState(new SpacecraftState(initialOrbit,
434 attitudeProv.getAttitude(initialOrbit,
435 initialOrbit.getDate(),
436 initialOrbit.getFrame()),
437 mass),
438 initialType);
439
440 }
441
442 /** {@inheritDoc}
443 * <p>The new initial state to consider
444 * must be defined with an osculating orbit.</p>
445 * @see #resetInitialState(SpacecraftState, PropagationType)
446 */
447 public void resetInitialState(final SpacecraftState state) {
448 resetInitialState(state, PropagationType.OSCULATING);
449 }
450
451 /** Reset the propagator initial state.
452 * @param state new initial state to consider
453 * @param stateType mean Eckstein-Hechler orbit or osculating orbit
454 * @since 10.2
455 */
456 public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
457 super.resetInitialState(state);
458 final CircularOrbit circular = (CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit());
459 this.initialModel = (stateType == PropagationType.MEAN) ?
460 new EHModel(circular, state.getMass(), referenceRadius, mu, ck0) :
461 computeMeanParameters(circular, state.getMass());
462 this.models = new TimeSpanMap<EHModel>(initialModel);
463 }
464
465 /** {@inheritDoc} */
466 protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
467 final EHModel newModel = computeMeanParameters((CircularOrbit) OrbitType.CIRCULAR.convertType(state.getOrbit()),
468 state.getMass());
469 if (forward) {
470 models.addValidAfter(newModel, state.getDate());
471 } else {
472 models.addValidBefore(newModel, state.getDate());
473 }
474 stateChanged(state);
475 }
476
477 /** Compute mean parameters according to the Eckstein-Hechler analytical model.
478 * @param osculating osculating orbit
479 * @param mass constant mass
480 * @return Eckstein-Hechler mean model
481 */
482 private EHModel computeMeanParameters(final CircularOrbit osculating, final double mass) {
483
484 // sanity check
485 if (osculating.getA() < referenceRadius) {
486 throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
487 osculating.getA());
488 }
489
490 // rough initialization of the mean parameters
491 EHModel current = new EHModel(osculating, mass, referenceRadius, mu, ck0);
492
493 // threshold for each parameter
494 final double epsilon = 1.0e-13;
495 final double thresholdA = epsilon * (1 + FastMath.abs(current.mean.getA()));
496 final double thresholdE = epsilon * (1 + current.mean.getE());
497 final double thresholdAngles = epsilon * FastMath.PI;
498
499 int i = 0;
500 while (i++ < 100) {
501
502 // recompute the osculating parameters from the current mean parameters
503 final UnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate());
504
505 // adapted parameters residuals
506 final double deltaA = osculating.getA() - parameters[0].getValue();
507 final double deltaEx = osculating.getCircularEx() - parameters[1].getValue();
508 final double deltaEy = osculating.getCircularEy() - parameters[2].getValue();
509 final double deltaI = osculating.getI() - parameters[3].getValue();
510 final double deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
511 parameters[4].getValue(),
512 0.0);
513 final double deltaAlphaM = MathUtils.normalizeAngle(osculating.getAlphaM() - parameters[5].getValue(), 0.0);
514
515 // update mean parameters
516 current = new EHModel(new CircularOrbit(current.mean.getA() + deltaA,
517 current.mean.getCircularEx() + deltaEx,
518 current.mean.getCircularEy() + deltaEy,
519 current.mean.getI() + deltaI,
520 current.mean.getRightAscensionOfAscendingNode() + deltaRAAN,
521 current.mean.getAlphaM() + deltaAlphaM,
522 PositionAngle.MEAN,
523 current.mean.getFrame(),
524 current.mean.getDate(), mu),
525 mass, referenceRadius, mu, ck0);
526
527 // check convergence
528 if (FastMath.abs(deltaA) < thresholdA &&
529 FastMath.abs(deltaEx) < thresholdE &&
530 FastMath.abs(deltaEy) < thresholdE &&
531 FastMath.abs(deltaI) < thresholdAngles &&
532 FastMath.abs(deltaRAAN) < thresholdAngles &&
533 FastMath.abs(deltaAlphaM) < thresholdAngles) {
534 return current;
535 }
536
537 }
538
539 throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_ECKSTEIN_HECHLER_MEAN_PARAMETERS, i);
540
541 }
542
543 /** {@inheritDoc} */
544 public CartesianOrbit propagateOrbit(final AbsoluteDate date) {
545 // compute Cartesian parameters, taking derivatives into account
546 // to make sure velocity and acceleration are consistent
547 final EHModel current = models.get(date);
548 return new CartesianOrbit(toCartesian(date, current.propagateParameters(date)),
549 current.mean.getFrame(), mu);
550 }
551
552 /** Local class for Eckstein-Hechler model, with fixed mean parameters. */
553 private static class EHModel implements Serializable {
554
555 /** Serializable UID. */
556 private static final long serialVersionUID = 20160115L;
557
558 /** Mean orbit. */
559 private final CircularOrbit mean;
560
561 /** Constant mass. */
562 private final double mass;
563
564 // CHECKSTYLE: stop JavadocVariable check
565
566 // preprocessed values
567 private final double xnotDot;
568 private final double rdpom;
569 private final double rdpomp;
570 private final double eps1;
571 private final double eps2;
572 private final double xim;
573 private final double ommD;
574 private final double rdl;
575 private final double aMD;
576
577 private final double kh;
578 private final double kl;
579
580 private final double ax1;
581 private final double ay1;
582 private final double as1;
583 private final double ac2;
584 private final double axy3;
585 private final double as3;
586 private final double ac4;
587 private final double as5;
588 private final double ac6;
589
590 private final double ex1;
591 private final double exx2;
592 private final double exy2;
593 private final double ex3;
594 private final double ex4;
595
596 private final double ey1;
597 private final double eyx2;
598 private final double eyy2;
599 private final double ey3;
600 private final double ey4;
601
602 private final double rx1;
603 private final double ry1;
604 private final double r2;
605 private final double r3;
606 private final double rl;
607
608 private final double iy1;
609 private final double ix1;
610 private final double i2;
611 private final double i3;
612 private final double ih;
613
614 private final double lx1;
615 private final double ly1;
616 private final double l2;
617 private final double l3;
618 private final double ll;
619
620 // CHECKSTYLE: resume JavadocVariable check
621
622 /** Create a model for specified mean orbit.
623 * @param mean mean orbit
624 * @param mass constant mass
625 * @param referenceRadius reference radius of the central body attraction model (m)
626 * @param mu central attraction coefficient (m³/s²)
627 * @param ck0 un-normalized zonal coefficients
628 */
629 EHModel(final CircularOrbit mean, final double mass,
630 final double referenceRadius, final double mu, final double[] ck0) {
631
632 this.mean = mean;
633 this.mass = mass;
634
635 // preliminary processing
636 double q = referenceRadius / mean.getA();
637 double ql = q * q;
638 final double g2 = ck0[2] * ql;
639 ql *= q;
640 final double g3 = ck0[3] * ql;
641 ql *= q;
642 final double g4 = ck0[4] * ql;
643 ql *= q;
644 final double g5 = ck0[5] * ql;
645 ql *= q;
646 final double g6 = ck0[6] * ql;
647
648 final SinCos sc = FastMath.sinCos(mean.getI());
649 final double cosI1 = sc.cos();
650 final double sinI1 = sc.sin();
651 final double sinI2 = sinI1 * sinI1;
652 final double sinI4 = sinI2 * sinI2;
653 final double sinI6 = sinI2 * sinI4;
654
655 if (sinI2 < 1.0e-10) {
656 throw new OrekitException(OrekitMessages.ALMOST_EQUATORIAL_ORBIT,
657 FastMath.toDegrees(mean.getI()));
658 }
659
660 if (FastMath.abs(sinI2 - 4.0 / 5.0) < 1.0e-3) {
661 throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT,
662 FastMath.toDegrees(mean.getI()));
663 }
664
665 if (mean.getE() > 0.1) {
666 // if 0.005 < e < 0.1 no error is triggered, but accuracy is poor
667 throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL,
668 mean.getE());
669 }
670
671 xnotDot = FastMath.sqrt(mu / mean.getA()) / mean.getA();
672
673 rdpom = -0.75 * g2 * (4.0 - 5.0 * sinI2);
674 rdpomp = 7.5 * g4 * (1.0 - 31.0 / 8.0 * sinI2 + 49.0 / 16.0 * sinI4) -
675 13.125 * g6 * (1.0 - 8.0 * sinI2 + 129.0 / 8.0 * sinI4 - 297.0 / 32.0 * sinI6);
676
677 q = 3.0 / (32.0 * rdpom);
678 eps1 = q * g4 * sinI2 * (30.0 - 35.0 * sinI2) -
679 175.0 * q * g6 * sinI2 * (1.0 - 3.0 * sinI2 + 2.0625 * sinI4);
680 q = 3.0 * sinI1 / (8.0 * rdpom);
681 eps2 = q * g3 * (4.0 - 5.0 * sinI2) - q * g5 * (10.0 - 35.0 * sinI2 + 26.25 * sinI4);
682
683 xim = mean.getI();
684 ommD = cosI1 * (1.50 * g2 - 2.25 * g2 * g2 * (2.5 - 19.0 / 6.0 * sinI2) +
685 0.9375 * g4 * (7.0 * sinI2 - 4.0) +
686 3.28125 * g6 * (2.0 - 9.0 * sinI2 + 8.25 * sinI4));
687
688 rdl = 1.0 - 1.50 * g2 * (3.0 - 4.0 * sinI2);
689 aMD = rdl +
690 2.25 * g2 * g2 * (9.0 - 263.0 / 12.0 * sinI2 + 341.0 / 24.0 * sinI4) +
691 15.0 / 16.0 * g4 * (8.0 - 31.0 * sinI2 + 24.5 * sinI4) +
692 105.0 / 32.0 * g6 * (-10.0 / 3.0 + 25.0 * sinI2 - 48.75 * sinI4 + 27.5 * sinI6);
693
694 final double qq = -1.5 * g2 / rdl;
695 final double qA = 0.75 * g2 * g2 * sinI2;
696 final double qB = 0.25 * g4 * sinI2;
697 final double qC = 105.0 / 16.0 * g6 * sinI2;
698 final double qD = -0.75 * g3 * sinI1;
699 final double qE = 3.75 * g5 * sinI1;
700 kh = 0.375 / rdpom;
701 kl = kh / sinI1;
702
703 ax1 = qq * (2.0 - 3.5 * sinI2);
704 ay1 = qq * (2.0 - 2.5 * sinI2);
705 as1 = qD * (4.0 - 5.0 * sinI2) +
706 qE * (2.625 * sinI4 - 3.5 * sinI2 + 1.0);
707 ac2 = qq * sinI2 +
708 qA * 7.0 * (2.0 - 3.0 * sinI2) +
709 qB * (15.0 - 17.5 * sinI2) +
710 qC * (3.0 * sinI2 - 1.0 - 33.0 / 16.0 * sinI4);
711 axy3 = qq * 3.5 * sinI2;
712 as3 = qD * 5.0 / 3.0 * sinI2 +
713 qE * 7.0 / 6.0 * sinI2 * (1.0 - 1.125 * sinI2);
714 ac4 = qA * sinI2 +
715 qB * 4.375 * sinI2 +
716 qC * 0.75 * (1.1 * sinI4 - sinI2);
717
718 as5 = qE * 21.0 / 80.0 * sinI4;
719
720 ac6 = qC * -11.0 / 80.0 * sinI4;
721
722 ex1 = qq * (1.0 - 1.25 * sinI2);
723 exx2 = qq * 0.5 * (3.0 - 5.0 * sinI2);
724 exy2 = qq * (2.0 - 1.5 * sinI2);
725 ex3 = qq * 7.0 / 12.0 * sinI2;
726 ex4 = qq * 17.0 / 8.0 * sinI2;
727
728 ey1 = qq * (1.0 - 1.75 * sinI2);
729 eyx2 = qq * (1.0 - 3.0 * sinI2);
730 eyy2 = qq * (2.0 * sinI2 - 1.5);
731 ey3 = qq * 7.0 / 12.0 * sinI2;
732 ey4 = qq * 17.0 / 8.0 * sinI2;
733
734 q = -qq * cosI1;
735 rx1 = 3.5 * q;
736 ry1 = -2.5 * q;
737 r2 = -0.5 * q;
738 r3 = 7.0 / 6.0 * q;
739 rl = g3 * cosI1 * (4.0 - 15.0 * sinI2) -
740 2.5 * g5 * cosI1 * (4.0 - 42.0 * sinI2 + 52.5 * sinI4);
741
742 q = 0.5 * qq * sinI1 * cosI1;
743 iy1 = q;
744 ix1 = -q;
745 i2 = q;
746 i3 = q * 7.0 / 3.0;
747 ih = -g3 * cosI1 * (4.0 - 5.0 * sinI2) +
748 2.5 * g5 * cosI1 * (4.0 - 14.0 * sinI2 + 10.5 * sinI4);
749
750 lx1 = qq * (7.0 - 77.0 / 8.0 * sinI2);
751 ly1 = qq * (55.0 / 8.0 * sinI2 - 7.50);
752 l2 = qq * (1.25 * sinI2 - 0.5);
753 l3 = qq * (77.0 / 24.0 * sinI2 - 7.0 / 6.0);
754 ll = g3 * (53.0 * sinI2 - 4.0 - 57.5 * sinI4) +
755 2.5 * g5 * (4.0 - 96.0 * sinI2 + 269.5 * sinI4 - 183.75 * sinI6);
756
757 }
758
759 /** Extrapolate an orbit up to a specific target date.
760 * @param date target date for the orbit
761 * @return propagated parameters
762 */
763 public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) {
764
765 // Keplerian evolution
766 final UnivariateDerivative2 dt = new UnivariateDerivative2(date.durationFrom(mean.getDate()), 1.0, 0.0);
767 final UnivariateDerivative2 xnot = dt.multiply(xnotDot);
768
769 // secular effects
770
771 // eccentricity
772 final UnivariateDerivative2 x = xnot.multiply(rdpom + rdpomp);
773 final UnivariateDerivative2 cx = x.cos();
774 final UnivariateDerivative2 sx = x.sin();
775 final UnivariateDerivative2 exm = cx.multiply(mean.getCircularEx()).
776 add(sx.multiply(eps2 - (1.0 - eps1) * mean.getCircularEy()));
777 final UnivariateDerivative2 eym = sx.multiply((1.0 + eps1) * mean.getCircularEx()).
778 add(cx.multiply(mean.getCircularEy() - eps2)).
779 add(eps2);
780
781 // no secular effect on inclination
782
783 // right ascension of ascending node
784 final UnivariateDerivative2 omm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getRightAscensionOfAscendingNode() + ommD * xnot.getValue(),
785 FastMath.PI),
786 ommD * xnotDot,
787 0.0);
788
789 // latitude argument
790 final UnivariateDerivative2 xlm = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getAlphaM() + aMD * xnot.getValue(),
791 FastMath.PI),
792 aMD * xnotDot,
793 0.0);
794
795 // periodical terms
796 final UnivariateDerivative2 cl1 = xlm.cos();
797 final UnivariateDerivative2 sl1 = xlm.sin();
798 final UnivariateDerivative2 cl2 = cl1.multiply(cl1).subtract(sl1.multiply(sl1));
799 final UnivariateDerivative2 sl2 = cl1.multiply(sl1).add(sl1.multiply(cl1));
800 final UnivariateDerivative2 cl3 = cl2.multiply(cl1).subtract(sl2.multiply(sl1));
801 final UnivariateDerivative2 sl3 = cl2.multiply(sl1).add(sl2.multiply(cl1));
802 final UnivariateDerivative2 cl4 = cl3.multiply(cl1).subtract(sl3.multiply(sl1));
803 final UnivariateDerivative2 sl4 = cl3.multiply(sl1).add(sl3.multiply(cl1));
804 final UnivariateDerivative2 cl5 = cl4.multiply(cl1).subtract(sl4.multiply(sl1));
805 final UnivariateDerivative2 sl5 = cl4.multiply(sl1).add(sl4.multiply(cl1));
806 final UnivariateDerivative2 cl6 = cl5.multiply(cl1).subtract(sl5.multiply(sl1));
807
808 final UnivariateDerivative2 qh = eym.subtract(eps2).multiply(kh);
809 final UnivariateDerivative2 ql = exm.multiply(kl);
810
811 final UnivariateDerivative2 exmCl1 = exm.multiply(cl1);
812 final UnivariateDerivative2 exmSl1 = exm.multiply(sl1);
813 final UnivariateDerivative2 eymCl1 = eym.multiply(cl1);
814 final UnivariateDerivative2 eymSl1 = eym.multiply(sl1);
815 final UnivariateDerivative2 exmCl2 = exm.multiply(cl2);
816 final UnivariateDerivative2 exmSl2 = exm.multiply(sl2);
817 final UnivariateDerivative2 eymCl2 = eym.multiply(cl2);
818 final UnivariateDerivative2 eymSl2 = eym.multiply(sl2);
819 final UnivariateDerivative2 exmCl3 = exm.multiply(cl3);
820 final UnivariateDerivative2 exmSl3 = exm.multiply(sl3);
821 final UnivariateDerivative2 eymCl3 = eym.multiply(cl3);
822 final UnivariateDerivative2 eymSl3 = eym.multiply(sl3);
823 final UnivariateDerivative2 exmCl4 = exm.multiply(cl4);
824 final UnivariateDerivative2 exmSl4 = exm.multiply(sl4);
825 final UnivariateDerivative2 eymCl4 = eym.multiply(cl4);
826 final UnivariateDerivative2 eymSl4 = eym.multiply(sl4);
827
828 // semi major axis
829 final UnivariateDerivative2 rda = exmCl1.multiply(ax1).
830 add(eymSl1.multiply(ay1)).
831 add(sl1.multiply(as1)).
832 add(cl2.multiply(ac2)).
833 add(exmCl3.add(eymSl3).multiply(axy3)).
834 add(sl3.multiply(as3)).
835 add(cl4.multiply(ac4)).
836 add(sl5.multiply(as5)).
837 add(cl6.multiply(ac6));
838
839 // eccentricity
840 final UnivariateDerivative2 rdex = cl1.multiply(ex1).
841 add(exmCl2.multiply(exx2)).
842 add(eymSl2.multiply(exy2)).
843 add(cl3.multiply(ex3)).
844 add(exmCl4.add(eymSl4).multiply(ex4));
845 final UnivariateDerivative2 rdey = sl1.multiply(ey1).
846 add(exmSl2.multiply(eyx2)).
847 add(eymCl2.multiply(eyy2)).
848 add(sl3.multiply(ey3)).
849 add(exmSl4.subtract(eymCl4).multiply(ey4));
850
851 // ascending node
852 final UnivariateDerivative2 rdom = exmSl1.multiply(rx1).
853 add(eymCl1.multiply(ry1)).
854 add(sl2.multiply(r2)).
855 add(eymCl3.subtract(exmSl3).multiply(r3)).
856 add(ql.multiply(rl));
857
858 // inclination
859 final UnivariateDerivative2 rdxi = eymSl1.multiply(iy1).
860 add(exmCl1.multiply(ix1)).
861 add(cl2.multiply(i2)).
862 add(exmCl3.add(eymSl3).multiply(i3)).
863 add(qh.multiply(ih));
864
865 // latitude argument
866 final UnivariateDerivative2 rdxl = exmSl1.multiply(lx1).
867 add(eymCl1.multiply(ly1)).
868 add(sl2.multiply(l2)).
869 add(exmSl3.subtract(eymCl3).multiply(l3)).
870 add(ql.multiply(ll));
871
872 // osculating parameters
873 return new UnivariateDerivative2[] {
874 rda.add(1.0).multiply(mean.getA()),
875 rdex.add(exm),
876 rdey.add(eym),
877 rdxi.add(xim),
878 rdom.add(omm),
879 rdxl.add(xlm)
880 };
881
882 }
883
884 }
885
886 /** Convert circular parameters <em>with derivatives</em> to Cartesian coordinates.
887 * @param date date of the orbital parameters
888 * @param parameters circular parameters (a, ex, ey, i, raan, alphaM)
889 * @return Cartesian coordinates consistent with values and derivatives
890 */
891 private TimeStampedPVCoordinates toCartesian(final AbsoluteDate date, final UnivariateDerivative2[] parameters) {
892
893 // evaluate coordinates in the orbit canonical reference frame
894 final UnivariateDerivative2 cosOmega = parameters[4].cos();
895 final UnivariateDerivative2 sinOmega = parameters[4].sin();
896 final UnivariateDerivative2 cosI = parameters[3].cos();
897 final UnivariateDerivative2 sinI = parameters[3].sin();
898 final UnivariateDerivative2 alphaE = meanToEccentric(parameters[5], parameters[1], parameters[2]);
899 final UnivariateDerivative2 cosAE = alphaE.cos();
900 final UnivariateDerivative2 sinAE = alphaE.sin();
901 final UnivariateDerivative2 ex2 = parameters[1].multiply(parameters[1]);
902 final UnivariateDerivative2 ey2 = parameters[2].multiply(parameters[2]);
903 final UnivariateDerivative2 exy = parameters[1].multiply(parameters[2]);
904 final UnivariateDerivative2 q = ex2.add(ey2).subtract(1).negate().sqrt();
905 final UnivariateDerivative2 beta = q.add(1).reciprocal();
906 final UnivariateDerivative2 bx2 = beta.multiply(ex2);
907 final UnivariateDerivative2 by2 = beta.multiply(ey2);
908 final UnivariateDerivative2 bxy = beta.multiply(exy);
909 final UnivariateDerivative2 u = bxy.multiply(sinAE).subtract(parameters[1].add(by2.subtract(1).multiply(cosAE)));
910 final UnivariateDerivative2 v = bxy.multiply(cosAE).subtract(parameters[2].add(bx2.subtract(1).multiply(sinAE)));
911 final UnivariateDerivative2 x = parameters[0].multiply(u);
912 final UnivariateDerivative2 y = parameters[0].multiply(v);
913
914 // canonical orbit reference frame
915 final FieldVector3D<UnivariateDerivative2> p =
916 new FieldVector3D<>(x.multiply(cosOmega).subtract(y.multiply(cosI.multiply(sinOmega))),
917 x.multiply(sinOmega).add(y.multiply(cosI.multiply(cosOmega))),
918 y.multiply(sinI));
919
920 // dispatch derivatives
921 final Vector3D p0 = new Vector3D(p.getX().getValue(),
922 p.getY().getValue(),
923 p.getZ().getValue());
924 final Vector3D p1 = new Vector3D(p.getX().getFirstDerivative(),
925 p.getY().getFirstDerivative(),
926 p.getZ().getFirstDerivative());
927 final Vector3D p2 = new Vector3D(p.getX().getSecondDerivative(),
928 p.getY().getSecondDerivative(),
929 p.getZ().getSecondDerivative());
930 return new TimeStampedPVCoordinates(date, p0, p1, p2);
931
932 }
933
934 /** Computes the eccentric latitude argument from the mean latitude argument.
935 * @param alphaM = M + Ω mean latitude argument (rad)
936 * @param ex e cos(Ω), first component of circular eccentricity vector
937 * @param ey e sin(Ω), second component of circular eccentricity vector
938 * @return the eccentric latitude argument.
939 */
940 private UnivariateDerivative2 meanToEccentric(final UnivariateDerivative2 alphaM,
941 final UnivariateDerivative2 ex,
942 final UnivariateDerivative2 ey) {
943 // Generalization of Kepler equation to circular parameters
944 // with alphaE = PA + E and
945 // alphaM = PA + M = alphaE - ex.sin(alphaE) + ey.cos(alphaE)
946 UnivariateDerivative2 alphaE = alphaM;
947 UnivariateDerivative2 shift = alphaM.getField().getZero();
948 UnivariateDerivative2 alphaEMalphaM = alphaM.getField().getZero();
949 UnivariateDerivative2 cosAlphaE = alphaE.cos();
950 UnivariateDerivative2 sinAlphaE = alphaE.sin();
951 int iter = 0;
952 do {
953 final UnivariateDerivative2 f2 = ex.multiply(sinAlphaE).subtract(ey.multiply(cosAlphaE));
954 final UnivariateDerivative2 f1 = alphaM.getField().getOne().subtract(ex.multiply(cosAlphaE)).subtract(ey.multiply(sinAlphaE));
955 final UnivariateDerivative2 f0 = alphaEMalphaM.subtract(f2);
956
957 final UnivariateDerivative2 f12 = f1.multiply(2);
958 shift = f0.multiply(f12).divide(f1.multiply(f12).subtract(f0.multiply(f2)));
959
960 alphaEMalphaM = alphaEMalphaM.subtract(shift);
961 alphaE = alphaM.add(alphaEMalphaM);
962 cosAlphaE = alphaE.cos();
963 sinAlphaE = alphaE.sin();
964
965 } while (++iter < 50 && FastMath.abs(shift.getValue()) > 1.0e-12);
966
967 return alphaE;
968
969 }
970
971 /** {@inheritDoc} */
972 protected double getMass(final AbsoluteDate date) {
973 return models.get(date).mass;
974 }
975
976 }