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