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.integration;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
21  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
22  import org.junit.jupiter.api.AfterEach;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.BeforeEach;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.forces.gravity.potential.GravityFieldFactory;
28  import org.orekit.forces.gravity.potential.SHMFormatReader;
29  import org.orekit.forces.maneuvers.ImpulseManeuver;
30  import org.orekit.frames.FramesFactory;
31  import org.orekit.orbits.EquinoctialOrbit;
32  import org.orekit.orbits.Orbit;
33  import org.orekit.orbits.OrbitType;
34  import org.orekit.propagation.SpacecraftState;
35  import org.orekit.propagation.ToleranceProvider;
36  import org.orekit.propagation.events.DateDetector;
37  import org.orekit.propagation.numerical.NumericalPropagator;
38  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
39  import org.orekit.time.AbsoluteDate;
40  import org.orekit.utils.PVCoordinates;
41  
42  public class AdditionalDerivativesProvidersTest {
43  
44      private double          mu;
45      private AbsoluteDate    initDate;
46      private SpacecraftState initialState;
47      private double[][]      tolerance;
48  
49      /** Test for issue #401
50       *  with a numerical propagator */
51      @Test
52      public void testInitNumerical() {
53  
54          // setup
55          final double reference = 1.25;
56          final double rate      = 1.5;
57          final double dt        = 600.0;
58          Linear linear = new Linear("linear", reference, rate);
59          Assertions.assertFalse(linear.wasCalled());
60  
61          // action
62          AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
63          integrator.setInitialStepSize(60);
64          NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
65          propagatorNumerical.setInitialState(initialState.addAdditionalData(linear.getName(), reference));
66          propagatorNumerical.addAdditionalDerivativesProvider(linear);
67          SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
68  
69          // verify
70          Assertions.assertTrue(linear.wasCalled());
71          Assertions.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0], 1.0e-10);
72  
73      }
74  
75      /** Test for issue #401
76       *  with a DSST propagator */
77      @Test
78      public void testInitDSST() {
79  
80          // setup
81          final double reference = 3.5;
82          final double rate      = 1.5;
83          final double dt        = 600.0;
84          Linear linear = new Linear("linear", reference, rate);
85          Assertions.assertFalse(linear.wasCalled());
86  
87          // action
88          AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
89          integrator.setInitialStepSize(60);
90          DSSTPropagator propagatorDSST = new DSSTPropagator(integrator);
91          propagatorDSST.setInitialState(initialState.addAdditionalData(linear.getName(), reference));
92          propagatorDSST.addAdditionalDerivativesProvider(linear);
93          SpacecraftState finalState = propagatorDSST.propagate(initDate.shiftedBy(dt));
94  
95          // verify
96          Assertions.assertTrue(linear.wasCalled());
97          Assertions.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0], 1.0e-10);
98  
99      }
100 
101     @Test
102     public void testResetState() {
103 
104         // setup
105         final double reference1 = 3.5;
106         final double rate1      = 1.5;
107         Linear linear1 = new Linear("linear-1", reference1, rate1);
108         Assertions.assertFalse(linear1.wasCalled());
109         final double reference2 = 4.5;
110         final double rate2      = 1.25;
111         Linear linear2 = new Linear("linear-2", reference2, rate2);
112         Assertions.assertFalse(linear2.wasCalled());
113         final double dt = 600;
114 
115         // action
116         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
117         integrator.setInitialStepSize(60);
118         NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
119         propagatorNumerical.setInitialState(initialState.
120                                             addAdditionalData(linear1.getName(), reference1).
121                                             addAdditionalData(linear2.getName(), reference2));
122         propagatorNumerical.addAdditionalDerivativesProvider(linear1);
123         propagatorNumerical.addAdditionalDerivativesProvider(linear2);
124         propagatorNumerical.addEventDetector(new ImpulseManeuver(new DateDetector(initDate.shiftedBy(dt / 2.0)),
125                                                                  new Vector3D(0.1, 0.2, 0.3), 350.0));
126         SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
127 
128         // verify
129         Assertions.assertTrue(linear1.wasCalled());
130         Assertions.assertTrue(linear2.wasCalled());
131         Assertions.assertEquals(reference1 + dt * rate1, finalState.getAdditionalState(linear1.getName())[0], 1.0e-10);
132         Assertions.assertEquals(reference2 + dt * rate2, finalState.getAdditionalState(linear2.getName())[0], 1.0e-10);
133 
134     }
135 
136     @Test
137     public void testYield() {
138 
139         // setup
140         final double init1 = 1.0;
141         final double init2 = 2.0;
142         final double rate  = 0.5;
143         final double dt    = 600;
144         Yield yield1 = new Yield(null, "yield-1", rate);
145         Yield yield2 = new Yield(yield1.getName(), "yield-2", Double.NaN);
146 
147         // action
148         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
149         integrator.setInitialStepSize(60);
150         NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
151         propagatorNumerical.setInitialState(initialState.
152                                             addAdditionalData(yield1.getName(), init1).
153                                             addAdditionalData(yield2.getName(), init2));
154         propagatorNumerical.addAdditionalDerivativesProvider(yield2); // we intentionally register yield2 before yield 1 to check reordering
155         propagatorNumerical.addAdditionalDerivativesProvider(yield1);
156         SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
157 
158         // verify
159         Assertions.assertEquals(init1 + dt * rate, finalState.getAdditionalState(yield1.getName())[0], 1.0e-10);
160         Assertions.assertEquals(init2 + dt * rate, finalState.getAdditionalState(yield2.getName())[0], 1.0e-10);
161         Assertions.assertEquals(rate, finalState.getAdditionalStateDerivative(yield1.getName())[0], 1.0e-10);
162         Assertions.assertEquals(rate, finalState.getAdditionalStateDerivative(yield2.getName())[0], 1.0e-10);
163 
164     }
165 
166     @Test
167     public void testCoupling() {
168 
169         // setup
170         final double   dt       = 600.0;
171         final Coupling coupling = new Coupling("coupling", 3.5, -2.0, 1.0);
172 
173         // action
174         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
175         integrator.setInitialStepSize(60);
176         NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
177         propagatorNumerical.setInitialState(initialState.
178                                             addAdditionalData(coupling.getName(),
179                                                                coupling.secondaryInit));
180         propagatorNumerical.addAdditionalDerivativesProvider(coupling);
181         SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
182 
183         // verify
184         Assertions.assertEquals(coupling.secondaryInit + dt * coupling.secondaryRate,
185                 finalState.getAdditionalState(coupling.getName())[0], 1.0e-10);
186         Assertions.assertEquals(initialState.getOrbit().getA() + dt * coupling.smaRate, finalState.getOrbit().getA(), 1.0e-10);
187 
188     }
189 
190     @BeforeEach
191     public void setUp() {
192         Utils.setDataRoot("regular-data:potential/shm-format");
193         GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
194         mu  = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
195         final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
196         final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
197         initDate = AbsoluteDate.J2000_EPOCH;
198         final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position,  velocity),
199                                                  FramesFactory.getEME2000(), initDate, mu);
200         initialState = new SpacecraftState(orbit);
201         tolerance = ToleranceProvider.getDefaultToleranceProvider(0.001).getTolerances(orbit, OrbitType.EQUINOCTIAL);
202     }
203 
204     @AfterEach
205     public void tearDown() {
206         initDate     = null;
207         initialState = null;
208         tolerance    = null;
209     }
210 
211     private static class Linear implements AdditionalDerivativesProvider {
212 
213         private String  name;
214         private double  expectedAtInit;
215         private double  rate;
216         private boolean called;
217 
218         public Linear(final String name, final double expectedAtInit, final double rate) {
219             this.name           = name;
220             this.expectedAtInit = expectedAtInit;
221             this.rate           = rate;
222             this.called         = false;
223         }
224 
225         @Override
226         public void init(SpacecraftState initiaState, AbsoluteDate target) {
227             Assertions.assertEquals(expectedAtInit, initiaState.getAdditionalState(getName())[0], 1.0e-15);
228             called = true;
229         }
230 
231         @Override
232         public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
233             return new CombinedDerivatives(new double[] { rate }, null);
234         }
235 
236         @Override
237         public int getDimension() {
238             return 1;
239         }
240 
241         @Override
242         public String getName() {
243             return name;
244         }
245 
246         public boolean wasCalled() {
247             return called;
248         }
249 
250     }
251 
252     private static class Yield implements AdditionalDerivativesProvider {
253 
254         private String dependency;
255         private String name;
256         private double rate;
257 
258         public Yield(final String dependency, final String name, final double rate) {
259             this.dependency = dependency;
260             this.name       = name;
261             this.rate       = rate;
262         }
263 
264         @Override
265         public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
266             return new CombinedDerivatives(dependency == null ?
267                                            new double[] { rate } :
268                                            s.getAdditionalStateDerivative(dependency),
269                                            null);
270         }
271 
272         @Override
273         public boolean yields(final SpacecraftState state) {
274             return dependency != null && !state.hasAdditionalStateDerivative(dependency);
275         }
276 
277         @Override
278         public int getDimension() {
279             return 1;
280         }
281 
282         @Override
283         public String getName() {
284             return name;
285         }
286 
287     }
288 
289     private static class Coupling implements AdditionalDerivativesProvider {
290 
291         private final String  name;
292         private final double  secondaryInit;
293         private final double  secondaryRate;
294         private final double  smaRate;
295 
296         public Coupling(final String name,
297                         final double secondaryInit, final double secondaryRate,
298                         final double smaRate) {
299             this.name          = name;
300             this.secondaryInit = secondaryInit;
301             this.secondaryRate = secondaryRate;
302             this.smaRate       = smaRate;
303         }
304 
305         @Override
306         public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
307             return new CombinedDerivatives(new double[] { secondaryRate },
308                                            new double[] { smaRate, 0.0, 0.0, 0.0, 0.0, 0.0 });
309         }
310 
311         @Override
312         public int getDimension() {
313             return 1;
314         }
315 
316         @Override
317         public String getName() {
318             return name;
319         }
320 
321     }
322 
323 }