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