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.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.After;
23  import org.junit.Assert;
24  import org.junit.Before;
25  import org.junit.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.events.DateDetector;
36  import org.orekit.propagation.numerical.NumericalPropagator;
37  import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
38  import org.orekit.time.AbsoluteDate;
39  import org.orekit.utils.PVCoordinates;
40  
41  public class AdditionalDerivativesProvidersTest {
42  
43      private double          mu;
44      private AbsoluteDate    initDate;
45      private SpacecraftState initialState;
46      private double[][]      tolerance;
47  
48      /** Test for issue #401
49       *  with a numerical propagator */
50      @Test
51      public void testInitNumerical() {
52  
53          // setup
54          final double reference = 1.25;
55          final double rate      = 1.5;
56          final double dt        = 600.0;
57          Linear linear = new Linear("linear", reference, rate);
58          Assert.assertFalse(linear.wasCalled());
59  
60          // action
61          AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
62          integrator.setInitialStepSize(60);
63          NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
64          propagatorNumerical.setInitialState(initialState.addAdditionalState(linear.getName(), reference));
65          propagatorNumerical.addAdditionalDerivativesProvider(linear);
66          SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
67  
68          // verify
69          Assert.assertTrue(linear.wasCalled());
70          Assert.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0], 1.0e-10);
71  
72      }
73  
74      /** Test for issue #401
75       *  with a DSST propagator */
76      @Test
77      public void testInitDSST() {
78  
79          // setup
80          final double reference = 3.5;
81          final double rate      = 1.5;
82          final double dt        = 600.0;
83          Linear linear = new Linear("linear", reference, rate);
84          Assert.assertFalse(linear.wasCalled());
85  
86          // action
87          AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
88          integrator.setInitialStepSize(60);
89          DSSTPropagator propagatorDSST = new DSSTPropagator(integrator);
90          propagatorDSST.setInitialState(initialState.addAdditionalState(linear.getName(), reference));
91          propagatorDSST.addAdditionalDerivativesProvider(linear);
92          SpacecraftState finalState = propagatorDSST.propagate(initDate.shiftedBy(dt));
93  
94          // verify
95          Assert.assertTrue(linear.wasCalled());
96          Assert.assertEquals(reference + dt * rate, finalState.getAdditionalState(linear.getName())[0], 1.0e-10);
97  
98      }
99  
100     @Test
101     public void testResetState() {
102 
103         // setup
104         final double reference1 = 3.5;
105         final double rate1      = 1.5;
106         Linear linear1 = new Linear("linear-1", reference1, rate1);
107         Assert.assertFalse(linear1.wasCalled());
108         final double reference2 = 4.5;
109         final double rate2      = 1.25;
110         Linear linear2 = new Linear("linear-2", reference2, rate2);
111         Assert.assertFalse(linear2.wasCalled());
112         final double dt = 600;
113 
114         // action
115         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
116         integrator.setInitialStepSize(60);
117         NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
118         propagatorNumerical.setInitialState(initialState.
119                                             addAdditionalState(linear1.getName(), reference1).
120                                             addAdditionalState(linear2.getName(), reference2));
121         propagatorNumerical.addAdditionalDerivativesProvider(linear1);
122         propagatorNumerical.addAdditionalDerivativesProvider(linear2);
123         propagatorNumerical.addEventDetector(new ImpulseManeuver<>(new DateDetector(initDate.shiftedBy(dt / 2.0)),
124                                                                    new Vector3D(0.1, 0.2, 0.3), 350.0));
125         SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
126 
127         // verify
128         Assert.assertTrue(linear1.wasCalled());
129         Assert.assertTrue(linear2.wasCalled());
130         Assert.assertEquals(reference1 + dt * rate1, finalState.getAdditionalState(linear1.getName())[0], 1.0e-10);
131         Assert.assertEquals(reference2 + dt * rate2, finalState.getAdditionalState(linear2.getName())[0], 1.0e-10);
132 
133     }
134 
135     @Test
136     public void testYield() {
137 
138         // setup
139         final double init1 = 1.0;
140         final double init2 = 2.0;
141         final double rate  = 0.5;
142         final double dt    = 600;
143         Yield yield1 = new Yield(null, "yield-1", rate);
144         Yield yield2 = new Yield(yield1.getName(), "yield-2", Double.NaN);
145 
146         // action
147         AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
148         integrator.setInitialStepSize(60);
149         NumericalPropagator propagatorNumerical = new NumericalPropagator(integrator);
150         propagatorNumerical.setInitialState(initialState.
151                                             addAdditionalState(yield1.getName(), init1).
152                                             addAdditionalState(yield2.getName(), init2));
153         propagatorNumerical.addAdditionalDerivativesProvider(yield2); // we intentionally register yield2 before yield 1 to check reordering
154         propagatorNumerical.addAdditionalDerivativesProvider(yield1);
155         SpacecraftState finalState = propagatorNumerical.propagate(initDate.shiftedBy(dt));
156 
157         // verify
158         Assert.assertEquals(init1 + dt * rate, finalState.getAdditionalState(yield1.getName())[0],           1.0e-10);
159         Assert.assertEquals(init2 + dt * rate, finalState.getAdditionalState(yield2.getName())[0],           1.0e-10);
160         Assert.assertEquals(rate,              finalState.getAdditionalStateDerivative(yield1.getName())[0], 1.0e-10);
161         Assert.assertEquals(rate,              finalState.getAdditionalStateDerivative(yield2.getName())[0], 1.0e-10);
162 
163     }
164 
165     @Before
166     public void setUp() {
167         Utils.setDataRoot("regular-data:potential/shm-format");
168         GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
169         mu  = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
170         final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
171         final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
172         initDate = AbsoluteDate.J2000_EPOCH;
173         final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position,  velocity),
174                                                  FramesFactory.getEME2000(), initDate, mu);
175         initialState = new SpacecraftState(orbit);
176         tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL);
177     }
178 
179     @After
180     public void tearDown() {
181         initDate     = null;
182         initialState = null;
183         tolerance    = null;
184     }
185 
186     private static class Linear implements AdditionalDerivativesProvider {
187 
188         private String  name;
189         private double  expectedAtInit;
190         private double  rate;
191         private boolean called;
192 
193         public Linear(final String name, final double expectedAtInit, final double rate) {
194             this.name           = name;
195             this.expectedAtInit = expectedAtInit;
196             this.rate           = rate;
197             this.called         = false;
198         }
199 
200         @Override
201         public void init(SpacecraftState initiaState, AbsoluteDate target) {
202             Assert.assertEquals(expectedAtInit, initiaState.getAdditionalState(getName())[0], 1.0e-15);
203             called = true;
204         }
205 
206         @Override
207         public double[] derivatives(SpacecraftState s) {
208             return new double[] { rate };
209         }
210 
211         @Override
212         public int getDimension() {
213             return 1;
214         }
215 
216         @Override
217         public String getName() {
218             return name;
219         }
220 
221         public boolean wasCalled() {
222             return called;
223         }
224 
225     }
226 
227     private static class Yield implements AdditionalDerivativesProvider {
228 
229         private String dependency;
230         private String name;
231         private double rate;
232 
233         public Yield(final String dependency, final String name, final double rate) {
234             this.dependency = dependency;
235             this.name       = name;
236             this.rate       = rate;
237         }
238 
239         @Override
240         public double[] derivatives(final SpacecraftState s) {
241             return dependency == null ? new double[] { rate } : s.getAdditionalStateDerivative(dependency);
242         }
243 
244         @Override
245         public boolean yield(final SpacecraftState state) {
246             return dependency != null && !state.hasAdditionalStateDerivative(dependency);
247         }
248 
249         @Override
250         public int getDimension() {
251             return 1;
252         }
253 
254         @Override
255         public String getName() {
256             return name;
257         }
258 
259     }
260 
261 }