1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.analytical;
18
19
20 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.Collection;
23 import java.util.List;
24
25 import org.hipparchus.Field;
26 import org.hipparchus.CalculusFieldElement;
27 import org.hipparchus.exception.DummyLocalizable;
28 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
30 import org.hipparchus.util.Decimal64Field;
31 import org.hipparchus.util.FastMath;
32 import org.hipparchus.util.MathUtils;
33 import org.junit.After;
34 import org.junit.Assert;
35 import org.junit.Before;
36 import org.junit.Test;
37 import org.orekit.Utils;
38 import org.orekit.attitudes.Attitude;
39 import org.orekit.attitudes.AttitudeProvider;
40 import org.orekit.attitudes.FieldAttitude;
41 import org.orekit.attitudes.LofOffset;
42 import org.orekit.bodies.GeodeticPoint;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.errors.OrekitException;
45 import org.orekit.errors.OrekitMessages;
46 import org.orekit.forces.gravity.potential.GravityFieldFactory;
47 import org.orekit.forces.gravity.potential.TideSystem;
48 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
49 import org.orekit.frames.Frame;
50 import org.orekit.frames.FramesFactory;
51 import org.orekit.frames.LOFType;
52 import org.orekit.frames.TopocentricFrame;
53 import org.orekit.orbits.FieldCircularOrbit;
54 import org.orekit.orbits.FieldEquinoctialOrbit;
55 import org.orekit.orbits.FieldKeplerianOrbit;
56 import org.orekit.orbits.FieldOrbit;
57 import org.orekit.orbits.OrbitType;
58 import org.orekit.orbits.PositionAngle;
59 import org.orekit.propagation.FieldSpacecraftState;
60 import org.orekit.propagation.PropagationType;
61 import org.orekit.propagation.Propagator;
62 import org.orekit.propagation.events.FieldApsideDetector;
63 import org.orekit.propagation.events.FieldDateDetector;
64 import org.orekit.propagation.events.FieldElevationDetector;
65 import org.orekit.propagation.events.FieldEventDetector;
66 import org.orekit.propagation.events.FieldNodeDetector;
67 import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
68 import org.orekit.propagation.sampling.FieldOrekitFixedStepHandler;
69 import org.orekit.propagation.semianalytical.dsst.FieldDSSTPropagator;
70 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
71 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
72 import org.orekit.time.AbsoluteDate;
73 import org.orekit.time.DateComponents;
74 import org.orekit.time.FieldAbsoluteDate;
75 import org.orekit.time.TimeComponents;
76 import org.orekit.time.TimeScalesFactory;
77 import org.orekit.utils.CartesianDerivativesFilter;
78 import org.orekit.utils.FieldPVCoordinates;
79 import org.orekit.utils.FieldPVCoordinatesProvider;
80 import org.orekit.utils.IERSConventions;
81 import org.orekit.utils.PVCoordinatesProvider;
82 import org.orekit.utils.TimeStampedFieldPVCoordinates;
83
84
85 public class FieldEcksteinHechlerPropagatorTest {
86
87 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
88
89 private double mu;
90 private double ae;
91
92 @Before
93 public void setUp() {
94 Utils.setDataRoot("regular-data");
95 mu = 3.9860047e14;
96 ae = 6.378137e6;
97 double[][] cnm = new double[][] {
98 { 0 }, { 0 }, { -1.08263e-3 }, { 2.54e-6 }, { 1.62e-6 }, { 2.3e-7 }, { -5.5e-7 }
99 };
100 double[][] snm = new double[][] {
101 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
102 };
103
104 provider = GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN, cnm, snm);
105 }
106
107 @Test
108 public void sameDateCartesian() {
109 doSameDateCartesian(Decimal64Field.getInstance());
110 }
111
112 private <T extends CalculusFieldElement<T>> void doSameDateCartesian(Field<T> field) {
113 T zero = field.getZero();
114 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
115
116
117
118 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
119 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
120
121 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
122 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
123 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
124
125
126
127 FieldEcksteinHechlerPropagator<T> extrapolator =
128 new FieldEcksteinHechlerPropagator<>(initialOrbit, provider);
129
130
131
132 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(initDate);
133
134
135 Assert.assertEquals(0.0,
136 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
137 finalOrbit.getPVCoordinates().getPosition()).getReal(),
138 1.0e-8);
139
140
141
142
143
144
145
146
147
148
149
150 Assert.assertEquals(0.137,
151 FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
152 finalOrbit.getPVCoordinates().getVelocity()).getReal(),
153 1.0e-3);
154 Assert.assertEquals(125.2, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.1);
155
156 }
157
158 @Test
159 public void sameDateKeplerian() {
160 doSameDateKeplerian(Decimal64Field.getInstance());
161 }
162
163 private <T extends CalculusFieldElement<T>> void doSameDateKeplerian(Field<T> field) {
164
165 T zero = field.getZero();
166 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
167
168
169
170 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
171 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(7209668.0), zero.add(0.5e-4), zero.add(1.7), zero.add( 2.1), zero.add( 2.9),
172 zero.add(6.2), PositionAngle.TRUE,
173 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
174
175
176
177 FieldEcksteinHechlerPropagator<T> extrapolator =
178 new FieldEcksteinHechlerPropagator<>(initialOrbit, zero.add(Propagator.DEFAULT_MASS), provider);
179
180
181
182 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(initDate);
183
184
185 Assert.assertEquals(0.0,
186 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
187 finalOrbit.getPVCoordinates().getPosition()).getReal(),
188 3.0e-8);
189
190
191
192
193
194
195
196
197
198
199 Assert.assertEquals(0.137,
200 FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
201 finalOrbit.getPVCoordinates().getVelocity()).getReal(),
202 1.0e-3);
203 Assert.assertEquals(126.8, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.1);
204
205 }
206
207 @Test
208 public void almostSphericalBody() {
209 doAlmostSphericalBody(Decimal64Field.getInstance());
210 }
211
212 private <T extends CalculusFieldElement<T>> void doAlmostSphericalBody(Field<T> field) {
213 T zero = field.getZero();
214 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
215
216
217
218
219 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
220 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
221
222 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
223 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
224 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
225
226
227
228
229
230
231
232 UnnormalizedSphericalHarmonicsProvider kepProvider =
233 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
234 TideSystem.UNKNOWN,
235 new double[][] {
236 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
237 }, new double[][] {
238 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
239 });
240
241
242
243 FieldEcksteinHechlerPropagator<T> extrapolatorAna =
244 new FieldEcksteinHechlerPropagator<>(initialOrbit, kepProvider);
245 FieldKeplerianPropagator<T> extrapolatorKep = new FieldKeplerianPropagator<>(initialOrbit);
246
247
248
249 double delta_t = 100.0;
250 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
251
252 FieldSpacecraftState<T> finalOrbitAna = extrapolatorAna.propagate(extrapDate);
253 FieldSpacecraftState<T> finalOrbitKep = extrapolatorKep.propagate(extrapDate);
254
255 Assert.assertEquals(finalOrbitAna.getDate().durationFrom(extrapDate).getReal(), 0.0,
256 Utils.epsilonTest);
257
258 Assert.assertEquals(finalOrbitAna.getA().getReal(), finalOrbitKep.getA().getReal(), 10
259 * Utils.epsilonTest * finalOrbitKep.getA().getReal());
260 Assert.assertEquals(finalOrbitAna.getEquinoctialEx().getReal(), finalOrbitKep.getEquinoctialEx().getReal(), Utils.epsilonE
261 * finalOrbitKep.getE().getReal());
262 Assert.assertEquals(finalOrbitAna.getEquinoctialEy().getReal(), finalOrbitKep.getEquinoctialEy().getReal(), Utils.epsilonE
263 * finalOrbitKep.getE().getReal());
264 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx().getReal(), finalOrbitKep.getHx().getReal()),
265 finalOrbitKep.getHx().getReal(), Utils.epsilonAngle
266 * FastMath.abs(finalOrbitKep.getI().getReal()));
267 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy().getReal(), finalOrbitKep.getHy().getReal()),
268 finalOrbitKep.getHy().getReal(), Utils.epsilonAngle
269 * FastMath.abs(finalOrbitKep.getI().getReal()));
270 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv().getReal(), finalOrbitKep.getLv().getReal()),
271 finalOrbitKep.getLv().getReal(), Utils.epsilonAngle
272 * FastMath.abs(finalOrbitKep.getLv().getReal()));
273 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE().getReal(), finalOrbitKep.getLE().getReal()),
274 finalOrbitKep.getLE().getReal(), Utils.epsilonAngle
275 * FastMath.abs(finalOrbitKep.getLE().getReal()));
276 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM().getReal(), finalOrbitKep.getLM().getReal()),
277 finalOrbitKep.getLM().getReal(), Utils.epsilonAngle
278 * FastMath.abs(finalOrbitKep.getLM().getReal()));
279
280 }
281
282 @Test
283 public void propagatedCartesian() {
284 doPropagatedCartesian(Decimal64Field.getInstance());
285 }
286
287 private <T extends CalculusFieldElement<T>> void doPropagatedCartesian(Field<T> field) {
288 T zero = field.getZero();
289 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
290
291
292
293 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
294 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
295
296 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
297 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
298 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
299
300
301
302 FieldEcksteinHechlerPropagator<T> extrapolator =
303 new FieldEcksteinHechlerPropagator<>(initialOrbit,
304 new LofOffset(initialOrbit.getFrame(),
305 LOFType.VNC, RotationOrder.XYZ, 0, 0, 0),
306 provider);
307
308
309
310 double delta_t = 100000.0;
311 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
312
313 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(extrapDate);
314
315 Assert.assertEquals(0.0, finalOrbit.getDate().durationFrom(extrapDate).getReal(), 1.0e-9);
316
317
318 T LM = finalOrbit.getLE().subtract(finalOrbit.getEquinoctialEx().multiply(
319 finalOrbit.getLE().sin())).add(finalOrbit.getEquinoctialEy()
320 .multiply(finalOrbit.getLE().cos()));
321
322 Assert.assertEquals(LM.getReal(), finalOrbit.getLM().getReal(), Utils.epsilonAngle
323 * FastMath.abs(finalOrbit.getLM().getReal()));
324
325
326 Assert.assertEquals(FastMath.tan((finalOrbit.getLE().getReal() - finalOrbit.getLv().getReal()) / 2.),
327 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit
328 .getEquinoctialEy()).getReal(), Utils.epsilonAngle);
329
330
331 T deltaM = finalOrbit.getLM().subtract(initialOrbit.getLM());
332 T deltaE = finalOrbit.getLE().subtract(initialOrbit.getLE());
333 T delta = finalOrbit.getEquinoctialEx().multiply(finalOrbit.getLE().sin()).subtract(
334 initialOrbit.getEquinoctialEx().multiply(initialOrbit.getLE().sin())).subtract(
335 finalOrbit.getEquinoctialEy().multiply(finalOrbit.getLE().cos())).add(
336 initialOrbit.getEquinoctialEy().multiply(initialOrbit.getLE().cos()));
337
338 Assert.assertEquals(deltaM.getReal(), deltaE.getReal() - delta.getReal(), Utils.epsilonAngle
339 * FastMath.abs(deltaE.getReal() - delta.getReal()));
340
341
342 T ex = finalOrbit.getEquinoctialEx();
343 T ey = finalOrbit.getEquinoctialEy();
344 T hx = finalOrbit.getHx();
345 T hy = finalOrbit.getHy();
346 T LE = finalOrbit.getLE();
347
348 T ex2 = ex.multiply(ex);
349 T ey2 = ey.multiply(ey);
350 T hx2 = hx.multiply(hx);
351 T hy2 = hy.multiply(hy);
352 T h2p1 = hx2.add(1.).add(hy2);
353 T beta = ex2.negate().add(1.).subtract(ey2).sqrt().add(1.).reciprocal();
354
355 T x3 = ex.negate().add(ey2.multiply(beta).negate().add(1.).multiply(LE.cos())).add(beta.multiply(ex).multiply(ey).multiply(
356 LE.sin()));
357 T y3 = ey.negate().add(ex2.negate().multiply(beta).add(1).multiply(LE.sin())).add(beta.multiply(ex).multiply(ey).multiply(LE.cos()));
358
359 FieldVector3D<T> U = new FieldVector3D<>(hx2.add(1).subtract(hy2).divide(h2p1),
360 hx.multiply(hy).multiply(2).divide(h2p1),
361 hy.multiply(-2).divide(h2p1));
362
363 FieldVector3D<T> V = new FieldVector3D<>(hx.multiply(2).multiply(hy).divide(h2p1),
364 hy2.add(1).subtract(hx2).divide(h2p1),
365 hx.multiply(2).divide(h2p1));
366
367 FieldVector3D<T> r = new FieldVector3D<>(finalOrbit.getA(), new FieldVector3D<>(x3, U, y3, V));
368
369 Assert.assertEquals(finalOrbit.getPVCoordinates().getPosition().getNorm().getReal(),
370 r.getNorm().getReal(),
371 Utils.epsilonTest * r.getNorm().getReal());
372
373 }
374
375 @Test
376 public void propagatedKeplerian() {
377
378 doPropagatedKeplerian(Decimal64Field.getInstance());
379
380 }
381
382 private <T extends CalculusFieldElement<T>> void doPropagatedKeplerian(Field<T> field) {
383 T zero = field.getZero();
384 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
385
386
387 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
388 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(7209668.0), zero.add(0.5e-4), zero.add(1.7), zero.add(2.1), zero.add(2.9),
389 zero.add(6.2), PositionAngle.TRUE,
390 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
391
392
393
394 FieldEcksteinHechlerPropagator<T> extrapolator =
395 new FieldEcksteinHechlerPropagator<>(initialOrbit,
396 new LofOffset(initialOrbit.getFrame(), LOFType.VNC),
397 zero.add(2000.0), provider);
398
399
400
401 double delta_t = 100000.0;
402 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
403
404 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(extrapDate);
405
406 Assert.assertEquals(0.0, finalOrbit.getDate().durationFrom(extrapDate).getReal(), 1.0e-9);
407
408
409 T LM = finalOrbit.getLE().subtract(finalOrbit.getEquinoctialEx().multiply(
410 finalOrbit.getLE().sin())).add(finalOrbit.getEquinoctialEy().multiply(
411 finalOrbit.getLE().cos()));
412
413 Assert.assertEquals(LM.getReal(), finalOrbit.getLM().getReal(), Utils.epsilonAngle);
414
415
416 Assert.assertEquals(FastMath.tan((finalOrbit.getLE().getReal() - finalOrbit.getLv().getReal()) / 2.),
417 tangLEmLv(finalOrbit.getLv(), finalOrbit.getEquinoctialEx(), finalOrbit
418 .getEquinoctialEy()).getReal(), Utils.epsilonAngle);
419
420
421
422 T deltaM = finalOrbit.getLM().subtract(initialOrbit.getLM());
423 T deltaE = finalOrbit.getLE().subtract(initialOrbit.getLE());
424 T delta = finalOrbit.getEquinoctialEx().multiply(finalOrbit.getLE().sin()).subtract(
425 initialOrbit.getEquinoctialEx().multiply(initialOrbit.getLE().sin())).subtract(
426 finalOrbit.getEquinoctialEy().multiply(finalOrbit.getLE().cos())).add(
427 initialOrbit.getEquinoctialEy().multiply(initialOrbit.getLE().cos()));
428
429 Assert.assertEquals(deltaM.getReal(), deltaE.getReal() - delta.getReal(), Utils.epsilonAngle
430 * FastMath.abs(deltaE.getReal() - delta.getReal()));
431
432
433 T ex = finalOrbit.getEquinoctialEx();
434 T ey = finalOrbit.getEquinoctialEy();
435 T hx = finalOrbit.getHx();
436 T hy = finalOrbit.getHy();
437 T LE = finalOrbit.getLE();
438
439 T ex2 = ex.multiply( ex);
440 T ey2 = ey.multiply( ey);
441 T hx2 = hx.multiply( hx);
442 T hy2 = hy.multiply( hy);
443 T h2p1 = hx2.add(1).add(hy2);
444 T beta = ex2.negate().add(1.).subtract(ey2).sqrt().add(1).reciprocal();
445
446 T x3 = ex.negate().add(beta.negate().multiply(ey2).add(1).multiply(LE.cos())).add(
447 beta.multiply(ex).multiply(ey).multiply(LE.sin()));
448 T y3 = ey.negate().add(beta.negate().multiply(ex2).add(1).multiply(LE.sin())).add(
449 beta.multiply(ex).multiply(ey).multiply(LE.cos()));
450
451 FieldVector3D<T> U = new FieldVector3D<>(hx2.subtract(hy2).add(1.).divide(h2p1),
452 hx.multiply(2).multiply(hy).divide(h2p1),
453 hy.multiply(-2.).divide(h2p1));
454 FieldVector3D<T> V = new FieldVector3D<>(hx.multiply(2.).multiply(hy).divide(h2p1),
455 hy2.subtract(hx2).add(1.).divide(h2p1),
456 hx.multiply(2).divide(h2p1));
457 FieldVector3D<T> r = new FieldVector3D<>(finalOrbit.getA(), new FieldVector3D<>(x3, U, y3, V));
458
459 Assert.assertEquals(finalOrbit.getPVCoordinates().getPosition().getNorm().getReal(), r.getNorm().getReal(),
460 Utils.epsilonTest * r.getNorm().getReal());
461
462 }
463
464 @Test
465 public void undergroundOrbit() {
466 doUndergroundOrbit(Decimal64Field.getInstance());
467 }
468
469 private <T extends CalculusFieldElement<T>> void doUndergroundOrbit(Field<T> field) {
470 T zero = field.getZero();
471 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
472
473 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add( 1.0e6), zero.add(4.0e6));
474 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(800.0), zero.add(100.0));
475 FieldAbsoluteDate<T> initDate = date;
476 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
477 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
478 try {
479
480
481 FieldEcksteinHechlerPropagator<T> extrapolator =
482 new FieldEcksteinHechlerPropagator<>(initialOrbit, provider);
483
484
485
486 double delta_t = 0.0;
487 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
488 extrapolator.propagate(extrapDate);
489 Assert.fail("an exception should have been thrown");
490 } catch (OrekitException oe) {
491 Assert.assertEquals(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, oe.getSpecifier());
492 }
493 }
494
495 @Test
496 public void equatorialOrbit() {
497 doEquatorialOrbit(Decimal64Field.getInstance());
498 }
499
500 private <T extends CalculusFieldElement<T>> void doEquatorialOrbit(Field<T> field) {
501 T zero = field.getZero();
502 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
503
504 FieldAbsoluteDate<T> initDate = date;
505 FieldOrbit<T> initialOrbit = new FieldCircularOrbit<>(zero.add(7000000), zero.add(1.0e-4), zero.add(-1.5e-4),
506 zero, zero.add(1.2), zero.add(2.3), PositionAngle.MEAN,
507 FramesFactory.getEME2000(),
508 initDate, zero.add(provider.getMu()));
509 try {
510
511
512 FieldEcksteinHechlerPropagator<T> extrapolator =
513 new FieldEcksteinHechlerPropagator<>(initialOrbit, provider);
514
515
516
517 double delta_t = 0.0;
518 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
519 extrapolator.propagate(extrapDate);
520 Assert.fail("an exception should have been thrown");
521 } catch (OrekitException oe) {
522 Assert.assertEquals(OrekitMessages.ALMOST_EQUATORIAL_ORBIT, oe.getSpecifier());
523 }
524 }
525
526 @Test
527 public void criticalInclination() {
528 doCriticalInclination(Decimal64Field.getInstance());
529 }
530
531 private <T extends CalculusFieldElement<T>> void doCriticalInclination(Field<T> field) {
532 T zero = field.getZero();
533 FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field);
534 FieldOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(zero.add(-3862363.8474653554),
535 zero.add(-3521533.9758022362),
536 zero.add(4647637.852558916)),
537 new FieldVector3D<>(zero.add(65.36170817232278),
538 zero.add(-6056.563439401233),
539 zero.add(-4511.1247889782757))),
540 FramesFactory.getEME2000(),
541 initDate, zero.add(provider.getMu()));
542
543 try {
544
545
546 FieldEcksteinHechlerPropagator<T> extrapolator =
547 new FieldEcksteinHechlerPropagator<>(initialOrbit, provider);
548
549
550
551 double delta_t = 0.0;
552 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
553 extrapolator.propagate(extrapDate);
554 Assert.fail("an exception should have been thrown");
555 } catch (OrekitException oe) {
556 Assert.assertEquals(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT, oe.getSpecifier());
557 }
558 }
559
560 @Test
561 public void tooEllipticalOrbit() {
562 doTooEllipticalOrbit(Decimal64Field.getInstance());
563 }
564
565 private <T extends CalculusFieldElement<T>> void doTooEllipticalOrbit(Field<T> field) {
566 T zero = field.getZero();
567 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
568
569 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add( 4.0e6));
570 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add( 8000.0), zero.add(1000.0));
571 FieldAbsoluteDate<T> initDate = date;
572 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
573 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
574 try {
575
576
577 FieldEcksteinHechlerPropagator<T> extrapolator =
578 new FieldEcksteinHechlerPropagator<>(initialOrbit, provider);
579
580
581
582 double delta_t = 0.0;
583 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
584 extrapolator.propagate(extrapDate);
585 Assert.fail("an exception should have been thrown");
586 } catch (OrekitException oe) {
587 Assert.assertEquals(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL, oe.getSpecifier());
588 }
589 }
590
591 @Test
592 public void hyperbolic() {
593 doHyperbolic(Decimal64Field.getInstance());
594 }
595
596 private <T extends CalculusFieldElement<T>> void doHyperbolic(Field<T> field) {
597 T zero = field.getZero();
598 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
599 FieldKeplerianOrbit<T> hyperbolic =
600 new FieldKeplerianOrbit<>(zero.add(-1.0e10), zero.add(2), zero, zero, zero, zero, PositionAngle.TRUE,
601 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
602 try {
603 FieldEcksteinHechlerPropagator<T> propagator =
604 new FieldEcksteinHechlerPropagator<>(hyperbolic, provider);
605 propagator.propagate(date.shiftedBy(10.0));
606 Assert.fail("an exception should have been thrown");
607 } catch (OrekitException oe) {
608 Assert.assertEquals(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, oe.getSpecifier());
609 }
610 }
611
612 @Test
613 public void wrongAttitude() {
614 doWrongAttitude(Decimal64Field.getInstance());
615 }
616
617 private <T extends CalculusFieldElement<T>> void doWrongAttitude(Field<T> field) {
618 T zero = field.getZero();
619 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
620 FieldKeplerianOrbit<T> orbit =
621 new FieldKeplerianOrbit<>(zero.add(1.0e10), zero.add(1.0e-4), zero.add(1.0e-2), zero, zero, zero, PositionAngle.TRUE,
622 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
623 final DummyLocalizable gasp = new DummyLocalizable("gasp");
624 AttitudeProvider wrongLaw = new AttitudeProvider() {
625
626 @Override
627 public Attitude getAttitude(PVCoordinatesProvider pvProv,
628 AbsoluteDate date, Frame frame)
629 {
630 throw new OrekitException(gasp, new RuntimeException());
631 }
632
633 @Override
634 public <Q extends CalculusFieldElement<Q>> FieldAttitude<Q>
635 getAttitude(FieldPVCoordinatesProvider<Q> pvProv,
636 FieldAbsoluteDate<Q> date, Frame frame)
637 {
638 throw new OrekitException(gasp, new RuntimeException());
639 }
640 };
641 try {
642 FieldEcksteinHechlerPropagator<T> propagator =
643 new FieldEcksteinHechlerPropagator<>(orbit, wrongLaw, provider);
644 propagator.propagate(date.shiftedBy(10.0));
645 Assert.fail("an exception should have been thrown");
646 } catch (OrekitException oe) {
647 Assert.assertSame(gasp, oe.getSpecifier());
648 }
649 }
650
651 @Test
652 public void testAcceleration() {
653 doTestAcceleration(Decimal64Field.getInstance());
654 }
655
656 private <T extends CalculusFieldElement<T>> void doTestAcceleration(Field<T> field) {
657 T zero = field.getZero();
658 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
659 final FieldKeplerianOrbit<T> orbit =
660 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
661 FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
662 FieldEcksteinHechlerPropagator<T> propagator =
663 new FieldEcksteinHechlerPropagator<>(orbit, provider);
664 FieldAbsoluteDate<T> target = date.shiftedBy(10000.0);
665 List<TimeStampedFieldPVCoordinates<T>> sample = new ArrayList<TimeStampedFieldPVCoordinates<T>>();
666 for (double dt : Arrays.asList(-0.5, 0.0, 0.5)) {
667 sample.add(propagator.propagate(target.shiftedBy(dt)).getPVCoordinates());
668 }
669 TimeStampedFieldPVCoordinates<T> interpolated =
670 TimeStampedFieldPVCoordinates.interpolate(target, CartesianDerivativesFilter.USE_P, sample);
671 FieldVector3D<T> computedP = sample.get(1).getPosition();
672 FieldVector3D<T> computedV = sample.get(1).getVelocity();
673 FieldVector3D<T> referenceP = interpolated.getPosition();
674 FieldVector3D<T> referenceV = interpolated.getVelocity();
675 FieldVector3D<T> computedA = sample.get(1).getAcceleration();
676 FieldVector3D<T> referenceA = interpolated.getAcceleration();
677 final FieldCircularOrbit<T> propagated = (FieldCircularOrbit<T>) OrbitType.CIRCULAR.convertType(propagator.propagateOrbit(target, null));
678 final FieldCircularOrbit<T> keplerian =
679 new FieldCircularOrbit<>(propagated.getA(),
680 propagated.getCircularEx(),
681 propagated.getCircularEy(),
682 propagated.getI(),
683 propagated.getRightAscensionOfAscendingNode(),
684 propagated.getAlphaM(), PositionAngle.MEAN,
685 propagated.getFrame(),
686 propagated.getDate(),
687 propagated.getMu());
688 FieldVector3D<T> keplerianP = keplerian.getPVCoordinates().getPosition();
689 FieldVector3D<T> keplerianV = keplerian.getPVCoordinates().getVelocity();
690 FieldVector3D<T> keplerianA = keplerian.getPVCoordinates().getAcceleration();
691
692
693 Assert.assertEquals(0.0, FieldVector3D.distance(referenceP, computedP).getReal(), 1.0e-15);
694 Assert.assertEquals(0.0, FieldVector3D.distance(referenceP, keplerianP).getReal(), 4.0e-9);
695
696
697
698 T computationErrorV = FieldVector3D.distance(referenceV, computedV);
699 T nonKeplerianEffectV = FieldVector3D.distance(referenceV, keplerianV);
700 Assert.assertEquals(0.0, nonKeplerianEffectV.getReal() - computationErrorV.getReal(), 9.0e-12);
701 Assert.assertEquals(2.2e-4, computationErrorV.getReal(), 3.0e-6);
702
703
704
705
706 T computationErrorA = FieldVector3D.distance(referenceA, computedA);
707 T nonKeplerianEffectA = FieldVector3D.distance(referenceA, keplerianA);
708 Assert.assertEquals(1.0e-7, computationErrorA.getReal(), 6.0e-9);
709 Assert.assertEquals(6.37e-3, nonKeplerianEffectA.getReal(), 7.0e-6);
710 Assert.assertTrue(computationErrorA.getReal() < nonKeplerianEffectA.getReal() / 60000);
711
712 }
713
714 @Test
715 public void ascendingNode() {
716 doAscendingNode(Decimal64Field.getInstance());
717 }
718
719 private <T extends CalculusFieldElement<T>> void doAscendingNode(Field<T> field) {
720 T zero = field.getZero();
721 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
722 final FieldKeplerianOrbit<T> orbit =
723 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
724 FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
725 FieldEcksteinHechlerPropagator<T> propagator =
726 new FieldEcksteinHechlerPropagator<>(orbit, provider);
727 FieldNodeDetector<T> detector = new FieldNodeDetector<>(orbit, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
728 Assert.assertTrue(FramesFactory.getITRF(IERSConventions.IERS_2010, true) == detector.getFrame());
729 propagator.addEventDetector(detector);
730
731 FieldAbsoluteDate<T> farTarget = date.shiftedBy(10000.0);
732
733 FieldSpacecraftState<T> propagated = propagator.propagate(farTarget);
734
735 FieldPVCoordinates<T> pv = propagated.getPVCoordinates(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
736 Assert.assertTrue(farTarget.durationFrom(propagated.getDate()).getReal() > 3500.0);
737 Assert.assertTrue(farTarget.durationFrom(propagated.getDate()).getReal() < 4000.0);
738
739 Assert.assertEquals(0, pv.getPosition().getZ().getReal(), 1.0e-6);
740 Assert.assertTrue(pv.getVelocity().getZ().getReal() > 0);
741 Collection<FieldEventDetector<T>> detectors = propagator.getEventsDetectors();
742 Assert.assertEquals(1, detectors.size());
743 propagator.clearEventsDetectors();
744 Assert.assertEquals(0, propagator.getEventsDetectors().size());
745 }
746
747 @Test
748 public void stopAtTargetDate() {
749 doStopAtTargetDate(Decimal64Field.getInstance());
750 }
751
752 private <T extends CalculusFieldElement<T>> void doStopAtTargetDate(Field<T> field) {
753 T zero = field.getZero();
754 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
755 final FieldKeplerianOrbit<T> orbit =
756 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
757 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
758 FieldEcksteinHechlerPropagator<T> propagator =
759 new FieldEcksteinHechlerPropagator<>(orbit, provider);
760 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
761 propagator.addEventDetector(new FieldNodeDetector<>(orbit, itrf).
762 withHandler(new FieldContinueOnEvent<FieldNodeDetector<T>, T>()));
763 FieldAbsoluteDate<T> farTarget = orbit.getDate().shiftedBy(10000.0);
764 FieldSpacecraftState<T> propagated = propagator.propagate(farTarget);
765 Assert.assertEquals(0.0, FastMath.abs(farTarget.durationFrom(propagated.getDate()).getReal()), 1.0e-3);
766 }
767
768 @Test
769 public void perigee() {
770 doPerigee(Decimal64Field.getInstance());
771 }
772
773 private <T extends CalculusFieldElement<T>> void doPerigee(Field<T> field) {
774 T zero = field.getZero();
775 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
776 final FieldKeplerianOrbit<T> orbit =
777 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
778 FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
779 FieldEcksteinHechlerPropagator<T> propagator =
780 new FieldEcksteinHechlerPropagator<>(orbit, provider);
781 propagator.addEventDetector(new FieldApsideDetector<>(orbit));
782 FieldAbsoluteDate<T> farTarget = date.shiftedBy(10000.0);
783 FieldSpacecraftState<T> propagated = propagator.propagate(farTarget);
784 FieldPVCoordinates<T> pv = propagated.getPVCoordinates(FramesFactory.getITRF(IERSConventions.IERS_2010, true));
785 Assert.assertTrue(farTarget.durationFrom(propagated.getDate()).getReal() > 3000.0);
786 Assert.assertTrue(farTarget.durationFrom(propagated.getDate()).getReal() < 3500.0);
787 Assert.assertEquals(orbit.getA().getReal() * (1.0 - orbit.getE().getReal()), pv.getPosition().getNorm().getReal(), 410);
788 }
789
790 @Test
791 public void date() {
792 doDate(Decimal64Field.getInstance());
793
794 }
795
796 private <T extends CalculusFieldElement<T>> void doDate(Field<T> field) {
797 T zero = field.getZero();
798 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
799 final FieldKeplerianOrbit<T> orbit =
800 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
801 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
802 FieldEcksteinHechlerPropagator<T> propagator =
803 new FieldEcksteinHechlerPropagator<>(orbit, provider);
804 final FieldAbsoluteDate<T> stopDate = date.shiftedBy(500.0);
805 propagator.addEventDetector(new FieldDateDetector<>(stopDate));
806 FieldAbsoluteDate<T> farTarget = date.shiftedBy(10000.0);
807 FieldSpacecraftState<T> propagated = propagator.propagate(farTarget);
808 Assert.assertEquals(0, stopDate.durationFrom(propagated.getDate()).getReal(), 1.0e-10);
809 }
810
811 @Test
812 public void fixedStep() {
813
814 doFixedStep(Decimal64Field.getInstance());
815
816
817 }
818
819 private <T extends CalculusFieldElement<T>> void doFixedStep(Field<T> field) {
820 T zero = field.getZero();
821 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
822 final FieldKeplerianOrbit<T> orbit =
823 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
824 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
825 FieldEcksteinHechlerPropagator<T> propagator =
826 new FieldEcksteinHechlerPropagator<>(orbit, provider);
827 final T step = zero.add(100.0);
828 propagator.setStepHandler(step, new FieldOrekitFixedStepHandler<T>() {
829 private FieldAbsoluteDate<T> previous;
830 @Override
831 public void handleStep(FieldSpacecraftState<T> currentState) {
832 if (previous != null) {
833 Assert.assertEquals(step.getReal(), currentState.getDate().durationFrom(previous).getReal(), 1.0e-10);
834 }
835 previous = currentState.getDate();
836 }
837 });
838 FieldAbsoluteDate<T> farTarget = date.shiftedBy(10000.0);
839 propagator.propagate(farTarget);
840 }
841
842 @Test
843 public void setting() {
844
845 doSetting(Decimal64Field.getInstance());
846
847 }
848
849 private <T extends CalculusFieldElement<T>> void doSetting(Field<T> field) {
850 T zero = field.getZero();
851 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
852 final FieldKeplerianOrbit<T> orbit =
853 new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
854 FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
855 FieldEcksteinHechlerPropagator<T> propagator =
856 new FieldEcksteinHechlerPropagator<>(orbit, provider);
857 final OneAxisEllipsoid earthShape =
858 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
859 final TopocentricFrame topo =
860 new TopocentricFrame(earthShape, new GeodeticPoint(0.389, -2.962, 0), null);
861 FieldElevationDetector<T> detector = new FieldElevationDetector<>(zero.add(60), zero.add(1.0e-9), topo).withConstantElevation(0.09);
862 Assert.assertEquals(0.09, detector.getMinElevation(), 1.0e-12);
863 Assert.assertTrue(topo == detector.getTopocentricFrame());
864 propagator.addEventDetector(detector);
865
866 FieldAbsoluteDate<T> farTarget = date.shiftedBy(10000.0);
867 FieldSpacecraftState<T> propagated = propagator.propagate(farTarget);
868 final double elevation = topo.getElevation(propagated.getPVCoordinates().getPosition().toVector3D(),
869 propagated.getFrame(),
870 propagated.getDate().toAbsoluteDate());
871 final double zVelocity = propagated.getPVCoordinates(topo).getVelocity().getZ().getReal();
872 Assert.assertTrue(farTarget.durationFrom(propagated.getDate()).getReal() > 7800.0);
873 Assert.assertTrue("Incorrect value " + farTarget.durationFrom(propagated.getDate()) + " !< 7900",
874 farTarget.durationFrom(propagated.getDate()).getReal() < 7900.0);
875 Assert.assertEquals(0.09, elevation, 1.0e-11);
876 Assert.assertTrue(zVelocity < 0);
877 }
878
879 @Test
880 public void testIssue504() {
881 doTestIssue504(Decimal64Field.getInstance());
882 }
883
884 private <T extends CalculusFieldElement<T>> void doTestIssue504(Field<T> field) {
885 final T zero = field.getZero();
886
887 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
888 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
889 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2018, 07, 15), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
890 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
891 FramesFactory.getEME2000(),
892 initDate,
893 zero.add(3.986004415E14)));
894
895
896 final List<DSSTForceModel> models = new ArrayList<>();
897 models.add(new DSSTZonal(provider));
898 final FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(initialState, DEFAULT_LAW, models);
899
900
901 final FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(meanState.getOrbit(), provider, PropagationType.MEAN);
902 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate);
903
904
905 Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 18.0);
906 Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-6);
907 Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 5.0e-6);
908 Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-6);
909 Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 2.0e-6);
910 Assert.assertEquals(0.0,
911 FieldVector3D.distance(initialState.getPVCoordinates().getPosition(),
912 finalState.getPVCoordinates().getPosition()).getReal(),
913 11.4);
914 Assert.assertEquals(0.0,
915 FieldVector3D.distance(initialState.getPVCoordinates().getVelocity(),
916 finalState.getPVCoordinates().getVelocity()).getReal(),
917 4.2e-2);
918 }
919
920 @Test
921 public void testIssue504Bis() {
922 doTestIssue504Bis(Decimal64Field.getInstance());
923 }
924
925 private <T extends CalculusFieldElement<T>> void doTestIssue504Bis(Field<T> field) {
926 final T zero = field.getZero();
927
928 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
929 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
930 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2018, 07, 15), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
931 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
932 FramesFactory.getEME2000(),
933 initDate,
934 zero.add(3.986004415E14)));
935
936
937 final List<DSSTForceModel> models = new ArrayList<>();
938 models.add(new DSSTZonal(provider));
939 final FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(initialState, DEFAULT_LAW, models);
940
941
942 final FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(meanState.getOrbit(), DEFAULT_LAW, zero.add(498.5), provider, PropagationType.MEAN);
943 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate);
944
945
946 Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 18.0);
947 Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-6);
948 Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 5.0e-6);
949 Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-6);
950 Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 2.0e-6);
951 Assert.assertEquals(0.0,
952 FieldVector3D.distance(initialState.getPVCoordinates().getPosition(),
953 finalState.getPVCoordinates().getPosition()).getReal(),
954 11.4);
955 Assert.assertEquals(0.0,
956 FieldVector3D.distance(initialState.getPVCoordinates().getVelocity(),
957 finalState.getPVCoordinates().getVelocity()).getReal(),
958 4.2e-2);
959 }
960
961 private <T extends CalculusFieldElement<T>> T tangLEmLv(T Lv, T ex, T ey) {
962
963 return ey.multiply(Lv.cos()).subtract(ex.multiply(Lv.sin())).divide(
964 ex.multiply(Lv.cos()).add(1.0).add(ey.multiply(Lv.sin()).add(ex.negate().multiply(ex).add(1.0).subtract(
965 ey.multiply(ey)).sqrt())));
966
967 }
968
969 @After
970 public void tearDown() {
971 provider = null;
972 }
973
974 private UnnormalizedSphericalHarmonicsProvider provider;
975
976 }