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