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