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.ode.events.Action;
27 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
28 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
29 import org.hipparchus.util.Decimal64Field;
30 import org.hipparchus.util.MathArrays;
31 import org.junit.Assert;
32 import org.junit.Before;
33 import org.junit.Test;
34 import org.orekit.Utils;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.errors.OrekitMessages;
37 import org.orekit.forces.gravity.potential.GravityFieldFactory;
38 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
39 import org.orekit.frames.FramesFactory;
40 import org.orekit.orbits.FieldEquinoctialOrbit;
41 import org.orekit.orbits.FieldOrbit;
42 import org.orekit.orbits.OrbitType;
43 import org.orekit.propagation.FieldAdditionalStateProvider;
44 import org.orekit.propagation.FieldBoundedPropagator;
45 import org.orekit.propagation.FieldEphemerisGenerator;
46 import org.orekit.propagation.FieldSpacecraftState;
47 import org.orekit.propagation.analytical.FieldKeplerianPropagator;
48 import org.orekit.propagation.events.FieldDateDetector;
49 import org.orekit.propagation.numerical.FieldNumericalPropagator;
50 import org.orekit.time.FieldAbsoluteDate;
51 import org.orekit.utils.Constants;
52 import org.orekit.utils.FieldPVCoordinates;
53
54
55 public class FieldIntegratedEphemerisTest {
56
57 @Test
58 public void testNormalKeplerIntegration() {
59 doTestNormalKeplerIntegration(Decimal64Field.getInstance());
60 }
61
62 @Test
63 public void testGetFrame() {
64 doTestGetFrame(Decimal64Field.getInstance());
65 }
66
67 @Test
68 public void testAdditionalState() {
69 doTestAdditionalState(Decimal64Field.getInstance());
70 }
71
72 @Test
73 public void testNoReset() {
74 doTestNoReset(Decimal64Field.getInstance());
75 }
76
77 private <T extends CalculusFieldElement<T>> void doTestNormalKeplerIntegration(Field<T> field) {
78 FieldOrbit<T> initialOrbit = createOrbit(field);
79 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
80
81 FieldKeplerianPropagator<T> keplerEx = new FieldKeplerianPropagator<>(initialOrbit);
82
83
84
85
86 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(3600.0);
87 final FieldEphemerisGenerator<T> generator1 = numericalPropagator.getEphemerisGenerator();
88 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
89 numericalPropagator.propagate(finalDate);
90 Assert.assertTrue(numericalPropagator.getCalls() < 3200);
91 FieldBoundedPropagator<T> ephemeris = generator1.getGeneratedEphemeris();
92
93
94 for (double dt = 1; dt <= 3600.0; dt += 1) {
95 FieldAbsoluteDate<T> intermediateDate = initialOrbit.getDate().shiftedBy(dt);
96 FieldSpacecraftState<T> keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
97 FieldSpacecraftState<T> numericIntermediateOrbit = ephemeris.propagate(intermediateDate);
98 FieldVector3D<T> kepPosition = keplerIntermediateOrbit.getPVCoordinates().getPosition();
99 FieldVector3D<T> numPosition = numericIntermediateOrbit.getPVCoordinates().getPosition();
100
101 Assert.assertEquals(0, kepPosition.subtract(numPosition).getNorm().getReal(), 5.0e-2);
102 }
103
104
105 FieldAbsoluteDate<T> intermediateDate = initialOrbit.getDate().shiftedBy(1589);
106 FieldSpacecraftState<T> keplerIntermediateOrbit = keplerEx.propagate(intermediateDate);
107 FieldSpacecraftState<T> state = keplerEx.propagate(finalDate);
108 numericalPropagator.setInitialState(state);
109 final FieldEphemerisGenerator<T> generator2 = numericalPropagator.getEphemerisGenerator();
110 numericalPropagator.propagate(initialOrbit.getDate());
111 FieldBoundedPropagator<T> invEphemeris = generator2.getGeneratedEphemeris();
112 FieldSpacecraftState<T> numericIntermediateOrbit = invEphemeris.propagate(intermediateDate);
113 FieldVector3D<T> kepPosition = keplerIntermediateOrbit.getPVCoordinates().getPosition();
114 FieldVector3D<T> numPosition = numericIntermediateOrbit.getPVCoordinates().getPosition();
115
116 Assert.assertEquals(0, kepPosition.subtract(numPosition).getNorm().getReal(), 3.0e-3);
117
118 }
119
120 private <T extends CalculusFieldElement<T>> void doTestGetFrame(Field<T> field) {
121 FieldOrbit<T> initialOrbit = createOrbit(field);
122 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
123
124 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
125 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
126 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
127 numericalPropagator.propagate(finalDate);
128 Assert.assertTrue(numericalPropagator.getCalls() < 3200);
129 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
130
131
132 Assert.assertNotNull(ephemeris.getFrame());
133 Assert.assertSame(ephemeris.getFrame(), numericalPropagator.getFrame());
134 }
135
136 private <T extends CalculusFieldElement<T>> void doTestAdditionalState(Field<T> field) {
137 final FieldOrbit<T> initialOrbit = createOrbit(field);
138 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
139
140
141 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
142 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
143 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
144 numericalPropagator.propagate(finalDate);
145 Assert.assertTrue(numericalPropagator.getCalls() < 3200);
146 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
147 ephemeris.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {
148
149 @Override
150 public String getName() {
151 return "time-since-start";
152 }
153
154 @Override
155 public T[] getAdditionalState(FieldSpacecraftState<T> state) {
156 T[] array = MathArrays.buildArray(state.getDate().getField(), 1);
157 array[0] = state.getDate().durationFrom(initialOrbit.getDate());
158 return array;
159 }
160 });
161
162
163 FieldSpacecraftState<T> s = ephemeris.propagate(initialOrbit.getDate().shiftedBy(20.0));
164 Assert.assertEquals(20.0, s.getAdditionalState("time-since-start")[0].getReal(), 1.0e-10);
165
166
167 try {
168
169 Method getDrivers = ephemeris.getClass().getDeclaredMethod("getParametersDrivers", (Class[]) null);
170 getDrivers.setAccessible(true);
171 Assert.assertTrue(((List<?>) getDrivers.invoke(ephemeris, (Object[]) null)).isEmpty());
172
173 Method getMass = ephemeris.getClass().getDeclaredMethod("getMass", FieldAbsoluteDate.class);
174 getMass.setAccessible(true);
175 @SuppressWarnings("unchecked")
176 T mass = (T) getMass.invoke(ephemeris, finalDate);
177 Assert.assertEquals(1000.0, mass.getReal(), 1.0e-10);
178
179 Method propagateOrbit = ephemeris.getClass().getDeclaredMethod("propagateOrbit",
180 FieldAbsoluteDate.class,
181 CalculusFieldElement[].class);
182 propagateOrbit.setAccessible(true);
183 @SuppressWarnings("unchecked")
184 FieldOrbit<T> orbit = (FieldOrbit<T>) propagateOrbit.invoke(ephemeris, finalDate, null);
185 Assert.assertEquals(initialOrbit.getA().getReal(), orbit.getA().getReal(), 1.0e-10);
186
187 } catch (IllegalAccessException | NoSuchMethodException | SecurityException |
188 IllegalArgumentException | InvocationTargetException e) {
189 Assert.fail(e.getLocalizedMessage());
190 }
191
192 }
193
194 private <T extends CalculusFieldElement<T>> void doTestNoReset(Field<T> field) {
195 final FieldOrbit<T> initialOrbit = createOrbit(field);
196 FieldNumericalPropagator<T> numericalPropagator = createPropagator(field);
197
198
199 FieldAbsoluteDate<T> finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
200 final FieldEphemerisGenerator<T> generator = numericalPropagator.getEphemerisGenerator();
201 numericalPropagator.setInitialState(new FieldSpacecraftState<>(initialOrbit));
202 numericalPropagator.propagate(finalDate);
203 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
204 ephemeris.addEventDetector(new FieldDateDetector<>(initialOrbit.getDate().shiftedBy(10)).
205 withHandler((FieldSpacecraftState<T> s,
206 FieldDateDetector<T> detector,
207 boolean increasing) -> Action.RESET_STATE));
208
209 try {
210 ephemeris.propagate(initialOrbit.getDate(), finalDate);
211 Assert.fail("an exception should habe been thrown");
212 } catch (OrekitException oe) {
213 Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
214 }
215
216 }
217
218 @Before
219 public void setUp() {
220 mu = 3.9860047e14;
221 Utils.setDataRoot("regular-data:potential/icgem-format");
222 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
223 }
224
225 private <T extends CalculusFieldElement<T>> FieldNumericalPropagator<T> createPropagator(Field<T> field) {
226 double[] absTolerance= {
227 0.0001, 1.0e-11, 1.0e-11, 1.0e-8, 1.0e-8, 1.0e-8, 0.001
228 };
229 double[] relTolerance = {
230 1.0e-8, 1.0e-8, 1.0e-8, 1.0e-9, 1.0e-9, 1.0e-9, 1.0e-7
231 };
232 OrbitType type = OrbitType.EQUINOCTIAL;
233 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 500, absTolerance, relTolerance);
234 integrator.setInitialStepSize(100);
235 FieldNumericalPropagator<T> numericalPropagator = new FieldNumericalPropagator<>(field, integrator);
236 numericalPropagator.setOrbitType(type);
237 return numericalPropagator;
238 }
239
240 private <T extends CalculusFieldElement<T>> FieldOrbit<T> createOrbit(Field<T> field) {
241 T zero = field.getZero();
242 FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
243 FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
244
245 double mu = 3.9860047e14;
246 FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
247 return new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
248 FramesFactory.getEME2000(), initDate, zero.add(mu));
249 }
250
251 double mu;
252 }