1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.conversion;
18
19 import static org.junit.Assert.assertFalse;
20 import static org.junit.Assert.assertTrue;
21
22 import java.io.IOException;
23 import java.text.ParseException;
24 import java.util.List;
25
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28 import org.hipparchus.util.FastMath;
29 import org.junit.Assert;
30 import org.junit.Before;
31 import org.junit.Test;
32 import org.orekit.Utils;
33 import org.orekit.bodies.OneAxisEllipsoid;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.errors.OrekitMessages;
36 import org.orekit.forces.ForceModel;
37 import org.orekit.forces.drag.DragForce;
38 import org.orekit.forces.drag.DragSensitive;
39 import org.orekit.forces.drag.IsotropicDrag;
40 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
41 import org.orekit.forces.gravity.NewtonianAttraction;
42 import org.orekit.forces.gravity.potential.GravityFieldFactory;
43 import org.orekit.frames.FramesFactory;
44 import org.orekit.models.earth.atmosphere.Atmosphere;
45 import org.orekit.models.earth.atmosphere.SimpleExponentialAtmosphere;
46 import org.orekit.orbits.EquinoctialOrbit;
47 import org.orekit.orbits.Orbit;
48 import org.orekit.orbits.OrbitType;
49 import org.orekit.orbits.PositionAngle;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.propagation.integration.AdditionalDerivativesProvider;
52 import org.orekit.propagation.numerical.NumericalPropagator;
53 import org.orekit.time.AbsoluteDate;
54 import org.orekit.time.TimeScalesFactory;
55 import org.orekit.utils.Constants;
56 import org.orekit.utils.IERSConventions;
57 import org.orekit.utils.PVCoordinates;
58 import org.orekit.utils.ParameterDriver;
59
60 public class NumericalConverterTest {
61
62 private double mu;
63 private double minStep;
64 private double maxStep;
65 private double dP;
66
67 private Orbit orbit;
68 private NumericalPropagator propagator;
69 private ForceModel gravity;
70 private ForceModel drag;
71 private Atmosphere atmosphere;
72 private double crossSection;
73
74 @Test
75 public void testIssue598() {
76
77 final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
78
79 final NumericalPropagatorBuilder builder =
80 new NumericalPropagatorBuilder(OrbitType.CIRCULAR.convertType(orbit),
81 dp54Builder,
82 PositionAngle.TRUE, 1.0);
83 builder.addForceModel(gravity);
84
85 assertFalse(hasNewtonianAttraction(builder.getAllForceModels()));
86
87 builder.buildPropagator(builder.getSelectedNormalizedParameters());
88
89 assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
90
91 builder.addForceModel(drag);
92 assertTrue(hasNewtonianAttraction(builder.getAllForceModels()));
93 }
94
95 @Test
96 public void testOnlyCartesianAllowed() {
97 NumericalPropagatorBuilder builder =
98 new NumericalPropagatorBuilder(OrbitType.CIRCULAR.convertType(orbit),
99 new LutherIntegratorBuilder(100.0),
100 PositionAngle.TRUE, 1.0);
101 builder.addForceModel(drag);
102 builder.addForceModel(gravity);
103 try {
104 new JacobianPropagatorConverter(builder, 1.0e-3, 5000);
105 Assert.fail("an exception should have been thrown");
106 } catch (OrekitException oe) {
107 Assert.assertEquals(OrekitMessages.ORBIT_TYPE_NOT_ALLOWED, oe.getSpecifier());
108 Assert.assertEquals(OrbitType.CIRCULAR, oe.getParts()[0]);
109 Assert.assertEquals(OrbitType.CARTESIAN, oe.getParts()[1]);
110 }
111 }
112
113 @Test
114 public void testConversionWithoutParameters() throws IOException, ParseException {
115 checkFit(orbit, 6000, 300, 1.0e-3, 0.855);
116 }
117
118 @Test
119 public void testConversionWithFreeParameter() throws IOException, ParseException {
120 checkFit(orbit, 6000, 300, 1.0e-3, 0.826,
121 DragSensitive.DRAG_COEFFICIENT, NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT);
122 }
123
124 @Test
125 public void testIntegrators01() {
126
127 ODEIntegratorBuilder abBuilder = new AdamsBashforthIntegratorBuilder(2, minStep, maxStep, dP);
128 checkFit(abBuilder);
129 }
130
131 @Test
132 public void testIntegrators02() {
133
134 ODEIntegratorBuilder amBuilder = new AdamsMoultonIntegratorBuilder(2, minStep, maxStep, dP);
135 checkFit(amBuilder);
136 }
137
138 @Test
139 public void testIntegrators03() {
140
141 final double stepSize = 100.;
142
143 ODEIntegratorBuilder crkBuilder = new ClassicalRungeKuttaIntegratorBuilder(stepSize);
144 checkFit(crkBuilder);
145 }
146
147 @Test
148 public void testIntegrators04() {
149
150 final double stepSize = 100.;
151
152 ODEIntegratorBuilder lBuilder = new LutherIntegratorBuilder(stepSize);
153 checkFit(lBuilder);
154 }
155
156 @Test
157 public void testIntegrators05() {
158
159 ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
160 checkFit(dp54Builder);
161 }
162
163 @Test
164 public void testIntegrators06() {
165
166 final double stepSize = 100.;
167
168 ODEIntegratorBuilder eBuilder = new EulerIntegratorBuilder(stepSize);
169 checkFit(eBuilder);
170 }
171
172 @Test
173 public void testIntegrators07() {
174
175 final double stepSize = 100.;
176
177 ODEIntegratorBuilder gBuilder = new GillIntegratorBuilder(stepSize);
178 checkFit(gBuilder);
179 }
180
181 @Test
182 public void testIntegrators08() {
183
184 ODEIntegratorBuilder gbsBuilder = new GraggBulirschStoerIntegratorBuilder(minStep, maxStep, dP);
185 checkFit(gbsBuilder);
186 }
187
188 @Test
189 public void testIntegrators09() {
190
191 ODEIntegratorBuilder hh54Builder = new HighamHall54IntegratorBuilder(minStep, maxStep, dP);
192 checkFit(hh54Builder);
193 }
194
195 @Test
196 public void testIntegrators10() {
197
198 final double stepSize = 100.;
199
200 ODEIntegratorBuilder mBuilder = new MidpointIntegratorBuilder(stepSize);
201 checkFit(mBuilder);
202 }
203
204 @Test
205 public void testIntegrators11() {
206
207 final double stepSize = 100.;
208
209 ODEIntegratorBuilder teBuilder = new ThreeEighthesIntegratorBuilder(stepSize);
210 checkFit(teBuilder);
211 }
212
213 @Test
214 public void testAdditionalEquations() {
215
216 final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
217
218 final NumericalPropagatorBuilder builder =
219 new NumericalPropagatorBuilder(OrbitType.CIRCULAR.convertType(orbit),
220 dp54Builder,
221 PositionAngle.TRUE, 1.0);
222 builder.addForceModel(drag);
223 builder.addForceModel(gravity);
224
225
226 builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
227
228 public String getName() {
229 return "linear";
230 }
231
232 public int getDimension() {
233 return 1;
234 }
235
236 public double[] derivatives(SpacecraftState s) {
237 return new double[] { 1.0 };
238 }
239
240 });
241
242 builder.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
243
244 public String getName() {
245 return "linear";
246 }
247
248 public int getDimension() {
249 return 1;
250 }
251
252 public double[] derivatives(SpacecraftState s) {
253 return new double[] { 1.0 };
254 }
255
256 });
257
258 try {
259
260 builder.buildPropagator(builder.getSelectedNormalizedParameters());
261 Assert.fail("an exception should have been thrown");
262 } catch (OrekitException oe) {
263 Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
264 }
265 }
266
267 @Test
268 public void testDeselectOrbitals() {
269
270 final ODEIntegratorBuilder dp54Builder = new DormandPrince54IntegratorBuilder(minStep, maxStep, dP);
271
272 final NumericalPropagatorBuilder builder =
273 new NumericalPropagatorBuilder(OrbitType.CIRCULAR.convertType(orbit),
274 dp54Builder,
275 PositionAngle.TRUE, 1.0);
276 for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
277 Assert.assertTrue(driver.isSelected());
278 }
279 builder.deselectDynamicParameters();
280 for (ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
281 Assert.assertFalse(driver.isSelected());
282 }
283 }
284
285 protected void checkFit(final Orbit orbit, final double duration,
286 final double stepSize, final double threshold,
287 final double expectedRMS,
288 final String... freeParameters)
289 throws IOException, ParseException {
290
291 NumericalPropagatorBuilder builder =
292 new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit),
293 new DormandPrince853IntegratorBuilder(minStep, maxStep, dP),
294 PositionAngle.TRUE, dP);
295
296 ForceModel guessedDrag = drag;
297 ForceModel guessedGravity = gravity;
298 for (String param: freeParameters) {
299 if (DragSensitive.DRAG_COEFFICIENT.equals(param)) {
300
301 ParameterDriver driver = drag.getParameterDriver(param);
302 double coeff = driver.getReferenceValue() - driver.getScale();
303 guessedDrag = new DragForce(atmosphere, new IsotropicDrag(crossSection, coeff));
304 } else if (NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT.equals(param)) {
305
306 guessedGravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
307 GravityFieldFactory.getNormalizedProvider(2, 0));
308 ParameterDriver driver = guessedGravity.getParameterDriver(param);
309 driver.setValue(driver.getReferenceValue() + driver.getScale());
310 }
311 }
312 builder.addForceModel(guessedDrag);
313 builder.addForceModel(guessedGravity);
314
315 JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder,
316 threshold,
317 5000);
318
319 fitter.convert(propagator, duration, 1 + (int) (duration / stepSize), freeParameters);
320
321 NumericalPropagator prop = (NumericalPropagator)fitter.getAdaptedPropagator();
322 Orbit fitted = prop.getInitialState().getOrbit();
323
324 for (String param: freeParameters) {
325 for (ForceModel force: propagator.getAllForceModels()) {
326 if (force.isSupported(param)) {
327 for (ForceModel model: prop.getAllForceModels()) {
328 if (model.isSupported(param)) {
329 Assert.assertEquals(force.getParameterDriver(param).getValue(),
330 model.getParameterDriver(param).getValue(),
331 3.0e-4 * FastMath.abs(force.getParameterDriver(param).getValue()));
332 }
333 }
334 }
335 }
336 }
337
338 Assert.assertEquals(expectedRMS, fitter.getRMS(), 0.01 * expectedRMS);
339
340 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getX(),
341 fitted.getPVCoordinates().getPosition().getX(),
342 1.1);
343 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getY(),
344 fitted.getPVCoordinates().getPosition().getY(),
345 1.1);
346 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getZ(),
347 fitted.getPVCoordinates().getPosition().getZ(),
348 1.1);
349
350 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getX(),
351 fitted.getPVCoordinates().getVelocity().getX(),
352 0.0005);
353 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getY(),
354 fitted.getPVCoordinates().getVelocity().getY(),
355 0.0005);
356 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getZ(),
357 fitted.getPVCoordinates().getVelocity().getZ(),
358 0.0005);
359 }
360
361 protected void checkFit(final ODEIntegratorBuilder foiBuilder) {
362
363 NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit),
364 foiBuilder,
365 PositionAngle.TRUE,
366 1.0);
367
368 builder.addForceModel(drag);
369 builder.addForceModel(gravity);
370 builder.setAttitudeProvider(Utils.defaultLaw());
371 builder.setMass(1000.0);
372
373 JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, 1.0, 500);
374
375 fitter.convert(propagator, 1000., 11);
376
377 NumericalPropagator prop = (NumericalPropagator)fitter.getAdaptedPropagator();
378 Orbit fitted = prop.getInitialState().getOrbit();
379
380 final double peps = 1.e-1;
381 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getX(),
382 fitted.getPVCoordinates().getPosition().getX(),
383 peps * FastMath.abs(orbit.getPVCoordinates().getPosition().getX()));
384 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getY(),
385 fitted.getPVCoordinates().getPosition().getY(),
386 peps * FastMath.abs(orbit.getPVCoordinates().getPosition().getY()));
387 Assert.assertEquals(orbit.getPVCoordinates().getPosition().getZ(),
388 fitted.getPVCoordinates().getPosition().getZ(),
389 peps * FastMath.abs(orbit.getPVCoordinates().getPosition().getZ()));
390
391 final double veps = 5.e-1;
392 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getX(),
393 fitted.getPVCoordinates().getVelocity().getX(),
394 veps * FastMath.abs(orbit.getPVCoordinates().getVelocity().getX()));
395 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getY(),
396 fitted.getPVCoordinates().getVelocity().getY(),
397 veps * FastMath.abs(orbit.getPVCoordinates().getVelocity().getY()));
398 Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getZ(),
399 fitted.getPVCoordinates().getVelocity().getZ(),
400 veps * FastMath.abs(orbit.getPVCoordinates().getVelocity().getZ()));
401 }
402
403 @Before
404 public void setUp() throws IOException, ParseException {
405
406 Utils.setDataRoot("regular-data:potential/shm-format");
407 gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
408 GravityFieldFactory.getNormalizedProvider(2, 0));
409 mu = gravity.getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).getValue();
410 minStep = 1.0;
411 maxStep = 600.0;
412 dP = 10.0;
413
414
415 final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6).normalize()
416 .scalarMultiply(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 300e3);
417 final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
418 final AbsoluteDate initDate = new AbsoluteDate(2010, 10, 10, 10, 10, 10.0, TimeScalesFactory.getUTC());
419 orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
420 FramesFactory.getEME2000(), initDate, mu);
421
422 final double[][] tol = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN);
423 propagator = new NumericalPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]));
424 propagator.setInitialState(new SpacecraftState(orbit));
425 propagator.setOrbitType(OrbitType.CARTESIAN);
426
427
428 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
429 Constants.WGS84_EARTH_FLATTENING,
430 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
431 earth.setAngularThreshold(1.e-7);
432 atmosphere = new SimpleExponentialAtmosphere(earth, 0.0004, 42000.0, 7500.0);
433 final double dragCoef = 2.0;
434 crossSection = 10.0;
435 drag = new DragForce(atmosphere, new IsotropicDrag(crossSection, dragCoef));
436
437 propagator.addForceModel(gravity);
438 propagator.addForceModel(drag);
439 }
440
441 private boolean hasNewtonianAttraction(final List<ForceModel> forceModels) {
442 final int last = forceModels.size() - 1;
443 return last >= 0 && forceModels.get(last) instanceof NewtonianAttraction;
444 }
445
446 }
447