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