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.drag;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.analysis.differentiation.DSFactory;
21  import org.hipparchus.analysis.differentiation.DerivativeStructure;
22  import org.hipparchus.analysis.differentiation.Gradient;
23  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
26  import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
27  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.BeforeEach;
31  import org.junit.jupiter.api.Test;
32  import org.orekit.Utils;
33  import org.orekit.attitudes.AttitudeProvider;
34  import org.orekit.attitudes.LofOffset;
35  import org.orekit.bodies.CelestialBody;
36  import org.orekit.bodies.CelestialBodyFactory;
37  import org.orekit.bodies.OneAxisEllipsoid;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  import org.orekit.forces.AbstractLegacyForceModelTest;
41  import org.orekit.forces.BoxAndSolarArraySpacecraft;
42  import org.orekit.forces.ForceModel;
43  import org.orekit.frames.Frame;
44  import org.orekit.frames.FramesFactory;
45  import org.orekit.frames.LOFType;
46  import org.orekit.frames.Transform;
47  import org.orekit.models.earth.atmosphere.Atmosphere;
48  import org.orekit.models.earth.atmosphere.HarrisPriester;
49  import org.orekit.orbits.CartesianOrbit;
50  import org.orekit.orbits.FieldKeplerianOrbit;
51  import org.orekit.orbits.KeplerianOrbit;
52  import org.orekit.orbits.Orbit;
53  import org.orekit.orbits.OrbitType;
54  import org.orekit.orbits.PositionAngleType;
55  import org.orekit.propagation.FieldSpacecraftState;
56  import org.orekit.propagation.SpacecraftState;
57  import org.orekit.propagation.ToleranceProvider;
58  import org.orekit.propagation.numerical.FieldNumericalPropagator;
59  import org.orekit.propagation.numerical.NumericalPropagator;
60  import org.orekit.time.AbsoluteDate;
61  import org.orekit.time.DateComponents;
62  import org.orekit.time.FieldAbsoluteDate;
63  import org.orekit.time.TimeComponents;
64  import org.orekit.time.TimeScale;
65  import org.orekit.time.TimeScalesFactory;
66  import org.orekit.utils.Constants;
67  import org.orekit.utils.FieldPVCoordinates;
68  import org.orekit.utils.IERSConventions;
69  import org.orekit.utils.PVCoordinates;
70  import org.orekit.utils.ParameterDriver;
71  import org.orekit.utils.TimeSpanMap;
72  
73  import java.util.List;
74  
75  class TimeSpanDragForceTest extends AbstractLegacyForceModelTest {
76  
77      /** UTC time scale. */
78      private TimeScale utc;
79  
80      /** Compute acceleration derivatives around input position at input date.
81       *  Using finite differences in position.
82       */
83      @Override
84      protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
85                                                                           final FieldSpacecraftState<DerivativeStructure> state) {
86          try {
87  
88              final AbsoluteDate                       date     = state.getDate().toAbsoluteDate();
89              final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
90              final FieldVector3D<DerivativeStructure> velocity = state.getPVCoordinates().getVelocity();
91              java.lang.reflect.Field atmosphereField = AbstractDragForceModel.class.getDeclaredField("atmosphere");
92              atmosphereField.setAccessible(true);
93              Atmosphere atmosphere = (Atmosphere) atmosphereField.get(forceModel);
94              java.lang.reflect.Field spacecraftField = DragForce.class.getDeclaredField("spacecraft");
95              spacecraftField.setAccessible(true);
96              // Get the DragSensitive model at date
97              DragSensitive spacecraft = ((TimeSpanDragForce) (forceModel)).getDragSensitive(date);
98  
99              // retrieve derivation properties
100             final DSFactory factory = state.getMass().getFactory();
101 
102             // get atmosphere properties in atmosphere own frame
103             final Frame      atmFrame  = atmosphere.getFrame();
104             final Transform  toBody    = state.getFrame().getTransformTo(atmFrame, date);
105             final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position);
106             final Vector3D   posBody   = posBodyDS.toVector3D();
107             final Vector3D   vAtmBody  = atmosphere.getVelocity(date, posBody, atmFrame);
108 
109             // estimate density model
110             if (factory.getCompiler().getOrder() > 1) {
111                 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_DERIVATION_ORDER, factory.getCompiler().getOrder());
112             }
113 
114             final DerivativeStructure rho = atmosphere.getDensity(new FieldAbsoluteDate<>(factory.getDerivativeField(), date),
115                     posBodyDS, atmFrame);
116 
117             // we consider that at first order the atmosphere velocity in atmosphere frame
118             // does not depend on local position; however atmosphere velocity in inertial
119             // frame DOES depend on position since the transform between the frames depends
120             // on it, due to central body rotation rate and velocity composition.
121             // So we use the transform to get the correct partial derivatives on vAtm
122             final FieldVector3D<DerivativeStructure> vAtmBodyDS =
123                             new FieldVector3D<>(factory.constant(vAtmBody.getX()),
124                                             factory.constant(vAtmBody.getY()),
125                                             factory.constant(vAtmBody.getZ()));
126             final FieldPVCoordinates<DerivativeStructure> pvAtmBody = new FieldPVCoordinates<>(posBodyDS, vAtmBodyDS);
127             final FieldPVCoordinates<DerivativeStructure> pvAtm     = toBody.getInverse().transformPVCoordinates(pvAtmBody);
128 
129             // now we can compute relative velocity, it takes into account partial derivatives with respect to position
130             final FieldVector3D<DerivativeStructure> relativeVelocity = pvAtm.getVelocity().subtract(velocity);
131 
132 
133             // Extract drag parameters of the proper model
134             DerivativeStructure[] allParameters = forceModel.getParameters(factory.getDerivativeField());
135             DerivativeStructure[] parameters = ((TimeSpanDragForce) (forceModel)).extractParameters(allParameters, state.getDate());
136 
137             // compute acceleration with all its partial derivatives
138             return spacecraft.dragAcceleration(state, rho, relativeVelocity, parameters);
139         } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
140             return null;
141         }
142     }
143 
144     /** Compute acceleration derivatives around input position at input date.
145      *  Using finite differences in position.
146      */
147     @Override
148     protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
149                                                                       final FieldSpacecraftState<Gradient> state) {
150         try {
151 
152             final AbsoluteDate                       date     = state.getDate().toAbsoluteDate();
153             final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
154             final FieldVector3D<Gradient> velocity = state.getPVCoordinates().getVelocity();
155             java.lang.reflect.Field atmosphereField = AbstractDragForceModel.class.getDeclaredField("atmosphere");
156             atmosphereField.setAccessible(true);
157             Atmosphere atmosphere = (Atmosphere) atmosphereField.get(forceModel);
158             java.lang.reflect.Field spacecraftField = DragForce.class.getDeclaredField("spacecraft");
159             spacecraftField.setAccessible(true);
160             // Get the DragSensitive model at date
161             DragSensitive spacecraft = ((TimeSpanDragForce) (forceModel)).getDragSensitive(date);
162 
163             final int freeParameters = state.getMass().getFreeParameters();
164 
165             // get atmosphere properties in atmosphere own frame
166             final Frame      atmFrame  = atmosphere.getFrame();
167             final Transform  toBody    = state.getFrame().getTransformTo(atmFrame, date);
168             final FieldVector3D<Gradient> posBodyG = toBody.transformPosition(position);
169             final Vector3D   posBody   = posBodyG.toVector3D();
170             final Vector3D   vAtmBody  = atmosphere.getVelocity(date, posBody, atmFrame);
171 
172             // estimate density model
173             final Gradient rho = atmosphere.getDensity(new FieldAbsoluteDate<>(posBodyG.getX().getField(), date),
174                     posBodyG, atmFrame);
175 
176             // we consider that at first order the atmosphere velocity in atmosphere frame
177             // does not depend on local position; however atmosphere velocity in inertial
178             // frame DOES depend on position since the transform between the frames depends
179             // on it, due to central body rotation rate and velocity composition.
180             // So we use the transform to get the correct partial derivatives on vAtm
181             final FieldVector3D<Gradient> vAtmBodyG =
182                             new FieldVector3D<>(Gradient.constant(freeParameters, vAtmBody.getX()),
183                                             Gradient.constant(freeParameters, vAtmBody.getY()),
184                                             Gradient.constant(freeParameters, vAtmBody.getZ()));
185             final FieldPVCoordinates<Gradient> pvAtmBody = new FieldPVCoordinates<>(posBodyG, vAtmBodyG);
186             final FieldPVCoordinates<Gradient> pvAtm     = toBody.getInverse().transformPVCoordinates(pvAtmBody);
187 
188             // now we can compute relative velocity, it takes into account partial derivatives with respect to position
189             final FieldVector3D<Gradient> relativeVelocity = pvAtm.getVelocity().subtract(velocity);
190 
191 
192             // Extract drag parameters of the proper model
193             Gradient[] allParameters = forceModel.getParameters(state.getDate().getField());
194             Gradient[] parameters = ((TimeSpanDragForce) (forceModel)).extractParameters(allParameters, state.getDate());
195 
196             // compute acceleration with all its partial derivatives
197             return spacecraft.dragAcceleration(state, rho, relativeVelocity, parameters);
198         } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
199             return null;
200         }
201     }
202 
203     /** Test that the getParameterDrivers method is working as expected
204      * on an IsotropicDrag-based (ie. spherical) TimeSpanDragForce model.
205      */
206     @Test
207     void testGetParameterDriversSphere() {
208 
209         // Atmosphere
210         final Atmosphere atmosphere = getAtmosphere();
211 
212         // A date
213         AbsoluteDate date = new AbsoluteDate("2000-01-01T00:00:00.000", TimeScalesFactory.getUTC());
214 
215         // One IsotropicDrag added, only one driver should be in the drivers' array
216         // its name should be the default name: IsotropicDrag.DRAG_COEFFICIENT
217         // -----------------------
218         double dragArea = 2.;
219         double dragCd0 = 0.;
220         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, new IsotropicDrag(dragArea, dragCd0));
221         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
222         List<ParameterDriver> drivers = forceModel.getParametersDrivers();
223         Assertions.assertEquals(2,  drivers.size());
224         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
225         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
226         Assertions.assertEquals(dragCd0,  drivers.get(1).getValue(), 0.);
227         Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT,  drivers.get(1).getName());
228 
229         // Extract drag model at an arbitrary epoch and check it is the one added
230         IsotropicDrag isoDrag = (IsotropicDrag) forceModel.getDragSensitive(date);
231         drivers = isoDrag.getDragParametersDrivers();
232         Assertions.assertEquals(2, drivers.size());
233         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
234         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
235         Assertions.assertEquals(dragCd0,  drivers.get(1).getValue(), 0.);
236         Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT,  drivers.get(1).getName());
237 
238         // 3 IsotropicDrag models added, with one default
239         // ----------------------------------------------
240         double dragCd1 = 1.;
241         double dragCd2 = 2.;
242         double dt = 120.;
243         // Build the force model
244         isoDrag = new IsotropicDrag(dragArea, dragCd0);
245         IsotropicDrag isoDrag1 = new IsotropicDrag(dragArea, dragCd1);
246         IsotropicDrag isoDrag2 = new IsotropicDrag(dragArea, dragCd2);
247         forceModel = new TimeSpanDragForce(atmosphere, isoDrag);
248         forceModel.addDragSensitiveValidAfter(isoDrag1, date.shiftedBy(dt));
249         forceModel.addDragSensitiveValidBefore(isoDrag2, date.shiftedBy(-dt));
250 
251         // Extract the drivers and check their values and names
252         drivers = forceModel.getParametersDrivers();
253         Assertions.assertEquals(6,  drivers.size());
254         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
255         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_BEFORE + date.shiftedBy(-dt).toString(utc),
256                             drivers.get(0).getName());
257         Assertions.assertEquals(dragCd2,  drivers.get(1).getValue(), 0.);
258         Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT + TimeSpanDragForce.DATE_BEFORE + date.shiftedBy(-dt).toString(utc),
259                             drivers.get(1).getName());
260         Assertions.assertEquals(1.0,  drivers.get(2).getValue(), 0.);
261         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(2).getName());
262         Assertions.assertEquals(dragCd0,  drivers.get(3).getValue(), 0.);
263         Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT,  drivers.get(3).getName());
264         Assertions.assertEquals(1.0,  drivers.get(4).getValue(), 0.);
265         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_AFTER + date.shiftedBy(+dt).toString(utc),
266                             drivers.get(4).getName());
267         Assertions.assertEquals(dragCd1,  drivers.get(5).getValue(), 0.);
268         Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT + TimeSpanDragForce.DATE_AFTER + date.shiftedBy(+dt).toString(utc),
269                             drivers.get(5).getName());
270 
271         // Check that proper models are returned at significant test dates
272         // Cd0 model
273         double eps = 1.e-14;
274         Assertions.assertEquals(isoDrag, forceModel.getDragSensitive(date));
275         Assertions.assertEquals(isoDrag, forceModel.getDragSensitive(date.shiftedBy(-dt)));
276         Assertions.assertEquals(isoDrag, forceModel.getDragSensitive(date.shiftedBy(+dt - eps)));
277         // Cd2 model
278         Assertions.assertEquals(isoDrag2, forceModel.getDragSensitive(date.shiftedBy(-dt - eps)));
279         Assertions.assertEquals(isoDrag2, forceModel.getDragSensitive(date.shiftedBy(-dt - 86400.)));
280         // Cd1 model
281         Assertions.assertEquals(isoDrag1, forceModel.getDragSensitive(date.shiftedBy(+dt)));
282         Assertions.assertEquals(isoDrag1, forceModel.getDragSensitive(date.shiftedBy(+dt + 86400.)));
283 
284         // Add a custom-named driver
285         // ------------
286         double dragCd3 = 3.;
287         IsotropicDrag isoDrag3 = new IsotropicDrag(dragArea, dragCd3);
288         isoDrag3.getDragParametersDrivers().get(0).setName("custom-fact");
289         isoDrag3.getDragParametersDrivers().get(1).setName("custom-Cd");
290         forceModel.addDragSensitiveValidAfter(isoDrag3, date.shiftedBy(2. * dt));
291         drivers = forceModel.getParametersDrivers();
292         Assertions.assertEquals(8,  drivers.size());
293         Assertions.assertEquals(1.0,  drivers.get(6).getValue(), 0.);
294         Assertions.assertEquals("custom-fact", drivers.get(6).getName());
295         Assertions.assertEquals(dragCd3,  drivers.get(7).getValue(), 0.);
296         Assertions.assertEquals("custom-Cd", drivers.get(7).getName());
297 
298 
299         // Test #getDragSensitiveSpan method
300         Assertions.assertEquals(isoDrag, forceModel.getDragSensitiveSpan(date).getData());
301         Assertions.assertEquals(isoDrag2, forceModel.getDragSensitiveSpan(date.shiftedBy(-dt - 86400.)).getData());
302         Assertions.assertEquals(isoDrag1, forceModel.getDragSensitiveSpan(date.shiftedBy(+dt + 1.)).getData());
303         Assertions.assertEquals(isoDrag3, forceModel.getDragSensitiveSpan(date.shiftedBy(2 * dt + 1.)).getData());
304 
305         // Test #extractDragSensitiveRange
306         TimeSpanMap<DragSensitive> dragMap = forceModel.extractDragSensitiveRange(date, date.shiftedBy(dt + 1.));
307         Assertions.assertEquals(isoDrag, dragMap.getSpan(date).getData());
308         Assertions.assertEquals(isoDrag1, dragMap.getSpan(date.shiftedBy(dt + 86400.)).getData());
309         Assertions.assertEquals(isoDrag, dragMap.getSpan(date.shiftedBy(-dt - 86400.)).getData());
310     }
311 
312 
313     /** Test parameter derivatives for an IsotropicDrag TimeSpanDragForce.
314      *  This test is more or less a copy of the same one in DragForceTest class
315      *  with addition of several IsotropicDrag models valid at different dates
316      *  to test that the different parameters' derivatives are computed correctly.
317      */
318     @Test
319     void testParameterDerivativeSphere() {
320 
321         // Low Earth orbit definition (about 360km altitude)
322         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
323         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
324         final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
325         final SpacecraftState state =
326                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
327                                                        FramesFactory.getGCRF(),
328                                                        date,
329                                                        Constants.EIGEN5C_EARTH_MU));
330 
331         // Atmosphere
332         final Atmosphere atmosphere = getAtmosphere();
333 
334         // Constant area for the different tests
335         final double dragArea = 2.5;
336 
337         // Initialize force model (first coef is valid at all epochs)
338         final double dragCd  = 1.2;
339         final IsotropicDrag isotropicDrag = new IsotropicDrag(dragArea, dragCd);
340         isotropicDrag.getDragParametersDrivers().get(0).setName("Cd");
341         final TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, isotropicDrag, TimeScalesFactory.getUTC());
342 
343 
344         // After t2 = t + 4h
345         final double dragCd2 = 3.;
346         final double dt2 = 4 * 3600.;
347         final AbsoluteDate date2 = date.shiftedBy(dt2);
348         final IsotropicDrag isotropicDrag2 = new IsotropicDrag(dragArea, dragCd2);
349         isotropicDrag2.getDragParametersDrivers().get(0).setName("Cd2");
350         forceModel.addDragSensitiveValidAfter(isotropicDrag2, date2);
351 
352         // Before t3 = t - 1day
353         final double dragCd3 = 3.;
354         final double dt3 = -86400.;
355         final AbsoluteDate date3 = date.shiftedBy(dt3);
356         final IsotropicDrag isotropicDrag3 = new IsotropicDrag(dragArea, dragCd3);
357         isotropicDrag3.getDragParametersDrivers().get(0).setName("Cd3");
358         forceModel.addDragSensitiveValidBefore(isotropicDrag3, date3);
359 
360 
361         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
362 
363         // Check parameter derivatives at initial date: only "Cd" shouldn't be 0.
364         checkParameterDerivative(state, forceModel, "Cd" , 1.0e-4, 2.0e-12);
365         checkParameterDerivative(state, forceModel, "Cd2", 1.0e-4, 0.);
366         checkParameterDerivative(state, forceModel, "Cd3", 1.0e-4, 0.);
367 
368         // Check parameter derivatives after date2: only "Cd2" shouldn't be 0.
369         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "Cd", 1.0e-4, 0.);
370         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "Cd2", 1.0e-4, 2.0e-12);
371         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "Cd3", 1.0e-4, 0.);
372 
373         // Check parameter derivatives after date3: only "Cd3" shouldn't be 0.
374         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "Cd", 1.0e-4, 0.);
375         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "Cd2", 1.0e-4, 0.);
376         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "Cd3", 1.0e-4, 2.0e-12);
377     }
378 
379     /** Test parameter derivatives for an IsotropicDrag TimeSpanDragForce.
380      *  This test is more or less a copy of the same one in DragForceTest class
381      *  with addition of several IsotropicDrag models valid at different dates
382      *  to test that the different parameters' derivatives are computed correctly.
383      */
384     @Test
385     void testParameterDerivativeSphereGradient() {
386 
387         // Low Earth orbit definition (about 360km altitude)
388         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
389         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
390         final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
391         final SpacecraftState state =
392                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
393                                                        FramesFactory.getGCRF(),
394                                                        date,
395                                                        Constants.EIGEN5C_EARTH_MU));
396 
397         // Atmosphere
398         final Atmosphere atmosphere = getAtmosphere();
399 
400         // Constant area for the different tests
401         final double dragArea = 2.5;
402 
403         // Initialize force model (first coef is valid at all epochs)
404         final double dragCd  = 1.2;
405         final IsotropicDrag isotropicDrag = new IsotropicDrag(dragArea, dragCd);
406         isotropicDrag.getDragParametersDrivers().get(0).setName("Cd");
407         final TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, isotropicDrag, TimeScalesFactory.getUTC());
408 
409 
410         // After t2 = t + 4h
411         final double dragCd2 = 3.;
412         final double dt2 = 4 * 3600.;
413         final AbsoluteDate date2 = date.shiftedBy(dt2);
414         final IsotropicDrag isotropicDrag2 = new IsotropicDrag(dragArea, dragCd2);
415         isotropicDrag2.getDragParametersDrivers().get(0).setName("Cd2");
416         forceModel.addDragSensitiveValidAfter(isotropicDrag2, date2);
417 
418         // Before t3 = t - 1day
419         final double dragCd3 = 3.;
420         final double dt3 = -86400.;
421         final AbsoluteDate date3 = date.shiftedBy(dt3);
422         final IsotropicDrag isotropicDrag3 = new IsotropicDrag(dragArea, dragCd3);
423         isotropicDrag3.getDragParametersDrivers().get(0).setName("Cd3");
424         forceModel.addDragSensitiveValidBefore(isotropicDrag3, date3);
425 
426 
427         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
428 
429         // Check parameter derivatives at initial date: only "Cd" shouldn't be 0.
430         checkParameterDerivativeGradient(state, forceModel, "Cd" , 1.0e-4, 2.0e-12);
431         checkParameterDerivativeGradient(state, forceModel, "Cd2", 1.0e-4, 0.);
432         checkParameterDerivativeGradient(state, forceModel, "Cd3", 1.0e-4, 0.);
433 
434         // Check parameter derivatives after date2: only "Cd2" shouldn't be 0.
435         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, "Cd", 1.0e-4, 0.);
436         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, "Cd2", 1.0e-4, 2.0e-12);
437         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, "Cd3", 1.0e-4, 0.);
438 
439         // Check parameter derivatives after date3: only "Cd3" shouldn't be 0.
440         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, "Cd", 1.0e-4, 0.);
441         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, "Cd2", 1.0e-4, 0.);
442         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, "Cd3", 1.0e-4, 2.0e-12);
443     }
444 
445     /** Test state Jacobian computation. */
446     @Test
447     void testStateJacobianSphere()
448         {
449 
450         // Initialization
451         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
452                                              new TimeComponents(13, 59, 27.816),
453                                              TimeScalesFactory.getUTC());
454         double i     = FastMath.toRadians(98.7);
455         double omega = FastMath.toRadians(93.0);
456         double OMEGA = FastMath.toRadians(15.0 * 22.5);
457         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
458                                          0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
459                                          Constants.EIGEN5C_EARTH_MU);
460         OrbitType integrationType = OrbitType.CARTESIAN;
461         double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
462 
463 
464 
465         // Atmosphere
466         final Atmosphere atmosphere = getAtmosphere();
467 
468         // Time span drag force model init
469         double dragArea = 2.;
470         double dragCd0 = 1.;
471         double dragCd1 = 2.;
472         double dragCd2 = 3.;
473         double dt = 1. * 3600.;
474         // Build the force model
475         IsotropicDrag isoDrag0 = new IsotropicDrag(dragArea, dragCd0);
476         IsotropicDrag isoDrag1 = new IsotropicDrag(dragArea, dragCd1);
477         IsotropicDrag isoDrag2 = new IsotropicDrag(dragArea, dragCd2);
478         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, isoDrag0);
479         forceModel.addDragSensitiveValidAfter(isoDrag1, date.shiftedBy(dt));
480         forceModel.addDragSensitiveValidBefore(isoDrag2, date.shiftedBy(-dt));
481 
482         // Check state derivatives inside first IsotropicDrag model
483         NumericalPropagator propagator =
484                         new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
485                                                                                tolerances[0], tolerances[1]));
486         propagator.setOrbitType(integrationType);
487         propagator.addForceModel(forceModel);
488         SpacecraftState state0 = new SpacecraftState(orbit);
489         // Set target date to 0.5*dt to be inside 1st drag model
490         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
491         checkStateJacobian(propagator, state0, date.shiftedBy(0.5 * dt),
492                            1e3, tolerances[0], 9.2e-10);
493 
494         // Check state derivatives inside 2nd IsotropicDrag model
495         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
496                                                                             tolerances[0], tolerances[1]));
497         propagator.setOrbitType(integrationType);
498         propagator.addForceModel(forceModel);
499         // Set target date to 1.5*dt to be inside 2nd drag model
500         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
501         checkStateJacobian(propagator, state0, date.shiftedBy(1.5 * dt),
502                            1e3, tolerances[0], 6.7e-9);
503 
504         // Check state derivatives inside 3rd IsotropicDrag model
505         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
506                                                                             tolerances[0], tolerances[1]));
507         propagator.setOrbitType(integrationType);
508         propagator.addForceModel(forceModel);
509         // Set target date to *1.5*dt to be inside 3rd drag model
510         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
511         checkStateJacobian(propagator, state0, date.shiftedBy(-1.5 * dt),
512                            1e3, tolerances[0], 6.0e-9);
513     }
514 
515     /** Test that the getParameterDrivers method is working as expected
516      * on an BoxAndSolarArraySpacecraft-based TimeSpanDragForce model.
517      * Here only the drag coefficient is modeled.
518      */
519     @Test
520     void testGetParameterDriversBoxOnlyDrag() {
521 
522         // Atmosphere
523         final CelestialBody sun = CelestialBodyFactory.getSun();
524         final Atmosphere atmosphere = getAtmosphere();
525 
526         // A date
527         AbsoluteDate date = new AbsoluteDate("2000-01-01T00:00:00.000", TimeScalesFactory.getUTC());
528 
529         // One BoxAndSolarArraySpacecraft added, test with one or two "default" drivers
530         // -----------------------
531 
532         double dragCd0 = 0.;
533         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
534                                                                          sun, 20.0, Vector3D.PLUS_J,
535                                                                          dragCd0, 0.0,
536                                                                          0.7, 0.2);
537         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
538         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
539         List<ParameterDriver> drivers = forceModel.getParametersDrivers();
540         Assertions.assertEquals(1,  drivers.size());
541         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
542         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
543 
544         // Extract drag model at an arbitrary epoch and check it is the one added
545         BoxAndSolarArraySpacecraft box = (BoxAndSolarArraySpacecraft) forceModel.getDragSensitive(date);
546         drivers = box.getDragParametersDrivers();
547         Assertions.assertEquals(1, drivers.size());
548         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
549         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
550 
551         // 3 BoxAndSolarArraySpacecraft models added, with one "default" in the middle
552         // ----------------------------------------------
553         double dragCd1 = 1.;
554         double dragCd2 = 2.;
555         double dt = 120.;
556         // Build the force model
557         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
558                                                                          sun, 20.0, Vector3D.PLUS_J,
559                                                                          dragCd1, 0.0,
560                                                                          0.7, 0.2);
561         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
562                                                                          sun, 20.0, Vector3D.PLUS_J,
563                                                                          dragCd2, 0.0,
564                                                                          0.7, 0.2);
565         forceModel = new TimeSpanDragForce(atmosphere, box0);
566         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
567         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
568 
569         // Extract the drivers and check their values and names
570         drivers = forceModel.getParametersDrivers();
571         Assertions.assertEquals(3,  drivers.size());
572         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
573         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_BEFORE + date.shiftedBy(-dt).toString(utc),
574                                 drivers.get(0).getName());
575 
576         Assertions.assertEquals(1.0,  drivers.get(1).getValue(), 0.);
577         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(1).getName());
578 
579         Assertions.assertEquals(1.0,  drivers.get(2).getValue(), 0.);
580         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_AFTER + date.shiftedBy(+dt).toString(utc),
581                                 drivers.get(2).getName());
582 
583         // Check the models at dates
584         // Cd0 model
585         double eps = 1.e-14;
586         Assertions.assertEquals(box0, forceModel.getDragSensitive(date));
587         Assertions.assertEquals(box0, forceModel.getDragSensitive(date.shiftedBy(-dt)));
588         Assertions.assertEquals(box0, forceModel.getDragSensitive(date.shiftedBy(+dt - eps)));
589         // Cd2 model
590         Assertions.assertEquals(box2, forceModel.getDragSensitive(date.shiftedBy(-dt - eps)));
591         Assertions.assertEquals(box2, forceModel.getDragSensitive(date.shiftedBy(-dt - 86400.)));
592         // Cd1 model
593         Assertions.assertEquals(box1, forceModel.getDragSensitive(date.shiftedBy(+dt)));
594         Assertions.assertEquals(box1, forceModel.getDragSensitive(date.shiftedBy(+dt + 86400.)));
595 
596         // Add a custom-named driver
597         // ----------------
598         double dragCd3 = 3.;
599         BoxAndSolarArraySpacecraft box3 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
600                                                                          sun, 20.0, Vector3D.PLUS_J,
601                                                                          dragCd3, 0.0,
602                                                                          0.7, 0.2);
603         box3.getDragParametersDrivers().get(0).setName("custom-Cd");
604         forceModel.addDragSensitiveValidAfter(box3, date.shiftedBy(2. * dt));
605         drivers = forceModel.getParametersDrivers();
606         Assertions.assertEquals(4,  drivers.size());
607         Assertions.assertEquals(1.0,  drivers.get(3).getValue(), 0.);
608         Assertions.assertEquals("custom-Cd", drivers.get(3).getName());
609     }
610 
611     /** Test that the getParameterDrivers method is working as expected
612      * on an BoxAndSolarArraySpacecraft-based TimeSpanDragForce model.
613      * Here both drag and lift coefficients are modeled.
614      */
615     @Test
616     void testGetParameterDriversBox() {
617 
618         // Atmosphere
619         final CelestialBody sun = CelestialBodyFactory.getSun();
620         final Atmosphere atmosphere = getAtmosphere();
621 
622         // A date
623         AbsoluteDate date = new AbsoluteDate("2000-01-01T00:00:00.000", TimeScalesFactory.getUTC());
624 
625         // One BoxAndSolarArraySpacecraft added, test with one or two "default" drivers
626         // -----------------------
627 
628         double dragCd0 = 0.;
629         double dragCl0 = 0.;
630         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
631                                                                          sun, 20.0, Vector3D.PLUS_J,
632                                                                          dragCd0, dragCl0,
633                                                                          0.7, 0.2);
634         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
635         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
636         List<ParameterDriver> drivers = forceModel.getParametersDrivers();
637         Assertions.assertEquals(1,  drivers.size());
638         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
639         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
640 
641         // Extract drag model at an arbitrary epoch and check it is the one added
642         BoxAndSolarArraySpacecraft box = (BoxAndSolarArraySpacecraft) forceModel.getDragSensitive(date);
643         drivers = box.getDragParametersDrivers();
644         Assertions.assertEquals(1, drivers.size());
645         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
646         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(0).getName());
647 
648         // 3 BoxAndSolarArraySpacecraft models added, with one "default" in the middle
649         // ----------------------------------------------
650         double dragCd1 = 1.;
651         double dragCl1 = 0.1;
652         double dragCd2 = 2.;
653         double dragCl2 = 0.2;
654         double dt = 120.;
655         // Build the force model
656         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
657                                                                          sun, 20.0, Vector3D.PLUS_J,
658                                                                          dragCd1, dragCl1,
659                                                                          0.7, 0.2);
660         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
661                                                                          sun, 20.0, Vector3D.PLUS_J,
662                                                                          dragCd2, dragCl2,
663                                                                          0.7, 0.2);
664         forceModel = new TimeSpanDragForce(atmosphere, box0);
665         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
666         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
667 
668         // Extract the drivers and check their values and names
669         drivers = forceModel.getParametersDrivers();
670         Assertions.assertEquals(3,  drivers.size());
671         Assertions.assertEquals(1.0,  drivers.get(0).getValue(), 0.);
672         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_BEFORE + date.shiftedBy(-dt).toString(utc),
673                                 drivers.get(0).getName());
674 
675         Assertions.assertEquals(1.0,  drivers.get(1).getValue(), 0.);
676         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR,  drivers.get(1).getName());
677 
678         Assertions.assertEquals(1.0,  drivers.get(2).getValue(), 0.);
679         Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR + TimeSpanDragForce.DATE_AFTER + date.shiftedBy(+dt).toString(utc),
680                                 drivers.get(2).getName());
681 
682         // Check the models at dates
683         // Cd0 model
684         double eps = 1.e-14;
685         Assertions.assertEquals(box0, forceModel.getDragSensitive(date));
686         Assertions.assertEquals(box0, forceModel.getDragSensitive(date.shiftedBy(-dt)));
687         Assertions.assertEquals(box0, forceModel.getDragSensitive(date.shiftedBy(+dt - eps)));
688         // Cd2 model
689         Assertions.assertEquals(box2, forceModel.getDragSensitive(date.shiftedBy(-dt - eps)));
690         Assertions.assertEquals(box2, forceModel.getDragSensitive(date.shiftedBy(-dt - 86400.)));
691         // Cd1 model
692         Assertions.assertEquals(box1, forceModel.getDragSensitive(date.shiftedBy(+dt)));
693         Assertions.assertEquals(box1, forceModel.getDragSensitive(date.shiftedBy(+dt + 86400.)));
694 
695 
696         // Add a custom-named driver
697         // ----------------
698         double dragCd3 = 3.;
699         double dragCl3 = 0.3;
700         BoxAndSolarArraySpacecraft box3 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
701                                                                          sun, 20.0, Vector3D.PLUS_J,
702                                                                          dragCd3, dragCl3,
703                                                                          0.7, 0.2);
704         box3.getDragParametersDrivers().get(0).setName("custom-Cd");
705         forceModel.addDragSensitiveValidAfter(box3, date.shiftedBy(2. * dt));
706         drivers = forceModel.getParametersDrivers();
707         Assertions.assertEquals(4,  drivers.size());
708         Assertions.assertEquals(1.0,  drivers.get(3).getValue(), 0.);
709         Assertions.assertEquals("custom-Cd", drivers.get(3).getName());
710     }
711 
712     /** Test parameter derivatives for an BoxAndSolarArraySpacecraft TimeSpanDragForce.
713      *  This test is more or less a copy of the same one in DragForceTest class
714      *  with addition of several BoxAndSolarArraySpacecraft models valid at different dates
715      *  to test that the different parameters' derivatives are computed correctly.
716      */
717     @Test
718     void testParametersDerivativesBox() {
719 
720         // Low Earth orbit definition (about 360km altitude)
721         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
722         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
723         final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
724         final SpacecraftState state =
725                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
726                                                        FramesFactory.getGCRF(),
727                                                        date,
728                                                        Constants.EIGEN5C_EARTH_MU));
729 
730         // Atmosphere
731         final CelestialBody sun = CelestialBodyFactory.getSun();
732         final Atmosphere atmosphere = getAtmosphere();
733 
734         // Initialize force model (first coef is valid at all epochs)
735         final double dragCd  = 1.;
736         final double dragCl  = 0.1;
737         final BoxAndSolarArraySpacecraft box = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
738                                                                                sun, 20.0, Vector3D.PLUS_J,
739                                                                                dragCd, dragCl, 0.7, 0.2);
740         final TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box);
741 
742 
743         // After t2 = 4h
744         final double dragCd2 = 2.;
745         final double dragCl2 = 0.2;
746         final double dt2 = 4 * 3600.;
747         final AbsoluteDate date2 = date.shiftedBy(dt2);
748         final BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
749                                                                                sun, 20.0, Vector3D.PLUS_J,
750                                                                                dragCd2, dragCl2, 0.7, 0.2);
751         box2.getDragParametersDrivers().get(0).setName("factor2");
752         forceModel.addDragSensitiveValidAfter(box2, date2);
753 
754         // Before t3 = 1day
755         final double dragCd3 = 3.;
756         final double dragCl3 = 0.3;
757         final double dt3 = -86400.;
758         final AbsoluteDate date3 = date.shiftedBy(dt3);
759         final BoxAndSolarArraySpacecraft box3 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
760                                                                                sun, 20.0, Vector3D.PLUS_J,
761                                                                                dragCd3, dragCl3, 0.7, 0.2);
762         box3.getDragParametersDrivers().get(0).setName("factor3");
763         forceModel.addDragSensitiveValidBefore(box3, date3);
764 
765         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
766 
767         // Check parameter derivatives at initial date: only 1st model parameter derivatives shouldn't be 0.
768         checkParameterDerivative(state, forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 2.2e-11);
769         checkParameterDerivative(state, forceModel, "factor2", 1.0e-4, 0.);
770         checkParameterDerivative(state, forceModel, "factor3", 1.0e-4, 0.);
771 
772         // Check parameter derivatives after date2: only 2nd model parameter derivatives shouldn't be 0.
773         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 0.);
774         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "factor2", 1.0e-4, 2.2e-12);
775         checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "factor3", 1.0e-4, 0.);
776 
777         // Check parameter derivatives before date3: only 3nd model parameter derivatives shouldn't be 0.
778         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 0.);
779         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "factor2", 1.0e-4, 0.);
780         checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "factor3", 1.0e-4, 2.0e-12);
781 
782     }
783 
784     /** Test parameter derivatives for an BoxAndSolarArraySpacecraft TimeSpanDragForce.
785      *  This test is more or less a copy of the same one in DragForceTest class
786      *  with addition of several BoxAndSolarArraySpacecraft models valid at different dates
787      *  to test that the different parameters' derivatives are computed correctly.
788      */
789     @Test
790     void testParametersDerivativesBoxGradient() {
791 
792         // Low Earth orbit definition (about 360km altitude)
793         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
794         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
795         final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
796         final SpacecraftState state =
797                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
798                                                        FramesFactory.getGCRF(),
799                                                        date,
800                                                        Constants.EIGEN5C_EARTH_MU));
801 
802         // Atmosphere
803         final CelestialBody sun = CelestialBodyFactory.getSun();
804         final Atmosphere atmosphere = getAtmosphere();
805 
806         // Initialize force model (first coef is valid at all epochs)
807         final double dragCd  = 1.;
808         final double dragCl  = 0.1;
809         final BoxAndSolarArraySpacecraft box = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
810                                                                                sun, 20.0, Vector3D.PLUS_J,
811                                                                                dragCd, dragCl, 0.7, 0.2);
812         final TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box);
813 
814 
815         // After t2 = 4h
816         final double dragCd2 = 2.;
817         final double dragCl2 = 0.2;
818         final double dt2 = 4 * 3600.;
819         final AbsoluteDate date2 = date.shiftedBy(dt2);
820         final BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
821                                                                                sun, 20.0, Vector3D.PLUS_J,
822                                                                                dragCd2, dragCl2, 0.7, 0.2);
823         box2.getDragParametersDrivers().get(0).setName("factor2");
824         forceModel.addDragSensitiveValidAfter(box2, date2);
825 
826         // Before t3 = 1day
827         final double dragCd3 = 3.;
828         final double dragCl3 = 0.3;
829         final double dt3 = -86400.;
830         final AbsoluteDate date3 = date.shiftedBy(dt3);
831         final BoxAndSolarArraySpacecraft box3 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
832                                                                                sun, 20.0, Vector3D.PLUS_J,
833                                                                                dragCd3, dragCl3, 0.7, 0.2);
834         box3.getDragParametersDrivers().get(0).setName("factor3");
835         forceModel.addDragSensitiveValidBefore(box3, date3);
836 
837         Assertions.assertFalse(forceModel.dependsOnPositionOnly());
838 
839         // Check parameter derivatives at initial date: only 1st model parameter derivatives shouldn't be 0.
840         checkParameterDerivativeGradient(state, forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 2.2e-11);
841         checkParameterDerivativeGradient(state, forceModel, "factor2", 1.0e-4, 0.);
842         checkParameterDerivativeGradient(state, forceModel, "factor3", 1.0e-4, 0.);
843 
844         // Check parameter derivatives after date2: only 2nd model parameter derivatives shouldn't be 0.
845         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 0.);
846         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, "factor2", 1.0e-4, 2.2e-12);
847         checkParameterDerivativeGradient(state.shiftedBy(dt2 * 1.1), forceModel, "factor3", 1.0e-4, 0.);
848 
849         // Check parameter derivatives before date3: only 3nd model parameter derivatives shouldn't be 0.
850         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 0.);
851         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, "factor2", 1.0e-4, 0.);
852         checkParameterDerivativeGradient(state.shiftedBy(dt3 * 1.1), forceModel, "factor3", 1.0e-4, 2.0e-12);
853 
854     }
855 
856     /** Test state Jacobian computation using finite differences in position
857      * and method {@link #accelerationDerivatives}
858      */
859     @Test
860     void testJacobianBoxVs80Implementation()
861         {
862 
863         // initialization
864         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
865                                              new TimeComponents(13, 59, 27.816),
866                                              TimeScalesFactory.getUTC());
867         double i     = FastMath.toRadians(98.7);
868         double omega = FastMath.toRadians(93.0);
869         double OMEGA = FastMath.toRadians(15.0 * 22.5);
870         Orbit refOrbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
871                                             0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
872                                             Constants.EIGEN5C_EARTH_MU);
873         CelestialBody sun = CelestialBodyFactory.getSun();
874 
875         // Atmosphere
876         final Atmosphere atmosphere = getAtmosphere();
877         // Time span drag force model initialization
878         double dragCd0 = 1.;
879         double dragCl0 = 0.1;
880         double dragCd1 = 2.;
881         double dragCl1 = 0.2;
882         double dragCd2 = 3.;
883         double dragCl2 = 0.3;
884         double dt = 3. * 3600.;
885 
886         // Build the force model
887         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
888                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
889         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
890                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
891         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
892                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
893         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
894         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
895         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
896 
897         // Check state derivatives inside first box model
898         Orbit orbit = refOrbit.shiftedBy(0.);
899         SpacecraftState state = new SpacecraftState(orbit,
900                                                     Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
901         checkStateJacobianVs80Implementation(state, forceModel,
902                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
903                                              5e-6, false);
904 
905         // Check state derivatives inside 2nd box model
906         orbit = refOrbit.shiftedBy(1.1 * dt);
907         state = new SpacecraftState(orbit,
908                                     Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
909         checkStateJacobianVs80Implementation(state, forceModel,
910                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
911                                              5e-6, false);
912 
913         // Check state derivatives inside 3rd box model
914         orbit = refOrbit.shiftedBy(-1.1 * dt);
915         state = new SpacecraftState(orbit,
916                                     Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
917         checkStateJacobianVs80Implementation(state, forceModel,
918                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
919                                              5e-6, false);
920 
921     }
922 
923     /** Test state Jacobian computation using finite differences in position
924      * and method {@link #accelerationDerivatives}
925      */
926     @Test
927     void testJacobianBoxVs80ImplementationGradient()
928         {
929 
930         // initialization
931         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
932                                              new TimeComponents(13, 59, 27.816),
933                                              TimeScalesFactory.getUTC());
934         double i     = FastMath.toRadians(98.7);
935         double omega = FastMath.toRadians(93.0);
936         double OMEGA = FastMath.toRadians(15.0 * 22.5);
937         Orbit refOrbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
938                                             0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
939                                             Constants.EIGEN5C_EARTH_MU);
940         CelestialBody sun = CelestialBodyFactory.getSun();
941         AttitudeProvider defaultLaw = Utils.defaultLaw();
942 
943         // Atmosphere
944         final Atmosphere atmosphere = getAtmosphere();
945         // Time span drag force model initialization
946         double dragCd0 = 1.;
947         double dragCl0 = 0.1;
948         double dragCd1 = 2.;
949         double dragCl1 = 0.2;
950         double dragCd2 = 3.;
951         double dragCl2 = 0.3;
952         double dt = 3. * 3600.;
953 
954         // Build the force model
955         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
956                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
957         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
958                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
959         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
960                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
961         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
962         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
963         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
964 
965         // Check state derivatives inside first box model
966         Orbit orbit = refOrbit.shiftedBy(0.);
967         SpacecraftState state = new SpacecraftState(orbit,
968                                                     Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
969         checkStateJacobianVs80ImplementationGradient(state, forceModel,
970                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
971                                              5e-6, false);
972 
973         // Check state derivatives inside 2nd box model
974         orbit = refOrbit.shiftedBy(1.1 * dt);
975         state = new SpacecraftState(orbit,
976                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
977         checkStateJacobianVs80ImplementationGradient(state, forceModel,
978                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
979                                              5e-6, false);
980 
981         // Check state derivatives inside 3rd box model
982         orbit = refOrbit.shiftedBy(-1.1 * dt);
983         state = new SpacecraftState(orbit,
984                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
985         checkStateJacobianVs80ImplementationGradient(state, forceModel,
986                                              new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
987                                              5e-6, false);
988 
989     }
990 
991     /** Test state Jacobian computation using finite differences once again.
992      * This time the differentiation is made using built-in numerical differentiation method.
993      */
994     @Test
995     void testJacobianBoxVsFiniteDifferences()
996         {
997 
998         // initialization
999         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
1000                                              new TimeComponents(13, 59, 27.816),
1001                                              TimeScalesFactory.getUTC());
1002         double i     = FastMath.toRadians(98.7);
1003         double omega = FastMath.toRadians(93.0);
1004         double OMEGA = FastMath.toRadians(15.0 * 22.5);
1005         Orbit refOrbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1006                                             0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
1007                                             Constants.EIGEN5C_EARTH_MU);
1008         CelestialBody sun = CelestialBodyFactory.getSun();
1009         AttitudeProvider defaultLaw = Utils.defaultLaw();
1010 
1011         // Atmosphere
1012         final Atmosphere atmosphere = getAtmosphere();
1013 
1014         // Time span drag force model initialization
1015         double dragCd0 = 1.;
1016         double dragCl0 = 0.1;
1017         double dragCd1 = 2.;
1018         double dragCl1 = 0.2;
1019         double dragCd2 = 3.;
1020         double dragCl2 = 0.3;
1021         double dt = 3. * 3600.;
1022 
1023         // Build the force model
1024         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1025                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1026         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1027                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1028         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1029                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1030         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1031         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
1032         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
1033 
1034         // Check state derivatives inside first box model
1035         Orbit orbit = refOrbit.shiftedBy(0.);
1036         SpacecraftState state = new SpacecraftState(orbit,
1037                                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1038         checkStateJacobianVsFiniteDifferences(state, forceModel, defaultLaw, 1.0, 5.0e-6, false);
1039 
1040         // Check state derivatives inside 2nd box model
1041         orbit = refOrbit.shiftedBy(1.1 * dt);
1042         state = new SpacecraftState(orbit,
1043                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1044         checkStateJacobianVsFiniteDifferences(state, forceModel, defaultLaw, 1.0, 5.0e-6, false);
1045 
1046         // Check state derivatives inside 3rd box model
1047         orbit = refOrbit.shiftedBy(-1.1 * dt);
1048         state = new SpacecraftState(orbit,
1049                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1050         checkStateJacobianVsFiniteDifferences(state, forceModel, defaultLaw, 1.0, 6.0e-6, false);
1051     }
1052 
1053     /** Test state Jacobian computation using finite differences once again.
1054      * This time the differentiation is made using built-in numerical differentiation method.
1055      */
1056     @Test
1057     void testJacobianBoxVsFiniteDifferencesGradient()
1058         {
1059 
1060         // initialization
1061         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
1062                                              new TimeComponents(13, 59, 27.816),
1063                                              TimeScalesFactory.getUTC());
1064         double i     = FastMath.toRadians(98.7);
1065         double omega = FastMath.toRadians(93.0);
1066         double OMEGA = FastMath.toRadians(15.0 * 22.5);
1067         Orbit refOrbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1068                                             0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
1069                                             Constants.EIGEN5C_EARTH_MU);
1070         CelestialBody sun = CelestialBodyFactory.getSun();
1071         AttitudeProvider defaultLaw = Utils.defaultLaw();
1072 
1073         // Atmosphere
1074         final Atmosphere atmosphere = getAtmosphere();
1075 
1076         // Time span drag force model initialization
1077         double dragCd0 = 1.;
1078         double dragCl0 = 0.1;
1079         double dragCd1 = 2.;
1080         double dragCl1 = 0.2;
1081         double dragCd2 = 3.;
1082         double dragCl2 = 0.3;
1083         double dt = 3. * 3600.;
1084 
1085         // Build the force model
1086         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1087                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1088         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1089                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1090         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1091                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1092         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1093         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
1094         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
1095 
1096         // Check state derivatives inside first box model
1097         Orbit orbit = refOrbit.shiftedBy(0.);
1098         final double dP = 1.;
1099         SpacecraftState state = new SpacecraftState(orbit,
1100                                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1101         checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, defaultLaw, dP, 5.0e-6, false);
1102 
1103         // Check state derivatives inside 2nd box model
1104         orbit = refOrbit.shiftedBy(1.1 * dt);
1105         state = new SpacecraftState(orbit,
1106                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1107         checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, defaultLaw, dP, 5.0e-6, false);
1108 
1109         // Check state derivatives inside 3rd box model
1110         orbit = refOrbit.shiftedBy(-1.1 * dt);
1111         state = new SpacecraftState(orbit,
1112                                     defaultLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
1113         checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, defaultLaw, dP, 6.0e-6, false);
1114     }
1115 
1116     /** Test state Jacobian computation. */
1117     @Test
1118     void testGlobalStateJacobianBox()
1119         {
1120 
1121         // Initialization
1122         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
1123                                              new TimeComponents(13, 59, 27.816),
1124                                              TimeScalesFactory.getUTC());
1125         double i     = FastMath.toRadians(98.7);
1126         double omega = FastMath.toRadians(93.0);
1127         double OMEGA = FastMath.toRadians(15.0 * 22.5);
1128         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
1129                                          0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
1130                                          Constants.EIGEN5C_EARTH_MU);
1131         OrbitType integrationType = OrbitType.CARTESIAN;
1132         double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
1133         CelestialBody sun = CelestialBodyFactory.getSun();
1134 
1135         // Atmosphere
1136         final Atmosphere atmosphere = getAtmosphere();
1137         // Time span drag force model init
1138         double dragCd0 = 1.;
1139         double dragCl0 = 0.1;
1140         double dragCd1 = 2.;
1141         double dragCl1 = 0.2;
1142         double dragCd2 = 3.;
1143         double dragCl2 = 0.3;
1144         double dt = 1. * 3600.;
1145 
1146         // Build the force model
1147         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1148                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1149         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1150                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1151         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1152                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1153         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1154         forceModel.addDragSensitiveValidAfter(box1, date.shiftedBy(dt));
1155         forceModel.addDragSensitiveValidBefore(box2, date.shiftedBy(-dt));
1156 
1157         // Check state derivatives inside first box model
1158         NumericalPropagator propagator =
1159                         new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
1160                                                                                tolerances[0], tolerances[1]));
1161         propagator.setOrbitType(integrationType);
1162         propagator.addForceModel(forceModel);
1163         SpacecraftState state0 = new SpacecraftState(orbit);
1164         // Set target date to 0.5*dt to be inside 1st drag model
1165         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
1166         checkStateJacobian(propagator, state0, date.shiftedBy(0.5 * dt),
1167                            1e3, tolerances[0], 1.1e-9);
1168 
1169         // Check state derivatives inside 2nd box model
1170         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
1171                                                                             tolerances[0], tolerances[1]));
1172         propagator.setOrbitType(integrationType);
1173         propagator.addForceModel(forceModel);
1174         // Set target date to 1.5*dt to be inside 2nd drag model
1175         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
1176         checkStateJacobian(propagator, state0, date.shiftedBy(1.5 * dt),
1177                            1e3, tolerances[0], 9.7e-9);
1178 
1179         // Check state derivatives inside 3rd box model
1180         propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
1181                                                                             tolerances[0], tolerances[1]));
1182         propagator.setOrbitType(integrationType);
1183         propagator.addForceModel(forceModel);
1184         // Set target date to *1.5*dt to be inside 3rd drag model
1185         // The further away we are from the initial date, the greater the checkTolerance parameter must be set
1186         checkStateJacobian(propagator, state0, date.shiftedBy(-1.5 * dt),
1187                            1e3, tolerances[0], 5.6e-9);
1188     }
1189 
1190     /** Testing if the propagation between the FieldPropagation and the propagation is equivalent.
1191      * Also testing if propagating X+dX with the propagation is equivalent to
1192      * propagation X with the FieldPropagation and then applying the Taylor
1193      * expansion of dX to the result.*/
1194     @Test
1195     public void RealFieldTest() {
1196         // Initial field Keplerian orbit
1197         // The variables are the six orbital parameters
1198         DSFactory factory = new DSFactory(6, 4);
1199         DerivativeStructure a_0 = factory.variable(0, 7e6);
1200         DerivativeStructure e_0 = factory.variable(1, 0.01);
1201         DerivativeStructure i_0 = factory.variable(2, 1.2);
1202         DerivativeStructure R_0 = factory.variable(3, 0.7);
1203         DerivativeStructure O_0 = factory.variable(4, 0.5);
1204         DerivativeStructure n_0 = factory.variable(5, 0.1);
1205 
1206         Field<DerivativeStructure> field = a_0.getField();
1207         DerivativeStructure zero = field.getZero();
1208 
1209         // Initial date = J2000 epoch
1210         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
1211 
1212         // J2000 frame
1213         Frame EME = FramesFactory.getEME2000();
1214 
1215         // Create initial field Keplerian orbit
1216         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
1217                                                                                  PositionAngleType.MEAN,
1218                                                                                  EME,
1219                                                                                  J2000,
1220                                                                                  zero.add(Constants.EIGEN5C_EARTH_MU));
1221 
1222         // Initial field and classical S/Cs
1223         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
1224         SpacecraftState iSR = initialState.toSpacecraftState();
1225 
1226         // Field integrator and classical integrator
1227         ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
1228                         new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
1229         ClassicalRungeKuttaIntegrator RIntegrator =
1230                         new ClassicalRungeKuttaIntegrator(6);
1231         OrbitType type = OrbitType.EQUINOCTIAL;
1232 
1233         // Field and classical numerical propagators
1234         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
1235         FNP.setOrbitType(type);
1236         FNP.setInitialState(initialState);
1237 
1238         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
1239         NP.setOrbitType(type);
1240         NP.setInitialState(iSR);
1241 
1242 
1243         // Set up force model
1244         CelestialBody sun = CelestialBodyFactory.getSun();
1245 
1246         // Atmosphere
1247         final Atmosphere atmosphere = getAtmosphere();
1248         // Time span drag force model init
1249         double dragCd0 = 1.;
1250         double dragCl0 = 0.1;
1251         double dragCd1 = 2.;
1252         double dragCl1 = 0.2;
1253         double dragCd2 = 3.;
1254         double dragCl2 = 0.3;
1255         double dt = 1000.;
1256 
1257         // Build the force model
1258         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1259                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1260         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1261                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1262         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1263                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1264         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1265         forceModel.addDragSensitiveValidAfter(box1, J2000.toAbsoluteDate().shiftedBy(dt));
1266         forceModel.addDragSensitiveValidBefore(box2, J2000.toAbsoluteDate().shiftedBy(-dt));
1267         FNP.addForceModel(forceModel);
1268         NP.addForceModel(forceModel);
1269 
1270         // Do the test
1271         // -----------
1272 
1273         // Propagate inside 1st drag model
1274         checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 0.9 * dt, NP, FNP,
1275                                   1.0e-30, 6.0e-09, 2.0e-10, 5.0e-11,
1276                                   1, false);
1277 
1278         // Propagate to 2nd drag model (reset propagator first)
1279         FNP.resetInitialState(initialState);
1280         NP.resetInitialState(iSR);
1281         checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1.1 * dt, NP, FNP,
1282                                   1.0e-30, 2.0e-08, 8.0e-11, 1.0e-10,
1283                                   1, false);
1284 
1285         // Propagate to 3rd drag model  (reset propagator first)
1286         FNP.resetInitialState(initialState);
1287         NP.resetInitialState(iSR);
1288         checkRealFieldPropagation(FKO, PositionAngleType.MEAN, -1.1 * dt, NP, FNP,
1289                                   1.0e-15, 2.0e-08, 2.0e-09, 2.0e-09,
1290                                   1, false);
1291     }
1292 
1293     /** Testing if the propagation between the FieldPropagation and the propagation is equivalent.
1294      * Also testing if propagating X+dX with the propagation is equivalent to
1295      * propagation X with the FieldPropagation and then applying the Taylor
1296      * expansion of dX to the result.*/
1297     @Test
1298     public void RealFieldGradientTest() {
1299         // Initial field Keplerian orbit
1300         // The variables are the six orbital parameters
1301         final int freeParameters = 6;
1302         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e6);
1303         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.01);
1304         Gradient i_0 = Gradient.variable(freeParameters, 2, 1.2);
1305         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
1306         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
1307         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
1308 
1309         Field<Gradient> field = a_0.getField();
1310         Gradient zero = field.getZero();
1311 
1312         // Initial date = J2000 epoch
1313         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
1314 
1315         // J2000 frame
1316         Frame EME = FramesFactory.getEME2000();
1317 
1318         // Create initial field Keplerian orbit
1319         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
1320                                                                       PositionAngleType.MEAN,
1321                                                                       EME,
1322                                                                       J2000,
1323                                                                       zero.add(Constants.EIGEN5C_EARTH_MU));
1324 
1325         // Initial field and classical S/Cs
1326         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
1327         SpacecraftState iSR = initialState.toSpacecraftState();
1328 
1329         // Field integrator and classical integrator
1330         ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
1331                         new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
1332         ClassicalRungeKuttaIntegrator RIntegrator =
1333                         new ClassicalRungeKuttaIntegrator(6);
1334         OrbitType type = OrbitType.EQUINOCTIAL;
1335 
1336         // Field and classical numerical propagators
1337         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
1338         FNP.setOrbitType(type);
1339         FNP.setInitialState(initialState);
1340         FNP.setPositionAngleType(PositionAngleType.TRUE);
1341 
1342         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
1343         NP.setOrbitType(FNP.getOrbitType());
1344         NP.setPositionAngleType(FNP.getPositionAngleType());
1345         NP.setInitialState(iSR);
1346 
1347 
1348         // Set up force model
1349         CelestialBody sun = CelestialBodyFactory.getSun();
1350 
1351         // Atmosphere
1352         final Atmosphere atmosphere = getAtmosphere();
1353         // Time span drag force model init
1354         double dragCd0 = 1.;
1355         double dragCl0 = 0.1;
1356         double dragCd1 = 2.;
1357         double dragCl1 = 0.2;
1358         double dragCd2 = 3.;
1359         double dragCl2 = 0.3;
1360         double dt = 1000.;
1361 
1362         // Build the force model
1363         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1364                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1365         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1366                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1367         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1368                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1369         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1370         forceModel.addDragSensitiveValidAfter(box1, J2000.toAbsoluteDate().shiftedBy(dt));
1371         forceModel.addDragSensitiveValidBefore(box2, J2000.toAbsoluteDate().shiftedBy(-dt));
1372         FNP.addForceModel(forceModel);
1373         NP.addForceModel(forceModel);
1374 
1375         // Do the test
1376         // -----------
1377 
1378         // Propagate inside 1st drag model
1379         checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 0.9 * dt, NP, FNP,
1380                                   1.0e-30, 2.5e-02, 7.7e-2, 1.9e-4,
1381                                   1, false);
1382 
1383         // Propagate to 2nd drag model (reset propagator first)
1384         FNP.resetInitialState(initialState);
1385         NP.resetInitialState(iSR);
1386         checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 1.1 * dt, NP, FNP,
1387                                   1.0e-30, 4.4e-02, 7.6e-5, 4.1e-4,
1388                                   1, false);
1389 
1390         // Propagate to 3rd drag model  (reset propagator first)
1391         FNP.resetInitialState(initialState);
1392         NP.resetInitialState(iSR);
1393         checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, -1.1 * dt, NP, FNP,
1394                                   1.0e-8, 2.4e-02, 2.3e-04, 3.9e-04,
1395                                   1, false);
1396     }
1397 
1398     /**Same test as the previous one but not adding the ForceModel to the NumericalPropagator
1399     it is a test to validate the previous test.
1400     (to test if the ForceModel it's actually
1401     doing something in the Propagator and the FieldPropagator)*/
1402     @Test
1403     public void RealFieldExpectErrorTest() {
1404         // Initial field Keplerian orbit
1405         // The variables are the six orbital parameters
1406         DSFactory factory = new DSFactory(6, 4);
1407         DerivativeStructure a_0 = factory.variable(0, 7e6);
1408         DerivativeStructure e_0 = factory.variable(1, 0.01);
1409         DerivativeStructure i_0 = factory.variable(2, 1.2);
1410         DerivativeStructure R_0 = factory.variable(3, 0.7);
1411         DerivativeStructure O_0 = factory.variable(4, 0.5);
1412         DerivativeStructure n_0 = factory.variable(5, 0.1);
1413 
1414         Field<DerivativeStructure> field = a_0.getField();
1415         DerivativeStructure zero = field.getZero();
1416 
1417         // Initial date = J2000 epoch
1418         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
1419 
1420         // J2000 frame
1421         Frame EME = FramesFactory.getEME2000();
1422 
1423         // Create initial field Keplerian orbit
1424         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
1425                                                                                  PositionAngleType.MEAN,
1426                                                                                  EME,
1427                                                                                  J2000,
1428                                                                                  zero.add(Constants.EIGEN5C_EARTH_MU));
1429 
1430         // Initial field and classical S/Cs
1431         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
1432         SpacecraftState iSR = initialState.toSpacecraftState();
1433 
1434         // Field integrator and classical integrator
1435         ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
1436                         new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
1437         ClassicalRungeKuttaIntegrator RIntegrator =
1438                         new ClassicalRungeKuttaIntegrator(6);
1439         OrbitType type = OrbitType.EQUINOCTIAL;
1440 
1441         // Field and classical numerical propagators
1442         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
1443         FNP.setOrbitType(type);
1444         FNP.setInitialState(initialState);
1445 
1446         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
1447         NP.setOrbitType(type);
1448         NP.setInitialState(iSR);
1449 
1450 
1451         // Set up force model
1452         CelestialBody sun = CelestialBodyFactory.getSun();
1453 
1454         // Atmosphere
1455         final Atmosphere atmosphere = getAtmosphere();
1456         // Time span drag force model init
1457         double dragCd0 = 1.;
1458         double dragCl0 = 0.1;
1459         double dragCd1 = 2.;
1460         double dragCl1 = 0.2;
1461         double dragCd2 = 3.;
1462         double dragCl2 = 0.3;
1463         double dt = 1. * 1100.;
1464 
1465         // Build the force model
1466         BoxAndSolarArraySpacecraft box0 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1467                                                                          Vector3D.PLUS_J, dragCd0, dragCl0, 0.7, 0.2);
1468         BoxAndSolarArraySpacecraft box1 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1469                                                                          Vector3D.PLUS_J, dragCd1, dragCl1, 0.7, 0.2);
1470         BoxAndSolarArraySpacecraft box2 = new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, sun, 20.0,
1471                                                                          Vector3D.PLUS_J, dragCd2, dragCl2, 0.7, 0.2);
1472         TimeSpanDragForce forceModel = new TimeSpanDragForce(atmosphere, box0);
1473         forceModel.addDragSensitiveValidAfter(box1, J2000.toAbsoluteDate().shiftedBy(dt));
1474         forceModel.addDragSensitiveValidBefore(box2, J2000.toAbsoluteDate().shiftedBy(-dt));
1475         FNP.addForceModel(forceModel);
1476         // NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR   NP.addForceModel(forceModel);
1477 
1478         FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(100.);
1479         FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
1480         SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
1481         FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
1482         PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
1483 
1484         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
1485         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
1486         Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
1487     }
1488 
1489     private Atmosphere getAtmosphere() {
1490         return new HarrisPriester(CelestialBodyFactory.getSun(),
1491                 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1492                         Constants.WGS84_EARTH_FLATTENING,
1493                         FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
1494     }
1495 
1496     @BeforeEach
1497     public void setUp() {
1498         Utils.setDataRoot("regular-data");
1499         utc = TimeScalesFactory.getUTC();
1500     }
1501 }