1   /* Copyright 2002-2021 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         // Propagation of the initial state at t + dt
265         final double dt = 3200.;
266         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
267 
268         // Check results
269         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
270         Assert.assertEquals(state.getA(), finalState.getA(), 0.);
271         Assert.assertEquals(state.getEquinoctialEx(), finalState.getEquinoctialEx(), 0.);
272         Assert.assertEquals(state.getEquinoctialEy(), finalState.getEquinoctialEy(), 0.);
273         Assert.assertEquals(state.getHx(), finalState.getHx(), 0.);
274         Assert.assertEquals(state.getHy(), finalState.getHy(), 0.);
275         Assert.assertEquals(state.getLM() + n * dt, finalState.getLM(), 1.e-14);
276 
277     }
278 
279     @Test
280     public void testEphemeris() {
281         SpacecraftState state = getGEOState();
282         setDSSTProp(state);
283 
284         // Set ephemeris mode
285         EphemerisGenerator generator = dsstProp.getEphemerisGenerator();
286 
287         // Propagation of the initial state at t + 10 days
288         final double dt = 2. * Constants.JULIAN_DAY;
289         dsstProp.propagate(state.getDate().shiftedBy(5. * dt));
290 
291         // Get ephemeris
292         BoundedPropagator ephem = generator.getGeneratedEphemeris();
293 
294         // Propagation of the initial state with ephemeris at t + 2 days
295         final SpacecraftState s = ephem.propagate(state.getDate().shiftedBy(dt));
296 
297         // Check results
298         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
299         Assert.assertEquals(state.getA(), s.getA(), 0.);
300         Assert.assertEquals(state.getEquinoctialEx(), s.getEquinoctialEx(), 0.);
301         Assert.assertEquals(state.getEquinoctialEy(), s.getEquinoctialEy(), 0.);
302         Assert.assertEquals(state.getHx(), s.getHx(), 0.);
303         Assert.assertEquals(state.getHy(), s.getHy(), 0.);
304         Assert.assertEquals(state.getLM() + n * dt, s.getLM(), 1.5e-14);
305 
306     }
307 
308     @Test
309     public void testImpulseManeuver() {
310         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);
311         final double a = initialOrbit.getA();
312         final double e = initialOrbit.getE();
313         final double i = initialOrbit.getI();
314         final double mu = initialOrbit.getMu();
315         final double vApo = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
316         double dv = 0.99 * FastMath.tan(i) * vApo;
317 
318         // Set propagator with state
319         setDSSTProp(new SpacecraftState(initialOrbit));
320 
321         // Add impulse maneuver
322         dsstProp.setAttitudeProvider(new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS));
323         dsstProp.addEventDetector(new ImpulseManeuver<NodeDetector>(new NodeDetector(initialOrbit, FramesFactory.getEME2000()), new Vector3D(dv, Vector3D.PLUS_J), 400.0));
324         SpacecraftState propagated = dsstProp.propagate(initialOrbit.getDate().shiftedBy(8000));
325 
326         Assert.assertEquals(0.0028257, propagated.getI(), 1.0e-6);
327     }
328 
329     @Test
330     public void testPropagationWithCentralBody() throws Exception {
331 
332         // Central Body geopotential 4x4
333         final UnnormalizedSphericalHarmonicsProvider provider =
334                 GravityFieldFactory.getUnnormalizedProvider(4, 4);
335         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
336 
337         // GPS Orbit
338         final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400,
339                                                        TimeScalesFactory.getUTC());
340         final Orbit orbit = new KeplerianOrbit(26559890.,
341                                                0.0041632,
342                                                FastMath.toRadians(55.2),
343                                                FastMath.toRadians(315.4985),
344                                                FastMath.toRadians(130.7562),
345                                                FastMath.toRadians(44.2377),
346                                                PositionAngle.MEAN,
347                                                FramesFactory.getEME2000(),
348                                                initDate,
349                                                provider.getMu());
350 
351         // Set propagator with state and force model
352         setDSSTProp(new SpacecraftState(orbit));
353         dsstProp.addForceModel(new DSSTZonal(provider, 4, 3, 9));
354         dsstProp.addForceModel(new DSSTTesseral(earthFrame,
355                                                 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
356                                                 4, 4, 4, 8, 4, 4, 2));
357 
358         // 5 days propagation
359         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
360 
361         // Ref GTDS_DSST:
362         // a    = 26559.92081 km
363         // h/ey =   0.2731622444E-03
364         // k/ex =   0.4164167597E-02
365         // p/hy =  -0.3399607878
366         // q/hx =   0.3971568634
367         // lM   = 140.6375352°
368         Assert.assertEquals(26559920.81, state.getA(), 1.e-1);
369         Assert.assertEquals(0.2731622444E-03, state.getEquinoctialEx(), 2.e-8);
370         Assert.assertEquals(0.4164167597E-02, state.getEquinoctialEy(), 2.e-8);
371         Assert.assertEquals(-0.3399607878, state.getHx(), 5.e-8);
372         Assert.assertEquals(0.3971568634, state.getHy(), 2.e-6);
373         Assert.assertEquals(140.6375352,
374                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
375                             5.e-3);
376     }
377 
378     @Test
379     public void testPropagationWithThirdBody() throws IOException {
380 
381         // Central Body geopotential 2x0
382         final UnnormalizedSphericalHarmonicsProvider provider =
383                 GravityFieldFactory.getUnnormalizedProvider(2, 0);
384         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
385         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
386         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
387                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
388                                                    provider, 2, 0, 0, 2, 2, 0, 0);
389 
390         // Third Bodies Force Model (Moon + Sun)
391         DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
392         DSSTForceModel sun  = new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu());
393 
394         // SIRIUS Orbit
395         final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
396                                                        TimeScalesFactory.getUTC());
397         final Orbit orbit = new KeplerianOrbit(42163393.,
398                                                0.2684,
399                                                FastMath.toRadians(63.435),
400                                                FastMath.toRadians(270.0),
401                                                FastMath.toRadians(285.0),
402                                                FastMath.toRadians(344.0),
403                                                PositionAngle.MEAN,
404                                                FramesFactory.getEME2000(),
405                                                initDate,
406                                                provider.getMu());
407 
408         // Set propagator with state and force model
409         setDSSTProp(new SpacecraftState(orbit));
410         dsstProp.addForceModel(zonal);
411         dsstProp.addForceModel(tesseral);
412         dsstProp.addForceModel(moon);
413         dsstProp.addForceModel(sun);
414 
415         // 5 days propagation
416         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
417 
418         // Ref Standalone_DSST:
419         // a    = 42163393.0 m
420         // h/ey =  -0.06893353670734315
421         // k/ex =  -0.2592789733084587
422         // p/hy =  -0.5968524904937771
423         // q/hx =   0.1595005111738418
424         // lM   = 183°9386620425922
425         Assert.assertEquals(42163393.0, state.getA(), 1.e-1);
426         Assert.assertEquals(-0.2592789733084587, state.getEquinoctialEx(), 5.e-7);
427         Assert.assertEquals(-0.06893353670734315, state.getEquinoctialEy(), 2.e-7);
428         Assert.assertEquals( 0.1595005111738418, state.getHx(), 2.e-7);
429         Assert.assertEquals(-0.5968524904937771, state.getHy(), 5.e-8);
430         Assert.assertEquals(183.9386620425922,
431                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
432                             3.e-2);
433     }
434 
435     @Test(expected=OrekitException.class)
436     public void testTooSmallMaxDegree() {
437         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 1, 0, 3);
438     }
439 
440     @Test(expected=OrekitException.class)
441     public void testTooLargeMaxDegree() {
442         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 8, 0, 8);
443     }
444 
445     @Test(expected=OrekitException.class)
446     public void testWrongMaxPower() {
447         new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(8, 8), 4, 4, 4);
448     }
449 
450     @Test
451     public void testPropagationWithDrag() {
452 
453         // Central Body geopotential 2x0
454         final UnnormalizedSphericalHarmonicsProvider provider =
455                 GravityFieldFactory.getUnnormalizedProvider(2, 0);
456         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
457         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 0, 5);
458         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
459                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
460                                                    provider, 2, 0, 0, 2, 2, 0, 0);
461 
462         // Drag Force Model
463         final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
464                                                             Constants.WGS84_EARTH_FLATTENING,
465                                                             earthFrame);
466         final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
467 
468         final double cd = 2.0;
469         final double area = 25.0;
470         DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area, provider.getMu());
471 
472         // LEO Orbit
473         final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
474                                                        TimeScalesFactory.getUTC());
475         final Orbit orbit = new KeplerianOrbit(7204535.848109440,
476                                                0.0012402238462686,
477                                                FastMath.toRadians(98.74341600466740),
478                                                FastMath.toRadians(111.1990175076630),
479                                                FastMath.toRadians(43.32990110790340),
480                                                FastMath.toRadians(68.66852509725620),
481                                                PositionAngle.MEAN,
482                                                FramesFactory.getEME2000(),
483                                                initDate,
484                                                provider.getMu());
485 
486         // Set propagator with state and force model
487         setDSSTProp(new SpacecraftState(orbit));
488         dsstProp.addForceModel(zonal);
489         dsstProp.addForceModel(tesseral);
490         dsstProp.addForceModel(drag);
491 
492         // 5 days propagation
493         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
494 
495         // Ref Standalone_DSST:
496         // a    = 7204521.657141485 m
497         // h/ey =  0.0007093755541595772
498         // k/ex = -0.001016800430994036
499         // p/hy =  0.8698955648709271
500         // q/hx =  0.7757573478894775
501         // lM   = 193°0939742953394
502         Assert.assertEquals(7204521.657141485, state.getA(), 6.e-1);
503         Assert.assertEquals(-0.001016800430994036, state.getEquinoctialEx(), 5.e-8);
504         Assert.assertEquals(0.0007093755541595772, state.getEquinoctialEy(), 2.e-8);
505         Assert.assertEquals(0.7757573478894775, state.getHx(), 5.e-8);
506         Assert.assertEquals(0.8698955648709271, state.getHy(), 5.e-8);
507         Assert.assertEquals(193.0939742953394,
508                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
509                             2.e-3);
510         //Assert.assertEquals(((DSSTAtmosphericDrag)drag).getCd(), cd, 1e-9);
511         //Assert.assertEquals(((DSSTAtmosphericDrag)drag).getArea(), area, 1e-9);
512         Assert.assertEquals(((DSSTAtmosphericDrag)drag).getAtmosphere(), atm);
513 
514         final double atmosphericMaxConstant = 1000000.0; //DSSTAtmosphericDrag.ATMOSPHERE_ALTITUDE_MAX
515         Assert.assertEquals(((DSSTAtmosphericDrag)drag).getRbar(), atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 1e-9);
516     }
517 
518     @Test
519     public void testPropagationWithSolarRadiationPressure() {
520 
521         // Central Body geopotential 2x0
522         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
523         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
524         DSSTForceModel tesseral = new DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(),
525                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
526                                                    provider, 2, 0, 0, 2, 2, 0, 0);
527 
528         // SRP Force Model
529         DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
530                                                             Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
531                                                             provider.getMu());
532 
533         // GEO Orbit
534         final AbsoluteDate initDate = new AbsoluteDate(2003, 9, 16, 0, 0, 00.000,
535                                                        TimeScalesFactory.getUTC());
536         final Orbit orbit = new KeplerianOrbit(42166258.,
537                                                0.0001,
538                                                FastMath.toRadians(0.001),
539                                                FastMath.toRadians(315.4985),
540                                                FastMath.toRadians(130.7562),
541                                                FastMath.toRadians(44.2377),
542                                                PositionAngle.MEAN,
543                                                FramesFactory.getGCRF(),
544                                                initDate,
545                                                provider.getMu());
546 
547         // Set propagator with state and force model
548         dsstProp = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(86400.));
549         dsstProp.setInitialState(new SpacecraftState(orbit), PropagationType.MEAN);
550         dsstProp.addForceModel(zonal);
551         dsstProp.addForceModel(tesseral);
552         dsstProp.addForceModel(srp);
553 
554         // 10 days propagation
555         final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(10. * 86400.));
556 
557         // Ref Standalone_DSST:
558         // a    = 42166257.99807995 m
559         // h/ey = -0.1191876027555493D-03
560         // k/ex = -0.1781865038201885D-05
561         // p/hy =  0.6618387121369373D-05
562         // q/hx = -0.5624363171289686D-05
563         // lM   = 140°3496229467104
564         Assert.assertEquals(42166257.99807995, state.getA(), 0.8);
565         Assert.assertEquals(-0.1781865038201885e-05, state.getEquinoctialEx(), 3.e-7);
566         Assert.assertEquals(-0.1191876027555493e-03, state.getEquinoctialEy(), 4.e-6);
567         Assert.assertEquals(-0.5624363171289686e-05, state.getHx(), 4.e-9);
568         Assert.assertEquals( 0.6618387121369373e-05, state.getHy(), 3.e-10);
569         Assert.assertEquals(140.3496229467104,
570                             FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)),
571                             2.e-4);
572     }
573 
574     @Test
575     public void testStopEvent() {
576         SpacecraftState state = getLEOState();
577         setDSSTProp(state);
578 
579         final AbsoluteDate stopDate = state.getDate().shiftedBy(1000);
580         CheckingHandler<DateDetector> checking = new CheckingHandler<DateDetector>(Action.STOP);
581         dsstProp.addEventDetector(new DateDetector(stopDate).withHandler(checking));
582         checking.assertEvent(false);
583         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(3200));
584         checking.assertEvent(true);
585         Assert.assertEquals(0, finalState.getDate().durationFrom(stopDate), 1.0e-10);
586     }
587 
588     @Test
589     public void testContinueEvent() {
590         SpacecraftState state = getLEOState();
591         setDSSTProp(state);
592 
593         final AbsoluteDate resetDate = state.getDate().shiftedBy(1000);
594         CheckingHandler<DateDetector> checking = new CheckingHandler<DateDetector>(Action.CONTINUE);
595         dsstProp.addEventDetector(new DateDetector(resetDate).withHandler(checking));
596         final double dt = 3200;
597         checking.assertEvent(false);
598         final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
599         checking.assertEvent(true);
600         final double n = FastMath.sqrt(state.getMu() / state.getA()) / state.getA();
601         Assert.assertEquals(state.getA(), finalState.getA(), 1.0e-10);
602         Assert.assertEquals(state.getEquinoctialEx(), finalState.getEquinoctialEx(), 1.0e-10);
603         Assert.assertEquals(state.getEquinoctialEy(), finalState.getEquinoctialEy(), 1.0e-10);
604         Assert.assertEquals(state.getHx(), finalState.getHx(), 1.0e-10);
605         Assert.assertEquals(state.getHy(), finalState.getHy(), 1.0e-10);
606         Assert.assertEquals(state.getLM() + n * dt, finalState.getLM(), 6.0e-10);
607     }
608 
609     @Test
610     public void testIssue157() {
611         Utils.setDataRoot("regular-data:potential/icgem-format");
612         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
613         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
614         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
615                                          FramesFactory.getTOD(false),
616                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
617                                          nshp.getMu());
618         double period = orbit.getKeplerianPeriod();
619         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
620         AdaptiveStepsizeIntegrator integrator =
621                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
622         integrator.setInitialStepSize(10 * period);
623         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.MEAN);
624         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
625                                                       Constants.WGS84_EARTH_FLATTENING,
626                                                       FramesFactory.getGTOD(false));
627         CelestialBody sun = CelestialBodyFactory.getSun();
628         CelestialBody moon = CelestialBodyFactory.getMoon();
629         propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
630         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
631                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
632                                                   nshp, 8, 8, 4, 12, 8, 8, 4));
633         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
634         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
635         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
636         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
637 
638 
639         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.OSCULATING);
640         SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
641         // the following comparison is in fact meaningless
642         // the initial orbit is osculating the final orbit is a mean orbit
643         // and they are not considered at the same epoch
644         // we keep it only as is was an historical test
645         Assert.assertEquals(2189.4, orbit.getA() - finalState.getA(), 1.0);
646 
647         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
648         finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
649         // the following comparison is realistic
650         // both the initial orbit and final orbit are mean orbits
651         Assert.assertEquals(1478.05, orbit.getA() - finalState.getA(), 1.0);
652 
653     }
654 
655     @Test
656     public void testEphemerisGeneration() {
657         Utils.setDataRoot("regular-data:potential/icgem-format");
658         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
659         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
660         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
661                                          FramesFactory.getTOD(false),
662                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
663                                          nshp.getMu());
664         double period = orbit.getKeplerianPeriod();
665         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
666         AdaptiveStepsizeIntegrator integrator =
667                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
668         integrator.setInitialStepSize(10 * period);
669         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
670         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
671                                                       Constants.WGS84_EARTH_FLATTENING,
672                                                       FramesFactory.getGTOD(false));
673         CelestialBody sun = CelestialBodyFactory.getSun();
674         CelestialBody moon = CelestialBodyFactory.getMoon();
675         propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
676         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
677                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
678                                                   nshp, 8, 8, 4, 12, 8, 8, 4));
679         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
680         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
681         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
682         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
683         propagator.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
684 
685         // direct generation of states
686         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
687         final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
688         propagator.setStepHandler(600, currentState -> states.add(currentState));
689         propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
690 
691         // ephemeris generation
692         propagator.setInitialState(new SpacecraftState(orbit, 45.0), PropagationType.MEAN);
693         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
694         propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
695         BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
696 
697         double maxError = 0;
698         for (final SpacecraftState state : states) {
699             final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
700             final double error = Vector3D.distance(state.getPVCoordinates().getPosition(),
701                                                    fromEphemeris.getPVCoordinates().getPosition());
702             maxError = FastMath.max(maxError, error);
703         }
704         Assert.assertEquals(0.0, maxError, 1.0e-10);
705     }
706 
707     @Test
708     public void testGetInitialOsculatingState() throws IllegalArgumentException, OrekitException {
709         final SpacecraftState initialState = getGEOState();
710 
711         // build integrator
712         final double minStep = initialState.getKeplerianPeriod() * 0.1;
713         final double maxStep = initialState.getKeplerianPeriod() * 10.0;
714         final double[][] tol = DSSTPropagator.tolerances(0.1, initialState.getOrbit());
715         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
716 
717         // build the propagator for the propagation of the mean elements
718         DSSTPropagator prop = new DSSTPropagator(integrator, PropagationType.MEAN);
719 
720         final UnnormalizedSphericalHarmonicsProvider provider =
721                 GravityFieldFactory.getUnnormalizedProvider(4, 0);
722         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
723         DSSTForceModel zonal    = new DSSTZonal(provider, 4, 3, 9);
724         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
725                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
726                                                    provider, 4, 0, 4, 8, 4, 0, 2);
727         prop.addForceModel(zonal);
728         prop.addForceModel(tesseral);
729 
730         // Set the initial state
731         prop.setInitialState(initialState, PropagationType.MEAN);
732         // Check the stored initial state is the osculating one
733         Assert.assertEquals(initialState, prop.getInitialState());
734         // Check that no propagation, i.e. propagation to the initial date, provides the initial
735         // osculating state although the propagator is configured to propagate mean elements !!!
736         Assert.assertEquals(initialState, prop.propagate(initialState.getDate()));
737     }
738 
739     @Test
740     public void testMeanToOsculatingState() throws IllegalArgumentException, OrekitException {
741         final SpacecraftState meanState = getGEOState();
742 
743         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
744         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
745         final DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
746         final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
747                                                          Constants.WGS84_EARTH_ANGULAR_VELOCITY,
748                                                          provider, 2, 0, 0, 2, 2, 0, 0);
749 
750         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
751         forces.add(zonal);
752         forces.add(tesseral);
753 
754         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
755         Assert.assertEquals(1559.1,
756                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
757                                               osculatingState.getPVCoordinates().getPosition()),
758                             1.0);
759     }
760 
761     @Test
762     public void testOsculatingToMeanState() throws IllegalArgumentException, OrekitException {
763         final SpacecraftState meanState = getGEOState();
764 
765         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
766         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
767 
768         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
769         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
770                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
771                                                    provider, 2, 0, 0, 2, 2, 0, 0);
772 
773         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
774         forces.add(zonal);
775         forces.add(tesseral);
776 
777         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
778 
779         // there are no Gaussian force models, we don't need an attitude provider
780         final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
781 
782         Assert.assertEquals(meanState.getA(), computedMeanState.getA(), 2.0e-8);
783         Assert.assertEquals(0.0,
784                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
785                                              computedMeanState.getPVCoordinates().getPosition()),
786                             2.0e-8);
787     }
788 
789     @Test
790     public void testShortPeriodCoefficients() {
791         Utils.setDataRoot("regular-data:potential/icgem-format");
792         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
793         UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(4, 4);
794         Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN,
795                                          FramesFactory.getTOD(false),
796                                          new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
797                                          nshp.getMu());
798         double period = orbit.getKeplerianPeriod();
799         double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
800         AdaptiveStepsizeIntegrator integrator =
801                 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
802         integrator.setInitialStepSize(10 * period);
803         DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
804         OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
805                                                       Constants.WGS84_EARTH_FLATTENING,
806                                                       FramesFactory.getGTOD(false));
807         CelestialBody sun = CelestialBodyFactory.getSun();
808         CelestialBody moon = CelestialBodyFactory.getMoon();
809         propagator.addForceModel(new DSSTZonal(nshp, 4, 3, 9));
810         propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
811                                                   Constants.WGS84_EARTH_ANGULAR_VELOCITY,
812                                                   nshp, 4, 4, 4, 8, 4, 4, 2));
813         propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
814         propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
815         propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
816         propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius(), nshp.getMu()));
817 
818         final AbsoluteDate finalDate = orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY);
819         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
820         final SpacecraftState stateNoConfig = propagator.propagate(finalDate);
821         Assert.assertEquals(0, stateNoConfig.getAdditionalStates().size());
822 
823         propagator.setSelectedCoefficients(new HashSet<String>());
824         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
825         final SpacecraftState stateConfigEmpty = propagator.propagate(finalDate);
826         Assert.assertEquals(234, stateConfigEmpty.getAdditionalStates().size());
827 
828         final Set<String> selected = new HashSet<String>();
829         selected.add("DSST-3rd-body-Moon-s[7]");
830         selected.add("DSST-central-body-tesseral-c[-2][3]");
831         propagator.setSelectedCoefficients(selected);
832         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
833         final SpacecraftState stateConfigeSelected = propagator.propagate(finalDate);
834         Assert.assertEquals(selected.size(), stateConfigeSelected.getAdditionalStates().size());
835 
836         propagator.setSelectedCoefficients(null);
837         propagator.resetInitialState(new SpacecraftState(orbit, 45.0));
838         final SpacecraftState stateConfigNull = propagator.propagate(finalDate);
839         Assert.assertEquals(0, stateConfigNull.getAdditionalStates().size());
840 
841     }
842 
843     @Test
844     public void testIssueMeanInclination() {
845 
846         final double earthAe = 6378137.0;
847         final double earthMu = 3.9860044E14;
848         final double earthJ2 = 0.0010826;
849 
850         // Initialize the DSST propagator with only J2 perturbation
851         Orbit orb = new KeplerianOrbit(new TimeStampedPVCoordinates(new AbsoluteDate("1992-10-08T15:20:38.821",
852                                                                                      TimeScalesFactory.getUTC()),
853                                                                     new Vector3D(5392808.809823, -4187618.3357927715, -44206.638015847195),
854                                                                     new Vector3D(2337.4472786270794, 2474.0146611860464, 6778.507766114648)),
855                                        FramesFactory.getTOD(false), earthMu);
856         final SpacecraftState ss = new SpacecraftState(orb);
857         final UnnormalizedSphericalHarmonicsProvider provider =
858               GravityFieldFactory.getUnnormalizedProvider(earthAe, earthMu, TideSystem.UNKNOWN,
859                                                           new double[][] { { 0.0 }, { 0.0 }, { -earthJ2 } },
860                                                           new double[][] { { 0.0 }, { 0.0 }, { 0.0 } });
861         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
862         DSSTForceModel zonal    = new DSSTZonal(provider, 2, 1, 5);
863         DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
864                                                    Constants.WGS84_EARTH_ANGULAR_VELOCITY,
865                                                    provider, 2, 0, 0, 2, 2, 0, 0);
866         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
867         forces.add(zonal);
868         forces.add(tesseral);
869         // Computes J2 mean elements using the DSST osculating to mean converter
870         final Orbit meanOrb = DSSTPropagator.computeMeanState(ss, null, forces).getOrbit();
871         Assert.assertEquals(0.0164196, FastMath.toDegrees(orb.getI() - meanOrb.getI()), 1.0e-7);
872     }
873 
874     @Test
875     public void testIssue257() {
876         final SpacecraftState meanState = getGEOState();
877 
878         // Third Bodies Force Model (Moon + Sun)
879         final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), meanState.getMu());
880         final DSSTForceModel sun  = new DSSTThirdBody(CelestialBodyFactory.getSun(), meanState.getMu());
881 
882         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
883         forces.add(moon);
884         forces.add(sun);
885 
886         final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
887         Assert.assertEquals(734.3,
888                             Vector3D.distance(meanState.getPVCoordinates().getPosition(),
889                                               osculatingState.getPVCoordinates().getPosition()),
890                             1.0);
891 
892         final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
893         Assert.assertEquals(734.3,
894                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
895                                               computedMeanState.getPVCoordinates().getPosition()),
896                             1.0);
897 
898         Assert.assertEquals(0.0,
899                             Vector3D.distance(computedMeanState.getPVCoordinates().getPosition(),
900                                               meanState.getPVCoordinates().getPosition()),
901                             5.0e-6);
902 
903     }
904 
905     @Test
906     public void testIssue339() {
907 
908         final SpacecraftState osculatingState = getLEOState();
909 
910         final CelestialBody    sun   = CelestialBodyFactory.getSun();
911         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
912                                                             Constants.WGS84_EARTH_FLATTENING,
913                                                             FramesFactory.getITRF(IERSConventions.IERS_2010,
914                                                                                   true));
915         final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
916                                                                                      sun,
917                                                                                      50.0, Vector3D.PLUS_J,
918                                                                                      2.0, 0.1,
919                                                                                      0.2, 0.6);
920         final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
921         final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(),
922                                                                 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
923                                                                 0.0, 0.0, 0.0);
924 
925         // Surface force models that require an attitude provider
926         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
927         forces.add(new DSSTSolarRadiationPressure(sun,
928                                                   Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
929                                                   boxAndWing,
930                                                   osculatingState.getMu()));
931         forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getMu()));
932 
933         final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
934         Assert.assertEquals(0.522,
935                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
936                                               meanState.getPVCoordinates().getPosition()),
937                             0.001);
938 
939         final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
940         Assert.assertEquals(0.0,
941                             Vector3D.distance(osculatingState.getPVCoordinates().getPosition(),
942                                               computedOsculatingState.getPVCoordinates().getPosition()),
943                             5.0e-6);
944 
945     }
946 
947     @Test
948     public void testIssue613() {
949         // Spacecraft state
950         final SpacecraftState state = getLEOState();
951 
952         // Body frame
953         final Frame itrf = FramesFactory .getITRF(IERSConventions.IERS_2010, true);
954         
955         // Earth
956         final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 4);
957         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
958                                                             Constants.WGS84_EARTH_FLATTENING, itrf);
959         // Detectors
960         final List<EventDetector> events = new ArrayList<>();
961         events.add(new AltitudeDetector(85.5, earth));
962         events.add(new LatitudeCrossingDetector(earth, 0.0));
963 
964         // Force models
965         final List<DSSTForceModel> forceModels = new ArrayList<>();
966         forceModels.add(new DSSTZonal(provider));
967         forceModels.add(new DSSTTesseral(itrf, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
968         forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu()));
969         forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu()));
970 
971         // Set up DSST propagator
972         final double[][] tol = DSSTPropagator.tolerances(10.0, state.getOrbit());
973         final ODEIntegrator integrator = new DormandPrince54Integrator(60.0, 3600.0, tol[0], tol[1]);
974         final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
975         for (DSSTForceModel force : forceModels) {
976             propagator.addForceModel(force);
977         }
978         for (EventDetector event : events) {
979             propagator.addEventDetector(event);
980         }
981         propagator.setInitialState(state);
982 
983         // Propagation
984         final SpacecraftState finalState = propagator.propagate(state.getDate().shiftedBy(86400.0));
985 
986         // Verify is the propagation is correctly performed
987         Assert.assertEquals(finalState.getMu(), 3.986004415E14, Double.MIN_VALUE);
988     }
989 
990     @Test
991     public void testIssue339WithAccelerations() {
992         final SpacecraftState osculatingState = getLEOStatePropagatedBy30Minutes();
993         final CelestialBody sun = CelestialBodyFactory.getSun();
994         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
995                                                             FramesFactory.getITRF(IERSConventions.IERS_2010, true));
996         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);
997         final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
998         final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(), LOFType.LVLH_CCSDS, RotationOrder.XYZ, 0.0, 0.0, 0.0);
999         // Surface force models that require an attitude provider
1000         final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
1001         forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getMu()));
1002         final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1003         final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1004         Assert.assertEquals(0.0, Vector3D.distance(osculatingState.getPVCoordinates().getPosition(), computedOsculatingState.getPVCoordinates().getPosition()),
1005                             5.0e-6);
1006     }
1007 
1008     @Test
1009     public void testIssue704() {
1010 
1011         // Coordinates
1012         final Orbit         orbit = getLEOState().getOrbit();
1013         final PVCoordinates pv    = orbit.getPVCoordinates();
1014 
1015         // dP
1016         final double dP = 10.0;
1017 
1018         // Computes dV
1019         final double r2 = pv.getPosition().getNormSq();
1020         final double v  = pv.getVelocity().getNorm();
1021         final double dV = orbit.getMu() * dP / (v * r2);
1022 
1023         // Verify
1024         final double[][] tol1 = DSSTPropagator.tolerances(dP, orbit);
1025         final double[][] tol2 = DSSTPropagator.tolerances(dP, dV, orbit);
1026         for (int i = 0; i < tol1.length; i++) {
1027             Assert.assertArrayEquals(tol1[i], tol2[i], Double.MIN_VALUE);
1028         }
1029 
1030     }
1031 
1032     /** This test is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1033     @Test
1034     public void testIssue670() {
1035 
1036         final NumericalForce force     = new NumericalForce();
1037         final DSSTForce      dsstForce = new DSSTForce(force, Constants.WGS84_EARTH_MU);
1038 
1039         SpacecraftState state = getLEOState();
1040         setDSSTProp(state);
1041         dsstProp.addForceModel(dsstForce);
1042 
1043         // Verify flag are false
1044         Assert.assertFalse(force.initialized);
1045         Assert.assertFalse(force.accComputed);
1046 
1047         // Propagation of the initial state at t + dt
1048         final double dt = 3200.;
1049         final AbsoluteDate target = state.getDate().shiftedBy(dt);
1050 
1051         dsstProp.propagate(target);
1052 
1053         // Flag must be true
1054         Assert.assertTrue(force.initialized);
1055         Assert.assertTrue(force.accComputed);
1056 
1057     }
1058 
1059 
1060     private SpacecraftState getGEOState() throws IllegalArgumentException, OrekitException {
1061         // No shadow at this date
1062         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
1063                                                        TimeScalesFactory.getUTC());
1064         final Orbit orbit = new EquinoctialOrbit(42164000,
1065                                                  10e-3,
1066                                                  10e-3,
1067                                                  FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
1068                                                  FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
1069                                                  PositionAngle.TRUE,
1070                                                  FramesFactory.getEME2000(),
1071                                                  initDate,
1072                                                  3.986004415E14);
1073         return new SpacecraftState(orbit);
1074     }
1075 
1076     private SpacecraftState getLEOState() throws IllegalArgumentException, OrekitException {
1077         final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1078         final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1079         // Spring equinoxe 21st mars 2003 1h00m
1080         final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1081         return new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
1082                                                         FramesFactory.getEME2000(),
1083                                                         initDate,
1084                                                         3.986004415E14));
1085     }
1086 
1087     private void setDSSTProp(SpacecraftState initialState) {
1088         initialState.getDate();
1089         final double minStep = initialState.getKeplerianPeriod();
1090         final double maxStep = 100. * minStep;
1091         final double[][] tol = DSSTPropagator.tolerances(1.0, initialState.getOrbit());
1092         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1093         dsstProp = new DSSTPropagator(integrator);
1094         dsstProp.setInitialState(initialState, PropagationType.MEAN);
1095     }
1096 
1097     private static class CheckingHandler<T extends EventDetector> implements EventHandler<T> {
1098 
1099         private final Action actionOnEvent;
1100         private boolean gotHere;
1101 
1102         public CheckingHandler(final Action actionOnEvent) {
1103             this.actionOnEvent = actionOnEvent;
1104             this.gotHere       = false;
1105         }
1106 
1107         public void assertEvent(boolean expected) {
1108             Assert.assertEquals(expected, gotHere);
1109         }
1110 
1111         @Override
1112         public Action eventOccurred(SpacecraftState s, T detector, boolean increasing) {
1113             gotHere = true;
1114             return actionOnEvent;
1115         }
1116 
1117     }
1118 
1119     /** This class is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1120     private class DSSTForce extends AbstractGaussianContribution {
1121 
1122         DSSTForce(ForceModel contribution, double mu) {
1123             super("DSST mock -", 6.0e-10, contribution, mu);
1124         }
1125 
1126         /** {@inheritDoc} */
1127         @Override
1128         public EventDetector[] getEventsDetectors() {
1129             return null;
1130         }
1131 
1132         /** {@inheritDoc} */
1133         @Override
1134         public <T extends CalculusFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
1135             return null;
1136         }
1137 
1138         /** {@inheritDoc} */
1139         @Override
1140         protected List<ParameterDriver> getParametersDriversWithoutMu() {
1141             return Collections.emptyList();
1142         }
1143 
1144         /** {@inheritDoc} */
1145         @Override
1146         protected double[] getLLimits(SpacecraftState state,
1147                                       AuxiliaryElements auxiliaryElements) {
1148             return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0),
1149                                    FastMath.PI + MathUtils.normalizeAngle(state.getLv(), 0) };
1150         }
1151 
1152         /** {@inheritDoc} */
1153         @Override
1154         protected <T extends CalculusFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
1155                                                                  FieldAuxiliaryElements<T> auxiliaryElements) {
1156             final Field<T> field = state.getDate().getField();
1157             final T zero = field.getZero();
1158             final T[] tab = MathArrays.buildArray(field, 2);
1159             tab[0] = MathUtils.normalizeAngle(state.getLv(), zero).subtract(FastMath.PI);
1160             tab[1] = MathUtils.normalizeAngle(state.getLv(), zero).add(FastMath.PI);
1161             return tab;
1162         }
1163         
1164     }
1165 
1166     /** This class is based on the example given by Orekit user kris06 in https://gitlab.orekit.org/orekit/orekit/-/issues/670. */
1167     private class NumericalForce extends AbstractForceModel {
1168 
1169         private boolean initialized;
1170         private boolean accComputed;
1171 
1172         NumericalForce() {
1173             this.initialized = false;
1174         }
1175 
1176         /** {@inheritDoc} */
1177         @Override
1178         public void init(final SpacecraftState s0, final AbsoluteDate t) {
1179             this.initialized = true;
1180             this.accComputed = false;
1181         }
1182 
1183         /** {@inheritDoc} */
1184         @Override
1185         public boolean dependsOnPositionOnly() {
1186             return false;
1187         }
1188 
1189         /** {@inheritDoc} */
1190         @Override
1191         public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1192             this.accComputed = true;
1193             return Vector3D.ZERO;
1194         }
1195 
1196         /** {@inheritDoc} */
1197         @Override
1198         public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1199             return FieldVector3D.getZero(s.getDate().getField());
1200         }
1201 
1202         /** {@inheritDoc} */
1203         @Override
1204         public Stream<EventDetector> getEventsDetectors() {
1205             return Stream.empty();
1206         }
1207 
1208         /** {@inheritDoc} */
1209         @Override
1210         public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
1211             return Stream.empty();
1212         }
1213 
1214 
1215         @Override
1216         public List<ParameterDriver> getParametersDrivers() {
1217             return Collections.emptyList();
1218         }
1219         
1220     }
1221 
1222     @Before
1223     public void setUp() throws IOException, ParseException {
1224         Utils.setDataRoot("regular-data:potential/shm-format");
1225     }
1226 
1227     @After
1228     public void tearDown() {
1229         dsstProp = null;
1230     }
1231 
1232     private SpacecraftState getLEOStatePropagatedBy30Minutes() throws IllegalArgumentException, OrekitException {
1233 
1234         final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1235         final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1236         // Spring equinoxe 21st mars 2003 1h00m
1237         final AbsoluteDate initialDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1238         final CartesianOrbit osculatingOrbit = new CartesianOrbit(new PVCoordinates(position, velocity), FramesFactory.getTOD(IERSConventions.IERS_1996, false),
1239                                                                   initialDate, Constants.WGS84_EARTH_MU);
1240         // Adaptive step integrator
1241         // with a minimum step of 0.001 and a maximum step of 1000
1242         double minStep = 0.001;
1243         double maxstep = 1000.0;
1244         double positionTolerance = 10.0;
1245         OrbitType propagationType = OrbitType.EQUINOCTIAL;
1246         double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, osculatingOrbit, propagationType);
1247         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
1248         NumericalPropagator propagator = new NumericalPropagator(integrator);
1249         propagator.setOrbitType(propagationType);
1250 
1251         NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
1252         ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
1253         propagator.addForceModel(holmesFeatherstone);
1254         propagator.setInitialState(new SpacecraftState(osculatingOrbit));
1255 
1256         return propagator.propagate(new AbsoluteDate(initialDate, 1800.));
1257     }
1258 
1259 }