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 static org.hamcrest.CoreMatchers.is;
20
21 import java.io.FileNotFoundException;
22 import java.io.IOException;
23 import java.io.UnsupportedEncodingException;
24 import java.text.ParseException;
25 import java.util.ArrayList;
26 import java.util.Collections;
27 import java.util.List;
28
29 import org.hamcrest.MatcherAssert;
30 import org.hipparchus.CalculusFieldElement;
31 import org.hipparchus.Field;
32 import org.hipparchus.analysis.differentiation.Gradient;
33 import org.hipparchus.linear.RealMatrix;
34 import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
35 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
36 import org.hipparchus.util.FastMath;
37 import org.hipparchus.util.MathArrays;
38 import org.junit.Assert;
39 import org.junit.Before;
40 import org.junit.Test;
41 import org.orekit.Utils;
42 import org.orekit.attitudes.Attitude;
43 import org.orekit.attitudes.AttitudeProvider;
44 import org.orekit.bodies.CelestialBodyFactory;
45 import org.orekit.bodies.OneAxisEllipsoid;
46 import org.orekit.errors.OrekitException;
47 import org.orekit.errors.OrekitMessages;
48 import org.orekit.forces.drag.DragSensitive;
49 import org.orekit.forces.drag.IsotropicDrag;
50 import org.orekit.forces.gravity.potential.GravityFieldFactory;
51 import org.orekit.forces.gravity.potential.SHMFormatReader;
52 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
53 import org.orekit.frames.Frame;
54 import org.orekit.frames.FramesFactory;
55 import org.orekit.models.earth.atmosphere.HarrisPriester;
56 import org.orekit.orbits.EquinoctialOrbit;
57 import org.orekit.orbits.KeplerianOrbit;
58 import org.orekit.orbits.Orbit;
59 import org.orekit.orbits.OrbitType;
60 import org.orekit.orbits.PositionAngle;
61 import org.orekit.propagation.FieldSpacecraftState;
62 import org.orekit.propagation.PropagationType;
63 import org.orekit.propagation.SpacecraftState;
64 import org.orekit.propagation.events.EventDetector;
65 import org.orekit.propagation.events.FieldEventDetector;
66 import org.orekit.propagation.numerical.NumericalPropagator;
67 import org.orekit.propagation.sampling.OrekitStepHandler;
68 import org.orekit.propagation.sampling.OrekitStepInterpolator;
69 import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
70 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
71 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
72 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
73 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
74 import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
75 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
76 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
77 import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
78 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
79 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
80 import org.orekit.time.AbsoluteDate;
81 import org.orekit.utils.Constants;
82 import org.orekit.utils.IERSConventions;
83 import org.orekit.utils.ParameterDriver;
84 import org.orekit.utils.ParameterDriversList;
85
86
87 @Deprecated
88 public class DSSTPartialDerivativesEquationsTest {
89
90
91 private static final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
92
93 private static final double gm = Constants.EIGEN5C_EARTH_MU;
94
95 private static final Frame eci = FramesFactory.getGCRF();
96
97
98 private DSSTPropagator propagator;
99
100 private MockForceModel forceModel;
101
102 private SpacecraftState state;
103
104 private DSSTPartialDerivativesEquations pde;
105
106
107
108
109
110 @Before
111 public void setUp() {
112 Utils.setDataRoot("regular-data:potential/shm-format");
113 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
114 propagator = new DSSTPropagator(new DormandPrince54Integrator(1, 500, 0.001, 0.001));
115 forceModel = new MockForceModel();
116 propagator.addForceModel(forceModel);
117 pde = new DSSTPartialDerivativesEquations("pde", propagator, PropagationType.MEAN);
118 final Orbit orbit = new EquinoctialOrbit(4.2163393E7,
119 -0.25925449177598586,
120 -0.06946703170551687,
121 0.15995912655021305,
122 -0.5969755874197339,
123 15.47576793123677,
124 PositionAngle.MEAN,
125 eci, date, gm);
126 state = new SpacecraftState(orbit).addAdditionalState("pde", new double[2 * 3 * 6]);
127 pde.setInitialJacobians(state);
128
129 }
130
131 @Test
132 public void testDragParametersDerivatives() throws ParseException, IOException {
133 doTestParametersDerivatives(DragSensitive.DRAG_COEFFICIENT,
134 2.4e-3,
135 PropagationType.MEAN,
136 OrbitType.EQUINOCTIAL);
137 }
138
139 @Test
140 public void testMuParametersDerivatives() throws ParseException, IOException {
141 doTestParametersDerivatives(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
142 5.e-3,
143 PropagationType.MEAN,
144 OrbitType.EQUINOCTIAL);
145 }
146
147 private void doTestParametersDerivatives(String parameterName, double tolerance,
148 PropagationType type,
149 OrbitType... orbitTypes) {
150
151 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
152 Constants.WGS84_EARTH_FLATTENING,
153 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
154
155 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
156
157 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
158
159 DSSTForceModel drag = new DSSTAtmosphericDrag(new HarrisPriester(CelestialBodyFactory.getSun(), earth),
160 new IsotropicDrag(2.5, 1.2),
161 provider.getMu());
162
163 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
164 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
165 4, 4, 4, 8, 4, 4, 2);
166 final DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
167
168 Orbit baseOrbit =
169 new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
170 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
171 provider.getMu());
172
173 double dt = 900;
174 double dP = 1.0;
175 for (OrbitType orbitType : orbitTypes) {
176 final Orbit initialOrbit = orbitType.convertType(baseOrbit);
177
178 DSSTPropagator propagator =
179 setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
180 propagator.setMu(provider.getMu());
181 for (final DSSTForceModel forceModel : propagator.getAllForceModels()) {
182 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
183 driver.setValue(driver.getReferenceValue());
184 driver.setSelected(driver.getName().equals(parameterName));
185 }
186 }
187
188 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, type);
189 final SpacecraftState initialState =
190 partials.setInitialJacobians(new SpacecraftState(initialOrbit));
191 propagator.setInitialState(initialState, PropagationType.MEAN);
192 final DSSTJacobiansMapper mapper = partials.getMapper();
193 PickUpHandler pickUp = new PickUpHandler(mapper, null);
194 propagator.setStepHandler(pickUp);
195 propagator.propagate(initialState.getDate().shiftedBy(dt));
196
197
198 double[][] dYdPRef = new double[6][1];
199 DSSTPropagator propagator2 = setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
200 propagator2.setMu(provider.getMu());
201 ParameterDriversList bound = new ParameterDriversList();
202 for (final DSSTForceModel forceModel : propagator2.getAllForceModels()) {
203 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
204 if (driver.getName().equals(parameterName)) {
205 driver.setSelected(true);
206 bound.add(driver);
207 } else {
208 driver.setSelected(false);
209 }
210 }
211 }
212 ParameterDriver selected = bound.getDrivers().get(0);
213 double p0 = selected.getReferenceValue();
214 double h = selected.getScale();
215 selected.setValue(p0 - 4 * h);
216 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
217 orbitType,
218 initialState.getFrame(), initialState.getDate(),
219 propagator2.getMu(),
220 initialState.getAttitude()));
221 SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
222 selected.setValue(p0 - 3 * h);
223 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
224 orbitType,
225 initialState.getFrame(), initialState.getDate(),
226 propagator2.getMu(),
227 initialState.getAttitude()));
228 SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
229 selected.setValue(p0 - 2 * h);
230 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
231 orbitType,
232 initialState.getFrame(), initialState.getDate(),
233 propagator2.getMu(),
234 initialState.getAttitude()));
235 SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
236 selected.setValue(p0 - 1 * h);
237 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
238 orbitType,
239 initialState.getFrame(), initialState.getDate(),
240 propagator2.getMu(),
241 initialState.getAttitude()));
242 SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
243 selected.setValue(p0 + 1 * h);
244 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
245 orbitType,
246 initialState.getFrame(), initialState.getDate(),
247 propagator2.getMu(),
248 initialState.getAttitude()));
249 SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
250 selected.setValue(p0 + 2 * h);
251 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
252 orbitType,
253 initialState.getFrame(), initialState.getDate(),
254 propagator2.getMu(),
255 initialState.getAttitude()));
256 SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
257 selected.setValue(p0 + 3 * h);
258 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
259 orbitType,
260 initialState.getFrame(), initialState.getDate(),
261 propagator2.getMu(),
262 initialState.getAttitude()));
263 SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
264 selected.setValue(p0 + 4 * h);
265 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
266 orbitType,
267 initialState.getFrame(), initialState.getDate(),
268 propagator2.getMu(),
269 initialState.getAttitude()));
270 SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
271 fillJacobianColumn(dYdPRef, 0, orbitType, h,
272 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
273
274 for (int i = 0; i < 6; ++i) {
275 Assert.assertEquals(dYdPRef[i][0], pickUp.dYdP.getEntry(i, 0), FastMath.abs(dYdPRef[i][0] * tolerance));
276 }
277
278 }
279
280 }
281
282 @Test
283 public void testPropagationTypesElliptical() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
284 doTestPropagation(PropagationType.MEAN, 7.0e-16);
285 }
286
287 @Test
288 public void testPropagationTypesEllipticalWithShortPeriod() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
289 doTestPropagation(PropagationType.OSCULATING, 3.3e-4);
290 }
291
292 private void doTestPropagation(PropagationType type,
293 double tolerance)
294 throws FileNotFoundException, UnsupportedEncodingException {
295
296 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
297
298 Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
299
300 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
301 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
302 4, 4, 4, 8, 4, 4, 2);
303
304 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
305 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
306 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
307 provider.getMu());
308
309 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
310
311 Orbit initialOrbit =
312 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
313 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
314 provider.getMu());
315 final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit);
316
317 double dt = 900;
318 double dP = 0.001;
319 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
320
321
322 DSSTPropagator propagator = setUpPropagator(type, orbit, dP, orbitType, srp, tesseral, zonal, moon);
323 propagator.setMu(provider.getMu());
324 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, type);
325 final SpacecraftState initialState =
326 partials.setInitialJacobians(new SpacecraftState(orbit));
327 final double[] stateVector = new double[6];
328 OrbitType.EQUINOCTIAL.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
329 final AbsoluteDate target = initialState.getDate().shiftedBy(dt);
330 propagator.setInitialState(initialState, type);
331 final DSSTJacobiansMapper mapper = partials.getMapper();
332 PickUpHandler pickUp = new PickUpHandler(mapper, null);
333 propagator.setStepHandler(pickUp);
334 propagator.propagate(target);
335
336
337 double[][] dYdY0Ref = new double[6][6];
338 DSSTPropagator propagator2 = setUpPropagator(type, orbit, dP, orbitType, srp, tesseral, zonal, moon);
339 propagator2.setMu(provider.getMu());
340 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
341 for (int i = 0; i < 6; ++i) {
342 propagator2.setInitialState(shiftState(initialState, orbitType, -4 * steps[i], i), type);
343 SpacecraftState sM4h = propagator2.propagate(target);
344 propagator2.setInitialState(shiftState(initialState, orbitType, -3 * steps[i], i), type);
345 SpacecraftState sM3h = propagator2.propagate(target);
346 propagator2.setInitialState(shiftState(initialState, orbitType, -2 * steps[i], i), type);
347 SpacecraftState sM2h = propagator2.propagate(target);
348 propagator2.setInitialState(shiftState(initialState, orbitType, -1 * steps[i], i), type);
349 SpacecraftState sM1h = propagator2.propagate(target);
350 propagator2.setInitialState(shiftState(initialState, orbitType, 1 * steps[i], i), type);
351 SpacecraftState sP1h = propagator2.propagate(target);
352 propagator2.setInitialState(shiftState(initialState, orbitType, 2 * steps[i], i), type);
353 SpacecraftState sP2h = propagator2.propagate(target);
354 propagator2.setInitialState(shiftState(initialState, orbitType, 3 * steps[i], i), type);
355 SpacecraftState sP3h = propagator2.propagate(target);
356 propagator2.setInitialState(shiftState(initialState, orbitType, 4 * steps[i], i), type);
357 SpacecraftState sP4h = propagator2.propagate(target);
358 fillJacobianColumn(dYdY0Ref, i, orbitType, steps[i],
359 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
360 }
361
362 for (int i = 0; i < 6; ++i) {
363 for (int j = 0; j < 6; ++j) {
364 if (stateVector[i] != 0) {
365 double error = FastMath.abs((pickUp.dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
366 Assert.assertEquals(0, error, tolerance);
367 }
368 }
369 }
370 }
371
372 @Test(expected=OrekitException.class)
373 public void testNotInitialized() {
374 Orbit initialOrbit =
375 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
376 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
377 Constants.EIGEN5C_EARTH_MU);
378 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
379
380 double dP = 0.001;
381 DSSTPropagator propagator =
382 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
383 new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN).getMapper();
384 }
385
386 @Test(expected=OrekitException.class)
387 public void testTooSmallDimension() {
388 Orbit initialOrbit =
389 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
390 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
391 Constants.EIGEN5C_EARTH_MU);
392 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
393
394 double dP = 0.001;
395 DSSTPropagator propagator =
396 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
397 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
398 partials.setInitialJacobians(new SpacecraftState(orbit),
399 new double[5][6], new double[6][2]);
400 }
401
402 @Test(expected=OrekitException.class)
403 public void testTooLargeDimension() {
404 Orbit initialOrbit =
405 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
406 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
407 Constants.EIGEN5C_EARTH_MU);
408 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
409
410 double dP = 0.001;
411 DSSTPropagator propagator =
412 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
413 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
414 partials.setInitialJacobians(new SpacecraftState(orbit),
415 new double[8][6], new double[6][2]);
416 }
417
418 @Test(expected=OrekitException.class)
419 public void testMismatchedDimensions() {
420 Orbit initialOrbit =
421 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
422 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
423 Constants.EIGEN5C_EARTH_MU);
424 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
425
426 double dP = 0.001;
427 DSSTPropagator propagator =
428 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL);
429 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
430 partials.setInitialJacobians(new SpacecraftState(orbit),
431 new double[6][6], new double[7][2]);
432 }
433
434 @Test
435 public void testWrongParametersDimension() {
436 Orbit initialOrbit =
437 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
438 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
439 Constants.EIGEN5C_EARTH_MU);
440 final Orbit orbit = OrbitType.EQUINOCTIAL.convertType(initialOrbit);
441
442 double dP = 0.001;
443 DSSTForceModel sunAttraction = new DSSTThirdBody(CelestialBodyFactory.getSun(), Constants.EIGEN5C_EARTH_MU);
444 DSSTForceModel moonAttraction = new DSSTThirdBody(CelestialBodyFactory.getMoon(), Constants.EIGEN5C_EARTH_MU);
445 DSSTPropagator propagator =
446 setUpPropagator(PropagationType.MEAN, orbit, dP, OrbitType.EQUINOCTIAL,
447 sunAttraction, moonAttraction);
448 DSSTPartialDerivativesEquations partials = new DSSTPartialDerivativesEquations("partials", propagator, PropagationType.MEAN);
449 try {
450 partials.setInitialJacobians(new SpacecraftState(orbit),
451 new double[6][6], new double[6][3]);
452 partials.derivatives(new SpacecraftState(orbit));
453 Assert.fail("an exception should have been thrown");
454 } catch (OrekitException oe) {
455 Assert.assertEquals(OrekitMessages.INITIAL_MATRIX_AND_PARAMETERS_NUMBER_MISMATCH,
456 oe.getSpecifier());
457 }
458 }
459
460 private void fillJacobianColumn(double[][] jacobian, int column,
461 OrbitType orbitType, double h,
462 SpacecraftState sM4h, SpacecraftState sM3h,
463 SpacecraftState sM2h, SpacecraftState sM1h,
464 SpacecraftState sP1h, SpacecraftState sP2h,
465 SpacecraftState sP3h, SpacecraftState sP4h) {
466 double[] aM4h = stateToArray(sM4h, orbitType)[0];
467 double[] aM3h = stateToArray(sM3h, orbitType)[0];
468 double[] aM2h = stateToArray(sM2h, orbitType)[0];
469 double[] aM1h = stateToArray(sM1h, orbitType)[0];
470 double[] aP1h = stateToArray(sP1h, orbitType)[0];
471 double[] aP2h = stateToArray(sP2h, orbitType)[0];
472 double[] aP3h = stateToArray(sP3h, orbitType)[0];
473 double[] aP4h = stateToArray(sP4h, orbitType)[0];
474 for (int i = 0; i < jacobian.length; ++i) {
475 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
476 32 * (aP3h[i] - aM3h[i]) -
477 168 * (aP2h[i] - aM2h[i]) +
478 672 * (aP1h[i] - aM1h[i])) / (840 * h);
479 }
480 }
481
482 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
483 double delta, int column) {
484
485 double[][] array = stateToArray(state, orbitType);
486 array[0][column] += delta;
487
488 return arrayToState(array, orbitType, state.getFrame(), state.getDate(),
489 state.getMu(), state.getAttitude());
490
491 }
492
493 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
494 double[][] array = new double[2][6];
495
496 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
497 return array;
498 }
499
500 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
501 Frame frame, AbsoluteDate date, double mu,
502 Attitude attitude) {
503 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngle.MEAN, date, mu, frame);
504 return new SpacecraftState(orbit, attitude);
505 }
506
507 private DSSTPropagator setUpPropagator(PropagationType type, Orbit orbit, double dP,
508 OrbitType orbitType,
509 DSSTForceModel... models) {
510
511 final double minStep = 6000.0;
512 final double maxStep = 86400.0;
513
514 double[][] tol = NumericalPropagator.tolerances(dP, orbit, orbitType);
515 DSSTPropagator propagator =
516 new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]), type);
517 for (DSSTForceModel model : models) {
518 propagator.addForceModel(model);
519 }
520 return propagator;
521 }
522
523 private static class PickUpHandler implements OrekitStepHandler {
524
525 private final DSSTJacobiansMapper mapper;
526 private final AbsoluteDate pickUpDate;
527 private RealMatrix dYdY0;
528 private RealMatrix dYdP;
529
530 public PickUpHandler(DSSTJacobiansMapper mapper, AbsoluteDate pickUpDate) {
531 this.mapper = mapper;
532 this.pickUpDate = pickUpDate;
533 }
534
535 public void handleStep(OrekitStepInterpolator interpolator) {
536 if (pickUpDate != null) {
537
538 double dt0 = pickUpDate.durationFrom(interpolator.getPreviousState().getDate());
539 double dt1 = pickUpDate.durationFrom(interpolator.getCurrentState().getDate());
540 if (dt0 * dt1 > 0) {
541
542 return;
543 } else {
544 checkState(interpolator.getInterpolatedState(pickUpDate));
545 }
546 }
547 }
548
549 public void finish(SpacecraftState finalState) {
550 checkState(finalState);
551 }
552
553 private void checkState(final SpacecraftState state) {
554 Assert.assertEquals(1, state.getAdditionalStates().size());
555 Assert.assertTrue(state.getAdditionalStates().containsKey(mapper.getName()));
556 mapper.setReferenceState(state);
557 dYdY0 = mapper.getStateTransitionMatrix(state);
558 dYdP = mapper.getParametersJacobian(state);
559
560 }
561
562 }
563
564
565
566
567
568 @Test
569 public void testIssue713() {
570 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
571
572 Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
573
574 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
575 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
576 4, 4, 4, 8, 4, 4, 2);
577
578 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
579 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
580 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
581 provider.getMu());
582
583 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
584
585 Orbit initialOrbit =
586 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
587 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
588 provider.getMu());
589 final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit);
590
591 double dP = 0.001;
592 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
593
594
595 DSSTPropagator propagatorMEAN = setUpPropagator(PropagationType.MEAN, orbit, dP, orbitType, srp, tesseral, zonal, moon);
596 propagatorMEAN.setMu(provider.getMu());
597 DSSTPartialDerivativesEquations partialsMEAN = new DSSTPartialDerivativesEquations("partials", propagatorMEAN, PropagationType.MEAN);
598 final SpacecraftState initialStateMEAN =
599 partialsMEAN.setInitialJacobians(new SpacecraftState(orbit));
600 final DSSTJacobiansMapper mapperMEAN = partialsMEAN.getMapper();
601 mapperMEAN.setReferenceState(initialStateMEAN);
602 RealMatrix dYdY0MEAN = mapperMEAN.getStateTransitionMatrix(initialStateMEAN);
603 for (int i = 0; i < 6; ++i) {
604 for (int j = 0; j < 6; ++j) {
605 Assert.assertEquals(i == j ? 1.0 : 0.0, dYdY0MEAN.getEntry(i, j), 1e-9);
606 }
607 }
608
609
610 DSSTPropagator propagatorOSC = setUpPropagator(PropagationType.OSCULATING, orbit, dP, orbitType, srp, tesseral, zonal, moon);
611 propagatorOSC.setMu(provider.getMu());
612 DSSTPartialDerivativesEquations partialsOSC = new DSSTPartialDerivativesEquations("partials", propagatorOSC, PropagationType.OSCULATING);
613 final SpacecraftState initialStateOSC =
614 partialsOSC.setInitialJacobians(new SpacecraftState(orbit));
615 final DSSTJacobiansMapper mapperOSC = partialsOSC.getMapper();
616 mapperOSC.setReferenceState(initialStateOSC);
617 RealMatrix dYdY0OSC = mapperOSC.getStateTransitionMatrix(initialStateOSC);
618 final double[] refLine1 = new double[] {1.0000, -5750.3478, 15270.6488, -2707.1208, -2165.0148, -178.3653};
619 final double[] refLine6 = new double[] {0.0000, 0.0035, 0.0013, -0.0005, 0.0005, 1.0000};
620 for (int i = 0; i < 6; ++i) {
621 Assert.assertEquals(refLine1[i], dYdY0OSC.getEntry(0, i), 1e-4);
622 Assert.assertEquals(refLine6[i], dYdY0OSC.getEntry(5, i), 1e-4);
623 }
624
625 }
626
627
628
629
630
631 @Test
632 public void testDerivatives() {
633
634
635 pde.derivatives(state);
636
637
638 MatcherAssert.assertThat(forceModel.sma.getReal(), is(state.getA()));
639 MatcherAssert.assertThat(forceModel.ex.getReal(), is(state.getEquinoctialEx()));
640 MatcherAssert.assertThat(forceModel.ey.getReal(), is(state.getEquinoctialEy()));
641 MatcherAssert.assertThat(forceModel.hx.getReal(), is(state.getHx()));
642 MatcherAssert.assertThat(forceModel.hy.getReal(), is(state.getHy()));
643 MatcherAssert.assertThat(forceModel.l.getReal(), is(state.getLv()));
644
645 }
646
647
648 private static class MockForceModel implements DSSTForceModel {
649
650
651 public Gradient sma;
652
653
654 public Gradient ex;
655
656
657 public Gradient ey;
658
659
660 public Gradient hx;
661
662
663 public Gradient hy;
664
665
666 public Gradient l;
667
668 @Override
669 public List<ShortPeriodTerms> initializeShortPeriodTerms(AuxiliaryElements auxiliaryElements,
670 PropagationType type,
671 double[] parameters) {
672 return new ArrayList<ShortPeriodTerms>();
673 }
674
675 @Override
676 public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(FieldAuxiliaryElements<T> auxiliaryElements,
677 PropagationType type,
678 T[] parameters) {
679 return new ArrayList<FieldShortPeriodTerms<T>>();
680 }
681
682 @Override
683 public double[] getMeanElementRate(SpacecraftState state,
684 AuxiliaryElements auxiliaryElements,
685 double[] parameters) {
686 return new double[] {state.getA(),
687 state.getEquinoctialEx(),
688 state.getEquinoctialEy(),
689 state.getHx(),
690 state.getHy(),
691 state.getLv()};
692 }
693
694 @Override
695 public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(FieldSpacecraftState<T> state,
696 FieldAuxiliaryElements<T> auxiliaryElements,
697 T[] parameters) {
698
699 final Field<T> field = state.getDate().getField();
700
701 this.sma = (Gradient) state.getA();
702 this.ex = (Gradient) state.getEquinoctialEx();
703 this.ey = (Gradient) state.getEquinoctialEy();
704 this.hx = (Gradient) state.getHx();
705 this.hy = (Gradient) state.getHy();
706 this.l = (Gradient) state.getLv();
707
708 final T[] elements = MathArrays.buildArray(field, 6);
709 elements[0] = state.getA();
710 elements[1] = state.getEquinoctialEx();
711 elements[2] = state.getEquinoctialEy();
712 elements[3] = state.getHx();
713 elements[4] = state.getHy();
714 elements[5] = state.getLv();
715
716 return elements;
717
718 }
719
720 @Override
721 public EventDetector[] getEventsDetectors() {
722 return new EventDetector[0];
723 }
724
725 @Override
726 public <T extends CalculusFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(Field<T> field) {
727 return null;
728 }
729
730 @Override
731 public void registerAttitudeProvider(AttitudeProvider provider) {
732 }
733
734 @Override
735 public void updateShortPeriodTerms(double[] parameters, SpacecraftState... meanStates) {
736 }
737
738 @Override
739 @SuppressWarnings("unchecked")
740 public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(T[] parameters,
741 FieldSpacecraftState<T>... meanStates) {
742 }
743
744 @Override
745 public List<ParameterDriver> getParametersDrivers() {
746 return Collections.emptyList();
747 }
748
749 }
750 }