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