1 package org.orekit.propagation.analytical;
2
3 import org.hipparchus.geometry.euclidean.threed.Vector3D;
4 import org.hipparchus.linear.RealMatrix;
5 import org.junit.Assert;
6 import org.junit.Before;
7 import org.junit.Test;
8 import org.orekit.Utils;
9 import org.orekit.forces.gravity.potential.GravityFieldFactory;
10 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
11 import org.orekit.frames.FramesFactory;
12 import org.orekit.orbits.CartesianOrbit;
13 import org.orekit.orbits.Orbit;
14 import org.orekit.orbits.OrbitType;
15 import org.orekit.orbits.PositionAngle;
16 import org.orekit.propagation.SpacecraftState;
17 import org.orekit.time.AbsoluteDate;
18 import org.orekit.utils.PVCoordinates;
19 import org.orekit.utils.ParameterDriver;
20 import org.orekit.utils.ParameterDriversList;
21
22 public class BrouwerLyddaneParametersDerivativesTest {
23
24
25 private UnnormalizedSphericalHarmonicsProvider provider;
26
27 @Before
28 public void setUp() {
29 Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
30 provider = GravityFieldFactory.getUnnormalizedProvider(5, 0);
31 }
32
33 @Test
34 public void testNoEstimatedParameters() {
35
36
37
38
39 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
40 Vector3D position = new Vector3D(3220103., 69623., 6149822.);
41 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
42 Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
43 FramesFactory.getEME2000(), initDate, provider.getMu());
44
45 final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(initialOrbit, provider, 1.0e-14);
46 final SpacecraftState initialState = propagator.getInitialState();
47 final double[] stateVector = new double[6];
48 OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
49 BrouwerLyddaneHarvester harvester = (BrouwerLyddaneHarvester) propagator.setupMatricesComputation("stm", null, null);
50 harvester.freezeColumnsNames();
51 RealMatrix dYdP = harvester.getParametersJacobian(initialState);
52 Assert.assertNull(dYdP);
53 }
54
55 @Test
56 public void testParametersDerivatives() {
57
58
59
60
61
62 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
63 Vector3D position = new Vector3D(3220103., 69623., 6149822.);
64 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
65 Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
66 FramesFactory.getEME2000(), initDate, provider.getMu());
67
68
69 final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(initialOrbit, provider, 1.0e-14);
70
71
72 ParameterDriversList bound = new ParameterDriversList();
73 final ParameterDriver M2 = propagator.getParametersDrivers().get(0);
74 M2.setSelected(true);
75 bound.add(M2);
76
77
78 final SpacecraftState initialState = propagator.getInitialState();
79 final double[] stateVector = new double[6];
80 OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
81 final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getKeplerianPeriod());
82 BrouwerLyddaneHarvester harvester = (BrouwerLyddaneHarvester) propagator.setupMatricesComputation("stm", null, null);
83 harvester.freezeColumnsNames();
84 final SpacecraftState finalState = propagator.propagate(target);
85 RealMatrix dYdP = harvester.getParametersJacobian(finalState);
86
87
88 OrbitType orbitType = OrbitType.CARTESIAN;
89 BrouwerLyddanePropagator propagator2;
90 double[][] dYdPRef = new double[6][1];
91
92 ParameterDriver selected = bound.getDrivers().get(0);
93 double p0 = selected.getReferenceValue();
94 double h = selected.getScale();
95 selected.setValue(p0 - 4 * h);
96 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
97 SpacecraftState sM4h = propagator2.propagate(target);
98 selected.setValue(p0 - 3 * h);
99 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
100 SpacecraftState sM3h = propagator2.propagate(target);
101 selected.setValue(p0 - 2 * h);
102 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
103 SpacecraftState sM2h = propagator2.propagate(target);
104 selected.setValue(p0 - 1 * h);
105 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
106 SpacecraftState sM1h = propagator2.propagate(target);
107 selected.setValue(p0 + 1 * h);
108 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
109 SpacecraftState sP1h = propagator2.propagate(target);
110 selected.setValue(p0 + 2 * h);
111 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
112 SpacecraftState sP2h = propagator2.propagate(target);
113 selected.setValue(p0 + 3 * h);
114 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
115 SpacecraftState sP3h = propagator2.propagate(target);
116 selected.setValue(p0 + 4 * h);
117 propagator2 = new BrouwerLyddanePropagator(initialOrbit, provider, selected.getValue());
118 SpacecraftState sP4h = propagator2.propagate(target);
119 fillJacobianColumn(dYdPRef, 0, orbitType, h,
120 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
121
122 for (int i = 0; i < 6; ++i) {
123 Assert.assertEquals(0.0, (dYdPRef[i][0] - dYdP.getEntry(i, 0)) / dYdPRef[i][0], 8.49e-13);
124 }
125
126 }
127
128 private void fillJacobianColumn(double[][] jacobian, int column,
129 OrbitType orbitType, double h,
130 SpacecraftState sM4h, SpacecraftState sM3h,
131 SpacecraftState sM2h, SpacecraftState sM1h,
132 SpacecraftState sP1h, SpacecraftState sP2h,
133 SpacecraftState sP3h, SpacecraftState sP4h) {
134 double[] aM4h = stateToArray(sM4h, orbitType)[0];
135 double[] aM3h = stateToArray(sM3h, orbitType)[0];
136 double[] aM2h = stateToArray(sM2h, orbitType)[0];
137 double[] aM1h = stateToArray(sM1h, orbitType)[0];
138 double[] aP1h = stateToArray(sP1h, orbitType)[0];
139 double[] aP2h = stateToArray(sP2h, orbitType)[0];
140 double[] aP3h = stateToArray(sP3h, orbitType)[0];
141 double[] aP4h = stateToArray(sP4h, orbitType)[0];
142 for (int i = 0; i < jacobian.length; ++i) {
143 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
144 32 * (aP3h[i] - aM3h[i]) -
145 168 * (aP2h[i] - aM2h[i]) +
146 672 * (aP1h[i] - aM1h[i])) / (840 * h);
147 }
148 }
149
150 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
151 double[][] array = new double[2][6];
152
153 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
154 return array;
155 }
156
157 }