1   /* Copyright 2002-2025 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.conversion;
18  
19  import static org.orekit.propagation.conversion.AbstractPropagatorBuilderTest.assertPropagatorBuilderIsACopy;
20  
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
23  import org.junit.jupiter.api.*;
24  import org.orekit.Utils;
25  import org.orekit.bodies.CelestialBodyFactory;
26  import org.orekit.errors.OrekitException;
27  import org.orekit.errors.OrekitMessages;
28  import org.orekit.forces.gravity.potential.GravityFieldFactory;
29  import org.orekit.frames.Frame;
30  import org.orekit.frames.FramesFactory;
31  import org.orekit.orbits.CartesianOrbit;
32  import org.orekit.orbits.EquinoctialOrbit;
33  import org.orekit.orbits.Orbit;
34  import org.orekit.orbits.OrbitType;
35  import org.orekit.orbits.PositionAngleType;
36  import org.orekit.propagation.CartesianToleranceProvider;
37  import org.orekit.propagation.PropagationType;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.ToleranceProvider;
40  import org.orekit.propagation.integration.AdditionalDerivativesProvider;
41  import org.orekit.propagation.integration.CombinedDerivatives;
42  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
43  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
44  import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
45  import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
46  import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
47  import org.orekit.time.AbsoluteDate;
48  import org.orekit.time.TimeScalesFactory;
49  import org.orekit.utils.Constants;
50  import org.orekit.utils.PVCoordinates;
51  import org.orekit.utils.ParameterDriver;
52  
53  import java.io.IOException;
54  import java.text.ParseException;
55  import java.util.List;
56  
57  public class DSSTPropagatorBuilderTest {
58  
59      private static final double eps  = 2.0e-10;
60  
61      private double minStep;
62      private double maxStep;
63      private ToleranceProvider toleranceProvider;
64      private double[][] tolerance;
65  
66      private AbsoluteDate initDate;
67      private EquinoctialOrbit orbit;
68      private DSSTPropagator propagator;
69      private DSSTForceModel moon;
70      private DSSTForceModel sun;
71  
72      @Test
73      public void testIntegrators01() {
74  
75          ODEIntegratorBuilder abBuilder = new AdamsBashforthIntegratorBuilder(2, minStep, maxStep, toleranceProvider);
76          doTestBuildPropagator(abBuilder);
77      }
78  
79      @Test
80      public void testIntegrators02() {
81  
82          ODEIntegratorBuilder amBuilder = new AdamsMoultonIntegratorBuilder(2, minStep, maxStep, toleranceProvider);
83          doTestBuildPropagator(amBuilder);
84      }
85  
86      @Test
87      public void testIntegrators03() {
88  
89          final double stepSize = 100.;
90  
91          ODEIntegratorBuilder crkBuilder = new ClassicalRungeKuttaIntegratorBuilder(stepSize);
92          doTestBuildPropagator(crkBuilder);
93      }
94  
95      @Test
96      public void testIntegrators04() {
97  
98          final double stepSize = 100.;
99  
100         ODEIntegratorBuilder lBuilder = new LutherIntegratorBuilder(stepSize);
101         doTestBuildPropagator(lBuilder);
102     }
103 
104     @Test
105     public void testIntegrators05() {
106 
107         ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
108         doTestBuildPropagator(dp54Builder);
109     }
110 
111     @Test
112     public void testIntegrators06() {
113 
114         final double stepSize = 100.;
115 
116         ODEIntegratorBuilder eBuilder = new EulerIntegratorBuilder(stepSize);
117         doTestBuildPropagator(eBuilder);
118     }
119 
120     @Test
121     public void testIntegrators07() {
122 
123         final double stepSize = 100.;
124 
125         ODEIntegratorBuilder gBuilder = new GillIntegratorBuilder(stepSize);
126         doTestBuildPropagator(gBuilder);
127     }
128 
129     @Test
130     public void testIntegrators08() {
131 
132         ODEIntegratorBuilder gbsBuilder = new GraggBulirschStoerIntegratorBuilder(minStep, maxStep, toleranceProvider);
133         doTestBuildPropagator(gbsBuilder);
134     }
135 
136     @Test
137     public void testIntegrators09() {
138 
139         ODEIntegratorBuilder hh54Builder = new HighamHall54IntegratorBuilder(minStep, maxStep, toleranceProvider);
140         doTestBuildPropagator(hh54Builder);
141     }
142 
143     @Test
144     public void testIntegrators10() {
145 
146         final double stepSize = 100.;
147 
148         ODEIntegratorBuilder mBuilder = new MidpointIntegratorBuilder(stepSize);
149         doTestBuildPropagator(mBuilder);
150     }
151 
152     @Test
153     public void testIntegrators11() {
154 
155         final double stepSize = 100.;
156 
157         ODEIntegratorBuilder teBuilder = new ThreeEighthesIntegratorBuilder(stepSize);
158         doTestBuildPropagator(teBuilder);
159     }
160 
161     @Test
162     void testClone() {
163 
164         // Given
165         final ODEIntegratorBuilder integratorBuilder = new ClassicalRungeKuttaIntegratorBuilder(3600.0);
166         final Orbit orbit = new CartesianOrbit(new PVCoordinates(
167                 new Vector3D(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS + 400000, 0, 0),
168                 new Vector3D(0, 7668.6, 0)), FramesFactory.getGCRF(),
169                 new AbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
170 
171         final DSSTPropagatorBuilder builder =
172                 new DSSTPropagatorBuilder(orbit, integratorBuilder, 1.0, PropagationType.OSCULATING, PropagationType.OSCULATING);
173 
174         builder.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0)));
175 
176         // When
177         final DSSTPropagatorBuilder copyBuilder = (DSSTPropagatorBuilder) builder.clone();
178 
179         // Then
180         assertDSSTPropagatorBuilderIsACopy(builder, copyBuilder);
181 
182     }
183 
184     private void assertDSSTPropagatorBuilderIsACopy(final DSSTPropagatorBuilder expected,
185                                                     final DSSTPropagatorBuilder actual) {
186         assertPropagatorBuilderIsACopy(expected, actual);
187 
188         Assertions.assertEquals(expected.getIntegratorBuilder().getClass(), actual.getIntegratorBuilder().getClass());
189         Assertions.assertEquals(expected.getPropagationType(), actual.getPropagationType());
190         Assertions.assertEquals(expected.getStateType(), actual.getStateType());
191         Assertions.assertEquals(expected.getMass(), actual.getMass());
192         final List<DSSTForceModel> expectedForceModelList = expected.getAllForceModels();
193         final List<DSSTForceModel> actualForceModelList   = actual.getAllForceModels();
194         Assertions.assertEquals(expectedForceModelList.size(), actualForceModelList.size());
195         for (int i = 0; i < expectedForceModelList.size(); i++) {
196             Assertions.assertEquals(expectedForceModelList.get(i).getClass(), actualForceModelList.get(i).getClass());
197         }
198 
199     }
200 
201     private void doTestBuildPropagator(final ODEIntegratorBuilder foiBuilder) {
202 
203         // We propagate using directly the propagator of the set up
204         final Orbit orbitWithPropagator = propagator.propagate(initDate.shiftedBy(600)).getOrbit();
205 
206         // We propagate using a build version of the propagator
207         // We shall have the same results than before
208         DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
209                                                                   foiBuilder,
210                                                                   1.0,
211                                                                   PropagationType.MEAN,
212                                                                   PropagationType.MEAN);
213 
214         builder.addForceModel(moon);
215         builder.setMass(1000.);
216 
217         final DSSTPropagator prop = (DSSTPropagator) builder.buildPropagator();
218 
219         final Orbit orbitWithBuilder = prop.propagate(initDate.shiftedBy(600)).getOrbit();
220 
221         // Verify
222         Assertions.assertEquals(orbitWithPropagator.getA(),             orbitWithBuilder.getA(), 1.e-1);
223         Assertions.assertEquals(orbitWithPropagator.getEquinoctialEx(), orbitWithBuilder.getEquinoctialEx(), eps);
224         Assertions.assertEquals(orbitWithPropagator.getEquinoctialEy(), orbitWithBuilder.getEquinoctialEy(), eps);
225         Assertions.assertEquals(orbitWithPropagator.getHx(),            orbitWithBuilder.getHx(), eps);
226         Assertions.assertEquals(orbitWithPropagator.getHy(),            orbitWithBuilder.getHy(), eps);
227         Assertions.assertEquals(orbitWithPropagator.getLM(),            orbitWithBuilder.getLM(), 8.0e-10);
228 
229     }
230 
231     @Test
232     public void testIssue598() {
233         // Integrator builder
234         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
235         // Propagator builder
236         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
237                                                                   dp54Builder,
238                                                                   1.0,
239                                                                   PropagationType.MEAN,
240                                                                   PropagationType.MEAN);
241         builder.addForceModel(moon);
242         // Verify that there is no Newtonian attraction force model
243         Assertions.assertFalse(hasNewtonianAttraction(builder.getAllForceModels()));
244         // Build the DSST propagator (not used here)
245         builder.buildPropagator();
246         // Verify the addition of the Newtonian attraction force model
247         Assertions.assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
248         // Add a new force model to ensure the Newtonian attraction stay at the last position
249         builder.addForceModel(sun);
250         Assertions.assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
251     }
252 
253     @Test
254     public void testAdditionalEquations() {
255         // Integrator builder
256         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
257         // Propagator builder
258         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
259                                                                   dp54Builder,
260                                                                   1.0,
261                                                                   PropagationType.MEAN,
262                                                                   PropagationType.MEAN);
263         builder.addForceModel(moon);
264         builder.addForceModel(sun);
265 
266         // Add additional equations
267         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
268 
269             public String getName() {
270                 return "linear";
271             }
272 
273             public int getDimension() {
274                 return 1;
275             }
276 
277             public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
278                 return new CombinedDerivatives(new double[] { 1.0 }, null);
279             }
280 
281         });
282 
283         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
284 
285             public String getName() {
286                 return "linear";
287             }
288 
289             public int getDimension() {
290                 return 1;
291             }
292 
293             public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
294                 return new CombinedDerivatives(new double[] { 1.0 }, null);
295             }
296 
297         });
298 
299         try {
300             // Build the propagator
301             builder.buildPropagator();
302             Assertions.fail("an exception should have been thrown");
303         } catch (OrekitException oe) {
304             Assertions.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
305         }
306     }
307 
308     @Test
309     public void testDeselectOrbitals() {
310         // Integrator builder
311         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
312         // Propagator builder
313         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
314                                                                   dp54Builder,
315                                                                   1.0,
316                                                                   PropagationType.MEAN,
317                                                                   PropagationType.MEAN);
318         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
319             Assertions.assertTrue(driver.isSelected());
320         }
321         builder.deselectDynamicParameters();
322         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
323             Assertions.assertFalse(driver.isSelected());
324         }
325     }
326 
327     @BeforeEach
328     public void setUp() throws IOException, ParseException {
329 
330         Utils.setDataRoot("regular-data:potential");
331 
332         minStep = 1.0;
333         maxStep = 600.0;
334         toleranceProvider = ToleranceProvider.of(CartesianToleranceProvider.of(10., 0.0001, 
335                 CartesianToleranceProvider.DEFAULT_ABSOLUTE_MASS_TOLERANCE));
336 
337         final Frame earthFrame = FramesFactory.getEME2000();
338         initDate = new AbsoluteDate(2003, 07, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
339 
340         final double mu = 3.986004415E14;
341         // a    = 42163393.0 m
342         // ex =  -0.25925449177598586
343         // ey =  -0.06946703170551687
344         // hx =   0.15995912655021305
345         // hy =  -0.5969755874197339
346         // lM   = 15.47576793123677 rad
347         orbit = new EquinoctialOrbit(4.2163393E7,
348                                      -0.25925449177598586,
349                                      -0.06946703170551687,
350                                      0.15995912655021305,
351                                      -0.5969755874197339,
352                                      15.47576793123677,
353                                      PositionAngleType.MEAN,
354                                      earthFrame,
355                                      initDate,
356                                      mu);
357 
358         moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), mu);
359         sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), mu);
360 
361         tolerance  = toleranceProvider.getTolerances(orbit, OrbitType.EQUINOCTIAL);
362         propagator = new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]));
363         propagator.setInitialState(new SpacecraftState(orbit, 1000.), PropagationType.MEAN);
364         propagator.addForceModel(moon);
365 
366     }
367 
368     private boolean hasNewtonianAttraction(final List<DSSTForceModel> forceModels) {
369         final int last = forceModels.size() - 1;
370         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
371     }
372 
373 }