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