1   package org.orekit.propagation.analytical;
2   
3   import org.hipparchus.geometry.euclidean.threed.Vector3D;
4   import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
5   import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
6   import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
7   import org.hipparchus.stat.descriptive.rank.Max;
8   import org.hipparchus.stat.descriptive.rank.Min;
9   import org.hipparchus.util.FastMath;
10  import org.hipparchus.util.MathUtils;
11  import org.junit.jupiter.api.AfterEach;
12  import org.junit.jupiter.api.Assertions;
13  import org.junit.jupiter.api.BeforeEach;
14  import org.junit.jupiter.api.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.forces.ForceModel;
22  import org.orekit.forces.drag.DragForce;
23  import org.orekit.forces.drag.IsotropicDrag;
24  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
25  import org.orekit.forces.gravity.potential.GravityFieldFactory;
26  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
27  import org.orekit.forces.gravity.potential.TideSystem;
28  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
29  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
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.EquinoctialOrbit;
35  import org.orekit.orbits.KeplerianOrbit;
36  import org.orekit.orbits.Orbit;
37  import org.orekit.orbits.OrbitType;
38  import org.orekit.orbits.PositionAngle;
39  import org.orekit.propagation.PropagationType;
40  import org.orekit.propagation.Propagator;
41  import org.orekit.propagation.SpacecraftState;
42  import org.orekit.propagation.numerical.NumericalPropagator;
43  import org.orekit.time.AbsoluteDate;
44  import org.orekit.time.TimeScale;
45  import org.orekit.time.TimeScalesFactory;
46  import org.orekit.utils.Constants;
47  import org.orekit.utils.IERSConventions;
48  import org.orekit.utils.PVCoordinates;
49  
50  public class BrouwerLyddanePropagatorTest {
51  
52      private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
53  
54      @Test
55      public void sameDateCartesian() {
56  
57          // Definition of initial conditions with position and velocity
58          // ------------------------------------------------------------
59          // e = 0.04152500499523033   and   i = 1.705015527659039
60  
61          AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
62          Vector3D position = new Vector3D(3220103., 69623., 6149822.);
63          Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
64  
65          Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
66                                                    FramesFactory.getEME2000(), initDate, provider.getMu());
67  
68          // Extrapolation at the initial date
69          // ---------------------------------
70          BrouwerLyddanePropagator extrapolator =
71                  new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
72          SpacecraftState finalOrbit = extrapolator.propagate(initDate);
73  
74          // positions  velocity and semi major axis match perfectly
75          Assertions.assertEquals(0.0,
76                              Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
77                                                finalOrbit.getPVCoordinates().getPosition()),
78                              3.6e-9);
79  
80          Assertions.assertEquals(0.0,
81                              Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
82                                                finalOrbit.getPVCoordinates().getVelocity()),
83                              3.8e-12);
84          Assertions.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0);
85  
86      }
87  
88      @Test
89      public void sameDateKeplerian() {
90  
91          // Definition of initial conditions with position and velocity
92          // ------------------------------------------------------------
93          AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
94          Orbit initialOrbit = new KeplerianOrbit(6767924.41, .0005,  1.7, 2.1, 2.9,
95                  6.2, PositionAngle.TRUE,
96                  FramesFactory.getEME2000(), initDate, provider.getMu());
97  
98          BrouwerLyddanePropagator extrapolator =
99                  new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
100 
101         SpacecraftState finalOrbit = extrapolator.propagate(initDate);
102 
103         // positions  velocity and semi major axis match perfectly
104         Assertions.assertEquals(0.0,
105                             Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
106                                               finalOrbit.getPVCoordinates().getPosition()),
107                             7.0e-9);
108 
109         Assertions.assertEquals(0.0,
110                             Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
111                                               finalOrbit.getPVCoordinates().getVelocity()),
112                             7.0e-12);
113         Assertions.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0);
114     }
115 
116 
117     @Test
118     public void almostSphericalBody() {
119 
120         // Definition of initial conditions
121         // ---------------------------------
122         // with e around e = 1.4e-4 and i = 1.7 rad
123         Vector3D position = new Vector3D(3220103., 69623., 8449822.);
124         Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
125 
126         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
127         Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
128                                                   FramesFactory.getEME2000(), initDate, provider.getMu());
129 
130         // Initialisation to simulate a Keplerian extrapolation
131         // To be noticed: in order to simulate a Keplerian extrapolation with the
132         // analytical
133         // extrapolator, one should put the zonal coefficients to 0. But due to
134         // numerical pbs
135         // one must put a non 0 value.
136         UnnormalizedSphericalHarmonicsProvider kepProvider =
137                 GravityFieldFactory.getUnnormalizedProvider(6.378137e6, 3.9860047e14,
138                                                             TideSystem.UNKNOWN,
139                                                             new double[][] {
140                                                                 { 0 }, { 0 }, { 0.1e-10 }, { 0.1e-13 }, { 0.1e-13 }, { 0.1e-14 }, { 0.1e-14 }
141                                                             }, new double[][] {
142                                                                 { 0 }, { 0 },  { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
143                                                             });
144 
145         // Extrapolators definitions
146         // -------------------------
147         BrouwerLyddanePropagator extrapolatorAna =
148             new BrouwerLyddanePropagator(initialOrbit, 1000.0, kepProvider, BrouwerLyddanePropagator.M2);
149         KeplerianPropagator extrapolatorKep = new KeplerianPropagator(initialOrbit);
150 
151         // Extrapolation at a final date different from initial date
152         // ---------------------------------------------------------
153         double delta_t = 100.0; // extrapolation duration in seconds
154         AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
155 
156         SpacecraftState finalOrbitAna = extrapolatorAna.propagate(extrapDate);
157         SpacecraftState finalOrbitKep = extrapolatorKep.propagate(extrapDate);
158 
159         Assertions.assertEquals(finalOrbitAna.getDate().durationFrom(extrapDate), 0.0,
160                      Utils.epsilonTest);
161         // comparison of each orbital parameters
162         Assertions.assertEquals(finalOrbitAna.getA(), finalOrbitKep.getA(), 10
163                      * Utils.epsilonTest * finalOrbitKep.getA());
164         Assertions.assertEquals(finalOrbitAna.getEquinoctialEx(), finalOrbitKep.getEquinoctialEx(), Utils.epsilonE
165                      * finalOrbitKep.getE());
166         Assertions.assertEquals(finalOrbitAna.getEquinoctialEy(), finalOrbitKep.getEquinoctialEy(), Utils.epsilonE
167                      * finalOrbitKep.getE());
168         Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHx(), finalOrbitKep.getHx()),
169                      finalOrbitKep.getHx(), Utils.epsilonAngle
170                      * FastMath.abs(finalOrbitKep.getI()));
171         Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getHy(), finalOrbitKep.getHy()),
172                      finalOrbitKep.getHy(), Utils.epsilonAngle
173                      * FastMath.abs(finalOrbitKep.getI()));
174         Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLv(), finalOrbitKep.getLv()),
175                      finalOrbitKep.getLv(), Utils.epsilonAngle
176                      * FastMath.abs(finalOrbitKep.getLv()));
177         Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLE(), finalOrbitKep.getLE()),
178                      finalOrbitKep.getLE(), Utils.epsilonAngle
179                      * FastMath.abs(finalOrbitKep.getLE()));
180         Assertions.assertEquals(MathUtils.normalizeAngle(finalOrbitAna.getLM(), finalOrbitKep.getLM()),
181                      finalOrbitKep.getLM(), Utils.epsilonAngle
182                      * FastMath.abs(finalOrbitKep.getLM()));
183 
184     }
185 
186 
187 
188 
189     @Test
190     public void compareToNumericalPropagation() {
191 
192         final Frame inertialFrame = FramesFactory.getEME2000();
193         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
194         double timeshift = 60000. ;
195 
196         // Initial orbit
197         final double a = 24396159; // semi major axis in meters
198         final double e = 0.01; // eccentricity
199         final double i = FastMath.toRadians(7); // inclination
200         final double omega = FastMath.toRadians(180); // perigee argument
201         final double raan = FastMath.toRadians(261); // right ascention of ascending node
202         final double lM = 0; // mean anomaly
203         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
204                                                       inertialFrame, initDate, provider.getMu());
205         // Initial state definition
206         final SpacecraftState initialState = new SpacecraftState(initialOrbit);
207 
208         //_______________________________________________________________________________________________
209         // SET UP A REFERENCE NUMERICAL PROPAGATION
210         //_______________________________________________________________________________________________
211 
212         // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
213         final double minStep = 0.001;
214         final double maxstep = 1000.0;
215         final double positionTolerance = 10.0;
216         final OrbitType propagationType = OrbitType.KEPLERIAN;
217         final double[][] tolerances =
218                 NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
219         final AdaptiveStepsizeIntegrator integrator =
220                 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
221 
222         // Numerical Propagator
223         final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
224         NumPropagator.setOrbitType(propagationType);
225 
226         final ForceModel holmesFeatherstone =
227                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
228         NumPropagator.addForceModel(holmesFeatherstone);
229 
230         // Set up initial state in the propagator
231         NumPropagator.setInitialState(initialState);
232 
233         // Extrapolate from the initial to the final date
234         final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
235         final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
236 
237         //_______________________________________________________________________________________________
238         // SET UP A BROUWER LYDDANE PROPAGATION
239         //_______________________________________________________________________________________________
240 
241         BrouwerLyddanePropagator BLextrapolator =
242                 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
243 
244         SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
245         final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
246 
247         Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.070);
248         Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
249         Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000000053);
250         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
251                 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0021);
252         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
253                 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0000013);
254         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
255                 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0021);
256     }
257 
258     @Test
259     public void compareToNumericalPropagationWithDrag() {
260 
261         final Frame inertialFrame = FramesFactory.getEME2000();
262         final TimeScale utc = TimeScalesFactory.getUTC();
263         final AbsoluteDate initDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
264         double timeshift = 60000. ;
265 
266         // Initial orbit
267         final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3; // semi major axis in meters
268         final double e = 0.01; // eccentricity
269         final double i = FastMath.toRadians(7); // inclination
270         final double omega = FastMath.toRadians(180); // perigee argument
271         final double raan = FastMath.toRadians(261); // right ascention of ascending node
272         final double lM = 0; // mean anomaly
273         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
274                                                       inertialFrame, initDate, provider.getMu());
275         // Initial state definition
276         final SpacecraftState initialState = new SpacecraftState(initialOrbit);
277 
278         //_______________________________________________________________________________________________
279         // SET UP A REFERENCE NUMERICAL PROPAGATION
280         //_______________________________________________________________________________________________
281 
282         // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
283         final double minStep = 0.001;
284         final double maxstep = 1000.0;
285         final double positionTolerance = 10.0;
286         final OrbitType propagationType = OrbitType.KEPLERIAN;
287         final double[][] tolerances =
288                 NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
289         final AdaptiveStepsizeIntegrator integrator =
290                 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
291 
292         // Numerical Propagator
293         final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
294         NumPropagator.setOrbitType(propagationType);
295 
296         // Atmosphere
297         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
298                                                             FramesFactory.getITRF(IERSConventions.IERS_2010, true));
299         MarshallSolarActivityFutureEstimation msafe =
300                         new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
301                                                                   MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
302         DataContext.getDefault().getDataProvidersManager().feed(msafe.getSupportedNames(), msafe);
303         DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
304 
305         // Force model
306         final ForceModel holmesFeatherstone =
307                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
308         final ForceModel drag =
309                         new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0));
310         NumPropagator.addForceModel(holmesFeatherstone);
311         NumPropagator.addForceModel(drag);
312 
313         // Set up initial state in the propagator
314         NumPropagator.setInitialState(initialState);
315 
316         // Extrapolate from the initial to the final date
317         final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
318         final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
319 
320         //_______________________________________________________________________________________________
321         // SET UP A BROUWER LYDDANE PROPAGATION WITHOUT DRAG
322         //_______________________________________________________________________________________________
323 
324         BrouwerLyddanePropagator BLextrapolator =
325                         new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
326 
327         SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
328         KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
329 
330         // Verify a and e differences without the drag effect on Brouwer-Lyddane
331         final double deltaSmaBefore = 20.44;
332         final double deltaEccBefore = 1.0125e-4;
333         Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore);
334         Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore);
335 
336         //_______________________________________________________________________________________________
337         // SET UP A BROUWER LYDDANE PROPAGATION WITH DRAG
338         //_______________________________________________________________________________________________
339 
340         double M2 = 1.0e-14;
341         BLextrapolator = new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2);
342         BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
343         BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
344 
345         // Verify a and e differences without the drag effect on Brouwer-Lyddane
346         final double deltaSmaAfter = 15.66;
347         final double deltaEccAfter = 1.0121e-4;
348         Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter);
349         Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter);
350         Assertions.assertTrue(deltaSmaAfter < deltaSmaBefore);
351         Assertions.assertTrue(deltaEccAfter < deltaEccBefore);
352 
353     }
354 
355 
356     @Test
357     public void compareToNumericalPropagationMeanInitialOrbit() {
358 
359         final Frame inertialFrame = FramesFactory.getEME2000();
360         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
361         double timeshift = 60000. ;
362 
363         // Initial orbit
364         final double a = 24396159; // semi major axis in meters
365         final double e = 0.01; // eccentricity
366         final double i = FastMath.toRadians(7); // inclination
367         final double omega = FastMath.toRadians(180); // perigee argument
368         final double raan = FastMath.toRadians(261); // right ascention of ascending node
369         final double lM = 0; // mean anomaly
370         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
371                                                       inertialFrame, initDate, provider.getMu());
372 
373 
374         BrouwerLyddanePropagator BLextrapolator =
375                 new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider),
376                                              PropagationType.MEAN, BrouwerLyddanePropagator.M2);
377         SpacecraftState initialOsculatingState = BLextrapolator.propagate(initDate);
378         final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit());
379 
380         SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
381         final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
382 
383         //_______________________________________________________________________________________________
384         // SET UP A REFERENCE NUMERICAL PROPAGATION
385         //_______________________________________________________________________________________________
386 
387 
388         // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
389         final double minStep = 0.001;
390         final double maxstep = 1000.0;
391         final double positionTolerance = 10.0;
392         final OrbitType propagationType = OrbitType.KEPLERIAN;
393         final double[][] tolerances =
394                 NumericalPropagator.tolerances(positionTolerance, InitOrbit, propagationType);
395         final AdaptiveStepsizeIntegrator integrator =
396                 new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
397 
398         // Numerical Propagator
399         final NumericalPropagator NumPropagator = new NumericalPropagator(integrator);
400         NumPropagator.setOrbitType(propagationType);
401 
402         final ForceModel holmesFeatherstone =
403                 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
404         NumPropagator.addForceModel(holmesFeatherstone);
405 
406         // Set up initial state in the propagator
407         NumPropagator.setInitialState(initialOsculatingState);
408 
409         // Extrapolate from the initial to the final date
410         final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift));
411         final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit());
412 
413         //_______________________________________________________________________________________________
414         // SET UP A BROUWER LYDDANE PROPAGATION
415         //_______________________________________________________________________________________________
416 
417         Assertions.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.17);
418         Assertions.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
419         Assertions.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
420         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
421                 MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.197);
422         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
423                 MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.00072);
424         Assertions.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
425                 MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12);
426 
427     }
428 
429 
430     @Test
431     public void compareToNumericalPropagationResetInitialIntermediate() {
432 
433         final Frame inertialFrame = FramesFactory.getEME2000();
434         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
435         double timeshift = 60000.;
436 
437         // Initial orbit
438         final double a = 24396159; // semi major axis in meters
439         final double e = 0.01; // eccentricity
440         final double i = FastMath.toRadians(7); // inclination
441         final double omega = FastMath.toRadians(180); // perigee argument
442         final double raan = FastMath.toRadians(261); // right ascention of ascending node
443         final double lM = 0; // mean anomaly
444         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
445                                                       inertialFrame, initDate, provider.getMu());
446         // Initial state definition
447         final SpacecraftState initialState = new SpacecraftState(initialOrbit);
448 
449         //_______________________________________________________________________________________________
450         // SET UP A BROUWER LYDDANE PROPAGATOR
451         //_______________________________________________________________________________________________
452 
453         BrouwerLyddanePropagator BLextrapolator1 =
454                 new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider),
455                                              PropagationType.OSCULATING, BrouwerLyddanePropagator.M2);
456         //_______________________________________________________________________________________________
457         // SET UP ANOTHER BROUWER LYDDANE PROPAGATOR
458         //_______________________________________________________________________________________________
459 
460         BrouwerLyddanePropagator BLextrapolator2 =
461                 new BrouwerLyddanePropagator( new KeplerianOrbit(a + 3000, e + 0.001, i - FastMath.toRadians(12.0), omega, raan, lM, PositionAngle.TRUE,
462                         inertialFrame, initDate, provider.getMu()),DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
463         // Reset BL2 with BL1 initial state
464         BLextrapolator2.resetInitialState(initialState);
465 
466         SpacecraftState BLFinalState1 = BLextrapolator1.propagate(initDate.shiftedBy(timeshift));
467         final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit());
468         SpacecraftState BLFinalState2 = BLextrapolator2.propagate(initDate.shiftedBy(timeshift));
469         BLextrapolator2.resetIntermediateState(BLFinalState1, true);
470         final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit());
471 
472         Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
473         Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
474         Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
475         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
476                 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
477         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
478                 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
479         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
480                 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
481 
482     }
483 
484     @Test
485     public void compareConstructors() {
486 
487         final Frame inertialFrame = FramesFactory.getEME2000();
488         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
489         double timeshift = 600. ;
490 
491         // Initial orbit
492         final double a = 24396159; // semi major axis in meters
493         final double e = 0.01; // eccentricity
494         final double i = FastMath.toRadians(7); // inclination
495         final double omega = FastMath.toRadians(180); // perigee argument
496         final double raan = FastMath.toRadians(261); // right ascention of ascending node
497         final double lM = 0; // mean anomaly
498         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
499                                                       inertialFrame, initDate, provider.getMu());
500 
501 
502         BrouwerLyddanePropagator BLPropagator1 = new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW,
503                 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
504         BrouwerLyddanePropagator BLPropagator2 = new BrouwerLyddanePropagator(initialOrbit,
505                 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
506         BrouwerLyddanePropagator BLPropagator3 = new BrouwerLyddanePropagator(initialOrbit, Propagator.DEFAULT_MASS,
507                 provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
508 
509         SpacecraftState BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift));
510         final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit());
511         SpacecraftState BLFinalState2 = BLPropagator2.propagate(initDate.shiftedBy(timeshift));
512         final KeplerianOrbit BLOrbit2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState2.getOrbit());
513         SpacecraftState BLFinalState3 = BLPropagator3.propagate(initDate.shiftedBy(timeshift));
514         final KeplerianOrbit BLOrbit3 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState3.getOrbit());
515 
516 
517         Assertions.assertEquals(BLOrbit1.getA(), BLOrbit2.getA(), 0.0);
518         Assertions.assertEquals(BLOrbit1.getE(), BLOrbit2.getE(), 0.0);
519         Assertions.assertEquals(BLOrbit1.getI(), BLOrbit2.getI(), 0.0);
520         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
521                 MathUtils.normalizeAngle(BLOrbit2.getPerigeeArgument(), FastMath.PI), 0.0);
522         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
523                 MathUtils.normalizeAngle(BLOrbit2.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
524         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
525                 MathUtils.normalizeAngle(BLOrbit2.getTrueAnomaly(), FastMath.PI), 0.0);
526         Assertions.assertEquals(BLOrbit1.getA(), BLOrbit3.getA(), 0.0);
527         Assertions.assertEquals(BLOrbit1.getE(), BLOrbit3.getE(), 0.0);
528         Assertions.assertEquals(BLOrbit1.getI(), BLOrbit3.getI(), 0.0);
529         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getPerigeeArgument(), FastMath.PI),
530                 MathUtils.normalizeAngle(BLOrbit3.getPerigeeArgument(), FastMath.PI), 0.0);
531         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getRightAscensionOfAscendingNode(), FastMath.PI),
532                 MathUtils.normalizeAngle(BLOrbit3.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0);
533         Assertions.assertEquals(MathUtils.normalizeAngle(BLOrbit1.getTrueAnomaly(), FastMath.PI),
534                 MathUtils.normalizeAngle(BLOrbit3.getTrueAnomaly(), FastMath.PI), 0.0);
535 
536     }
537 
538 
539     @Test
540     public void undergroundOrbit() {
541         Assertions.assertThrows(OrekitException.class, () -> {
542             // for a semi major axis < equatorial radius
543             Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
544             Vector3D velocity = new Vector3D(-500.0, 800.0, 100.0);
545             AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
546             Orbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
547                     FramesFactory.getEME2000(), initDate, provider.getMu());
548             // Extrapolator definition
549             // -----------------------
550             BrouwerLyddanePropagator extrapolator =
551                     new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, provider.getAe(), provider.getMu(),
552                             -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
553 
554             // Extrapolation at the initial date
555             // ---------------------------------
556             double delta_t = 0.0;
557             AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
558             extrapolator.propagate(extrapDate);
559         });
560     }
561 
562 
563     @Test
564     public void tooEllipticalOrbit() {
565         Assertions.assertThrows(OrekitException.class, () -> {
566             // for an eccentricity too big for the model
567             AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
568             Orbit initialOrbit = new KeplerianOrbit(67679244.0, 1.0,  1.85850, 2.1, 2.9,
569                     6.2, PositionAngle.TRUE, FramesFactory.getEME2000(),
570                     initDate, provider.getMu());
571             // Extrapolator definition
572             // -----------------------
573             BrouwerLyddanePropagator extrapolator =
574                     new BrouwerLyddanePropagator(initialOrbit, provider.getAe(), provider.getMu(),
575                             -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2);
576 
577             // Extrapolation at the initial date
578             // ---------------------------------
579             double delta_t = 0.0;
580             AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
581             extrapolator.propagate(extrapDate);
582         });
583     }
584 
585 
586     @Test
587     public void criticalInclination() {
588 
589         final Frame inertialFrame = FramesFactory.getEME2000();
590 
591         final double a = 24396159; // semi major axis in meters
592         final double e = 0.01; // eccentricity
593         final double i = FastMath.acos(1.0 / FastMath.sqrt(5.0)); // critical inclination
594         final double omega = FastMath.toRadians(180); // perigee argument
595         final double raan = FastMath.toRadians(261); // right ascention of ascending node
596         final double lM = 0; // mean anomaly
597 
598         AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
599         final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
600                                                       inertialFrame, initDate, provider.getMu());
601 
602         // Extrapolator definition
603         // -----------------------
604         BrouwerLyddanePropagator extrapolator =
605             new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
606 
607         // Extrapolation at the initial date
608         // ---------------------------------
609         SpacecraftState finalOrbit = extrapolator.propagate(initDate);
610 
611         // Verify
612         Assertions.assertEquals(0.0,
613                             Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
614                                               finalOrbit.getPVCoordinates().getPosition()),
615                             1.9e-8);
616 
617         Assertions.assertEquals(0.0,
618                             Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
619                                               finalOrbit.getPVCoordinates().getVelocity()),
620                             3.0e-12);
621         Assertions.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0);
622 
623     }
624 
625 
626     @Test
627     public void insideEarth() {
628         Assertions.assertThrows(OrekitException.class, () -> {
629             // for an eccentricity too big for the model
630             AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
631             Orbit initialOrbit = new KeplerianOrbit( provider.getAe()-1000.0, 0.01, FastMath.toRadians(10.0), 2.1, 2.9,
632                     6.2, PositionAngle.TRUE, FramesFactory.getEME2000(),
633                     initDate, provider.getMu());
634             // Extrapolator definition
635             // -----------------------
636             BrouwerLyddanePropagator extrapolator =
637                     new BrouwerLyddanePropagator(initialOrbit, Propagator.DEFAULT_MASS, provider.getAe(), provider.getMu(),
638                             -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7);
639 
640             // Extrapolation at the initial date
641             // ---------------------------------
642             double delta_t = 0.0;
643             AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
644             extrapolator.propagate(extrapDate);
645         });
646     }
647 
648     @Test
649     public void testUnableToComputeBLMeanParameters() {
650         Assertions.assertThrows(OrekitException.class, () -> {
651 
652             final Frame inertialFrame = FramesFactory.getEME2000();
653             AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
654 
655             // Initial orbit
656             final double a = 24396159; // semi major axis in meters
657             final double e = 0.9; // eccentricity
658             final double i = FastMath.toRadians(7); // inclination
659             final double omega = FastMath.toRadians(180); // perigee argument
660             final double raan = FastMath.toRadians(261); // right ascention of ascending node
661             final double lM = FastMath.toRadians(0); // mean anomaly
662             final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE,
663                     inertialFrame, initDate, provider.getMu());
664             // Extrapolator definition
665             // -----------------------
666             BrouwerLyddanePropagator extrapolator =
667                     new BrouwerLyddanePropagator(initialOrbit,  GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2);
668 
669             // Extrapolation at the initial date
670             // ---------------------------------
671             double delta_t = 0.0;
672             AbsoluteDate extrapDate = initDate.shiftedBy(delta_t);
673             extrapolator.propagate(extrapDate);
674         });
675     }
676 
677     @Test
678     public void testMeanOrbit() {
679         final KeplerianOrbit initialOsculating =
680                         new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
681                                            FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
682                                            provider.getMu());
683         final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
684         final UnnormalizedSphericalHarmonics ush = ushp.onDate(initialOsculating.getDate());
685 
686         // set up a reference numerical propagator starting for the specified start orbit
687         // using the same force models (i.e. the first few zonal terms)
688         double[][] tol = NumericalPropagator.tolerances(0.1, initialOsculating, OrbitType.KEPLERIAN);
689         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
690         integrator.setInitialStepSize(60);
691         NumericalPropagator num = new NumericalPropagator(integrator);
692         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
693         num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
694         num.setInitialState(new SpacecraftState(initialOsculating));
695         num.setOrbitType(OrbitType.KEPLERIAN);
696         final StorelessUnivariateStatistic oscMin  = new Min();
697         final StorelessUnivariateStatistic oscMax  = new Max();
698         final StorelessUnivariateStatistic meanMin = new Min();
699         final StorelessUnivariateStatistic meanMax = new Max();
700         num.getMultiplexer().add(60, state -> {
701             final Orbit osc = state.getOrbit();
702             oscMin.increment(osc.getA());
703             oscMax.increment(osc.getA());
704             // compute mean orbit at current date (this is what we test)
705             final Orbit mean = BrouwerLyddanePropagator.computeMeanOrbit(state.getOrbit(), ushp, ush, BrouwerLyddanePropagator.M2);
706             meanMin.increment(mean.getA());
707             meanMax.increment(mean.getA());
708         });
709         num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
710 
711         Assertions.assertEquals(3188.347, oscMax.getResult()  - oscMin.getResult(),  1.0e-3);
712         Assertions.assertEquals(  55.540, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
713 
714     }
715 
716     @BeforeEach
717     public void setUp() {
718         Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
719         provider = GravityFieldFactory.getNormalizedProvider(5, 0);
720     }
721 
722     @AfterEach
723     public void tearDown() {
724         provider = null;
725     }
726 
727     private NormalizedSphericalHarmonicsProvider provider;
728 
729 }
730 
731