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