1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.semianalytical.dsst;
18  
19  import java.io.IOException;
20  import java.text.ParseException;
21  import java.util.ArrayList;
22  import java.util.Collection;
23  import java.util.Collections;
24  import java.util.HashSet;
25  import java.util.List;
26  import java.util.Set;
27  import java.util.stream.Stream;
28  
29  import org.hamcrest.MatcherAssert;
30  import org.hipparchus.Field;
31  import org.hipparchus.CalculusFieldElement;
32  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
33  import org.hipparchus.geometry.euclidean.threed.RotationOrder;
34  import org.hipparchus.geometry.euclidean.threed.Vector3D;
35  import org.hipparchus.ode.ODEIntegrator;
36  import org.hipparchus.ode.events.Action;
37  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
38  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
39  import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
40  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
41  import org.hipparchus.util.FastMath;
42  import org.hipparchus.util.MathArrays;
43  import org.hipparchus.util.MathUtils;
44  import org.junit.After;
45  import org.junit.Assert;
46  import org.junit.Before;
47  import org.junit.Test;
48  import org.orekit.OrekitMatchers;
49  import org.orekit.Utils;
50  import org.orekit.attitudes.AttitudeProvider;
51  import org.orekit.attitudes.LofOffset;
52  import org.orekit.bodies.CelestialBody;
53  import org.orekit.bodies.CelestialBodyFactory;
54  import org.orekit.bodies.OneAxisEllipsoid;
55  import org.orekit.errors.OrekitException;
56  import org.orekit.forces.AbstractForceModel;
57  import org.orekit.forces.BoxAndSolarArraySpacecraft;
58  import org.orekit.forces.ForceModel;
59  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
60  import org.orekit.forces.gravity.potential.GRGSFormatReader;
61  import org.orekit.forces.gravity.potential.GravityFieldFactory;
62  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
63  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
64  import org.orekit.forces.gravity.potential.TideSystem;
65  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
66  import org.orekit.forces.maneuvers.ImpulseManeuver;
67  import org.orekit.frames.Frame;
68  import org.orekit.frames.FramesFactory;
69  import org.orekit.frames.LOFType;
70  import org.orekit.models.earth.atmosphere.Atmosphere;
71  import org.orekit.models.earth.atmosphere.HarrisPriester;
72  import org.orekit.orbits.CartesianOrbit;
73  import org.orekit.orbits.CircularOrbit;
74  import org.orekit.orbits.EquinoctialOrbit;
75  import org.orekit.orbits.KeplerianOrbit;
76  import org.orekit.orbits.Orbit;
77  import org.orekit.orbits.OrbitType;
78  import org.orekit.orbits.PositionAngle;
79  import org.orekit.propagation.BoundedPropagator;
80  import org.orekit.propagation.EphemerisGenerator;
81  import org.orekit.propagation.FieldSpacecraftState;
82  import org.orekit.propagation.PropagationType;
83  import org.orekit.propagation.Propagator;
84  import org.orekit.propagation.SpacecraftState;
85  import org.orekit.propagation.events.AltitudeDetector;
86  import org.orekit.propagation.events.DateDetector;
87  import org.orekit.propagation.events.EventDetector;
88  import org.orekit.propagation.events.FieldEventDetector;
89  import org.orekit.propagation.events.LatitudeCrossingDetector;
90  import org.orekit.propagation.events.NodeDetector;
91  import org.orekit.propagation.events.handlers.EventHandler;
92  import org.orekit.propagation.numerical.NumericalPropagator;
93  import org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution;
94  import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
95  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
96  import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
97  import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
98  import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
99  import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
100 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
101 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
102 import org.orekit.time.AbsoluteDate;
103 import org.orekit.time.DateComponents;
104 import org.orekit.time.TimeComponents;
105 import org.orekit.time.TimeScale;
106 import org.orekit.time.TimeScalesFactory;
107 import org.orekit.utils.Constants;
108 import org.orekit.utils.IERSConventions;
109 import org.orekit.utils.PVCoordinates;
110 import org.orekit.utils.ParameterDriver;
111 import org.orekit.utils.TimeStampedPVCoordinates;
112 
113 public class DSSTPropagatorTest {
114 
115     private DSSTPropagator dsstProp;
116 
117     @Test
118     public void testIssue363() {
119         Utils.setDataRoot("regular-data");
120         AbsoluteDate date = new AbsoluteDate("2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
121         CircularOrbit orbit = new CircularOrbit(7389068.5, 1.0e-15, 1.0e-15, 1.709573, 1.308398, 0, PositionAngle.MEAN,
122                                                 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
123                                                 date, Constants.WGS84_EARTH_MU);
124         SpacecraftState osculatingState = new SpacecraftState(orbit, 1116.2829);
125 
126         List<DSSTForceModel> dsstForceModels = new ArrayList<DSSTForceModel>();
127 
128         dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu()));
129         dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu()));
130 
131         SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
132         Assert.assertEquals( 0.421,   osculatingState.getA()             - meanState.getA(),             1.0e-3);
133         Assert.assertEquals(-5.23e-8, osculatingState.getEquinoctialEx() - meanState.getEquinoctialEx(), 1.0e-10);
134         Assert.assertEquals(15.22e-8, osculatingState.getEquinoctialEy() - meanState.getEquinoctialEy(), 1.0e-10);
135         Assert.assertEquals(-3.15e-8, osculatingState.getHx()            - meanState.getHx(),            1.0e-10);
136         Assert.assertEquals( 2.83e-8, osculatingState.getHy()            - meanState.getHy(),            1.0e-10);
137         Assert.assertEquals(15.96e-8, osculatingState.getLM()            - meanState.getLM(),            1.0e-10);
138 
139     }
140 
141     @Test
142     public void testIssue364() {
143         Utils.setDataRoot("regular-data");
144         AbsoluteDate date = new AbsoluteDate("2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
145         CircularOrbit orbit = new CircularOrbit(7389068.5, 0.0, 0.0, 1.709573, 1.308398, 0, PositionAngle.MEAN,
146                                                 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
147                                                 date, Constants.WGS84_EARTH_MU);
148         SpacecraftState osculatingState = new SpacecraftState(orbit, 1116.2829);
149 
150         List<DSSTForceModel> dsstForceModels = new ArrayList<DSSTForceModel>();
151 
152         dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu()));
153         dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu()));
154 
155         SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
156         Assert.assertEquals( 0.421,   osculatingState.getA()             - meanState.getA(),             1.0e-3);
157         Assert.assertEquals(-5.23e-8, osculatingState.getEquinoctialEx() - meanState.getEquinoctialEx(), 1.0e-10);
158         Assert.assertEquals(15.22e-8, osculatingState.getEquinoctialEy() - meanState.getEquinoctialEy(), 1.0e-10);
159         Assert.assertEquals(-3.15e-8, osculatingState.getHx()            - meanState.getHx(),            1.0e-10);
160         Assert.assertEquals( 2.83e-8, osculatingState.getHy()            - meanState.getHy(),            1.0e-10);
161         Assert.assertEquals(15.96e-8, osculatingState.getLM()            - meanState.getLM(),            1.0e-10);
162 
163     }
164 
165     @Test
166     public void testHighDegreesSetting() {
167 
168         Utils.setDataRoot("regular-data:potential/grgs-format");
169         GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
170         int earthDegree = 36;
171         int earthOrder  = 36;
172         int eccPower    = 4;
173         final UnnormalizedSphericalHarmonicsProvider provider =
174                 GravityFieldFactory.getUnnormalizedProvider(earthDegree, earthOrder);
175         final org.orekit.frames.Frame earthFrame =
176                 FramesFactory.getITRF(IERSConventions.IERS_2010, true); // terrestrial frame
177         final DSSTForceModel force =
178                 new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
179                                  earthDegree, earthOrder, eccPower, earthDegree + eccPower,
180                                  earthDegree, earthOrder, eccPower);
181         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
182         forces.add(force);
183         TimeScale tai = TimeScalesFactory.getTAI();
184         AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
185         Frame eci = FramesFactory.getGCRF();
186         KeplerianOrbit orbit = new KeplerianOrbit(
187                 7120000.0, 1.0e-3, FastMath.toRadians(60.0),
188                 FastMath.toRadians(120.0), FastMath.toRadians(47.0),
189                 FastMath.toRadians(12.0),
190                 PositionAngle.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
191         SpacecraftState oscuState = DSSTPropagator.computeOsculatingState(new SpacecraftState(orbit), null, forces);
192         Assert.assertEquals(7119927.097122, oscuState.getA(), 0.001);
193     }
194 
195     @Test
196     public void testEphemerisDates() {
197         //setup
198         TimeScale tai = TimeScalesFactory.getTAI();
199         AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
200         AbsoluteDate startDate = new AbsoluteDate("2015-07-03", tai).shiftedBy(-0.1);
201         AbsoluteDate endDate = new AbsoluteDate("2015-07-04", tai);
202         Frame eci = FramesFactory.getGCRF();
203         KeplerianOrbit orbit = new KeplerianOrbit(
204                 600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0, 0, 0, 0,
205                 PositionAngle.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
206         double[][] tol = DSSTPropagator
207                 .tolerances(1, orbit);
208         Propagator prop = new DSSTPropagator(
209                 new DormandPrince853Integrator(0.1, 500, tol[0], tol[1]));
210         prop.resetInitialState(new SpacecraftState(new CartesianOrbit(orbit)));
211 
212         //action
213         EphemerisGenerator generator = prop.getEphemerisGenerator();
214         prop.propagate(startDate, endDate);
215         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
216 
217         //verify
218         TimeStampedPVCoordinates actualPV = ephemeris.getPVCoordinates(startDate, eci);
219         TimeStampedPVCoordinates expectedPV = orbit.getPVCoordinates(startDate, eci);
220         MatcherAssert.assertThat(actualPV.getPosition(),
221                 OrekitMatchers.vectorCloseTo(expectedPV.getPosition(), 1.0));
222         MatcherAssert.assertThat(actualPV.getVelocity(),
223                 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity(), 1.0));
224         MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate),
225                 OrekitMatchers.closeTo(0, 0));
226         MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate),
227                 OrekitMatchers.closeTo(0, 0));
228         //test date
229         AbsoluteDate date = endDate.shiftedBy(-0.11);
230         Assert.assertEquals(
231                 ephemeris.propagate(date).getDate().durationFrom(date), 0, 0);
232     }
233 
234     @Test
235     public void testNoExtrapolation() {
236         SpacecraftState state = getLEOState();
237         setDSSTProp(state);
238 
239         // Propagation of the initial state at the initial date
240         final SpacecraftState finalState = dsstProp.propagate(state.getDate());
241 
242         // Initial orbit definition
243         final Vector3D initialPosition = state.getPVCoordinates().getPosition();
244         final Vector3D initialVelocity = state.getPVCoordinates().getVelocity();
245 
246         // Final orbit definition
247         final Vector3D finalPosition = finalState.getPVCoordinates().getPosition();
248         final Vector3D finalVelocity = finalState.getPVCoordinates().getVelocity();
249 
250         // Check results
251         Assert.assertEquals(initialPosition.getX(), finalPosition.getX(), 0.0);
252         Assert.assertEquals(initialPosition.getY(), finalPosition.getY(), 0.0);
253         Assert.assertEquals(initialPosition.getZ(), finalPosition.getZ(), 0.0);
254         Assert.assertEquals(initialVelocity.getX(), finalVelocity.getX(), 0.0);
255         Assert.assertEquals(initialVelocity.getY(), finalVelocity.getY(), 0.0);
256         Assert.assertEquals(initialVelocity.getZ(), finalVelocity.getZ(), 0.0);
257     }
258 
259     @Test
260     public void testKepler() {
261         SpacecraftState state = getLEOState();
262         setDSSTProp(state);
263 
264         Assert.assertEquals(2, dsstProp.getSatelliteRevolution());
265         dsstProp.setSatelliteRevolution(17);
266         Assert.assertEquals(17, dsstProp.getSatelliteRevolution());
267 
268         // Propagation of the initial state at t + dt
269         final double dt = 3200.;
270         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
271 
272         // Check results
273         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
274         Assert.assertEquals(state.getA(), finalState.getA(), 0.);
275         Assert.assertEquals(state.getEquinoctialEx(), finalState.getEquinoctialEx(), 0.);
276         Assert.assertEquals(state.getEquinoctialEy(), finalState.getEquinoctialEy(), 0.);
277         Assert.assertEquals(state.getHx(), finalState.getHx(), 0.);
278         Assert.assertEquals(state.getHy(), finalState.getHy(), 0.);
279         Assert.assertEquals(state.getLM() + n * dt, finalState.getLM(), 1.e-14);
280 
281     }
282 
283     @Test
284     public void testEphemeris() {
285         SpacecraftState state = getGEOState();
286         setDSSTProp(state);
287 
288         // Set ephemeris mode
289         EphemerisGenerator generator = dsstProp.getEphemerisGenerator();
290 
291         // Propagation of the initial state at t + 10 days
292         final double dt = 2. * Constants.JULIAN_DAY;
293         dsstProp.propagate(state.getDate().shiftedBy(5. * dt));
294 
295         // Get ephemeris
296         BoundedPropagator ephem = generator.getGeneratedEphemeris();
297 
298         // Propagation of the initial state with ephemeris at t + 2 days
299         final SpacecraftState s = ephem.propagate(state.getDate().shiftedBy(dt));
300 
301         // Check results
302         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
303         Assert.assertEquals(state.getA(), s.getA(), 0.);
304         Assert.assertEquals(state.getEquinoctialEx(), s.getEquinoctialEx(), 0.);
305         Assert.assertEquals(state.getEquinoctialEy(), s.getEquinoctialEy(), 0.);
306         Assert.assertEquals(state.getHx(), s.getHx(), 0.);
307         Assert.assertEquals(state.getHy(), s.getHy(), 0.);
308         Assert.assertEquals(state.getLM() + n * dt, s.getLM(), 1.5e-14);
309 
310     }
311 
312     @Test
313     public void testImpulseManeuver() {
314         final Orbit initialOrbit = new KeplerianOrbit(24532000.0, 0.72, 0.3, FastMath.PI, 0.4, 2.0, PositionAngle.MEAN, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2008, 06, 23), new TimeComponents(14, 18, 37), TimeScalesFactory.getUTC()), 3.986004415e14);
315         final double a = initialOrbit.getA();
316         final double e = initialOrbit.getE();
317         final double i = initialOrbit.getI();
318         final double mu = initialOrbit.getMu();
319         final double vApo = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
320         double dv = 0.99 * FastMath.tan(i) * vApo;
321 
322         // Set propagator with state
323         setDSSTProp(new SpacecraftState(initialOrbit));
324 
325         // Add impulse maneuver
326         dsstProp.setAttitudeProvider(new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS));
327         dsstProp.addEventDetector(new ImpulseManeuver<NodeDetector>(new NodeDetector(initialOrbit, FramesFactory.getEME2000()), new Vector3D(dv, Vector3D.PLUS_J), 400.0));
328         SpacecraftState propagated = dsstProp.propagate(initialOrbit.getDate().shiftedBy(8000));
329 
330         Assert.assertEquals(0.0028257, propagated.getI(), 1.0e-6);
331     }
332 
333     @Test
334     public void testPropagationWithCentralBody() throws Exception {
335 
336         // Central Body geopotential 4x4
337         final UnnormalizedSphericalHarmonicsProvider provider =
338                 GravityFieldFactory.getUnnormalizedProvider(4, 4);
339         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
340 
341         // GPS Orbit
342         final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400,
343                                                        TimeScalesFactory.getUTC());
344         final Orbit orbit = new KeplerianOrbit(26559890.,
345                                                0.0041632,
346                                                FastMath.toRadians(55.2),
347                                                FastMath.toRadians(315.4985),
348                                                FastMath.toRadians(130.7562),
349                                                FastMath.toRadians(44.2377),
350                                                PositionAngle.MEAN,
351                                                FramesFactory.getEME2000(),
352                                                initDate,
353                                                provider.getMu());
354 
355         // Set propagator with state and force model
356         setDSSTProp(new SpacecraftState(orbit));
357         dsstProp.addForceModel(new DSSTZonal(provider, 4, 3, 9));
358         dsstProp.addForceModel(new DSSTTesseral(earthFrame,
359                                                 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
360                                                 4, 4, 4, 8, 4, 4, 2));
361 
362         // 5 days propagation
363         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
364 
365         // Ref GTDS_DSST:
366         // a    = 26559.92081 km
367         // h/ey =   0.2731622444E-03
368         // k/ex =   0.4164167597E-02
369         // p/hy =  -0.3399607878
370         // q/hx =   0.3971568634
371         // lM   = 140.6375352°
372         Assert.assertEquals(26559920.81, state.getA(), 1.e-1);
373         Assert.assertEquals(0.2731622444E-03, state.getEquinoctialEx(), 2.e-8);
374         Assert.assertEquals(0.4164167597E-02, state.getEquinoctialEy(), 2.e-8);
375         Assert.assertEquals(-0.3399607878, state.getHx(), 5.e-8);
376         Assert.assertEquals(0.3971568634, state.getHy(), 2.e-6);
377         Assert.assertEquals(140.6375352,
378                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
379                             5.e-3);
380     }
381 
382     @Test
383     public void testPropagationWithThirdBody() throws IOException {
384 
385         // Central Body geopotential 2x0
386         final UnnormalizedSphericalHarmonicsProvider provider =
387                 GravityFieldFactory.getUnnormalizedProvider(2, 0);
388         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
389         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
390         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
391                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
392                                                    provider, 2, 0, 0, 2, 2, 0, 0);
393 
394         // Third Bodies Force Model (Moon + Sun)
395         DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
396         DSSTForceModel sun  = new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu());
397 
398         // SIRIUS Orbit
399         final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
400                                                        TimeScalesFactory.getUTC());
401         final Orbit orbit = new KeplerianOrbit(42163393.,
402                                                0.2684,
403                                                FastMath.toRadians(63.435),
404                                                FastMath.toRadians(270.0),
405                                                FastMath.toRadians(285.0),
406                                                FastMath.toRadians(344.0),
407                                                PositionAngle.MEAN,
408                                                FramesFactory.getEME2000(),
409                                                initDate,
410                                                provider.getMu());
411 
412         // Set propagator with state and force model
413         setDSSTProp(new SpacecraftState(orbit));
414         dsstProp.addForceModel(zonal);
415         dsstProp.addForceModel(tesseral);
416         dsstProp.addForceModel(moon);
417         dsstProp.addForceModel(sun);
418 
419         // 5 days propagation
420         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
421 
422         // Ref Standalone_DSST:
423         // a    = 42163393.0 m
424         // h/ey =  -0.06893353670734315
425         // k/ex =  -0.2592789733084587
426         // p/hy =  -0.5968524904937771
427         // q/hx =   0.1595005111738418
428         // lM   = 183°9386620425922
429         Assert.assertEquals(42163393.0, state.getA(), 1.e-1);
430         Assert.assertEquals(-0.2592789733084587, state.getEquinoctialEx(), 5.e-7);
431         Assert.assertEquals(-0.06893353670734315, state.getEquinoctialEy(), 2.e-7);
432         Assert.assertEquals( 0.1595005111738418, state.getHx(), 2.e-7);
433         Assert.assertEquals(-0.5968524904937771, state.getHy(), 5.e-8);
434         Assert.assertEquals(183.9386620425922,
435                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
436                             3.e-2);
437     }
438 
439     @Test(expected=OrekitException.class)
440     public void testTooSmallMaxDegree() {
441         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 1, 0, 3);
442     }
443 
444     @Test(expected=OrekitException.class)
445     public void testTooLargeMaxDegree() {
446         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 8, 0, 8);
447     }
448 
449     @Test(expected=OrekitException.class)
450     public void testWrongMaxPower() {
451         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(8, 8), 4, 4, 4);
452     }
453 
454     @Test
455     public void testPropagationWithDrag() {
456 
457         // Central Body geopotential 2x0
458         final UnnormalizedSphericalHarmonicsProvider provider =
459                 GravityFieldFactory.getUnnormalizedProvider(2, 0);
460         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
461         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 0, 5);
462         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
463                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
464                                                    provider, 2, 0, 0, 2, 2, 0, 0);
465 
466         // Drag Force Model
467         final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
468                                                             Constants.WGS84_EARTH_FLATTENING,
469                                                             earthFrame);
470         final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
471 
472         final double cd = 2.0;
473         final double area = 25.0;
474         DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area, provider.getMu());
475 
476         // LEO Orbit
477         final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
478                                                        TimeScalesFactory.getUTC());
479         final Orbit orbit = new KeplerianOrbit(7204535.848109440,
480                                                0.0012402238462686,
481                                                FastMath.toRadians(98.74341600466740),
482                                                FastMath.toRadians(111.1990175076630),
483                                                FastMath.toRadians(43.32990110790340),
484                                                FastMath.toRadians(68.66852509725620),
485                                                PositionAngle.MEAN,
486                                                FramesFactory.getEME2000(),
487                                                initDate,
488                                                provider.getMu());
489 
490         // Set propagator with state and force model
491         setDSSTProp(new SpacecraftState(orbit));
492         dsstProp.addForceModel(zonal);
493         dsstProp.addForceModel(tesseral);
494         dsstProp.addForceModel(drag);
495 
496         // 5 days propagation
497         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
498 
499         // Ref Standalone_DSST:
500         // a    = 7204521.657141485 m
501         // h/ey =  0.0007093755541595772
502         // k/ex = -0.001016800430994036
503         // p/hy =  0.8698955648709271
504         // q/hx =  0.7757573478894775
505         // lM   = 193°0939742953394
506         Assert.assertEquals(7204521.657141485, state.getA(), 6.e-1);
507         Assert.assertEquals(-0.001016800430994036, state.getEquinoctialEx(), 5.e-8);
508         Assert.assertEquals(0.0007093755541595772, state.getEquinoctialEy(), 2.e-8);
509         Assert.assertEquals(0.7757573478894775, state.getHx(), 5.e-8);
510         Assert.assertEquals(0.8698955648709271, state.getHy(), 5.e-8);
511         Assert.assertEquals(193.0939742953394,
512                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
513                             2.e-3);
514         //Assert.assertEquals(((DSSTAtmosphericDrag)drag).getCd(), cd, 1e-9);
515         //Assert.assertEquals(((DSSTAtmosphericDrag)drag).getArea(), area, 1e-9);
516         Assert.assertEquals(((DSSTAtmosphericDrag)drag).getAtmosphere(), atm);
517 
518         final double atmosphericMaxConstant = 1000000.0; //DSSTAtmosphericDrag.ATMOSPHERE_ALTITUDE_MAX
519         Assert.assertEquals(((DSSTAtmosphericDrag)drag).getRbar(), atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 1e-9);
520     }
521 
522     @Test
523     public void testPropagationWithSolarRadiationPressure() {
524 
525         // Central Body geopotential 2x0
526         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
527         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
528         DSSTForceModel tesseral = new DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(),
529                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
530                                                    provider, 2, 0, 0, 2, 2, 0, 0);
531 
532         // SRP Force Model
533         DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
534                                                             Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
535                                                             provider.getMu());
536 
537         // GEO Orbit
538         final AbsoluteDate initDate = new AbsoluteDate(2003, 9, 16, 0, 0, 00.000,
539                                                        TimeScalesFactory.getUTC());
540         final Orbit orbit = new KeplerianOrbit(42166258.,
541                                                0.0001,
542                                                FastMath.toRadians(0.001),
543                                                FastMath.toRadians(315.4985),
544                                                FastMath.toRadians(130.7562),
545                                                FastMath.toRadians(44.2377),
546                                                PositionAngle.MEAN,
547                                                FramesFactory.getGCRF(),
548                                                initDate,
549                                                provider.getMu());
550 
551         // Set propagator with state and force model
552         dsstProp = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(86400.));
553         dsstProp.setInitialState(new SpacecraftState(orbit), PropagationType.MEAN);
554         dsstProp.addForceModel(zonal);
555         dsstProp.addForceModel(tesseral);
556         dsstProp.addForceModel(srp);
557 
558         // 10 days propagation
559         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(10. * 86400.));
560 
561         // Ref Standalone_DSST:
562         // a    = 42166257.99807995 m
563         // h/ey = -0.1191876027555493D-03
564         // k/ex = -0.1781865038201885D-05
565         // p/hy =  0.6618387121369373D-05
566         // q/hx = -0.5624363171289686D-05
567         // lM   = 140°3496229467104
568         Assert.assertEquals(42166257.99807995, state.getA(), 0.8);
569         Assert.assertEquals(-0.1781865038201885e-05, state.getEquinoctialEx(), 3.e-7);
570         Assert.assertEquals(-0.1191876027555493e-03, state.getEquinoctialEy(), 4.e-6);
571         Assert.assertEquals(-0.5624363171289686e-05, state.getHx(), 4.e-9);
572         Assert.assertEquals( 0.6618387121369373e-05, state.getHy(), 3.e-10);
573         Assert.assertEquals(140.3496229467104,
574                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
575                             2.e-4);
576     }
577 
578     @Test
579     public void testStopEvent() {
580         SpacecraftState state = getLEOState();
581         setDSSTProp(state);
582 
583         final AbsoluteDate stopDate = state.getDate().shiftedBy(1000);
584         CheckingHandler<DateDetector> checking = new CheckingHandler<DateDetector>(Action.STOP);
585         dsstProp.addEventDetector(new DateDetector(stopDate).withHandler(checking));
586         checking.assertEvent(false);
587         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(3200));
588         checking.assertEvent(true);
589         Assert.assertEquals(0, finalState.getDate().durationFrom(stopDate), 1.0e-10);
590     }
591 
592     @Test
593     public void testContinueEvent() {
594         SpacecraftState state = getLEOState();
595         setDSSTProp(state);
596 
597         final AbsoluteDate resetDate = state.getDate().shiftedBy(1000);
598         CheckingHandler<DateDetector> checking = new CheckingHandler<DateDetector>(Action.CONTINUE);
599         dsstProp.addEventDetector(new DateDetector(resetDate).withHandler(checking));
600         final double dt = 3200;
601         checking.assertEvent(false);
602         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
603         checking.assertEvent(true);
604         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
605         Assert.assertEquals(state.getA(), finalState.getA(), 1.0e-10);
606         Assert.assertEquals(state.getEquinoctialEx(), finalState.getEquinoctialEx(), 1.0e-10);
607         Assert.assertEquals(state.getEquinoctialEy(), finalState.getEquinoctialEy(), 1.0e-10);
608         Assert.assertEquals(state.getHx(), finalState.getHx(), 1.0e-10);
609         Assert.assertEquals(state.getHy(), finalState.getHy(), 1.0e-10);
610         Assert.assertEquals(state.getLM() + n * dt, finalState.getLM(), 6.0e-10);
611     }
612 
613     @Test
614     public void testIssue157() {
615         Utils.setDataRoot("regular-data:potential/icgem-format");
616         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
617         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
618         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
619                                          FramesFactory.getTOD(false),
620                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
621                                          nshp.getMu());
622         double period = orbit.getKeplerianPeriod();
623         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
624         AdaptiveStepsizeIntegrator integrator =
625                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
626         integrator.setInitialStepSize(10 * period);
627         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.MEAN);
628         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
629                                                       Constants.WGS84_EARTH_FLATTENING,
630                                                       FramesFactory.getGTOD(false));
631         CelestialBody sun = CelestialBodyFactory.getSun();
632         CelestialBody moon = CelestialBodyFactory.getMoon();
633         propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
634         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
635                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
636                                                   nshp, 8, 8, 4, 12, 8, 8, 4));
637         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
638         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
639         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
640         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
641 
642 
643         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.OSCULATING);
644         SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
645         // the following comparison is in fact meaningless
646         // the initial orbit is osculating the final orbit is a mean orbit
647         // and they are not considered at the same epoch
648         // we keep it only as is was an historical test
649         Assert.assertEquals(2189.4, orbit.getA() - finalState.getA(), 1.0);
650 
651         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
652         finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
653         // the following comparison is realistic
654         // both the initial orbit and final orbit are mean orbits
655         Assert.assertEquals(1478.05, orbit.getA() - finalState.getA(), 1.0);
656 
657     }
658 
659     @Test
660     public void testEphemerisGeneration() {
661         Utils.setDataRoot("regular-data:potential/icgem-format");
662         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
663         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
664         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
665                                          FramesFactory.getTOD(false),
666                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
667                                          nshp.getMu());
668         double period = orbit.getKeplerianPeriod();
669         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
670         AdaptiveStepsizeIntegrator integrator =
671                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
672         integrator.setInitialStepSize(10 * period);
673         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
674         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
675                                                       Constants.WGS84_EARTH_FLATTENING,
676                                                       FramesFactory.getGTOD(false));
677         CelestialBody sun = CelestialBodyFactory.getSun();
678         CelestialBody moon = CelestialBodyFactory.getMoon();
679         propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
680         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
681                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
682                                                   nshp, 8, 8, 4, 12, 8, 8, 4));
683         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
684         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
685         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
686         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
687         propagator.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
688 
689         // direct generation of states
690         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
691         final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
692         propagator.setStepHandler(600, currentState -> states.add(currentState));
693         propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
694 
695         // ephemeris generation
696         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
697         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
698         propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
699         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
700 
701         double maxError = 0;
702         for (final SpacecraftState state : states) {
703             final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
704             final double error = Vector3D.distance(state.getPVCoordinates().getPosition(),
705                                                    fromEphemeris.getPVCoordinates().getPosition());
706             maxError = FastMath.max(maxError, error);
707         }
708         Assert.assertEquals(0.0, maxError, 1.0e-10);
709     }
710 
711     @Test
712     public void testGetInitialOsculatingState() throws IllegalArgumentException, OrekitException {
713         final SpacecraftState initialState = getGEOState();
714 
715         // build integrator
716         final double minStep = initialState.getKeplerianPeriod() * 0.1;
717         final double maxStep = initialState.getKeplerianPeriod() * 10.0;
718         final double[][] tol = DSSTPropagator.tolerances(0.1, initialState.getOrbit());
719         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
720 
721         // build the propagator for the propagation of the mean elements
722         DSSTPropagator prop = new DSSTPropagator(integrator, PropagationType.MEAN);
723 
724         final UnnormalizedSphericalHarmonicsProvider provider =
725                 GravityFieldFactory.getUnnormalizedProvider(4, 0);
726         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
727         DSSTForceModel zonal    = new DSSTZonal(provider, 4, 3, 9);
728         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
729                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
730                                                    provider, 4, 0, 4, 8, 4, 0, 2);
731         prop.addForceModel(zonal);
732         prop.addForceModel(tesseral);
733 
734         // Set the initial state
735         prop.setInitialState(initialState, PropagationType.MEAN);
736         // Check the stored initial state is the osculating one
737         Assert.assertEquals(initialState, prop.getInitialState());
738         // Check that no propagation, i.e. propagation to the initial date, provides the initial
739         // osculating state although the propagator is configured to propagate mean elements !!!
740         Assert.assertEquals(initialState, prop.propagate(initialState.getDate()));
741     }
742 
743     @Test
744     public void testMeanToOsculatingState() throws IllegalArgumentException, OrekitException {
745         final SpacecraftState meanState = getGEOState();
746 
747         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
748         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
749         final DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
750         final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
751                                                          Constants.WGS84_EARTH_ANGULAR_VELOCITY,
752                                                          provider, 2, 0, 0, 2, 2, 0, 0);
753 
754         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
755         forces.add(zonal);
756         forces.add(tesseral);
757 
758         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
759         Assert.assertEquals(1559.1,
760                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
761                                               osculatingState.getPVCoordinates().getPosition()),
762                             1.0);
763     }
764 
765     @Test
766     public void testOsculatingToMeanState() throws IllegalArgumentException, OrekitException {
767         final SpacecraftState meanState = getGEOState();
768 
769         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
770         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
771 
772         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
773         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
774                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
775                                                    provider, 2, 0, 0, 2, 2, 0, 0);
776 
777         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
778         forces.add(zonal);
779         forces.add(tesseral);
780 
781         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
782 
783         // there are no Gaussian force models, we don't need an attitude provider
784         final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
785 
786         Assert.assertEquals(meanState.getA(), computedMeanState.getA(), 2.0e-8);
787         Assert.assertEquals(0.0,
788                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
789                                              computedMeanState.getPVCoordinates().getPosition()),
790                             2.0e-8);
791     }
792 
793     @Test
794     public void testShortPeriodCoefficients() {
795         Utils.setDataRoot("regular-data:potential/icgem-format");
796         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
797         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(4, 4);
798         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
799                                          FramesFactory.getTOD(false),
800                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
801                                          nshp.getMu());
802         double period = orbit.getKeplerianPeriod();
803         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
804         AdaptiveStepsizeIntegrator integrator =
805                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
806         integrator.setInitialStepSize(10 * period);
807         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
808         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
809                                                       Constants.WGS84_EARTH_FLATTENING,
810                                                       FramesFactory.getGTOD(false));
811         CelestialBody sun = CelestialBodyFactory.getSun();
812         CelestialBody moon = CelestialBodyFactory.getMoon();
813         propagator.addForceModel(new DSSTZonal(nshp, 4, 3, 9));
814         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
815                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
816                                                   nshp, 4, 4, 4, 8, 4, 4, 2));
817         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
818         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
819         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
820         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
821 
822         final AbsoluteDate finalDate = orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY);
823         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
824         final SpacecraftState stateNoConfig = propagator.propagate(finalDate);
825         Assert.assertEquals(0, stateNoConfig.getAdditionalStatesValues().size());
826 
827         Assert.assertNull(propagator.getSelectedCoefficients());
828         propagator.setSelectedCoefficients(new HashSet<String>());
829         Assert.assertNotNull(propagator.getSelectedCoefficients());
830         Assert.assertTrue(propagator.getSelectedCoefficients().isEmpty());
831         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
832         final SpacecraftState stateConfigEmpty = propagator.propagate(finalDate);
833         Assert.assertEquals(234, stateConfigEmpty.getAdditionalStatesValues().size());
834 
835         final Set<String> selected = new HashSet<String>();
836         selected.add("DSST-3rd-body-Moon-s[7]");
837         selected.add("DSST-central-body-tesseral-c[-2][3]");
838         propagator.setSelectedCoefficients(selected);
839         Assert.assertEquals(2, propagator.getSelectedCoefficients().size());
840         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
841         final SpacecraftState stateConfigeSelected = propagator.propagate(finalDate);
842         Assert.assertEquals(selected.size(), stateConfigeSelected.getAdditionalStatesValues().size());
843 
844         propagator.setSelectedCoefficients(null);
845         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
846         final SpacecraftState stateConfigNull = propagator.propagate(finalDate);
847         Assert.assertEquals(0, stateConfigNull.getAdditionalStatesValues().size());
848 
849     }
850 
851     @Test
852     public void testIssueMeanInclination() {
853 
854         final double earthAe = 6378137.0;
855         final double earthMu = 3.9860044E14;
856         final double earthJ2 = 0.0010826;
857 
858         // Initialize the DSST propagator with only J2 perturbation
859         Orbit orb = new KeplerianOrbit(new TimeStampedPVCoordinates(new AbsoluteDate("1992-10-08T15:20:38.821",
860                                                                                      TimeScalesFactory.getUTC()),
861                                                                     new Vector3D(5392808.809823, -4187618.3357927715, -44206.638015847195),
862                                                                     new Vector3D(2337.4472786270794, 2474.0146611860464, 6778.507766114648)),
863                                        FramesFactory.getTOD(false), earthMu);
864         final SpacecraftState ss = new SpacecraftState(orb);
865         final UnnormalizedSphericalHarmonicsProvider provider =
866               GravityFieldFactory.getUnnormalizedProvider(earthAe, earthMu, TideSystem.UNKNOWN,
867                                                           new double[][] { { 0.0 }, { 0.0 }, { -earthJ2 } },
868                                                           new double[][] { { 0.0 }, { 0.0 }, { 0.0 } });
869         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
870         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
871         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
872                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
873                                                    provider, 2, 0, 0, 2, 2, 0, 0);
874         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
875         forces.add(zonal);
876         forces.add(tesseral);
877         // Computes J2 mean elements using the DSST osculating to mean converter
878         final Orbit meanOrb = DSSTPropagator.computeMeanState(ss, null, forces).getOrbit();
879         Assert.assertEquals(0.0164196, FastMath.toDegrees(orb.getI() - meanOrb.getI()), 1.0e-7);
880     }
881 
882     @Test
883     public void testIssue257() {
884         final SpacecraftState meanState = getGEOState();
885 
886         // Third Bodies Force Model (Moon + Sun)
887         final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), meanState.getMu());
888         final DSSTForceModel sun  = new DSSTThirdBody(CelestialBodyFactory.getSun(), meanState.getMu());
889 
890         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
891         forces.add(moon);
892         forces.add(sun);
893 
894         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
895         Assert.assertEquals(734.3,
896                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
897                                               osculatingState.getPVCoordinates().getPosition()),
898                             1.0);
899 
900         final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
901         Assert.assertEquals(734.3,
902                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
903                                               computedMeanState.getPVCoordinates().getPosition()),
904                             1.0);
905 
906         Assert.assertEquals(0.0,
907                             Vector3D.distance(computedMeanState.getPVCoordinates().getPosition(),
908                                               meanState.getPVCoordinates().getPosition()),
909                             5.0e-6);
910 
911     }
912 
913     @Test
914     public void testIssue339() {
915 
916         final SpacecraftState osculatingState = getLEOState();
917 
918         final CelestialBody    sun   = CelestialBodyFactory.getSun();
919         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
920                                                             Constants.WGS84_EARTH_FLATTENING,
921                                                             FramesFactory.getITRF(IERSConventions.IERS_2010,
922                                                                                   true));
923         final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
924                                                                                      sun,
925                                                                                      50.0, Vector3D.PLUS_J,
926                                                                                      2.0, 0.1,
927                                                                                      0.2, 0.6);
928         final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
929         final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(),
930                                                                 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
931                                                                 0.0, 0.0, 0.0);
932 
933         // Surface force models that require an attitude provider
934         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
935         forces.add(new DSSTSolarRadiationPressure(sun,
936                                                   Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
937                                                   boxAndWing,
938                                                   osculatingState.getMu()));
939         forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getMu()));
940 
941         final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
942         Assert.assertEquals(0.522,
943                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
944                                               meanState.getPVCoordinates().getPosition()),
945                             0.001);
946 
947         final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
948         Assert.assertEquals(0.0,
949                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
950                                               computedOsculatingState.getPVCoordinates().getPosition()),
951                             5.0e-6);
952 
953     }
954 
955     @Test
956     public void testIssue613() {
957         // Spacecraft state
958         final SpacecraftState state = getLEOState();
959 
960         // Body frame
961         final Frame itrf = FramesFactory .getITRF(IERSConventions.IERS_2010, true);
962         
963         // Earth
964         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 4);
965         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
966                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
967         // Detectors
968         final List<EventDetector> events = new ArrayList<>();
969         events.add(new AltitudeDetector(85.5, earth));
970         events.add(new LatitudeCrossingDetector(earth, 0.0));
971 
972         // Force models
973         final List<DSSTForceModel> forceModels = new ArrayList<>();
974         forceModels.add(new DSSTZonal(provider));
975         forceModels.add(new DSSTTesseral(itrf, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
976         forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu()));
977         forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu()));
978 
979         // Set up DSST propagator
980         final double[][] tol = DSSTPropagator.tolerances(10.0, state.getOrbit());
981         final ODEIntegrator integrator = new DormandPrince54Integrator(60.0, 3600.0, tol[0], tol[1]);
982         final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
983         for (DSSTForceModel force : forceModels) {
984             propagator.addForceModel(force);
985         }
986         for (EventDetector event : events) {
987             propagator.addEventDetector(event);
988         }
989         propagator.setInitialState(state);
990 
991         // Propagation
992         final SpacecraftState finalState = propagator.propagate(state.getDate().shiftedBy(86400.0));
993 
994         // Verify is the propagation is correctly performed
995         Assert.assertEquals(finalState.getMu(), 3.986004415E14, Double.MIN_VALUE);
996     }
997 
998     @Test
999     public void testIssue339WithAccelerations() {
1000         final SpacecraftState osculatingState = getLEOStatePropagatedBy30Minutes();
1001         final CelestialBody sun = CelestialBodyFactory.getSun();
1002         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
1003                                                             FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1004         final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0, sun, 50.0, Vector3D.PLUS_J, 2.0, 0.1, 0.2, 0.6);
1005         final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
1006         final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(), LOFType.LVLH_CCSDS, RotationOrder.XYZ, 0.0, 0.0, 0.0);
1007         // Surface force models that require an attitude provider
1008         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
1009         forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getMu()));
1010         final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1011         final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1012         Assert.assertEquals(0.0, Vector3D.distance(osculatingState.getPVCoordinates().getPosition(), computedOsculatingState.getPVCoordinates().getPosition()),
1013                             5.0e-6);
1014     }
1015 
1016     @Test
1017     public void testIssue704() {
1018 
1019         // Coordinates
1020         final Orbit         orbit = getLEOState().getOrbit();
1021         final PVCoordinates pv    = orbit.getPVCoordinates();
1022 
1023         // dP
1024         final double dP = 10.0;
1025 
1026         // Computes dV
1027         final double r2 = pv.getPosition().getNormSq();
1028         final double v  = pv.getVelocity().getNorm();
1029         final double dV = orbit.getMu() * dP / (v * r2);
1030 
1031         // Verify
1032         final double[][] tol1 = DSSTPropagator.tolerances(dP, orbit);
1033         final double[][] tol2 = DSSTPropagator.tolerances(dP, dV, orbit);
1034         for (int i = 0; i < tol1.length; i++) {
1035             Assert.assertArrayEquals(tol1[i], tol2[i], Double.MIN_VALUE);
1036         }
1037 
1038     }
1039 
1040     /** This test is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1041     @Test
1042     public void testIssue670() {
1043 
1044         final NumericalForce force     = new NumericalForce();
1045         final DSSTForce      dsstForce = new DSSTForce(force, Constants.WGS84_EARTH_MU);
1046 
1047         SpacecraftState state = getLEOState();
1048         setDSSTProp(state);
1049         dsstProp.addForceModel(dsstForce);
1050 
1051         // Verify flag are false
1052         Assert.assertFalse(force.initialized);
1053         Assert.assertFalse(force.accComputed);
1054 
1055         // Propagation of the initial state at t + dt
1056         final double dt = 3200.;
1057         final AbsoluteDate target = state.getDate().shiftedBy(dt);
1058 
1059         dsstProp.propagate(target);
1060 
1061         // Flag must be true
1062         Assert.assertTrue(force.initialized);
1063         Assert.assertTrue(force.accComputed);
1064 
1065     }
1066 
1067 
1068     private SpacecraftState getGEOState() throws IllegalArgumentException, OrekitException {
1069         // No shadow at this date
1070         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
1071                                                        TimeScalesFactory.getUTC());
1072         final Orbit orbit = new EquinoctialOrbit(42164000,
1073                                                  10e-3,
1074                                                  10e-3,
1075                                                  FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
1076                                                  FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
1077                                                  PositionAngle.TRUE,
1078                                                  FramesFactory.getEME2000(),
1079                                                  initDate,
1080                                                  3.986004415E14);
1081         return new SpacecraftState(orbit);
1082     }
1083 
1084     private SpacecraftState getLEOState() throws IllegalArgumentException, OrekitException {
1085         final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1086         final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1087         // Spring equinoxe 21st mars 2003 1h00m
1088         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1089         return new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
1090                                                         FramesFactory.getEME2000(),
1091                                                         initDate,
1092                                                         3.986004415E14));
1093     }
1094 
1095     private void setDSSTProp(SpacecraftState initialState) {
1096         initialState.getDate();
1097         final double minStep = initialState.getKeplerianPeriod();
1098         final double maxStep = 100. * minStep;
1099         final double[][] tol = DSSTPropagator.tolerances(1.0, initialState.getOrbit());
1100         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1101         dsstProp = new DSSTPropagator(integrator);
1102         dsstProp.setInitialState(initialState, PropagationType.MEAN);
1103     }
1104 
1105     private static class CheckingHandler<T extends EventDetector> implements EventHandler<T> {
1106 
1107         private final Action actionOnEvent;
1108         private boolean gotHere;
1109 
1110         public CheckingHandler(final Action actionOnEvent) {
1111             this.actionOnEvent = actionOnEvent;
1112             this.gotHere       = false;
1113         }
1114 
1115         public void assertEvent(boolean expected) {
1116             Assert.assertEquals(expected, gotHere);
1117         }
1118 
1119         @Override
1120         public Action eventOccurred(SpacecraftState s, T detector, boolean increasing) {
1121             gotHere = true;
1122             return actionOnEvent;
1123         }
1124 
1125     }
1126 
1127     /** This class is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1128     private class DSSTForce extends AbstractGaussianContribution {
1129 
1130         DSSTForce(ForceModel contribution, double mu) {
1131             super("DSST mock -", 6.0e-10, contribution, mu);
1132         }
1133 
1134         /** {@inheritDoc} */
1135         @Override
1136         public EventDetector[] getEventsDetectors() {
1137             return null;
1138         }
1139 
1140         /** {@inheritDoc} */
1141         @Override
1142         public <T extends CalculusFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
1143             return null;
1144         }
1145 
1146         /** {@inheritDoc} */
1147         @Override
1148         protected List<ParameterDriver> getParametersDriversWithoutMu() {
1149             return Collections.emptyList();
1150         }
1151 
1152         /** {@inheritDoc} */
1153         @Override
1154         protected double[] getLLimits(SpacecraftState state,
1155                                       AuxiliaryElements auxiliaryElements) {
1156             return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
1157                                    FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0) };
1158         }
1159 
1160         /** {@inheritDoc} */
1161         @Override
1162         protected <T extends CalculusFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
1163                                                                  FieldAuxiliaryElements<T> auxiliaryElements) {
1164             final Field<T> field = state.getDate().getField();
1165             final T zero = field.getZero();
1166             final T[] tab = MathArrays.buildArray(field, 2);
1167             tab[0] = MathUtils.normalizeAngle(state.getLv(), zero).subtract(FastMath.PI);
1168             tab[1] = MathUtils.normalizeAngle(state.getLv(), zero).add(FastMath.PI);
1169             return tab;
1170         }
1171         
1172     }
1173 
1174     /** This class is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1175     private class NumericalForce extends AbstractForceModel {
1176 
1177         private boolean initialized;
1178         private boolean accComputed;
1179 
1180         NumericalForce() {
1181             this.initialized = false;
1182         }
1183 
1184         /** {@inheritDoc} */
1185         @Override
1186         public void init(final SpacecraftState s0, final AbsoluteDate t) {
1187             this.initialized = true;
1188             this.accComputed = false;
1189         }
1190 
1191         /** {@inheritDoc} */
1192         @Override
1193         public boolean dependsOnPositionOnly() {
1194             return false;
1195         }
1196 
1197         /** {@inheritDoc} */
1198         @Override
1199         public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1200             this.accComputed = true;
1201             return Vector3D.ZERO;
1202         }
1203 
1204         /** {@inheritDoc} */
1205         @Override
1206         public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1207             return FieldVector3D.getZero(s.getDate().getField());
1208         }
1209 
1210         /** {@inheritDoc} */
1211         @Override
1212         public Stream<EventDetector> getEventsDetectors() {
1213             return Stream.empty();
1214         }
1215 
1216         /** {@inheritDoc} */
1217         @Override
1218         public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
1219             return Stream.empty();
1220         }
1221 
1222 
1223         @Override
1224         public List<ParameterDriver> getParametersDrivers() {
1225             return Collections.emptyList();
1226         }
1227         
1228     }
1229 
1230     @Before
1231     public void setUp() throws IOException, ParseException {
1232         Utils.setDataRoot("regular-data:potential/shm-format");
1233     }
1234 
1235     @After
1236     public void tearDown() {
1237         dsstProp = null;
1238     }
1239 
1240     private SpacecraftState getLEOStatePropagatedBy30Minutes() throws IllegalArgumentException, OrekitException {
1241 
1242         final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1243         final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1244         // Spring equinoxe 21st mars 2003 1h00m
1245         final AbsoluteDate initialDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1246         final CartesianOrbit osculatingOrbit = new CartesianOrbit(new PVCoordinates(position, velocity), FramesFactory.getTOD(IERSConventions.IERS_1996, false),
1247                                                                   initialDate, Constants.WGS84_EARTH_MU);
1248         // Adaptive step integrator
1249         // with a minimum step of 0.001 and a maximum step of 1000
1250         double minStep = 0.001;
1251         double maxstep = 1000.0;
1252         double positionTolerance = 10.0;
1253         OrbitType propagationType = OrbitType.EQUINOCTIAL;
1254         double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, osculatingOrbit, propagationType);
1255         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
1256         NumericalPropagator propagator = new NumericalPropagator(integrator);
1257         propagator.setOrbitType(propagationType);
1258 
1259         NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
1260         ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
1261         propagator.addForceModel(holmesFeatherstone);
1262         propagator.setInitialState(new SpacecraftState(osculatingOrbit));
1263 
1264         return propagator.propagate(new AbsoluteDate(initialDate, 1800.));
1265     }
1266 
1267 }