1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.propagation.conversion;
18  
19  import static org.junit.Assert.assertFalse;
20  import static org.junit.Assert.assertTrue;
21  
22  import java.io.IOException;
23  import java.text.ParseException;
24  import java.util.List;
25  
26  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
27  import org.junit.Assert;
28  import org.junit.Before;
29  import org.junit.Test;
30  import org.orekit.Utils;
31  import org.orekit.bodies.CelestialBodyFactory;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.errors.OrekitMessages;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.orbits.EquinoctialOrbit;
37  import org.orekit.orbits.Orbit;
38  import org.orekit.orbits.OrbitType;
39  import org.orekit.orbits.PositionAngle;
40  import org.orekit.propagation.PropagationType;
41  import org.orekit.propagation.SpacecraftState;
42  import org.orekit.propagation.integration.AdditionalDerivativesProvider;
43  import org.orekit.propagation.numerical.NumericalPropagator;
44  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
45  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
46  import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
47  import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
48  import org.orekit.time.AbsoluteDate;
49  import org.orekit.time.TimeScalesFactory;
50  import org.orekit.utils.ParameterDriver;
51  
52  public class DSSTPropagatorBuilderTest {
53  
54      private static final double eps  = 2.0e-10;
55  
56      private double minStep;
57      private double maxStep;
58      private double dP;
59      private double[][] tolerance;
60  
61      private AbsoluteDate initDate;
62      private EquinoctialOrbit orbit;
63      private DSSTPropagator propagator;
64      private DSSTForceModel moon;
65      private DSSTForceModel sun;
66  
67      @Test
68      public void testIntegrators01() {
69  
70          ODEIntegratorBuilder abBuilder = new AdamsBashforthIntegratorBuilder(2, minStep, maxStep, dP);
71          doTestBuildPropagator(abBuilder);
72      }
73  
74      @Test
75      public void testIntegrators02() {
76  
77          ODEIntegratorBuilder amBuilder = new AdamsMoultonIntegratorBuilder(2, minStep, maxStep, dP);
78          doTestBuildPropagator(amBuilder);
79      }
80  
81      @Test
82      public void testIntegrators03() {
83  
84          final double stepSize = 100.;
85  
86          ODEIntegratorBuilder crkBuilder = new ClassicalRungeKuttaIntegratorBuilder(stepSize);
87          doTestBuildPropagator(crkBuilder);
88      }
89  
90      @Test
91      public void testIntegrators04() {
92  
93          final double stepSize = 100.;
94  
95          ODEIntegratorBuilder lBuilder = new LutherIntegratorBuilder(stepSize);
96          doTestBuildPropagator(lBuilder);
97      }
98  
99      @Test
100     public void testIntegrators05() {
101 
102         ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
103         doTestBuildPropagator(dp54Builder);
104     }
105 
106     @Test
107     public void testIntegrators06() {
108 
109         final double stepSize = 100.;
110 
111         ODEIntegratorBuilder eBuilder = new EulerIntegratorBuilder(stepSize);
112         doTestBuildPropagator(eBuilder);
113     }
114 
115     @Test
116     public void testIntegrators07() {
117 
118         final double stepSize = 100.;
119 
120         ODEIntegratorBuilder gBuilder = new GillIntegratorBuilder(stepSize);
121         doTestBuildPropagator(gBuilder);
122     }
123 
124     @Test
125     public void testIntegrators08() {
126 
127         ODEIntegratorBuilder gbsBuilder = new GraggBulirschStoerIntegratorBuilder(minStep, maxStep, dP);
128         doTestBuildPropagator(gbsBuilder);
129     }
130 
131     @Test
132     public void testIntegrators09() {
133 
134         ODEIntegratorBuilder hh54Builder = new HighamHall54IntegratorBuilder(minStep, maxStep, dP);
135         doTestBuildPropagator(hh54Builder);
136     }
137 
138     @Test
139     public void testIntegrators10() {
140 
141         final double stepSize = 100.;
142 
143         ODEIntegratorBuilder mBuilder = new MidpointIntegratorBuilder(stepSize);
144         doTestBuildPropagator(mBuilder);
145     }
146 
147     @Test
148     public void testIntegrators11() {
149 
150         final double stepSize = 100.;
151 
152         ODEIntegratorBuilder teBuilder = new ThreeEighthesIntegratorBuilder(stepSize);
153         doTestBuildPropagator(teBuilder);
154     }
155 
156     private void doTestBuildPropagator(final ODEIntegratorBuilder foiBuilder) {
157         
158         // We propagate using directly the propagator of the set up
159         final Orbit orbitWithPropagator = propagator.propagate(initDate.shiftedBy(600)).getOrbit();
160         
161         // We propagate using a build version of the propagator
162         // We shall have the same results than before
163         DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
164                                                                   foiBuilder,
165                                                                   1.0,
166                                                                   PropagationType.MEAN,
167                                                                   PropagationType.MEAN);
168         
169         builder.addForceModel(moon);
170         builder.setMass(1000.);
171         
172         final DSSTPropagator prop = builder.buildPropagator(builder.getSelectedNormalizedParameters());
173         
174         final Orbit orbitWithBuilder = prop.propagate(initDate.shiftedBy(600)).getOrbit();
175         
176         // Verify
177         Assert.assertEquals(orbitWithPropagator.getA(),             orbitWithBuilder.getA(), 1.e-1);
178         Assert.assertEquals(orbitWithPropagator.getEquinoctialEx(), orbitWithBuilder.getEquinoctialEx(), eps);
179         Assert.assertEquals(orbitWithPropagator.getEquinoctialEy(), orbitWithBuilder.getEquinoctialEy(), eps);
180         Assert.assertEquals(orbitWithPropagator.getHx(),            orbitWithBuilder.getHx(), eps);
181         Assert.assertEquals(orbitWithPropagator.getHy(),            orbitWithBuilder.getHy(), eps);
182         Assert.assertEquals(orbitWithPropagator.getLM(),            orbitWithBuilder.getLM(), 8.0e-10);
183         
184     }
185 
186     @Test
187     public void testIssue598() {
188         // Integrator builder
189         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
190         // Propagator builder
191         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
192                                                                   dp54Builder,
193                                                                   1.0,
194                                                                   PropagationType.MEAN,
195                                                                   PropagationType.MEAN);
196         builder.addForceModel(moon);
197         // Verify that there is no Newtonian attration force model
198         assertFalse(hasNewtonianAttraction(builder.getAllForceModels()));
199         // Build the DSST propagator (not used here)
200         builder.buildPropagator(builder.getSelectedNormalizedParameters());
201         // Verify the addition of the Newtonian attraction force model
202         assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
203         // Add a new force model to ensure the Newtonian attraction stay at the last position
204         builder.addForceModel(sun);
205         assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
206     }
207 
208     @Test
209     public void testAdditionalEquations() {
210         // Integrator builder
211         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
212         // Propagator builder
213         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
214                                                                   dp54Builder,
215                                                                   1.0,
216                                                                   PropagationType.MEAN,
217                                                                   PropagationType.MEAN);
218         builder.addForceModel(moon);
219         builder.addForceModel(sun);
220 
221         // Add additional equations
222         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
223 
224             public String getName() {
225                 return "linear";
226             }
227 
228             public int getDimension() {
229                 return 1;
230             }
231 
232             public double[] derivatives(SpacecraftState s) {
233                 return new double[] { 1.0 };
234             }
235         });
236 
237         builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
238 
239     	    public String getName() {
240     	        return "linear";
241     	    }
242 
243     	    public int getDimension() {
244     	        return 1;
245     	    }
246 
247             public double[] derivatives(SpacecraftState s) {
248                 return new double[] { 1.0 };
249             }
250         });
251 
252         try {
253     	    // Build the propagator
254     	    builder.buildPropagator(builder.getSelectedNormalizedParameters());
255             Assert.fail("an exception should have been thrown");
256         } catch (OrekitException oe) {
257             Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
258         }
259     }
260     
261     @Test
262     public void testDeselectOrbitals() {
263         // Integrator builder
264         final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
265         // Propagator builder
266         final DSSTPropagatorBuilder builder = new DSSTPropagatorBuilder(orbit,
267                                                                   dp54Builder,
268                                                                   1.0,
269                                                                   PropagationType.MEAN,
270                                                                   PropagationType.MEAN);
271         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
272             Assert.assertTrue(driver.isSelected());
273         }
274         builder.deselectDynamicParameters();
275         for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
276             Assert.assertFalse(driver.isSelected());
277         }
278     }
279     
280     @Before
281     public void setUp() throws IOException, ParseException {
282         
283         Utils.setDataRoot("regular-data");
284 
285         minStep = 1.0;
286         maxStep = 600.0;
287         dP      = 10.0;
288 
289         final Frame earthFrame = FramesFactory.getEME2000();
290         initDate = new AbsoluteDate(2003, 07, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
291         
292         final double mu = 3.986004415E14;
293         // a    = 42163393.0 m
294         // ex =  -0.25925449177598586
295         // ey =  -0.06946703170551687
296         // hx =   0.15995912655021305
297         // hy =  -0.5969755874197339
298         // lM   = 15.47576793123677 rad
299         orbit = new EquinoctialOrbit(4.2163393E7,
300                                      -0.25925449177598586,
301                                      -0.06946703170551687,
302                                      0.15995912655021305,
303                                      -0.5969755874197339,
304                                      15.47576793123677,
305                                      PositionAngle.MEAN,
306                                      earthFrame,
307                                      initDate,
308                                      mu);
309 
310         moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), mu);
311         sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), mu);
312 
313         tolerance  = NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
314         propagator = new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]));
315         propagator.setInitialState(new SpacecraftState(orbit, 1000.), PropagationType.MEAN);
316         propagator.addForceModel(moon);
317 
318     }
319 
320     private boolean hasNewtonianAttraction(final List<DSSTForceModel> forceModels) {
321         final int last = forceModels.size() - 1;
322         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
323     }
324 
325 }