1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.integration;
18
19 import org.hipparchus.geometry.euclidean.threed.Rotation;
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.MatrixUtils;
22 import org.hipparchus.linear.RealMatrix;
23 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.Binary64Field;
26 import org.junit.jupiter.api.Assertions;
27 import org.junit.jupiter.api.BeforeEach;
28 import org.junit.jupiter.api.Test;
29 import org.orekit.TestUtils;
30 import org.orekit.Utils;
31 import org.orekit.attitudes.*;
32 import org.orekit.bodies.CelestialBody;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.forces.gravity.potential.GravityFieldFactory;
36 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.orbits.EquinoctialOrbit;
39 import org.orekit.orbits.Orbit;
40 import org.orekit.orbits.OrbitType;
41 import org.orekit.propagation.*;
42 import org.orekit.propagation.analytical.KeplerianPropagator;
43 import org.orekit.propagation.events.DateDetector;
44 import org.orekit.propagation.events.EventDetector;
45 import org.orekit.propagation.numerical.NumericalPropagator;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.utils.AbsolutePVCoordinatesTest;
48 import org.orekit.utils.Constants;
49 import org.orekit.utils.PVCoordinates;
50
51 import java.util.stream.Stream;
52
53 class IntegratedEphemerisTest {
54
55 @Test
56 void testNormalKeplerIntegration() {
57
58
59 KeplerianPropagator keplerEx = new KeplerianPropagator(initialOrbit);
60
61
62
63
64 AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
65 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
66 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
67 numericalPropagator.propagate(finalDate);
68 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
69 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
70
71
72 for (int i = 1; i <= Constants.JULIAN_DAY; i++) {
73 AbsoluteDate intermediateDate = initialOrbit.getDate().shiftedBy(i);
74 SpacecraftState keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
75 SpacecraftState numericIntermediateOrbit = ephemeris.propagate(intermediateDate);
76 Vector3D kepPosition = keplerIntermediateOrbit.getPosition();
77 Vector3D numPosition = numericIntermediateOrbit.getPosition();
78 Assertions.assertEquals(0, kepPosition.subtract(numPosition).getNorm(), 0.06);
79 }
80
81
82 AbsoluteDate intermediateDate = initialOrbit.getDate().shiftedBy(41589);
83 SpacecraftState keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
84 SpacecraftState state = keplerEx.propagate(finalDate);
85 numericalPropagator.setInitialState(state);
86 final EphemerisGenerator generator2 = numericalPropagator.getEphemerisGenerator();
87 numericalPropagator.propagate(initialOrbit.getDate());
88 BoundedPropagator invEphemeris = generator2.getGeneratedEphemeris();
89 SpacecraftState numericIntermediateOrbit = invEphemeris.propagate(intermediateDate);
90 Vector3D kepPosition = keplerIntermediateOrbit.getPosition();
91 Vector3D numPosition = numericIntermediateOrbit.getPosition();
92 Assertions.assertEquals(0, kepPosition.subtract(numPosition).getNorm(), 10e-2);
93
94 }
95
96 @Test
97 void testGetAttitudeProvider() {
98
99 final EventDetector detector = new DateDetector();
100 final AttitudeProvider attitudeProvider = new FrameAlignedProvider(initialOrbit.getFrame()) {
101 @Override
102 public Stream<EventDetector> getEventDetectors() {
103 return Stream.of(detector);
104 }
105 };
106 numericalPropagator.setAttitudeProvider(attitudeProvider);
107 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
108 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
109 numericalPropagator.propagate(initialOrbit.getDate().shiftedBy(10));
110
111 final BoundedPropagator boundedPropagator = generator.getGeneratedEphemeris();
112
113 final AttitudeProvider ephemerisAttitudeProvider = boundedPropagator.getAttitudeProvider();
114 Assertions.assertEquals(0, ephemerisAttitudeProvider.getEventDetectors().count());
115 Assertions.assertEquals(0, ephemerisAttitudeProvider.getFieldEventDetectors(Binary64Field.getInstance()).count());
116 Assertions.assertInstanceOf(AttitudeProviderModifier.class, ephemerisAttitudeProvider);
117 final AttitudeProviderModifier providerModifier = (AttitudeProviderModifier) ephemerisAttitudeProvider;
118 Assertions.assertEquals(attitudeProvider, providerModifier.getUnderlyingAttitudeProvider());
119 final Attitude expectedAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
120 final Attitude actualAttitude = ephemerisAttitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
121 Assertions.assertEquals(expectedAttitude.getSpin(), actualAttitude.getSpin());
122 Assertions.assertEquals(expectedAttitude.getRotationAcceleration(), actualAttitude.getRotationAcceleration());
123 Assertions.assertEquals(0., Rotation.distance(expectedAttitude.getRotation(), actualAttitude.getRotation()));
124 }
125
126 @Test
127 void testPartialDerivativesIssue16() {
128
129 final String eqName = "derivatives";
130 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
131 numericalPropagator.setOrbitType(OrbitType.CARTESIAN);
132 final MatricesHarvester harvester = numericalPropagator.setupMatricesComputation(eqName, null, null);
133 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
134 numericalPropagator.propagate(initialOrbit.getDate().shiftedBy(3600.0));
135 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
136 ephemeris.setStepHandler(interpolator -> {
137 SpacecraftState state = interpolator.getCurrentState();
138 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
139 harvester.getParametersJacobian(state);
140 RealMatrix deltaId = dYdY0.subtract(MatrixUtils.createRealIdentityMatrix(6));
141 Assertions.assertTrue(deltaId.getNorm1() > 100);
142 Assertions.assertTrue(deltaId.getNorm1() < 3100);
143 });
144
145 ephemeris.propagate(initialOrbit.getDate().shiftedBy(1800.0));
146
147 }
148
149 @Test
150 void testGetFrame() {
151
152 AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
153 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
154 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
155 numericalPropagator.propagate(finalDate);
156 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
157 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
158
159
160 Assertions.assertNotNull(ephemeris.getFrame());
161 Assertions.assertSame(ephemeris.getFrame(), numericalPropagator.getFrame());
162 }
163
164 @Test
165 void testIssue766() {
166
167
168 AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
169 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
170 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
171 numericalPropagator.propagate(finalDate);
172 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
173 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
174
175
176 Assertions.assertInstanceOf(AttitudeProviderModifier.class, ephemeris.getAttitudeProvider());
177
178
179 CelestialBody sun = CelestialBodyFactory.getSun();
180 ephemeris.setAttitudeProvider(new CelestialBodyPointed(FramesFactory.getEME2000(), sun, Vector3D.PLUS_K,
181 Vector3D.PLUS_I, Vector3D.PLUS_K));
182 Assertions.assertInstanceOf(CelestialBodyPointed.class, ephemeris.getAttitudeProvider());
183
184 }
185
186 @Test
187 void testAdditionalDerivatives() {
188
189 AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(10.0);
190 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(1e-3).getTolerances(initialOrbit, OrbitType.CARTESIAN);
191 DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-6, 10.0, tolerances[0], tolerances[1]);
192 integrator.setInitialStepSize(1.0e-3);
193 NumericalPropagator propagator = new NumericalPropagator(integrator);
194 final DerivativesProvider provider1 = new DerivativesProvider("provider-1", 3);
195 propagator.addAdditionalDerivativesProvider(provider1);
196 final DerivativesProvider provider2 = new DerivativesProvider("provider-2", 1);
197 propagator.addAdditionalDerivativesProvider(provider2);
198 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
199 propagator.setInitialState(new SpacecraftState(initialOrbit).
200 addAdditionalData(provider1.getName(), new double[provider1.getDimension()]).
201 addAdditionalData(provider2.getName(), new double[provider2.getDimension()]));
202 propagator.propagate(finalDate);
203 IntegratedEphemeris ephemeris = (IntegratedEphemeris) generator.getGeneratedEphemeris();
204
205 for (double dt = 0; dt < ephemeris.getMaxDate().durationFrom(ephemeris.getMinDate()); dt += 0.1) {
206 SpacecraftState state = ephemeris.propagate(ephemeris.getMinDate().shiftedBy(dt));
207 checkState(dt, state, provider1);
208 checkState(dt, state, provider2);
209 }
210
211
212 Assertions.assertEquals(SpacecraftState.DEFAULT_MASS, ephemeris.getMass(initialOrbit.getDate()));
213 assertOrbit(initialOrbit, ephemeris.propagateOrbit(initialOrbit.getDate()), 1e-8);
214
215
216 Exception thrown = Assertions.assertThrows(OrekitException.class,
217 () -> ephemeris.resetIntermediateState(new SpacecraftState(initialOrbit), true));
218 Assertions.assertEquals(
219 "reset state not allowed", thrown.getMessage());
220
221 }
222
223
224 @Test
225 void testIssue949() {
226
227 final AbsoluteDate initialDate = new AbsoluteDate();
228 numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
229 numericalPropagator.setOrbitType(OrbitType.CARTESIAN);
230
231
232 final AdditionalDataProvider<double[]> additionalDataProvider = TestUtils.getAdditionalProviderWithInit();
233 numericalPropagator.addAdditionalDataProvider(additionalDataProvider);
234
235
236 final EphemerisGenerator generator = numericalPropagator.getEphemerisGenerator();
237 numericalPropagator.propagate(initialDate.shiftedBy(1));
238
239 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
240
241
242 Assertions.assertDoesNotThrow(() -> ephemeris.propagate(ephemeris.getMaxDate()), "No error should have been thrown");
243
244 }
245
246 public static void assertOrbit(final Orbit expected, final Orbit actual, final double epsilon) {
247 Assertions.assertEquals(expected.getMu(), actual.getMu());
248 Assertions.assertEquals(expected.getType(), actual.getType());
249 AbsolutePVCoordinatesTest.assertPV(expected.getPVCoordinates(expected.getFrame()),
250 actual.getPVCoordinates(expected.getFrame()),
251 epsilon);
252 }
253
254 private void checkState(final double dt, final SpacecraftState state, final DerivativesProvider provider) {
255
256 Assertions.assertTrue(state.hasAdditionalData(provider.getName()));
257 Assertions.assertEquals(provider.getDimension(), state.getAdditionalState(provider.getName()).length);
258 for (int i = 0; i < provider.getDimension(); ++i) {
259 Assertions.assertEquals(i * dt, state.getAdditionalState(provider.getName())[i], 4.0e-15 * i * dt);
260 }
261
262 Assertions.assertTrue(state.hasAdditionalStateDerivative(provider.getName()));
263 Assertions.assertEquals(provider.getDimension(), state.getAdditionalStateDerivative(provider.getName()).length);
264 for (int i = 0; i < provider.getDimension(); ++i) {
265 Assertions.assertEquals(i, state.getAdditionalStateDerivative(provider.getName())[i], 2.0e-14 * i);
266 }
267
268 }
269
270 @BeforeEach
271 public void setUp() {
272
273 Utils.setDataRoot("regular-data:potential/icgem-format");
274 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
275
276
277 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
278 Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
279 double mu = 3.9860047e14;
280
281 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
282 initialOrbit =
283 new EquinoctialOrbit(new PVCoordinates(position, velocity),
284 FramesFactory.getEME2000(), initDate, mu);
285
286
287 double[] absTolerance = {
288 0.0001, 1.0e-11, 1.0e-11, 1.0e-8, 1.0e-8, 1.0e-8, 0.001
289 };
290 double[] relTolerance = {
291 1.0e-8, 1.0e-8, 1.0e-8, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-7
292 };
293 AdaptiveStepsizeIntegrator integrator =
294 new DormandPrince853Integrator(0.001, 500, absTolerance, relTolerance);
295 integrator.setInitialStepSize(100);
296 numericalPropagator = new NumericalPropagator(integrator);
297
298 }
299
300 private Orbit initialOrbit;
301 private NumericalPropagator numericalPropagator;
302
303 private static class DerivativesProvider implements AdditionalDerivativesProvider {
304 private final String name;
305 private final double[] derivatives;
306 DerivativesProvider(final String name, final int dimension) {
307 this.name = name;
308 this.derivatives = new double[dimension];
309 for (int i = 0; i < dimension; ++i) {
310 derivatives[i] = i;
311 }
312 }
313 public String getName() {
314 return name;
315 }
316 public int getDimension() {
317 return derivatives.length;
318 }
319 public CombinedDerivatives combinedDerivatives(final SpacecraftState s) {
320 return new CombinedDerivatives(derivatives, null);
321 }
322 }
323
324 }