1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.orbits;
18
19 import org.hamcrest.MatcherAssert;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.Field;
22 import org.hipparchus.analysis.UnivariateFunction;
23 import org.hipparchus.analysis.differentiation.DSFactory;
24 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
25 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
26 import org.hipparchus.complex.Complex;
27 import org.hipparchus.complex.ComplexField;
28 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.hipparchus.linear.FieldMatrixPreservingVisitor;
31 import org.hipparchus.linear.MatrixUtils;
32 import org.hipparchus.util.Binary64;
33 import org.hipparchus.util.Binary64Field;
34 import org.hipparchus.util.FastMath;
35 import org.hipparchus.util.MathArrays;
36 import org.hipparchus.util.MathUtils;
37 import org.junit.jupiter.api.Assertions;
38 import org.junit.jupiter.api.BeforeEach;
39 import org.junit.jupiter.api.Test;
40 import org.junit.jupiter.params.ParameterizedTest;
41 import org.junit.jupiter.params.provider.EnumSource;
42 import org.orekit.Utils;
43 import org.orekit.errors.OrekitIllegalArgumentException;
44 import org.orekit.errors.OrekitMessages;
45 import org.orekit.frames.Frame;
46 import org.orekit.frames.FramesFactory;
47 import org.orekit.frames.Transform;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.FieldAbsoluteDate;
50 import org.orekit.time.TimeScalesFactory;
51 import org.orekit.utils.Constants;
52 import org.orekit.utils.FieldPVCoordinates;
53 import org.orekit.utils.PVCoordinates;
54 import org.orekit.utils.TimeStampedFieldPVCoordinates;
55
56 import java.util.function.Function;
57
58 import static org.orekit.OrekitMatchers.relativelyCloseTo;
59
60
61 class FieldEquinoctialOrbitTest {
62
63
64 private double mu;
65
66 @BeforeEach
67 public void setUp() {
68
69 Utils.setDataRoot("regular-data");
70
71
72 mu = 3.9860047e14;
73 }
74
75 @ParameterizedTest
76 @EnumSource(PositionAngleType.class)
77 void testWithCachedPositionAngleType(final PositionAngleType positionAngleType) {
78
79 final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
80 final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
81 final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
82 final double muEarth = 3.9860047e14;
83 final CartesianOrbit cartesianOrbit = new CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(), AbsoluteDate.ARBITRARY_EPOCH, muEarth);
84 final EquinoctialOrbit equinoctialOrbit = new EquinoctialOrbit(cartesianOrbit);
85 final Binary64Field field = Binary64Field.getInstance();
86 final FieldEquinoctialOrbit<Binary64> fieldOrbit = new FieldEquinoctialOrbit<>(field, equinoctialOrbit);
87
88 final FieldEquinoctialOrbit<Binary64> orbit = fieldOrbit.withCachedPositionAngleType(positionAngleType);
89
90 Assertions.assertEquals(fieldOrbit.getFrame(), orbit.getFrame());
91 Assertions.assertEquals(fieldOrbit.getDate(), orbit.getDate());
92 Assertions.assertEquals(fieldOrbit.getMu(), orbit.getMu());
93 final Vector3D relativePosition = fieldOrbit.getPosition(orbit.getFrame()).subtract(
94 orbit.getPosition()).toVector3D();
95 Assertions.assertEquals(0., relativePosition.getNorm(), 1e-6);
96 Assertions.assertEquals(fieldOrbit.hasNonKeplerianAcceleration(),
97 orbit.hasNonKeplerianAcceleration());
98 }
99
100 @Test
101 void testInFrameNonKeplerian() {
102 testTemplateInFrame(Vector3D.MINUS_J, PositionAngleType.TRUE);
103 }
104
105 @ParameterizedTest
106 @EnumSource(PositionAngleType.class)
107 void testInFrameKeplerian(final PositionAngleType positionAngleType) {
108 testTemplateInFrame(Vector3D.ZERO, positionAngleType);
109 }
110
111 private void testTemplateInFrame(final Vector3D acceleration, final PositionAngleType positionAngleType) {
112
113 final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
114 final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
115 final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
116 final double muEarth = 3.9860047e14;
117 final CartesianOrbit cartesianOrbit = new CartesianOrbit(pvCoordinates, FramesFactory.getEME2000(),
118 AbsoluteDate.ARBITRARY_EPOCH, muEarth);
119 final FieldEquinoctialOrbit<Binary64> fieldOrbit = new FieldEquinoctialOrbit<>(Binary64Field.getInstance(),
120 cartesianOrbit).withCachedPositionAngleType(positionAngleType);
121
122 final FieldEquinoctialOrbit<Binary64> fieldOrbitWithOtherFrame = fieldOrbit.inFrame(FramesFactory.getGCRF());
123
124 Assertions.assertNotEquals(fieldOrbit.getFrame(), fieldOrbitWithOtherFrame.getFrame());
125 Assertions.assertEquals(fieldOrbit.getDate(), fieldOrbitWithOtherFrame.getDate());
126 Assertions.assertEquals(fieldOrbit.getMu(), fieldOrbitWithOtherFrame.getMu());
127 final FieldVector3D<Binary64> relativePosition = fieldOrbit.getPosition(fieldOrbitWithOtherFrame.getFrame()).subtract(
128 fieldOrbitWithOtherFrame.getPosition());
129 Assertions.assertEquals(0., relativePosition.getNorm().getReal(), 1e-6);
130 Assertions.assertEquals(fieldOrbit.hasNonKeplerianAcceleration(),
131 fieldOrbitWithOtherFrame.hasNonKeplerianAcceleration());
132 }
133
134 @Test
135 void testEquinoctialToEquinoctialEll() {
136 doTestEquinoctialToEquinoctialEll(Binary64Field.getInstance());
137 }
138
139 @Test
140 void testEquinoctialToEquinoctialCirc() {
141 doTestEquinoctialToEquinoctialCirc(Binary64Field.getInstance());
142 }
143
144 @Test
145 void testEquinoctialToCartesian() {
146 doTestEquinoctialToCartesian(Binary64Field.getInstance());
147 }
148
149 @Test
150 void testEquinoctialToKeplerian() {
151 doTestEquinoctialToKeplerian(Binary64Field.getInstance());
152 }
153
154 @Test
155 void testNumericalIssue25() {
156 doTestNumericalIssue25(Binary64Field.getInstance());
157 }
158
159 @Test
160 void testAnomaly() {
161 doTestAnomaly(Binary64Field.getInstance());
162 }
163
164 @Test
165 void testPositionVelocityNorms() {
166 doTestPositionVelocityNorms(Binary64Field.getInstance());
167 }
168
169 @Test
170 void testGeometry() {
171 doTestGeometry(Binary64Field.getInstance());
172 }
173
174 @Test
175 void testRadiusOfCurvature() {
176 doTestRadiusOfCurvature(Binary64Field.getInstance());
177 }
178
179 @Test
180 void testSymmetry() {
181 doTestSymmetry(Binary64Field.getInstance());
182 }
183
184 @Test
185 void testJacobianReference() {
186 doTestJacobianReference(Binary64Field.getInstance());
187 }
188
189 @Test
190 void testJacobianFinitedifferences() {
191 doTestJacobianFinitedifferences(Binary64Field.getInstance());
192 }
193
194 @Test
195 void testHyperbolic1() {
196 doTestHyperbolic1(Binary64Field.getInstance());
197 }
198
199 @Test
200 void testHyperbolic2() {
201 doTestHyperbolic2(Binary64Field.getInstance());
202 }
203
204 @Test
205 void testToOrbitWithoutDerivatives() {
206 doTestToOrbitWithoutDerivatives(Binary64Field.getInstance());
207 }
208
209 @Test
210 void testToOrbitWithDerivatives() {
211 doTestToOrbitWithDerivatives(Binary64Field.getInstance());
212 }
213
214 @Test
215 void testDerivativesConversionSymmetry() {
216 doTestDerivativesConversionSymmetry(Binary64Field.getInstance());
217 }
218
219 @Test
220 void testToString() {
221 doTestToString(Binary64Field.getInstance());
222 }
223
224 @Test
225 void testNonInertialFrame() {
226 Assertions.assertThrows(IllegalArgumentException.class,
227 () -> doTestNonInertialFrame(Binary64Field.getInstance()));
228 }
229
230 @Test
231 void testNonKeplerianDerivatives() {
232 doTestNonKeplerianDerivatives(Binary64Field.getInstance());
233 }
234
235 @Test
236 void testPositionAngleDerivatives() {
237 doTestPositionAngleDerivatives(Binary64Field.getInstance());
238 }
239
240 @Test
241 void testEquatorialRetrograde() {
242 doTestEquatorialRetrograde(Binary64Field.getInstance());
243 }
244
245 @Test
246 void testCopyNonKeplerianAcceleration() {
247 doTestCopyNonKeplerianAcceleration(Binary64Field.getInstance());
248 }
249
250 @Test
251 void testNormalize() {
252 doTestNormalize(Binary64Field.getInstance());
253 }
254
255 @Test
256 void testWithKeplerianRates() {
257
258 final ComplexField field = ComplexField.getInstance();
259 final EquinoctialOrbit expectedOrbit = createOrbitTestFromEquinoctialOrbit(true);
260 final FieldEquinoctialOrbit<Complex> fieldOrbit = new FieldEquinoctialOrbit<>(field, expectedOrbit);
261
262 final FieldEquinoctialOrbit<Complex> actualFieldOrbit = fieldOrbit.withKeplerianRates();
263
264 Assertions.assertFalse(actualFieldOrbit.hasNonKeplerianRates());
265 Assertions.assertEquals(fieldOrbit.getMu(), actualFieldOrbit.getMu());
266 Assertions.assertEquals(fieldOrbit.getDate(), actualFieldOrbit.getDate());
267 Assertions.assertEquals(fieldOrbit.getFrame(), actualFieldOrbit.getFrame());
268 Assertions.assertEquals(fieldOrbit.getPosition(), actualFieldOrbit.getPosition());
269 }
270
271 @Test
272 void testFromEquinoctialOrbitWithDerivatives() {
273
274 final ComplexField field = ComplexField.getInstance();
275 final EquinoctialOrbit expectedOrbit = createOrbitTestFromEquinoctialOrbit(true);
276
277 final FieldEquinoctialOrbit<Complex> fieldOrbit = new FieldEquinoctialOrbit<>(field, expectedOrbit);
278
279 compareFieldOrbitToOrbit(fieldOrbit, expectedOrbit);
280 }
281
282 @Test
283 void testFromEquinoctialOrbitWithoutDerivatives() {
284
285 final ComplexField field = ComplexField.getInstance();
286 final EquinoctialOrbit expectedOrbit = createOrbitTestFromEquinoctialOrbit(false);
287
288 final FieldEquinoctialOrbit<Complex> fieldOrbit = new FieldEquinoctialOrbit<>(field, expectedOrbit);
289
290 compareFieldOrbitToOrbit(fieldOrbit, expectedOrbit);
291 }
292
293 private EquinoctialOrbit createOrbitTestFromEquinoctialOrbit(final boolean withDerivatives) {
294 final PositionAngleType positionAngleType = PositionAngleType.TRUE;
295 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
296 final Frame frame = FramesFactory.getEME2000();
297 final double a = 10000.e3;
298 final double ex = 1e-4;
299 final double ey = 1e-3;
300 final double hx = 1.e-2;
301 final double hy = -2.e-3;
302 final double lv = 0.5;
303 if (withDerivatives) {
304 final double derivative = 0.;
305 return new EquinoctialOrbit(a, ex, ey, hx, hy, lv, derivative, derivative, derivative, derivative,
306 derivative, derivative, positionAngleType, frame, date, mu);
307 } else {
308 return new EquinoctialOrbit(a, ex, ey, hx, hy, lv, positionAngleType, frame, date, mu);
309 }
310 }
311
312 @Test
313 void testCoverageCachedPositionAngleTypeWithRates() {
314
315 final double semiMajorAxis = 1e4;
316 final double ex = 0.;
317 final double ey = 0.;
318 final double expectedL = 0.;
319 final double expectedLDot = 0.;
320 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
321 final Binary64Field field = Binary64Field.getInstance();
322 final Binary64 zero = field.getZero();
323
324 for (final PositionAngleType inputPositionAngleType: PositionAngleType.values()) {
325 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
326 final FieldEquinoctialOrbit<Binary64> fieldOrbit = new FieldEquinoctialOrbit<>(
327 zero.newInstance(semiMajorAxis), zero.newInstance(ex), zero.newInstance(ey), zero, zero,
328 zero.newInstance(expectedL), zero, zero, zero, zero, zero, zero.newInstance(expectedLDot),
329 inputPositionAngleType, cachedPositionAngleType,
330 FramesFactory.getGCRF(), new FieldAbsoluteDate<>(field, date), zero.newInstance(mu));
331 Assertions.assertEquals(cachedPositionAngleType, fieldOrbit.getCachedPositionAngleType());
332 Assertions.assertEquals(expectedL, fieldOrbit.getLv().getReal());
333 Assertions.assertEquals(expectedL, fieldOrbit.getLM().getReal());
334 Assertions.assertEquals(expectedL, fieldOrbit.getLE().getReal());
335 Assertions.assertEquals(expectedLDot, fieldOrbit.getLvDot().getReal());
336 Assertions.assertEquals(expectedLDot, fieldOrbit.getLMDot().getReal());
337 Assertions.assertEquals(expectedLDot, fieldOrbit.getLEDot().getReal());
338 }
339 }
340 }
341
342 private <T extends CalculusFieldElement<T>> void compareFieldOrbitToOrbit(final FieldEquinoctialOrbit<T> fieldOrbit,
343 final EquinoctialOrbit orbit) {
344 Assertions.assertEquals(orbit.getFrame(), fieldOrbit.getFrame());
345 Assertions.assertEquals(orbit.getMu(), fieldOrbit.getMu().getReal());
346 Assertions.assertEquals(orbit.getDate(), fieldOrbit.getDate().toAbsoluteDate());
347 Assertions.assertEquals(orbit.getA(), fieldOrbit.getA().getReal());
348 Assertions.assertEquals(orbit.getEquinoctialEx(), fieldOrbit.getEquinoctialEx().getReal());
349 Assertions.assertEquals(orbit.getEquinoctialEy(), fieldOrbit.getEquinoctialEy().getReal());
350 Assertions.assertEquals(orbit.getHx(), fieldOrbit.getHx().getReal());
351 Assertions.assertEquals(orbit.getHy(), fieldOrbit.getHy().getReal());
352 Assertions.assertEquals(orbit.getLv(), fieldOrbit.getLv().getReal());
353 Assertions.assertEquals(orbit.hasNonKeplerianAcceleration(), fieldOrbit.hasNonKeplerianAcceleration());
354 Assertions.assertEquals(orbit.getADot(), fieldOrbit.getADot().getReal());
355 Assertions.assertEquals(orbit.getEquinoctialExDot(), fieldOrbit.getEquinoctialExDot().getReal());
356 Assertions.assertEquals(orbit.getEquinoctialEyDot(), fieldOrbit.getEquinoctialEyDot().getReal());
357 Assertions.assertEquals(orbit.getHxDot(), fieldOrbit.getHxDot().getReal());
358 Assertions.assertEquals(orbit.getHyDot(), fieldOrbit.getHyDot().getReal());
359 Assertions.assertEquals(orbit.getLvDot(), fieldOrbit.getLvDot().getReal());
360
361 }
362
363 @Test
364 void testGetLVersusDouble() {
365
366 final double semiMajorAxis = 1e7;
367 final double ex = 1e-2;
368 final double ey = 1e-3;
369 final double expectedL = 2;
370 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
371 final Binary64Field field = Binary64Field.getInstance();
372 final Binary64 zero = field.getZero();
373
374 for (final PositionAngleType inputPositionAngleType: PositionAngleType.values()) {
375 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
376 final FieldEquinoctialOrbit<Binary64> fieldOrbit = new FieldEquinoctialOrbit<>(
377 zero.newInstance(semiMajorAxis), zero.newInstance(ex), zero.newInstance(ey), zero, zero,
378 zero.newInstance(expectedL), inputPositionAngleType, cachedPositionAngleType,
379 FramesFactory.getGCRF(), new FieldAbsoluteDate<>(field, date), zero.newInstance(mu));
380 final EquinoctialOrbit equinoctialOrbit = fieldOrbit.toOrbit();
381 Assertions.assertEquals(equinoctialOrbit.getLE(), fieldOrbit.getLE().getReal());
382 Assertions.assertEquals(equinoctialOrbit.getLv(), fieldOrbit.getLv().getReal());
383 Assertions.assertEquals(equinoctialOrbit.getLM(), fieldOrbit.getLM().getReal());
384 }
385 }
386 }
387
388 private <T extends CalculusFieldElement<T>> void doTestEquinoctialToEquinoctialEll(Field<T> field) {
389 T zero = field.getZero();
390 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
391
392 T ix = zero.add(1.200e-04);
393 T iy = zero.add(-1.16e-04);
394 T inc = ix.multiply(ix).add(iy.multiply(iy)).divide(4.).sqrt().asin().multiply(2);
395 T hx = inc.divide(2.).tan().multiply(ix).divide(inc.divide(2.).sin().multiply(2));
396 T hy = inc.divide(2.).tan().multiply(iy).divide(inc.divide(2.).sin().multiply(2));
397
398
399 FieldEquinoctialOrbit<T> equi =
400 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add( 0.5), zero.add(-0.5), hx, hy,
401 zero.add(5.300), PositionAngleType.MEAN,
402 FramesFactory.getEME2000(), date, zero.add(mu));
403 FieldVector3D<T> pos = equi.getPosition();
404 FieldVector3D<T> vit = equi.getPVCoordinates().getVelocity();
405
406 FieldPVCoordinates<T> FieldPVCoordinates = new FieldPVCoordinates<>(pos, vit);
407
408 FieldEquinoctialOrbit<T> param = new FieldEquinoctialOrbit<>(FieldPVCoordinates, FramesFactory.getEME2000(), date, zero.add(mu));
409 Assertions.assertEquals(param.getA().getReal(), equi.getA().getReal(), Utils.epsilonTest * equi.getA().getReal());
410 Assertions.assertEquals(param.getEquinoctialEx().getReal(), equi.getEquinoctialEx().getReal(),
411 Utils.epsilonE * FastMath.abs(equi.getE().getReal()));
412 Assertions.assertEquals(param.getEquinoctialEy().getReal(), equi.getEquinoctialEy().getReal(),
413 Utils.epsilonE * FastMath.abs(equi.getE().getReal()));
414 Assertions.assertEquals(param.getHx().getReal(), equi.getHx().getReal(), Utils.epsilonAngle
415 * FastMath.abs(equi.getI().getReal()));
416 Assertions.assertEquals(param.getHy().getReal(), equi.getHy().getReal(), Utils.epsilonAngle
417 * FastMath.abs(equi.getI().getReal()));
418 Assertions.assertEquals(MathUtils.normalizeAngle(param.getLv().getReal(), equi.getLv().getReal()), equi.getLv().getReal(),
419 Utils.epsilonAngle * FastMath.abs(equi.getLv().getReal()));
420
421 }
422
423 private <T extends CalculusFieldElement<T>> void doTestEquinoctialToEquinoctialCirc(Field<T> field) {
424 T zero = field.getZero();
425 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
426
427 T ix = zero.add(1.200e-04);
428 T iy = zero.add(-1.16e-04);
429 T inc = ix.multiply(ix).add(iy.multiply(iy)).divide(4.).sqrt().asin().multiply(2);
430 T hx = inc.divide(2.).tan().multiply(ix).divide(inc.divide(2.).sin().multiply(2));
431 T hy = inc.divide(2.).tan().multiply(iy).divide(inc.divide(2.).sin().multiply(2));
432
433
434 FieldEquinoctialOrbit<T> equiCir =
435 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.1e-10), zero.add(-0.1e-10), hx, hy,
436 zero.add(5.300), PositionAngleType.MEAN,
437 FramesFactory.getEME2000(), date, zero.add(mu));
438 FieldVector3D<T> posCir = equiCir.getPosition();
439 FieldVector3D<T> vitCir = equiCir.getPVCoordinates().getVelocity();
440
441 FieldPVCoordinates<T> FieldPVCoordinates = new FieldPVCoordinates<>(posCir, vitCir);
442
443 FieldEquinoctialOrbit<T> paramCir = new FieldEquinoctialOrbit<>(FieldPVCoordinates, FramesFactory.getEME2000(),
444 date, zero.add(mu));
445 Assertions.assertEquals(paramCir.getA().getReal(), equiCir.getA().getReal(), Utils.epsilonTest
446 * equiCir.getA().getReal());
447 Assertions.assertEquals(paramCir.getEquinoctialEx().getReal(), equiCir.getEquinoctialEx().getReal(),
448 Utils.epsilonEcir * FastMath.abs(equiCir.getE().getReal()));
449 Assertions.assertEquals(paramCir.getEquinoctialEy().getReal(), equiCir.getEquinoctialEy().getReal(),
450 Utils.epsilonEcir * FastMath.abs(equiCir.getE().getReal()));
451 Assertions.assertEquals(paramCir.getHx().getReal(), equiCir.getHx().getReal(), Utils.epsilonAngle
452 * FastMath.abs(equiCir.getI().getReal()));
453 Assertions.assertEquals(paramCir.getHy().getReal(), equiCir.getHy().getReal(), Utils.epsilonAngle
454 * FastMath.abs(equiCir.getI().getReal()));
455 Assertions.assertEquals(MathUtils.normalizeAngle(paramCir.getLv().getReal(), equiCir.getLv().getReal()), equiCir
456 .getLv().getReal(), Utils.epsilonAngle * FastMath.abs(equiCir.getLv().getReal()));
457
458 }
459
460 private <T extends CalculusFieldElement<T>> void doTestEquinoctialToCartesian(Field<T> field) {
461 T zero = field.getZero();
462 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
463
464 T ix = zero.add(1.200e-04);
465 T iy = zero.add(-1.16e-04);
466 T inc = ix.multiply(ix).add(iy.multiply(iy)).divide(4.).sqrt().asin().multiply(2);
467 T hx = inc.divide(2.).tan().multiply(ix).divide(inc.divide(2.).sin().multiply(2));
468 T hy = inc.divide(2.).tan().multiply(iy).divide(inc.divide(2.).sin().multiply(2));
469
470 FieldEquinoctialOrbit<T> equi =
471 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(-7.900e-06), zero.add(1.100e-04), hx, hy,
472 zero.add(5.300), PositionAngleType.MEAN,
473 FramesFactory.getEME2000(), date, zero.add(mu));
474 FieldVector3D<T> pos = equi.getPosition();
475 FieldVector3D<T> vit = equi.getPVCoordinates().getVelocity();
476
477
478 T oneovera = (pos.getNorm().reciprocal().multiply(2)).subtract(vit.getNorm().multiply(vit.getNorm()).divide(mu));
479 Assertions.assertEquals(oneovera.getReal(), 1. / equi.getA().getReal(), 1.0e-7);
480
481 Assertions.assertEquals(0.233745668678733e+08, pos.getX().getReal(), Utils.epsilonTest
482 * FastMath.abs(pos.getX().getReal()));
483 Assertions.assertEquals(-0.350998914352669e+08, pos.getY().getReal(), Utils.epsilonTest
484 * FastMath.abs(pos.getY().getReal()));
485 Assertions.assertEquals(-0.150053723123334e+04, pos.getZ().getReal(), Utils.epsilonTest
486 * FastMath.abs(pos.getZ().getReal()));
487
488 Assertions.assertEquals(2558.7096558809967, vit.getX().getReal(), Utils.epsilonTest
489 * FastMath.abs(vit.getX().getReal()));
490 Assertions.assertEquals(1704.1586039092576, vit.getY().getReal(), Utils.epsilonTest
491 * FastMath.abs(vit.getY().getReal()));
492 Assertions.assertEquals(0.5013093577879, vit.getZ().getReal(), Utils.epsilonTest
493 * FastMath.abs(vit.getZ().getReal()));
494 }
495
496 private <T extends CalculusFieldElement<T>> void doTestEquinoctialToKeplerian(Field<T> field) {
497 T zero = field.getZero();
498 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
499
500 T ix = zero.add(1.200e-04);
501 T iy = zero.add(-1.16e-04);
502 T inc = ix.multiply(ix).add(iy.multiply(iy)).divide(4.).sqrt().asin().multiply(2);
503 T hx = inc.divide(2.).tan().multiply(ix).divide(inc.divide(2.).sin().multiply(2));
504 T hy = inc.divide(2.).tan().multiply(iy).divide(inc.divide(2.).sin().multiply(2));
505
506 FieldEquinoctialOrbit<T> equi =
507 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(-7.900e-06), zero.add(1.100e-04), hx, hy,
508 zero.add(5.300), PositionAngleType.MEAN,
509 FramesFactory.getEME2000(), date, zero.add(mu));
510
511 FieldKeplerianOrbit<T> kep = new FieldKeplerianOrbit<>(equi);
512
513 Assertions.assertEquals(42166712.000, equi.getA().getReal(), Utils.epsilonTest * kep.getA().getReal());
514 Assertions.assertEquals(0.110283316961361e-03, kep.getE().getReal(), Utils.epsilonE
515 * FastMath.abs(kep.getE().getReal()));
516 Assertions.assertEquals(0.166901168553917e-03, kep.getI().getReal(), Utils.epsilonAngle
517 * FastMath.abs(kep.getI().getReal()));
518 Assertions.assertEquals(MathUtils.normalizeAngle(-3.87224326008837, kep.getPerigeeArgument().getReal()),
519 kep.getPerigeeArgument().getReal(), Utils.epsilonTest
520 * FastMath.abs(kep.getPerigeeArgument().getReal()));
521 Assertions.assertEquals(MathUtils.normalizeAngle(5.51473467358854, kep
522 .getRightAscensionOfAscendingNode().getReal()), kep
523 .getRightAscensionOfAscendingNode().getReal(), Utils.epsilonTest
524 * FastMath.abs(kep.getRightAscensionOfAscendingNode().getReal()));
525 Assertions.assertEquals(MathUtils.normalizeAngle(3.65750858649982, kep.getMeanAnomaly().getReal()), kep
526 .getMeanAnomaly().getReal(), Utils.epsilonTest * FastMath.abs(kep.getMeanAnomaly().getReal()));
527
528 }
529
530 private <T extends CalculusFieldElement<T>> void doTestHyperbolic1(Field<T> field) {
531 try {
532 T zero = field.getZero();
533 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
534 new FieldEquinoctialOrbit<>(zero.add(-42166712.0), zero.add(1.9), zero.add(0.5), zero.add(0.01), zero.add(-0.02), zero.add(5.300),
535 PositionAngleType.MEAN, FramesFactory.getEME2000(), date, zero.add(mu));
536 Assertions.fail("an exception should have been thrown");
537 } catch (OrekitIllegalArgumentException oe) {
538 Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
539 }
540 }
541
542 private <T extends CalculusFieldElement<T>> void doTestHyperbolic2(Field<T> field) {
543 T zero = field.getZero();
544 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
545 FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(-42166712.0), zero.add(1.9), zero.add(0.5),
546 zero.add(0.01), zero.add(-0.02), zero.add(5.300),
547 PositionAngleType.MEAN, FramesFactory.getEME2000(), date, zero.add(mu));
548 try {
549 new FieldEquinoctialOrbit<>(orbit.getPVCoordinates(), orbit.getFrame(), orbit.getMu());
550 Assertions.fail("an exception should have been thrown");
551 } catch (OrekitIllegalArgumentException oe) {
552 Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
553 }
554 }
555
556 private <T extends CalculusFieldElement<T>> void doTestToOrbitWithoutDerivatives(Field<T> field) {
557 T zero = field.getZero();
558 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
559
560 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
561 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
562 FieldPVCoordinates<T> pvCoordinates = new FieldPVCoordinates<>(position, velocity);
563 FieldEquinoctialOrbit<T> fieldOrbit = new FieldEquinoctialOrbit<>(pvCoordinates, FramesFactory.getEME2000(), date, zero.add(mu));
564 EquinoctialOrbit orbit = fieldOrbit.toOrbit();
565 Assertions.assertFalse(orbit.hasNonKeplerianAcceleration());
566 MatcherAssert.assertThat(orbit.getA(), relativelyCloseTo(fieldOrbit.getA().getReal(), 0));
567 MatcherAssert.assertThat(orbit.getEquinoctialEx(), relativelyCloseTo(fieldOrbit.getEquinoctialEx().getReal(), 0));
568 MatcherAssert.assertThat(orbit.getEquinoctialEy(), relativelyCloseTo(fieldOrbit.getEquinoctialEy().getReal(), 0));
569 MatcherAssert.assertThat(orbit.getHx(), relativelyCloseTo(fieldOrbit.getHx().getReal(), 0));
570 MatcherAssert.assertThat(orbit.getHy(), relativelyCloseTo(fieldOrbit.getHy().getReal(), 0));
571 MatcherAssert.assertThat(orbit.getLv(), relativelyCloseTo(fieldOrbit.getLv().getReal(), 0));
572 }
573
574 private <T extends CalculusFieldElement<T>> void doTestToOrbitWithDerivatives(Field<T> field) {
575 T zero = field.getZero();
576 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
577
578 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
579 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
580 T r2 = position.getNormSq();
581 T r = r2.sqrt();
582 final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(zero.add(mu).negate()),
583 position);
584 final FieldVector3D<T> nonKeplerianAcceleration = keplerianAcceleration.scalarMultiply(1.1);
585 final FieldPVCoordinates<T> pvCoordinates = new FieldPVCoordinates<>(position, velocity, nonKeplerianAcceleration);
586 final FieldEquinoctialOrbit<T> fieldOrbit = new FieldEquinoctialOrbit<>(pvCoordinates, FramesFactory.getEME2000(), date, zero.add(mu));
587 EquinoctialOrbit orbit = fieldOrbit.toOrbit();
588 Assertions.assertTrue(orbit.hasNonKeplerianAcceleration());
589 MatcherAssert.assertThat(orbit.getA(), relativelyCloseTo(fieldOrbit.getA().getReal(), 0));
590 MatcherAssert.assertThat(orbit.getEquinoctialEx(), relativelyCloseTo(fieldOrbit.getEquinoctialEx().getReal(), 0));
591 MatcherAssert.assertThat(orbit.getEquinoctialEy(), relativelyCloseTo(fieldOrbit.getEquinoctialEy().getReal(), 0));
592 MatcherAssert.assertThat(orbit.getHx(), relativelyCloseTo(fieldOrbit.getHx().getReal(), 0));
593 MatcherAssert.assertThat(orbit.getHy(), relativelyCloseTo(fieldOrbit.getHy().getReal(), 0));
594 MatcherAssert.assertThat(orbit.getLv(), relativelyCloseTo(fieldOrbit.getLv().getReal(), 0));
595 MatcherAssert.assertThat(orbit.getADot(), relativelyCloseTo(fieldOrbit.getADot().getReal(), 0));
596 MatcherAssert.assertThat(orbit.getEquinoctialExDot(), relativelyCloseTo(fieldOrbit.getEquinoctialExDot().getReal(), 0));
597 MatcherAssert.assertThat(orbit.getEquinoctialEyDot(), relativelyCloseTo(fieldOrbit.getEquinoctialEyDot().getReal(), 0));
598 MatcherAssert.assertThat(orbit.getHxDot(), relativelyCloseTo(fieldOrbit.getHxDot().getReal(), 0));
599 MatcherAssert.assertThat(orbit.getHyDot(), relativelyCloseTo(fieldOrbit.getHyDot().getReal(), 0));
600 MatcherAssert.assertThat(orbit.getLvDot(), relativelyCloseTo(fieldOrbit.getLvDot().getReal(), 0));
601 MatcherAssert.assertThat(orbit.getLEDot(), relativelyCloseTo(fieldOrbit.getLEDot().getReal(), 0));
602 MatcherAssert.assertThat(orbit.getLMDot(), relativelyCloseTo(fieldOrbit.getLMDot().getReal(), 0));
603 MatcherAssert.assertThat(orbit.getEDot(), relativelyCloseTo(fieldOrbit.getEDot().getReal(), 0));
604 MatcherAssert.assertThat(orbit.getIDot(), relativelyCloseTo(fieldOrbit.getIDot().getReal(), 0));
605 MatcherAssert.assertThat(orbit.getLDot(PositionAngleType.TRUE), relativelyCloseTo(fieldOrbit.getLDot(PositionAngleType.TRUE).getReal(), 0));
606 MatcherAssert.assertThat(orbit.getLDot(PositionAngleType.ECCENTRIC), relativelyCloseTo(fieldOrbit.getLDot(PositionAngleType.ECCENTRIC).getReal(), 0));
607 MatcherAssert.assertThat(orbit.getLDot(PositionAngleType.MEAN), relativelyCloseTo(fieldOrbit.getLDot(PositionAngleType.MEAN).getReal(), 0));
608 }
609
610 private <T extends CalculusFieldElement<T>> void doTestNumericalIssue25(Field<T> field) {
611 T zero = field.getZero();
612 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3782116.14107698), zero.add(416663.11924914), zero.add(5875541.62103057));
613 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-6349.7848910501), zero.add(288.4061811651), zero.add(4066.9366759691));
614 FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
615 FramesFactory.getEME2000(),
616 new FieldAbsoluteDate<>(field, "2004-01-01T23:00:00.000",
617 TimeScalesFactory.getUTC()),
618 zero.add(3.986004415E14));
619 Assertions.assertEquals(0.0, orbit.getE().getReal(), 2.0e-14);
620 }
621
622 private <T extends CalculusFieldElement<T>> void doTestAnomaly(Field<T> field) {
623
624 T one = field.getOne();
625 T zero = field.getZero();
626 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
627
628 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
629 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
630
631 FieldEquinoctialOrbit<T> p = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), date, zero.add(mu));
632 FieldKeplerianOrbit<T> kep = new FieldKeplerianOrbit<>(p);
633
634 T e = p.getE();
635 T eRatio = one.subtract(e).divide(e.add(1)).sqrt();
636 T paPraan = kep.getPerigeeArgument().add(
637 kep.getRightAscensionOfAscendingNode());
638
639 T lv = zero.add(1.1);
640
641 T lE = eRatio.multiply(lv.subtract(paPraan).divide(2).tan()).atan().multiply(2).add(paPraan);
642 T lM = lE.subtract(e.multiply(lE.subtract(paPraan).sin()));
643
644 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
645 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
646 p.getFrame(), p.getDate(), p.getMu());
647 Assertions.assertEquals(p.getLv().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
648 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
649 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
650 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
651 p.getEquinoctialEy() , p.getHx(), p.getHy() , zero , PositionAngleType.TRUE,
652 p.getFrame(), p.getDate(), p.getMu());
653
654 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
655 p.getEquinoctialEy() , p.getHx(), p.getHy() , lE , PositionAngleType.ECCENTRIC,
656 p.getFrame(), p.getDate(), p.getMu());
657 Assertions.assertEquals(p.getLv().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
658 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
659 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
660 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
661 p.getEquinoctialEy() , p.getHx(), p.getHy() , zero , PositionAngleType.TRUE,
662 p.getFrame(), p.getDate(), p.getMu());
663
664 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
665 p.getEquinoctialEy() , p.getHx(), p.getHy() , lM , PositionAngleType.MEAN,
666 p.getFrame(), p.getDate(), p.getMu());
667 Assertions.assertEquals(p.getLv ().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
668 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
669 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
670
671
672 p = new FieldEquinoctialOrbit<>(p.getA(), zero ,
673 zero, p.getHx(), p.getHy() , p.getLv() , PositionAngleType.TRUE,
674 p.getFrame(), p.getDate(), p.getMu());
675
676 lE = lv;
677 lM = lE;
678
679 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
680 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
681 p.getFrame(), p.getDate(), p.getMu());
682 Assertions.assertEquals(p.getLv().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
683 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
684 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
685 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
686 p.getEquinoctialEy() , p.getHx(), p.getHy() , zero , PositionAngleType.TRUE,
687 p.getFrame(), p.getDate(), p.getMu());
688
689 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
690 p.getEquinoctialEy() , p.getHx(), p.getHy() , lE , PositionAngleType.ECCENTRIC,
691 p.getFrame(), p.getDate(), p.getMu());
692 Assertions.assertEquals(p.getLv().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
693 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
694 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
695 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
696 p.getEquinoctialEy() , p.getHx(), p.getHy() , zero , PositionAngleType.TRUE,
697 p.getFrame(), p.getDate(), p.getMu());
698
699 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
700 p.getEquinoctialEy() , p.getHx(), p.getHy() , lM , PositionAngleType.MEAN, p.getFrame(), p.getDate(), p.getMu());
701 Assertions.assertEquals(p.getLv().getReal(), lv.getReal(), Utils.epsilonAngle * FastMath.abs(lv.getReal()));
702 Assertions.assertEquals(p.getLE().getReal(), lE.getReal(), Utils.epsilonAngle * FastMath.abs(lE.getReal()));
703 Assertions.assertEquals(p.getLM().getReal(), lM.getReal(), Utils.epsilonAngle * FastMath.abs(lM.getReal()));
704 }
705
706 private <T extends CalculusFieldElement<T>> void doTestPositionVelocityNorms(Field<T> field) {
707
708 T zero = field.getZero();
709 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
710
711
712 FieldEquinoctialOrbit<T> p =
713 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.5), zero.add(-0.5), zero.add(1.200), zero.add(2.1),
714 zero.add(0.67), PositionAngleType.TRUE,
715 FramesFactory.getEME2000(), date, zero.add(mu));
716
717 T ex = p.getEquinoctialEx();
718 T ey = p.getEquinoctialEy();
719 T lv = p.getLv();
720 T ksi = ex.multiply(lv.cos()).add(ey.multiply(lv.sin())).add(1);
721 T nu = ex.multiply(lv.sin()).subtract(ey.multiply(lv.cos()));
722 T epsilon = ex.negate().multiply(ex).subtract(ey.multiply(ey)).add(1.0).sqrt();
723
724 T a = p.getA();
725 T na = a.reciprocal().multiply(p.getMu()).sqrt();
726
727 Assertions.assertEquals(a.getReal() * epsilon.getReal() * epsilon.getReal() / ksi.getReal(),
728 p.getPosition().getNorm().getReal(),
729 Utils.epsilonTest * FastMath.abs(p.getPosition().getNorm().getReal()));
730 Assertions.assertEquals(na.getReal() * FastMath.sqrt(ksi.getReal() * ksi.getReal() + nu.getReal() * nu.getReal()) / epsilon.getReal(),
731 p.getPVCoordinates().getVelocity().getNorm().getReal(),
732 Utils.epsilonTest * FastMath.abs(p.getPVCoordinates().getVelocity().getNorm().getReal()));
733
734
735 FieldEquinoctialOrbit<T> pCirEqua =
736 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.1e-8), zero.add(0.1e-8), zero.add(0.1e-8), zero.add(0.1e-8),
737 zero.add(0.67), PositionAngleType.TRUE,
738 FramesFactory.getEME2000(), date, zero.add(mu));
739
740 ex = pCirEqua.getEquinoctialEx();
741 ey = pCirEqua.getEquinoctialEy();
742 lv = pCirEqua.getLv();
743 ksi = ex.multiply(lv.cos()).add(ey.multiply(lv.sin())).add(1);
744 nu = ex.multiply(lv.sin()).subtract(ey.multiply(lv.cos()));
745 epsilon = ex.negate().multiply(ex).subtract(ey.multiply(ey)).add(1.0).sqrt();
746
747 a = pCirEqua.getA();
748 na = a.reciprocal().multiply(p.getMu()).sqrt();
749
750 Assertions.assertEquals(a.getReal() * epsilon.getReal() * epsilon.getReal() / ksi.getReal(), pCirEqua.getPosition()
751 .getNorm().getReal(), Utils.epsilonTest
752 * FastMath.abs(pCirEqua.getPosition().getNorm().getReal()));
753 Assertions.assertEquals(na.getReal() * FastMath.sqrt(ksi.getReal() * ksi.getReal() + nu.getReal() * nu.getReal()) / epsilon.getReal(),
754 pCirEqua.getPVCoordinates().getVelocity().getNorm().getReal(), Utils.epsilonTest
755 * FastMath.abs(pCirEqua.getPVCoordinates().getVelocity().getNorm().getReal()));
756 }
757
758 private <T extends CalculusFieldElement<T>> void doTestGeometry(Field<T> field) {
759
760
761 T one = field.getOne();
762 T zero = field.getZero();
763 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
764
765
766
767 FieldEquinoctialOrbit<T> p =
768 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.5), zero.add(-0.5), zero.add(1.200), zero.add(2.1),
769 zero.add(0.67), PositionAngleType.TRUE,
770 FramesFactory.getEME2000(), date, zero.add(mu));
771
772 FieldVector3D<T> position = p.getPosition();
773 FieldVector3D<T> velocity = p.getPVCoordinates().getVelocity();
774 FieldVector3D<T> momentum = p.getPVCoordinates().getMomentum().normalize();
775
776 T apogeeRadius = p.getA().multiply(p.getE().add(1.0));
777 T perigeeRadius = p.getA().multiply(one.subtract(p.getE()));
778
779 for (T lv = zero; lv.getReal() <= 2 * FastMath.PI; lv = lv.add(2 * FastMath.PI / 100.)) {
780 p = new FieldEquinoctialOrbit<>(p.getA(), p.getEquinoctialEx(),
781 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
782 p.getFrame(), p.getDate(), p.getMu());
783 position = p.getPosition();
784
785
786
787
788 Assertions.assertTrue((position.getNorm().getReal() - apogeeRadius.getReal()) <= (apogeeRadius.getReal() * Utils.epsilonTest));
789 Assertions.assertTrue((position.getNorm().getReal() - perigeeRadius.getReal()) >= (-perigeeRadius.getReal() * Utils.epsilonTest));
790
791 position = position.normalize();
792 velocity = p.getPVCoordinates().getVelocity();
793 velocity = velocity.normalize();
794
795
796
797
798
799 Assertions.assertTrue(FastMath.abs(FieldVector3D.dotProduct(position, momentum).getReal()) < Utils.epsilonTest);
800
801 Assertions.assertTrue(FastMath.abs(FieldVector3D.dotProduct(velocity, momentum).getReal()) < Utils.epsilonTest);
802 }
803
804
805
806 FieldEquinoctialOrbit<T> pCirEqua =
807 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.1e-8), zero.add(0.1e-8), zero.add(0.1e-8), zero.add(0.1e-8),
808 zero.add(0.67), PositionAngleType.TRUE,
809 FramesFactory.getEME2000(), date, zero.add(mu));
810
811 position = pCirEqua.getPosition();
812 velocity = pCirEqua.getPVCoordinates().getVelocity();
813
814 momentum = FieldVector3D.crossProduct(position, velocity).normalize();
815
816 apogeeRadius = pCirEqua.getA().multiply(pCirEqua.getE().add(1));
817 perigeeRadius = pCirEqua.getA().multiply(one.subtract(pCirEqua.getE()));
818
819 Assertions.assertEquals(perigeeRadius.getReal(), apogeeRadius.getReal(), 1.e+4 * Utils.epsilonTest
820 * apogeeRadius.getReal());
821
822 for (T lv = zero; lv.getReal() <= 2 * FastMath.PI; lv =lv.add(2 * FastMath.PI / 100.)) {
823 pCirEqua = new FieldEquinoctialOrbit<>(pCirEqua.getA(), pCirEqua.getEquinoctialEx(),
824 pCirEqua.getEquinoctialEy() , pCirEqua.getHx(), pCirEqua.getHy() , lv , PositionAngleType.TRUE,
825 pCirEqua.getFrame(), p.getDate(), p.getMu());
826 position = pCirEqua.getPosition();
827
828
829
830 Assertions.assertTrue((position.getNorm().getReal() - apogeeRadius.getReal()) <= (apogeeRadius.getReal() * Utils.epsilonTest));
831 Assertions.assertTrue((position.getNorm().getReal() - perigeeRadius.getReal()) >= (-perigeeRadius.getReal() * Utils.epsilonTest));
832
833 position = position.normalize();
834 velocity = pCirEqua.getPVCoordinates().getVelocity();
835 velocity = velocity.normalize();
836
837
838
839
840
841 Assertions.assertTrue(FastMath.abs(FieldVector3D.dotProduct(position, momentum).getReal()) < Utils.epsilonTest);
842
843 Assertions.assertTrue(FastMath.abs(FieldVector3D.dotProduct(velocity, momentum).getReal()) < Utils.epsilonTest);
844 }
845 }
846
847 private <T extends CalculusFieldElement<T>> void doTestRadiusOfCurvature(Field<T> field) {
848 T one = field.getOne();
849 T zero = field.getZero();
850 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
851
852
853 FieldEquinoctialOrbit<T> p =
854 new FieldEquinoctialOrbit<>(zero.add(42166712.0), zero.add(0.5), zero.add(-0.5), zero.add(1.200), zero.add(2.1),
855 zero.add(0.67), PositionAngleType.TRUE,
856 FramesFactory.getEME2000(), date, zero.add(mu));
857
858
859 FieldVector3D<T> u = p.getPVCoordinates().getMomentum().orthogonal();
860 FieldVector3D<T> v = FieldVector3D.crossProduct(p.getPVCoordinates().getMomentum(), u).normalize();
861
862
863 T xDot = FieldVector3D.dotProduct(p.getPVCoordinates().getVelocity(), u);
864 T yDot = FieldVector3D.dotProduct(p.getPVCoordinates().getVelocity(), v);
865 T xDotDot = FieldVector3D.dotProduct(p.getPVCoordinates().getAcceleration(), u);
866 T yDotDot = FieldVector3D.dotProduct(p.getPVCoordinates().getAcceleration(), v);
867 T dot2 = xDot.multiply(xDot).add(yDot.multiply(yDot));
868 T rCart = dot2.multiply(dot2.sqrt()).divide(
869 xDot.multiply(yDotDot).subtract(yDot.multiply(xDotDot).abs()));
870
871
872 T ex = p.getEquinoctialEx();
873 T ey = p.getEquinoctialEy();
874 T f = ex.multiply(p.getLE().cos()).add( ey.multiply(p.getLE().sin()));
875
876 T oMf2 = one.subtract(f.multiply(f));
877 T rOrb = p.getA().multiply(oMf2).multiply(oMf2.divide(
878 one.subtract(ex.multiply(ex).add(ey.multiply(ey)))).sqrt());
879
880
881 Assertions.assertEquals(rCart.getReal(), rOrb.getReal(), 2.0e-15 * p.getA().getReal());
882
883
884
885 Assertions.assertEquals(0.8477 * p.getA().getReal(), rCart.getReal(), 1.0e-4 * p.getA().getReal());
886
887 }
888
889 private <T extends CalculusFieldElement<T>> void doTestSymmetry(Field<T> field) {
890 T zero = field.getZero();
891 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
892
893
894 FieldVector3D<T> position = new FieldVector3D<>(zero.add(4512.9), zero.add(18260.), zero.add( -5127.));
895 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(134664.6), zero.add(90066.8), zero.add(72047.6));
896
897 FieldEquinoctialOrbit<T> p = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
898 FramesFactory.getEME2000(), date, zero.add(mu));
899
900 FieldVector3D<T> positionOffset = p.getPosition().subtract(position);
901 FieldVector3D<T> velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
902
903 Assertions.assertEquals(0, positionOffset.getNorm().getReal(), 7.5e-12);
904 Assertions.assertEquals(0, velocityOffset.getNorm().getReal(), 1.0e-15);
905
906
907 position = new FieldVector3D<>(zero.add(33051.2), zero.add(26184.9), zero.add(-1.3E-5));
908 velocity = new FieldVector3D<>(zero.add(-60376.2), zero.add(76208.), zero.add(2.7E-4));
909
910 p = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
911 FramesFactory.getEME2000(), date, zero.add(mu));
912
913 positionOffset = p.getPosition().subtract(position);
914 velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
915
916 Assertions.assertEquals(0, positionOffset.getNorm().getReal(), 1.1e-11);
917 Assertions.assertEquals(0, velocityOffset.getNorm().getReal(), 1.0e-15);
918 }
919
920 private <T extends CalculusFieldElement<T>> void doTestNonInertialFrame(Field<T> field) throws IllegalArgumentException {
921
922 T zero = field.getZero();
923 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
924
925 FieldVector3D<T> position = new FieldVector3D<>(zero.add(4512.9), zero.add(18260.), zero.add(-5127.));
926 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(134664.6), zero.add(90066.8), zero.add(72047.6));
927 FieldPVCoordinates<T> FieldPVCoordinates = new FieldPVCoordinates<>( position, velocity);
928 new FieldEquinoctialOrbit<>(FieldPVCoordinates,
929 new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false),
930 date, zero.add(mu));
931 }
932
933 private <T extends CalculusFieldElement<T>> void doTestJacobianReference(Field<T> field) {
934
935 T zero = field.getZero();
936 FieldAbsoluteDate<T> dateTca = new FieldAbsoluteDate<>(field, 2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
937 double mu = 3.986004415e+14;
938 FieldEquinoctialOrbit<T> orbEqu = new FieldEquinoctialOrbit<>(zero.add(7000000.0), zero.add(0.01), zero.add(-0.02), zero.add(1.2), zero.add(2.1),
939 zero.add(FastMath.toRadians(40.)), PositionAngleType.MEAN,
940 FramesFactory.getEME2000(), dateTca, zero.add(mu));
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011 FieldVector3D<T> pRef = new FieldVector3D<>(zero.add(2004367.298657628707588), zero.add(6575317.978060320019722), zero.add(-1518024.843913963763043));
1012 FieldVector3D<T> vRef = new FieldVector3D<>(zero.add(5574.048661495634406), zero.add(-368.839015744295409), zero.add(5009.529487849066754));
1013 double[][] jRefR = {
1014 { 0.56305379787310628, 1.8470954710993663, -0.42643364527246025, 1370.4369387322224, -90.682848736736688 , 1231.6441195141242 },
1015 { 9.52434720041122055E-008, 9.49704503778007296E-008, 4.46607520107935678E-008, 1.69704446323098610E-004, 7.05603505855828105E-005, 1.14825140460141970E-004 },
1016 { -5.41784097802642701E-008, 9.54903765833015538E-008, -8.95815777332234450E-008, 1.01864980963344096E-004, -1.03194262242761416E-004, 1.40668700715197768E-004 },
1017 { 1.96680305426455816E-007, -1.12388745957974467E-007, -2.27118924123407353E-007, 2.06472886488132167E-004, -1.17984506564646906E-004, -2.38427023682723818E-004 },
1018 { -2.24382495052235118E-007, 1.28218568601277626E-007, 2.59108357381747656E-007, 1.89034327703662092E-004, -1.08019615830663994E-004, -2.18289640324466583E-004 },
1019 { -3.04001022071876804E-007, 1.22214683774559989E-007, 1.35141804810132761E-007, -1.34034616931480536E-004, -2.14283975204169379E-004, 1.29018773893081404E-004 }
1020 };
1021 T[][] jRef = MathArrays.buildArray(field, 6, 6);
1022 for (int ii = 0; ii<6; ii++){
1023 for (int jj = 0; jj< 6; jj++){
1024 jRef[ii][jj] = zero.add(jRefR[ii][jj]);
1025 }
1026 }
1027 FieldPVCoordinates<T> pv = orbEqu.getPVCoordinates();
1028 Assertions.assertEquals(0, pv.getPosition().subtract(pRef).getNorm().getReal(), 2.0e-16 * pRef.getNorm().getReal());
1029 Assertions.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm().getReal(), 2.0e-16 * vRef.getNorm().getReal());
1030
1031 T[][] jacobian = MathArrays.buildArray(field, 6, 6);
1032 orbEqu.getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);
1033
1034
1035 for (int i = 0; i < jacobian.length; i++) {
1036 T[] row = jacobian[i];
1037 T[] rowRef = jRef[i];
1038 for (int j = 0; j < row.length; j++) {
1039 Assertions.assertEquals(0, (row[j].getReal() - rowRef[j].getReal()) / rowRef[j].getReal(), 4.0e-15);
1040 }
1041 }
1042
1043 }
1044
1045 private <T extends CalculusFieldElement<T>> void doTestJacobianFinitedifferences(Field<T> field) {
1046
1047 T zero = field.getZero();
1048 FieldAbsoluteDate<T> dateTca = new FieldAbsoluteDate<>(field, 2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
1049 double mu = 3.986004415e+14;
1050 FieldEquinoctialOrbit<T> orbEqu = new FieldEquinoctialOrbit<>(zero.add(7000000.0), zero.add(0.01), zero.add(-0.02), zero.add(1.2), zero.add(2.1),
1051 zero.add(FastMath.toRadians(40.)), PositionAngleType.MEAN,
1052 FramesFactory.getEME2000(), dateTca, zero.add(mu));
1053
1054 for (PositionAngleType type : PositionAngleType.values()) {
1055 T hP = zero.add(2.0);
1056 T[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbEqu, hP);
1057 T[][] jacobian = MathArrays.buildArray(field, 6, 6);
1058 orbEqu.getJacobianWrtCartesian(type, jacobian);
1059
1060 for (int i = 0; i < jacobian.length; i++) {
1061 T[] row = jacobian[i];
1062 T[] rowRef = finiteDiffJacobian[i];
1063 for (int j = 0; j < row.length; j++) {
1064 Assertions.assertEquals(0, (row[j].getReal() - rowRef[j].getReal()) / rowRef[j].getReal(), 8.0e-9);
1065 }
1066 }
1067
1068 T[][] invJacobian = MathArrays.buildArray(field, 6, 6);
1069 orbEqu.getJacobianWrtParameters(type, invJacobian);
1070 MatrixUtils.createFieldMatrix(jacobian).
1071 multiply(MatrixUtils.createFieldMatrix(invJacobian)).
1072 walkInRowOrder(new FieldMatrixPreservingVisitor<T>() {
1073 public void start(int rows, int columns,
1074 int startRow, int endRow, int startColumn, int endColumn) {
1075 }
1076
1077 public void visit(int row, int column, T value) {
1078 Assertions.assertEquals(row == column ? 1.0 : 0.0, value.getReal(), 1.4e-9);
1079 }
1080
1081 public T end() {
1082 return null;
1083 }
1084 });
1085
1086 }
1087
1088 }
1089
1090 private <T extends CalculusFieldElement<T>> T[][] finiteDifferencesJacobian(PositionAngleType type, FieldEquinoctialOrbit<T> orbit, T hP)
1091 {
1092 Field<T> field = hP.getField();
1093 T[][] jacobian = MathArrays.buildArray(field, 6, 6);
1094 for (int i = 0; i < 6; ++i) {
1095 fillColumn(type, i, orbit, hP, jacobian);
1096 }
1097 return jacobian;
1098 }
1099
1100 private <T extends CalculusFieldElement<T>> void fillColumn(PositionAngleType type, int i, FieldEquinoctialOrbit<T> orbit, T hP, T[][] jacobian) {
1101 T zero = hP.getField().getZero();
1102
1103
1104 FieldVector3D<T> p = orbit.getPosition();
1105 FieldVector3D<T> v = orbit.getPVCoordinates().getVelocity();
1106 T hV = hP.multiply(orbit.getMu()).divide(v.getNorm().multiply(p.getNormSq()));
1107
1108 T h;
1109 FieldVector3D<T> dP = new FieldVector3D<>(hP.getField().getZero(), hP.getField().getZero(), hP.getField().getZero());
1110 FieldVector3D<T> dV = new FieldVector3D<>(hP.getField().getZero(), hP.getField().getZero(), hP.getField().getZero());
1111 switch (i) {
1112 case 0:
1113 h = hP;
1114 dP = new FieldVector3D<>(hP, zero, zero);
1115 break;
1116 case 1:
1117 h = hP;
1118 dP = new FieldVector3D<>(zero, hP, zero);
1119 break;
1120 case 2:
1121 h = hP;
1122 dP = new FieldVector3D<>(zero, zero, hP);
1123 break;
1124 case 3:
1125 h = hV;
1126 dV = new FieldVector3D<>(hV, zero, zero);
1127 break;
1128 case 4:
1129 h = hV;
1130 dV = new FieldVector3D<>(zero, hV, zero);
1131 break;
1132 default:
1133 h = hV;
1134 dV = new FieldVector3D<>(zero, zero, hV);
1135 break;
1136 }
1137
1138 FieldEquinoctialOrbit<T> oM4h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, -4, dP), new FieldVector3D<>(1, v, -4, dV)),
1139 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1140 FieldEquinoctialOrbit<T> oM3h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, -3, dP), new FieldVector3D<>(1, v, -3, dV)),
1141 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1142 FieldEquinoctialOrbit<T> oM2h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, -2, dP), new FieldVector3D<>(1, v, -2, dV)),
1143 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1144 FieldEquinoctialOrbit<T> oM1h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, -1, dP), new FieldVector3D<>(1, v, -1, dV)),
1145 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1146 FieldEquinoctialOrbit<T> oP1h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, +1, dP), new FieldVector3D<>(1, v, +1, dV)),
1147 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1148 FieldEquinoctialOrbit<T> oP2h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, +2, dP), new FieldVector3D<>(1, v, +2, dV)),
1149 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1150 FieldEquinoctialOrbit<T> oP3h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, +3, dP), new FieldVector3D<>(1, v, +3, dV)),
1151 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1152 FieldEquinoctialOrbit<T> oP4h = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(1, p, +4, dP), new FieldVector3D<>(1, v, +4, dV)),
1153 orbit.getFrame(), orbit.getDate(), orbit.getMu());
1154
1155 jacobian[0][i] =(zero.add( -3).multiply(oP4h.getA() .subtract(oM4h.getA())).add(
1156 zero.add( 32).multiply(oP3h.getA() .subtract(oM3h.getA()))).subtract(
1157 zero.add(168).multiply(oP2h.getA() .subtract(oM2h.getA()))).add(
1158 zero.add(672).multiply(oP1h.getA() .subtract(oM1h.getA())))).divide(h.multiply(840));
1159 jacobian[1][i] =(zero.add( -3).multiply(oP4h.getEquinoctialEx().subtract(oM4h.getEquinoctialEx())).add(
1160 zero.add( 32).multiply(oP3h.getEquinoctialEx().subtract(oM3h.getEquinoctialEx()))).subtract(
1161 zero.add(168).multiply(oP2h.getEquinoctialEx().subtract(oM2h.getEquinoctialEx()))).add(
1162 zero.add(672).multiply(oP1h.getEquinoctialEx().subtract(oM1h.getEquinoctialEx())))).divide(h.multiply(840));
1163 jacobian[2][i] =(zero.add( -3).multiply(oP4h.getEquinoctialEy().subtract(oM4h.getEquinoctialEy())).add(
1164 zero.add( 32).multiply(oP3h.getEquinoctialEy().subtract(oM3h.getEquinoctialEy()))).subtract(
1165 zero.add(168).multiply(oP2h.getEquinoctialEy().subtract(oM2h.getEquinoctialEy()))).add(
1166 zero.add(672).multiply(oP1h.getEquinoctialEy().subtract(oM1h.getEquinoctialEy())))).divide(h.multiply(840));
1167 jacobian[3][i] =(zero.add( -3).multiply(oP4h.getHx() .subtract(oM4h.getHx())).add(
1168 zero.add( 32).multiply(oP3h.getHx() .subtract(oM3h.getHx()))).subtract(
1169 zero.add(168).multiply(oP2h.getHx() .subtract(oM2h.getHx()))).add(
1170 zero.add(672).multiply(oP1h.getHx() .subtract(oM1h.getHx())))).divide(h.multiply(840));
1171 jacobian[4][i] =(zero.add( -3).multiply(oP4h.getHy() .subtract(oM4h.getHy())).add(
1172 zero.add( 32).multiply(oP3h.getHy() .subtract(oM3h.getHy()))).subtract(
1173 zero.add(168).multiply(oP2h.getHy() .subtract(oM2h.getHy()))).add(
1174 zero.add(672).multiply(oP1h.getHy() .subtract(oM1h.getHy())))).divide(h.multiply(840));
1175 jacobian[5][i] =(zero.add( -3).multiply(oP4h.getL(type) .subtract(oM4h.getL(type))).add(
1176 zero.add( 32).multiply(oP3h.getL(type) .subtract(oM3h.getL(type)))).subtract(
1177 zero.add(168).multiply(oP2h.getL(type) .subtract(oM2h.getL(type)))).add(
1178 zero.add(672).multiply(oP1h.getL(type) .subtract(oM1h.getL(type))))).divide(h.multiply(840));
1179
1180 }
1181
1182 private <T extends CalculusFieldElement<T>> void doTestNonKeplerianDerivatives(Field<T> field) {
1183 final T zero = field.getZero();
1184
1185 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1186 final FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(6896874.444705), field.getZero().add(1956581.072644), field.getZero().add(-147476.245054));
1187 final FieldVector3D<T> velocity = new FieldVector3D<>(field.getZero().add(166.816407662), field.getZero().add(-1106.783301861), field.getZero().add(-7372.745712770));
1188 final FieldVector3D <T> acceleration = new FieldVector3D<>(field.getZero().add(-7.466182457944), field.getZero().add(-2.118153357345), field.getZero().add(0.160004048437));
1189 final TimeStampedFieldPVCoordinates<T> pv = new TimeStampedFieldPVCoordinates<>(date, position, velocity, acceleration);
1190 final Frame frame = FramesFactory.getEME2000();
1191 final T mu = zero.add(Constants.EIGEN5C_EARTH_MU);
1192 final FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(pv, frame, mu);
1193
1194 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getA),
1195 orbit.getADot().getReal(),
1196 4.3e-8);
1197 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getEquinoctialEx),
1198 orbit.getEquinoctialExDot().getReal(),
1199 2.1e-15);
1200 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getEquinoctialEy),
1201 orbit.getEquinoctialEyDot().getReal(),
1202 5.3e-16);
1203 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getHx),
1204 orbit.getHxDot().getReal(),
1205 4.4e-15);
1206 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getHy),
1207 orbit.getHyDot().getReal(),
1208 1.5e-15);
1209 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getLv),
1210 orbit.getLvDot().getReal(),
1211 1.2e-15);
1212 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getLE),
1213 orbit.getLEDot().getReal(),
1214 7.7e-16);
1215 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getLM),
1216 orbit.getLMDot().getReal(),
1217 8.8e-16);
1218 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getE),
1219 orbit.getEDot().getReal(),
1220 6.9e-16);
1221 Assertions.assertEquals(differentiate(pv, frame, mu, FieldEquinoctialOrbit::getI),
1222 orbit.getIDot().getReal(),
1223 3.5e-15);
1224 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.TRUE)),
1225 orbit.getLDot(PositionAngleType.TRUE).getReal(),
1226 1.2e-15);
1227 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.ECCENTRIC)),
1228 orbit.getLDot(PositionAngleType.ECCENTRIC).getReal(),
1229 7.7e-16);
1230 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.MEAN)),
1231 orbit.getLDot(PositionAngleType.MEAN).getReal(),
1232 8.8e-16);
1233
1234 }
1235
1236 private <T extends CalculusFieldElement<T>, S extends Function<FieldEquinoctialOrbit<T>, T>>
1237 double differentiate(TimeStampedFieldPVCoordinates<T> pv, Frame frame, T mu, S picker) {
1238 final DSFactory factory = new DSFactory(1, 1);
1239 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
1240 UnivariateDifferentiableFunction diff =
1241 differentiator.differentiate((UnivariateFunction) dt -> picker.apply(new FieldEquinoctialOrbit<>(pv.shiftedBy(dt),
1242 frame, mu)).getReal());
1243 return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
1244 }
1245
1246 private <T extends CalculusFieldElement<T>> void doTestPositionAngleDerivatives(final Field<T> field) {
1247 final T zero = field.getZero();
1248
1249 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1250 final FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(6896874.444705), field.getZero().add(1956581.072644), field.getZero().add(-147476.245054));
1251 final FieldVector3D<T> velocity = new FieldVector3D<>(field.getZero().add(166.816407662), field.getZero().add(-1106.783301861), field.getZero().add(-7372.745712770));
1252 final FieldVector3D <T> acceleration = new FieldVector3D<>(field.getZero().add(-7.466182457944), field.getZero().add(-2.118153357345), field.getZero().add(0.160004048437));
1253 final TimeStampedFieldPVCoordinates<T> pv = new TimeStampedFieldPVCoordinates<>(date, position, velocity, acceleration);
1254 final Frame frame = FramesFactory.getEME2000();
1255 final double mu = Constants.EIGEN5C_EARTH_MU;
1256 final FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(pv, frame, zero.add(mu));
1257
1258 for (PositionAngleType type : PositionAngleType.values()) {
1259 final FieldEquinoctialOrbit<T> rebuilt = new FieldEquinoctialOrbit<>(orbit.getA(),
1260 orbit.getEquinoctialEx(),
1261 orbit.getEquinoctialEy(),
1262 orbit.getHx(),
1263 orbit.getHy(),
1264 orbit.getL(type),
1265 orbit.getADot(),
1266 orbit.getEquinoctialExDot(),
1267 orbit.getEquinoctialEyDot(),
1268 orbit.getHxDot(),
1269 orbit.getHyDot(),
1270 orbit.getLDot(type),
1271 type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
1272 MatcherAssert.assertThat(rebuilt.getA().getReal(), relativelyCloseTo(orbit.getA().getReal(), 1));
1273 MatcherAssert.assertThat(rebuilt.getEquinoctialEx().getReal(), relativelyCloseTo(orbit.getEquinoctialEx().getReal(), 1));
1274 MatcherAssert.assertThat(rebuilt.getEquinoctialEy().getReal(), relativelyCloseTo(orbit.getEquinoctialEy().getReal(), 1));
1275 MatcherAssert.assertThat(rebuilt.getHx().getReal(), relativelyCloseTo(orbit.getHx().getReal(), 1));
1276 MatcherAssert.assertThat(rebuilt.getHy().getReal(), relativelyCloseTo(orbit.getHy().getReal(), 1));
1277 MatcherAssert.assertThat(rebuilt.getADot().getReal(), relativelyCloseTo(orbit.getADot().getReal(), 1));
1278 MatcherAssert.assertThat(rebuilt.getEquinoctialExDot().getReal(), relativelyCloseTo(orbit.getEquinoctialExDot().getReal(), 1));
1279 MatcherAssert.assertThat(rebuilt.getEquinoctialEyDot().getReal(), relativelyCloseTo(orbit.getEquinoctialEyDot().getReal(), 1));
1280 MatcherAssert.assertThat(rebuilt.getHxDot().getReal(), relativelyCloseTo(orbit.getHxDot().getReal(), 1));
1281 MatcherAssert.assertThat(rebuilt.getHyDot().getReal(), relativelyCloseTo(orbit.getHyDot().getReal(), 1));
1282 for (PositionAngleType type2 : PositionAngleType.values()) {
1283 MatcherAssert.assertThat(rebuilt.getL(type2).getReal(), relativelyCloseTo(orbit.getL(type2).getReal(), 1));
1284 MatcherAssert.assertThat(rebuilt.getLDot(type2).getReal(), relativelyCloseTo(orbit.getLDot(type2).getReal(), 1));
1285 }
1286 }
1287
1288 }
1289
1290 private <T extends CalculusFieldElement<T>> void doTestEquatorialRetrograde(final Field<T> field) {
1291 final T zero = field.getZero();
1292 FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(10000000.0), field.getZero(), field.getZero());
1293 FieldVector3D<T> velocity = new FieldVector3D<>(field.getZero(), field.getZero().add(-6500.0), field.getZero());
1294 T r2 = position.getNormSq();
1295 T r = r2.sqrt();
1296 FieldVector3D<T> acceleration = new FieldVector3D<>(r.multiply(r2.reciprocal().multiply(zero.add(mu).negate())), position,
1297 field.getOne(), new FieldVector3D<>(field.getZero().add(-0.1),
1298 field.getZero().add(0.2),
1299 field.getZero().add(0.3)));
1300 FieldPVCoordinates<T> pvCoordinates = new FieldPVCoordinates<>(position, velocity, acceleration);
1301
1302
1303
1304 FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldKeplerianOrbit<>(pvCoordinates,
1305 FramesFactory.getEME2000(),
1306 FieldAbsoluteDate.getJ2000Epoch(field),
1307 zero.add(mu)));
1308 Assertions.assertEquals(10637829.465, orbit.getA().getReal(), 1.0e-3);
1309 Assertions.assertEquals(-738.145, orbit.getADot().getReal(), 1.0e-3);
1310 Assertions.assertEquals(0.05995861, orbit.getE().getReal(), 1.0e-8);
1311 Assertions.assertEquals(-6.523e-5, orbit.getEDot().getReal(), 1.0e-8);
1312 Assertions.assertTrue(Double.isNaN(orbit.getI().getReal()));
1313 Assertions.assertTrue(Double.isNaN(orbit.getHx().getReal()));
1314 Assertions.assertTrue(Double.isNaN(orbit.getHxDot().getReal()));
1315 Assertions.assertTrue(Double.isNaN(orbit.getHy().getReal()));
1316 Assertions.assertTrue(Double.isNaN(orbit.getHyDot().getReal()));
1317
1318 }
1319
1320 private <T extends CalculusFieldElement<T>> void doTestDerivativesConversionSymmetry(Field<T> field) {
1321 T zero = field.getZero();
1322 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-05-01T00:01:20.000", TimeScalesFactory.getUTC());
1323 FieldVector3D<T> position = new FieldVector3D<>(zero.add(6893443.400234382),
1324 zero.add(1886406.1073757345),
1325 zero.add(-589265.1150359757));
1326 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-281.1261461082365),
1327 zero.add(-1231.6165642450928),
1328 zero.add(-7348.756363469432));
1329 FieldVector3D<T> acceleration = new FieldVector3D<>(zero.add(-7.460341170581685),
1330 zero.add(-2.0415957334584527),
1331 zero.add(0.6393322823627762));
1332 FieldPVCoordinates<T> pvCoordinates = new FieldPVCoordinates<>( position, velocity, acceleration);
1333 FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(pvCoordinates, FramesFactory.getEME2000(),
1334 date, zero.add(Constants.EIGEN5C_EARTH_MU));
1335 Assertions.assertTrue(orbit.hasNonKeplerianAcceleration());
1336 T r2 = position.getNormSq();
1337 T r = r2.sqrt();
1338 FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(orbit.getMu().negate()),
1339 position);
1340 Assertions.assertEquals(0.0101, FieldVector3D.distance(keplerianAcceleration, acceleration).getReal(), 1.0e-4);
1341
1342 for (OrbitType type : OrbitType.values()) {
1343 FieldOrbit<T> converted = type.convertType(orbit);
1344 Assertions.assertTrue(converted.hasNonKeplerianAcceleration());
1345 FieldEquinoctialOrbit<T> rebuilt = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.convertType(converted);
1346 Assertions.assertTrue(rebuilt.hasNonKeplerianAcceleration());
1347 Assertions.assertEquals(orbit.getADot().getReal(), rebuilt.getADot().getReal(), 3.0e-13);
1348 Assertions.assertEquals(orbit.getEquinoctialExDot().getReal(), rebuilt.getEquinoctialExDot().getReal(), 1.0e-15);
1349 Assertions.assertEquals(orbit.getEquinoctialEyDot().getReal(), rebuilt.getEquinoctialEyDot().getReal(), 1.0e-15);
1350 Assertions.assertEquals(orbit.getHxDot().getReal(), rebuilt.getHxDot().getReal(), 1.0e-15);
1351 Assertions.assertEquals(orbit.getHyDot().getReal(), rebuilt.getHyDot().getReal(), 1.0e-15);
1352 Assertions.assertEquals(orbit.getLvDot().getReal(), rebuilt.getLvDot().getReal(), 1.0e-15);
1353 }
1354
1355 }
1356
1357 private <T extends CalculusFieldElement<T>> void doTestToString(Field<T> field) {
1358 FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(-29536113.0),
1359 field.getZero().add(30329259.0),
1360 field.getZero().add(-100125.0));
1361 FieldVector3D<T> velocity = new FieldVector3D<>(field.getZero().add(-2194.0),
1362 field.getZero().add(-2141.0),
1363 field.getZero().add(-8.0));
1364 FieldPVCoordinates<T> pvCoordinates = new FieldPVCoordinates<>(position, velocity);
1365 FieldEquinoctialOrbit<T> orbit = new FieldEquinoctialOrbit<>(pvCoordinates, FramesFactory.getEME2000(),
1366 FieldAbsoluteDate.getJ2000Epoch(field), field.getZero().add(mu));
1367 Assertions.assertEquals("equinoctial parameters: {a: 4.225517000282565E7; ex: 5.927324978565528E-4; ey: -0.002062743969643666; hx: 6.401103130239252E-5; hy: -0.0017606836670756732; lv: 134.24111947709974;}",
1368 orbit.toString());
1369 }
1370
1371 private <T extends CalculusFieldElement<T>> void doTestCopyNonKeplerianAcceleration(Field<T> field) {
1372
1373 final Frame eme2000 = FramesFactory.getEME2000();
1374
1375
1376 final FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(42164140),
1377 field.getZero(),
1378 field.getZero());
1379
1380 final FieldPVCoordinates<T> pv =
1381 new FieldPVCoordinates<>(position,
1382 new FieldVector3D<>(field.getZero(),
1383 position.getNorm().reciprocal().multiply(mu).sqrt(),
1384 field.getZero()));
1385
1386 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), field.getZero().add(mu));
1387
1388
1389 final FieldOrbit<T> orbitCopy = new FieldKeplerianOrbit<>(orbit);
1390
1391
1392 final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10);
1393 final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10);
1394
1395 Assertions.assertEquals(0.0,
1396 FieldVector3D.distance(shiftedOrbit.getPosition(),
1397 shiftedOrbitCopy.getPosition()).getReal(),
1398 1.0e-10);
1399 Assertions.assertEquals(0.0,
1400 FieldVector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
1401 shiftedOrbitCopy.getPVCoordinates().getVelocity()).getReal(),
1402 1.0e-10);
1403
1404 }
1405
1406 private <T extends CalculusFieldElement<T>> void doTestNormalize(Field<T> field) {
1407 final T zero = field.getZero();
1408 FieldEquinoctialOrbit<T> withoutDerivatives =
1409 new FieldEquinoctialOrbit<>(zero.newInstance(42166712.0), zero.newInstance(0.005),
1410 zero.newInstance(-0.025), zero.newInstance(0.17),
1411 zero.newInstance(0.34), zero.newInstance(0.4), PositionAngleType.MEAN,
1412 FramesFactory.getEME2000(), FieldAbsoluteDate.getJ2000Epoch(field),
1413 zero.newInstance(mu));
1414 FieldEquinoctialOrbit<T> ref =
1415 new FieldEquinoctialOrbit<>(zero.newInstance(24000000.0), zero.newInstance(-0.012),
1416 zero.newInstance(0.01), zero.newInstance(0.2),
1417 zero.newInstance(0.1), zero.newInstance(-6.28), PositionAngleType.MEAN,
1418 FramesFactory.getEME2000(), FieldAbsoluteDate.getJ2000Epoch(field),
1419 zero.newInstance(mu));
1420
1421 FieldEquinoctialOrbit<T> normalized1 = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.normalize(withoutDerivatives, ref);
1422 Assertions.assertFalse(normalized1.hasNonKeplerianAcceleration());
1423 Assertions.assertEquals(0.0, normalized1.getA().subtract(withoutDerivatives.getA()).getReal(), 1.0e-6);
1424 Assertions.assertEquals(0.0, normalized1.getEquinoctialEx().subtract(withoutDerivatives.getEquinoctialEx()).getReal(), 1.0e-10);
1425 Assertions.assertEquals(0.0, normalized1.getEquinoctialEy().subtract(withoutDerivatives.getEquinoctialEy()).getReal(), 1.0e-10);
1426 Assertions.assertEquals(0.0, normalized1.getHx().subtract(withoutDerivatives.getHx()).getReal(), 1.0e-10);
1427 Assertions.assertEquals(0.0, normalized1.getHy().subtract(withoutDerivatives.getHy()).getReal(), 1.0e-10);
1428 Assertions.assertEquals(-MathUtils.TWO_PI, normalized1.getLv().subtract(withoutDerivatives.getLv()).getReal(), 1.0e-10);
1429
1430 T[] p = MathArrays.buildArray(field, 6);
1431 T[] pDot = MathArrays.buildArray(field, 6);
1432 for (int i = 0; i < pDot.length; ++i) {
1433 pDot[i] = zero.newInstance(i);
1434 }
1435 OrbitType.EQUINOCTIAL.mapOrbitToArray(withoutDerivatives, PositionAngleType.TRUE, p, null);
1436 FieldEquinoctialOrbit<T> withDerivatives = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.mapArrayToOrbit(p, pDot,
1437 PositionAngleType.TRUE,
1438 withoutDerivatives.getDate(),
1439 withoutDerivatives.getMu(),
1440 withoutDerivatives.getFrame());
1441 FieldEquinoctialOrbit<T> normalized2 = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.normalize(withDerivatives, ref);
1442 Assertions.assertTrue(normalized2.hasNonKeplerianAcceleration());
1443 Assertions.assertEquals(0.0, normalized2.getA().subtract(withDerivatives.getA()).getReal(), 1.0e-6);
1444 Assertions.assertEquals(0.0, normalized2.getEquinoctialEx().subtract(withDerivatives.getEquinoctialEx()).getReal(), 1.0e-10);
1445 Assertions.assertEquals(0.0, normalized2.getEquinoctialEy().subtract(withDerivatives.getEquinoctialEy()).getReal(), 1.0e-10);
1446 Assertions.assertEquals(0.0, normalized2.getHx().subtract(withDerivatives.getHx()).getReal(), 1.0e-10);
1447 Assertions.assertEquals(0.0, normalized2.getHy().subtract(withDerivatives.getHy()).getReal(), 1.0e-10);
1448 Assertions.assertEquals(-MathUtils.TWO_PI, normalized2.getLv().subtract(withDerivatives.getLv()).getReal(), 1.0e-10);
1449 Assertions.assertEquals(0.0, normalized2.getADot().subtract(withDerivatives.getADot()).getReal(), 1.0e-6);
1450 Assertions.assertEquals(0.0, normalized2.getEquinoctialExDot().subtract(withDerivatives.getEquinoctialExDot()).getReal(), 1.0e-10);
1451 Assertions.assertEquals(0.0, normalized2.getEquinoctialEyDot().subtract(withDerivatives.getEquinoctialEyDot()).getReal(), 1.0e-10);
1452 Assertions.assertEquals(0.0, normalized2.getHxDot().subtract(withDerivatives.getHxDot()).getReal(), 1.0e-10);
1453 Assertions.assertEquals(0.0, normalized2.getHyDot().subtract(withDerivatives.getHyDot()).getReal(), 1.0e-10);
1454 Assertions.assertEquals(0.0, normalized2.getLvDot().subtract(withDerivatives.getLvDot()).getReal(), 1.0e-10);
1455
1456 }
1457
1458 }