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.util.ArrayList;
21 import java.util.List;
22 import java.util.concurrent.TimeUnit;
23
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
26 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
27 import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
28 import org.hipparchus.stat.descriptive.rank.Max;
29 import org.hipparchus.stat.descriptive.rank.Min;
30 import org.hipparchus.util.FastMath;
31 import org.hipparchus.util.MathUtils;
32 import org.junit.jupiter.api.AfterEach;
33 import org.junit.jupiter.api.Assertions;
34 import org.junit.jupiter.api.BeforeEach;
35 import org.junit.jupiter.api.Test;
36 import org.orekit.Utils;
37 import org.orekit.attitudes.AttitudeProvider;
38 import org.orekit.bodies.CelestialBodyFactory;
39 import org.orekit.bodies.OneAxisEllipsoid;
40 import org.orekit.errors.OrekitException;
41 import org.orekit.errors.OrekitIllegalArgumentException;
42 import org.orekit.forces.ForceModel;
43 import org.orekit.forces.drag.DragForce;
44 import org.orekit.forces.drag.IsotropicDrag;
45 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
46 import org.orekit.forces.gravity.potential.GravityFieldFactory;
47 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
48 import org.orekit.forces.gravity.potential.TideSystem;
49 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
50 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
51 import org.orekit.frames.Frame;
52 import org.orekit.frames.FramesFactory;
53 import org.orekit.models.earth.atmosphere.DTM2000;
54 import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
55 import org.orekit.orbits.CartesianOrbit;
56 import org.orekit.orbits.CircularOrbit;
57 import org.orekit.orbits.EquinoctialOrbit;
58 import org.orekit.orbits.KeplerianOrbit;
59 import org.orekit.orbits.Orbit;
60 import org.orekit.orbits.OrbitType;
61 import org.orekit.orbits.PositionAngleType;
62 import org.orekit.propagation.PropagationType;
63 import org.orekit.propagation.Propagator;
64 import org.orekit.propagation.SpacecraftState;
65 import org.orekit.propagation.ToleranceProvider;
66 import org.orekit.propagation.analytical.tle.TLE;
67 import org.orekit.propagation.analytical.tle.TLEPropagator;
68 import org.orekit.propagation.numerical.NumericalPropagator;
69 import org.orekit.time.AbsoluteDate;
70 import org.orekit.time.TimeScale;
71 import org.orekit.time.TimeScalesFactory;
72 import org.orekit.utils.Constants;
73 import org.orekit.utils.IERSConventions;
74 import org.orekit.utils.PVCoordinates;
75
76 public class BrouwerLyddanePropagatorTest {
77
78 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
79
80 @Test
81 public void testIssue947ComputeMeanOrbit() {
82
83 TLE tleOrbit = new TLE("1 43196U 18015E 21055.59816856 .00000894 00000-0 38966-4 0 9996",
84 "2 43196 97.4662 188.8169 0016935 299.6845 60.2706 15.24746686170319");
85 Propagator propagator = TLEPropagator.selectExtrapolator(tleOrbit);
86
87
88 SpacecraftState tleState = propagator.getInitialState();
89 SpacecraftState tleStateAtDate = propagator.propagate(propagator.getInitialState().getDate().shiftedBy(3, TimeUnit.DAYS));
90
91
92 double epsilon = 1.0e-12;
93 int maxIter = 500;
94 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 0);
95
96 Assertions.assertDoesNotThrow(() -> {
97 BrouwerLyddanePropagator.computeMeanOrbit(tleState.getOrbit(), provider,
98 provider.onDate(tleState.getDate()),
99 BrouwerLyddanePropagator.M2,
100 epsilon, maxIter);
101 });
102
103 Assertions.assertDoesNotThrow(() -> {
104 BrouwerLyddanePropagator.computeMeanOrbit(tleStateAtDate.getOrbit(), provider,
105 provider.onDate(tleStateAtDate.getDate()),
106 BrouwerLyddanePropagator.M2,
107 epsilon, maxIter);
108 });
109 }
110
111 @Test
112 public void testIssue947Propagate() {
113
114 final double x1 = -1772.619869591273 * 1000;
115 final double y1 = -3908.6521138424428 * 1000;
116 final double z1 = 5266.68093513367 * 1000;
117 final double dx1 = 6.35969327821623 * 1000;
118 final double dy1 = -4.165238186695803 * 1000;
119 final double dz1 = -0.9458311825913897 * 1000;
120
121 final Vector3D pos = new Vector3D(x1, y1, z1);
122 final Vector3D vel = new Vector3D(dx1, dy1, dz1);
123 final TimeScale tai = TimeScalesFactory.getTAI();
124 final AbsoluteDate date = new AbsoluteDate(2023, 3, 1, 14, 6, 29.639584, tai);
125 final CircularOrbit orbit = new CircularOrbit(new PVCoordinates(pos, vel),
126 FramesFactory.getGCRF(), date,
127 Constants.EGM96_EARTH_MU);
128
129 final BrouwerLyddanePropagator blPropagator = new BrouwerLyddanePropagator(orbit, 1000.,
130 Constants.IERS2010_EARTH_EQUATORIAL_RADIUS, Constants.IERS2010_EARTH_MU,
131 Constants.EIGEN5C_EARTH_C20, Constants.EIGEN5C_EARTH_C30,
132 Constants.EIGEN5C_EARTH_C40, Constants.EIGEN5C_EARTH_C50,
133 BrouwerLyddanePropagator.M2);
134 final List<SpacecraftState> states = new ArrayList<>();
135 blPropagator.setStepHandler(10, states::add);
136 final AbsoluteDate startDate = new AbsoluteDate(2023, 3, 1, 14, 6, 29.639584, tai);
137 final AbsoluteDate endDate = new AbsoluteDate(2023, 3, 1, 15, 6, 29.639584, tai);
138 Assertions.assertDoesNotThrow(() -> blPropagator.propagate(startDate, endDate));
139 }
140
141 @Test
142 public void testIssue947Constructor() {
143
144 final Frame eme2000 = FramesFactory.getEME2000();
145
146 final AbsoluteDate t = new AbsoluteDate(2023,06,26,12,0,0.0,TimeScalesFactory.getUTC());
147 final double x = 6313554.48504233;
148 final double y = 2775620.8433687;
149 final double z = -2111774.8221765;
150 final double vx = 1575.31305905025;
151 final double vy = 1786.43351611741;
152 final double vz = 7042.05214468662;
153
154 final PVCoordinates pv = new PVCoordinates(new Vector3D(x, y, z), new Vector3D(vx, vy, vz));
155 final CartesianOrbit orb = new CartesianOrbit(pv, eme2000, t, Constants.EGM96_EARTH_MU);
156
157 Assertions.assertDoesNotThrow(() -> {
158 new BrouwerLyddanePropagator(orb, Constants.EGM96_EARTH_EQUATORIAL_RADIUS,
159 Constants.EGM96_EARTH_MU,
160 Constants.EGM96_EARTH_C20, Constants.EGM96_EARTH_C30,
161 Constants.EGM96_EARTH_C40, Constants.EGM96_EARTH_C50,
162 BrouwerLyddanePropagator.M2);
163 });
164 }
165
166 @Test
167 public void sameDateCartesian() {
168
169
170
171
172
173 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
174 Vector3D position = new Vector3D(3220103., 69623., 6149822.);
175 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
176
177 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
178 FramesFactory.getEME2000(), initDate, provider.getMu());
179
180
181
182 BrouwerLyddanePropagator extrapolator =
183 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
184 SpacecraftState finalOrbit = extrapolator.propagate(initDate);
185
186
187 Assertions.assertEquals(0.0,
188 Vector3D.distance(initialOrbit.getPosition(),
189 finalOrbit.getPosition()),
190 4.7e-7);
191 Assertions.assertEquals(0.0,
192 Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
193 finalOrbit.getPVCoordinates().getVelocity()),
194 2.8e-10);
195 Assertions.assertEquals(0.0, finalOrbit.getOrbit().getA() - initialOrbit.getA(), 1.3e-8);
196
197 }
198
199 @Test
200 public void sameDateKeplerian() {
201
202
203
204 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
205 Orbit initialOrbit = new KeplerianOrbit(6767924.41, .0005, 1.7, 2.1, 2.9,
206 6.2, PositionAngleType.TRUE,
207 FramesFactory.getEME2000(), initDate, provider.getMu());
208
209 BrouwerLyddanePropagator extrapolator =
210 new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
211
212 SpacecraftState finalOrbit = extrapolator.propagate(initDate);
213
214
215 Assertions.assertEquals(0.0,
216 Vector3D.distance(initialOrbit.getPosition(),
217 finalOrbit.getPosition()),
218 4.0e-7);
219 Assertions.assertEquals(0.0,
220 Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
221 finalOrbit.getPVCoordinates().getVelocity()),
222 2.9e-10);
223 Assertions.assertEquals(0.0, finalOrbit.getOrbit().getA() - initialOrbit.getA(), 3.8e-9);
224 }
225
226
227 @Test
228 public void almostSphericalBody() {
229
230
231
232
233 Vector3D position = new Vector3D(3220103., 69623., 8449822.);
234 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
235
236 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
237 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
238 FramesFactory.getEME2000(), initDate, provider.getMu());
239
240
241
242
243
244
245
246 UnnormalizedSphericalHarmonicsProvider kepProvider =
247 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
248 TideSystem.UNKNOWN,
249 new double[][] {
250 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
251 }, new double[][] {
252 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
253 });
254
255
256
257 BrouwerLyddanePropagator extrapolatorAna =
258 new BrouwerLyddanePropagator(initialOrbit, 1000.0, kepProvider, BrouwerLyddanePropagator.M2);
259 KeplerianPropagator extrapolatorKep = new KeplerianPropagator(initialOrbit);
260
261
262
263 double delta_t = 100.0;
264 AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
265
266 Orbit finalOrbitAna = extrapolatorAna.propagate(extrapDate).getOrbit();
267 Orbit finalOrbitKep = extrapolatorKep.propagate(extrapDate).getOrbit();
268
269 Assertions.assertEquals(0.0, finalOrbitAna.getDate().durationFrom(extrapDate), 0.0);
270
271 Assertions.assertEquals(finalOrbitAna.getA(), finalOrbitKep.getA(), 10
272 * Utils.epsilonTest * finalOrbitKep.getA());
273 Assertions.assertEquals(finalOrbitAna.getEquinoctialEx(), finalOrbitKep.getEquinoctialEx(), Utils.epsilonE
274 * finalOrbitKep.getE());
275 Assertions.assertEquals(finalOrbitAna.getEquinoctialEy(), finalOrbitKep.getEquinoctialEy(), Utils.epsilonE
276 * finalOrbitKep.getE());
277 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx(), finalOrbitKep.getHx()),
278 finalOrbitKep.getHx(), Utils.epsilonAngle
279 * FastMath.abs(finalOrbitKep.getI()));
280 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy(), finalOrbitKep.getHy()),
281 finalOrbitKep.getHy(), Utils.epsilonAngle
282 * FastMath.abs(finalOrbitKep.getI()));
283 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv(), finalOrbitKep.getLv()),
284 finalOrbitKep.getLv(), Utils.epsilonAngle
285 * FastMath.abs(finalOrbitKep.getLv()));
286 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE(), finalOrbitKep.getLE()),
287 finalOrbitKep.getLE(), Utils.epsilonAngle
288 * FastMath.abs(finalOrbitKep.getLE()));
289 Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM(), finalOrbitKep.getLM()),
290 finalOrbitKep.getLM(), Utils.epsilonAngle
291 * FastMath.abs(finalOrbitKep.getLM()));
292 }
293
294 @Test
295 public void compareToNumericalPropagation() {
296
297 final Frame inertialFrame = FramesFactory.getEME2000();
298 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
299 double timeshift = 60000. ;
300
301
302 final double a = 24396159;
303 final double e = 0.01;
304 final double i = FastMath.toRadians(7);
305 final double omega = FastMath.toRadians(180);
306 final double raan = FastMath.toRadians(261);
307 final double lM = 0;
308 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.TRUE,
309 inertialFrame, initDate, provider.getMu());
310
311 final SpacecraftState initialState = new SpacecraftState(initialOrbit);
312
313
314
315
316
317
318 final double minStep = 0.001;
319 final double maxstep = 1000.0;
320 final double positionTolerance = 10.0;
321 final OrbitType propagationType = OrbitType.KEPLERIAN;
322 final double[][] tolerances =
323 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit, propagationType);
324 final AdaptiveStepsizeIntegrator integrator =
325 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
326
327
328 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
329 NumPropagator.setOrbitType(propagationType);
330
331 final ForceModel holmesFeatherstone =
332 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
333 NumPropagator.addForceModel(holmesFeatherstone);
334
335
336 NumPropagator.setInitialState(initialState);
337
338
339 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
340 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
341
342
343
344
345
346 BrouwerLyddanePropagator BLextrapolator =
347 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
348
349 SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
350 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
351
352 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.175);
353 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 3.2e-6);
354 Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 6.9e-8);
355 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
356 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0053);
357 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
358 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 1.2e-6);
359 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
360 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0052);
361 }
362
363 @Test
364 public void compareToNumericalPropagationWithDrag() {
365
366 final Frame inertialFrame = FramesFactory.getEME2000();
367 final TimeScale utc = TimeScalesFactory.getUTC();
368 final AbsoluteDate initDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
369 double timeshift = 60000. ;
370
371
372 final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3;
373 final double e = 0.01;
374 final double i = FastMath.toRadians(7);
375 final double omega = FastMath.toRadians(180);
376 final double raan = FastMath.toRadians(261);
377 final double lM = 0;
378 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.TRUE,
379 inertialFrame, initDate, provider.getMu());
380
381 final SpacecraftState initialState = new SpacecraftState(initialOrbit);
382
383
384
385
386
387
388 final double minStep = 0.001;
389 final double maxstep = 1000.0;
390 final double positionTolerance = 10.0;
391 final OrbitType propagationType = OrbitType.KEPLERIAN;
392 final double[][] tolerances =
393 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit, propagationType);
394 final AdaptiveStepsizeIntegrator integrator =
395 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
396
397
398 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
399 NumPropagator.setOrbitType(propagationType);
400
401
402 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
403 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
404 MarshallSolarActivityFutureEstimation msafe =
405 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
406 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
407 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
408
409
410 final ForceModel holmesFeatherstone =
411 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
412 final ForceModel drag = new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0));
413 NumPropagator.addForceModel(holmesFeatherstone);
414 NumPropagator.addForceModel(drag);
415
416
417 NumPropagator.setInitialState(initialState);
418
419
420 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
421 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
422
423
424
425
426
427 BrouwerLyddanePropagator BLextrapolator =
428 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
429
430 SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
431 KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
432
433
434 final double deltaSmaBefore = 15.81;
435 final double deltaEccBefore = 9.2681e-5;
436 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore);
437 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore);
438
439
440
441
442
443 double M2 = 1.0e-14;
444 BLextrapolator = new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2);
445 BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
446 BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
447
448
449 final double deltaSmaAfter = 11.04;
450 final double deltaEccAfter = 9.2639e-5;
451 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter);
452 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter);
453
454 Assertions.assertTrue(deltaSmaAfter < deltaSmaBefore);
455 Assertions.assertTrue(deltaEccAfter < deltaEccBefore);
456 }
457
458
459 @Test
460 public void compareToNumericalPropagationMeanInitialOrbit() {
461
462 final Frame inertialFrame = FramesFactory.getEME2000();
463 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
464 double timeshift = 60000. ;
465
466
467 final double a = 24396159;
468 final double e = 0.01;
469 final double i = FastMath.toRadians(7);
470 final double omega = FastMath.toRadians(180);
471 final double raan = FastMath.toRadians(261);
472 final double lM = 0;
473 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.TRUE,
474 inertialFrame, initDate, provider.getMu());
475
476 BrouwerLyddanePropagator BLextrapolator =
477 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider),
478 PropagationType.MEAN, BrouwerLyddanePropagator.M2);
479 SpacecraftState initialOsculatingState = BLextrapolator.propagate(initDate);
480 final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit());
481
482 SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
483 final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
484
485
486
487
488
489
490 final double minStep = 0.001;
491 final double maxstep = 1000.0;
492 final double positionTolerance = 10.0;
493 final OrbitType propagationType = OrbitType.KEPLERIAN;
494 final double[][] tolerances =
495 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(InitOrbit, propagationType);
496 final AdaptiveStepsizeIntegrator integrator =
497 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
498
499
500 final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
501 NumPropagator.setOrbitType(propagationType);
502
503 final ForceModel holmesFeatherstone =
504 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
505 NumPropagator.addForceModel(holmesFeatherstone);
506
507
508 NumPropagator.setInitialState(initialOsculatingState);
509
510
511 final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
512 final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
513
514
515 Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.174);
516 Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 3.2e-6);
517 Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 6.9e-8);
518 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
519 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0053);
520 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
521 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 1.2e-6);
522 Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
523 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0052);
524 }
525
526 @Test
527 public void compareToNumericalPropagationResetInitialIntermediate() {
528
529 final Frame inertialFrame = FramesFactory.getEME2000();
530 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
531 double timeshift = 60000.;
532
533
534 final double a = 24396159;
535 final double e = 0.01;
536 final double i = FastMath.toRadians(7);
537 final double omega = FastMath.toRadians(180);
538 final double raan = FastMath.toRadians(261);
539 final double lM = 0;
540 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.TRUE,
541 inertialFrame, initDate, provider.getMu());
542
543 final SpacecraftState initialState = new SpacecraftState(initialOrbit);
544
545
546
547
548
549 BrouwerLyddanePropagator BLextrapolator1 =
550 new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider),
551 PropagationType.OSCULATING, BrouwerLyddanePropagator.M2);
552
553
554
555
556 BrouwerLyddanePropagator BLextrapolator2 =
557 new BrouwerLyddanePropagator( new KeplerianOrbit(a + 3000, e + 0.001, i - FastMath.toRadians(12.0), omega, raan, lM, PositionAngleType.TRUE,
558 inertialFrame, initDate, provider.getMu()),DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
559
560 BLextrapolator2.resetInitialState(initialState);
561
562 SpacecraftState BLFinalState1 = BLextrapolator1.propagate(initDate.shiftedBy(timeshift));
563 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit());
564 SpacecraftState BLFinalState2 = BLextrapolator2.propagate(initDate.shiftedBy(timeshift));
565 BLextrapolator2.resetIntermediateState(BLFinalState1, true);
566 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit());
567
568
569 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
570 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
571 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
572 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
573 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
574 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
575 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
576 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
577 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
578 }
579
580 @Test
581 public void compareConstructors() {
582
583 final Frame inertialFrame = FramesFactory.getEME2000();
584 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
585 double timeshift = 600. ;
586
587
588 final double a = 24396159;
589 final double e = 0.01;
590 final double i = FastMath.toRadians(7);
591 final double omega = FastMath.toRadians(180);
592 final double raan = FastMath.toRadians(261);
593 final double lV = 0;
594 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lV, PositionAngleType.TRUE,
595 inertialFrame, initDate, provider.getMu());
596
597 BrouwerLyddanePropagator BLPropagator1 = new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW,
598 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
599 BrouwerLyddanePropagator BLPropagator2 = new BrouwerLyddanePropagator(initialOrbit,
600 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
601 BrouwerLyddanePropagator BLPropagator3 = new BrouwerLyddanePropagator(initialOrbit, Propagator.DEFAULT_MASS,
602 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
603
604 SpacecraftState BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift));
605 final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit());
606 SpacecraftState BLFinalState2 = BLPropagator2.propagate(initDate.shiftedBy(timeshift));
607 final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit());
608 SpacecraftState BLFinalState3 = BLPropagator3.propagate(initDate.shiftedBy(timeshift));
609 final KeplerianOrbit BLOrbit3 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState3.getOrbit());
610
611
612 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
613 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
614 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
615 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
616 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
617 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
618 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
619 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
620 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
621
622 Assertions.assertEquals(BLOrbit1.getA(), BLOrbit3.getA(), 0.0);
623 Assertions.assertEquals(BLOrbit1.getE(), BLOrbit3.getE(), 0.0);
624 Assertions.assertEquals(BLOrbit1.getI(), BLOrbit3.getI(), 0.0);
625 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
626 MathUtils.normalizeAngle(BLOrbit3.getPerigeeArgument(), FastMath.PI), 0.0);
627 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
628 MathUtils.normalizeAngle(BLOrbit3.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
629 Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
630 MathUtils.normalizeAngle(BLOrbit3.getTrueAnomaly(), FastMath.PI), 0.0);
631 }
632
633 @Test
634 public void undergroundOrbit() {
635
636 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
637 Vector3D velocity = new Vector3D(-500.0, 800.0, 100.0);
638 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
639 Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
640 FramesFactory.getEME2000(), initDate, provider.getMu());
641
642
643 OrekitException oe = Assertions.assertThrows(OrekitException.class, () -> {
644 new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, provider.getAe(), provider.getMu(),
645 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
646 });
647 Assertions.assertTrue(oe.getMessage().contains("trajectory inside the Brillouin sphere (r ="));
648 }
649
650
651 @Test
652 public void tooEllipticalOrbit() {
653
654 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
655 Orbit initialOrbit = new KeplerianOrbit(67679244.0, 1.0, 1.85850, 2.1, 2.9,
656 6.2, PositionAngleType.TRUE, FramesFactory.getEME2000(),
657 initDate, provider.getMu());
658
659
660 OrekitException oe = Assertions.assertThrows(OrekitException.class, () -> {
661 new BrouwerLyddanePropagator(initialOrbit, provider.getAe(), provider.getMu(),
662 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
663 });
664 Assertions.assertTrue(oe.getMessage().contains("too large eccentricity for propagation model: e ="));
665 }
666
667
668 @Test
669 public void criticalInclination() {
670
671 final Frame inertialFrame = FramesFactory.getEME2000();
672
673 final double a = 24396159;
674 final double e = 0.01;
675 final double i = FastMath.acos(1.0 / FastMath.sqrt(5.0));
676 final double omega = FastMath.toRadians(180);
677 final double raan = FastMath.toRadians(261);
678 final double lV = 0;
679
680 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
681 final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lV, PositionAngleType.TRUE,
682 inertialFrame, initDate, provider.getMu());
683
684
685
686 BrouwerLyddanePropagator extrapolator =
687 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
688
689
690
691 SpacecraftState finalOrbit = extrapolator.propagate(initDate);
692
693
694
695 Assertions.assertEquals(0.0,
696 Vector3D.distance(initialOrbit.getPosition(),
697 finalOrbit.getPosition()),
698 1.5e-8);
699 Assertions.assertEquals(0.0,
700 Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
701 finalOrbit.getPVCoordinates().getVelocity()),
702 2.7e-12);
703 Assertions.assertEquals(0.0, finalOrbit.getOrbit().getA() - initialOrbit.getA(), 7.5e-9);
704 }
705
706 @Test
707 public void insideEarth() {
708
709 final AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
710 final Orbit initialOrbit = new KeplerianOrbit(provider.getAe()-1000.0, 0.01,
711 FastMath.toRadians(10.0), 2.1, 2.9,
712 6.2, PositionAngleType.TRUE,
713 FramesFactory.getEME2000(),
714 initDate, provider.getMu());
715
716
717 OrekitIllegalArgumentException oe = Assertions.assertThrows(OrekitIllegalArgumentException.class, () -> {
718 new BrouwerLyddanePropagator(initialOrbit, Propagator.DEFAULT_MASS, provider.getAe(), provider.getMu(),
719 -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7);
720 });
721 Assertions.assertTrue(oe.getMessage().contains("orbit should be either elliptic with a > 0 and e < 1 or hyperbolic with a < 0 and e > 1"));
722 }
723
724 @Test
725 public void testUnableToComputeBLMeanParameters() {
726 final Frame eci = FramesFactory.getEME2000();
727 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
728
729
730 final double a = 24396159;
731 final double e = 0.9;
732 final double i = FastMath.toRadians(7);
733 final double pa = FastMath.toRadians(180);
734 final double raan = FastMath.toRadians(261);
735 final double lM = FastMath.toRadians(0);
736 final Orbit orbit = new KeplerianOrbit(a, e, i, pa, raan, lM, PositionAngleType.MEAN,
737 eci, date, provider.getMu());
738
739
740 final UnnormalizedSphericalHarmonicsProvider up = GravityFieldFactory.getUnnormalizedProvider(provider);
741 final UnnormalizedSphericalHarmonics uh = up.onDate(date);
742 final double eps = 1.e-13;
743 final int maxIter = 10;
744 OrekitException oe = Assertions.assertThrows(OrekitException.class, () -> {
745 BrouwerLyddanePropagator.computeMeanOrbit(orbit, up, uh, BrouwerLyddanePropagator.M2,
746 eps, maxIter);
747 });
748 Assertions.assertTrue(oe.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
749 }
750
751 @Test
752 public void testMeanOrbit() {
753 final KeplerianOrbit initialOsculating =
754 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
755 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
756 provider.getMu());
757 final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
758 final UnnormalizedSphericalHarmonics ush = ushp.onDate(initialOsculating.getDate());
759
760
761
762 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.1).getTolerances(initialOsculating, OrbitType.KEPLERIAN);
763 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
764 integrator.setInitialStepSize(60);
765 NumericalPropagator num = new NumericalPropagator(integrator);
766 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
767 num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
768 num.setInitialState(new SpacecraftState(initialOsculating));
769 num.setOrbitType(OrbitType.KEPLERIAN);
770 num.setPositionAngleType(initialOsculating.getCachedPositionAngleType());
771 final StorelessUnivariateStatistic oscMin = new Min();
772 final StorelessUnivariateStatistic oscMax = new Max();
773 final StorelessUnivariateStatistic meanMin = new Min();
774 final StorelessUnivariateStatistic meanMax = new Max();
775 num.getMultiplexer().add(60, state -> {
776 final Orbit osc = state.getOrbit();
777 oscMin.increment(osc.getA());
778 oscMax.increment(osc.getA());
779
780 final Orbit mean = BrouwerLyddanePropagator.computeMeanOrbit(state.getOrbit(), ushp, ush, BrouwerLyddanePropagator.M2);
781 meanMin.increment(mean.getA());
782 meanMax.increment(mean.getA());
783 });
784 num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
785
786
787 Assertions.assertEquals(3188.347, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
788 Assertions.assertEquals( 25.794, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
789
790 }
791
792 @Test
793 public void testGeostationaryOrbit() {
794
795 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
796 final Frame eci = FramesFactory.getEME2000();
797 final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
798
799
800 final double sma = FastMath.cbrt(Constants.IERS2010_EARTH_MU /
801 Constants.IERS2010_EARTH_ANGULAR_VELOCITY /
802 Constants.IERS2010_EARTH_ANGULAR_VELOCITY);
803 final double ecc = 0.0;
804 final double inc = 0.0;
805 final double pa = 0.0;
806 final double raan = 0.0;
807 final double lV = 0.0;
808 final Orbit orbit = new KeplerianOrbit(sma, ecc, inc, pa, raan, lV,
809 PositionAngleType.TRUE,
810 eci, date, provider.getMu());
811
812
813 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(orbit, ushp, PropagationType.MEAN,
814 BrouwerLyddanePropagator.M2);
815
816
817 final SpacecraftState state = bl.propagate(date.shiftedBy(Constants.JULIAN_DAY));
818 final KeplerianOrbit orbOsc = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit());
819
820
821 Assertions.assertTrue(Double.isFinite(orbOsc.getA()));
822 Assertions.assertTrue(Double.isFinite(orbOsc.getE()));
823 Assertions.assertTrue(Double.isFinite(orbOsc.getI()));
824 Assertions.assertTrue(Double.isFinite(orbOsc.getPerigeeArgument()));
825 Assertions.assertTrue(Double.isFinite(orbOsc.getRightAscensionOfAscendingNode()));
826 Assertions.assertTrue(Double.isFinite(orbOsc.getTrueAnomaly()));
827
828
829 final BrouwerLyddanePropagator bl2 = new BrouwerLyddanePropagator(orbit, ushp, PropagationType.OSCULATING,
830 BrouwerLyddanePropagator.M2);
831
832
833 final SpacecraftState state2 = bl2.propagate(date.shiftedBy(Constants.JULIAN_DAY));
834 final KeplerianOrbit orbOsc2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state2.getOrbit());
835
836
837 Assertions.assertTrue(Double.isFinite(orbOsc2.getA()));
838 Assertions.assertTrue(Double.isFinite(orbOsc2.getE()));
839 Assertions.assertTrue(Double.isFinite(orbOsc2.getI()));
840 Assertions.assertTrue(Double.isFinite(orbOsc2.getPerigeeArgument()));
841 Assertions.assertTrue(Double.isFinite(orbOsc2.getRightAscensionOfAscendingNode()));
842 Assertions.assertTrue(Double.isFinite(orbOsc2.getTrueAnomaly()));
843 }
844
845 @BeforeEach
846 public void setUp() {
847 Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
848 provider = GravityFieldFactory.getNormalizedProvider(5, 0);
849 }
850
851 @AfterEach
852 public void tearDown() {
853 provider = null;
854 }
855
856 private NormalizedSphericalHarmonicsProvider provider;
857
858 }