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