1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.semianalytical.dsst;
18
19 import java.io.FileNotFoundException;
20 import java.io.IOException;
21 import java.io.UnsupportedEncodingException;
22 import java.text.ParseException;
23
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.FastMath;
26 import org.junit.Assert;
27 import org.junit.Before;
28 import org.junit.Test;
29 import org.orekit.Utils;
30 import org.orekit.attitudes.Attitude;
31 import org.orekit.bodies.CelestialBodyFactory;
32 import org.orekit.bodies.OneAxisEllipsoid;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitMessages;
35 import org.orekit.forces.drag.DragSensitive;
36 import org.orekit.forces.drag.IsotropicDrag;
37 import org.orekit.forces.gravity.potential.GravityFieldFactory;
38 import org.orekit.forces.gravity.potential.SHMFormatReader;
39 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
40 import org.orekit.frames.Frame;
41 import org.orekit.frames.FramesFactory;
42 import org.orekit.models.earth.atmosphere.HarrisPriester;
43 import org.orekit.orbits.EquinoctialOrbit;
44 import org.orekit.orbits.KeplerianOrbit;
45 import org.orekit.orbits.Orbit;
46 import org.orekit.orbits.OrbitType;
47 import org.orekit.orbits.PositionAngle;
48 import org.orekit.propagation.PropagationType;
49 import org.orekit.propagation.SpacecraftState;
50 import org.orekit.propagation.numerical.NumericalPropagator;
51 import org.orekit.propagation.sampling.OrekitStepHandler;
52 import org.orekit.propagation.sampling.OrekitStepInterpolator;
53 import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
54 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
55 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
56 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
57 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
58 import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
59 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
60 import org.orekit.time.AbsoluteDate;
61 import org.orekit.utils.Constants;
62 import org.orekit.utils.IERSConventions;
63 import org.orekit.utils.ParameterDriver;
64 import org.orekit.utils.ParameterDriversList;
65
66 @Deprecated
67 public class DSSTPartialDerivativesTest {
68
69 @Test
70 public void testDragParametersDerivatives() throws ParseException, IOException {
71 doTestParametersDerivatives(DragSensitive.DRAG_COEFFICIENT,
72 2.4e-3,
73 PropagationType.MEAN,
74 OrbitType.EQUINOCTIAL);
75 }
76
77 @Test
78 public void testMuParametersDerivatives() throws ParseException, IOException {
79 doTestParametersDerivatives(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
80 5.e-3,
81 PropagationType.MEAN,
82 OrbitType.EQUINOCTIAL);
83 }
84
85 private void doTestParametersDerivatives(String parameterName, double tolerance,
86 PropagationType type,
87 OrbitType... orbitTypes) {
88
89 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
90 Constants.WGS84_EARTH_FLATTENING,
91 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
92
93 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
94
95 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
96
97 DSSTForceModel drag = new DSSTAtmosphericDrag(new HarrisPriester(CelestialBodyFactory.getSun(), earth),
98 new IsotropicDrag(2.5, 1.2),
99 provider.getMu());
100
101 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
102 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
103 4, 4, 4, 8, 4, 4, 2);
104 final DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
105
106 Orbit baseOrbit =
107 new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
108 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
109 provider.getMu());
110
111 double dt = 900;
112 double dP = 1.0;
113 for (OrbitType orbitType : orbitTypes) {
114 final Orbit initialOrbit = orbitType.convertType(baseOrbit);
115
116 DSSTPropagator propagator =
117 setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
118 propagator.setMu(provider.getMu());
119 for (final DSSTForceModel forceModel : propagator.getAllForceModels()) {
120 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
121 driver.setValue(driver.getReferenceValue());
122 driver.setSelected(driver.getName().equals(parameterName));
123 }
124 }
125
126 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, type);
127 final SpacecraftState initialState =
128 partials.setInitialJacobians(new SpacecraftState(initialOrbit));
129 propagator.setInitialState(initialState, PropagationType.MEAN);
130 final DSSTJacobiansMapper mapper = partials.getMapper();
131 PickUpHandler pickUp = new PickUpHandler(mapper, null);
132 propagator.setStepHandler(pickUp);
133 propagator.propagate(initialState.getDate().shiftedBy(dt));
134 double[][] dYdP = pickUp.getdYdP();
135
136
137 double[][] dYdPRef = new double[6][1];
138 DSSTPropagator propagator2 = setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
139 propagator2.setMu(provider.getMu());
140 ParameterDriversList bound = new ParameterDriversList();
141 for (final DSSTForceModel forceModel : propagator2.getAllForceModels()) {
142 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
143 if (driver.getName().equals(parameterName)) {
144 driver.setSelected(true);
145 bound.add(driver);
146 } else {
147 driver.setSelected(false);
148 }
149 }
150 }
151 ParameterDriver selected = bound.getDrivers().get(0);
152 double p0 = selected.getReferenceValue();
153 double h = selected.getScale();
154 selected.setValue(p0 - 4 * h);
155 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
156 orbitType,
157 initialState.getFrame(), initialState.getDate(),
158 propagator2.getMu(),
159 initialState.getAttitude()));
160 SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
161 selected.setValue(p0 - 3 * h);
162 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
163 orbitType,
164 initialState.getFrame(), initialState.getDate(),
165 propagator2.getMu(),
166 initialState.getAttitude()));
167 SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
168 selected.setValue(p0 - 2 * h);
169 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
170 orbitType,
171 initialState.getFrame(), initialState.getDate(),
172 propagator2.getMu(),
173 initialState.getAttitude()));
174 SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
175 selected.setValue(p0 - 1 * h);
176 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
177 orbitType,
178 initialState.getFrame(), initialState.getDate(),
179 propagator2.getMu(),
180 initialState.getAttitude()));
181 SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
182 selected.setValue(p0 + 1 * h);
183 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
184 orbitType,
185 initialState.getFrame(), initialState.getDate(),
186 propagator2.getMu(),
187 initialState.getAttitude()));
188 SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
189 selected.setValue(p0 + 2 * h);
190 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
191 orbitType,
192 initialState.getFrame(), initialState.getDate(),
193 propagator2.getMu(),
194 initialState.getAttitude()));
195 SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
196 selected.setValue(p0 + 3 * h);
197 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
198 orbitType,
199 initialState.getFrame(), initialState.getDate(),
200 propagator2.getMu(),
201 initialState.getAttitude()));
202 SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
203 selected.setValue(p0 + 4 * h);
204 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
205 orbitType,
206 initialState.getFrame(), initialState.getDate(),
207 propagator2.getMu(),
208 initialState.getAttitude()));
209 SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
210 fillJacobianColumn(dYdPRef, 0, orbitType, h,
211 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
212
213 for (int i = 0; i < 6; ++i) {
214 Assert.assertEquals(dYdPRef[i][0], dYdP[i][0], FastMath.abs(dYdPRef[i][0] * tolerance));
215 }
216
217 }
218
219 }
220
221 @Test
222 public void testPropagationTypesElliptical() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
223 doTestPropagation(PropagationType.MEAN, 7.0e-16);
224 }
225
226 @Test
227 public void testPropagationTypesEllipticalWithShortPeriod() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
228 doTestPropagation(PropagationType.OSCULATING, 3.3e-4);
229 }
230
231 private void doTestPropagation(PropagationType type,
232 double tolerance)
233 throws FileNotFoundException, UnsupportedEncodingException {
234
235 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
236
237 Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
238
239 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
240 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
241 4, 4, 4, 8, 4, 4, 2);
242
243 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
244 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
245 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
246 provider.getMu());
247
248 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
249
250 Orbit initialOrbit =
251 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
252 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
253 provider.getMu());
254 final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit);
255
256 double dt = 900;
257 double dP = 0.001;
258 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
259
260
261 DSSTPropagator propagator = setUpPropagator(type, orbit, dP, orbitType, srp, tesseral, zonal, moon);
262 propagator.setMu(provider.getMu());
263 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, type);
264 final SpacecraftState initialState =
265 partials.setInitialJacobians(new SpacecraftState(orbit));
266 final double[] stateVector = new double[6];
267 OrbitType.EQUINOCTIAL.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
268 final AbsoluteDate target = initialState.getDate().shiftedBy(dt);
269 propagator.setInitialState(initialState, type);
270 final DSSTJacobiansMapper mapper = partials.getMapper();
271 PickUpHandler pickUp = new PickUpHandler(mapper, null);
272 propagator.setStepHandler(pickUp);
273 propagator.propagate(target);
274 double[][] dYdY0 = pickUp.getdYdY0();
275
276
277 double[][] dYdY0Ref = new double[6][6];
278 DSSTPropagator propagator2 = setUpPropagator(type, orbit, dP, orbitType, srp, tesseral, zonal, moon);
279 propagator2.setMu(provider.getMu());
280 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
281 for (int i = 0; i < 6; ++i) {
282 propagator2.setInitialState(shiftState(initialState, orbitType, -4 * steps[i], i), type);
283 SpacecraftState sM4h = propagator2.propagate(target);
284 propagator2.setInitialState(shiftState(initialState, orbitType, -3 * steps[i], i), type);
285 SpacecraftState sM3h = propagator2.propagate(target);
286 propagator2.setInitialState(shiftState(initialState, orbitType, -2 * steps[i], i), type);
287 SpacecraftState sM2h = propagator2.propagate(target);
288 propagator2.setInitialState(shiftState(initialState, orbitType, -1 * steps[i], i), type);
289 SpacecraftState sM1h = propagator2.propagate(target);
290 propagator2.setInitialState(shiftState(initialState, orbitType, 1 * steps[i], i), type);
291 SpacecraftState sP1h = propagator2.propagate(target);
292 propagator2.setInitialState(shiftState(initialState, orbitType, 2 * steps[i], i), type);
293 SpacecraftState sP2h = propagator2.propagate(target);
294 propagator2.setInitialState(shiftState(initialState, orbitType, 3 * steps[i], i), type);
295 SpacecraftState sP3h = propagator2.propagate(target);
296 propagator2.setInitialState(shiftState(initialState, orbitType, 4 * steps[i], i), type);
297 SpacecraftState sP4h = propagator2.propagate(target);
298 fillJacobianColumn(dYdY0Ref, i, orbitType, steps[i],
299 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
300 }
301
302 for (int i = 0; i < 6; ++i) {
303 for (int j = 0; j < 6; ++j) {
304 if (stateVector[i] != 0) {
305 double error = FastMath.abs((dYdY0[i][j] - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
306 Assert.assertEquals(0, error, tolerance);
307 }
308 }
309 }
310 }
311
312 @Test(expected=OrekitException.class)
313 public void testNotInitialized() {
314 Orbit initialOrbit =
315 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
316 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
317 Constants.EIGEN5C_EARTH_MU);
318 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
319
320 double dP = 0.001;
321 DSSTPropagator propagator =
322 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
323 new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN).getMapper();
324 }
325
326 @Test(expected=OrekitException.class)
327 public void testTooSmallDimension() {
328 Orbit initialOrbit =
329 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
330 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
331 Constants.EIGEN5C_EARTH_MU);
332 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
333
334 double dP = 0.001;
335 DSSTPropagator propagator =
336 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
337 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
338 partials.setInitialJacobians(new SpacecraftState(orbit),
339 new double[5][6], new double[6][2]);
340 }
341
342 @Test(expected=OrekitException.class)
343 public void testTooLargeDimension() {
344 Orbit initialOrbit =
345 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
346 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
347 Constants.EIGEN5C_EARTH_MU);
348 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
349
350 double dP = 0.001;
351 DSSTPropagator propagator =
352 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
353 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
354 partials.setInitialJacobians(new SpacecraftState(orbit),
355 new double[8][6], new double[6][2]);
356 }
357
358 @Test(expected=OrekitException.class)
359 public void testMismatchedDimensions() {
360 Orbit initialOrbit =
361 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
362 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
363 Constants.EIGEN5C_EARTH_MU);
364 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
365
366 double dP = 0.001;
367 DSSTPropagator propagator =
368 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
369 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
370 partials.setInitialJacobians(new SpacecraftState(orbit),
371 new double[6][6], new double[7][2]);
372 }
373
374 @Test
375 public void testWrongParametersDimension() {
376 Orbit initialOrbit =
377 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
378 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
379 Constants.EIGEN5C_EARTH_MU);
380 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
381
382 double dP = 0.001;
383 DSSTForceModel sunAttraction = new DSSTThirdBody(CelestialBodyFactory.getSun(), Constants.EIGEN5C_EARTH_MU);
384 DSSTForceModel moonAttraction = new DSSTThirdBody(CelestialBodyFactory.getMoon(), Constants.EIGEN5C_EARTH_MU);
385 DSSTPropagator propagator =
386 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL,
387 sunAttraction, moonAttraction);
388 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
389 try {
390 partials.setInitialJacobians(new SpacecraftState(orbit),
391 new double[6][6], new double[6][3]);
392 partials.computeDerivatives(new SpacecraftState(orbit), new double[6]);
393 Assert.fail("an exception should have been thrown");
394 } catch (OrekitException oe) {
395 Assert.assertEquals(OrekitMessages.INITIAL_MATRIX_AND_PARAMETERS_NUMBER_MISMATCH,
396 oe.getSpecifier());
397 }
398 }
399
400 private void fillJacobianColumn(double[][] jacobian, int column,
401 OrbitType orbitType, double h,
402 SpacecraftState sM4h, SpacecraftState sM3h,
403 SpacecraftState sM2h, SpacecraftState sM1h,
404 SpacecraftState sP1h, SpacecraftState sP2h,
405 SpacecraftState sP3h, SpacecraftState sP4h) {
406 double[] aM4h = stateToArray(sM4h, orbitType)[0];
407 double[] aM3h = stateToArray(sM3h, orbitType)[0];
408 double[] aM2h = stateToArray(sM2h, orbitType)[0];
409 double[] aM1h = stateToArray(sM1h, orbitType)[0];
410 double[] aP1h = stateToArray(sP1h, orbitType)[0];
411 double[] aP2h = stateToArray(sP2h, orbitType)[0];
412 double[] aP3h = stateToArray(sP3h, orbitType)[0];
413 double[] aP4h = stateToArray(sP4h, orbitType)[0];
414 for (int i = 0; i < jacobian.length; ++i) {
415 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
416 32 * (aP3h[i] - aM3h[i]) -
417 168 * (aP2h[i] - aM2h[i]) +
418 672 * (aP1h[i] - aM1h[i])) / (840 * h);
419 }
420 }
421
422 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
423 double delta, int column) {
424
425 double[][] array = stateToArray(state, orbitType);
426 array[0][column] += delta;
427
428 return arrayToState(array, orbitType, state.getFrame(), state.getDate(),
429 state.getMu(), state.getAttitude());
430
431 }
432
433 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
434 double[][] array = new double[2][6];
435
436 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
437 return array;
438 }
439
440 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
441 Frame frame, AbsoluteDate date, double mu,
442 Attitude attitude) {
443 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngle.MEAN, date, mu, frame);
444 return new SpacecraftState(orbit, attitude);
445 }
446
447 private DSSTPropagator setUpPropagator(PropagationType type, Orbit orbit, double dP,
448 OrbitType orbitType,
449 DSSTForceModel... models) {
450
451 final double minStep = 6000.0;
452 final double maxStep = 86400.0;
453
454 double[][] tol = NumericalPropagator.tolerances(dP, orbit, orbitType);
455 DSSTPropagator propagator =
456 new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]), type);
457 for (DSSTForceModel model : models) {
458 propagator.addForceModel(model);
459 }
460 return propagator;
461 }
462
463 private static class PickUpHandler implements OrekitStepHandler {
464
465 private final DSSTJacobiansMapper mapper;
466 private final AbsoluteDate pickUpDate;
467 private final double[][] dYdY0;
468 private final double[][] dYdP;
469
470 public PickUpHandler(DSSTJacobiansMapper mapper, AbsoluteDate pickUpDate) {
471 this.mapper = mapper;
472 this.pickUpDate = pickUpDate;
473 dYdY0 = new double[DSSTJacobiansMapper.STATE_DIMENSION][DSSTJacobiansMapper.STATE_DIMENSION];
474 dYdP = new double[DSSTJacobiansMapper.STATE_DIMENSION][mapper.getParameters()];
475 }
476
477 public double[][] getdYdY0() {
478 return dYdY0;
479 }
480
481 public double[][] getdYdP() {
482 return dYdP;
483 }
484
485 public void handleStep(OrekitStepInterpolator interpolator) {
486 if (pickUpDate != null) {
487
488 double dt0 = pickUpDate.durationFrom(interpolator.getPreviousState().getDate());
489 double dt1 = pickUpDate.durationFrom(interpolator.getCurrentState().getDate());
490 if (dt0 * dt1 > 0) {
491
492 return;
493 } else {
494 checkState(interpolator.getInterpolatedState(pickUpDate));
495 }
496 }
497 }
498
499 public void finish(SpacecraftState finalState) {
500 checkState(finalState);
501 }
502
503 private void checkState(final SpacecraftState state) {
504 Assert.assertEquals(1, state.getAdditionalStatesValues().size());
505 Assert.assertEquals(mapper.getName(), state.getAdditionalStatesValues().getData().get(0).getKey());
506 mapper.setShortPeriodJacobians(state);
507 mapper.getStateJacobian(state, dYdY0);
508 mapper.getParametersJacobian(state, dYdP);
509
510 }
511
512 }
513
514
515
516
517
518 @Test
519 public void testIssue713() {
520 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
521
522 Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
523
524 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
525 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
526 4, 4, 4, 8, 4, 4, 2);
527
528 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
529 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
530 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
531 provider.getMu());
532
533 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
534
535 Orbit initialOrbit =
536 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
537 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
538 provider.getMu());
539 final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit);
540
541 double dP = 0.001;
542 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
543
544
545 DSSTPropagator propagatorMEAN = setUpPropagator(PropagationType.MEAN, orbit, dP, orbitType, srp, tesseral, zonal, moon);
546 propagatorMEAN.setMu(provider.getMu());
547 DSSTPartialDerivativesEquations partialsMEAN = new DSSTPartialDerivativesEquations("partials", propagatorMEAN, PropagationType.MEAN);
548 final SpacecraftState initialStateMEAN =
549 partialsMEAN.setInitialJacobians(new SpacecraftState(orbit));
550 final DSSTJacobiansMapper mapperMEAN = partialsMEAN.getMapper();
551 mapperMEAN.setShortPeriodJacobians(initialStateMEAN);
552 double[][] dYdY0MEAN = new double[DSSTJacobiansMapper.STATE_DIMENSION][DSSTJacobiansMapper.STATE_DIMENSION];
553 mapperMEAN.getStateJacobian(initialStateMEAN, dYdY0MEAN);
554 for (int i = 0; i < 6; ++i) {
555 for (int j = 0; j < 6; ++j) {
556 if (i == j) {
557 Assert.assertEquals(1.0, dYdY0MEAN[i][j], 1e-9);
558 }
559 else {
560 Assert.assertEquals(0.0, dYdY0MEAN[i][j], 1e-9);
561 }
562 }
563 }
564
565
566 DSSTPropagator propagatorOSC = setUpPropagator(PropagationType.OSCULATING, orbit, dP, orbitType, srp, tesseral, zonal, moon);
567 propagatorOSC.setMu(provider.getMu());
568 DSSTPartialDerivativesEquations partialsOSC = new DSSTPartialDerivativesEquations("partials", propagatorOSC, PropagationType.OSCULATING);
569 final SpacecraftState initialStateOSC =
570 partialsOSC.setInitialJacobians(new SpacecraftState(orbit));
571 final DSSTJacobiansMapper mapperOSC = partialsOSC.getMapper();
572 mapperOSC.setShortPeriodJacobians(initialStateOSC);
573 double[][] dYdY0OSC = new double[DSSTJacobiansMapper.STATE_DIMENSION][DSSTJacobiansMapper.STATE_DIMENSION];
574 mapperOSC.getStateJacobian(initialStateOSC, dYdY0OSC);
575 final double[] refLine1 = new double[] {1.0000, -5750.3478, 15270.6488, -2707.1208, -2165.0148, -178.3653};
576 final double[] refLine6 = new double[] {0.0000, 0.0035, 0.0013, -0.0005, 0.0005, 1.0000};
577 for (int i = 0; i < 6; ++i) {
578 Assert.assertEquals(refLine1[i], dYdY0OSC[0][i], 1e-4);
579 Assert.assertEquals(refLine6[i], dYdY0OSC[5][i], 1e-4);
580 }
581
582 }
583
584 @Before
585 public void setUp() {
586 Utils.setDataRoot("regular-data:potential/shm-format");
587 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
588 }
589
590 }
591