1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.gravity;
18
19 import org.hipparchus.Field;
20 import org.hipparchus.analysis.differentiation.DSFactory;
21 import org.hipparchus.analysis.differentiation.DerivativeStructure;
22 import org.hipparchus.analysis.differentiation.Gradient;
23 import org.hipparchus.analysis.differentiation.GradientField;
24 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.linear.MatrixUtils;
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.jupiter.api.Assertions;
34 import org.junit.jupiter.api.BeforeEach;
35 import org.junit.jupiter.api.DisplayName;
36 import org.junit.jupiter.api.Test;
37 import org.mockito.Mockito;
38 import org.orekit.Utils;
39 import org.orekit.attitudes.LofOffset;
40 import org.orekit.bodies.CelestialBody;
41 import org.orekit.bodies.CelestialBodyFactory;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.forces.AbstractLegacyForceModelTest;
44 import org.orekit.forces.ForceModel;
45 import org.orekit.frames.Frame;
46 import org.orekit.frames.FramesFactory;
47 import org.orekit.frames.LOFType;
48 import org.orekit.orbits.CartesianOrbit;
49 import org.orekit.orbits.EquinoctialOrbit;
50 import org.orekit.orbits.FieldKeplerianOrbit;
51 import org.orekit.orbits.KeplerianOrbit;
52 import org.orekit.orbits.Orbit;
53 import org.orekit.orbits.OrbitType;
54 import org.orekit.orbits.PositionAngleType;
55 import org.orekit.propagation.FieldSpacecraftState;
56 import org.orekit.propagation.SpacecraftState;
57 import org.orekit.propagation.numerical.FieldNumericalPropagator;
58 import org.orekit.propagation.numerical.NumericalPropagator;
59 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
60 import org.orekit.time.AbsoluteDate;
61 import org.orekit.time.DateComponents;
62 import org.orekit.time.FieldAbsoluteDate;
63 import org.orekit.time.TimeComponents;
64 import org.orekit.time.TimeScalesFactory;
65 import org.orekit.utils.Constants;
66 import org.orekit.utils.FieldPVCoordinates;
67 import org.orekit.utils.PVCoordinates;
68
69 public class ThirdBodyAttractionTest extends AbstractLegacyForceModelTest {
70
71 private double mu;
72
73 @Override
74 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
75 final FieldSpacecraftState<DerivativeStructure> state) {
76 try {
77 final AbsoluteDate date = state.getDate().toAbsoluteDate();
78 final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
79 java.lang.reflect.Field bodyField = ThirdBodyAttraction.class.getDeclaredField("body");
80 bodyField.setAccessible(true);
81 CelestialBody body = (CelestialBody) bodyField.get(forceModel);
82 double gm = forceModel.
83 getParameterDriver(body.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX).
84 getValue(date);
85
86
87 final Vector3D centralToBody = body.getPosition(date, state.getFrame());
88 final double r2Central = centralToBody.getNormSq();
89 final FieldVector3D<DerivativeStructure> satToBody = position.subtract(centralToBody).negate();
90 final DerivativeStructure r2Sat = satToBody.getNormSq();
91
92
93 final FieldVector3D<DerivativeStructure> satAcc =
94 new FieldVector3D<>(r2Sat.sqrt().multiply(r2Sat).reciprocal().multiply(gm), satToBody);
95 final Vector3D centralAcc =
96 new Vector3D(gm / (r2Central * FastMath.sqrt(r2Central)), centralToBody);
97 return satAcc.subtract(centralAcc);
98
99
100 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
101 return null;
102 }
103 }
104
105 @Override
106 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
107 final FieldSpacecraftState<Gradient> state) {
108 try {
109 final AbsoluteDate date = state.getDate().toAbsoluteDate();
110 final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
111 java.lang.reflect.Field bodyField = ThirdBodyAttraction.class.getDeclaredField("body");
112 bodyField.setAccessible(true);
113 CelestialBody body = (CelestialBody) bodyField.get(forceModel);
114 double gm = forceModel.
115 getParameterDriver(body.getName() + ThirdBodyAttraction.ATTRACTION_COEFFICIENT_SUFFIX).
116 getValue(date);
117
118
119 final Vector3D centralToBody = body.getPosition(date, state.getFrame());
120 final double r2Central = centralToBody.getNormSq();
121 final FieldVector3D<Gradient> satToBody = position.subtract(centralToBody).negate();
122 final Gradient r2Sat = satToBody.getNormSq();
123
124
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
139 public void testSunContrib() {
140 Assertions.assertThrows(OrekitException.class, () -> {
141
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, PositionAngleType.TRUE, FramesFactory.getEME2000(), date, mu);
149 double period = 2 * FastMath.PI * orbit.getA() * FastMath.sqrt(orbit.getA() / orbit.getMu());
150
151
152 NumericalPropagator calc =
153 new NumericalPropagator(new GraggBulirschStoerIntegrator(10.0, period, 0, 1.0e-5));
154 calc.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
155
156
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
174
175
176
177
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 PositionAngleType.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
229 checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
230 1.0e-16, 5.0e-10, 3.0e-11, 3.0e-10,
231 1, false);
232 }
233
234
235
236
237
238
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 PositionAngleType.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
290 checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 1005., NP, FNP,
291 1.0e-16, 1.3e-2, 2.9e-4, 1.4e-3,
292 1, false);
293 }
294
295
296
297
298
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 PositionAngleType.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
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 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
356 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
357 Assertions.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
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, PositionAngleType.TRUE, FramesFactory.getEME2000(), date, mu);
371 double period = 2 * FastMath.PI * orbit.getA() * FastMath.sqrt(orbit.getA() / orbit.getMu());
372
373
374 NumericalPropagator calc =
375 new NumericalPropagator(new GraggBulirschStoerIntegrator(10.0, period, 0, 1.0e-5));
376 calc.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
377
378
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 Assertions.assertEquals(hXRef(t), currentState.getHx(), 1e-4);
408 Assertions.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 Assertions.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 Assertions.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
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, PositionAngleType.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
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, PositionAngleType.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
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, PositionAngleType.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 @Test
525 @DisplayName("Test that acceleration derivatives with respect to absolute date are not equal to zero.")
526 public void testIssue1070() {
527
528
529 final int freeParameters = 1;
530 final GradientField field = GradientField.getField(freeParameters);
531 final Gradient zero = field.getZero();
532 final Gradient variable = Gradient.variable(freeParameters, 0, 0.);
533 final FieldAbsoluteDate<Gradient> fieldAbsoluteDate = new FieldAbsoluteDate<>(field, AbsoluteDate.ARBITRARY_EPOCH).
534 shiftedBy(variable);
535
536
537 @SuppressWarnings("unchecked")
538 final FieldSpacecraftState<Gradient> stateMock = Mockito.mock(FieldSpacecraftState.class);
539 Mockito.when(stateMock.getDate()).thenReturn(fieldAbsoluteDate);
540 Mockito.when(stateMock.getPosition()).thenReturn(new FieldVector3D<>(zero, zero));
541 Mockito.when(stateMock.getFrame()).thenReturn(FramesFactory.getGCRF());
542
543
544 final CelestialBody moon = CelestialBodyFactory.getMoon();
545 final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon);
546
547
548 final Gradient gm = zero.add(moon.getGM());
549 final FieldVector3D<Gradient> accelerationVector = forceModel.acceleration(stateMock, new Gradient[] { gm });
550
551
552 final double[] derivatives = accelerationVector.getNormSq().getGradient();
553 Assertions.assertNotEquals(0., MatrixUtils.createRealVector(derivatives).getNorm());
554 }
555
556 @BeforeEach
557 public void setUp() {
558 mu = 3.986e14;
559 Utils.setDataRoot("regular-data");
560 }
561
562 }