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  import org.orekit.utils.TimeStampedPVCoordinates;
53  
54  import java.io.IOException;
55  import java.text.ParseException;
56  import java.util.List;
57  
58  public class DSSTPropagatorBuilderTest {
59  
60      private static final double eps  = 2.0e-10;
61  
62      private double minStep;
63      private double maxStep;
64      private ToleranceProvider toleranceProvider;
65      private double[][] tolerance;
66  
67      private AbsoluteDate initDate;
68      private EquinoctialOrbit orbit;
69      private DSSTPropagator propagator;
70      private DSSTForceModel moon;
71      private DSSTForceModel sun;
72  
73      @Test
74      public void testIntegrators01() {
75  
76          ODEIntegratorBuilder abBuilder = new AdamsBashforthIntegratorBuilder(2, minStep, maxStep, toleranceProvider);
77          doTestBuildPropagator(abBuilder);
78      }
79  
80      @Test
81      public void testIntegrators02() {
82  
83          ODEIntegratorBuilder amBuilder = new AdamsMoultonIntegratorBuilder(2, minStep, maxStep, toleranceProvider);
84          doTestBuildPropagator(amBuilder);
85      }
86  
87      @Test
88      public void testIntegrators03() {
89  
90          final double stepSize = 100.;
91  
92          ODEIntegratorBuilder crkBuilder = new ClassicalRungeKuttaIntegratorBuilder(stepSize);
93          doTestBuildPropagator(crkBuilder);
94      }
95  
96      @Test
97      public void testIntegrators04() {
98  
99          final double stepSize = 100.;
100 
101         ODEIntegratorBuilder lBuilder = new LutherIntegratorBuilder(stepSize);
102         doTestBuildPropagator(lBuilder);
103     }
104 
105     @Test
106     public void testIntegrators05() {
107 
108         ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
109         doTestBuildPropagator(dp54Builder);
110     }
111 
112     @Test
113     public void testIntegrators06() {
114 
115         final double stepSize = 100.;
116 
117         ODEIntegratorBuilder eBuilder = new EulerIntegratorBuilder(stepSize);
118         doTestBuildPropagator(eBuilder);
119     }
120 
121     @Test
122     public void testIntegrators07() {
123 
124         final double stepSize = 100.;
125 
126         ODEIntegratorBuilder gBuilder = new GillIntegratorBuilder(stepSize);
127         doTestBuildPropagator(gBuilder);
128     }
129 
130     @Test
131     public void testIntegrators08() {
132 
133         ODEIntegratorBuilder gbsBuilder = new GraggBulirschStoerIntegratorBuilder(minStep, maxStep, toleranceProvider);
134         doTestBuildPropagator(gbsBuilder);
135     }
136 
137     @Test
138     public void testIntegrators09() {
139 
140         ODEIntegratorBuilder hh54Builder = new HighamHall54IntegratorBuilder(minStep, maxStep, toleranceProvider);
141         doTestBuildPropagator(hh54Builder);
142     }
143 
144     @Test
145     public void testIntegrators10() {
146 
147         final double stepSize = 100.;
148 
149         ODEIntegratorBuilder mBuilder = new MidpointIntegratorBuilder(stepSize);
150         doTestBuildPropagator(mBuilder);
151     }
152 
153     @Test
154     public void testIntegrators11() {
155 
156         final double stepSize = 100.;
157 
158         ODEIntegratorBuilder teBuilder = new ThreeEighthesIntegratorBuilder(stepSize);
159         doTestBuildPropagator(teBuilder);
160     }
161 
162     @Test
163     void testClone() {
164 
165         // Given
166         final ODEIntegratorBuilder integratorBuilder = new ClassicalRungeKuttaIntegratorBuilder(3600.0);
167         final Orbit orbit = new CartesianOrbit(new PVCoordinates(
168                 new Vector3D(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS + 400000, 0, 0),
169                 new Vector3D(0, 7668.6, 0)), FramesFactory.getGCRF(),
170                 new AbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
171 
172         final DSSTPropagatorBuilder builder =
173                 new DSSTPropagatorBuilder(orbit, integratorBuilder, 1.0, PropagationType.OSCULATING, PropagationType.OSCULATING);
174 
175         builder.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0)));
176 
177         // Select propagation parameters to test this point
178         builder.getPropagationParametersDrivers().getDrivers().forEach(d -> d.setSelected(true));
179 
180         // When
181         final DSSTPropagatorBuilder copyBuilder = builder.clone();
182 
183         // Then
184         assertDSSTPropagatorBuilderIsACopy(builder, copyBuilder);
185     }
186 
187     private void assertDSSTPropagatorBuilderIsACopy(final DSSTPropagatorBuilder expected,
188                                                     final DSSTPropagatorBuilder actual) {
189         assertPropagatorBuilderIsACopy(expected, actual);
190 
191         Assertions.assertEquals(expected.getIntegratorBuilder().getClass(), actual.getIntegratorBuilder().getClass());
192         Assertions.assertEquals(expected.getPropagationType(), actual.getPropagationType());
193         Assertions.assertEquals(expected.getStateType(), actual.getStateType());
194         Assertions.assertEquals(expected.getMass(), actual.getMass());
195         final List<DSSTForceModel> expectedForceModelList = expected.getAllForceModels();
196         final List<DSSTForceModel> actualForceModelList   = actual.getAllForceModels();
197         Assertions.assertEquals(expectedForceModelList.size(), actualForceModelList.size());
198         for (int i = 0; i < expectedForceModelList.size(); i++) {
199             Assertions.assertEquals(expectedForceModelList.get(i).getClass(), actualForceModelList.get(i).getClass());
200         }
201 
202     }
203 
204     private void doTestBuildPropagator(final ODEIntegratorBuilder foiBuilder) {
205 
206         // We propagate using directly the propagator of the set up
207         final Orbit orbitWithPropagator = propagator.propagate(initDate.shiftedBy(600)).getOrbit();
208 
209         // We propagate using a build version of the propagator
210         // We shall have the same results than before
211         DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
212                                                                   foiBuilder,
213                                                                   1.0,
214                                                                   PropagationType.MEAN,
215                                                                   PropagationType.MEAN);
216 
217         builder.addForceModel(moon);
218         builder.setMass(1000.);
219 
220         final DSSTPropagator prop = builder.buildPropagator();
221 
222         final Orbit orbitWithBuilder = prop.propagate(initDate.shiftedBy(600)).getOrbit();
223 
224         // Verify
225         Assertions.assertEquals(orbitWithPropagator.getA(),             orbitWithBuilder.getA(), 1.e-1);
226         Assertions.assertEquals(orbitWithPropagator.getEquinoctialEx(), orbitWithBuilder.getEquinoctialEx(), eps);
227         Assertions.assertEquals(orbitWithPropagator.getEquinoctialEy(), orbitWithBuilder.getEquinoctialEy(), eps);
228         Assertions.assertEquals(orbitWithPropagator.getHx(),            orbitWithBuilder.getHx(), eps);
229         Assertions.assertEquals(orbitWithPropagator.getHy(),            orbitWithBuilder.getHy(), eps);
230         Assertions.assertEquals(orbitWithPropagator.getLM(),            orbitWithBuilder.getLM(), 8.0e-10);
231 
232     }
233 
234     @Test
235     public void testIssue598() {
236         // Integrator builder
237         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
238         // Propagator builder
239         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
240                                                                   dp54Builder,
241                                                                   1.0,
242                                                                   PropagationType.MEAN,
243                                                                   PropagationType.MEAN);
244         builder.addForceModel(moon);
245         // Verify that there is no Newtonian attraction force model
246         Assertions.assertFalse(hasNewtonianAttraction(builder.getAllForceModels()));
247         // Build the DSST propagator (not used here)
248         builder.buildPropagator();
249         // Verify the addition of the Newtonian attraction force model
250         Assertions.assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
251         // Add a new force model to ensure the Newtonian attraction stay at the last position
252         builder.addForceModel(sun);
253         Assertions.assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
254     }
255 
256     @Test
257     public void testAdditionalEquations() {
258         // Integrator builder
259         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
260         // Propagator builder
261         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
262                                                                   dp54Builder,
263                                                                   1.0,
264                                                                   PropagationType.MEAN,
265                                                                   PropagationType.MEAN);
266         builder.addForceModel(moon);
267         builder.addForceModel(sun);
268 
269         // Add additional equations
270         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
271 
272             public String getName() {
273                 return "linear";
274             }
275 
276             public int getDimension() {
277                 return 1;
278             }
279 
280             public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
281                 return new CombinedDerivatives(new double[] { 1.0 }, null);
282             }
283 
284         });
285 
286         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
287 
288             public String getName() {
289                 return "linear";
290             }
291 
292             public int getDimension() {
293                 return 1;
294             }
295 
296             public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
297                 return new CombinedDerivatives(new double[] { 1.0 }, null);
298             }
299 
300         });
301 
302         try {
303             // Build the propagator
304             builder.buildPropagator();
305             Assertions.fail("an exception should have been thrown");
306         } catch (OrekitException oe) {
307             Assertions.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
308         }
309     }
310 
311     @Test
312     public void testDeselectOrbitals() {
313         // Integrator builder
314         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, toleranceProvider);
315         // Propagator builder
316         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
317                                                                   dp54Builder,
318                                                                   1.0,
319                                                                   PropagationType.MEAN,
320                                                                   PropagationType.MEAN);
321         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
322             Assertions.assertTrue(driver.isSelected());
323         }
324         builder.deselectDynamicParameters();
325         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
326             Assertions.assertFalse(driver.isSelected());
327         }
328     }
329 
330     /** Test for issue #1741.
331      * <p>This test checks that orbital drivers in cloned propagator builders
332      * ain't at the same physical address, i.e. that they're not linked anymore.</p>
333      */
334     @Test
335     void testIssue1741() {
336 
337         // Given
338         final ODEIntegratorBuilder integratorBuilder = new ClassicalRungeKuttaIntegratorBuilder(3600.0);
339         final Orbit orbit = new CartesianOrbit(new PVCoordinates(
340                         new Vector3D(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS + 400000, 0, 0),
341                         new Vector3D(0, 7668.6, 0)), FramesFactory.getGCRF(),
342                                                new AbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
343 
344         final DSSTPropagatorBuilder builder =
345                         new DSSTPropagatorBuilder(orbit, integratorBuilder, 1.0, PropagationType.OSCULATING, PropagationType.OSCULATING);
346 
347         builder.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0)));
348 
349         // When
350         final DSSTPropagatorBuilder copyBuilder = builder.clone();
351 
352         // Change orbit of the copied builder
353         final TimeStampedPVCoordinates modifiedPv = orbit.shiftedBy(3600.).getPVCoordinates();
354         copyBuilder.resetOrbit(new CartesianOrbit(modifiedPv, orbit.getFrame(), orbit.getDate(), orbit.getMu()));
355 
356 
357         // Then
358         // Original builder should still have original orbit
359         final PVCoordinates originalPv = orbit.getPVCoordinates();
360         final PVCoordinates initialPv = builder.createInitialOrbit().getPVCoordinates();
361         final double dP = originalPv.getPosition().distance(initialPv.getPosition());
362         final double dV = originalPv.getVelocity().distance(initialPv.getVelocity());
363         final double dA = originalPv.getAcceleration().distance(initialPv.getAcceleration());
364         Assertions.assertEquals(0., dP, 0.);
365         Assertions.assertEquals(0., dV, 0.);
366         Assertions.assertEquals(0., dA, 0.);
367     }
368 
369     @BeforeEach
370     public void setUp() throws IOException, ParseException {
371 
372         Utils.setDataRoot("regular-data:potential");
373 
374         minStep = 1.0;
375         maxStep = 600.0;
376         toleranceProvider = ToleranceProvider.of(CartesianToleranceProvider.of(10., 0.0001, 
377                 CartesianToleranceProvider.DEFAULT_ABSOLUTE_MASS_TOLERANCE));
378 
379         final Frame earthFrame = FramesFactory.getEME2000();
380         initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000, TimeScalesFactory.getUTC());
381 
382         final double mu = 3.986004415E14;
383         // a    = 42163393.0 m
384         // ex =  -0.25925449177598586
385         // ey =  -0.06946703170551687
386         // hx =   0.15995912655021305
387         // hy =  -0.5969755874197339
388         // lM   = 15.47576793123677 rad
389         orbit = new EquinoctialOrbit(4.2163393E7,
390                                      -0.25925449177598586,
391                                      -0.06946703170551687,
392                                      0.15995912655021305,
393                                      -0.5969755874197339,
394                                      15.47576793123677,
395                                      PositionAngleType.MEAN,
396                                      earthFrame,
397                                      initDate,
398                                      mu);
399 
400         moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), mu);
401         sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), mu);
402 
403         tolerance  = toleranceProvider.getTolerances(orbit, OrbitType.EQUINOCTIAL);
404         propagator = new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]));
405         propagator.setInitialState(new SpacecraftState(orbit).withMass(1000.), PropagationType.MEAN);
406         propagator.addForceModel(moon);
407 
408     }
409 
410     private boolean hasNewtonianAttraction(final List<DSSTForceModel> forceModels) {
411         final int last = forceModels.size() - 1;
412         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
413     }
414 
415 }