1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.maneuvers.jacobians;
18
19 import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
20 import org.hipparchus.geometry.euclidean.threed.Rotation;
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.hipparchus.linear.RealMatrix;
23 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.MathUtils;
27 import org.junit.jupiter.api.Assertions;
28 import org.junit.jupiter.api.BeforeEach;
29 import org.junit.jupiter.api.Test;
30 import org.orekit.Utils;
31 import org.orekit.attitudes.AttitudeProvider;
32 import org.orekit.attitudes.FrameAlignedProvider;
33 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
34 import org.orekit.forces.gravity.NewtonianAttraction;
35 import org.orekit.forces.gravity.potential.GravityFieldFactory;
36 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
37 import org.orekit.forces.maneuvers.Maneuver;
38 import org.orekit.forces.maneuvers.propulsion.BasicConstantThrustPropulsionModel;
39 import org.orekit.forces.maneuvers.propulsion.PropulsionModel;
40 import org.orekit.forces.maneuvers.trigger.DateBasedManeuverTriggers;
41 import org.orekit.frames.FramesFactory;
42 import org.orekit.orbits.KeplerianOrbit;
43 import org.orekit.orbits.Orbit;
44 import org.orekit.orbits.OrbitType;
45 import org.orekit.orbits.PositionAngleType;
46 import org.orekit.propagation.*;
47 import org.orekit.propagation.numerical.NumericalPropagator;
48 import org.orekit.propagation.sampling.MultiSatStepHandler;
49 import org.orekit.propagation.sampling.OrekitStepInterpolator;
50 import org.orekit.time.AbsoluteDate;
51 import org.orekit.time.DateComponents;
52 import org.orekit.time.TimeComponents;
53 import org.orekit.time.TimeScalesFactory;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.IERSConventions;
56
57 import java.util.ArrayList;
58 import java.util.Arrays;
59 import java.util.List;
60
61 public class TriggersDerivativesTest {
62
63 @Test
64 void testDerivativeWrtStartTimeCartesianForward() {
65 doTestDerivativeWrtStartStopTime(true, OrbitType.CARTESIAN, true,
66 4.3e-6, 2.8e-6, 4.3e-6, 1.5e-5, 8.1e-6, 1.5e-5);
67 }
68
69 @Test
70 void testDerivativeWrtStartTimeKeplerianForward() {
71 doTestDerivativeWrtStartStopTime(true, OrbitType.KEPLERIAN, true,
72 1.5e-5, 1.5e-5, 1.5e-5, 6.8e-6, 1.1e-4, 3.2e-5);
73 }
74
75 @Test
76 void testDerivativeWrtStartTimeCartesianBackward() {
77 doTestDerivativeWrtStartStopTime(true, OrbitType.CARTESIAN, false,
78 9.7e-5, 9.7e-5, 9.7e-5, 4.7e-8, 4.7e-8, 4.7e-8);
79 }
80
81 @Test
82 void testDerivativeWrtStartTimeKeplerianBackward() {
83 doTestDerivativeWrtStartStopTime(true, OrbitType.KEPLERIAN, false,
84 6.5e-8, 6.5e-8, 6.5e-8, 1.1e-6, 1.8e-6, 6.2e-7);
85 }
86
87 @Test
88 void testDerivativeWrtStopTimeCartesianForward() {
89 doTestDerivativeWrtStartStopTime(false, OrbitType.CARTESIAN, true,
90 5.2e-10, 6.0e-8, 6.6e-11, 6.8e-12, 3.9e-11, 7.5e-12);
91 }
92
93 @Test
94 void testDerivativeWrtStopTimeKeplerianForward() {
95 doTestDerivativeWrtStartStopTime(false, OrbitType.KEPLERIAN, true,
96 1.8e-11, 9.7e-12, 3.0e-12, 2.6e-9, 2.9e-9, 1.8e-9);
97 }
98
99 @Test
100 void testDerivativeWrtStopTimeCartesianBackward() {
101 doTestDerivativeWrtStartStopTime(false, OrbitType.CARTESIAN, false,
102 9.7e-5, 9.7e-5, 9.7e-5, 1.5e-5, 2.1e-5, 1.5e-5);
103 }
104
105 @Test
106 void testDerivativeWrtStopTimeKeplerianBackward() {
107 doTestDerivativeWrtStartStopTime(false, OrbitType.KEPLERIAN, false,
108 1.5e-5, 1.5e-5, 1.5e-5, 3.8e-6, 1.2e-4, 3.0e-4);
109 }
110
111 @Test
112 void testDerivativeWrtMedianCartesianForward() {
113 doTestDerivativeWrtMedianDuration(true, OrbitType.CARTESIAN, true,
114 1.8e-5, 1.2e-5, 1.8e-5, 2.0e-2, 1.7e-2, 2.0e-2);
115 }
116
117 @Test
118 void testDerivativeWrtMedianKeplerianForward() {
119 doTestDerivativeWrtMedianDuration(true, OrbitType.KEPLERIAN, true,
120 0.095, 0.13, 0.11, 9.0e-6, 4.8e-5, 3.2e-4);
121 }
122
123 @Test
124 void testDerivativeWrtMedianCartesianBackward() {
125 doTestDerivativeWrtMedianDuration(true, OrbitType.CARTESIAN, false,
126 9.7e-5, 9.7e-5, 9.7e-5, 1.9e-2, 2.1e-2, 2.0e-2);
127 }
128
129 @Test
130 void testDerivativeWrtMedianKeplerianBackward() {
131 doTestDerivativeWrtMedianDuration(true, OrbitType.KEPLERIAN, false,
132 0.091, 0.13, 0.11, 7.4e-6, 4.7e-5, 3.4e-4);
133 }
134
135 @Test
136 void testDerivativeWrtDurationCartesianForward() {
137 doTestDerivativeWrtMedianDuration(false, OrbitType.CARTESIAN, true,
138 2.5e-6, 1.6e-6, 2.5e-6, 7.3e-6, 4.1e-6, 7.2e-6);
139 }
140
141 @Test
142 void testDerivativeWrtDurationKeplerianForward() {
143 doTestDerivativeWrtMedianDuration(false, OrbitType.KEPLERIAN, true,
144 7.2e-6, 7.2e-6, 7.2e-6, 2.5e-6, 2.5e-5, 1.5e-5);
145 }
146
147 @Test
148 void testDerivativeWrtDurationCartesianBackward() {
149 doTestDerivativeWrtMedianDuration(false, OrbitType.CARTESIAN, false,
150 9.7e-5, 9.7e-5, 9.7e-5, 7.1e-6, 1.1e-5, 7.2e-6);
151 }
152
153 @Test
154 void testDerivativeWrtDurationKeplerianBackward() {
155 doTestDerivativeWrtMedianDuration(false, OrbitType.KEPLERIAN, false,
156 7.2e-6, 7.2e-6, 7.2e-6, 2.3e-6, 2.9e-5, 3.0e-4);
157 }
158
159 private void doTestDerivativeWrtStartStopTime(final boolean start, final OrbitType orbitType, final boolean forward,
160 final double...tolerance) {
161
162 final PositionAngleType positionAngleType = PositionAngleType.TRUE;
163 final int degree = 20;
164 final double duration = 200.0;
165 final double h = 1.0;
166 final double samplingtep = 2.0;
167
168 final KeplerianOrbit initial = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(buildInitialState(buildAttitudeProvider()).getOrbit());
169 final double firingM = MathUtils.normalizeAngle(FastMath.PI, initial.getMeanAnomaly() + (forward ? FastMath.PI : -FastMath.PI));
170 final AbsoluteDate apogeeDate = initial.getDate().shiftedBy((firingM - initial.getMeanAnomaly()) / initial.getKeplerianMeanMotion());
171 final AbsoluteDate firing = apogeeDate.shiftedBy(-0.5 * duration);
172
173 final List<Propagator> propagators = new ArrayList<>();
174
175
176 for (int k = -4; k <= 4; ++k) {
177 final DateBasedManeuverTriggers trigger = start ?
178 new DateBasedManeuverTriggers("MAN_0", firing.shiftedBy(k * h), duration - k * h) :
179 new DateBasedManeuverTriggers("MAN_0", firing, duration + k * h);
180 propagators.add(buildPropagator(orbitType, positionAngleType, degree, firing, duration, trigger));
181 }
182
183
184 final NumericalPropagator autonomous = (NumericalPropagator) propagators.get(4);
185 final MatricesHarvester harvester = autonomous.setupMatricesComputation("stm", null, null);
186 autonomous.
187 getAllForceModels().
188 forEach(fm -> fm.
189 getParametersDrivers().
190 stream().
191 filter(d -> d.getName().equals(start ? "MAN_0_START" : "MAN_0_STOP") ||
192 d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
193 forEach(d -> d.setSelected(true)));
194
195 DerivativesSampler sampler = start ?
196 new DerivativesSampler(harvester, 4, null, -1, null, -1, null, -1, orbitType, positionAngleType,
197 firing, duration, h, samplingtep) :
198 new DerivativesSampler(null, -1, harvester, 4, null, -1, null, -1, orbitType, positionAngleType,
199 firing, duration, h, samplingtep);
200
201 final PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(propagators, sampler);
202 if (forward) {
203 parallelizer.propagate(firing.shiftedBy(-30 * samplingtep), firing.shiftedBy(duration + 300 * samplingtep));
204 } else {
205 parallelizer.propagate(firing.shiftedBy(duration + 30 * samplingtep), firing.shiftedBy(-300 * samplingtep));
206 }
207
208 double[] maxRelativeError = new double[tolerance.length];
209 for (final Entry entry : sampler.sample) {
210 final Result result = start ? entry.start : entry.stop;
211 for (int i = 0; i < 6; ++i) {
212 double f = result.finiteDifferences[i].getFirstDerivative();
213 double c = result.closedForm[i].getFirstDerivative();
214 maxRelativeError[i] = FastMath.max(maxRelativeError[i], FastMath.abs(f - c) / FastMath.max(1.0e-10, FastMath.abs(f)));
215 }
216 }
217
218 for (int i = 0; i < tolerance.length; ++i) {
219 Assertions.assertEquals(0.0, maxRelativeError[i], tolerance[i]);
220 }
221
222 }
223
224 private void doTestDerivativeWrtMedianDuration(final boolean median, final OrbitType orbitType, final boolean forward,
225 final double...tolerance) {
226
227 final PositionAngleType positionAngleType = PositionAngleType.TRUE;
228 final int degree = 20;
229 final double duration = 200.0;
230 final double h = 1.0;
231 final double samplingtep = 2.0;
232
233 final KeplerianOrbit initial = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(buildInitialState(buildAttitudeProvider()).getOrbit());
234 final double firingM = MathUtils.normalizeAngle(FastMath.PI, initial.getMeanAnomaly() + (forward ? FastMath.PI : -FastMath.PI));
235 final AbsoluteDate apogeeDate = initial.getDate().shiftedBy((firingM - initial.getMeanAnomaly()) / initial.getKeplerianMeanMotion());
236 final AbsoluteDate firing = apogeeDate.shiftedBy(-0.5 * duration);
237
238 final List<Propagator> propagators = new ArrayList<>();
239
240
241 for (int k = -4; k <= 4; ++k) {
242 final DateBasedManeuverTriggers triggers = median ?
243 new DateBasedManeuverTriggers("MAN_0", firing.shiftedBy(k * h), duration) :
244 new DateBasedManeuverTriggers("MAN_0", firing.shiftedBy(-0.5 * k * h), duration + k * h);
245 propagators.add(buildPropagator(orbitType, positionAngleType, degree, firing, duration, triggers));
246 }
247 for (int k = -4; k <= 4; ++k) {
248 final DateBasedManeuverTriggers triggers =
249 new DateBasedManeuverTriggers("MAN_1", firing.shiftedBy(k * h), duration - k * h);
250 propagators.add(buildPropagator(orbitType, positionAngleType, degree, firing, duration, triggers));
251 }
252 for (int k = -4; k <= 4; ++k) {
253 final DateBasedManeuverTriggers triggers =
254 new DateBasedManeuverTriggers("MAN_2", firing, duration + k * h);
255 propagators.add(buildPropagator(orbitType, positionAngleType, degree, firing, duration, triggers));
256 }
257
258
259
260
261 final NumericalPropagator autonomous = (NumericalPropagator) propagators.get(4);
262 final MatricesHarvester harvester = autonomous.setupMatricesComputation("stm-0", null, null);
263 autonomous.
264 getAllForceModels().
265 forEach(fm -> fm.
266 getParametersDrivers().
267 stream().
268 filter(d -> d.getName().equals(median ? "MAN_0_MEDIAN" : "MAN_0_DURATION") ||
269 d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
270 forEach(d -> d.setSelected(true)));
271
272 final NumericalPropagator autonomousStart = (NumericalPropagator) propagators.get(13);
273 final MatricesHarvester harvesterStart = autonomousStart.setupMatricesComputation("stm-1", null, null);
274 autonomousStart.
275 getAllForceModels().
276 forEach(fm -> fm.
277 getParametersDrivers().
278 stream().
279 filter(d -> d.getName().equals("MAN_1_START") ||
280 d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
281 forEach(d -> d.setSelected(true)));
282
283 final NumericalPropagator autonomousStop = (NumericalPropagator) propagators.get(22);
284 final MatricesHarvester harvesterStop = autonomousStop.setupMatricesComputation("stm-2", null, null);
285 autonomousStop.
286 getAllForceModels().
287 forEach(fm -> fm.
288 getParametersDrivers().
289 stream().
290 filter(d -> d.getName().equals("MAN_2_STOP") ||
291 d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
292 forEach(d -> d.setSelected(true)));
293
294 DerivativesSampler sampler = median ?
295 new DerivativesSampler(harvesterStart, 13, harvesterStop, 22, harvester, 4, null, -1,
296 orbitType, positionAngleType,
297 firing, duration, h, samplingtep) :
298 new DerivativesSampler(harvesterStart, 13, harvesterStop, 22, null, -1, harvester, 4,
299 orbitType, positionAngleType,
300 firing, duration, h, samplingtep);
301
302 final PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(propagators, sampler);
303 if (forward) {
304 parallelizer.propagate(firing.shiftedBy(-30 * samplingtep), firing.shiftedBy(duration + 300 * samplingtep));
305 } else {
306 parallelizer.propagate(firing.shiftedBy(duration + 30 * samplingtep), firing.shiftedBy(-300 * samplingtep));
307 }
308
309 double[] maxRelativeError = new double[tolerance.length];
310 for (final Entry entry : sampler.sample) {
311 final Result result = median ? entry.median : entry.duration;
312 for (int i = 0; i < 6; ++i) {
313 double f = result.finiteDifferences[i].getFirstDerivative();
314 double c = result.closedForm[i].getFirstDerivative();
315 maxRelativeError[i] = FastMath.max(maxRelativeError[i], FastMath.abs(f - c) / FastMath.max(1.0e-10, FastMath.abs(f)));
316 }
317 }
318
319 for (int i = 0; i < tolerance.length; ++i) {
320 Assertions.assertEquals(0.0, maxRelativeError[i], tolerance[i]);
321 }
322
323 }
324
325 private NumericalPropagator buildPropagator(final OrbitType orbitType, final PositionAngleType positionAngleType,
326 final int degree, final AbsoluteDate firing, final double duration,
327 final DateBasedManeuverTriggers triggers) {
328
329 final AttitudeProvider attitudeProvider = buildAttitudeProvider();
330 SpacecraftState initialState = buildInitialState(attitudeProvider);
331
332 final double isp = 318;
333 final double f = 420;
334 PropulsionModel propulsionModel = new BasicConstantThrustPropulsionModel(f, isp, Vector3D.PLUS_I, "ABM");
335
336 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(initialState.getOrbit(), orbitType);
337 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 600, tol[0], tol[1]);
338 integrator.setInitialStepSize(60);
339 final NumericalPropagator propagator = new NumericalPropagator(integrator);
340
341 propagator.setOrbitType(orbitType);
342 propagator.setPositionAngleType(positionAngleType);
343 propagator.setAttitudeProvider(attitudeProvider);
344 if (degree > 0) {
345 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
346 GravityFieldFactory.getNormalizedProvider(degree, degree)));
347 }
348 final Maneuver maneuver = new Maneuver(null, triggers, propulsionModel);
349 propagator.addForceModel(maneuver);
350 propagator.addAdditionalDataProvider(new AdditionalDataProvider<double[]>() {
351 public String getName() { return triggers.getName().concat("-acc"); }
352 public double[] getAdditionalData(SpacecraftState state) {
353 double[] parameters = Arrays.copyOfRange(maneuver.getParameters(initialState.getDate()), 0, propulsionModel.getParametersDrivers().size());
354 return new double[] {
355 propulsionModel.getAcceleration(state, state.getAttitude(), parameters).getNorm()
356 };
357 }
358 });
359 propagator.setInitialState(initialState);
360 return propagator;
361
362 }
363
364 private SpacecraftState buildInitialState(final AttitudeProvider attitudeProvider) {
365 final double mass = 2500;
366 final double a = 24396159;
367 final double e = 0.72831215;
368 final double i = FastMath.toRadians(7);
369 final double omega = FastMath.toRadians(180);
370 final double OMEGA = FastMath.toRadians(261);
371 final double lv = 0;
372
373 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 1, 1), new TimeComponents(23, 30, 00.000),
374 TimeScalesFactory.getUTC());
375 final Orbit orbit = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
376 FramesFactory.getEME2000(), initDate, Constants.EIGEN5C_EARTH_MU);
377 return new SpacecraftState(orbit, attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
378 }
379
380 private AttitudeProvider buildAttitudeProvider() {
381 final double delta = FastMath.toRadians(-7.4978);
382 final double alpha = FastMath.toRadians(351);
383 return new FrameAlignedProvider(new Rotation(new Vector3D(alpha, delta), Vector3D.PLUS_I));
384 }
385
386 private class DerivativesSampler implements MultiSatStepHandler {
387
388 final MatricesHarvester harvesterStart;
389 final int indexStart;
390 final MatricesHarvester harvesterStop;
391 final int indexStop;
392 final MatricesHarvester harvesterMedian;
393 final int indexMedian;
394 final MatricesHarvester harvesterDuration;
395 final int indexDuration;
396 final OrbitType orbitType;
397 final PositionAngleType positionAngleType;
398 final AbsoluteDate firing;
399 final double duration;
400 final double h;
401 final double samplingtep;
402 final List<Entry> sample;
403 boolean forward;
404 AbsoluteDate next;
405
406 DerivativesSampler(final MatricesHarvester harvesterStart, final int indexStart,
407 final MatricesHarvester harvesterStop, final int indexStop,
408 final MatricesHarvester harvesterMedian, final int indexMedian,
409 final MatricesHarvester harvesterDuration, final int indexDuration,
410 final OrbitType orbitType, final PositionAngleType positionAngleType,
411 final AbsoluteDate firing, final double duration,
412 final double h, final double samplingtep) {
413 this.harvesterStart = harvesterStart;
414 this.indexStart = indexStart;
415 this.harvesterStop = harvesterStop;
416 this.indexStop = indexStop;
417 this.harvesterMedian = harvesterMedian;
418 this.indexMedian = indexMedian;
419 this.harvesterDuration = harvesterDuration;
420 this.indexDuration = indexDuration;
421 this.orbitType = orbitType;
422 this.positionAngleType = positionAngleType;
423 this.firing = firing;
424 this.duration = duration;
425 this.h = h;
426 this.samplingtep = samplingtep;
427 this.sample = new ArrayList<>();
428 this.next = null;
429 }
430
431 public void init(final List<SpacecraftState> states0, final AbsoluteDate t) {
432 final AbsoluteDate t0 = states0.get(0).getDate();
433 if (t.isAfterOrEqualTo(t0)) {
434 forward = true;
435 next = t0.shiftedBy(samplingtep);
436 } else {
437 forward = false;
438 next = t0.shiftedBy(-samplingtep);
439 }
440 }
441
442 public void handleStep(final List<OrekitStepInterpolator> interpolators) {
443 final AbsoluteDate previousDate = interpolators.get(0).getPreviousState().getDate();
444 final AbsoluteDate currentDate = interpolators.get(0).getCurrentState().getDate();
445 while ( forward && (next.isAfter(previousDate) && next.isBeforeOrEqualTo(currentDate)) ||
446 !forward && (next.isBefore(previousDate) && next.isAfterOrEqualTo(currentDate))) {
447
448 if (!(tooClose(next, firing) || tooClose(next, firing.shiftedBy(duration)))) {
449 final Entry entry = new Entry(next.getDate());
450 fill(interpolators, harvesterMedian, indexMedian, entry.median);
451 fill(interpolators, harvesterStart, indexStart, entry.start);
452 fill(interpolators, harvesterStop, indexStop, entry.stop);
453 fill(interpolators, harvesterDuration, indexDuration, entry.duration);
454 sample.add(entry);
455 }
456 next = next.shiftedBy(forward ? samplingtep : -samplingtep);
457 }
458 }
459
460 private void fill(final List<OrekitStepInterpolator> interpolators, final MatricesHarvester harvester, final int index, final Result result) {
461 if (index >= 0) {
462 final double[][] o = new double[9][6];
463 for (int i = 0; i < o.length; ++i) {
464 orbitType.mapOrbitToArray(interpolators.get(index + i - 4).getInterpolatedState(next).getOrbit(),
465 positionAngleType, o[i], null);
466 }
467 final SpacecraftState centralState = interpolators.get(index).getInterpolatedState(next);
468 final RealMatrix jacobian = harvester.getParametersJacobian(centralState);
469 for (int i = 0; i < result.closedForm.length; ++i) {
470 result.closedForm[i] = new UnivariateDerivative1(o[4][i], jacobian.getEntry(i, 0));
471 result.finiteDifferences[i] = new UnivariateDerivative1(o[4][i],
472 differential8(o[0][i], o[1][i], o[2][i], o[3][i],
473 o[5][i], o[6][i], o[7][i], o[8][i],
474 h));
475 }
476 }
477 }
478
479 private boolean tooClose(final AbsoluteDate date, final AbsoluteDate discontinuity) {
480 final double maxOffset = h * 4;
481 return date.shiftedBy(-maxOffset).isBeforeOrEqualTo(discontinuity) &&
482 date.shiftedBy(+maxOffset).isAfterOrEqualTo(discontinuity);
483 }
484
485 }
486
487 private static class Entry {
488 private Result start;
489 private Result stop;
490 private Result median;
491 private Result duration;
492 Entry(final AbsoluteDate date) {
493 this.start = new Result();
494 this.stop = new Result();
495 this.median = new Result();
496 this.duration = new Result();
497 }
498 }
499
500 private static class Result {
501 private UnivariateDerivative1[] finiteDifferences;
502 private UnivariateDerivative1[] closedForm;
503 Result() {
504 this.finiteDifferences = new UnivariateDerivative1[6];
505 this.closedForm = new UnivariateDerivative1[6];
506 }
507 }
508
509 private double differential8(final double fM4h, final double fM3h, final double fM2h, final double fM1h,
510 final double fP1h, final double fP2h, final double fP3h, final double fP4h,
511 final double h) {
512
513
514 return (-3 * (fP4h - fM4h) + 32 * (fP3h - fM3h) - 168 * (fP2h - fM2h) + 672 * (fP1h - fM1h)) / (840 * h);
515
516 }
517
518 @BeforeEach
519 public void setUp() {
520 Utils.setDataRoot("orbit-determination/february-2016:potential/icgem-format");
521 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
522 }
523
524 }