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 org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
20 import org.hipparchus.util.FastMath;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.BeforeEach;
23 import org.junit.jupiter.api.Test;
24 import org.orekit.Utils;
25 import org.orekit.attitudes.Attitude;
26 import org.orekit.bodies.CelestialBodyFactory;
27 import org.orekit.bodies.OneAxisEllipsoid;
28 import org.orekit.forces.drag.DragSensitive;
29 import org.orekit.forces.drag.IsotropicDrag;
30 import org.orekit.forces.gravity.potential.GravityFieldFactory;
31 import org.orekit.forces.gravity.potential.SHMFormatReader;
32 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
33 import org.orekit.frames.Frame;
34 import org.orekit.frames.FramesFactory;
35 import org.orekit.models.earth.atmosphere.HarrisPriester;
36 import org.orekit.orbits.EquinoctialOrbit;
37 import org.orekit.orbits.KeplerianOrbit;
38 import org.orekit.orbits.Orbit;
39 import org.orekit.orbits.OrbitType;
40 import org.orekit.orbits.PositionAngleType;
41 import org.orekit.propagation.PropagationType;
42 import org.orekit.propagation.SpacecraftState;
43 import org.orekit.propagation.ToleranceProvider;
44 import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
45 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
46 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
47 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
48 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
49 import org.orekit.time.AbsoluteDate;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.IERSConventions;
52 import org.orekit.utils.ParameterDriver;
53 import org.orekit.utils.ParameterDriversList;
54
55 import java.io.IOException;
56 import java.text.ParseException;
57
58
59 public class DSSTIntegrableJacobianColumnGeneratorTest {
60
61 @BeforeEach
62 public void setUp() {
63 Utils.setDataRoot("regular-data:potential/shm-format");
64 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
65 }
66
67 @Test
68 public void testDragParametersDerivatives() throws ParseException, IOException {
69 doTestParametersDerivatives(DragSensitive.DRAG_COEFFICIENT,
70 2.4e-3,
71 PropagationType.MEAN,
72 OrbitType.EQUINOCTIAL);
73 }
74
75 @Test
76 public void testMuParametersDerivatives() throws ParseException, IOException {
77 doTestParametersDerivatives(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
78 5.e-3,
79 PropagationType.MEAN,
80 OrbitType.EQUINOCTIAL);
81 }
82
83 private void doTestParametersDerivatives(String parameterName, double tolerance,
84 PropagationType type,
85 OrbitType... orbitTypes) {
86
87 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
88 Constants.WGS84_EARTH_FLATTENING,
89 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
90
91 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
92
93 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
94
95 DSSTForceModel drag = new DSSTAtmosphericDrag(new HarrisPriester(CelestialBodyFactory.getSun(), earth),
96 new IsotropicDrag(2.5, 1.2),
97 provider.getMu());
98
99 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
100 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
101 4, 4, 4, 8, 4, 4, 2);
102 final DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
103
104 Orbit baseOrbit =
105 new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngleType.MEAN,
106 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
107 provider.getMu());
108
109 double dt = 900;
110 double dP = 1.0;
111 for (OrbitType orbitType : orbitTypes) {
112 final Orbit initialOrbit = orbitType.convertType(baseOrbit);
113
114 DSSTPropagator propagator =
115 setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
116 propagator.setMu(provider.getMu());
117 for (final DSSTForceModel forceModel : propagator.getAllForceModels()) {
118 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
119 driver.setValue(driver.getReferenceValue());
120 driver.setSelected(driver.getName().equals(parameterName));
121 }
122 }
123
124 final SpacecraftState initialState = new SpacecraftState(initialOrbit);
125 propagator.setInitialState(initialState, PropagationType.MEAN);
126 PickUpHandler pickUp = new PickUpHandler(propagator, null, null, DragSensitive.DRAG_COEFFICIENT);
127 propagator.setStepHandler(pickUp);
128 propagator.propagate(initialState.getDate().shiftedBy(dt));
129
130
131 double[][] dYdPRef = new double[6][1];
132 DSSTPropagator propagator2 = setUpPropagator(type, initialOrbit, dP, orbitType, zonal, tesseral, drag);
133 propagator2.setMu(provider.getMu());
134 ParameterDriversList bound = new ParameterDriversList();
135 for (final DSSTForceModel forceModel : propagator2.getAllForceModels()) {
136 for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
137 if (driver.getName().equals(parameterName)) {
138 driver.setSelected(true);
139 bound.add(driver);
140 } else {
141 driver.setSelected(false);
142 }
143 }
144 }
145 ParameterDriver selected = bound.getDrivers().get(0);
146 double p0 = selected.getReferenceValue();
147 double h = selected.getScale();
148 selected.setValue(p0 - 4 * h);
149 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
150 orbitType,
151 initialState.getFrame(), initialState.getDate(),
152 propagator2.getMu(),
153 initialState.getAttitude()));
154 SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
155 selected.setValue(p0 - 3 * h);
156 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
157 orbitType,
158 initialState.getFrame(), initialState.getDate(),
159 propagator2.getMu(),
160 initialState.getAttitude()));
161 SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
162 selected.setValue(p0 - 2 * h);
163 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
164 orbitType,
165 initialState.getFrame(), initialState.getDate(),
166 propagator2.getMu(),
167 initialState.getAttitude()));
168 SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
169 selected.setValue(p0 - 1 * h);
170 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
171 orbitType,
172 initialState.getFrame(), initialState.getDate(),
173 propagator2.getMu(),
174 initialState.getAttitude()));
175 SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
176 selected.setValue(p0 + 1 * h);
177 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
178 orbitType,
179 initialState.getFrame(), initialState.getDate(),
180 propagator2.getMu(),
181 initialState.getAttitude()));
182 SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
183 selected.setValue(p0 + 2 * h);
184 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
185 orbitType,
186 initialState.getFrame(), initialState.getDate(),
187 propagator2.getMu(),
188 initialState.getAttitude()));
189 SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
190 selected.setValue(p0 + 3 * h);
191 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
192 orbitType,
193 initialState.getFrame(), initialState.getDate(),
194 propagator2.getMu(),
195 initialState.getAttitude()));
196 SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
197 selected.setValue(p0 + 4 * h);
198 propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType),
199 orbitType,
200 initialState.getFrame(), initialState.getDate(),
201 propagator2.getMu(),
202 initialState.getAttitude()));
203 SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
204 fillJacobianColumn(dYdPRef, 0, orbitType, h,
205 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
206
207 for (int i = 0; i < 6; ++i) {
208 Assertions.assertEquals(dYdPRef[i][0], pickUp.getdYdP().getEntry(i, 0), FastMath.abs(dYdPRef[i][0] * tolerance));
209 }
210
211 }
212
213 }
214
215 private void fillJacobianColumn(double[][] jacobian, int column,
216 OrbitType orbitType, double h,
217 SpacecraftState sM4h, SpacecraftState sM3h,
218 SpacecraftState sM2h, SpacecraftState sM1h,
219 SpacecraftState sP1h, SpacecraftState sP2h,
220 SpacecraftState sP3h, SpacecraftState sP4h) {
221 double[] aM4h = stateToArray(sM4h, orbitType)[0];
222 double[] aM3h = stateToArray(sM3h, orbitType)[0];
223 double[] aM2h = stateToArray(sM2h, orbitType)[0];
224 double[] aM1h = stateToArray(sM1h, orbitType)[0];
225 double[] aP1h = stateToArray(sP1h, orbitType)[0];
226 double[] aP2h = stateToArray(sP2h, orbitType)[0];
227 double[] aP3h = stateToArray(sP3h, orbitType)[0];
228 double[] aP4h = stateToArray(sP4h, orbitType)[0];
229 for (int i = 0; i < jacobian.length; ++i) {
230 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
231 32 * (aP3h[i] - aM3h[i]) -
232 168 * (aP2h[i] - aM2h[i]) +
233 672 * (aP1h[i] - aM1h[i])) / (840 * h);
234 }
235 }
236
237 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
238 double[][] array = new double[2][6];
239
240 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
241 return array;
242 }
243
244 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
245 Frame frame, AbsoluteDate date, double mu,
246 Attitude attitude) {
247 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngleType.MEAN, date, mu, frame);
248 return new SpacecraftState(orbit, attitude);
249 }
250
251 private DSSTPropagator setUpPropagator(PropagationType type, Orbit orbit, double dP,
252 OrbitType orbitType,
253 DSSTForceModel... models) {
254
255 final double minStep = 6000.0;
256 final double maxStep = 86400.0;
257
258 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, orbitType);
259 DSSTPropagator propagator =
260 new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]), type);
261 for (DSSTForceModel model : models) {
262 propagator.addForceModel(model);
263 }
264 return propagator;
265 }
266
267 }