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.forces.empirical;
18  
19  import java.util.List;
20  
21  import org.hipparchus.Field;
22  import org.hipparchus.analysis.differentiation.Gradient;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
25  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
26  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
27  import org.hipparchus.util.FastMath;
28  import org.junit.Assert;
29  import org.junit.Before;
30  import org.junit.Test;
31  import org.orekit.Utils;
32  import org.orekit.attitudes.AttitudeProvider;
33  import org.orekit.attitudes.LofOffset;
34  import org.orekit.forces.AbstractForceModelTest;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.frames.LOFType;
38  import org.orekit.orbits.CartesianOrbit;
39  import org.orekit.orbits.FieldKeplerianOrbit;
40  import org.orekit.orbits.KeplerianOrbit;
41  import org.orekit.orbits.Orbit;
42  import org.orekit.orbits.OrbitType;
43  import org.orekit.orbits.PositionAngle;
44  import org.orekit.propagation.FieldSpacecraftState;
45  import org.orekit.propagation.SpacecraftState;
46  import org.orekit.propagation.numerical.FieldNumericalPropagator;
47  import org.orekit.propagation.numerical.NumericalPropagator;
48  import org.orekit.time.AbsoluteDate;
49  import org.orekit.time.DateComponents;
50  import org.orekit.time.FieldAbsoluteDate;
51  import org.orekit.time.TimeComponents;
52  import org.orekit.time.TimeScale;
53  import org.orekit.time.TimeScalesFactory;
54  import org.orekit.utils.Constants;
55  import org.orekit.utils.PVCoordinates;
56  import org.orekit.utils.ParameterDriver;
57  import org.orekit.utils.TimeSpanMap;
58  
59  public class TimeSpanParametricAccelerationTest extends AbstractForceModelTest {
60  
61      @Test
62      public void testGetParameterDrivers() {
63  
64          // A date
65          final TimeScale utc = TimeScalesFactory.getUTC();
66          AbsoluteDate date = new AbsoluteDate("2000-01-01T00:00:00.000", utc);
67          
68          // One AccelerationModel added, only one driver should be in the drivers' array
69          TimeSpanParametricAcceleration forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_K, false,
70                                                                                         new PolynomialAccelerationModel("driver", date, 0));
71  
72          Assert.assertFalse(forceModel.dependsOnPositionOnly());
73          List<ParameterDriver> drivers = forceModel.getParametersDrivers();
74          Assert.assertEquals(1,  drivers.size());
75          Assert.assertEquals("driver[0]",  drivers.get(0).getName());
76  
77          // Extract acceleration model at an arbitrary epoch and check it is the one added
78          PolynomialAccelerationModel accModel = (PolynomialAccelerationModel) forceModel.getAccelerationModel(date);
79          drivers = accModel.getParametersDrivers();
80          Assert.assertEquals(1, drivers.size());
81          Assert.assertEquals(0.0,  drivers.get(0).getValue(), 0.);
82          Assert.assertEquals("driver[0]",  drivers.get(0).getName());
83          
84          // 3 AccelerationModel added, with one default
85          // ----------------------------------------------
86          String prefix = "C";
87          int order = 0;
88          double dt = 120.;
89          // Build the force model
90          accModel = new PolynomialAccelerationModel(prefix, date, order);
91          PolynomialAccelerationModel accModel1 = new PolynomialAccelerationModel(prefix + "1", date, order);
92          PolynomialAccelerationModel accModel2 = new PolynomialAccelerationModel(prefix + "2", date, order);
93          forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_K, false, accModel);
94          forceModel.addAccelerationModelValidAfter(accModel1, date.shiftedBy(dt));
95          forceModel.addAccelerationModelValidBefore(accModel2, date.shiftedBy(-dt));
96          
97          // Extract the drivers and check their values and names
98          drivers = forceModel.getParametersDrivers();
99          Assert.assertEquals(3,  drivers.size());
100         Assert.assertEquals(0.0,  drivers.get(0).getValue(), 0.);
101         Assert.assertEquals("C2[0]", drivers.get(0).getName());
102         Assert.assertEquals(0.0,  drivers.get(1).getValue(), 0.);
103         Assert.assertEquals("C[0]",  drivers.get(1).getName());
104         Assert.assertEquals(0.0,  drivers.get(2).getValue(), 0.);
105         Assert.assertEquals("C1[0]", drivers.get(2).getName());
106         
107         // Check that proper models are returned at significant test dates
108         double eps = 1.e-14;
109         Assert.assertEquals(accModel,  forceModel.getAccelerationModel(date));
110         Assert.assertEquals(accModel,  forceModel.getAccelerationModel(date.shiftedBy(-dt)));
111         Assert.assertEquals(accModel,  forceModel.getAccelerationModel(date.shiftedBy(+dt - eps)));
112         Assert.assertEquals(accModel2, forceModel.getAccelerationModel(date.shiftedBy(-dt - eps)));
113         Assert.assertEquals(accModel2, forceModel.getAccelerationModel(date.shiftedBy(-dt - 86400.)));
114         Assert.assertEquals(accModel1, forceModel.getAccelerationModel(date.shiftedBy(+dt)));
115         Assert.assertEquals(accModel1, forceModel.getAccelerationModel(date.shiftedBy(+dt + 86400.)));
116 
117         // Test #getAccelerationModelSpan method
118         Assert.assertEquals(accModel,  forceModel.getAccelerationModelSpan(date).getData());
119         Assert.assertEquals(accModel2, forceModel.getAccelerationModelSpan(date.shiftedBy(-dt - 86400.)).getData());
120         Assert.assertEquals(accModel1, forceModel.getAccelerationModelSpan(date.shiftedBy(+dt + 1.)).getData());
121         
122         // Test #extractAccelerationModelRange
123         TimeSpanMap<AccelerationModel> dragMap = forceModel.extractAccelerationModelRange(date, date.shiftedBy(dt + 1.));
124         Assert.assertEquals(accModel,  dragMap.getSpan(date).getData());
125         Assert.assertEquals(accModel1, dragMap.getSpan(date.shiftedBy(dt + 86400.)).getData());
126         Assert.assertEquals(accModel,  dragMap.getSpan(date.shiftedBy(-dt - 86400.)).getData());
127     }
128 
129     @Test
130     public void testParameterDerivative() {
131 
132         // Time scale
133         final TimeScale tai = TimeScalesFactory.getTAI();
134 
135         // Low Earth orbit definition (about 360km altitude)
136         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
137         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
138         final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, tai);
139         final SpacecraftState state =
140                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
141                                                        FramesFactory.getGCRF(),
142                                                        date,
143                                                        Constants.EIGEN5C_EARTH_MU));
144 
145         // AccelerationModel
146         final PolynomialAccelerationModel accelerationModel = new PolynomialAccelerationModel("C1", null, 0);
147         accelerationModel.getParametersDrivers().get(0).setValue(1.4);
148         final TimeSpanParametricAcceleration forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_I, true, accelerationModel);
149 
150         // After t2 = t + 4h
151         final double dt2 = 4 * 3600.;
152         final AbsoluteDate date2 = date.shiftedBy(dt2);
153         final PolynomialAccelerationModel accelerationModel2 = new PolynomialAccelerationModel("C2", null, 0);
154         accelerationModel2.getParametersDrivers().get(0).setValue(0.7);
155         forceModel.addAccelerationModelValidAfter(accelerationModel2, date2);
156         
157         // Before t3 = t - 1day
158         final double dt3 = -86400.;
159         final AbsoluteDate date3 = date.shiftedBy(dt3);
160         final PolynomialAccelerationModel accelerationModel3 = new PolynomialAccelerationModel("C3", null, 0);
161         accelerationModel3.getParametersDrivers().get(0).setValue(2.7);
162         forceModel.addAccelerationModelValidBefore(accelerationModel3, date3);
163 
164         // Initialize model
165         forceModel.init(state, null);
166 
167         Assert.assertTrue(forceModel.dependsOnPositionOnly());
168 
169         // Check parameter derivatives at initial date: only "C1" shouldn't be 0.
170         checkParameterDerivative(state, forceModel, "C1[0]" , 1.0e-4, 2.0e-12);
171         checkParameterDerivative(state, forceModel, "C2[0]", 1.0e-4, 0.);
172         checkParameterDerivative(state, forceModel, "C3[0]", 1.0e-4, 0.);
173         
174         // Check parameter derivatives after date2: only "C2" shouldn't be 0.
175         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "C1[0]", 1.0e-4, 0.);
176         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "C2[0]", 1.0e-4, 2.0e-12);
177         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "C3[0]", 1.0e-4, 0.);
178         
179         // Check parameter derivatives after date3: only "C3" shouldn't be 0.
180         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "C1[0]", 1.0e-4, 0.);
181         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "C2[0]", 1.0e-4, 0.);
182         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "C3[0]", 1.0e-4, 2.0e-12);
183     }
184 
185     @Test
186     public void testStateJacobian() {
187 
188         // Initialization
189         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
190                                              new TimeComponents(13, 59, 27.816),
191                                              TimeScalesFactory.getUTC());
192         double i     = FastMath.toRadians(98.7);
193         double omega = FastMath.toRadians(93.0);
194         double OMEGA = FastMath.toRadians(15.0 * 22.5);
195         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
196                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
197                                          Constants.EIGEN5C_EARTH_MU);
198         OrbitType integrationType = OrbitType.CARTESIAN;
199         double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
200         
201         // Time span acceleration force model init
202         double dt = 1. * 3600.;
203         // Build the force model
204         PolynomialAccelerationModel polyAcc0 = new PolynomialAccelerationModel("C0", null, 0);
205         polyAcc0.getParametersDrivers().get(0).setValue(2.7);
206         PolynomialAccelerationModel polyAcc1 = new PolynomialAccelerationModel("C1", null, 0);
207         polyAcc1.getParametersDrivers().get(0).setValue(0.9);
208         PolynomialAccelerationModel polyAcc2 = new PolynomialAccelerationModel("C2", null, 0);
209         polyAcc2.getParametersDrivers().get(0).setValue(1.4);
210         TimeSpanParametricAcceleration forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_J, null, polyAcc0);
211         forceModel.addAccelerationModelValidAfter(polyAcc1, date.shiftedBy(dt));
212         forceModel.addAccelerationModelValidBefore(polyAcc2, date.shiftedBy(-dt));
213         
214         // Check state derivatives inside first AccelerationModel
215         NumericalPropagator propagator =
216                         new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
217                                                                                tolerances[0], tolerances[1]));
218         propagator.setOrbitType(integrationType);
219         propagator.addForceModel(forceModel);
220         SpacecraftState state0 = new SpacecraftState(orbit);
221         // Set target date to 0.5*dt to be inside 1st AccelerationModel
222         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
223         checkStateJacobian(propagator, state0, date.shiftedBy(0.5 * dt),
224                            1e3, tolerances[0], 1.0e-9);
225         
226         // Check state derivatives inside 2nd AccelerationModel
227         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
228                                                                             tolerances[0], tolerances[1]));
229         propagator.setOrbitType(integrationType);
230         propagator.addForceModel(forceModel);
231         // Set target date to 1.5*dt to be inside 2nd AccelerationModel
232         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
233         checkStateJacobian(propagator, state0, date.shiftedBy(1.5 * dt),
234                            1e3, tolerances[0], 1.5e-2);
235         
236         // Check state derivatives inside 3rd AccelerationModel
237         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
238                                                                             tolerances[0], tolerances[1]));
239         propagator.setOrbitType(integrationType);
240         propagator.addForceModel(forceModel);
241         // Set target date to *1.5*dt to be inside 3rd AccelerationModel
242         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
243         checkStateJacobian(propagator, state0, date.shiftedBy(-1.5 * dt),
244                            1e3, tolerances[0], 6.0e-3);
245     }
246 
247     @Test
248     public void testStateJacobianAttitudeOverride() {
249 
250         // Initialization
251         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
252                                              new TimeComponents(13, 59, 27.816),
253                                              TimeScalesFactory.getUTC());
254         double i     = FastMath.toRadians(98.7);
255         double omega = FastMath.toRadians(93.0);
256         double OMEGA = FastMath.toRadians(15.0 * 22.5);
257         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
258                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
259                                          Constants.EIGEN5C_EARTH_MU);
260         OrbitType integrationType = OrbitType.CARTESIAN;
261         double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
262         
263         // Time span acceleration force model init
264         final AttitudeProvider attitudeOverride = new LofOffset(orbit.getFrame(), LOFType.VNC);
265         double dt = 1. * 3600.;
266         // Build the force model
267         PolynomialAccelerationModel polyAcc0 = new PolynomialAccelerationModel("C0", null, 0);
268         polyAcc0.getParametersDrivers().get(0).setValue(2.7);
269         PolynomialAccelerationModel polyAcc1 = new PolynomialAccelerationModel("C1", null, 0);
270         polyAcc1.getParametersDrivers().get(0).setValue(0.9);
271         PolynomialAccelerationModel polyAcc2 = new PolynomialAccelerationModel("C2", null, 0);
272         polyAcc2.getParametersDrivers().get(0).setValue(1.4);
273         TimeSpanParametricAcceleration forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_J, attitudeOverride, polyAcc0);
274         forceModel.addAccelerationModelValidAfter(polyAcc1, date.shiftedBy(dt));
275         forceModel.addAccelerationModelValidBefore(polyAcc2, date.shiftedBy(-dt));
276         
277         // Check state derivatives inside first AccelerationModel
278         NumericalPropagator propagator =
279                         new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
280                                                                                tolerances[0], tolerances[1]));
281         propagator.setOrbitType(integrationType);
282         propagator.addForceModel(forceModel);
283         SpacecraftState state0 = new SpacecraftState(orbit);
284         // Set target date to 0.5*dt to be inside 1st AccelerationModel
285         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
286         checkStateJacobian(propagator, state0, date.shiftedBy(0.5 * dt),
287                            1e3, tolerances[0], 1.7e-9);
288         
289         // Check state derivatives inside 2nd AccelerationModel
290         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
291                                                                             tolerances[0], tolerances[1]));
292         propagator.setOrbitType(integrationType);
293         propagator.addForceModel(forceModel);
294         // Set target date to 1.5*dt to be inside 2nd AccelerationModel
295         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
296         checkStateJacobian(propagator, state0, date.shiftedBy(1.5 * dt),
297                            1e3, tolerances[0], 1.8e-2);
298         
299         // Check state derivatives inside 3rd AccelerationModel
300         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
301                                                                             tolerances[0], tolerances[1]));
302         propagator.setOrbitType(integrationType);
303         propagator.addForceModel(forceModel);
304         // Set target date to *1.5*dt to be inside 3rd AccelerationModel
305         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
306         checkStateJacobian(propagator, state0, date.shiftedBy(-1.5 * dt),
307                            1e3, tolerances[0], 2.0e-2);
308     }
309 
310     @Test
311     public void RealFieldGradientTest() {
312         // Initial field Keplerian orbit
313         // The variables are the six orbital parameters
314         final int freeParameters = 6;
315         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e6);
316         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.01);
317         Gradient i_0 = Gradient.variable(freeParameters, 2, 1.2);
318         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
319         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
320         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
321 
322         Field<Gradient> field = a_0.getField();
323         Gradient zero = field.getZero();
324 
325         // Initial date = J2000 epoch
326         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
327         
328         // J2000 frame
329         Frame EME = FramesFactory.getEME2000();
330 
331         // Create initial field Keplerian orbit
332         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
333                                                                       PositionAngle.MEAN,
334                                                                       EME,
335                                                                       J2000,
336                                                                       zero.add(Constants.EIGEN5C_EARTH_MU));
337         
338         // Initial field and classical S/Cs
339         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
340         SpacecraftState iSR = initialState.toSpacecraftState();
341 
342         // Field integrator and classical integrator
343         ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
344                         new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
345         ClassicalRungeKuttaIntegrator RIntegrator =
346                         new ClassicalRungeKuttaIntegrator(6);
347         OrbitType type = OrbitType.EQUINOCTIAL;
348         
349         // Field and classical numerical propagators
350         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
351         FNP.setOrbitType(type);
352         FNP.setInitialState(initialState);
353 
354         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
355         NP.setOrbitType(type);
356         NP.setInitialState(iSR);
357         
358         // Set up force model
359          
360         // AccelerationModel init
361         double dt = 1000.;
362         // Build the force model
363         PolynomialAccelerationModel polyAcc0 = new PolynomialAccelerationModel("C0", null, 0);
364         polyAcc0.getParametersDrivers().get(0).setValue(2.7);
365         PolynomialAccelerationModel polyAcc1 = new PolynomialAccelerationModel("C1", null, 0);
366         polyAcc1.getParametersDrivers().get(0).setValue(0.9);
367         PolynomialAccelerationModel polyAcc2 = new PolynomialAccelerationModel("C2", null, 0);
368         polyAcc2.getParametersDrivers().get(0).setValue(1.4);
369         TimeSpanParametricAcceleration forceModel = new TimeSpanParametricAcceleration(Vector3D.PLUS_J, null, polyAcc0);
370         forceModel.addAccelerationModelValidAfter(polyAcc1, J2000.toAbsoluteDate().shiftedBy(dt));
371         forceModel.addAccelerationModelValidBefore(polyAcc2, J2000.toAbsoluteDate().shiftedBy(-dt));
372         FNP.addForceModel(forceModel);
373         NP.addForceModel(forceModel);
374         
375         // Do the test
376         // -----------
377         
378         // Propagate inside 1st AccelerationModel
379         checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 0.9 * dt, NP, FNP,
380                                   1.0e-9, 2.4e-02, 7.8e-5, 1.5e-3,
381                                   1, false);
382         
383         // Propagate to 2nd AccelerationModel (reset propagator first)
384         FNP.resetInitialState(initialState);
385         NP.resetInitialState(iSR);
386         checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 1.1 * dt, NP, FNP,
387                                   1.0e-9, 4.3e-02, 8.2e-5, 3.1e-4,
388                                   1, false);
389         
390         // Propagate to 3rd AccelerationModel (reset propagator first)
391         FNP.resetInitialState(initialState);
392         NP.resetInitialState(iSR);
393         checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, -1.1 * dt, NP, FNP,
394                                   1.0e-9, 2.0e-02, 3.9e-04, 4.5e-04,
395                                   1, false);
396     }
397 
398     @Before
399     public void setUp() {
400         Utils.setDataRoot("regular-data");
401     }
402 
403 }