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 java.lang.reflect.InvocationTargetException;
20 import java.lang.reflect.Method;
21 import java.util.List;
22
23 import org.hipparchus.CalculusFieldElement;
24 import org.hipparchus.Field;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.geometry.euclidean.threed.Rotation;
27 import org.hipparchus.ode.events.Action;
28 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
29 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
30 import org.hipparchus.util.Binary64;
31 import org.hipparchus.util.Binary64Field;
32 import org.hipparchus.util.MathArrays;
33 import org.junit.jupiter.api.Assertions;
34 import org.junit.jupiter.api.BeforeEach;
35 import org.junit.jupiter.api.Test;
36 import org.orekit.TestUtils;
37 import org.orekit.Utils;
38 import org.orekit.attitudes.Attitude;
39 import org.orekit.attitudes.AttitudeProvider;
40 import org.orekit.attitudes.AttitudeProviderModifier;
41 import org.orekit.attitudes.FrameAlignedProvider;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.errors.OrekitMessages;
44 import org.orekit.forces.gravity.potential.GravityFieldFactory;
45 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
46 import org.orekit.frames.FramesFactory;
47 import org.orekit.orbits.FieldEquinoctialOrbit;
48 import org.orekit.orbits.FieldOrbit;
49 import org.orekit.orbits.Orbit;
50 import org.orekit.orbits.OrbitType;
51 import org.orekit.propagation.*;
52 import org.orekit.propagation.analytical.FieldKeplerianPropagator;
53 import org.orekit.propagation.events.FieldDateDetector;
54 import org.orekit.propagation.numerical.FieldNumericalPropagator;
55 import org.orekit.time.FieldAbsoluteDate;
56 import org.orekit.utils.Constants;
57 import org.orekit.utils.FieldPVCoordinates;
58
59
60 public class FieldIntegratedEphemerisTest {
61
62 @Test
63 public void testNormalKeplerIntegration() {
64 doTestNormalKeplerIntegration(Binary64Field.getInstance());
65 }
66
67 @Test
68 public void testGetFrame() {
69 doTestGetFrame(Binary64Field.getInstance());
70 }
71
72 @Test
73 public void testAdditionalState() {
74 doTestAdditionalState(Binary64Field.getInstance());
75 }
76
77 @Test
78 public void testNoReset() {
79 doTestNoReset(Binary64Field.getInstance());
80 }
81
82 @Test
83 public void testAdditionalDerivatives() {
84 doTestAdditionalDerivatives(Binary64Field.getInstance());
85 }
86
87
88 @Test
89 public void testIssue949() {
90 doTestIssue949(Binary64Field.getInstance());
91 }
92
93
94 private <T extends CalculusFieldElement<T>> void doTestNormalKeplerIntegration(Field<T> field) {
95 FieldOrbit<T> initialOrbit = createOrbit(field);
96 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
97
98 FieldKeplerianPropagator<T> keplerEx = new FieldKeplerianPropagator<>(initialOrbit);
99
100
101
102
103 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(3600.0);
104 final FieldEphemerisGenerator<T> generator1 = numericalPropagator.getEphemerisGenerator();
105 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
106 numericalPropagator.propagate(finalDate);
107 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
108 FieldBoundedPropagator<T> ephemeris = generator1.getGeneratedEphemeris();
109
110
111 for (double dt = 1; dt <= 3600.0; dt += 1) {
112 FieldAbsoluteDate<T> intermediateDate = initialOrbit.getDate().shiftedBy(dt);
113 FieldSpacecraftState<T> keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
114 FieldSpacecraftState<T> numericIntermediateOrbit = ephemeris.propagate(intermediateDate);
115 FieldVector3D<T> kepPosition = keplerIntermediateOrbit.getPosition();
116 FieldVector3D<T> numPosition = numericIntermediateOrbit.getPosition();
117
118 Assertions.assertEquals(0, kepPosition.subtract(numPosition).getNorm().getReal(), 5.0e-2);
119 }
120
121
122 FieldAbsoluteDate<T> intermediateDate = initialOrbit.getDate().shiftedBy(1589);
123 FieldSpacecraftState<T> keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
124 FieldSpacecraftState<T> state = keplerEx.propagate(finalDate);
125 numericalPropagator.setInitialState(state);
126 final FieldEphemerisGenerator<T> generator2 = numericalPropagator.getEphemerisGenerator();
127 numericalPropagator.propagate(initialOrbit.getDate());
128 FieldBoundedPropagator<T> invEphemeris = generator2.getGeneratedEphemeris();
129 FieldSpacecraftState<T> numericIntermediateOrbit = invEphemeris.propagate(intermediateDate);
130 FieldVector3D<T> kepPosition = keplerIntermediateOrbit.getPosition();
131 FieldVector3D<T> numPosition = numericIntermediateOrbit.getPosition();
132
133 Assertions.assertEquals(0, kepPosition.subtract(numPosition).getNorm().getReal(), 3.0e-3);
134
135 }
136
137 @Test
138 void testGetAttitudeProvider() {
139
140 final Binary64Field field = Binary64Field.getInstance();
141 final AttitudeProvider attitudeProvider = new FrameAlignedProvider(FramesFactory.getGCRF());
142 final FieldNumericalPropagator<Binary64> fieldPropagator = createPropagator(field);
143 fieldPropagator.setAttitudeProvider(attitudeProvider);
144 final FieldOrbit<Binary64> fieldOrbit = createOrbit(field);
145 fieldPropagator.setInitialState(new FieldSpacecraftState<>(fieldOrbit));
146 final FieldEphemerisGenerator<Binary64> generator = fieldPropagator.getEphemerisGenerator();
147 fieldPropagator.propagate(fieldPropagator.getInitialState().getDate().shiftedBy(10));
148
149 final FieldBoundedPropagator<Binary64> boundedPropagator = generator.getGeneratedEphemeris();
150
151 final AttitudeProvider ephemerisAttitudeProvider = boundedPropagator.getAttitudeProvider();
152 Assertions.assertInstanceOf(AttitudeProviderModifier.class, ephemerisAttitudeProvider);
153 final AttitudeProviderModifier providerModifier = (AttitudeProviderModifier) ephemerisAttitudeProvider;
154 Assertions.assertEquals(attitudeProvider, providerModifier.getUnderlyingAttitudeProvider());
155 Assertions.assertEquals(0, ephemerisAttitudeProvider.getEventDetectors().count());
156 Assertions.assertEquals(0, ephemerisAttitudeProvider.getFieldEventDetectors(Binary64Field.getInstance()).count());
157 final Orbit orbit = fieldOrbit.toOrbit();
158 final Attitude expectedAttitude = attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
159 final Attitude actualAttitude = ephemerisAttitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
160 Assertions.assertEquals(expectedAttitude.getSpin(), actualAttitude.getSpin());
161 Assertions.assertEquals(expectedAttitude.getRotationAcceleration(), actualAttitude.getRotationAcceleration());
162 Assertions.assertEquals(0., Rotation.distance(expectedAttitude.getRotation(), actualAttitude.getRotation()));
163 }
164
165 private <T extends CalculusFieldElement<T>> void doTestGetFrame(Field<T> field) {
166 FieldOrbit<T> initialOrbit = createOrbit(field);
167 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
168
169 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
170 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
171 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
172 numericalPropagator.propagate(finalDate);
173 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
174 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
175
176
177 Assertions.assertNotNull(ephemeris.getFrame());
178 Assertions.assertSame(ephemeris.getFrame(), numericalPropagator.getFrame());
179 }
180
181 private <T extends CalculusFieldElement<T>> void doTestAdditionalState(Field<T> field) {
182 final FieldOrbit<T> initialOrbit = createOrbit(field);
183 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
184
185
186 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
187 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
188 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
189 numericalPropagator.propagate(finalDate);
190 Assertions.assertTrue(numericalPropagator.getCalls() < 3200);
191 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
192 ephemeris.addAdditionalDataProvider(new FieldAdditionalDataProvider<T[], T>() {
193
194 @Override
195 public String getName() {
196 return "time-since-start";
197 }
198
199 @Override
200 public T[] getAdditionalData(FieldSpacecraftState<T> state) {
201 T[] array = MathArrays.buildArray(state.getDate().getField(), 1);
202 array[0] = state.getDate().durationFrom(initialOrbit.getDate());
203 return array;
204 }
205 });
206
207
208 FieldSpacecraftState<T> s = ephemeris.propagate(initialOrbit.getDate().shiftedBy(20.0));
209 Assertions.assertEquals(20.0, s.getAdditionalState("time-since-start")[0].getReal(), 1.0e-10);
210
211
212 try {
213
214 Method getDrivers = ephemeris.getClass().getDeclaredMethod("getParametersDrivers", (Class[]) null);
215 getDrivers.setAccessible(true);
216 Assertions.assertTrue(((List<?>) getDrivers.invoke(ephemeris, (Object[]) null)).isEmpty());
217
218 Method getMass = ephemeris.getClass().getDeclaredMethod("getMass", FieldAbsoluteDate.class);
219 getMass.setAccessible(true);
220 @SuppressWarnings("unchecked")
221 T mass = (T) getMass.invoke(ephemeris, finalDate);
222 Assertions.assertEquals(1000.0, mass.getReal(), 1.0e-10);
223
224 Method propagateOrbit = ephemeris.getClass().getDeclaredMethod("propagateOrbit",
225 FieldAbsoluteDate.class,
226 CalculusFieldElement[].class);
227 propagateOrbit.setAccessible(true);
228 @SuppressWarnings("unchecked")
229 FieldOrbit<T> orbit = (FieldOrbit<T>) propagateOrbit.invoke(ephemeris, finalDate, null);
230 Assertions.assertEquals(initialOrbit.getA().getReal(), orbit.getA().getReal(), 1.0e-10);
231
232 } catch (IllegalAccessException | NoSuchMethodException | SecurityException |
233 IllegalArgumentException | InvocationTargetException e) {
234 Assertions.fail(e.getLocalizedMessage());
235 }
236
237 }
238
239 private <T extends CalculusFieldElement<T>> void doTestNoReset(Field<T> field) {
240 final FieldOrbit<T> initialOrbit = createOrbit(field);
241 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
242
243
244 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
245 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
246 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
247 numericalPropagator.propagate(finalDate);
248 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
249 FieldDateDetector<T> detector = new FieldDateDetector<>(initialOrbit.getDate().getField(),
250 initialOrbit.getDate().shiftedBy(10)).
251 withHandler((s, d, increasing) -> Action.RESET_STATE);
252 ephemeris.addEventDetector(detector);
253
254 try {
255 ephemeris.propagate(initialOrbit.getDate(), finalDate);
256 Assertions.fail("an exception should habe been thrown");
257 } catch (OrekitException oe) {
258 Assertions.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
259 }
260
261 }
262
263 private <T extends CalculusFieldElement<T>> void doTestAdditionalDerivatives(final Field<T> field) {
264
265 final FieldOrbit<T> initialOrbit = createOrbit(field);
266 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(10.0);
267 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(1e-3).getTolerances(
268 initialOrbit, OrbitType.CARTESIAN);
269 DormandPrince853FieldIntegrator<T> integrator =
270 new DormandPrince853FieldIntegrator<>(field, 1.0e-6, 10.0, tolerances[0], tolerances[1]);
271 integrator.setInitialStepSize(1.0e-3);
272 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
273 final DerivativesProvider<T> provider1 = new DerivativesProvider<>(field, "provider-1", 3);
274 propagator.addAdditionalDerivativesProvider(provider1);
275 final DerivativesProvider<T> provider2 = new DerivativesProvider<>(field, "provider-2", 1);
276 propagator.addAdditionalDerivativesProvider(provider2);
277 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
278 propagator.setInitialState(new FieldSpacecraftState<>(initialOrbit).
279 addAdditionalData(provider1.getName(), MathArrays.buildArray(field, provider1.getDimension())).
280 addAdditionalData(provider2.getName(), MathArrays.buildArray(field, provider2.getDimension())));
281 propagator.propagate(finalDate);
282 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
283
284 for (double dt = 0; dt < ephemeris.getMaxDate().durationFrom(ephemeris.getMinDate()).getReal(); dt += 0.1) {
285 FieldSpacecraftState<T> state = ephemeris.propagate(ephemeris.getMinDate().shiftedBy(dt));
286 checkState(dt, state, provider1);
287 checkState(dt, state, provider2);
288 }
289
290 }
291
292 private <T extends CalculusFieldElement<T>> void doTestIssue949(Field<T> field) {
293
294 final FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field);
295 final FieldOrbit<T> initialOrbit = createOrbit(field);
296 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
297 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
298 numericalPropagator.setOrbitType(OrbitType.CARTESIAN);
299
300
301 final FieldAdditionalDataProvider<T[], T> additionalDataProvider = TestUtils.getFieldAdditionalProviderWithInit();
302 numericalPropagator.addAdditionalDataProvider(additionalDataProvider);
303
304
305 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
306 numericalPropagator.propagate(initialDate.shiftedBy(1));
307
308 final FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
309
310
311 Assertions.assertDoesNotThrow(() -> ephemeris.propagate(ephemeris.getMaxDate()), "No error should have been thrown");
312
313 }
314
315 private <T extends CalculusFieldElement<T>> void checkState(final double dt, final FieldSpacecraftState<T> state,
316 final DerivativesProvider<T> provider) {
317
318 Assertions.assertTrue(state.hasAdditionalData(provider.getName()));
319 Assertions.assertEquals(provider.getDimension(), state.getAdditionalState(provider.getName()).length);
320 for (int i = 0; i < provider.getDimension(); ++i) {
321 Assertions.assertEquals(i * dt,
322 state.getAdditionalState(provider.getName())[i].getReal(),
323 4.0e-15 * i * dt);
324 }
325
326 Assertions.assertTrue(state.hasAdditionalStateDerivative(provider.getName()));
327 Assertions.assertEquals(provider.getDimension(), state.getAdditionalStateDerivative(provider.getName()).length);
328 for (int i = 0; i < provider.getDimension(); ++i) {
329 Assertions.assertEquals(i,
330 state.getAdditionalStateDerivative(provider.getName())[i].getReal(),
331 2.0e-14 * i);
332 }
333
334 }
335
336 @BeforeEach
337 public void setUp() {
338 mu = 3.9860047e14;
339 Utils.setDataRoot("regular-data:potential/icgem-format");
340 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
341 }
342
343 private <T extends CalculusFieldElement<T>> FieldNumericalPropagator<T> createPropagator(Field<T> field) {
344 double[] absTolerance= {
345 0.0001, 1.0e-11, 1.0e-11, 1.0e-8, 1.0e-8, 1.0e-8, 0.001
346 };
347 double[] relTolerance = {
348 1.0e-8, 1.0e-8, 1.0e-8, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-7
349 };
350 OrbitType type = OrbitType.EQUINOCTIAL;
351 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 500, absTolerance, relTolerance);
352 integrator.setInitialStepSize(100);
353 FieldNumericalPropagator<T> numericalPropagator = new FieldNumericalPropagator<>(field, integrator);
354 numericalPropagator.setOrbitType(type);
355 return numericalPropagator;
356 }
357
358 private <T extends CalculusFieldElement<T>> FieldOrbit<T> createOrbit(Field<T> field) {
359 T zero = field.getZero();
360 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
361 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
362
363 double mu = 3.9860047e14;
364 FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
365 return new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
366 FramesFactory.getEME2000(), initDate, zero.add(mu));
367 }
368
369 double mu;
370
371 private static class DerivativesProvider<T extends CalculusFieldElement<T>> implements FieldAdditionalDerivativesProvider<T> {
372 private final String name;
373 private final T[] derivatives;
374 DerivativesProvider(final Field<T> field, final String name, final int dimension) {
375 this.name = name;
376 this.derivatives = MathArrays.buildArray(field, dimension);
377 for (int i = 0; i < dimension; ++i) {
378 derivatives[i] = field.getZero().newInstance(i);
379 }
380 }
381 public String getName() {
382 return name;
383 }
384 public int getDimension() {
385 return derivatives.length;
386 }
387 public FieldCombinedDerivatives<T> combinedDerivatives(final FieldSpacecraftState<T> s) {
388 return new FieldCombinedDerivatives<>(derivatives, null);
389 }
390 }
391
392 }