1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.orekit.propagation.analytical;
19
20 import java.io.IOException;
21 import java.util.concurrent.TimeUnit;
22
23 import org.hipparchus.CalculusFieldElement;
24 import org.hipparchus.Field;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
27 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
28 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
29 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
30 import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
31 import org.hipparchus.stat.descriptive.rank.Max;
32 import org.hipparchus.stat.descriptive.rank.Min;
33 import org.hipparchus.util.Binary64Field;
34 import org.hipparchus.util.FastMath;
35 import org.hipparchus.util.MathUtils;
36 import org.junit.jupiter.api.AfterEach;
37 import org.junit.jupiter.api.Assertions;
38 import org.junit.jupiter.api.BeforeEach;
39 import org.junit.jupiter.api.Test;
40 import org.orekit.Utils;
41 import org.orekit.attitudes.AttitudeProvider;
42 import org.orekit.bodies.CelestialBodyFactory;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.errors.OrekitException;
45 import org.orekit.errors.OrekitMessages;
46 import org.orekit.forces.ForceModel;
47 import org.orekit.forces.drag.DragForce;
48 import org.orekit.forces.drag.IsotropicDrag;
49 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
50 import org.orekit.forces.gravity.potential.GravityFieldFactory;
51 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
52 import org.orekit.forces.gravity.potential.TideSystem;
53 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
54 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
55 import org.orekit.frames.Frame;
56 import org.orekit.frames.FramesFactory;
57 import org.orekit.models.earth.atmosphere.DTM2000;
58 import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
59 import org.orekit.orbits.FieldEquinoctialOrbit;
60 import org.orekit.orbits.FieldKeplerianOrbit;
61 import org.orekit.orbits.FieldOrbit;
62 import org.orekit.orbits.KeplerianOrbit;
63 import org.orekit.orbits.OrbitType;
64 import org.orekit.orbits.PositionAngleType;
65 import org.orekit.propagation.FieldSpacecraftState;
66 import org.orekit.propagation.PropagationType;
67 import org.orekit.propagation.Propagator;
68 import org.orekit.propagation.SpacecraftState;
69 import org.orekit.propagation.ToleranceProvider;
70 import org.orekit.propagation.analytical.tle.TLE;
71 import org.orekit.propagation.analytical.tle.TLEPropagator;
72 import org.orekit.propagation.numerical.FieldNumericalPropagator;
73 import org.orekit.propagation.numerical.NumericalPropagator;
74 import org.orekit.time.AbsoluteDate;
75 import org.orekit.time.FieldAbsoluteDate;
76 import org.orekit.time.TimeScale;
77 import org.orekit.time.TimeScalesFactory;
78 import org.orekit.utils.Constants;
79 import org.orekit.utils.FieldPVCoordinates;
80 import org.orekit.utils.IERSConventions;
81
82 public class FieldBrouwerLyddanePropagatorTest {
83
84 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
85
86 @Test
87 public void testIssue947() {
88 doTestIssue947(Binary64Field.getInstance());
89 }
90
91 private <T extends CalculusFieldElement<T>> void doTestIssue947(Field<T> field) {
92
93 TLE tleOrbit = new TLE("1 43196U 18015E 21055.59816856 .00000894 00000-0 38966-4 0 9996",
94 "2 43196 97.4662 188.8169 0016935 299.6845 60.2706 15.24746686170319");
95 Propagator propagator = TLEPropagator.selectExtrapolator(tleOrbit);
96
97
98 SpacecraftState tleState = propagator.getInitialState();
99 FieldOrbit<T> orbIni = OrbitType.KEPLERIAN.convertToFieldOrbit(field, tleState.getOrbit());
100 SpacecraftState tleStateAtDate = propagator.propagate(propagator.getInitialState().getDate().shiftedBy(3, TimeUnit.DAYS));
101 FieldOrbit<T> orbAtDate = OrbitType.KEPLERIAN.convertToFieldOrbit(field, tleStateAtDate.getOrbit());
102
103 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 0);
104
105 Assertions.assertDoesNotThrow(() -> {
106 FieldBrouwerLyddanePropagator.computeMeanOrbit(orbIni, provider,
107 provider.onDate(tleState.getDate()),
108 BrouwerLyddanePropagator.M2);
109 });
110
111 Assertions.assertDoesNotThrow(() -> {
112 FieldBrouwerLyddanePropagator.computeMeanOrbit(orbAtDate, provider,
113 provider.onDate(tleStateAtDate.getDate()),
114 BrouwerLyddanePropagator.M2);
115 });
116 }
117
118 @Test
119 public void sameDateCartesian() {
120 doSameDateCartesian(Binary64Field.getInstance());
121 }
122
123 private <T extends CalculusFieldElement<T>> void doSameDateCartesian(Field<T> field) {
124
125 T zero = field.getZero();
126 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
127
128
129
130
131 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
132 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6149822.));
133 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
134
135 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
136 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
137
138
139
140 FieldBrouwerLyddanePropagator<T> extrapolator =
141 new FieldBrouwerLyddanePropagator<T>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
142 FieldOrbit<T> finalOrbit = extrapolator.propagate(initDate).getOrbit();
143
144
145 Assertions.assertEquals(0.0,
146 FieldVector3D.distance(initialOrbit.getPosition(),
147 finalOrbit.getPosition()).getReal(),
148 4.8e-7);
149 Assertions.assertEquals(0.0,
150 FieldVector3D.distance(initialOrbit.getVelocity(),
151 finalOrbit.getVelocity()).getReal(),
152 2.9e-10);
153 Assertions.assertEquals(0.0,
154 finalOrbit.getA().getReal() - initialOrbit.getA().getReal(),
155 6.6e-9);
156
157 }
158
159 @Test
160 public void sameDateKeplerian() {
161 doSameDateKeplerian(Binary64Field.getInstance());
162 }
163
164 private <T extends CalculusFieldElement<T>> void doSameDateKeplerian(Field<T> field) {
165
166 T zero = field.getZero();
167 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
168
169
170 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
171 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(6767924.41), zero.add(.005), zero.add(1.7),zero.add( 2.1),
172 zero.add(2.9), zero.add(6.2), PositionAngleType.TRUE,
173 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
174
175 FieldBrouwerLyddanePropagator<T> extrapolator =
176 new FieldBrouwerLyddanePropagator<T>(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
177
178 FieldOrbit<T> finalOrbit = extrapolator.propagate(initDate).getOrbit();
179
180
181 Assertions.assertEquals(0.0,
182 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
183 finalOrbit.getPVCoordinates().getPosition()).getReal(),
184 1.2e-6);
185 Assertions.assertEquals(0.0,
186 FieldVector3D.distance(initialOrbit.getVelocity(),
187 finalOrbit.getVelocity()).getReal(),
188 7.0e-10);
189 Assertions.assertEquals(0.0,
190 finalOrbit.getA().getReal() - initialOrbit.getA().getReal(),
191 9.4e-9);
192 }
193
194 @Test
195 public void almostSphericalBody() {
196 doAlmostSphericalBody(Binary64Field.getInstance());
197 }
198
199 private <T extends CalculusFieldElement<T>> void doAlmostSphericalBody(Field<T> field) {
200
201 T zero = field.getZero();
202 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
203
204
205
206 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(8449822.));
207 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
208
209
210
211 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
212 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
213 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
214
215
216
217
218
219
220
221 UnnormalizedSphericalHarmonicsProvider kepProvider =
222 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
223 TideSystem.UNKNOWN,
224 new double[][] {
225 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
226 }, new double[][] {
227 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
228 });
229
230
231
232 FieldBrouwerLyddanePropagator<T> extrapolatorAna =
233 new FieldBrouwerLyddanePropagator<>(initialOrbit, zero.add(1000.0), kepProvider, BrouwerLyddanePropagator.M2);
234 FieldKeplerianPropagator<T> extrapolatorKep = new FieldKeplerianPropagator<>(initialOrbit);
235
236
237
238 double delta_t = 100.0;
239 FieldAbsoluteDate<T> extrapDate = date.shiftedBy(delta_t);
240
241 FieldOrbit<T> finalOrbitAna = extrapolatorAna.propagate(extrapDate).getOrbit();
242 FieldOrbit<T> finalOrbitKep = extrapolatorKep.propagate(extrapDate).getOrbit();
243
244 Assertions.assertEquals(finalOrbitAna.getDate().durationFrom(extrapDate).getReal(), 0.0,
245 Utils.epsilonTest);
246
247 Assertions.assertEquals(finalOrbitAna.getA().getReal(), finalOrbitKep.getA().getReal(), 10
248 * Utils.epsilonTest * finalOrbitKep.getA().getReal());
249 Assertions.assertEquals(finalOrbitAna.getEquinoctialEx().getReal(), finalOrbitKep.getEquinoctialEx().getReal(), Utils.epsilonE
250 * finalOrbitKep.getE().getReal());
251 Assertions.assertEquals(finalOrbitAna.getEquinoctialEy().getReal(), finalOrbitKep.getEquinoctialEy().getReal(), Utils.epsilonE
252 * finalOrbitKep.getE().getReal());
253 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx().getReal(), finalOrbitKep.getHx().getReal()),
254 finalOrbitKep.getHx().getReal(), Utils.epsilonAngle
255 * FastMath.abs(finalOrbitKep.getI().getReal()));
256 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy().getReal(), finalOrbitKep.getHy().getReal()),
257 finalOrbitKep.getHy().getReal(), Utils.epsilonAngle
258 * FastMath.abs(finalOrbitKep.getI().getReal()));
259 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv().getReal(), finalOrbitKep.getLv().getReal()),
260 finalOrbitKep.getLv().getReal(), Utils.epsilonAngle
261 * FastMath.abs(finalOrbitKep.getLv().getReal()));
262 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE().getReal(), finalOrbitKep.getLE().getReal()),
263 finalOrbitKep.getLE().getReal(), Utils.epsilonAngle
264 * FastMath.abs(finalOrbitKep.getLE().getReal()));
265 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM().getReal(), finalOrbitKep.getLM().getReal()),
266 finalOrbitKep.getLM().getReal(), Utils.epsilonAngle
267 * FastMath.abs(finalOrbitKep.getLM().getReal()));
268
269 }
270
271 @Test
272 public void compareToNumericalPropagation() {
273 doCompareToNumericalPropagation(Binary64Field.getInstance());
274 }
275
276 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagation(Field<T> field) {
277
278 T zero = field.getZero();
279 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
280 final Frame inertialFrame = FramesFactory.getEME2000();
281 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
282 double timeshift = 60000. ;
283
284
285 final double a = 24396159;
286 final double e = 0.01;
287 final double i = FastMath.toRadians(7);
288 final double omega = FastMath.toRadians(180);
289 final double raan = FastMath.toRadians(261);
290 final double lM = 0;
291 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
292 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
293 inertialFrame, initDate, zero.add(provider.getMu()));
294
295
296 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
297
298
299
300
301
302
303 final double minStep = 0.001;
304 final double maxstep = 1000.0;
305 final double positionTolerance = 10.0;
306 final OrbitType propagationType = OrbitType.KEPLERIAN;
307 final double[][] tolerances =
308 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit.toOrbit(), propagationType);
309 final AdaptiveStepsizeIntegrator integrator =
310 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
311
312
313 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
314 NumPropagator.setOrbitType(propagationType);
315
316 final ForceModel holmesFeatherstone =
317 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
318 NumPropagator.addForceModel(holmesFeatherstone);
319
320
321 NumPropagator.setInitialState(initialState.toSpacecraftState());
322
323
324 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
325 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
326
327
328
329
330
331 FieldBrouwerLyddanePropagator<T> BLextrapolator =
332 new FieldBrouwerLyddanePropagator<T>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
333
334 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
335 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
336
337 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.175);
338 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 3.2e-6);
339 Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 6.9e-8);
340 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
341 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0053);
342 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
343 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 1.2e-6);
344 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
345 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0052);
346 }
347
348 @Test
349 public void compareToNumericalPropagationWithDrag() {
350 doCompareToNumericalPropagationWithDrag(Binary64Field.getInstance());
351 }
352
353 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationWithDrag(Field<T> field) {
354
355 T zero = field.getZero();
356 final Frame inertialFrame = FramesFactory.getEME2000();
357 final TimeScale utc = TimeScalesFactory.getUTC();
358 final AbsoluteDate date = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
359 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(date, zero);
360 double timeshift = 60000. ;
361
362
363 final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3;
364 final double e = 0.01;
365 final double i = FastMath.toRadians(7);
366 final double omega = FastMath.toRadians(180);
367 final double raan = FastMath.toRadians(261);
368 final double lM = 0;
369 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
370 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
371 inertialFrame, initDate, zero.add(provider.getMu()));
372
373 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
374
375
376
377
378
379
380 final double minStep = 0.001;
381 final double maxstep = 1000.0;
382 final double positionTolerance = 10.0;
383 final OrbitType propagationType = OrbitType.KEPLERIAN;
384 final double[][] tolerances =
385 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit.toOrbit(), propagationType);
386 final AdaptiveStepsizeIntegrator integrator =
387 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
388
389
390 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
391 NumPropagator.setOrbitType(propagationType);
392
393
394 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
395 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
396 MarshallSolarActivityFutureEstimation msafe =
397 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
398 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
399 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
400
401
402 final ForceModel holmesFeatherstone =
403 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
404 final ForceModel drag = new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0));
405 NumPropagator.addForceModel(holmesFeatherstone);
406 NumPropagator.addForceModel(drag);
407
408
409 NumPropagator.setInitialState(initialState.toSpacecraftState());
410
411
412 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
413 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
414
415
416
417
418
419 FieldBrouwerLyddanePropagator<T> BLextrapolator =
420 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
421
422 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
423 KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
424
425
426 final double deltaSmaBefore = 15.81;
427 final double deltaEccBefore = 9.2681e-5;
428 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore);
429 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore);
430
431
432
433
434
435 double M2 = 1.0e-14;
436 BLextrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2);
437 BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
438 BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
439
440
441 final double deltaSmaAfter = 11.04;
442 final double deltaEccAfter = 9.2639e-5;
443 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter);
444 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter);
445 Assertions.assertTrue(deltaSmaAfter < deltaSmaBefore);
446 Assertions.assertTrue(deltaEccAfter < deltaEccBefore);
447 Assertions.assertEquals(M2, BLextrapolator.getM2(), Double.MIN_VALUE);
448
449 }
450
451 @Test
452 public void compareToNumericalPropagationMeanInitialOrbit() {
453 doCompareToNumericalPropagationMeanInitialOrbit(Binary64Field.getInstance());
454 }
455
456 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationMeanInitialOrbit(Field<T> field) {
457
458 T zero = field.getZero();
459 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
460 final Frame inertialFrame = FramesFactory.getEME2000();
461 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
462 double timeshift = 60000. ;
463
464
465 final double a = 24396159;
466 final double e = 0.01;
467 final double i = FastMath.toRadians(7);
468 final double omega = FastMath.toRadians(180);
469 final double raan = FastMath.toRadians(261);
470 final double lM = 0;
471 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
472 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
473 inertialFrame, initDate, zero.add(provider.getMu()));
474
475 FieldBrouwerLyddanePropagator<T> BLextrapolator =
476 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider),
477 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
478 FieldSpacecraftState<T> initialOsculatingState = BLextrapolator.propagate(initDate);
479 final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit().toOrbit());
480
481 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
482 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
483
484
485
486
487
488
489 final double minStep = 0.001;
490 final double maxstep = 1000.0;
491 final double positionTolerance = 10.0;
492 final OrbitType propagationType = OrbitType.KEPLERIAN;
493 final double[][] tolerances =
494 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(InitOrbit, propagationType);
495 final AdaptiveStepsizeIntegrator integrator =
496 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
497
498
499 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
500 NumPropagator.setOrbitType(propagationType);
501
502 final ForceModel holmesFeatherstone =
503 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
504 NumPropagator.addForceModel(holmesFeatherstone);
505
506
507 NumPropagator.setInitialState(initialOsculatingState.toSpacecraftState());
508
509
510 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
511 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
512
513
514 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.174);
515 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 3.2e-6);
516 Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 6.9e-8);
517 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
518 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0053);
519 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
520 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 1.2e-6);
521 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
522 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0052);
523 }
524
525 @Test
526 public void compareToNumericalPropagationResetInitialIntermediate() {
527 doCompareToNumericalPropagationResetInitialIntermediate(Binary64Field.getInstance());
528 }
529
530 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationResetInitialIntermediate(Field<T> field) {
531
532 T zero = field.getZero();
533 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
534 final Frame inertialFrame = FramesFactory.getEME2000();
535 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
536 double timeshift = 60000. ;
537
538
539 final double a = 24396159;
540 final double e = 0.01;
541 final double i = FastMath.toRadians(7);
542 final double omega = FastMath.toRadians(180);
543 final double raan = FastMath.toRadians(261);
544 final double lM = 0;
545 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
546 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
547 inertialFrame, initDate, zero.add(provider.getMu()));
548
549 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
550
551
552
553
554
555 FieldBrouwerLyddanePropagator<T> BLextrapolator1 =
556 new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, zero.add(Propagator.DEFAULT_MASS), GravityFieldFactory.getUnnormalizedProvider(provider),
557 PropagationType.OSCULATING, BrouwerLyddanePropagator.M2);
558
559
560
561
562 FieldBrouwerLyddanePropagator<T> BLextrapolator2 =
563 new FieldBrouwerLyddanePropagator<>( new FieldKeplerianOrbit<>(zero.add(a + 3000), zero.add(e + 0.001),
564 zero.add(i - FastMath.toRadians(12.0)), zero.add(omega),
565 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
566 inertialFrame, initDate, zero.add(provider.getMu())),DEFAULT_LAW,
567 zero.add(Propagator.DEFAULT_MASS),
568 GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
569
570 BLextrapolator2.resetInitialState(initialState);
571
572 FieldSpacecraftState<T> BLFinalState1 = BLextrapolator1.propagate(initDate.shiftedBy(timeshift));
573 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit().toOrbit());
574 FieldSpacecraftState<T> BLFinalState2 = BLextrapolator2.propagate(initDate.shiftedBy(timeshift));
575 BLextrapolator2.resetIntermediateState(BLFinalState1, true);
576 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit().toOrbit());
577
578 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
579 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
580 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
581 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
582 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
583 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
584 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
585 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
586 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
587
588 }
589
590 @Test
591 public void compareConstructors() {
592 doCompareConstructors(Binary64Field.getInstance());
593 }
594
595 private <T extends CalculusFieldElement<T>> void doCompareConstructors(Field<T> field) {
596
597 T zero = field.getZero();
598 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
599 final Frame inertialFrame = FramesFactory.getEME2000();
600 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
601 double timeshift = 600. ;
602
603
604 final double a = 24396159;
605 final double e = 0.01;
606 final double i = FastMath.toRadians(7);
607 final double omega = FastMath.toRadians(180);
608 final double raan = FastMath.toRadians(261);
609 final double lM = 0;
610 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
611 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
612 inertialFrame, initDate, zero.add(provider.getMu()));
613
614
615 FieldBrouwerLyddanePropagator<T> BLPropagator1 = new FieldBrouwerLyddanePropagator<T>(initialOrbit, DEFAULT_LAW,
616 provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
617 FieldBrouwerLyddanePropagator<T> BLPropagator2 = new FieldBrouwerLyddanePropagator<>(initialOrbit,
618 provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
619 FieldBrouwerLyddanePropagator<T> BLPropagator3 = new FieldBrouwerLyddanePropagator<>(initialOrbit,
620 zero.add(Propagator.DEFAULT_MASS), provider.getAe(), zero.add(provider.getMu()), -1.08263e-3,
621 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
622
623 FieldSpacecraftState<T> BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift));
624 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit().toOrbit());
625 FieldSpacecraftState<T> BLFinalState2 = BLPropagator2.propagate(initDate.shiftedBy(timeshift));
626 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit().toOrbit());
627 FieldSpacecraftState<T> BLFinalState3 = BLPropagator3.propagate(initDate.shiftedBy(timeshift));
628 final KeplerianOrbit BLOrbit3 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState3.getOrbit().toOrbit());
629
630 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
631 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
632 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
633 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
634 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
635 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
636 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
637 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
638 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
639
640 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit3.getA(), 0.0);
641 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit3.getE(), 0.0);
642 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit3.getI(), 0.0);
643 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
644 MathUtils.normalizeAngle(BLOrbit3.getPerigeeArgument(), FastMath.PI), 0.0);
645 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
646 MathUtils.normalizeAngle(BLOrbit3.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
647 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
648 MathUtils.normalizeAngle(BLOrbit3.getTrueAnomaly(), FastMath.PI), 0.0);
649
650 }
651
652 @Test
653 public void undergroundOrbit() {
654 doUndergroundOrbit(Binary64Field.getInstance());
655 }
656
657 private <T extends CalculusFieldElement<T>> void doUndergroundOrbit(Field<T> field) {
658
659 T zero = field.getZero();
660 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
661 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
662
663
664 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
665 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(800.0), zero.add(100.0));
666
667 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<T>(position, velocity),
668 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
669
670
671 try {
672
673 FieldBrouwerLyddanePropagator<T> extrapolator =
674 new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, provider.getAe(), zero.add(provider.getMu()),
675 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
676
677
678
679 double delta_t = 0.0;
680 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
681 extrapolator.propagate(extrapDate);
682
683 } catch (OrekitException oe) {
684 Assertions.assertEquals(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, oe.getSpecifier());
685 }
686 }
687
688 @Test
689 public void tooEllipticalOrbit() {
690 doTooEllipticalOrbit(Binary64Field.getInstance());
691 }
692
693 private <T extends CalculusFieldElement<T>> void doTooEllipticalOrbit(Field<T> field) {
694
695 T zero = field.getZero();
696 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
697 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
698 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(67679244.0), zero.add(1.0), zero.add(1.85850),
699 zero.add(2.1), zero.add(2.9), zero.add(6.2), PositionAngleType.TRUE,
700 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
701 try {
702
703
704 FieldBrouwerLyddanePropagator<T> extrapolator =
705 new FieldBrouwerLyddanePropagator<>(initialOrbit, provider.getAe(), zero.add(provider.getMu()),
706 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
707
708
709
710 double delta_t = 0.0;
711 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
712 extrapolator.propagate(extrapDate);
713
714 } catch (OrekitException oe) {
715 Assertions.assertEquals(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL, oe.getSpecifier());
716 }
717 }
718
719 @Test
720 public void criticalInclination() {
721 doCriticalInclination(Binary64Field.getInstance());
722 }
723
724 private <T extends CalculusFieldElement<T>> void doCriticalInclination(Field<T> field) {
725
726 final Frame inertialFrame = FramesFactory.getEME2000();
727
728
729 final double a = 24396159;
730 final double e = 0.01;
731 final double i = FastMath.acos(1.0 / FastMath.sqrt(5.0));
732 final double omega = FastMath.toRadians(180);
733 final double raan = FastMath.toRadians(261);
734 final double lV = 0;
735
736 T zero = field.getZero();
737 FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field);
738 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
739 zero.add(raan), zero.add(lV), PositionAngleType.TRUE,
740 inertialFrame, initDate, zero.add(provider.getMu()));
741
742
743
744 FieldBrouwerLyddanePropagator<T> extrapolator =
745 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
746
747
748
749 final FieldOrbit<T> finalOrbit = extrapolator.propagate(initDate).getOrbit();
750
751
752
753 Assertions.assertEquals(0.0,
754 FieldVector3D.distance(initialOrbit.getPosition(),
755 finalOrbit.getPosition()).getReal(),
756 1.5e-8);
757 Assertions.assertEquals(0.0,
758 FieldVector3D.distance(initialOrbit.getVelocity(),
759 finalOrbit.getVelocity()).getReal(),
760 2.7e-12);
761 Assertions.assertEquals(0.0,
762 finalOrbit.getA().getReal() - initialOrbit.getA().getReal(),
763 7.5e-9);
764 }
765
766 @Test
767 public void testUnableToComputeBLMeanParameters() {
768 OrekitException oe = Assertions.assertThrows(OrekitException.class, () -> {
769 doTestUnableToComputeBLMeanParameters(Binary64Field.getInstance());
770 });
771 Assertions.assertTrue(oe.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
772 }
773
774 private <T extends CalculusFieldElement<T>> void doTestUnableToComputeBLMeanParameters(Field<T> field) {
775 final Frame eci = FramesFactory.getEME2000();
776 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
777
778
779 final double a = 24396159;
780 final double e = 0.9;
781 final double i = FastMath.toRadians(7);
782 final double pa = FastMath.toRadians(180);
783 final double raan = FastMath.toRadians(261);
784 final double lM = FastMath.toRadians(0);
785 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(field,
786 new KeplerianOrbit(a, e, i, pa, raan, lM,
787 PositionAngleType.MEAN,
788 eci, date, provider.getMu()));
789
790
791 final UnnormalizedSphericalHarmonicsProvider up = GravityFieldFactory.getUnnormalizedProvider(provider);
792 final UnnormalizedSphericalHarmonics uh = up.onDate(date);
793 final double eps = 1.e-13;
794 final int maxIter = 10;
795 FieldBrouwerLyddanePropagator.computeMeanOrbit(orbit, up, uh, BrouwerLyddanePropagator.M2,
796 eps, maxIter);
797 }
798
799 @Test
800 public void testMeanComparisonWithNonField() {
801 doTestMeanComparisonWithNonField(Binary64Field.getInstance());
802 }
803
804 private <T extends CalculusFieldElement<T>> void doTestMeanComparisonWithNonField(Field<T> field) {
805
806 T zero = field.getZero();
807 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
808 final Frame inertialFrame = FramesFactory.getEME2000();
809 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
810 double timeshift = -59400.0;
811
812
813 final double a = 24396159;
814 final double e = 0.9;
815 final double i = FastMath.toRadians(7);
816 final double omega = FastMath.toRadians(180);
817 final double raan = FastMath.toRadians(261);
818 final double lM = FastMath.toRadians(0);
819 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
820 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
821 inertialFrame, initDate, zero.add(provider.getMu()));
822
823
824 final FieldSpacecraftState<T> initialStateField = new FieldSpacecraftState<>(initialOrbit);
825 final SpacecraftState initialState = initialStateField.toSpacecraftState();
826
827
828 final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(),
829 GravityFieldFactory.getUnnormalizedProvider(provider),
830 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
831 final FieldSpacecraftState<T> finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift));
832 final FieldKeplerianOrbit<T> finalOrbitField = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit());
833 final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit();
834
835
836 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(),
837 GravityFieldFactory.getUnnormalizedProvider(provider),
838 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
839 final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift));
840 final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
841
842 Assertions.assertEquals(finalOrbitFieldReal.getA(), finalOrbit.getA(), Double.MIN_VALUE);
843 Assertions.assertEquals(finalOrbitFieldReal.getE(), finalOrbit.getE(), Double.MIN_VALUE);
844 Assertions.assertEquals(finalOrbitFieldReal.getI(), finalOrbit.getI(), Double.MIN_VALUE);
845 Assertions.assertEquals(finalOrbitFieldReal.getRightAscensionOfAscendingNode(), finalOrbit.getRightAscensionOfAscendingNode(), Double.MIN_VALUE);
846 Assertions.assertEquals(finalOrbitFieldReal.getPerigeeArgument(), finalOrbit.getPerigeeArgument(), Double.MIN_VALUE);
847 Assertions.assertEquals(finalOrbitFieldReal.getMeanAnomaly(), finalOrbit.getMeanAnomaly(), Double.MIN_VALUE);
848 Assertions.assertEquals(0.0, finalOrbitFieldReal.getPosition().distance(finalOrbit.getPosition()), Double.MIN_VALUE);
849 Assertions.assertEquals(0.0, finalOrbitFieldReal.getVelocity().distance(finalOrbit.getVelocity()), Double.MIN_VALUE);
850 }
851
852 @Test
853 public void testOsculatingComparisonWithNonField() {
854 doTestOsculatingComparisonWithNonField(Binary64Field.getInstance());
855 }
856
857 private <T extends CalculusFieldElement<T>> void doTestOsculatingComparisonWithNonField(Field<T> field) {
858
859 T zero = field.getZero();
860 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
861 final Frame inertialFrame = FramesFactory.getEME2000();
862 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
863 double timeshift = 60000.0;
864
865
866 final double a = 24396159;
867 final double e = 0.01;
868 final double i = FastMath.toRadians(7);
869 final double omega = FastMath.toRadians(180);
870 final double raan = FastMath.toRadians(261);
871 final double lM = FastMath.toRadians(0);
872 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
873 zero.add(raan), zero.add(lM), PositionAngleType.TRUE,
874 inertialFrame, initDate, zero.add(provider.getMu()));
875
876
877 final FieldSpacecraftState<T> initialStateField = new FieldSpacecraftState<>(initialOrbit);
878 final SpacecraftState initialState = initialStateField.toSpacecraftState();
879
880
881 final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(),
882 GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
883 final FieldSpacecraftState<T> finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift));
884 final FieldKeplerianOrbit<T> finalOrbitField = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit());
885 final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit();
886
887
888 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(),
889 GravityFieldFactory.getUnnormalizedProvider(provider),
890 BrouwerLyddanePropagator.M2);
891 final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift));
892 final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
893
894 Assertions.assertEquals(finalOrbitFieldReal.getA(), finalOrbit.getA(), Double.MIN_VALUE);
895 Assertions.assertEquals(finalOrbitFieldReal.getE(), finalOrbit.getE(), Double.MIN_VALUE);
896 Assertions.assertEquals(finalOrbitFieldReal.getI(), finalOrbit.getI(), Double.MIN_VALUE);
897 Assertions.assertEquals(finalOrbitFieldReal.getRightAscensionOfAscendingNode(), finalOrbit.getRightAscensionOfAscendingNode(), Double.MIN_VALUE);
898 Assertions.assertEquals(finalOrbitFieldReal.getPerigeeArgument(), finalOrbit.getPerigeeArgument(), Double.MIN_VALUE);
899 Assertions.assertEquals(finalOrbitFieldReal.getMeanAnomaly(), finalOrbit.getMeanAnomaly(), Double.MIN_VALUE);
900 Assertions.assertEquals(0.0, finalOrbitFieldReal.getPosition().distance(finalOrbit.getPosition()), Double.MIN_VALUE);
901 Assertions.assertEquals(0.0, finalOrbitFieldReal.getVelocity().distance(finalOrbit.getVelocity()), Double.MIN_VALUE);
902 }
903
904 @Test
905 public void testMeanOrbit() throws IOException {
906 doTestMeanOrbit(Binary64Field.getInstance());
907 }
908
909 private <T extends CalculusFieldElement<T>> void doTestMeanOrbit(Field<T> field) {
910 T zero = field.getZero();
911 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
912 final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
913 final FieldKeplerianOrbit<T> initialOsculating =
914 new FieldKeplerianOrbit<>(zero.newInstance(7.8e6), zero.newInstance(0.032), zero.newInstance(0.4),
915 zero.newInstance(0.1), zero.newInstance(0.2), zero.newInstance(0.3),
916 PositionAngleType.TRUE,
917 FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
918 final UnnormalizedSphericalHarmonics ush = ushp.onDate(initialOsculating.getDate().toAbsoluteDate());
919
920
921
922 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.1).getTolerances(initialOsculating, OrbitType.KEPLERIAN);
923 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, tol[0], tol[1]);
924 integrator.setInitialStepSize(60);
925 FieldNumericalPropagator<T> num = new FieldNumericalPropagator<>(field, integrator);
926 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
927 num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
928 num.setInitialState(new FieldSpacecraftState<>(initialOsculating));
929 num.setOrbitType(OrbitType.KEPLERIAN);
930 num.setPositionAngleType(initialOsculating.getCachedPositionAngleType());
931 final StorelessUnivariateStatistic oscMin = new Min();
932 final StorelessUnivariateStatistic oscMax = new Max();
933 final StorelessUnivariateStatistic meanMin = new Min();
934 final StorelessUnivariateStatistic meanMax = new Max();
935 num.getMultiplexer().add(zero.newInstance(60), state -> {
936 final FieldOrbit<T> osc = state.getOrbit();
937 oscMin.increment(osc.getA().getReal());
938 oscMax.increment(osc.getA().getReal());
939
940 final FieldOrbit<T> mean = FieldBrouwerLyddanePropagator.computeMeanOrbit(state.getOrbit(),
941 ushp, ush,
942 BrouwerLyddanePropagator.M2);
943 meanMin.increment(mean.getA().getReal());
944 meanMax.increment(mean.getA().getReal());
945 });
946 num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
947
948 Assertions.assertEquals(3188.347, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
949 Assertions.assertEquals( 25.794, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
950 }
951
952 @Test
953 public void testGeostationaryOrbit() {
954 doTestGeostationaryOrbit(Binary64Field.getInstance());
955 }
956
957 private <T extends CalculusFieldElement<T>> void doTestGeostationaryOrbit(Field<T> field) {
958
959 final Frame eci = FramesFactory.getEME2000();
960 final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
961
962 T zero = field.getZero();
963 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
964
965
966 final double sma = FastMath.cbrt(Constants.IERS2010_EARTH_MU /
967 Constants.IERS2010_EARTH_ANGULAR_VELOCITY /
968 Constants.IERS2010_EARTH_ANGULAR_VELOCITY);
969 final double ecc = 0.0;
970 final double inc = 0.0;
971 final double pa = 0.0;
972 final double raan = 0.0;
973 final double lV = 0.0;
974 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(sma), zero.add(ecc), zero.add(inc),
975 zero.add(pa), zero.add(raan), zero.add(lV),
976 PositionAngleType.TRUE,
977 eci, date, zero.add(provider.getMu()));
978
979
980 FieldBrouwerLyddanePropagator<T> fbl =
981 new FieldBrouwerLyddanePropagator<>(orbit, ushp, PropagationType.MEAN,
982 BrouwerLyddanePropagator.M2);
983
984
985 final FieldSpacecraftState<T> state = fbl.propagate(date.shiftedBy(Constants.JULIAN_DAY));
986 final KeplerianOrbit orbOsc = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit().toOrbit());
987
988
989 Assertions.assertTrue(Double.isFinite(orbOsc.getA()));
990 Assertions.assertTrue(Double.isFinite(orbOsc.getE()));
991 Assertions.assertTrue(Double.isFinite(orbOsc.getI()));
992 Assertions.assertTrue(Double.isFinite(orbOsc.getPerigeeArgument()));
993 Assertions.assertTrue(Double.isFinite(orbOsc.getRightAscensionOfAscendingNode()));
994 Assertions.assertTrue(Double.isFinite(orbOsc.getTrueAnomaly()));
995
996
997 FieldBrouwerLyddanePropagator<T> fbl2 =
998 new FieldBrouwerLyddanePropagator<>(orbit, ushp, PropagationType.OSCULATING,
999 BrouwerLyddanePropagator.M2);
1000
1001
1002 final FieldSpacecraftState<T> state2 = fbl2.propagate(date.shiftedBy(Constants.JULIAN_DAY));
1003 final KeplerianOrbit orbOsc2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state2.getOrbit().toOrbit());
1004
1005
1006 Assertions.assertTrue(Double.isFinite(orbOsc2.getA()));
1007 Assertions.assertTrue(Double.isFinite(orbOsc2.getE()));
1008 Assertions.assertTrue(Double.isFinite(orbOsc2.getI()));
1009 Assertions.assertTrue(Double.isFinite(orbOsc2.getPerigeeArgument()));
1010 Assertions.assertTrue(Double.isFinite(orbOsc2.getRightAscensionOfAscendingNode()));
1011 Assertions.assertTrue(Double.isFinite(orbOsc2.getTrueAnomaly()));
1012 }
1013
1014 @BeforeEach
1015 public void setUp() {
1016 Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
1017 provider = GravityFieldFactory.getNormalizedProvider(5, 0);
1018 }
1019
1020 @AfterEach
1021 public void tearDown() {
1022 provider = null;
1023 }
1024
1025 private NormalizedSphericalHarmonicsProvider provider;
1026 }