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