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