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