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.gravity;
18  
19  
20  import org.hipparchus.Field;
21  import org.hipparchus.analysis.differentiation.DSFactory;
22  import org.hipparchus.analysis.differentiation.DerivativeStructure;
23  import org.hipparchus.analysis.differentiation.Gradient;
24  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
25  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
28  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
29  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
30  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
31  import org.hipparchus.ode.nonstiff.GraggBulirschStoerIntegrator;
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.LofOffset;
38  import org.orekit.bodies.CelestialBody;
39  import org.orekit.bodies.CelestialBodyFactory;
40  import org.orekit.errors.OrekitException;
41  import org.orekit.forces.AbstractLegacyForceModelTest;
42  import org.orekit.forces.ForceModel;
43  import org.orekit.frames.Frame;
44  import org.orekit.frames.FramesFactory;
45  import org.orekit.frames.LOFType;
46  import org.orekit.orbits.CartesianOrbit;
47  import org.orekit.orbits.EquinoctialOrbit;
48  import org.orekit.orbits.FieldKeplerianOrbit;
49  import org.orekit.orbits.KeplerianOrbit;
50  import org.orekit.orbits.Orbit;
51  import org.orekit.orbits.OrbitType;
52  import org.orekit.orbits.PositionAngle;
53  import org.orekit.propagation.FieldSpacecraftState;
54  import org.orekit.propagation.SpacecraftState;
55  import org.orekit.propagation.numerical.FieldNumericalPropagator;
56  import org.orekit.propagation.numerical.NumericalPropagator;
57  import org.orekit.propagation.sampling.OrekitFixedStepHandler;
58  import org.orekit.time.AbsoluteDate;
59  import org.orekit.time.DateComponents;
60  import org.orekit.time.FieldAbsoluteDate;
61  import org.orekit.time.TimeComponents;
62  import org.orekit.time.TimeScalesFactory;
63  import org.orekit.utils.Constants;
64  import org.orekit.utils.FieldPVCoordinates;
65  import org.orekit.utils.PVCoordinates;
66  
67  public class ThirdBodyAttractionTest extends AbstractLegacyForceModelTest {
68  
69      private double mu;
70  
71      @Override
72      protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
73                                                                           final AbsoluteDate date, final  Frame frame,
74                                                                           final FieldVector3D<DerivativeStructure> position,
75                                                                           final FieldVector3D<DerivativeStructure> velocity,
76                                                                           final FieldRotation<DerivativeStructure> rotation,
77                                                                           final DerivativeStructure mass)
78          {
79          try {
80              java.lang.reflect.Field bodyField = ThirdBodyAttraction.class.getDeclaredField("body");
81              bodyField.setAccessible(true);
82              CelestialBody body = (CelestialBody) bodyField.get(forceModel);
83              double gm = forceModel.getParameterDriver(body.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX).getValue();
84  
85              // compute bodies separation vectors and squared norm
86              final Vector3D centralToBody    = body.getPVCoordinates(date, frame).getPosition();
87              final double r2Central          = centralToBody.getNormSq();
88              final FieldVector3D<DerivativeStructure> satToBody = position.subtract(centralToBody).negate();
89              final DerivativeStructure r2Sat = satToBody.getNormSq();
90  
91              // compute relative acceleration
92              final FieldVector3D<DerivativeStructure> satAcc =
93                      new FieldVector3D<>(r2Sat.sqrt().multiply(r2Sat).reciprocal().multiply(gm), satToBody);
94              final Vector3D centralAcc =
95                      new Vector3D(gm / (r2Central * FastMath.sqrt(r2Central)), centralToBody);
96              return satAcc.subtract(centralAcc);
97  
98  
99          } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
100             return null;
101         }
102     }
103 
104     @Override
105     protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
106                                                                       final AbsoluteDate date, final  Frame frame,
107                                                                       final FieldVector3D<Gradient> position,
108                                                                       final FieldVector3D<Gradient> velocity,
109                                                                       final FieldRotation<Gradient> rotation,
110                                                                       final Gradient mass)
111         {
112         try {
113             java.lang.reflect.Field bodyField = ThirdBodyAttraction.class.getDeclaredField("body");
114             bodyField.setAccessible(true);
115             CelestialBody body = (CelestialBody) bodyField.get(forceModel);
116             double gm = forceModel.getParameterDriver(body.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX).getValue();
117 
118             // compute bodies separation vectors and squared norm
119             final Vector3D centralToBody    = body.getPVCoordinates(date, frame).getPosition();
120             final double r2Central          = centralToBody.getNormSq();
121             final FieldVector3D<Gradient> satToBody = position.subtract(centralToBody).negate();
122             final Gradient r2Sat = satToBody.getNormSq();
123 
124             // compute relative acceleration
125             final FieldVector3D<Gradient> satAcc =
126                     new FieldVector3D<>(r2Sat.sqrt().multiply(r2Sat).reciprocal().multiply(gm), satToBody);
127             final Vector3D centralAcc =
128                     new Vector3D(gm / (r2Central * FastMath.sqrt(r2Central)), centralToBody);
129             return satAcc.subtract(centralAcc);
130 
131 
132         } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
133             return null;
134         }
135     }
136 
137 
138     @Test(expected= OrekitException.class)
139     public void testSunContrib() {
140 
141         // initialization
142         AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 07, 01),
143                                              new TimeComponents(13, 59, 27.816),
144                                              TimeScalesFactory.getUTC());
145         Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
146                                            FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
147                                            FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3),
148                                            0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
149         double period = 2 * FastMath.PI * orbit.getA() * FastMath.sqrt(orbit.getA() / orbit.getMu());
150 
151         // set up propagator
152         NumericalPropagator calc =
153             new NumericalPropagator(new GraggBulirschStoerIntegrator(10.0, period, 0, 1.0e-5));
154         calc.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
155 
156         // set up step handler to perform checks
157         calc.setStepHandler(FastMath.floor(period), new ReferenceChecker(date) {
158             protected double hXRef(double t) {
159                 return -1.06757e-3 + 0.221415e-11 * t + 18.9421e-5 *
160                 FastMath.cos(3.9820426e-7*t) - 7.59983e-5 * FastMath.sin(3.9820426e-7*t);
161             }
162             protected double hYRef(double t) {
163                 return 1.43526e-3 + 7.49765e-11 * t + 6.9448e-5 *
164                 FastMath.cos(3.9820426e-7*t) + 17.6083e-5 * FastMath.sin(3.9820426e-7*t);
165             }
166         });
167         AbsoluteDate finalDate = date.shiftedBy(365 * period);
168         calc.setInitialState(new SpacecraftState(orbit));
169         calc.propagate(finalDate);
170 
171     }
172 
173     /**Testing if the propagation between the FieldPropagation and the propagation
174      * is equivalent.
175      * Also testing if propagating X+dX with the propagation is equivalent to
176      * propagation X with the FieldPropagation and then applying the taylor
177      * expansion of dX to the result.*/
178     @Test
179     public void RealFieldTest() {
180         DSFactory factory = new DSFactory(6, 5);
181         DerivativeStructure a_0 = factory.variable(0, 7e7);
182         DerivativeStructure e_0 = factory.variable(1, 0.4);
183         DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
184         DerivativeStructure R_0 = factory.variable(3, 0.7);
185         DerivativeStructure O_0 = factory.variable(4, 0.5);
186         DerivativeStructure n_0 = factory.variable(5, 0.1);
187         DerivativeStructure mu  = factory.constant(Constants.EIGEN5C_EARTH_MU);
188 
189         Field<DerivativeStructure> field = a_0.getField();
190 
191         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
192 
193         Frame EME = FramesFactory.getEME2000();
194 
195         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
196                                                                                  PositionAngle.MEAN,
197                                                                                  EME,
198                                                                                  J2000,
199                                                                                  mu);
200 
201         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
202 
203         SpacecraftState iSR = initialState.toSpacecraftState();
204         OrbitType type = OrbitType.KEPLERIAN;
205         double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
206 
207 
208         AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
209                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
210         integrator.setInitialStepSize(60);
211         AdaptiveStepsizeIntegrator RIntegrator =
212                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
213         RIntegrator.setInitialStepSize(60);
214 
215         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
216         FNP.setOrbitType(type);
217         FNP.setInitialState(initialState);
218 
219         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
220         NP.setOrbitType(type);
221         NP.setInitialState(iSR);
222 
223         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
224 
225         FNP.addForceModel(forceModel);
226         NP.addForceModel(forceModel);
227         
228         // Do the test
229         checkRealFieldPropagation(FKO, PositionAngle.MEAN, 1005., NP, FNP,
230                                   1.0e-16, 5.0e-10, 3.0e-11, 3.0e-10,
231                                   1, false);
232     }
233 
234     /**Testing if the propagation between the FieldPropagation and the propagation
235      * is equivalent.
236      * Also testing if propagating X+dX with the propagation is equivalent to
237      * propagation X with the FieldPropagation and then applying the taylor
238      * expansion of dX to the result.*/
239     @Test
240     public void RealFieldGradientTest() {
241         int freeParameters = 6;
242         Gradient a_0 = Gradient.variable(freeParameters, 0, 7e7);
243         Gradient e_0 = Gradient.variable(freeParameters, 1, 0.4);
244         Gradient i_0 = Gradient.variable(freeParameters, 2, 85 * FastMath.PI / 180);
245         Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
246         Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
247         Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
248         Gradient mu  = Gradient.constant(freeParameters, Constants.EIGEN5C_EARTH_MU);
249 
250         Field<Gradient> field = a_0.getField();
251 
252         FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
253 
254         Frame EME = FramesFactory.getEME2000();
255 
256         FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
257                                                                       PositionAngle.MEAN,
258                                                                       EME,
259                                                                       J2000,
260                                                                       mu);
261 
262         FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
263 
264         SpacecraftState iSR = initialState.toSpacecraftState();
265         OrbitType type = OrbitType.KEPLERIAN;
266         double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
267 
268 
269         AdaptiveStepsizeFieldIntegrator<Gradient> integrator =
270                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
271         integrator.setInitialStepSize(60);
272         AdaptiveStepsizeIntegrator RIntegrator =
273                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
274         RIntegrator.setInitialStepSize(60);
275 
276         FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
277         FNP.setOrbitType(type);
278         FNP.setInitialState(initialState);
279 
280         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
281         NP.setOrbitType(type);
282         NP.setInitialState(iSR);
283 
284         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
285 
286         FNP.addForceModel(forceModel);
287         NP.addForceModel(forceModel);
288         
289         // Do the test
290         checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 1005., NP, FNP,
291                                   1.0e-16, 1.3e-2, 2.9e-4, 1.4e-3,
292                                   1, false);
293     }
294 
295     /**Same test as the previous one but not adding the ForceModel to the NumericalPropagator
296     it is a test to validate the previous test.
297     (to test if the ForceModel it's actually
298     doing something in the Propagator and the FieldPropagator)*/
299     @Test
300     public void RealFieldExpectErrorTest() {
301         DSFactory factory = new DSFactory(6, 5);
302         DerivativeStructure a_0 = factory.variable(0, 7e7);
303         DerivativeStructure e_0 = factory.variable(1, 0.4);
304         DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
305         DerivativeStructure R_0 = factory.variable(3, 0.7);
306         DerivativeStructure O_0 = factory.variable(4, 0.5);
307         DerivativeStructure n_0 = factory.variable(5, 0.1);
308 
309         Field<DerivativeStructure> field = a_0.getField();
310         DerivativeStructure zero = field.getZero();
311 
312         FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
313 
314         Frame EME = FramesFactory.getEME2000();
315 
316         FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
317                                                                                  PositionAngle.MEAN,
318                                                                                  EME,
319                                                                                  J2000,
320                                                                                  zero.add(Constants.EIGEN5C_EARTH_MU));
321 
322         FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
323 
324         SpacecraftState iSR = initialState.toSpacecraftState();
325         OrbitType type = OrbitType.KEPLERIAN;
326         double[][] tolerance = NumericalPropagator.tolerances(0.001, FKO.toOrbit(), type);
327 
328 
329         AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
330                         new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
331         integrator.setInitialStepSize(60);
332         AdaptiveStepsizeIntegrator RIntegrator =
333                         new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
334         RIntegrator.setInitialStepSize(60);
335 
336         FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
337         FNP.setOrbitType(type);
338         FNP.setInitialState(initialState);
339 
340         NumericalPropagator NP = new NumericalPropagator(RIntegrator);
341         NP.setOrbitType(type);
342         NP.setInitialState(iSR);
343 
344         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
345 
346         FNP.addForceModel(forceModel);
347      //NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR   NP.addForceModel(forceModel);
348 
349         FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
350         FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
351         SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
352         FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
353         PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
354 
355         Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
356         Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
357         Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
358     }
359     @Test
360     public void testMoonContrib() {
361 
362         // initialization
363         AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 07, 01),
364                                              new TimeComponents(13, 59, 27.816),
365                                              TimeScalesFactory.getUTC());
366         Orbit orbit =
367             new EquinoctialOrbit(42164000, 10e-3, 10e-3,
368                                  FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
369                                  FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3),
370                                  0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
371         double period = 2 * FastMath.PI * orbit.getA() * FastMath.sqrt(orbit.getA() / orbit.getMu());
372 
373         // set up propagator
374         NumericalPropagator calc =
375             new NumericalPropagator(new GraggBulirschStoerIntegrator(10.0, period, 0, 1.0e-5));
376         calc.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
377 
378         // set up step handler to perform checks
379         calc.setStepHandler(FastMath.floor(period), new ReferenceChecker(date) {
380             protected double hXRef(double t) {
381                 return  -0.000906173 + 1.93933e-11 * t +
382                          1.0856e-06  * FastMath.cos(5.30637e-05 * t) -
383                          1.22574e-06 * FastMath.sin(5.30637e-05 * t);
384             }
385             protected double hYRef(double t) {
386                 return 0.00151973 + 1.88991e-10 * t -
387                        1.25972e-06  * FastMath.cos(5.30637e-05 * t) -
388                        1.00581e-06 * FastMath.sin(5.30637e-05 * t);
389             }
390         });
391         AbsoluteDate finalDate = date.shiftedBy(31 * period);
392         calc.setInitialState(new SpacecraftState(orbit));
393         calc.propagate(finalDate);
394 
395     }
396 
397     private static abstract class ReferenceChecker implements OrekitFixedStepHandler {
398 
399         private final AbsoluteDate reference;
400 
401         protected ReferenceChecker(AbsoluteDate reference) {
402             this.reference = reference;
403         }
404 
405         public void handleStep(SpacecraftState currentState) {
406             double t = currentState.getDate().durationFrom(reference);
407             Assert.assertEquals(hXRef(t), currentState.getHx(), 1e-4);
408             Assert.assertEquals(hYRef(t), currentState.getHy(), 1e-4);
409         }
410 
411         protected abstract double hXRef(double t);
412 
413         protected abstract double hYRef(double t);
414 
415     }
416 
417     @Test
418     public void testParameterDerivative() {
419 
420         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
421         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
422         final SpacecraftState state =
423                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
424                                                        FramesFactory.getGCRF(),
425                                                        new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
426                                                        Constants.EIGEN5C_EARTH_MU));
427 
428         final CelestialBody moon = CelestialBodyFactory.getMoon();
429         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
430         Assert.assertTrue(forceModel.dependsOnPositionOnly());
431         final String name = moon.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX;
432         checkParameterDerivative(state, forceModel, name, 1.0, 7.0e-15);
433 
434     }
435 
436     @Test
437     public void testParameterDerivativeGradient() {
438 
439         final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
440         final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
441         final SpacecraftState state =
442                 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
443                                                        FramesFactory.getGCRF(),
444                                                        new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
445                                                        Constants.EIGEN5C_EARTH_MU));
446 
447         final CelestialBody moon = CelestialBodyFactory.getMoon();
448         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
449         Assert.assertTrue(forceModel.dependsOnPositionOnly());
450         final String name = moon.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX;
451         checkParameterDerivativeGradient(state, forceModel, name, 1.0, 7.0e-15);
452 
453     }
454 
455     @Test
456     public void testJacobianVs80Implementation() {
457         // initialization
458         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
459                                              new TimeComponents(13, 59, 27.816),
460                                              TimeScalesFactory.getUTC());
461         double i     = FastMath.toRadians(98.7);
462         double omega = FastMath.toRadians(93.0);
463         double OMEGA = FastMath.toRadians(15.0 * 22.5);
464         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
465                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
466                                          Constants.EIGEN5C_EARTH_MU);
467         final CelestialBody moon = CelestialBodyFactory.getMoon();
468         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
469         checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
470                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
471                                              1.0e-50, false);
472     }
473 
474     @Test
475     public void testJacobianVs80ImplementationGradient() {
476         // initialization
477         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
478                                              new TimeComponents(13, 59, 27.816),
479                                              TimeScalesFactory.getUTC());
480         double i     = FastMath.toRadians(98.7);
481         double omega = FastMath.toRadians(93.0);
482         double OMEGA = FastMath.toRadians(15.0 * 22.5);
483         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
484                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
485                                          Constants.EIGEN5C_EARTH_MU);
486         final CelestialBody moon = CelestialBodyFactory.getMoon();
487         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
488         checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
489                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
490                                              1.0e-15, false);
491     }
492 
493     @Test
494     public void testGlobalStateJacobian()
495         {
496 
497         // initialization
498         AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
499                                              new TimeComponents(13, 59, 27.816),
500                                              TimeScalesFactory.getUTC());
501         double i     = FastMath.toRadians(98.7);
502         double omega = FastMath.toRadians(93.0);
503         double OMEGA = FastMath.toRadians(15.0 * 22.5);
504         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
505                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
506                                          Constants.EIGEN5C_EARTH_MU);
507         OrbitType integrationType = OrbitType.CARTESIAN;
508         double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
509 
510         NumericalPropagator propagator =
511                 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
512                                                                        tolerances[0], tolerances[1]));
513         propagator.setOrbitType(integrationType);
514         final CelestialBody moon = CelestialBodyFactory.getMoon();
515         final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
516         propagator.addForceModel(forceModel);
517         SpacecraftState state0 = new SpacecraftState(orbit);
518 
519         checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
520                            1e4, tolerances[0], 2.0e-9);
521 
522     }
523 
524     @Before
525     public void setUp() {
526         mu = 3.986e14;
527         Utils.setDataRoot("regular-data");
528     }
529 
530 }