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