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.hamcrest.MatcherAssert;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.exception.DummyLocalizable;
22 import org.hipparchus.exception.LocalizedCoreFormats;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.linear.MatrixUtils;
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.MathUtils;
27 import org.junit.jupiter.api.AfterEach;
28 import org.junit.jupiter.api.Assertions;
29 import org.junit.jupiter.api.BeforeEach;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.OrekitMatchers;
32 import org.orekit.TestUtils;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.Attitude;
35 import org.orekit.attitudes.AttitudeProvider;
36 import org.orekit.attitudes.FieldAttitude;
37 import org.orekit.attitudes.LofOffset;
38 import org.orekit.bodies.BodyShape;
39 import org.orekit.bodies.GeodeticPoint;
40 import org.orekit.bodies.OneAxisEllipsoid;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.forces.maneuvers.ImpulseManeuver;
43 import org.orekit.frames.Frame;
44 import org.orekit.frames.FramesFactory;
45 import org.orekit.frames.LOFType;
46 import org.orekit.frames.TopocentricFrame;
47 import org.orekit.orbits.*;
48 import org.orekit.propagation.*;
49 import org.orekit.propagation.events.*;
50 import org.orekit.propagation.events.handlers.ContinueOnEvent;
51 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
52 import org.orekit.propagation.sampling.OrekitStepHandler;
53 import org.orekit.propagation.sampling.OrekitStepInterpolator;
54 import org.orekit.time.*;
55 import org.orekit.utils.*;
56
57 import java.util.ArrayList;
58 import java.util.List;
59
60
61 public class KeplerianPropagatorTest {
62
63
64 private double mu;
65
66 @Test
67 void testClearMatricesComputation() {
68
69 final KeplerianPropagator propagator = new KeplerianPropagator(TestUtils.getDefaultOrbit(AbsoluteDate.ARBITRARY_EPOCH));
70 propagator.setupMatricesComputation("stm", MatrixUtils.createRealIdentityMatrix(6), null);
71
72 propagator.clearMatricesComputation();
73
74 Assertions.assertTrue(propagator.getAdditionalDataProviders().isEmpty());
75 }
76
77
78
79
80
81 @Test
82 public void testPropagationDate() {
83
84 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
85
86 AbsoluteDate target =
87 initDate.shiftedBy(20.0).shiftedBy(FastMath.ulp(20.0) / 4);
88 Orbit ic = new KeplerianOrbit(6378137 + 500e3, 1e-3, 0, 0, 0, 0,
89 PositionAngleType.TRUE, FramesFactory.getGCRF(), initDate, mu);
90 Propagator propagator = new KeplerianPropagator(ic);
91
92
93 SpacecraftState actual = propagator.propagate(target);
94
95
96 Assertions.assertEquals(target, actual.getDate());
97 }
98
99 @Test
100 public void testEphemerisModeWithHandler() {
101
102 AbsoluteDate initDate = AbsoluteDate.GPS_EPOCH;
103 Orbit ic = new KeplerianOrbit(6378137 + 500e3, 1e-3, 0, 0, 0, 0,
104 PositionAngleType.TRUE, FramesFactory.getGCRF(), initDate, mu);
105 Propagator propagator = new KeplerianPropagator(ic);
106 AbsoluteDate end = initDate.shiftedBy(90 * 60);
107
108
109 final List<SpacecraftState> states = new ArrayList<>();
110 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
111 propagator.setStepHandler(interpolator -> {
112 states.add(interpolator.getCurrentState());
113 states.add(interpolator.getPreviousState());
114 });
115 propagator.propagate(end);
116 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
117
118
119 Assertions.assertTrue(states.size() > 1);
120 for (SpacecraftState state : states) {
121 PVCoordinates actual =
122 ephemeris.propagate(state.getDate()).getPVCoordinates();
123 MatcherAssert.assertThat(actual, OrekitMatchers.pvIs(state.getPVCoordinates()));
124 }
125 }
126
127 @Test
128 public void testAdditionalState() {
129 AbsoluteDate initDate = AbsoluteDate.GPS_EPOCH;
130 Orbit ic = new KeplerianOrbit(6378137 + 500e3, 1e-3, 0, 0, 0, 0, PositionAngleType.TRUE, FramesFactory.getGCRF(), initDate, mu);
131 Propagator propagator = new KeplerianPropagator(ic);
132 SpacecraftState initialState = propagator.getInitialState().addAdditionalData("myState", 4.2);
133 propagator.resetInitialState(initialState);
134 AbsoluteDate end = initDate.shiftedBy(90 * 60);
135 EphemerisGenerator generator = propagator.getEphemerisGenerator();
136 SpacecraftState finalStateKeplerianPropagator = propagator.propagate(end);
137 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
138 SpacecraftState ephemerisInitialState = ephemeris.getInitialState();
139 SpacecraftState finalStateBoundedPropagator = ephemeris.propagate(end);
140 Assertions.assertEquals(4.2, finalStateKeplerianPropagator.getAdditionalState("myState")[0], 1.0e-15);
141 Assertions.assertEquals(4.2, ephemerisInitialState.getAdditionalState("myState")[0], 1.0e-15);
142 Assertions.assertEquals(4.2, finalStateBoundedPropagator.getAdditionalState("myState")[0], 1.0e-15);
143 }
144
145 @Test
146 public void sameDateCartesian() {
147
148
149
150 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
151 Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
152
153 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
154 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
155 FramesFactory.getEME2000(), initDate, mu);
156
157
158
159 KeplerianPropagator extrapolator = new KeplerianPropagator(initialOrbit);
160
161
162
163 double delta_t = 0.0;
164 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
165
166 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
167
168 double a = finalOrbit.getA();
169
170 double n = FastMath.sqrt(finalOrbit.getMu()/FastMath.pow(a, 3));
171
172 Assertions.assertEquals(n*delta_t,
173 finalOrbit.getLM() - initialOrbit.getLM(),
174 Utils.epsilonTest * FastMath.abs(n*delta_t));
175 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbit.getLM(), initialOrbit.getLM()), initialOrbit.getLM(), Utils.epsilonAngle * FastMath.abs(initialOrbit.getLM()));
176
177 Assertions.assertEquals(finalOrbit.getA(), initialOrbit.getA(), Utils.epsilonTest * initialOrbit.getA());
178 Assertions.assertEquals(finalOrbit.getE(), initialOrbit.getE(), Utils.epsilonE * initialOrbit.getE());
179 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbit.getI(), initialOrbit.getI()), initialOrbit.getI(), Utils.epsilonAngle * FastMath.abs(initialOrbit.getI()));
180
181 }
182
183 @Test
184 public void sameDateKeplerian() {
185
186
187 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
188 Orbit initialOrbit = new KeplerianOrbit(7209668.0, 0.5e-4, 1.7, 2.1, 2.9,
189 6.2, PositionAngleType.TRUE,
190 FramesFactory.getEME2000(), initDate, mu);
191
192
193
194 KeplerianPropagator extrapolator = new KeplerianPropagator(initialOrbit);
195
196
197
198 double delta_t = 0.0;
199 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
200
201 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
202
203 double a = finalOrbit.getA();
204
205 double n = FastMath.sqrt(finalOrbit.getMu()/FastMath.pow(a, 3));
206
207 Assertions.assertEquals(n*delta_t,
208 finalOrbit.getLM() - initialOrbit.getLM(),
209 Utils.epsilonTest * FastMath.max(100., FastMath.abs(n*delta_t)));
210 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbit.getLM(), initialOrbit.getLM()), initialOrbit.getLM(), Utils.epsilonAngle * FastMath.abs(initialOrbit.getLM()));
211
212 Assertions.assertEquals(finalOrbit.getA(), initialOrbit.getA(), Utils.epsilonTest * initialOrbit.getA());
213 Assertions.assertEquals(finalOrbit.getE(), initialOrbit.getE(), Utils.epsilonE * initialOrbit.getE());
214 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbit.getI(), initialOrbit.getI()), initialOrbit.getI(), Utils.epsilonAngle * FastMath.abs(initialOrbit.getI()));
215
216 }
217
218 @Test
219 public void propagatedCartesian() {
220
221
222
223 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
224 Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
225 double mu = 3.9860047e14;
226
227 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
228 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
229 FramesFactory.getEME2000(), initDate, mu);
230
231
232
233 KeplerianPropagator extrapolator = new KeplerianPropagator(initialOrbit);
234
235
236
237 double delta_t = 100000.0;
238 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
239
240 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
241
242
243
244 double a = finalOrbit.getA();
245
246 double n = FastMath.sqrt(finalOrbit.getMu()/FastMath.pow(a, 3));
247
248 Assertions.assertEquals(n * delta_t,
249 finalOrbit.getLM() - initialOrbit.getLM(),
250 Utils.epsilonAngle);
251
252
253 double LM = finalOrbit.getLE()
254 - finalOrbit.getEquinoctialEx()*FastMath.sin(finalOrbit.getLE())
255 + finalOrbit.getEquinoctialEy()*FastMath.cos(finalOrbit.getLE());
256
257 Assertions.assertEquals(LM , finalOrbit.getLM() , Utils.epsilonAngle);
258
259
260 Assertions.assertEquals(FastMath.tan((finalOrbit.getLE() - finalOrbit.getLv())/2.),
261 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit.getEquinoctialEy()),
262 Utils.epsilonAngle);
263
264
265
266 double deltaM = finalOrbit.getLM() - initialOrbit.getLM();
267 double deltaE = finalOrbit.getLE() - initialOrbit.getLE();
268 double delta = finalOrbit.getEquinoctialEx() * (FastMath.sin(finalOrbit.getLE()) - FastMath.sin(initialOrbit.getLE()))
269 - finalOrbit.getEquinoctialEy() * (FastMath.cos(finalOrbit.getLE()) - FastMath.cos(initialOrbit.getLE()));
270
271 Assertions.assertEquals(deltaM, deltaE - delta, Utils.epsilonAngle);
272
273
274 Assertions.assertEquals(finalOrbit.getA(), initialOrbit.getA(), Utils.epsilonTest * initialOrbit.getA());
275 Assertions.assertEquals(finalOrbit.getEquinoctialEx(), initialOrbit.getEquinoctialEx(), Utils.epsilonE);
276 Assertions.assertEquals(finalOrbit.getEquinoctialEy(), initialOrbit.getEquinoctialEy(), Utils.epsilonE);
277 Assertions.assertEquals(finalOrbit.getHx(), initialOrbit.getHx(), Utils.epsilonAngle);
278 Assertions.assertEquals(finalOrbit.getHy(), initialOrbit.getHy(), Utils.epsilonAngle);
279
280
281 double ex = finalOrbit.getEquinoctialEx();
282 double ey = finalOrbit.getEquinoctialEy();
283 double hx = finalOrbit.getHx();
284 double hy = finalOrbit.getHy();
285 double LE = finalOrbit.getLE();
286
287 double ex2 = ex*ex;
288 double ey2 = ey*ey;
289 double hx2 = hx*hx;
290 double hy2 = hy*hy;
291 double h2p1 = 1. + hx2 + hy2;
292 double beta = 1. / (1. + FastMath.sqrt(1. - ex2 - ey2));
293
294 double x3 = -ex + (1.- beta*ey2)*FastMath.cos(LE) + beta*ex*ey*FastMath.sin(LE);
295 double y3 = -ey + (1. -beta*ex2)*FastMath.sin(LE) + beta*ex*ey*FastMath.cos(LE);
296
297 Vector3D U = new Vector3D((1. + hx2 - hy2)/ h2p1,
298 (2.*hx*hy)/h2p1,
299 (-2.*hy)/h2p1);
300
301 Vector3D V = new Vector3D((2.*hx*hy)/ h2p1,
302 (1.- hx2+ hy2)/h2p1,
303 (2.*hx)/h2p1);
304
305 Vector3D r = new Vector3D(finalOrbit.getA(), new Vector3D(x3, U, y3, V));
306
307 Assertions.assertEquals(finalOrbit.getPosition().getNorm(), r.getNorm(), Utils.epsilonTest * r.getNorm());
308
309 }
310
311 @Test
312 public void propagatedKeplerian() {
313
314
315
316 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
317 Orbit initialOrbit = new KeplerianOrbit(7209668.0, 0.5e-4, 1.7, 2.1, 2.9,
318 6.2, PositionAngleType.TRUE,
319 FramesFactory.getEME2000(), initDate, mu);
320
321
322
323 KeplerianPropagator extrapolator = new KeplerianPropagator(initialOrbit);
324
325
326
327 double delta_t = 100000.0;
328 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
329
330 Orbit finalOrbit = extrapolator.propagate(extrapDate).getOrbit();
331 Assertions.assertEquals(6092.3362422560844633, finalOrbit.getKeplerianPeriod(), 1.0e-12);
332 Assertions.assertEquals(0.001031326088602888358, finalOrbit.getKeplerianMeanMotion(), 1.0e-16);
333
334
335 double a = finalOrbit.getA();
336
337 double n = FastMath.sqrt(finalOrbit.getMu()/FastMath.pow(a, 3));
338
339 Assertions.assertEquals(n * delta_t,
340 finalOrbit.getLM() - initialOrbit.getLM(),
341 Utils.epsilonAngle);
342
343
344 double LM = finalOrbit.getLE()
345 - finalOrbit.getEquinoctialEx()*FastMath.sin(finalOrbit.getLE())
346 + finalOrbit.getEquinoctialEy()*FastMath.cos(finalOrbit.getLE());
347
348 Assertions.assertEquals(LM , finalOrbit.getLM() , Utils.epsilonAngle);
349
350
351 Assertions.assertEquals(FastMath.tan((finalOrbit.getLE() - finalOrbit.getLv())/2.),
352 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit.getEquinoctialEy()),
353 Utils.epsilonAngle);
354
355
356
357 double deltaM = finalOrbit.getLM() - initialOrbit.getLM();
358 double deltaE = finalOrbit.getLE() - initialOrbit.getLE();
359 double delta = finalOrbit.getEquinoctialEx() * (FastMath.sin(finalOrbit.getLE()) - FastMath.sin(initialOrbit.getLE())) - finalOrbit.getEquinoctialEy() * (FastMath.cos(finalOrbit.getLE()) - FastMath.cos(initialOrbit.getLE()));
360
361 Assertions.assertEquals(deltaM, deltaE - delta, Utils.epsilonAngle);
362
363
364 Assertions.assertEquals(finalOrbit.getA(), initialOrbit.getA(), Utils.epsilonTest * initialOrbit.getA());
365 Assertions.assertEquals(finalOrbit.getEquinoctialEx(), initialOrbit.getEquinoctialEx(), Utils.epsilonE);
366 Assertions.assertEquals(finalOrbit.getEquinoctialEy(), initialOrbit.getEquinoctialEy(), Utils.epsilonE);
367 Assertions.assertEquals(finalOrbit.getHx(), initialOrbit.getHx(), Utils.epsilonAngle);
368 Assertions.assertEquals(finalOrbit.getHy(), initialOrbit.getHy(), Utils.epsilonAngle);
369
370
371 double ex = finalOrbit.getEquinoctialEx();
372 double ey = finalOrbit.getEquinoctialEy();
373 double hx = finalOrbit.getHx();
374 double hy = finalOrbit.getHy();
375 double LE = finalOrbit.getLE();
376
377 double ex2 = ex*ex;
378 double ey2 = ey*ey;
379 double hx2 = hx*hx;
380 double hy2 = hy*hy;
381 double h2p1 = 1. + hx2 + hy2;
382 double beta = 1. / (1. + FastMath.sqrt(1. - ex2 - ey2));
383
384 double x3 = -ex + (1.- beta*ey2)*FastMath.cos(LE) + beta*ex*ey*FastMath.sin(LE);
385 double y3 = -ey + (1. -beta*ex2)*FastMath.sin(LE) + beta*ex*ey*FastMath.cos(LE);
386
387 Vector3D U = new Vector3D((1. + hx2 - hy2)/ h2p1,
388 (2.*hx*hy)/h2p1,
389 (-2.*hy)/h2p1);
390
391 Vector3D V = new Vector3D((2.*hx*hy)/ h2p1,
392 (1.- hx2+ hy2)/h2p1,
393 (2.*hx)/h2p1);
394
395 Vector3D r = new Vector3D(finalOrbit.getA(), new Vector3D(x3, U, y3, V));
396
397 Assertions.assertEquals(finalOrbit.getPosition().getNorm(), r.getNorm(), Utils.epsilonTest * r.getNorm());
398
399 }
400
401 @Test
402 public void wrongAttitude() {
403 Assertions.assertThrows(OrekitException.class, () -> {
404 KeplerianOrbit orbit =
405 new KeplerianOrbit(1.0e10, 1.0e-4, 1.0e-2, 0, 0, 0, PositionAngleType.TRUE,
406 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
407 AttitudeProvider wrongLaw = new AttitudeProvider() {
408 public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
409 throw new OrekitException(new DummyLocalizable("gasp"), new RuntimeException());
410 }
411 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(FieldPVCoordinatesProvider<T> pvProv,
412 FieldAbsoluteDate<T> date, Frame frame)
413 {
414 throw new OrekitException(new DummyLocalizable("gasp"), new RuntimeException());
415 }
416 };
417 KeplerianPropagator propagator = new KeplerianPropagator(orbit, wrongLaw);
418 propagator.propagate(AbsoluteDate.J2000_EPOCH.shiftedBy(10.0));
419 });
420 }
421
422 @Test
423 public void testStepException() {
424 Assertions.assertThrows(OrekitException.class, () -> {
425 final KeplerianOrbit orbit =
426 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
427 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
428 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
429 propagator.getMultiplexer().add(new OrekitStepHandler() {
430 public void handleStep(OrekitStepInterpolator interpolator) {
431 }
432 public void finish(SpacecraftState finalState) {
433 throw new OrekitException((Throwable) null, new DummyLocalizable("dummy error"));
434 }
435 });
436
437 propagator.propagate(orbit.getDate().shiftedBy(-3600));
438 });
439 }
440
441 @Test
442 public void tesWrapedAttitudeException() {
443 Assertions.assertThrows(OrekitException.class, () -> {
444 final KeplerianOrbit orbit =
445 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
446 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
447 KeplerianPropagator propagator = new KeplerianPropagator(orbit,
448 new AttitudeProvider() {
449 public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date,
450 Frame frame)
451 {
452 throw new OrekitException((Throwable) null,
453 new DummyLocalizable("dummy error"));
454 }
455 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(FieldPVCoordinatesProvider<T> pvProv,
456 FieldAbsoluteDate<T> date, Frame frame)
457 {
458 throw new OrekitException((Throwable) null,
459 new DummyLocalizable("dummy error"));
460 }
461 });
462 propagator.propagate(orbit.getDate().shiftedBy(10.09));
463 });
464 }
465
466 @Test
467 public void ascendingNode() {
468 final KeplerianOrbit orbit =
469 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
470 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
471 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
472 propagator.addEventDetector(new NodeDetector(orbit, FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
473 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
474 SpacecraftState propagated = propagator.propagate(farTarget);
475 PVCoordinates pv = propagated.getPVCoordinates(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
476 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 3500.0);
477 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 4000.0);
478 Assertions.assertEquals(0, pv.getPosition().getZ(), 2.0e-6);
479 Assertions.assertTrue(pv.getVelocity().getZ() > 0);
480 }
481
482 @Test
483 public void stopAtTargetDate() {
484 final KeplerianOrbit orbit =
485 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
486 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
487 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
488 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
489 propagator.addEventDetector(new NodeDetector(orbit, itrf).withHandler(new ContinueOnEvent()));
490 AbsoluteDate farTarget = orbit.getDate().shiftedBy(10000.0);
491 SpacecraftState propagated = propagator.propagate(farTarget);
492 Assertions.assertEquals(0.0, FastMath.abs(farTarget.durationFrom(propagated.getDate())), 1.0e-3);
493 }
494
495 @Test
496 public void perigee() {
497 final KeplerianOrbit orbit =
498 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
499 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
500 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
501 propagator.addEventDetector(new ApsideDetector(orbit));
502 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
503 SpacecraftState propagated = propagator.propagate(farTarget);
504 Vector3D pos = propagated.getPosition(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
505 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 3000.0);
506 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 3500.0);
507 Assertions.assertEquals(orbit.getA() * (1.0 - orbit.getE()), pos.getNorm(), 1.0e-6);
508 }
509
510 @Test
511 public void altitude() {
512 final KeplerianOrbit orbit =
513 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
514 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
515 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
516 BodyShape bodyShape =
517 new OneAxisEllipsoid(6378137.0, 1.0 / 298.257222101, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
518 AltitudeDetector detector =
519 new AltitudeDetector(0.05 * orbit.getKeplerianPeriod(),
520 1500000, bodyShape);
521 Assertions.assertEquals(1500000, detector.getAltitude(), 1.0e-12);
522 propagator.addEventDetector(detector);
523 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
524 SpacecraftState propagated = propagator.propagate(farTarget);
525 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 5400.0);
526 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 5500.0);
527 GeodeticPoint gp = bodyShape.transform(propagated.getPosition(),
528 propagated.getFrame(), propagated.getDate());
529 Assertions.assertEquals(1500000, gp.getAltitude(), 0.1);
530 }
531
532 @Test
533 public void date() {
534 final KeplerianOrbit orbit =
535 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
536 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
537 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
538 final AbsoluteDate stopDate = AbsoluteDate.J2000_EPOCH.shiftedBy(500.0);
539 propagator.addEventDetector(new DateDetector(stopDate));
540 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
541 SpacecraftState propagated = propagator.propagate(farTarget);
542 Assertions.assertEquals(0, stopDate.durationFrom(propagated.getDate()), 1.0e-10);
543 }
544
545 @Test
546 public void setting() {
547 final KeplerianOrbit orbit =
548 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
549 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
550 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
551 final OneAxisEllipsoid earthShape =
552 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
553 final TopocentricFrame topo =
554 new TopocentricFrame(earthShape, new GeodeticPoint(0.389, -2.962, 0), null);
555 propagator.addEventDetector(new ElevationDetector(60, AbstractDetector.DEFAULT_THRESHOLD, topo).withConstantElevation(0.09));
556 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
557 SpacecraftState propagated = propagator.propagate(farTarget);
558 final double elevation = topo.
559 getTrackingCoordinates(propagated.getPosition(),
560 propagated.getFrame(),
561 propagated.getDate()).
562 getElevation();
563 final double zVelocity = propagated.getPVCoordinates(topo).getVelocity().getZ();
564 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) > 7800.0);
565 Assertions.assertTrue(farTarget.durationFrom(propagated.getDate()) < 7900.0);
566 Assertions.assertEquals(0.09, elevation, 1.0e-9);
567 Assertions.assertTrue(zVelocity < 0);
568 }
569
570 @Test
571 public void fixedStep() {
572 final KeplerianOrbit orbit =
573 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
574 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
575 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
576 final double step = 100.0;
577 final int[] counter = new int[] {0};
578 propagator.setStepHandler(step, new OrekitFixedStepHandler() {
579 private AbsoluteDate previous;
580 public void handleStep(SpacecraftState currentState) {
581 if (previous != null) {
582 Assertions.assertEquals(step, currentState.getDate().durationFrom(previous), 1.0e-10);
583 }
584
585 PVCoordinates expected = new KeplerianPropagator(orbit)
586 .propagate(currentState.getDate()).getPVCoordinates();
587 MatcherAssert.assertThat(
588 currentState.getPVCoordinates(),
589 OrekitMatchers.pvIs(expected));
590 previous = currentState.getDate();
591 counter[0]++;
592 }
593 });
594 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
595 propagator.propagate(farTarget);
596
597 Assertions.assertEquals(
598 counter[0],
599 (int) (farTarget.durationFrom(orbit.getDate()) / step) + 1);
600 }
601
602 @Test
603 public void variableStep() {
604 final KeplerianOrbit orbit =
605 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
606 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
607 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
608 final double step = orbit.getKeplerianPeriod() / 100;
609 final int[] counter = new int[] {0};
610 propagator.setStepHandler(new OrekitStepHandler() {
611 private AbsoluteDate t = orbit.getDate();
612 @Override
613 public void handleStep(OrekitStepInterpolator interpolator) {
614
615 do {
616 PVCoordinates expected = new KeplerianPropagator(orbit)
617 .propagate(t).getPVCoordinates();
618 MatcherAssert.assertThat(
619 interpolator.getInterpolatedState(t).getPVCoordinates(),
620 OrekitMatchers.pvIs(expected));
621 t = t.shiftedBy(step);
622 counter[0]++;
623 } while (t.compareTo(interpolator.getCurrentState().getDate()) <= 0);
624 }
625 });
626 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
627 propagator.propagate(farTarget);
628
629 Assertions.assertEquals(
630 counter[0],
631 (int) (farTarget.durationFrom(orbit.getDate()) / step) + 1);
632 }
633
634 @Test
635 public void ephemeris() {
636 final KeplerianOrbit orbit =
637 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
638 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
639 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
640 AbsoluteDate farTarget = AbsoluteDate.J2000_EPOCH.shiftedBy(10000.0);
641 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
642 propagator.propagate(farTarget);
643 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
644 Assertions.assertEquals(0.0, ephemeris.getMinDate().durationFrom(orbit.getDate()), 1.0e10);
645 Assertions.assertEquals(0.0, ephemeris.getMaxDate().durationFrom(farTarget), 1.0e10);
646 }
647
648 @Test
649 public void testIssue14() {
650 AbsoluteDate initialDate = AbsoluteDate.J2000_EPOCH;
651 final KeplerianOrbit initialOrbit =
652 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
653 FramesFactory.getEME2000(), initialDate, 3.986004415e14);
654 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);
655
656 propagator.getEphemerisGenerator();
657 propagator.propagate(initialDate.shiftedBy(initialOrbit.getKeplerianPeriod()));
658 PVCoordinates pv1 = propagator.getPVCoordinates(initialDate, FramesFactory.getEME2000());
659
660 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
661 propagator.propagate(initialDate.shiftedBy(initialOrbit.getKeplerianPeriod()));
662 PVCoordinates pv2 = generator.getGeneratedEphemeris().getPVCoordinates(initialDate, FramesFactory.getEME2000());
663
664 Assertions.assertEquals(0.0, pv1.getPosition().subtract(pv2.getPosition()).getNorm(), 1.0e-15);
665 Assertions.assertEquals(0.0, pv1.getVelocity().subtract(pv2.getVelocity()).getNorm(), 1.0e-15);
666
667 }
668
669 @Test
670 public void testIssue107() {
671 final TimeScale utc = TimeScalesFactory.getUTC();
672 final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
673 final Vector3D velocity = new Vector3D(505.848, 942.781, 7435.922);
674 final AbsoluteDate date = new AbsoluteDate(2003, 9, 16, utc);
675 final Orbit orbit = new CircularOrbit(new PVCoordinates(position, velocity),
676 FramesFactory.getEME2000(), date, mu);
677
678 Propagator propagator = new KeplerianPropagator(orbit) {
679 AbsoluteDate lastDate = AbsoluteDate.PAST_INFINITY;
680
681 public SpacecraftState basicPropagate(final AbsoluteDate date) {
682 if (date.compareTo(lastDate) < 0) {
683 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
684 "no backward propagation allowed");
685 }
686 lastDate = date;
687 return super.basicPropagate(date);
688 }
689 };
690
691 SpacecraftState finalState = propagator.propagate(date.shiftedBy(3600.0));
692 Assertions.assertEquals(3600.0, finalState.getDate().durationFrom(date), 1.0e-15);
693
694 }
695
696 @Test
697 public void testMu() {
698 final KeplerianOrbit orbit1 =
699 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
700 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
701 Constants.WGS84_EARTH_MU);
702 final KeplerianOrbit orbit2 =
703 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
704 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
705 Constants.EIGEN5C_EARTH_MU);
706 final AbsoluteDate target = orbit1.getDate().shiftedBy(10000.0);
707 PVCoordinates pv1 = new KeplerianPropagator(orbit1).propagate(target).getPVCoordinates();
708 PVCoordinates pv2 = new KeplerianPropagator(orbit2).propagate(target).getPVCoordinates();
709 PVCoordinates pvWithMu1 = new KeplerianPropagator(orbit2, orbit1.getMu()).propagate(target).getPVCoordinates();
710 Assertions.assertEquals(0.026054, Vector3D.distance(pv1.getPosition(), pv2.getPosition()), 1.0e-6);
711 Assertions.assertEquals(0.0, Vector3D.distance(pv1.getPosition(), pvWithMu1.getPosition()), 1.0e-15);
712 }
713
714 @Test
715 public void testResetStateForward() {
716 final Frame eme2000 = FramesFactory.getEME2000();
717 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
718 new TimeComponents(14, 0, 0),
719 TimeScalesFactory.getUTC());
720 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
721 eme2000,
722 date, Constants.EIGEN5C_EARTH_MU);
723 final KeplerianPropagator propagator = new KeplerianPropagator(orbit, new LofOffset(eme2000, LOFType.LVLH));
724
725
726 final AbsoluteDate maneuverDate = date.shiftedBy(1000.0);
727 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
728 new Vector3D(0.0, 0.0, -100.0),
729 350.0));
730
731 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
732 propagator.setStepHandler(60.0, state -> {
733 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
734 if (state.getDate().isBefore(maneuverDate)) {
735 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
736 } else {
737 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
738 }
739 });
740
741 propagator.propagate(orbit.getDate().shiftedBy(1500.0));
742
743 }
744
745 @Test
746 public void testResetStateBackward() {
747 final Frame eme2000 = FramesFactory.getEME2000();
748 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
749 new TimeComponents(14, 0, 0),
750 TimeScalesFactory.getUTC());
751 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
752 eme2000,
753 date, Constants.EIGEN5C_EARTH_MU);
754 final KeplerianPropagator propagator = new KeplerianPropagator(orbit, new LofOffset(eme2000, LOFType.LVLH));
755
756
757 final AbsoluteDate maneuverDate = date.shiftedBy(-1000.0);
758 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
759 new Vector3D(0.0, 0.0, -100.0),
760 350.0));
761
762 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
763 propagator.setStepHandler(60.0, state -> {
764 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
765 if (state.getDate().isAfter(maneuverDate)) {
766 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
767 } else {
768 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
769 }
770 });
771
772 propagator.propagate(orbit.getDate().shiftedBy(-1500.0));
773
774 }
775
776 @Test
777 public void testNoDerivatives() {
778 for (OrbitType type : OrbitType.values()) {
779
780
781 final AbsoluteDate date = new AbsoluteDate(2003, 9, 16, TimeScalesFactory.getUTC());
782 final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
783 final Vector3D velocity = new Vector3D(505.848, 942.781, 7435.922);
784 final Vector3D keplerAcceleration = new Vector3D(-mu / position.getNormSq(), position.normalize());
785 final Vector3D nonKeplerAcceleration = new Vector3D(0.001, 0.002, 0.003);
786 final Vector3D acceleration = keplerAcceleration.add(nonKeplerAcceleration);
787 final TimeStampedPVCoordinates pva = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
788 final Orbit initial = type.convertType(new CartesianOrbit(pva, FramesFactory.getEME2000(), mu));
789 Assertions.assertEquals(type, initial.getType());
790
791
792 checkDerivatives(initial, true);
793
794 KeplerianPropagator propagator = new KeplerianPropagator(initial);
795 Assertions.assertEquals(type, propagator.getInitialState().getOrbit().getType());
796
797
798 checkDerivatives(propagator.getInitialState().getOrbit(), false);
799
800 PVCoordinates initPV = propagator.getInitialState().getOrbit().getPVCoordinates();
801 Assertions.assertEquals(nonKeplerAcceleration.getNorm(), Vector3D.distance(acceleration, initPV.getAcceleration()), 2.0e-15);
802 Assertions.assertEquals(0.0,
803 Vector3D.distance(keplerAcceleration, initPV.getAcceleration()),
804 4.0e-15);
805
806 double dt = 0.2 * initial.getKeplerianPeriod();
807 Orbit orbit = propagator.propagateOrbit(initial.getDate().shiftedBy(dt));
808 Assertions.assertEquals(type, orbit.getType());
809
810
811 checkDerivatives(orbit, false);
812
813
814 checkDerivatives(initial.shiftedBy(dt), true);
815
816 }
817 }
818
819 @Test
820 public void testStackableGenerators() {
821 final Frame eme2000 = FramesFactory.getEME2000();
822 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
823 new TimeComponents(14, 0, 0),
824 TimeScalesFactory.getUTC());
825 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
826 eme2000,
827 date, Constants.EIGEN5C_EARTH_MU);
828 final KeplerianPropagator propagator = new KeplerianPropagator(orbit, new LofOffset(eme2000, LOFType.LVLH));
829
830
831 propagator.addAdditionalDataProvider(new DependentGenerator("F", "E"));
832 propagator.addAdditionalDataProvider(new DependentGenerator("B", "A"));
833 propagator.addAdditionalDataProvider(new DependentGenerator("E", "D"));
834 propagator.addAdditionalDataProvider(new DependentGenerator("C", "B"));
835 propagator.addAdditionalDataProvider(new DependentGenerator("A", null));
836 propagator.addAdditionalDataProvider(new DependentGenerator("D", "C"));
837
838 final SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(3600.0));
839 Assertions.assertEquals(1, finalState.getAdditionalState("A").length);
840 Assertions.assertEquals(1.0, finalState.getAdditionalState("A")[0], 1.0e-15);
841 Assertions.assertEquals(1, finalState.getAdditionalState("B").length);
842 Assertions.assertEquals(2.0, finalState.getAdditionalState("B")[0], 1.0e-15);
843 Assertions.assertEquals(1, finalState.getAdditionalState("C").length);
844 Assertions.assertEquals(3.0, finalState.getAdditionalState("C")[0], 1.0e-15);
845 Assertions.assertEquals(1, finalState.getAdditionalState("D").length);
846 Assertions.assertEquals(4.0, finalState.getAdditionalState("D")[0], 1.0e-15);
847 Assertions.assertEquals(1, finalState.getAdditionalState("E").length);
848 Assertions.assertEquals(5.0, finalState.getAdditionalState("E")[0], 1.0e-15);
849 Assertions.assertEquals(1, finalState.getAdditionalState("F").length);
850 Assertions.assertEquals(6.0, finalState.getAdditionalState("F")[0], 1.0e-15);
851
852 }
853
854 @Test
855 public void testCircularDependency() {
856 final Frame eme2000 = FramesFactory.getEME2000();
857 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
858 new TimeComponents(14, 0, 0),
859 TimeScalesFactory.getUTC());
860 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
861 eme2000,
862 date, Constants.EIGEN5C_EARTH_MU);
863 final KeplerianPropagator propagator = new KeplerianPropagator(orbit, new LofOffset(eme2000, LOFType.LVLH));
864
865
866 propagator.addAdditionalDataProvider(new DependentGenerator("F", "E"));
867 propagator.addAdditionalDataProvider(new DependentGenerator("B", "A"));
868 propagator.addAdditionalDataProvider(new DependentGenerator("E", "D"));
869 propagator.addAdditionalDataProvider(new DependentGenerator("C", "B"));
870 propagator.addAdditionalDataProvider(new DependentGenerator("A", null));
871 propagator.addAdditionalDataProvider(new DependentGenerator("D", "F"));
872
873 final SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(3600.0));
874 Assertions.assertEquals(1, finalState.getAdditionalState("A").length);
875 Assertions.assertEquals(1.0, finalState.getAdditionalState("A")[0], 1.0e-15);
876 Assertions.assertEquals(1, finalState.getAdditionalState("B").length);
877 Assertions.assertEquals(2.0, finalState.getAdditionalState("B")[0], 1.0e-15);
878 Assertions.assertEquals(1, finalState.getAdditionalState("C").length);
879 Assertions.assertEquals(3.0, finalState.getAdditionalState("C")[0], 1.0e-15);
880 Assertions.assertFalse(finalState.hasAdditionalData("D"));
881 Assertions.assertFalse(finalState.hasAdditionalData("E"));
882 Assertions.assertFalse(finalState.hasAdditionalData("F"));
883
884 }
885
886 private void checkDerivatives(final Orbit orbit, final boolean expectedDerivatives) {
887 Assertions.assertEquals(expectedDerivatives, orbit.hasNonKeplerianAcceleration());
888 if (!expectedDerivatives) {
889 Assertions.assertEquals(0., orbit.getADot());
890 Assertions.assertEquals(0., orbit.getEDot());
891 Assertions.assertEquals(0., orbit.getIDot());
892 Assertions.assertEquals(0., orbit.getEquinoctialExDot());
893 Assertions.assertEquals(0., orbit.getEquinoctialEyDot());
894 Assertions.assertEquals(0., orbit.getHxDot());
895 Assertions.assertEquals(0., orbit.getHyDot());
896 }
897 }
898
899 private static double tangLEmLv(double Lv, double ex, double ey){
900
901 return (ey*FastMath.cos(Lv) - ex*FastMath.sin(Lv)) /
902 (1 + ex*FastMath.cos(Lv) + ey*FastMath.sin(Lv) + FastMath.sqrt(1 - ex*ex - ey*ey));
903 }
904
905 @BeforeEach
906 public void setUp() {
907 Utils.setDataRoot("regular-data");
908 mu = 3.9860047e14;
909 }
910
911 @AfterEach
912 public void tearDown() {
913 mu = Double.NaN;
914 }
915
916 private static class DependentGenerator implements AdditionalDataProvider<double[]> {
917
918 private final String name;
919 private final String dependency;
920
921 DependentGenerator(final String name, final String dependency) {
922 this.name = name;
923 this.dependency = dependency;
924 }
925
926 public String getName() {
927 return name;
928 }
929
930 public boolean yields(final SpacecraftState state) {
931 return dependency != null && state.getAdditionalDataValues().getEntry(dependency) == null;
932 }
933
934 public double[] getAdditionalData(final SpacecraftState state) {
935 return new double[] {
936 dependency == null ? 1.0 : state.getAdditionalState(dependency)[0] + 1.0
937 };
938 }
939
940 }
941
942 }
943