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