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