1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.analytical;
18
19 import java.io.IOException;
20 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.Collection;
23 import java.util.List;
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.exception.DummyLocalizable;
27 import org.hipparchus.geometry.euclidean.threed.Rotation;
28 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
31 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
32 import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
33 import org.hipparchus.stat.descriptive.rank.Max;
34 import org.hipparchus.stat.descriptive.rank.Min;
35 import org.hipparchus.util.FastMath;
36 import org.hipparchus.util.MathUtils;
37 import org.junit.jupiter.api.AfterEach;
38 import org.junit.jupiter.api.Assertions;
39 import org.junit.jupiter.api.BeforeEach;
40 import org.junit.jupiter.api.Test;
41 import org.orekit.Utils;
42 import org.orekit.attitudes.Attitude;
43 import org.orekit.attitudes.AttitudeProvider;
44 import org.orekit.attitudes.FieldAttitude;
45 import org.orekit.attitudes.LofOffset;
46 import org.orekit.bodies.GeodeticPoint;
47 import org.orekit.bodies.OneAxisEllipsoid;
48 import org.orekit.errors.OrekitException;
49 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
50 import org.orekit.forces.gravity.potential.GravityFieldFactory;
51 import org.orekit.forces.gravity.potential.TideSystem;
52 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
53 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
54 import org.orekit.frames.Frame;
55 import org.orekit.frames.FramesFactory;
56 import org.orekit.frames.LOFType;
57 import org.orekit.frames.TopocentricFrame;
58 import org.orekit.frames.Transform;
59 import org.orekit.orbits.CircularOrbit;
60 import org.orekit.orbits.EquinoctialOrbit;
61 import org.orekit.orbits.KeplerianOrbit;
62 import org.orekit.orbits.Orbit;
63 import org.orekit.orbits.OrbitType;
64 import org.orekit.orbits.PositionAngleType;
65 import org.orekit.propagation.PropagationType;
66 import org.orekit.propagation.Propagator;
67 import org.orekit.propagation.SpacecraftState;
68 import org.orekit.propagation.ToleranceProvider;
69 import org.orekit.propagation.conversion.EcksteinHechlerPropagatorBuilder;
70 import org.orekit.propagation.conversion.FiniteDifferencePropagatorConverter;
71 import org.orekit.propagation.conversion.PropagatorConverter;
72 import org.orekit.propagation.events.ApsideDetector;
73 import org.orekit.propagation.events.DateDetector;
74 import org.orekit.propagation.events.ElevationDetector;
75 import org.orekit.propagation.events.EventDetector;
76 import org.orekit.propagation.events.NodeDetector;
77 import org.orekit.propagation.events.handlers.ContinueOnEvent;
78 import org.orekit.propagation.numerical.NumericalPropagator;
79 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
80 import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
81 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
82 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
83 import org.orekit.time.AbsoluteDate;
84 import org.orekit.time.DateComponents;
85 import org.orekit.time.FieldAbsoluteDate;
86 import org.orekit.time.TimeComponents;
87 import org.orekit.time.TimeInterpolator;
88 import org.orekit.time.TimeScalesFactory;
89 import org.orekit.utils.CartesianDerivativesFilter;
90 import org.orekit.utils.Constants;
91 import org.orekit.utils.FieldPVCoordinatesProvider;
92 import org.orekit.utils.IERSConventions;
93 import org.orekit.utils.PVCoordinates;
94 import org.orekit.utils.PVCoordinatesProvider;
95 import org.orekit.utils.TimeStampedPVCoordinates;
96 import org.orekit.utils.TimeStampedPVCoordinatesHermiteInterpolator;
97
98
99 public class EcksteinHechlerPropagatorTest {
100
101 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
102
103 @Test
104 public void sameDateCartesian() {
105
106
107
108
109 Vector3D position = new Vector3D(3220103., 69623., 6449822.);
110 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
111
112 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
113 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
114 FramesFactory.getEME2000(), initDate, provider.getMu());
115
116
117
118 EcksteinHechlerPropagator extrapolator =
119 new EcksteinHechlerPropagator(initialOrbit, provider);
120
121
122
123 SpacecraftState finalOrbit = extrapolator.propagate(initDate);
124
125
126 Assertions.assertEquals(0.0,
127 Vector3D.distance(initialOrbit.getPosition(),
128 finalOrbit.getPosition()),
129 1.3e-6);
130
131
132
133
134
135
136
137
138
139
140 Assertions.assertEquals(0.137,
141 Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
142 finalOrbit.getPVCoordinates().getVelocity()),
143 1.0e-3);
144 Assertions.assertEquals(125.2, finalOrbit.getOrbit().getA() - initialOrbit.getA(), 0.1);
145
146 }
147
148 @Test
149 public void sameDateKeplerian() {
150
151
152
153 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
154 Orbit initialOrbit = new KeplerianOrbit(7209668.0, 0.5e-4, 1.7, 2.1, 2.9,
155 6.2, PositionAngleType.TRUE,
156 FramesFactory.getEME2000(), initDate, provider.getMu());
157
158
159
160 EcksteinHechlerPropagator extrapolator =
161 new EcksteinHechlerPropagator(initialOrbit, Propagator.DEFAULT_MASS, provider);
162
163
164
165 SpacecraftState finalOrbit = extrapolator.propagate(initDate);
166
167
168 Assertions.assertEquals(0.0,
169 Vector3D.distance(initialOrbit.getPosition(),
170 finalOrbit.getPosition()),
171 1.3e-6);
172
173
174
175
176
177
178
179
180
181
182 Assertions.assertEquals(0.137,
183 Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
184 finalOrbit.getPVCoordinates().getVelocity()),
185 1.0e-3);
186 Assertions.assertEquals(126.8, finalOrbit.getOrbit().getA() - initialOrbit.getA(), 0.1);
187
188 }
189
190 @Test
191 public void almostSphericalBody() {
192
193
194
195
196 Vector3D position = new Vector3D(3220103., 69623., 6449822.);
197 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
198
199 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
200 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
201 FramesFactory.getEME2000(), initDate, provider.getMu());
202
203
204
205
206
207
208
209 UnnormalizedSphericalHarmonicsProvider kepProvider =
210 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
211 TideSystem.UNKNOWN,
212 new double[][] {
213 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
214 }, new double[][] {
215 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
216 });
217
218
219
220 EcksteinHechlerPropagator extrapolatorAna =
221 new EcksteinHechlerPropagator(initialOrbit, 1000.0, kepProvider);
222 KeplerianPropagator extrapolatorKep = new KeplerianPropagator(initialOrbit);
223
224
225
226 double delta_t = 100.0;
227 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
228
229 Orbit finalOrbitAna = extrapolatorAna.propagate(extrapDate).getOrbit();
230 Orbit finalOrbitKep = extrapolatorKep.propagate(extrapDate).getOrbit();
231
232 Assertions.assertEquals(finalOrbitAna.getDate().durationFrom(extrapDate), 0.0,
233 Utils.epsilonTest);
234
235 Assertions.assertEquals(finalOrbitAna.getA(), finalOrbitKep.getA(), 10
236 * Utils.epsilonTest * finalOrbitKep.getA());
237 Assertions.assertEquals(finalOrbitAna.getEquinoctialEx(), finalOrbitKep.getEquinoctialEx(), Utils.epsilonE
238 * finalOrbitKep.getE());
239 Assertions.assertEquals(finalOrbitAna.getEquinoctialEy(), finalOrbitKep.getEquinoctialEy(), Utils.epsilonE
240 * finalOrbitKep.getE());
241 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx(), finalOrbitKep.getHx()),
242 finalOrbitKep.getHx(), Utils.epsilonAngle
243 * FastMath.abs(finalOrbitKep.getI()));
244 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy(), finalOrbitKep.getHy()),
245 finalOrbitKep.getHy(), Utils.epsilonAngle
246 * FastMath.abs(finalOrbitKep.getI()));
247 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv(), finalOrbitKep.getLv()),
248 finalOrbitKep.getLv(), Utils.epsilonAngle
249 * FastMath.abs(finalOrbitKep.getLv()));
250 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE(), finalOrbitKep.getLE()),
251 finalOrbitKep.getLE(), Utils.epsilonAngle
252 * FastMath.abs(finalOrbitKep.getLE()));
253 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM(), finalOrbitKep.getLM()),
254 finalOrbitKep.getLM(), Utils.epsilonAngle
255 * FastMath.abs(finalOrbitKep.getLM()));
256
257 }
258
259 @Test
260 public void propagatedCartesian() {
261
262
263
264 Vector3D position = new Vector3D(3220103., 69623., 6449822.);
265 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
266
267 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
268 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
269 FramesFactory.getEME2000(), initDate, provider.getMu());
270
271
272
273 EcksteinHechlerPropagator extrapolator =
274 new EcksteinHechlerPropagator(initialOrbit,
275 new LofOffset(initialOrbit.getFrame(),
276 LOFType.VNC, RotationOrder.XYZ, 0, 0, 0),
277 provider);
278
279
280
281 double delta_t = 100000.0;
282 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
283
284 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
285
286 Assertions.assertEquals(0.0, finalOrbit.getDate().durationFrom(extrapDate), 1.0e-9);
287
288
289 double LM = finalOrbit.getLE() - finalOrbit.getEquinoctialEx()
290 * FastMath.sin(finalOrbit.getLE()) + finalOrbit.getEquinoctialEy()
291 * FastMath.cos(finalOrbit.getLE());
292
293 Assertions.assertEquals(LM, finalOrbit.getLM(), Utils.epsilonAngle
294 * FastMath.abs(finalOrbit.getLM()));
295
296
297 Assertions.assertEquals(FastMath.tan((finalOrbit.getLE() - finalOrbit.getLv()) / 2.),
298 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit
299 .getEquinoctialEy()), Utils.epsilonAngle);
300
301
302 double deltaM = finalOrbit.getLM() - initialOrbit.getLM();
303 double deltaE = finalOrbit.getLE() - initialOrbit.getLE();
304 double delta = finalOrbit.getEquinoctialEx() * FastMath.sin(finalOrbit.getLE())
305 - initialOrbit.getEquinoctialEx() * FastMath.sin(initialOrbit.getLE())
306 - finalOrbit.getEquinoctialEy() * FastMath.cos(finalOrbit.getLE())
307 + initialOrbit.getEquinoctialEy() * FastMath.cos(initialOrbit.getLE());
308
309 Assertions.assertEquals(deltaM, deltaE - delta, Utils.epsilonAngle
310 * FastMath.abs(deltaE - delta));
311
312
313 double ex = finalOrbit.getEquinoctialEx();
314 double ey = finalOrbit.getEquinoctialEy();
315 double hx = finalOrbit.getHx();
316 double hy = finalOrbit.getHy();
317 double LE = finalOrbit.getLE();
318
319 double ex2 = ex * ex;
320 double ey2 = ey * ey;
321 double hx2 = hx * hx;
322 double hy2 = hy * hy;
323 double h2p1 = 1. + hx2 + hy2;
324 double beta = 1. / (1. + FastMath.sqrt(1. - ex2 - ey2));
325
326 double x3 = -ex + (1. - beta * ey2) * FastMath.cos(LE) + beta * ex * ey
327 * FastMath.sin(LE);
328 double y3 = -ey + (1. - beta * ex2) * FastMath.sin(LE) + beta * ex * ey
329 * FastMath.cos(LE);
330
331 Vector3D U = new Vector3D((1. + hx2 - hy2) / h2p1, (2. * hx * hy) / h2p1,
332 (-2. * hy) / h2p1);
333
334 Vector3D V = new Vector3D((2. * hx * hy) / h2p1, (1. - hx2 + hy2) / h2p1,
335 (2. * hx) / h2p1);
336
337 Vector3D r = new Vector3D(finalOrbit.getA(), (new Vector3D(x3, U, y3, V)));
338
339 Assertions.assertEquals(finalOrbit.getPosition().getNorm(), r.getNorm(),
340 Utils.epsilonTest * r.getNorm());
341
342 }
343
344 @Test
345 public void propagatedKeplerian() {
346
347
348 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
349 Orbit initialOrbit = new KeplerianOrbit(7209668.0, 0.5e-4, 1.7, 2.1, 2.9,
350 6.2, PositionAngleType.TRUE,
351 FramesFactory.getEME2000(), initDate, provider.getMu());
352
353
354
355 EcksteinHechlerPropagator extrapolator =
356 new EcksteinHechlerPropagator(initialOrbit,
357 new LofOffset(initialOrbit.getFrame(),
358 LOFType.VNC, RotationOrder.XYZ, 0, 0, 0),
359 2000.0, provider);
360
361
362
363 double delta_t = 100000.0;
364 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
365
366 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
367
368 Assertions.assertEquals(0.0, finalOrbit.getDate().durationFrom(extrapDate), 1.0e-9);
369
370
371 double LM = finalOrbit.getLE() - finalOrbit.getEquinoctialEx()
372 * FastMath.sin(finalOrbit.getLE()) + finalOrbit.getEquinoctialEy()
373 * FastMath.cos(finalOrbit.getLE());
374
375 Assertions.assertEquals(LM, finalOrbit.getLM(), Utils.epsilonAngle);
376
377
378 Assertions.assertEquals(FastMath.tan((finalOrbit.getLE() - finalOrbit.getLv()) / 2.),
379 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit
380 .getEquinoctialEy()), Utils.epsilonAngle);
381
382
383
384 double deltaM = finalOrbit.getLM() - initialOrbit.getLM();
385 double deltaE = finalOrbit.getLE() - initialOrbit.getLE();
386 double delta = finalOrbit.getEquinoctialEx() * FastMath.sin(finalOrbit.getLE())
387 - initialOrbit.getEquinoctialEx() * FastMath.sin(initialOrbit.getLE())
388 - finalOrbit.getEquinoctialEy() * FastMath.cos(finalOrbit.getLE())
389 + initialOrbit.getEquinoctialEy() * FastMath.cos(initialOrbit.getLE());
390
391 Assertions.assertEquals(deltaM, deltaE - delta, Utils.epsilonAngle
392 * FastMath.abs(deltaE - delta));
393
394
395 double ex = finalOrbit.getEquinoctialEx();
396 double ey = finalOrbit.getEquinoctialEy();
397 double hx = finalOrbit.getHx();
398 double hy = finalOrbit.getHy();
399 double LE = finalOrbit.getLE();
400
401 double ex2 = ex * ex;
402 double ey2 = ey * ey;
403 double hx2 = hx * hx;
404 double hy2 = hy * hy;
405 double h2p1 = 1. + hx2 + hy2;
406 double beta = 1. / (1. + FastMath.sqrt(1. - ex2 - ey2));
407
408 double x3 = -ex + (1. - beta * ey2) * FastMath.cos(LE) + beta * ex * ey
409 * FastMath.sin(LE);
410 double y3 = -ey + (1. - beta * ex2) * FastMath.sin(LE) + beta * ex * ey
411 * FastMath.cos(LE);
412
413 Vector3D U = new Vector3D((1. + hx2 - hy2) / h2p1, (2. * hx * hy) / h2p1,
414 (-2. * hy) / h2p1);
415
416 Vector3D V = new Vector3D((2. * hx * hy) / h2p1, (1. - hx2 + hy2) / h2p1,
417 (2. * hx) / h2p1);
418
419 Vector3D r = new Vector3D(finalOrbit.getA(), (new Vector3D(x3, U, y3, V)));
420
421 Assertions.assertEquals(finalOrbit.getPosition().getNorm(), r.getNorm(),
422 Utils.epsilonTest * r.getNorm());
423
424 }
425
426 @Test
427 public void undergroundOrbit() {
428 Assertions.assertThrows(OrekitException.class, () -> {
429
430 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
431 Vector3D velocity = new Vector3D(-500.0, 800.0, 100.0);
432 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
433 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
434 FramesFactory.getEME2000(), initDate, provider.getMu());
435
436
437 EcksteinHechlerPropagator extrapolator =
438 new EcksteinHechlerPropagator(initialOrbit, provider);
439
440
441
442 double delta_t = 0.0;
443 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
444 extrapolator.propagate(extrapDate);
445 });
446 }
447
448 @Test
449 public void equatorialOrbit() {
450 Assertions.assertThrows(OrekitException.class, () -> {
451 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
452 Orbit initialOrbit = new CircularOrbit(7000000, 1.0e-4, -1.5e-4,
453 0.0, 1.2, 2.3, PositionAngleType.MEAN,
454 FramesFactory.getEME2000(),
455 initDate, provider.getMu());
456
457
458 EcksteinHechlerPropagator extrapolator =
459 new EcksteinHechlerPropagator(initialOrbit, provider);
460
461
462
463 double delta_t = 0.0;
464 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
465 extrapolator.propagate(extrapDate);
466 });
467 }
468
469 @Test
470 public void criticalInclination() {
471 Assertions.assertThrows(OrekitException.class, () -> {
472 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
473 Orbit initialOrbit = new CircularOrbit(new PVCoordinates(new Vector3D(-3862363.8474653554,
474 -3521533.9758022362,
475 4647637.852558916),
476 new Vector3D(65.36170817232278,
477 -6056.563439401233,
478 -4511.1247889782757)),
479 FramesFactory.getEME2000(),
480 initDate, provider.getMu());
481
482
483
484 EcksteinHechlerPropagator extrapolator =
485 new EcksteinHechlerPropagator(initialOrbit, provider);
486
487
488
489 double delta_t = 0.0;
490 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
491 extrapolator.propagate(extrapDate);
492 });
493 }
494
495 @Test
496 public void tooEllipticalOrbit() {
497 Assertions.assertThrows(OrekitException.class, () -> {
498
499 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
500 Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
501 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
502 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
503 FramesFactory.getEME2000(), initDate, provider.getMu());
504
505
506 EcksteinHechlerPropagator extrapolator =
507 new EcksteinHechlerPropagator(initialOrbit, provider);
508
509
510
511 double delta_t = 0.0;
512 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
513 extrapolator.propagate(extrapDate);
514 });
515 }
516
517 @Test
518 public void hyperbolic() {
519 Assertions.assertThrows(OrekitException.class, () -> {
520 KeplerianOrbit hyperbolic =
521 new KeplerianOrbit(-1.0e10, 2, 0, 0, 0, 0, PositionAngleType.TRUE,
522 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
523 EcksteinHechlerPropagator propagator =
524 new EcksteinHechlerPropagator(hyperbolic, provider);
525 propagator.propagate(AbsoluteDate.J2000_EPOCH.shiftedBy(10.0));
526 });
527 }
528
529 @Test
530 public void wrongAttitude() {
531 Assertions.assertThrows(OrekitException.class, () -> {
532 KeplerianOrbit orbit =
533 new KeplerianOrbit(1.0e10, 1.0e-4, 1.0e-2, 0, 0, 0, PositionAngleType.TRUE,
534 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
535 AttitudeProvider wrongLaw = new AttitudeProvider() {
536 public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
537 throw new OrekitException(new DummyLocalizable("gasp"), new RuntimeException());
538 }
539 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(FieldPVCoordinatesProvider<T> pvProv,
540 FieldAbsoluteDate<T> date, Frame frame)
541 {
542 throw new OrekitException(new DummyLocalizable("gasp"), new RuntimeException());
543 }
544 };
545 EcksteinHechlerPropagator propagator =
546 new EcksteinHechlerPropagator(orbit, wrongLaw, provider);
547 propagator.propagate(AbsoluteDate.J2000_EPOCH.shiftedBy(10.0));
548 });
549 }
550
551 @Test
552 public void testAcceleration() {
553 final KeplerianOrbit orbit =
554 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
555 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
556 EcksteinHechlerPropagator propagator =
557 new EcksteinHechlerPropagator(orbit, provider);
558 AbsoluteDate target = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
559 List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
560 for (double dt : Arrays.asList(-0.5, 0.0, 0.5)) {
561 sample.add(propagator.propagate(target.shiftedBy(dt)).getPVCoordinates());
562 }
563
564
565 final TimeInterpolator<TimeStampedPVCoordinates> interpolator =
566 new TimeStampedPVCoordinatesHermiteInterpolator(sample.size(), CartesianDerivativesFilter.USE_P);
567
568 TimeStampedPVCoordinates interpolated = interpolator.interpolate(target, sample);
569 Vector3D computedP = sample.get(1).getPosition();
570 Vector3D computedV = sample.get(1).getVelocity();
571 Vector3D referenceP = interpolated.getPosition();
572 Vector3D referenceV = interpolated.getVelocity();
573 Vector3D computedA = sample.get(1).getAcceleration();
574 Vector3D referenceA = interpolated.getAcceleration();
575 final CircularOrbit propagated = (CircularOrbit) OrbitType.CIRCULAR.convertType(propagator.propagateOrbit(target));
576 final CircularOrbit keplerian =
577 new CircularOrbit(propagated.getA(),
578 propagated.getCircularEx(),
579 propagated.getCircularEy(),
580 propagated.getI(),
581 propagated.getRightAscensionOfAscendingNode(),
582 propagated.getAlphaM(), PositionAngleType.MEAN,
583 propagated.getFrame(),
584 propagated.getDate(),
585 propagated.getMu());
586 Vector3D keplerianP = keplerian.getPosition();
587 Vector3D keplerianV = keplerian.getPVCoordinates().getVelocity();
588 Vector3D keplerianA = keplerian.getPVCoordinates().getAcceleration();
589
590
591 Assertions.assertEquals(0.0, Vector3D.distance(referenceP, computedP), 1.0e-15);
592 Assertions.assertEquals(0.0, Vector3D.distance(referenceP, keplerianP), 4.0e-9);
593
594
595
596 double computationErrorV = Vector3D.distance(referenceV, computedV);
597 double nonKeplerianEffectV = Vector3D.distance(referenceV, keplerianV);
598 Assertions.assertEquals(nonKeplerianEffectV, computationErrorV, 2.0e-12);
599 Assertions.assertEquals(2.2e-4, computationErrorV, 3.0e-6);
600
601
602
603
604 double computationErrorA = Vector3D.distance(referenceA, computedA);
605 double nonKeplerianEffectA = Vector3D.distance(referenceA, keplerianA);
606 Assertions.assertEquals(1.36e-7, computationErrorA, 1.0e-9);
607 Assertions.assertEquals(6.37e-3, nonKeplerianEffectA, 7.0e-6);
608 Assertions.assertTrue(computationErrorA < nonKeplerianEffectA / 40000);
609
610 }
611
612 @Test
613 public void ascendingNode() {
614 final KeplerianOrbit orbit =
615 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
616 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
617 EcksteinHechlerPropagator propagator =
618 new EcksteinHechlerPropagator(orbit, provider);
619 NodeDetector detector = new NodeDetector(orbit, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
620 Assertions.assertTrue(FramesFactory.getITRF(IERSConventions.IERS_2010, true) == detector.getFrame());
621 propagator.addEventDetector(detector);
622 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
623 SpacecraftState propagated = propagator.propagate(farTarget);
624 PVCoordinates pv = propagated.getPVCoordinates(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
625 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 3500.0);
626 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 4000.0);
627 Assertions.assertEquals(0, pv.getPosition().getZ(), 1.0e-6);
628 Assertions.assertTrue(pv.getVelocity().getZ() > 0);
629 Collection<EventDetector> detectors = propagator.getEventDetectors();
630 Assertions.assertEquals(1, detectors.size());
631 propagator.clearEventsDetectors();
632 Assertions.assertEquals(0, propagator.getEventDetectors().size());
633 }
634
635 @Test
636 public void stopAtTargetDate() {
637 final KeplerianOrbit orbit =
638 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
639 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
640 EcksteinHechlerPropagator propagator =
641 new EcksteinHechlerPropagator(orbit, provider);
642 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
643 propagator.addEventDetector(new NodeDetector(orbit, itrf).withHandler(new ContinueOnEvent()));
644 AbsoluteDate farTarget = orbit.getDate().shiftedBy(10000.0);
645 SpacecraftState propagated = propagator.propagate(farTarget);
646 Assertions.assertEquals(0.0, FastMath.abs(farTarget.durationFrom(propagated.getDate())), 1.0e-3);
647 }
648
649 @Test
650 public void perigee() {
651 final KeplerianOrbit orbit =
652 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
653 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
654 EcksteinHechlerPropagator propagator =
655 new EcksteinHechlerPropagator(orbit, provider);
656 propagator.addEventDetector(new ApsideDetector(orbit));
657 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
658 SpacecraftState propagated = propagator.propagate(farTarget);
659 Vector3D pos = propagated.getPosition(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
660 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 3000.0);
661 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 3500.0);
662 Assertions.assertEquals(orbit.getA() * (1.0 - orbit.getE()), pos.getNorm(), 410);
663 }
664
665 @Test
666 public void date() {
667 final KeplerianOrbit orbit =
668 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
669 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
670 EcksteinHechlerPropagator propagator =
671 new EcksteinHechlerPropagator(orbit, provider);
672 final AbsoluteDate stopDate = AbsoluteDate.J2000_EPOCH.shiftedBy(500.0);
673 propagator.addEventDetector(new DateDetector(stopDate));
674 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
675 SpacecraftState propagated = propagator.propagate(farTarget);
676 Assertions.assertEquals(0, stopDate.durationFrom(propagated.getDate()), 1.0e-10);
677 }
678
679 @Test
680 public void fixedStep() {
681 final KeplerianOrbit orbit =
682 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
683 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
684 EcksteinHechlerPropagator propagator =
685 new EcksteinHechlerPropagator(orbit, provider);
686 final double step = 100.0;
687 propagator.setStepHandler(step, new OrekitFixedStepHandler() {
688 private AbsoluteDate previous;
689 public void handleStep(SpacecraftState currentState) {
690 if (previous != null) {
691 Assertions.assertEquals(step, currentState.getDate().durationFrom(previous), 1.0e-10);
692 }
693 previous = currentState.getDate();
694 }
695 });
696 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
697 propagator.propagate(farTarget);
698 }
699
700 @Test
701 public void setting() {
702 final KeplerianOrbit orbit =
703 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
704 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
705 EcksteinHechlerPropagator propagator =
706 new EcksteinHechlerPropagator(orbit, provider);
707 final OneAxisEllipsoid earthShape =
708 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
709 final TopocentricFrame topo =
710 new TopocentricFrame(earthShape, new GeodeticPoint(0.389, -2.962, 0), null);
711 ElevationDetector detector = new ElevationDetector(60, 1.0e-9, topo).withConstantElevation(0.09);
712 Assertions.assertEquals(0.09, detector.getMinElevation(), 1.0e-12);
713 Assertions.assertTrue(topo == detector.getTopocentricFrame());
714 propagator.addEventDetector(detector);
715 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
716 SpacecraftState propagated = propagator.propagate(farTarget);
717 final double elevation = topo.getTrackingCoordinates(propagated.getPosition(),
718 propagated.getFrame(),
719 propagated.getDate()).
720 getElevation();
721 final double zVelocity = propagated.getPVCoordinates(topo).getVelocity().getZ();
722 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 7800.0);
723 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 7900.0,"Incorrect value " + farTarget.durationFrom(propagated.getDate()) + " !< 7900");
724 Assertions.assertEquals(0.09, elevation, 1.0e-11);
725 Assertions.assertTrue(zVelocity < 0);
726 }
727
728 @Test
729 public void testInitializationCorrectness()
730 throws IOException {
731
732
733 AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(154.);
734 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
735 Frame eme2000 = FramesFactory.getEME2000();
736 Vector3D pole = itrf.getStaticTransformTo(eme2000, date).transformVector(Vector3D.PLUS_K);
737 Frame poleAligned = new Frame(FramesFactory.getEME2000(),
738 new Transform(date, new Rotation(pole, Vector3D.PLUS_K)),
739 "pole aligned", true);
740 CircularOrbit initial = new CircularOrbit(7208669.8179538045, 1.3740461966386876E-4, -3.2364250248363356E-5,
741 FastMath.toRadians(97.40236024565775),
742 FastMath.toRadians(166.15873160992115),
743 FastMath.toRadians(90.1282370098961), PositionAngleType.MEAN,
744 poleAligned, date, provider.getMu());
745
746
747 EcksteinHechlerPropagator defaultEH = new EcksteinHechlerPropagator(initial, provider);
748
749
750
751 CircularOrbit defaultOrbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(defaultEH.propagateOrbit(initial.getDate()));
752 Assertions.assertEquals(267.4, defaultOrbit.getA() - initial.getA(), 0.1);
753
754
755 Assertions.assertEquals(0.0,
756 Vector3D.distance(defaultOrbit.getPosition(),
757 initial.getPosition()),
758 2.0e-6);
759
760
761
762 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.1).getTolerances(initial, OrbitType.CIRCULAR);
763 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
764 integrator.setInitialStepSize(60);
765 NumericalPropagator num = new NumericalPropagator(integrator);
766 num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider)));
767 num.setInitialState(new SpacecraftState(initial));
768 num.setOrbitType(OrbitType.CIRCULAR);
769
770
771 PropagatorConverter converter =
772 new FiniteDifferencePropagatorConverter(new EcksteinHechlerPropagatorBuilder(initial,
773 provider,
774 PositionAngleType.TRUE,
775 1.0),
776 1.0e-6, 100);
777 EcksteinHechlerPropagator fittedEH =
778 (EcksteinHechlerPropagator) converter.convert(num, 3 * initial.getKeplerianPeriod(), 300);
779
780
781
782 CircularOrbit fittedOrbit = (CircularOrbit) OrbitType.CIRCULAR.convertType(fittedEH.propagateOrbit(initial.getDate()));
783 Assertions.assertEquals(0.623, defaultOrbit.getA() - fittedOrbit.getA(), 0.1);
784
785
786
787
788 Assertions.assertEquals(58.0,
789 Vector3D.distance(defaultOrbit.getPosition(),
790 fittedOrbit.getPosition()),
791 0.1);
792
793 }
794
795 @Test
796 public void testIssue504() {
797
798 final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
799 final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
800 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2018, 07, 15), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
801 final SpacecraftState initialState = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
802 FramesFactory.getEME2000(),
803 initDate,
804 provider.getMu()));
805
806
807 final List<DSSTForceModel> models = new ArrayList<>();
808 models.add(new DSSTZonal(provider));
809 final SpacecraftState meanState = DSSTPropagator.computeMeanState(initialState, DEFAULT_LAW, models);
810
811
812 final EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(meanState.getOrbit(), provider, PropagationType.MEAN);
813 final SpacecraftState finalState = propagator.propagate(initDate);
814
815
816 final Orbit initialOrbit = initialState.getOrbit();
817 final Orbit finalOrbit = finalState.getOrbit();
818 Assertions.assertEquals(initialOrbit.getA(), finalOrbit.getA(), 18.0);
819 Assertions.assertEquals(initialOrbit.getEquinoctialEx(), finalOrbit.getEquinoctialEx(), 1.0e-6);
820 Assertions.assertEquals(initialOrbit.getEquinoctialEy(), finalOrbit.getEquinoctialEy(), 5.0e-6);
821 Assertions.assertEquals(initialOrbit.getHx(), finalOrbit.getHx(), 1.0e-6);
822 Assertions.assertEquals(initialOrbit.getHy(), finalOrbit.getHy(), 2.0e-6);
823 Assertions.assertEquals(0.0,
824 Vector3D.distance(initialState.getPosition(),
825 finalState.getPosition()),
826 11.4);
827 Assertions.assertEquals(0.0,
828 Vector3D.distance(initialState.getPVCoordinates().getVelocity(),
829 finalState.getPVCoordinates().getVelocity()),
830 4.2e-2);
831 }
832
833 @Test
834 public void testIssue504Bis() {
835
836 final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
837 final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
838 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2018, 07, 15), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
839 final SpacecraftState initialState = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
840 FramesFactory.getEME2000(),
841 initDate,
842 provider.getMu()));
843
844
845 final List<DSSTForceModel> models = new ArrayList<>();
846 models.add(new DSSTZonal(provider));
847 final SpacecraftState meanState = DSSTPropagator.computeMeanState(initialState, DEFAULT_LAW, models);
848
849
850 final EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(meanState.getOrbit(), DEFAULT_LAW, 458.6, provider, PropagationType.MEAN);
851 final SpacecraftState finalState = propagator.propagate(initDate);
852
853
854 final Orbit initialOrbit = initialState.getOrbit();
855 final Orbit finalOrbit = finalState.getOrbit();
856 Assertions.assertEquals(initialOrbit.getA(), finalOrbit.getA(), 18.0);
857 Assertions.assertEquals(initialOrbit.getEquinoctialEx(), finalOrbit.getEquinoctialEx(), 1.0e-6);
858 Assertions.assertEquals(initialOrbit.getEquinoctialEy(), finalOrbit.getEquinoctialEy(), 5.0e-6);
859 Assertions.assertEquals(initialOrbit.getHx(), finalOrbit.getHx(), 1.0e-6);
860 Assertions.assertEquals(initialOrbit.getHy(), finalOrbit.getHy(), 2.0e-6);
861 Assertions.assertEquals(0.0,
862 Vector3D.distance(initialState.getPosition(),
863 finalState.getPosition()),
864 11.4);
865 Assertions.assertEquals(0.0,
866 Vector3D.distance(initialState.getPVCoordinates().getVelocity(),
867 finalState.getPVCoordinates().getVelocity()),
868 4.2e-2);
869 }
870
871 @Test
872 public void testMeanOrbit() {
873 final KeplerianOrbit initialOsculating =
874 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
875 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
876 provider.getMu());
877 final UnnormalizedSphericalHarmonics ush = provider.onDate(initialOsculating.getDate());
878
879
880
881 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.1).getTolerances(initialOsculating, OrbitType.CIRCULAR);
882 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
883 integrator.setInitialStepSize(60);
884 NumericalPropagator num = new NumericalPropagator(integrator);
885 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
886 num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider)));
887 num.setInitialState(new SpacecraftState(initialOsculating));
888 num.setOrbitType(OrbitType.CIRCULAR);
889 num.setPositionAngleType(initialOsculating.getCachedPositionAngleType());
890 final StorelessUnivariateStatistic oscMin = new Min();
891 final StorelessUnivariateStatistic oscMax = new Max();
892 final StorelessUnivariateStatistic meanMin = new Min();
893 final StorelessUnivariateStatistic meanMax = new Max();
894 num.getMultiplexer().add(60, state -> {
895 final Orbit osc = state.getOrbit();
896 oscMin.increment(osc.getA());
897 oscMax.increment(osc.getA());
898
899 final Orbit mean = EcksteinHechlerPropagator.computeMeanOrbit(state.getOrbit(), provider, ush);
900 meanMin.increment(mean.getA());
901 meanMax.increment(mean.getA());
902 });
903 num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
904
905 Assertions.assertEquals(3190.029, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
906 Assertions.assertEquals( 49.638, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
907
908 }
909
910 private static double tangLEmLv(double Lv, double ex, double ey) {
911
912 return (ey * FastMath.cos(Lv) - ex * FastMath.sin(Lv))
913 / (1 + ex * FastMath.cos(Lv) + ey * FastMath.sin(Lv) + FastMath.sqrt(1 - ex * ex
914 - ey * ey));
915 }
916
917 @BeforeEach
918 public void setUp() {
919 Utils.setDataRoot("regular-data");
920 double mu = 3.9860047e14;
921 double ae = 6.378137e6;
922 double[][] cnm = new double[][] {
923 { 0 }, { 0 }, { -1.08263e-3 }, { 2.54e-6 }, { 1.62e-6 }, { 2.3e-7 }, { -5.5e-7 }
924 };
925 double[][] snm = new double[][] {
926 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
927 };
928 provider = GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN, cnm, snm);
929 }
930
931 @AfterEach
932 public void tearDown() {
933 provider = null;
934 }
935
936 private UnnormalizedSphericalHarmonicsProvider provider;
937
938 }