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