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