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