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