1 package org.orekit.propagation.analytical;
2
3 import org.hipparchus.CalculusFieldElement;
4 import org.hipparchus.Field;
5 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
6 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
7 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
8 import org.hipparchus.util.Decimal64Field;
9 import org.hipparchus.util.FastMath;
10 import org.hipparchus.util.MathUtils;
11 import org.junit.After;
12 import org.junit.Assert;
13 import org.junit.Before;
14 import org.junit.Test;
15 import org.orekit.Utils;
16 import org.orekit.attitudes.AttitudeProvider;
17 import org.orekit.bodies.CelestialBodyFactory;
18 import org.orekit.bodies.OneAxisEllipsoid;
19 import org.orekit.data.DataContext;
20 import org.orekit.errors.OrekitException;
21 import org.orekit.errors.OrekitMessages;
22 import org.orekit.forces.ForceModel;
23 import org.orekit.forces.drag.DragForce;
24 import org.orekit.forces.drag.IsotropicDrag;
25 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
26 import org.orekit.forces.gravity.potential.GravityFieldFactory;
27 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
28 import org.orekit.forces.gravity.potential.TideSystem;
29 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
30 import org.orekit.frames.Frame;
31 import org.orekit.frames.FramesFactory;
32 import org.orekit.models.earth.atmosphere.DTM2000;
33 import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
34 import org.orekit.orbits.FieldEquinoctialOrbit;
35 import org.orekit.orbits.FieldKeplerianOrbit;
36 import org.orekit.orbits.FieldOrbit;
37 import org.orekit.orbits.KeplerianOrbit;
38 import org.orekit.orbits.OrbitType;
39 import org.orekit.orbits.PositionAngle;
40 import org.orekit.propagation.FieldSpacecraftState;
41 import org.orekit.propagation.PropagationType;
42 import org.orekit.propagation.Propagator;
43 import org.orekit.propagation.SpacecraftState;
44 import org.orekit.propagation.numerical.NumericalPropagator;
45 import org.orekit.time.AbsoluteDate;
46 import org.orekit.time.FieldAbsoluteDate;
47 import org.orekit.time.TimeScale;
48 import org.orekit.time.TimeScalesFactory;
49 import org.orekit.utils.Constants;
50 import org.orekit.utils.FieldPVCoordinates;
51 import org.orekit.utils.IERSConventions;
52
53 public class FieldBrouwerLyddanePropagatorTest {
54 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
55
56 @Test
57 public void sameDateCartesian() {
58 doSameDateCartesian(Decimal64Field.getInstance());
59 }
60 private <T extends CalculusFieldElement<T>> void doSameDateCartesian(Field<T> field) {
61
62 T zero = field.getZero();
63 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
64
65
66
67
68 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
69 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6149822.));
70 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
71
72 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
73 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
74
75
76
77 FieldBrouwerLyddanePropagator<T> extrapolator =
78 new FieldBrouwerLyddanePropagator<T>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
79 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(initDate);
80
81
82 Assert.assertEquals(0.0,
83 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
84 finalOrbit.getPVCoordinates().getPosition()).getReal(),
85 1.0e-8);
86
87 Assert.assertEquals(0.0,
88 FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
89 finalOrbit.getPVCoordinates().getVelocity()).getReal(),
90 1.0e-11);
91 Assert.assertEquals(0.0, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.0);
92
93 }
94
95
96 @Test
97 public void sameDateKeplerian() {
98 doSameDateKeplerian(Decimal64Field.getInstance());
99 }
100 private <T extends CalculusFieldElement<T>> void doSameDateKeplerian(Field<T> field) {
101
102 T zero = field.getZero();
103 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
104
105
106 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
107 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(6767924.41), zero.add(.005), zero.add(1.7),zero.add( 2.1),
108 zero.add(2.9), zero.add(6.2), PositionAngle.TRUE,
109 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
110
111 FieldBrouwerLyddanePropagator<T> extrapolator =
112 new FieldBrouwerLyddanePropagator<T>(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
113
114 FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(initDate);
115
116
117 Assert.assertEquals(0.0,
118 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
119 finalOrbit.getPVCoordinates().getPosition()).getReal(),
120 1.1e-8);
121
122 Assert.assertEquals(0.0,
123 FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
124 finalOrbit.getPVCoordinates().getVelocity()).getReal(),
125 1.4e-11);
126 Assert.assertEquals(0.0, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.0);
127 }
128
129
130 @Test
131 public void almostSphericalBody() {
132 doAlmostSphericalBody(Decimal64Field.getInstance());
133 }
134 private <T extends CalculusFieldElement<T>> void doAlmostSphericalBody(Field<T> field) {
135
136 T zero = field.getZero();
137 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
138
139
140
141 FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(8449822.));
142 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
143
144
145
146 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
147 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
148 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
149
150
151
152
153
154
155
156 UnnormalizedSphericalHarmonicsProvider kepProvider =
157 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
158 TideSystem.UNKNOWN,
159 new double[][] {
160 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
161 }, new double[][] {
162 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
163 });
164
165
166
167 FieldBrouwerLyddanePropagator<T> extrapolatorAna =
168 new FieldBrouwerLyddanePropagator<>(initialOrbit, zero.add(1000.0), kepProvider, BrouwerLyddanePropagator.M2);
169 FieldKeplerianPropagator<T> extrapolatorKep = new FieldKeplerianPropagator<>(initialOrbit);
170
171
172
173 double delta_t = 100.0;
174 FieldAbsoluteDate<T> extrapDate = date.shiftedBy(delta_t);
175
176 FieldSpacecraftState<T> finalOrbitAna = extrapolatorAna.propagate(extrapDate);
177 FieldSpacecraftState<T> finalOrbitKep = extrapolatorKep.propagate(extrapDate);
178
179 Assert.assertEquals(finalOrbitAna.getDate().durationFrom(extrapDate).getReal(), 0.0,
180 Utils.epsilonTest);
181
182 Assert.assertEquals(finalOrbitAna.getA().getReal(), finalOrbitKep.getA().getReal(), 10
183 * Utils.epsilonTest * finalOrbitKep.getA().getReal());
184 Assert.assertEquals(finalOrbitAna.getEquinoctialEx().getReal(), finalOrbitKep.getEquinoctialEx().getReal(), Utils.epsilonE
185 * finalOrbitKep.getE().getReal());
186 Assert.assertEquals(finalOrbitAna.getEquinoctialEy().getReal(), finalOrbitKep.getEquinoctialEy().getReal(), Utils.epsilonE
187 * finalOrbitKep.getE().getReal());
188 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx().getReal(), finalOrbitKep.getHx().getReal()),
189 finalOrbitKep.getHx().getReal(), Utils.epsilonAngle
190 * FastMath.abs(finalOrbitKep.getI().getReal()));
191 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy().getReal(), finalOrbitKep.getHy().getReal()),
192 finalOrbitKep.getHy().getReal(), Utils.epsilonAngle
193 * FastMath.abs(finalOrbitKep.getI().getReal()));
194 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv().getReal(), finalOrbitKep.getLv().getReal()),
195 finalOrbitKep.getLv().getReal(), Utils.epsilonAngle
196 * FastMath.abs(finalOrbitKep.getLv().getReal()));
197 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE().getReal(), finalOrbitKep.getLE().getReal()),
198 finalOrbitKep.getLE().getReal(), Utils.epsilonAngle
199 * FastMath.abs(finalOrbitKep.getLE().getReal()));
200 Assert.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM().getReal(), finalOrbitKep.getLM().getReal()),
201 finalOrbitKep.getLM().getReal(), Utils.epsilonAngle
202 * FastMath.abs(finalOrbitKep.getLM().getReal()));
203
204 }
205
206
207 @Test
208 public void compareToNumericalPropagation() {
209 doCompareToNumericalPropagation(Decimal64Field.getInstance());
210 }
211 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagation(Field<T> field) {
212
213 T zero = field.getZero();
214 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
215 final Frame inertialFrame = FramesFactory.getEME2000();
216 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
217 double timeshift = 60000. ;
218
219
220 final double a = 24396159;
221 final double e = 0.01;
222 final double i = FastMath.toRadians(7);
223 final double omega = FastMath.toRadians(180);
224 final double raan = FastMath.toRadians(261);
225 final double lM = 0;
226 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
227 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
228 inertialFrame, initDate, zero.add(provider.getMu()));
229
230
231 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
232
233
234
235
236
237
238 final double minStep = 0.001;
239 final double maxstep = 1000.0;
240 final double positionTolerance = 10.0;
241 final OrbitType propagationType = OrbitType.KEPLERIAN;
242 final double[][] tolerances =
243 NumericalPropagator.tolerances(positionTolerance, initialOrbit.toOrbit(), propagationType);
244 final AdaptiveStepsizeIntegrator integrator =
245 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
246
247
248 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
249 NumPropagator.setOrbitType(propagationType);
250
251 final ForceModel holmesFeatherstone =
252 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
253 NumPropagator.addForceModel(holmesFeatherstone);
254
255
256 NumPropagator.setInitialState(initialState.toSpacecraftState());
257
258
259 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
260 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
261
262
263
264
265
266 FieldBrouwerLyddanePropagator<T> BLextrapolator =
267 new FieldBrouwerLyddanePropagator<T>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
268
269 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
270 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
271
272
273 Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.072);
274 Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
275 Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
276 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
277 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.119);
278 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
279 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.000072);
280 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
281 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12);
282 }
283
284 @Test
285 public void compareToNumericalPropagationWithDrag() {
286 doCompareToNumericalPropagationWithDrag(Decimal64Field.getInstance());
287 }
288
289 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationWithDrag(Field<T> field) {
290
291 T zero = field.getZero();
292 final Frame inertialFrame = FramesFactory.getEME2000();
293 final TimeScale utc = TimeScalesFactory.getUTC();
294 final AbsoluteDate date = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
295 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(date, zero);
296 double timeshift = 60000. ;
297
298
299 final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3;
300 final double e = 0.01;
301 final double i = FastMath.toRadians(7);
302 final double omega = FastMath.toRadians(180);
303 final double raan = FastMath.toRadians(261);
304 final double lM = 0;
305 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
306 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
307 inertialFrame, initDate, zero.add(provider.getMu()));
308
309 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
310
311
312
313
314
315
316 final double minStep = 0.001;
317 final double maxstep = 1000.0;
318 final double positionTolerance = 10.0;
319 final OrbitType propagationType = OrbitType.KEPLERIAN;
320 final double[][] tolerances =
321 NumericalPropagator.tolerances(positionTolerance, initialOrbit.toOrbit(), propagationType);
322 final AdaptiveStepsizeIntegrator integrator =
323 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
324
325
326 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
327 NumPropagator.setOrbitType(propagationType);
328
329
330 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
331 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
332 MarshallSolarActivityFutureEstimation msafe =
333 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
334 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
335 DataContext.getDefault().getDataProvidersManager().feed(msafe.getSupportedNames(), msafe);
336 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
337
338
339 final ForceModel holmesFeatherstone =
340 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
341 final ForceModel drag =
342 new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0));
343 NumPropagator.addForceModel(holmesFeatherstone);
344 NumPropagator.addForceModel(drag);
345
346
347 NumPropagator.setInitialState(initialState.toSpacecraftState());
348
349
350 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
351 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
352
353
354
355
356
357 FieldBrouwerLyddanePropagator<T> BLextrapolator =
358 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
359
360 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
361 KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
362
363
364 final double deltaSmaBefore = 20.44;
365 final double deltaEccBefore = 1.0125e-4;
366 Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore);
367 Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore);
368
369
370
371
372
373 double M2 = 1.0e-14;
374 BLextrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2);
375 BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
376 BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
377
378
379 final double deltaSmaAfter = 15.66;
380 final double deltaEccAfter = 1.0121e-4;
381 Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter);
382 Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter);
383 Assert.assertTrue(deltaSmaAfter < deltaSmaBefore);
384 Assert.assertTrue(deltaEccAfter < deltaEccBefore);
385 Assert.assertEquals(M2, BLextrapolator.getM2(), Double.MIN_VALUE);
386
387 }
388
389 @Test
390 public void compareToNumericalPropagationMeanInitialOrbit() {
391 doCompareToNumericalPropagationMeanInitialOrbit(Decimal64Field.getInstance());
392 }
393 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationMeanInitialOrbit(Field<T> field) {
394
395 T zero = field.getZero();
396 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
397 final Frame inertialFrame = FramesFactory.getEME2000();
398 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
399 double timeshift = 60000. ;
400
401
402
403 final double a = 24396159;
404 final double e = 0.01;
405 final double i = FastMath.toRadians(7);
406 final double omega = FastMath.toRadians(180);
407 final double raan = FastMath.toRadians(261);
408 final double lM = 0;
409 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
410 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
411 inertialFrame, initDate, zero.add(provider.getMu()));
412
413 FieldBrouwerLyddanePropagator<T> BLextrapolator =
414 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider),
415 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
416 FieldSpacecraftState<T> initialOsculatingState = BLextrapolator.propagate(initDate);
417 final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit().toOrbit());
418
419 FieldSpacecraftState<T> BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
420 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
421
422
423
424
425
426
427
428 final double minStep = 0.001;
429 final double maxstep = 1000.0;
430 final double positionTolerance = 10.0;
431 final OrbitType propagationType = OrbitType.KEPLERIAN;
432 final double[][] tolerances =
433 NumericalPropagator.tolerances(positionTolerance, InitOrbit, propagationType);
434 final AdaptiveStepsizeIntegrator integrator =
435 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
436
437
438 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
439 NumPropagator.setOrbitType(propagationType);
440
441 final ForceModel holmesFeatherstone =
442 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
443 NumPropagator.addForceModel(holmesFeatherstone);
444
445
446 NumPropagator.setInitialState(initialOsculatingState.toSpacecraftState());
447
448
449 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift));
450 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
451
452
453
454
455
456 Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.17);
457 Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
458 Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
459 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
460 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.197);
461 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
462 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.00072);
463 Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
464 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12);
465 }
466
467
468 @Test
469 public void compareToNumericalPropagationResetInitialIntermediate() {
470 doCompareToNumericalPropagationResetInitialIntermediate(Decimal64Field.getInstance());
471 }
472 private <T extends CalculusFieldElement<T>> void doCompareToNumericalPropagationResetInitialIntermediate(Field<T> field) {
473
474 T zero = field.getZero();
475 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
476 final Frame inertialFrame = FramesFactory.getEME2000();
477 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
478 double timeshift = 60000. ;
479
480
481 final double a = 24396159;
482 final double e = 0.01;
483 final double i = FastMath.toRadians(7);
484 final double omega = FastMath.toRadians(180);
485 final double raan = FastMath.toRadians(261);
486 final double lM = 0;
487 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
488 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
489 inertialFrame, initDate, zero.add(provider.getMu()));
490
491 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
492
493
494
495
496
497 FieldBrouwerLyddanePropagator<T> BLextrapolator1 =
498 new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, zero.add(Propagator.DEFAULT_MASS), GravityFieldFactory.getUnnormalizedProvider(provider),
499 PropagationType.OSCULATING, BrouwerLyddanePropagator.M2);
500
501
502
503
504 FieldBrouwerLyddanePropagator<T> BLextrapolator2 =
505 new FieldBrouwerLyddanePropagator<>( new FieldKeplerianOrbit<>(zero.add(a + 3000), zero.add(e + 0.001),
506 zero.add(i - FastMath.toRadians(12.0)), zero.add(omega),
507 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
508 inertialFrame, initDate, zero.add(provider.getMu())),DEFAULT_LAW,
509 zero.add(Propagator.DEFAULT_MASS),
510 GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
511
512 BLextrapolator2.resetInitialState(initialState);
513
514 FieldSpacecraftState<T> BLFinalState1 = BLextrapolator1.propagate(initDate.shiftedBy(timeshift));
515 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit().toOrbit());
516 FieldSpacecraftState<T> BLFinalState2 = BLextrapolator2.propagate(initDate.shiftedBy(timeshift));
517 BLextrapolator2.resetIntermediateState(BLFinalState1, true);
518 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit().toOrbit());
519
520 Assert.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
521 Assert.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
522 Assert.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
523 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
524 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
525 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
526 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
527 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
528 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
529
530 }
531
532
533 @Test
534 public void compareConstructors() {
535 doCompareConstructors(Decimal64Field.getInstance());
536 }
537 private <T extends CalculusFieldElement<T>> void doCompareConstructors(Field<T> field) {
538
539 T zero = field.getZero();
540 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
541 final Frame inertialFrame = FramesFactory.getEME2000();
542 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
543 double timeshift = 600. ;
544
545
546 final double a = 24396159;
547 final double e = 0.01;
548 final double i = FastMath.toRadians(7);
549 final double omega = FastMath.toRadians(180);
550 final double raan = FastMath.toRadians(261);
551 final double lM = 0;
552 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
553 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
554 inertialFrame, initDate, zero.add(provider.getMu()));
555
556
557 FieldBrouwerLyddanePropagator<T> BLPropagator1 = new FieldBrouwerLyddanePropagator<T>(initialOrbit, DEFAULT_LAW,
558 provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
559 FieldBrouwerLyddanePropagator<T> BLPropagator2 = new FieldBrouwerLyddanePropagator<>(initialOrbit,
560 provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
561 FieldBrouwerLyddanePropagator<T> BLPropagator3 = new FieldBrouwerLyddanePropagator<>(initialOrbit,
562 zero.add(Propagator.DEFAULT_MASS), provider.getAe(), zero.add(provider.getMu()), -1.08263e-3,
563 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
564
565 FieldSpacecraftState<T> BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift));
566 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit().toOrbit());
567 FieldSpacecraftState<T> BLFinalState2 = BLPropagator2.propagate(initDate.shiftedBy(timeshift));
568 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit().toOrbit());
569 FieldSpacecraftState<T> BLFinalState3 = BLPropagator3.propagate(initDate.shiftedBy(timeshift));
570 final KeplerianOrbit BLOrbit3 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState3.getOrbit().toOrbit());
571
572
573 Assert.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
574 Assert.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
575 Assert.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
576 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
577 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
578 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
579 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
580 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
581 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
582 Assert.assertEquals(BLOrbit1.getA(), BLOrbit3.getA(), 0.0);
583 Assert.assertEquals(BLOrbit1.getE(), BLOrbit3.getE(), 0.0);
584 Assert.assertEquals(BLOrbit1.getI(), BLOrbit3.getI(), 0.0);
585 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
586 MathUtils.normalizeAngle(BLOrbit3.getPerigeeArgument(), FastMath.PI), 0.0);
587 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
588 MathUtils.normalizeAngle(BLOrbit3.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
589 Assert.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
590 MathUtils.normalizeAngle(BLOrbit3.getTrueAnomaly(), FastMath.PI), 0.0);
591
592 }
593
594
595 @Test
596 public void undergroundOrbit() {
597 doUndergroundOrbit(Decimal64Field.getInstance());
598 }
599 private <T extends CalculusFieldElement<T>> void doUndergroundOrbit(Field<T> field) {
600
601 T zero = field.getZero();
602 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
603 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
604
605
606 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
607 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(800.0), zero.add(100.0));
608
609 FieldOrbit<T> initialOrbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<T>(position, velocity),
610 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
611
612
613 try {
614
615 FieldBrouwerLyddanePropagator<T> extrapolator =
616 new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, provider.getAe(), zero.add(provider.getMu()),
617 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
618
619
620
621 double delta_t = 0.0;
622 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
623 extrapolator.propagate(extrapDate);
624
625 } catch (OrekitException oe) {
626 Assert.assertEquals(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE, oe.getSpecifier());
627 }
628 }
629
630
631 @Test
632 public void tooEllipticalOrbit() {
633 doTooEllipticalOrbit(Decimal64Field.getInstance());
634 }
635 private <T extends CalculusFieldElement<T>> void doTooEllipticalOrbit(Field<T> field) {
636
637 T zero = field.getZero();
638 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
639 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
640 FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(67679244.0), zero.add(1.0), zero.add(1.85850),
641 zero.add(2.1), zero.add(2.9), zero.add(6.2), PositionAngle.TRUE,
642 FramesFactory.getEME2000(), initDate, zero.add(provider.getMu()));
643 try {
644
645
646 FieldBrouwerLyddanePropagator<T> extrapolator =
647 new FieldBrouwerLyddanePropagator<>(initialOrbit, provider.getAe(), zero.add(provider.getMu()),
648 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
649
650
651
652 double delta_t = 0.0;
653 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
654 extrapolator.propagate(extrapDate);
655
656 } catch (OrekitException oe) {
657 Assert.assertEquals(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL, oe.getSpecifier());
658 }
659 }
660
661
662 @Test
663 public void criticalInclination() {
664 doCriticalInclination(Decimal64Field.getInstance());
665 }
666 private <T extends CalculusFieldElement<T>> void doCriticalInclination(Field<T> field) {
667
668 final Frame inertialFrame = FramesFactory.getEME2000();
669
670
671 final double a = 24396159;
672 final double e = 0.01;
673 final double i = FastMath.toRadians(7);
674 final double omega = FastMath.toRadians(180);
675 final double raan = FastMath.toRadians(261);
676 final double lM = 0;
677
678 T zero = field.getZero();
679 FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field);
680 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
681 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
682 inertialFrame, initDate, zero.add(provider.getMu()));
683
684
685
686 FieldBrouwerLyddanePropagator<T> extrapolator =
687 new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
688
689
690
691 final FieldSpacecraftState<T> finalOrbit = extrapolator.propagate(initDate);
692
693
694 Assert.assertEquals(0.0,
695 FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
696 finalOrbit.getPVCoordinates().getPosition()).getReal(),
697 7.0e-8);
698
699 Assert.assertEquals(0.0,
700 FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
701 finalOrbit.getPVCoordinates().getVelocity()).getReal(),
702 1.2e-11);
703
704 Assert.assertEquals(0.0, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.0);
705
706 }
707
708 @Test(expected = OrekitException.class)
709 public void testUnableToComputeBLMeanParameters() {
710 doTestUnableToComputeBLMeanParameters(Decimal64Field.getInstance());
711 }
712
713 private <T extends CalculusFieldElement<T>> void doTestUnableToComputeBLMeanParameters(Field<T> field) {
714
715 T zero = field.getZero();
716 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
717 final Frame inertialFrame = FramesFactory.getEME2000();
718 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
719
720
721 final double a = 24396159;
722 final double e = 0.9;
723 final double i = FastMath.toRadians(7);
724 final double omega = FastMath.toRadians(180);
725 final double raan = FastMath.toRadians(261);
726 final double lM = FastMath.toRadians(0);
727 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
728 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
729 inertialFrame, initDate, zero.add(provider.getMu()));
730
731
732
733 final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
734
735
736
737 T delta_t = zero;
738 FieldAbsoluteDate<T> extrapDate = initDate.shiftedBy(delta_t);
739 blField.propagate(extrapDate);
740
741 }
742
743 @Test
744 public void testMeanComparisonWithNonField() {
745 doTestMeanComparisonWithNonField(Decimal64Field.getInstance());
746 }
747
748 private <T extends CalculusFieldElement<T>> void doTestMeanComparisonWithNonField(Field<T> field) {
749
750 T zero = field.getZero();
751 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
752 final Frame inertialFrame = FramesFactory.getEME2000();
753 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
754 double timeshift = -60000.0;
755
756
757 final double a = 24396159;
758 final double e = 0.9;
759 final double i = FastMath.toRadians(7);
760 final double omega = FastMath.toRadians(180);
761 final double raan = FastMath.toRadians(261);
762 final double lM = FastMath.toRadians(0);
763 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
764 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
765 inertialFrame, initDate, zero.add(provider.getMu()));
766
767
768 final FieldSpacecraftState<T> initialStateField = new FieldSpacecraftState<>(initialOrbit);
769 final SpacecraftState initialState = initialStateField.toSpacecraftState();
770
771
772 final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(),
773 GravityFieldFactory.getUnnormalizedProvider(provider),
774 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
775 final FieldSpacecraftState<T> finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift));
776 final FieldKeplerianOrbit<T> finalOrbitField = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit());
777 final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit();
778
779
780 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(),
781 GravityFieldFactory.getUnnormalizedProvider(provider),
782 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
783 final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift));
784 final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
785
786 Assert.assertEquals(finalOrbitFieldReal.getA(), finalOrbit.getA(), Double.MIN_VALUE);
787 Assert.assertEquals(finalOrbitFieldReal.getE(), finalOrbit.getE(), Double.MIN_VALUE);
788 Assert.assertEquals(finalOrbitFieldReal.getI(), finalOrbit.getI(), Double.MIN_VALUE);
789 Assert.assertEquals(finalOrbitFieldReal.getRightAscensionOfAscendingNode(), + finalOrbit.getRightAscensionOfAscendingNode(), Double.MIN_VALUE);
790 Assert.assertEquals(finalOrbitFieldReal.getPerigeeArgument(), finalOrbit.getPerigeeArgument(), Double.MIN_VALUE);
791 Assert.assertEquals(finalOrbitFieldReal.getMeanAnomaly(), finalOrbit.getMeanAnomaly(), Double.MIN_VALUE);
792 Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getPosition().distance(finalOrbit.getPVCoordinates().getPosition()), Double.MIN_VALUE);
793 Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getVelocity().distance(finalOrbit.getPVCoordinates().getVelocity()), Double.MIN_VALUE);
794 }
795
796 @Test
797 public void testOsculatingComparisonWithNonField() {
798 doTestOsculatingComparisonWithNonField(Decimal64Field.getInstance());
799 }
800
801 private <T extends CalculusFieldElement<T>> void doTestOsculatingComparisonWithNonField(Field<T> field) {
802
803 T zero = field.getZero();
804 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
805 final Frame inertialFrame = FramesFactory.getEME2000();
806 FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
807 double timeshift = 60000.0;
808
809
810 final double a = 24396159;
811 final double e = 0.01;
812 final double i = FastMath.toRadians(7);
813 final double omega = FastMath.toRadians(180);
814 final double raan = FastMath.toRadians(261);
815 final double lM = FastMath.toRadians(0);
816 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
817 zero.add(raan), zero.add(lM), PositionAngle.TRUE,
818 inertialFrame, initDate, zero.add(provider.getMu()));
819
820
821 final FieldSpacecraftState<T> initialStateField = new FieldSpacecraftState<>(initialOrbit);
822 final SpacecraftState initialState = initialStateField.toSpacecraftState();
823
824
825 final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(),
826 GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
827 final FieldSpacecraftState<T> finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift));
828 final FieldKeplerianOrbit<T> finalOrbitField = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit());
829 final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit();
830
831
832 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(),
833 GravityFieldFactory.getUnnormalizedProvider(provider),
834 BrouwerLyddanePropagator.M2);
835 final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift));
836 final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
837
838 Assert.assertEquals(finalOrbitFieldReal.getA(), finalOrbit.getA(), Double.MIN_VALUE);
839 Assert.assertEquals(finalOrbitFieldReal.getE(), finalOrbit.getE(), Double.MIN_VALUE);
840 Assert.assertEquals(finalOrbitFieldReal.getI(), finalOrbit.getI(), Double.MIN_VALUE);
841 Assert.assertEquals(finalOrbitFieldReal.getRightAscensionOfAscendingNode(), + finalOrbit.getRightAscensionOfAscendingNode(), Double.MIN_VALUE);
842 Assert.assertEquals(finalOrbitFieldReal.getPerigeeArgument(), finalOrbit.getPerigeeArgument(), Double.MIN_VALUE);
843 Assert.assertEquals(finalOrbitFieldReal.getMeanAnomaly(), finalOrbit.getMeanAnomaly(), Double.MIN_VALUE);
844 Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getPosition().distance(finalOrbit.getPVCoordinates().getPosition()), Double.MIN_VALUE);
845 Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getVelocity().distance(finalOrbit.getPVCoordinates().getVelocity()), Double.MIN_VALUE);
846 }
847
848 @Before
849 public void setUp() {
850 Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
851 provider = GravityFieldFactory.getNormalizedProvider(5, 0);
852 }
853
854 @After
855 public void tearDown() {
856 provider = null;
857 }
858
859 private NormalizedSphericalHarmonicsProvider provider;
860 }