1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.numerical;
18
19 import java.lang.reflect.Array;
20
21 import org.hamcrest.MatcherAssert;
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
26 import org.hipparchus.ode.FieldODEIntegrator;
27 import org.hipparchus.ode.events.Action;
28 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
29 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
30 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
31 import org.hipparchus.util.Decimal64Field;
32 import org.hipparchus.util.FastMath;
33 import org.hipparchus.util.MathArrays;
34 import org.junit.Assert;
35 import org.junit.Before;
36 import org.junit.Test;
37 import org.orekit.OrekitMatchers;
38 import org.orekit.Utils;
39 import org.orekit.bodies.CelestialBodyFactory;
40 import org.orekit.bodies.OneAxisEllipsoid;
41 import org.orekit.data.DataContext;
42 import org.orekit.errors.OrekitException;
43 import org.orekit.errors.OrekitMessages;
44 import org.orekit.forces.ForceModel;
45 import org.orekit.forces.drag.DragForce;
46 import org.orekit.forces.drag.IsotropicDrag;
47 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
48 import org.orekit.forces.gravity.ThirdBodyAttraction;
49 import org.orekit.forces.gravity.potential.GRGSFormatReader;
50 import org.orekit.forces.gravity.potential.GravityFieldFactory;
51 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
52 import org.orekit.forces.gravity.potential.SHMFormatReader;
53 import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
54 import org.orekit.forces.radiation.SolarRadiationPressure;
55 import org.orekit.frames.Frame;
56 import org.orekit.frames.FramesFactory;
57 import org.orekit.models.earth.atmosphere.DTM2000;
58 import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
59 import org.orekit.orbits.FieldCartesianOrbit;
60 import org.orekit.orbits.FieldEquinoctialOrbit;
61 import org.orekit.orbits.FieldKeplerianOrbit;
62 import org.orekit.orbits.FieldOrbit;
63 import org.orekit.orbits.OrbitType;
64 import org.orekit.orbits.PositionAngle;
65 import org.orekit.propagation.FieldAdditionalStateProvider;
66 import org.orekit.propagation.FieldBoundedPropagator;
67 import org.orekit.propagation.FieldEphemerisGenerator;
68 import org.orekit.propagation.FieldSpacecraftState;
69 import org.orekit.propagation.events.FieldAbstractDetector;
70 import org.orekit.propagation.events.FieldApsideDetector;
71 import org.orekit.propagation.events.FieldDateDetector;
72 import org.orekit.propagation.events.FieldEventDetector;
73 import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
74 import org.orekit.propagation.events.handlers.FieldEventHandler;
75 import org.orekit.propagation.events.handlers.FieldStopOnEvent;
76 import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator;
77 import org.orekit.propagation.integration.FieldAdditionalDerivativesProvider;
78 import org.orekit.propagation.sampling.FieldOrekitStepHandler;
79 import org.orekit.propagation.sampling.FieldOrekitStepInterpolator;
80 import org.orekit.time.FieldAbsoluteDate;
81 import org.orekit.time.FieldTimeStamped;
82 import org.orekit.time.TimeScale;
83 import org.orekit.time.TimeScalesFactory;
84 import org.orekit.utils.Constants;
85 import org.orekit.utils.FieldPVCoordinates;
86 import org.orekit.utils.IERSConventions;
87 import org.orekit.utils.TimeStampedFieldPVCoordinates;
88
89
90 public class FieldNumericalPropagatorTest {
91
92 private double mu;
93
94 @Test(expected=OrekitException.class)
95 public void testNotInitialised1() {
96 doTestNotInitialised1(Decimal64Field.getInstance());
97 }
98
99 private <T extends CalculusFieldElement<T>> void doTestNotInitialised1(Field<T> field) {
100
101 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
102
103 final FieldAbstractIntegratedPropagator<T> notInitialised =
104 new FieldNumericalPropagator<>(field, new ClassicalRungeKuttaFieldIntegrator<>(field, field.getZero().add(10.0)));
105 notInitialised.propagate(initDate);
106 }
107
108 @Test(expected=OrekitException.class)
109 public void testNotInitialised2() {
110 doTestNotInitialised2(Decimal64Field.getInstance());
111 }
112
113 private <T extends CalculusFieldElement<T>> void doTestNotInitialised2(Field<T> field) {
114
115 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
116 final FieldAbstractIntegratedPropagator<T> notInitialised =
117 new FieldNumericalPropagator<>(field, new ClassicalRungeKuttaFieldIntegrator<>(field, field.getZero().add((10.0))));
118 notInitialised.propagate(initDate, initDate.shiftedBy(3600));
119 }
120
121 @Test
122 public void testEventAtEndOfEphemeris() {
123 doTestEventAtEndOfEphemeris(Decimal64Field.getInstance());
124 }
125
126 private <T extends CalculusFieldElement<T>, D extends FieldEventDetector<T>> void doTestEventAtEndOfEphemeris(Field<T> field) {
127
128 T zero = field.getZero();
129 FieldNumericalPropagator<T> propagator = createPropagator(field);
130 FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
131
132
133 FieldAbsoluteDate<T> end = initDate.shiftedBy(100)
134 .shiftedBy(3 * FastMath.ulp(100.0) / 4);
135 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
136 propagator.propagate(end);
137 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
138 CountingHandler<D, T> handler = new CountingHandler<D, T>();
139 FieldDateDetector<T> detector = new FieldDateDetector<>(zero.add(10), zero.add(1e-9), toArray(end)).withHandler(handler);
140
141 ephemeris.addEventDetector(detector);
142
143
144
145 FieldSpacecraftState<T> actual = ephemeris.propagate(end);
146
147
148 Assert.assertEquals(actual.getDate().durationFrom(end).getReal(), 0.0, 0.0);
149 Assert.assertEquals(1, handler.eventCount);
150 }
151
152 @Test
153 public void testEventAtBeginningOfEphemeris() {
154 doTestEventAtBeginningOfEphemeris(Decimal64Field.getInstance());
155 }
156
157 private <T extends CalculusFieldElement<T>, D extends FieldEventDetector<T>> void doTestEventAtBeginningOfEphemeris(Field<T> field) {
158
159 T zero = field.getZero();
160 FieldNumericalPropagator<T> propagator = createPropagator(field);
161 FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
162
163
164
165 FieldAbsoluteDate<T> end = initDate.shiftedBy(100)
166 .shiftedBy(3 * FastMath.ulp(100.0) / 4);
167 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
168 propagator.propagate(end);
169 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
170 CountingHandler<D, T> handler = new CountingHandler<D, T>();
171
172
173 FieldAbsoluteDate<T> eventDate = initDate.shiftedBy(FastMath.ulp(100.0) / 10.0);
174 FieldDateDetector<T> detector = new FieldDateDetector<>(zero.add(10), zero.add(1e-9), toArray(eventDate))
175 .withHandler(handler);
176
177 ephemeris.addEventDetector(detector);
178
179
180
181 Assert.assertEquals(ephemeris.propagate(end).getDate().durationFrom(end).getReal(), 0.0, 0.0);
182
183 Assert.assertEquals(ephemeris.propagate(initDate).getDate().durationFrom(initDate).getReal(), 0.0, 0.0);
184 Assert.assertEquals(2, handler.eventCount);
185 }
186
187 public class CountingHandler <D extends FieldEventDetector<T>, T extends CalculusFieldElement<T>>
188 implements FieldEventHandler<FieldEventDetector<T>, T> {
189
190
191
192
193
194 private int eventCount = 0;
195
196 @Override
197 public Action eventOccurred(FieldSpacecraftState<T> s,
198 FieldEventDetector<T> detector,
199 boolean increasing) {
200 eventCount++;
201 return Action.CONTINUE;
202 }
203
204 @Override
205 public FieldSpacecraftState<T> resetState(FieldEventDetector<T> detector,
206 FieldSpacecraftState<T> oldState) {
207 return null;
208 }
209
210
211 }
212
213 @Test
214 public void testCloseEventDates() {
215 doTestCloseEventDates(Decimal64Field.getInstance());
216 }
217
218 private <T extends CalculusFieldElement<T>> void doTestCloseEventDates(Field<T> field) {
219
220 T zero = field.getZero();
221 FieldNumericalPropagator<T> propagator = createPropagator(field);
222 FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
223
224
225 FieldDateDetector<T> d1 = new FieldDateDetector<>(zero.add(10), zero.add(1), toArray(initDate.shiftedBy(15)))
226 .withHandler(new FieldContinueOnEvent<FieldDateDetector<T>, T>());
227 FieldDateDetector<T> d2 = new FieldDateDetector<>(zero.add(10), zero.add(1), toArray(initDate.shiftedBy(15.5)))
228 .withHandler(new FieldContinueOnEvent<FieldDateDetector<T>, T>());
229 propagator.addEventDetector(d1);
230 propagator.addEventDetector(d2);
231
232
233 FieldAbsoluteDate<T> end = initDate.shiftedBy(30);
234 FieldSpacecraftState<T> actual = propagator.propagate(end);
235
236
237 Assert.assertEquals(actual.getDate().durationFrom(end).getReal(), 0.0, 0.0);
238 }
239
240 @Test
241 public void testEphemerisDates() {
242 doTestEphemerisDates(Decimal64Field.getInstance());
243 }
244
245 private <T extends CalculusFieldElement<T>> void doTestEphemerisDates(Field<T> field) {
246
247 T zero = field.getZero();
248
249
250
251 TimeScale tai = TimeScalesFactory.getTAI();
252 FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, "2015-07-01", tai);
253 FieldAbsoluteDate<T> startDate = new FieldAbsoluteDate<>(field, "2015-07-03", tai).shiftedBy(-0.1);
254 FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(field, "2015-07-04", tai);
255 Frame eci = FramesFactory.getGCRF();
256 FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS), zero, zero, zero, zero, zero,
257 PositionAngle.TRUE, eci, initialDate, zero.add(mu));
258 OrbitType type = OrbitType.CARTESIAN;
259 double[][] tol = NumericalPropagator.tolerances(1e-3, orbit.toOrbit(), type);
260 FieldNumericalPropagator<T> prop = new FieldNumericalPropagator<>(field,
261 new DormandPrince853FieldIntegrator<>(field, 0.1, 500, tol[0], tol[1]));
262 prop.setOrbitType(type);
263 prop.resetInitialState(new FieldSpacecraftState<>(new FieldCartesianOrbit<>(orbit)));
264
265
266 final FieldEphemerisGenerator<T> generator = prop.getEphemerisGenerator();
267 prop.propagate(startDate, endDate);
268 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
269
270
271 TimeStampedFieldPVCoordinates<T> actualPV = ephemeris.getPVCoordinates(startDate, eci);
272 TimeStampedFieldPVCoordinates<T> expectedPV = orbit.getPVCoordinates(startDate, eci);
273 MatcherAssert.assertThat(actualPV.getPosition().toVector3D(),
274 OrekitMatchers.vectorCloseTo(expectedPV.getPosition().toVector3D(), 1.0));
275 MatcherAssert.assertThat(actualPV.getVelocity().toVector3D(),
276 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity().toVector3D(), 1.0));
277 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate).getReal(),
278 OrekitMatchers.closeTo(0, 0));
279 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate).getReal(),
280 OrekitMatchers.closeTo(0, 0));
281
282 FieldAbsoluteDate<T> date = endDate.shiftedBy(-0.11);
283 Assert.assertEquals(
284 ephemeris.propagate(date).getDate().durationFrom(date).getReal(), 0, 0);
285
286 Assert.assertTrue(prop.getAdditionalDerivativesProviders().isEmpty());
287
288 }
289
290 @Test
291 public void testEphemerisDatesBackward() {
292 doTestEphemerisDatesBackward(Decimal64Field.getInstance());
293 }
294
295 private <T extends CalculusFieldElement<T>> void doTestEphemerisDatesBackward(Field<T> field) {
296
297 T zero = field.getZero();
298
299 TimeScale tai = TimeScalesFactory.getTAI();
300 FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, "2015-07-05", tai);
301 FieldAbsoluteDate<T> startDate = new FieldAbsoluteDate<>(field, "2015-07-03", tai).shiftedBy(-0.1);
302 FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(field, "2015-07-04", tai);
303 Frame eci = FramesFactory.getGCRF();
304 FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS), zero, zero, zero, zero, zero,
305 PositionAngle.TRUE, eci, initialDate, zero.add(mu));
306 OrbitType type = OrbitType.CARTESIAN;
307 double[][] tol = NumericalPropagator.tolerances(1e-3, orbit.toOrbit(), type);
308 FieldNumericalPropagator<T> prop = new FieldNumericalPropagator<>(field,
309 new DormandPrince853FieldIntegrator<>(field, 0.1, 500, tol[0], tol[1]));
310 prop.setOrbitType(type);
311 prop.resetInitialState(new FieldSpacecraftState<>(new FieldCartesianOrbit<>(orbit)));
312
313
314 final FieldEphemerisGenerator<T> generator = prop.getEphemerisGenerator();
315 prop.propagate(endDate, startDate);
316 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
317
318
319 TimeStampedFieldPVCoordinates<T> actualPV = ephemeris.getPVCoordinates(startDate, eci);
320 TimeStampedFieldPVCoordinates<T> expectedPV = orbit.getPVCoordinates(startDate, eci);
321 MatcherAssert.assertThat(actualPV.getPosition().toVector3D(),
322 OrekitMatchers.vectorCloseTo(expectedPV.getPosition().toVector3D(), 1.0));
323 MatcherAssert.assertThat(actualPV.getVelocity().toVector3D(),
324 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity().toVector3D(), 1.0));
325 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate).getReal(),
326 OrekitMatchers.closeTo(0, 0));
327 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate).getReal(),
328 OrekitMatchers.closeTo(0, 0));
329
330 FieldAbsoluteDate<T> date = endDate.shiftedBy(-0.11);
331 Assert.assertEquals(
332 ephemeris.propagate(date).getDate().durationFrom(date).getReal(), 0, 0);
333 }
334
335 @Test
336 public void testNoExtrapolation() {
337 doTestNoExtrapolation(Decimal64Field.getInstance());
338 }
339
340 private <T extends CalculusFieldElement<T>> void doTestNoExtrapolation(Field<T> field) {
341
342 T zero = field.getZero();
343
344 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
345 FieldSpacecraftState<T> initialState;
346 FieldNumericalPropagator<T> propagator;
347 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
348 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
349
350 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
351 FramesFactory.getEME2000(), initDate, zero.add(mu));
352 initialState = new FieldSpacecraftState<>(orbit);
353 OrbitType type = OrbitType.EQUINOCTIAL;
354 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
355 AdaptiveStepsizeFieldIntegrator<T>integrator =
356 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
357 integrator.setInitialStepSize(60);
358 propagator = new FieldNumericalPropagator<>(field, integrator);
359 propagator.setOrbitType(type);
360 propagator.setInitialState(initialState);
361
362
363
364
365
366
367 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate);
368
369 final FieldVector3D<T> initialPosition = initialState.getPVCoordinates().getPosition();
370 final FieldVector3D<T> initialVelocity = initialState.getPVCoordinates().getVelocity();
371
372
373 final FieldVector3D<T> finalPosition = finalState.getPVCoordinates().getPosition();
374 final FieldVector3D<T> finalVelocity = finalState.getPVCoordinates().getVelocity();
375
376
377 Assert.assertEquals(initialPosition.getX().getReal(), finalPosition.getX().getReal(), 1.0e-10);
378 Assert.assertEquals(initialPosition.getY().getReal(), finalPosition.getY().getReal(), 1.0e-10);
379 Assert.assertEquals(initialPosition.getZ().getReal(), finalPosition.getZ().getReal(), 1.0e-10);
380 Assert.assertEquals(initialVelocity.getX().getReal(), finalVelocity.getX().getReal(), 1.0e-10);
381 Assert.assertEquals(initialVelocity.getY().getReal(), finalVelocity.getY().getReal(), 1.0e-10);
382 Assert.assertEquals(initialVelocity.getZ().getReal(), finalVelocity.getZ().getReal(), 1.0e-10);
383
384 }
385
386 @Test
387 public void testKepler() {
388 doTestKepler(Decimal64Field.getInstance());
389 }
390
391 private <T extends CalculusFieldElement<T>> void doTestKepler(Field<T> field) {
392 T zero = field.getZero();
393
394 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
395 FieldSpacecraftState<T> initialState;
396 FieldNumericalPropagator<T> propagator;
397 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
398 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
399
400 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
401 FramesFactory.getEME2000(), initDate, zero.add(mu));
402 initialState = new FieldSpacecraftState<>(orbit);
403 OrbitType type = OrbitType.EQUINOCTIAL;
404 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
405 AdaptiveStepsizeFieldIntegrator<T>integrator =
406 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
407 integrator.setInitialStepSize(60);
408 propagator = new FieldNumericalPropagator<>(field, integrator);
409 propagator.setOrbitType(type);
410 propagator.setInitialState(initialState);
411
412
413 final double dt = 3200;
414 final FieldSpacecraftState<T> finalState =
415 propagator.propagate(initDate.shiftedBy(-60), initDate.shiftedBy(dt));
416
417
418 final double n = FastMath.sqrt(initialState.getMu().divide(initialState.getA())).getReal() / initialState.getA().getReal();
419 Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 1.0e-10);
420 Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-10);
421 Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 1.0e-10);
422 Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-10);
423 Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 1.0e-10);
424 Assert.assertEquals(initialState.getLM().getReal() + n * dt, finalState.getLM().getReal(), 2.0e-9);
425
426 }
427
428 @Test
429 public void testCartesian() {
430 doTestCartesian(Decimal64Field.getInstance());
431 }
432
433 private <T extends CalculusFieldElement<T>> void doTestCartesian(Field<T> field) {
434 T zero = field.getZero();
435
436 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
437 FieldSpacecraftState<T> initialState;
438 FieldNumericalPropagator<T> propagator;
439 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
440 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
441
442 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
443 FramesFactory.getEME2000(), initDate, zero.add(mu));
444 initialState = new FieldSpacecraftState<>(orbit);
445 OrbitType type = OrbitType.EQUINOCTIAL;
446 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
447 AdaptiveStepsizeFieldIntegrator<T>integrator =
448 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
449 integrator.setInitialStepSize(60);
450 propagator = new FieldNumericalPropagator<>(field, integrator);
451 propagator.setOrbitType(type);
452 propagator.setInitialState(initialState);
453
454
455 final T dt = zero.add(3200);
456 propagator.setOrbitType(OrbitType.CARTESIAN);
457 final FieldPVCoordinates<T> finalState =
458 propagator.propagate(initDate.shiftedBy(dt)).getPVCoordinates();
459 final FieldVector3D<T> pFin = finalState.getPosition();
460 final FieldVector3D<T> vFin = finalState.getVelocity();
461
462
463 final FieldPVCoordinates<T> reference = initialState.shiftedBy(dt).getPVCoordinates();
464 final FieldVector3D<T> pRef = reference.getPosition();
465 final FieldVector3D<T> vRef = reference.getVelocity();
466 Assert.assertEquals(0, pRef.subtract(pFin).getNorm().getReal(), 2e-4);
467 Assert.assertEquals(0, vRef.subtract(vFin).getNorm().getReal(), 7e-8);
468
469 }
470
471 @Test
472 public void testPropagationTypesElliptical() {
473 doTestPropagationTypesElliptical(Decimal64Field.getInstance());
474 }
475
476 private <T extends CalculusFieldElement<T>> void doTestPropagationTypesElliptical(Field<T> field) {
477
478 T zero = field.getZero();
479
480 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
481 FieldSpacecraftState<T> initialState;
482 FieldNumericalPropagator<T> propagator;
483 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
484 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
485
486 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
487 FramesFactory.getEME2000(), initDate, zero.add(mu));
488 initialState = new FieldSpacecraftState<>(orbit);
489 OrbitType type = OrbitType.EQUINOCTIAL;
490 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
491 AdaptiveStepsizeFieldIntegrator<T>integrator =
492 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
493 integrator.setInitialStepSize(60);
494 propagator = new FieldNumericalPropagator<>(field, integrator);
495 propagator.setOrbitType(type);
496 propagator.setInitialState(initialState);
497
498 ForceModel gravityField =
499 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
500 GravityFieldFactory.getNormalizedProvider(5, 5));
501 propagator.addForceModel(gravityField);
502
503
504 final FieldPVCoordinates<T> pv = initialState.getPVCoordinates();
505 final T dP = zero.add(0.001);
506 final T dV = pv.getPosition().getNormSq().multiply(pv.getVelocity().getNorm()).reciprocal().multiply(dP.multiply(initialState.getMu()));
507
508 final FieldPVCoordinates<T> pvcM = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.MEAN, propagator);
509 final FieldPVCoordinates<T> pviM = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.MEAN, propagator);
510 final FieldPVCoordinates<T> pveM = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.MEAN, propagator);
511 final FieldPVCoordinates<T> pvkM = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.MEAN, propagator);
512
513 final FieldPVCoordinates<T> pvcE = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, propagator);
514 final FieldPVCoordinates<T> pviE = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.ECCENTRIC, propagator);
515 final FieldPVCoordinates<T> pveE = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC, propagator);
516 final FieldPVCoordinates<T> pvkE = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, propagator);
517
518 final FieldPVCoordinates<T> pvcT = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.TRUE, propagator);
519 final FieldPVCoordinates<T> pviT = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.TRUE, propagator);
520 final FieldPVCoordinates<T> pveT = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.TRUE, propagator);
521 final FieldPVCoordinates<T> pvkT = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.TRUE, propagator);
522 Assert.assertEquals(0, pvcM.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 3.0);
523 Assert.assertEquals(0, pvcM.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 2.0);
524 Assert.assertEquals(0, pviM.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.6);
525 Assert.assertEquals(0, pviM.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
526 Assert.assertEquals(0, pvkM.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.5);
527 Assert.assertEquals(0, pvkM.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.3);
528 Assert.assertEquals(0, pveM.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.2);
529 Assert.assertEquals(0, pveM.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.2);
530
531 Assert.assertEquals(0, pvcE.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 3.0);
532 Assert.assertEquals(0, pvcE.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 2.0);
533
534 Assert.assertEquals(0, pviE.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.03);
535 Assert.assertEquals(0, pviE.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.04);
536 Assert.assertEquals(0, pvkE.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.4);
537 Assert.assertEquals(0, pvkE.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.3);
538 Assert.assertEquals(0, pveE.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.2);
539 Assert.assertEquals(0, pveE.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.07);
540
541 Assert.assertEquals(0, pvcT.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 3.0);
542 Assert.assertEquals(0, pvcT.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 2.0);
543 Assert.assertEquals(0, pviT.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
544 Assert.assertEquals(0, pviT.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.2);
545 Assert.assertEquals(0, pvkT.getPosition().subtract(pveT.getPosition()).getNorm().getReal() / dP.getReal(), 0.4);
546 Assert.assertEquals(0, pvkT.getVelocity().subtract(pveT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.2);
547
548 }
549
550 @Test
551 public void testPropagationTypesHyperbolic() {
552 doTestPropagationTypesHyperbolic(Decimal64Field.getInstance());
553 }
554
555 private <T extends CalculusFieldElement<T>> void doTestPropagationTypesHyperbolic(Field<T> field) {
556 T zero = field.getZero();
557
558 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
559 FieldSpacecraftState<T> initialState;
560 FieldNumericalPropagator<T> propagator;
561 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
562 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
563
564 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
565 FramesFactory.getEME2000(), initDate, zero.add(mu));
566 initialState = new FieldSpacecraftState<>(orbit);
567 OrbitType type = OrbitType.EQUINOCTIAL;
568 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
569 AdaptiveStepsizeFieldIntegrator<T>integrator =
570 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
571 integrator.setInitialStepSize(60);
572 propagator = new FieldNumericalPropagator<>(field, integrator);
573 propagator.setOrbitType(type);
574 propagator.setInitialState(initialState);
575
576 FieldSpacecraftState<T> state =
577 new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(zero.add(-10000000.0), zero.add(2.5), zero.add(0.3), zero, zero, zero,
578 PositionAngle.TRUE,
579 FramesFactory.getEME2000(), initDate,
580 zero.add(mu)));
581
582 ForceModel gravityField =
583 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
584 GravityFieldFactory.getNormalizedProvider(5, 5));
585 propagator.addForceModel(gravityField);
586
587
588 final FieldPVCoordinates<T> pv = state.getPVCoordinates();
589 final T dP = zero.add(0.001);
590 final T dV = dP.multiply(state.getMu()).divide(
591 pv.getPosition().getNormSq().multiply(pv.getVelocity().getNorm()));
592
593 final FieldPVCoordinates<T> pvcM = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.MEAN, propagator);
594 final FieldPVCoordinates<T> pvkM = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.MEAN, propagator);
595
596 final FieldPVCoordinates<T> pvcE = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, propagator);
597 final FieldPVCoordinates<T> pvkE = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, propagator);
598
599 final FieldPVCoordinates<T> pvcT = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.TRUE, propagator);
600 final FieldPVCoordinates<T> pvkT = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.TRUE, propagator);
601
602 Assert.assertEquals(0, pvcM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
603 Assert.assertEquals(0, pvcM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
604 Assert.assertEquals(0, pvkM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.2);
605 Assert.assertEquals(0, pvkM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.3);
606 Assert.assertEquals(0, pvcE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
607 Assert.assertEquals(0, pvcE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
608 Assert.assertEquals(0, pvkE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.009);
609 Assert.assertEquals(0, pvkE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.006);
610 Assert.assertEquals(0, pvcT.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
611 Assert.assertEquals(0, pvcT.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
612
613 }
614
615 private <T extends CalculusFieldElement<T>> FieldPVCoordinates<T> propagateInType(FieldSpacecraftState<T> state, T dP,
616 OrbitType type , PositionAngle angle, FieldNumericalPropagator<T> propagator)
617 {
618 T zero = dP.getField().getZero();
619 final T dt = zero.add(3200);
620 final double minStep = 0.001;
621 final double maxStep = 1000;
622 double[][] tol = NumericalPropagator.tolerances(dP.getReal(), state.getOrbit().toOrbit(), type);
623 AdaptiveStepsizeFieldIntegrator<T> integrator =
624 new DormandPrince853FieldIntegrator<>(zero.getField(), minStep, maxStep, tol[0], tol[1]);
625 FieldNumericalPropagator<T> newPropagator = new FieldNumericalPropagator<>(zero.getField(), integrator);
626 newPropagator.setOrbitType(type);
627 newPropagator.setPositionAngleType(angle);
628 newPropagator.setInitialState(state);
629 for (ForceModel force: propagator.getAllForceModels()) {
630 newPropagator.addForceModel(force);
631 }
632 return newPropagator.propagate(state.getDate().shiftedBy(dt)).getPVCoordinates();
633
634 }
635
636 @Test(expected=OrekitException.class)
637 public void testException() {
638 doTestException(Decimal64Field.getInstance());
639 }
640
641 private <T extends CalculusFieldElement<T>> void doTestException(Field<T> field) {
642 T zero = field.getZero();
643
644 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
645 FieldSpacecraftState<T> initialState;
646 FieldNumericalPropagator<T> propagator;
647 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
648 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
649
650 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
651 FramesFactory.getEME2000(), initDate, zero.add(mu));
652 initialState = new FieldSpacecraftState<>(orbit);
653 OrbitType type = OrbitType.EQUINOCTIAL;
654 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
655 AdaptiveStepsizeFieldIntegrator<T>integrator =
656 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
657 integrator.setInitialStepSize(60);
658 propagator = new FieldNumericalPropagator<>(field, integrator);
659 propagator.setOrbitType(type);
660 propagator.setInitialState(initialState);
661
662 propagator.setStepHandler(new FieldOrekitStepHandler<T>() {
663 private int countDown = 3;
664 private FieldAbsoluteDate<T> previousCall = null;
665 public void init(FieldSpacecraftState<T> s0, FieldAbsoluteDate<T> t) {
666 }
667 public void handleStep(FieldOrekitStepInterpolator<T> interpolator) {
668 if (previousCall != null) {
669 System.out.println(interpolator.getCurrentState().getDate().compareTo(previousCall) < 0);
670 }
671 if (--countDown == 0) {
672 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE, "dummy error");
673 }
674 }
675 });
676 propagator.propagate(initDate.shiftedBy(-3600));
677
678 }
679
680 @Test
681 public void testStopEvent() {
682 doTestStopEvent(Decimal64Field.getInstance());
683 }
684
685 private <T extends CalculusFieldElement<T>> void doTestStopEvent(Field<T> field) {
686 T zero = field.getZero();
687
688 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
689 FieldSpacecraftState<T> initialState;
690 FieldNumericalPropagator<T> propagator;
691 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
692 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
693
694 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
695 FramesFactory.getEME2000(), initDate, zero.add(mu));
696 initialState = new FieldSpacecraftState<>(orbit);
697 OrbitType type = OrbitType.EQUINOCTIAL;
698 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
699 AdaptiveStepsizeFieldIntegrator<T>integrator =
700 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
701 integrator.setInitialStepSize(60);
702 propagator = new FieldNumericalPropagator<>(field, integrator);
703 propagator.setOrbitType(type);
704 propagator.setInitialState(initialState);
705
706 final FieldAbsoluteDate<T> stopDate = initDate.shiftedBy(1000);
707 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.STOP);
708 propagator.addEventDetector(new FieldDateDetector<>(stopDate).withHandler(checking));
709 Assert.assertEquals(1, propagator.getEventsDetectors().size());
710 checking.assertEvent(false);
711 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(3200));
712 checking.assertEvent(true);
713 Assert.assertEquals(0, finalState.getDate().durationFrom(stopDate).getReal(), 1.0e-10);
714 propagator.clearEventsDetectors();
715 Assert.assertEquals(0, propagator.getEventsDetectors().size());
716
717 }
718
719 @Test
720 public void testResetStateEvent() {
721 doTestResetStateEvent(Decimal64Field.getInstance());
722 }
723
724 private <T extends CalculusFieldElement<T>> void doTestResetStateEvent(Field<T> field) {
725 T zero = field.getZero();
726
727 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
728 FieldSpacecraftState<T> initialState;
729 FieldNumericalPropagator<T> propagator;
730 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
731 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
732
733 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
734 FramesFactory.getEME2000(), initDate, zero.add(mu));
735 initialState = new FieldSpacecraftState<>(orbit);
736 OrbitType type = OrbitType.EQUINOCTIAL;
737 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
738 AdaptiveStepsizeFieldIntegrator<T>integrator =
739 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
740 integrator.setInitialStepSize(60);
741 propagator = new FieldNumericalPropagator<>(field, integrator);
742 propagator.setOrbitType(type);
743 propagator.setInitialState(initialState);
744 final FieldAbsoluteDate<T> resetDate = initDate.shiftedBy(1000);
745 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.RESET_STATE) {
746 public FieldSpacecraftState<T> resetState(FieldDateDetector<T> detector, FieldSpacecraftState<T> oldState) {
747 return new FieldSpacecraftState<>(oldState.getOrbit(), oldState.getAttitude(), oldState.getMass().subtract(200.0));
748 }
749 };
750 propagator.addEventDetector(new FieldDateDetector<>(resetDate).withHandler(checking));
751 checking.assertEvent(false);
752 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(3200));
753 checking.assertEvent(true);
754 Assert.assertEquals(initialState.getMass().getReal() - 200, finalState.getMass().getReal(), 1.0e-10);
755 }
756
757 @Test
758 public void testResetDerivativesEvent() {
759 doTestResetDerivativesEvent(Decimal64Field.getInstance());
760 }
761
762 private <T extends CalculusFieldElement<T>> void doTestResetDerivativesEvent(Field<T> field) {
763 T zero = field.getZero();
764
765 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
766 FieldSpacecraftState<T> initialState;
767 FieldNumericalPropagator<T> propagator;
768 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
769 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
770
771 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
772 FramesFactory.getEME2000(), initDate, zero.add(mu));
773 initialState = new FieldSpacecraftState<>(orbit);
774 OrbitType type = OrbitType.EQUINOCTIAL;
775 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
776 AdaptiveStepsizeFieldIntegrator<T>integrator =
777 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
778 integrator.setInitialStepSize(60);
779 propagator = new FieldNumericalPropagator<>(field, integrator);
780 propagator.setOrbitType(type);
781 propagator.setInitialState(initialState);
782 final FieldAbsoluteDate<T> resetDate = initDate.shiftedBy(1000);
783 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.RESET_DERIVATIVES);
784 propagator.addEventDetector(new FieldDateDetector<>(resetDate).withHandler(checking));
785 final double dt = 3200;
786 checking.assertEvent(false);
787 Assert.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
788 propagator.setResetAtEnd(true);
789 final FieldSpacecraftState<T> finalState =
790 propagator.propagate(initDate.shiftedBy(dt));
791 Assert.assertEquals(dt, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
792 checking.assertEvent(true);
793 final double n = FastMath.sqrt(initialState.getMu().getReal() / initialState.getA().getReal()) / initialState.getA().getReal();
794 Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 1.0e-10);
795 Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-10);
796 Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 1.0e-10);
797 Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-10);
798 Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 1.0e-10);
799 Assert.assertEquals(initialState.getLM().getReal() + n * dt, finalState.getLM().getReal(), 6.0e-10);
800 }
801
802 @Test
803 public void testContinueEvent() {
804 doTestContinueEvent(Decimal64Field.getInstance());
805 }
806
807 private <T extends CalculusFieldElement<T>> void doTestContinueEvent(Field<T> field) {
808 T zero = field.getZero();
809
810 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
811 FieldSpacecraftState<T> initialState;
812 FieldNumericalPropagator<T> propagator;
813 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
814 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
815
816 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
817 FramesFactory.getEME2000(), initDate, zero.add(mu));
818 initialState = new FieldSpacecraftState<>(orbit);
819 OrbitType type = OrbitType.EQUINOCTIAL;
820 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
821 AdaptiveStepsizeFieldIntegrator<T>integrator =
822 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
823 integrator.setInitialStepSize(60);
824 propagator = new FieldNumericalPropagator<>(field, integrator);
825 propagator.setOrbitType(type);
826 propagator.setInitialState(initialState);
827
828
829
830
831 final FieldAbsoluteDate<T> resetDate = initDate.shiftedBy(1000);
832 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.CONTINUE);
833 propagator.addEventDetector(new FieldDateDetector<>(resetDate).withHandler(checking));
834 final double dt = 3200;
835 checking.assertEvent(false);
836 Assert.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
837 propagator.setResetAtEnd(false);
838 final FieldSpacecraftState<T> finalState =
839 propagator.propagate(initDate.shiftedBy(dt));
840 Assert.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
841 checking.assertEvent(true);
842 final double n = FastMath.sqrt(initialState.getMu().getReal() / initialState.getA().getReal()) / initialState.getA().getReal();
843 Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 1.0e-10);
844 Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-10);
845 Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 1.0e-10);
846 Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-10);
847 Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 1.0e-10);
848 Assert.assertEquals(initialState.getLM().getReal() + n * dt, finalState.getLM().getReal(), 6.0e-10);
849 }
850
851 @Test
852 public void testAdditionalStateEvent() {
853 doTestAdditionalStateEvent(Decimal64Field.getInstance());
854 }
855
856 private <T extends CalculusFieldElement<T>> void doTestAdditionalStateEvent(Field<T> field) {
857 T zero = field.getZero();
858
859 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
860 FieldSpacecraftState<T> initialState;
861 FieldNumericalPropagator<T> propagator;
862 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
863 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
864
865 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
866 FramesFactory.getEME2000(), initDate, zero.add(mu));
867 initialState = new FieldSpacecraftState<>(orbit);
868 OrbitType type = OrbitType.EQUINOCTIAL;
869 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
870 AdaptiveStepsizeFieldIntegrator<T>integrator =
871 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
872 integrator.setInitialStepSize(60);
873 propagator = new FieldNumericalPropagator<>(field, integrator);
874 propagator.setOrbitType(type);
875 propagator.setInitialState(initialState);
876
877 propagator.addAdditionalDerivativesProvider(new FieldAdditionalDerivativesProvider<T>() {
878
879 public String getName() {
880 return "linear";
881 }
882
883 public int getDimension() {
884 return 1;
885 }
886
887 public T[] derivatives(FieldSpacecraftState<T> s) {
888 T[] pDot = MathArrays.buildArray(field, 1);
889 pDot[0] = field.getOne();
890 return pDot;
891 }
892 });
893 try {
894 propagator.addAdditionalDerivativesProvider(new FieldAdditionalDerivativesProvider<T>() {
895
896 public String getName() {
897 return "linear";
898 }
899
900 public int getDimension() {
901 return 1;
902 }
903
904 public T[] derivatives(FieldSpacecraftState<T> s) {
905 T[] pDot = MathArrays.buildArray(field, 1);
906 pDot[0] = field.getOne();
907 return pDot;
908 }
909 });
910 Assert.fail("an exception should have been thrown");
911 } catch (OrekitException oe) {
912 Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
913 }
914 try {
915 propagator.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {
916 public String getName() {
917 return "linear";
918 }
919
920 public T[] getAdditionalState(FieldSpacecraftState<T> state) {
921 return null;
922 }
923 });
924 Assert.fail("an exception should have been thrown");
925 } catch (OrekitException oe) {
926 Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
927 }
928 propagator.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {
929 public String getName() {
930 return "constant";
931 }
932
933 public T[] getAdditionalState(FieldSpacecraftState<T> state) {
934 T[] ret = MathArrays.buildArray(field, 1);
935 ret[0] = zero.add(1.0);
936 return ret;
937 }
938 });
939 Assert.assertTrue(propagator.isAdditionalStateManaged("linear"));
940 Assert.assertTrue(propagator.isAdditionalStateManaged("constant"));
941 Assert.assertFalse(propagator.isAdditionalStateManaged("non-managed"));
942 Assert.assertEquals(2, propagator.getManagedAdditionalStates().length);
943 propagator.setInitialState(propagator.getInitialState().addAdditionalState("linear", zero.add(1.5)));
944
945 CheckingHandler<AdditionalStateLinearDetector<T>, T> checking =
946 new CheckingHandler<AdditionalStateLinearDetector<T>, T>(Action.STOP);
947 propagator.addEventDetector(new AdditionalStateLinearDetector<T>(zero.add(10.0), zero.add(1.0e-8)).withHandler(checking));
948
949 final double dt = 3200;
950 checking.assertEvent(false);
951 final FieldSpacecraftState<T> finalState =
952 propagator.propagate(initDate.shiftedBy(dt));
953 checking.assertEvent(true);
954 Assert.assertEquals(3.0, finalState.getAdditionalState("linear")[0].getReal(), 1.0e-8);
955 Assert.assertEquals(1.5, finalState.getDate().durationFrom(initDate).getReal(), 1.0e-8);
956
957 }
958
959 private static class AdditionalStateLinearDetector<T extends CalculusFieldElement<T>>
960 extends FieldAbstractDetector<AdditionalStateLinearDetector<T>, T> {
961
962 public AdditionalStateLinearDetector(T maxCheck, T threshold) {
963 this(maxCheck, threshold, DEFAULT_MAX_ITER, new FieldStopOnEvent<AdditionalStateLinearDetector<T>, T>());
964 }
965
966 private AdditionalStateLinearDetector(T maxCheck, T threshold, int maxIter,
967 FieldEventHandler<? super AdditionalStateLinearDetector<T>, T> handler) {
968 super(maxCheck, threshold, maxIter, handler);
969 }
970
971 protected AdditionalStateLinearDetector<T> create(final T newMaxCheck, final T newThreshold,
972 final int newMaxIter,
973 final FieldEventHandler<? super AdditionalStateLinearDetector<T>, T> newHandler) {
974 return new AdditionalStateLinearDetector<T>(newMaxCheck, newThreshold, newMaxIter, newHandler);
975 }
976
977 public T g(FieldSpacecraftState<T> s) {
978 return s.getAdditionalState("linear")[0].subtract(3.0);
979 }
980
981 }
982
983
984 @Test
985 public void testResetAdditionalStateEvent() {
986 doTestResetAdditionalStateEvent(Decimal64Field.getInstance());
987 }
988
989 private <T extends CalculusFieldElement<T>> void doTestResetAdditionalStateEvent(final Field<T> field) {
990 FieldNumericalPropagator<T> propagator = createPropagator(field);
991
992
993 propagator.addAdditionalDerivativesProvider(new FieldAdditionalDerivativesProvider<T>() {
994
995 public String getName() {
996 return "linear";
997 }
998
999 public int getDimension() {
1000 return 1;
1001 }
1002
1003 public T[] derivatives(FieldSpacecraftState<T> s) {
1004 T[] pDot = MathArrays.buildArray(field, 1);
1005 pDot[0] = field.getOne();
1006 return pDot;
1007 }
1008 });
1009 propagator.setInitialState(propagator.getInitialState().addAdditionalState("linear",
1010 field.getZero().add(1.5)));
1011
1012 CheckingHandler<AdditionalStateLinearDetector<T>, T> checking =
1013 new CheckingHandler<AdditionalStateLinearDetector<T>, T>(Action.RESET_STATE) {
1014 public FieldSpacecraftState<T> resetState(AdditionalStateLinearDetector<T> detector, FieldSpacecraftState<T> oldState)
1015 {
1016 return oldState.addAdditionalState("linear", oldState.getAdditionalState("linear")[0].multiply(2));
1017 }
1018 };
1019
1020 propagator.addEventDetector(new AdditionalStateLinearDetector<T>(field.getZero().add(10.0),
1021 field.getZero().add(1.0e-8)).withHandler(checking));
1022
1023 final double dt = 3200;
1024 checking.assertEvent(false);
1025 final FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
1026 final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(dt));
1027
1028 Assert.assertEquals(dt + 4.5, finalState.getAdditionalState("linear")[0].getReal(), 1.0e-8);
1029 Assert.assertEquals(dt, finalState.getDate().durationFrom(initDate).getReal(), 1.0e-8);
1030
1031 }
1032
1033 @Test
1034 public void testEventDetectionBug() {
1035 doTestEventDetectionBug(Decimal64Field.getInstance());
1036 }
1037
1038 private <T extends CalculusFieldElement<T>> void doTestEventDetectionBug(final Field<T> field) {
1039
1040 T zero = field.getZero();
1041 TimeScale utc = TimeScalesFactory.getUTC();
1042 FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, 2005, 1, 1, 0, 0, 0.0, utc);
1043 T duration = zero.add(100000.0);
1044 FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(initialDate, duration);
1045
1046
1047 Frame EME2000 = FramesFactory.getEME2000();
1048
1049
1050
1051 double a = 35786000. + 6378137.0;
1052 double e = 0.70;
1053 double rApogee = a*(1+e);
1054 double vApogee = FastMath.sqrt(mu*(1-e)/(a*(1+e)));
1055 FieldOrbit<T> geo = new FieldCartesianOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(zero.add(rApogee), zero, zero),
1056 new FieldVector3D<>(zero, zero.add(vApogee), zero)),
1057 EME2000, initialDate, zero.add(mu));
1058
1059
1060 duration = geo.getKeplerianPeriod();
1061 endDate = new FieldAbsoluteDate<>(initialDate, duration);
1062
1063
1064 final double minStep = 0.001;
1065 final double maxStep = 1000;
1066 final double initStep = 60;
1067 final OrbitType type = OrbitType.EQUINOCTIAL;
1068 final double[] absTolerance = {
1069 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001};
1070 final double[] relTolerance = {
1071 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7};
1072
1073 AdaptiveStepsizeFieldIntegrator<T> integrator =
1074 new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, absTolerance, relTolerance);
1075 integrator.setInitialStepSize(initStep);
1076
1077
1078 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
1079 propagator.setOrbitType(type);
1080 T mass = field.getZero().add(1000.0);
1081 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(geo, mass);
1082 propagator.setInitialState(initialState);
1083 propagator.setOrbitType(OrbitType.CARTESIAN);
1084
1085
1086
1087 FieldApsideDetector<T> event1 = new FieldApsideDetector<>(geo);
1088 propagator.addEventDetector(event1);
1089
1090
1091 propagator.clearStepHandlers();
1092
1093
1094 FieldSpacecraftState<T> finalState = propagator.propagate(endDate);
1095
1096
1097 Assert.assertTrue(endDate.durationFrom(finalState.getDate()).getReal() > 40000.0);
1098 }
1099
1100 @Test
1101 public void testEphemerisGenerationIssue14() {
1102 doTestEphemerisGenerationIssue14(Decimal64Field.getInstance());
1103 }
1104
1105 private <T extends CalculusFieldElement<T>> void doTestEphemerisGenerationIssue14(Field<T> field)
1106 {
1107
1108
1109 FieldNumericalPropagator<T> propagator = createPropagator(field);
1110 final double dt = 3200;
1111 final FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
1112
1113 propagator.setOrbitType(OrbitType.CARTESIAN);
1114 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
1115 propagator.propagate(initDate.shiftedBy(dt));
1116 final FieldBoundedPropagator<T> ephemeris1 = generator.getGeneratedEphemeris();
1117 Assert.assertEquals(initDate, ephemeris1.getMinDate());
1118 Assert.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
1119
1120 propagator.getPVCoordinates(initDate.shiftedBy( 2 * dt), FramesFactory.getEME2000());
1121 propagator.getPVCoordinates(initDate.shiftedBy(-2 * dt), FramesFactory.getEME2000());
1122
1123
1124 Assert.assertEquals(initDate, ephemeris1.getMinDate());
1125 Assert.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
1126
1127 final FieldBoundedPropagator<T> ephemeris2 = generator.getGeneratedEphemeris();
1128 Assert.assertEquals(initDate.shiftedBy(-2 * dt), ephemeris2.getMinDate());
1129 Assert.assertEquals(initDate.shiftedBy( 2 * dt), ephemeris2.getMaxDate());
1130
1131
1132 Assert.assertEquals(initDate, ephemeris1.getMinDate());
1133 Assert.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
1134
1135 }
1136
1137 @Test
1138 public void testEphemerisAdditionalState() {
1139 doTestEphemerisAdditionalState(Decimal64Field.getInstance());
1140 }
1141
1142 private <T extends CalculusFieldElement<T>> void doTestEphemerisAdditionalState(final Field<T> field) {
1143
1144
1145 final double dt = -3200;
1146 final double rate = 2.0;
1147 FieldNumericalPropagator<T> propagator = createPropagator(field);
1148 FieldAbsoluteDate<T> initDate = propagator.getInitialState().getDate();
1149
1150 propagator.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {
1151 public String getName() {
1152 return "squaredA";
1153 }
1154 public T[] getAdditionalState(FieldSpacecraftState<T> state) {
1155 T[] a = MathArrays.buildArray(field, 1);
1156 a[0] = state.getA().multiply(state.getA());
1157 return a;
1158 }
1159 });
1160 propagator.addAdditionalDerivativesProvider(new FieldAdditionalDerivativesProvider<T>() {
1161 public String getName() {
1162 return "extra";
1163 }
1164 public int getDimension() {
1165 return 1;
1166 }
1167 public T[] derivatives(FieldSpacecraftState<T> s) {
1168 T[] pDot = MathArrays.buildArray(field, 1);
1169 pDot[0] = field.getZero().newInstance(rate);
1170 return pDot;
1171 }
1172 });
1173 propagator.setInitialState(propagator.getInitialState().addAdditionalState("extra", field.getZero().add(1.5)));
1174
1175 propagator.setOrbitType(OrbitType.CARTESIAN);
1176 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
1177 propagator.propagate(initDate.shiftedBy(dt));
1178 final FieldBoundedPropagator<T> ephemeris1 = generator.getGeneratedEphemeris();
1179 Assert.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMinDate());
1180 Assert.assertEquals(initDate, ephemeris1.getMaxDate());
1181 try {
1182 ephemeris1.propagate(ephemeris1.getMinDate().shiftedBy(-10.0));
1183 Assert.fail("an exception should have been thrown");
1184 } catch (OrekitException pe) {
1185 Assert.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE, pe.getSpecifier());
1186 }
1187 try {
1188 ephemeris1.propagate(ephemeris1.getMaxDate().shiftedBy(+10.0));
1189 Assert.fail("an exception should have been thrown");
1190 } catch (OrekitException pe) {
1191 Assert.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER, pe.getSpecifier());
1192 }
1193
1194 double shift = -60;
1195 FieldSpacecraftState<T> s = ephemeris1.propagate(initDate.shiftedBy(shift));
1196 Assert.assertEquals(2, s.getAdditionalStatesValues().size());
1197 Assert.assertTrue(s.hasAdditionalState("squaredA"));
1198 Assert.assertTrue(s.hasAdditionalState("extra"));
1199 Assert.assertEquals(s.getA().multiply(s.getA()).getReal(), s.getAdditionalState("squaredA")[0].getReal(), 1.0e-10);
1200 Assert.assertEquals(1.5 + shift * rate, s.getAdditionalState("extra")[0].getReal(), 1.0e-10);
1201
1202 try {
1203 ephemeris1.resetInitialState(s);
1204 Assert.fail("an exception should have been thrown");
1205 } catch (OrekitException oe) {
1206 Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
1207 }
1208
1209 }
1210
1211 @Test
1212 public void testIssue157() {
1213 doTestIssue157(Decimal64Field.getInstance());
1214 }
1215
1216 private <T extends CalculusFieldElement<T>> void doTestIssue157(final Field<T> field) {
1217 try {
1218 FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(field.getZero().add(13378000),
1219 field.getZero().add(0.05),
1220 field.getZero().add(0),
1221 field.getZero().add(0),
1222 field.getZero().add(FastMath.PI),
1223 field.getZero().add(0),
1224 PositionAngle.MEAN,
1225 FramesFactory.getTOD(false),
1226 new FieldAbsoluteDate<>(field, 2003, 5, 6, TimeScalesFactory.getUTC()),
1227 field.getZero().add(Constants.EIGEN5C_EARTH_MU));
1228 FieldNumericalPropagator.tolerances(field.getZero().add(1.0), orbit, OrbitType.KEPLERIAN);
1229 Assert.fail("an exception should have been thrown");
1230 } catch (OrekitException oe) {
1231 Assert.assertEquals(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, oe.getSpecifier());
1232 }
1233 }
1234
1235 @Test
1236 public void testIssue704() {
1237 doTestIssue704(Decimal64Field.getInstance());
1238 }
1239
1240 private <T extends CalculusFieldElement<T>> void doTestIssue704(final Field<T> field) {
1241
1242 T zero = field.getZero();
1243
1244
1245 final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
1246 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
1247 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
1248 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
1249 FramesFactory.getEME2000(), initDate, zero.add(mu));
1250 final FieldPVCoordinates<T> pv = orbit.getPVCoordinates();
1251
1252
1253 final T dP = zero.add(10.0);
1254
1255
1256 final T r2 = pv.getPosition().getNormSq();
1257 final T v = pv.getVelocity().getNorm();
1258 final T dV = dP.multiply(orbit.getMu()).divide(v.multiply(r2));
1259
1260
1261 final double[][] tolCart1 = FieldNumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN);
1262 final double[][] tolCart2 = FieldNumericalPropagator.tolerances(dP, dV, orbit, OrbitType.CARTESIAN);
1263 for (int i = 0; i < tolCart1.length; i++) {
1264 Assert.assertArrayEquals(tolCart1[i], tolCart2[i], Double.MIN_VALUE);
1265 }
1266
1267
1268 final double[][] tolKep1 = FieldNumericalPropagator.tolerances(dP, orbit, OrbitType.KEPLERIAN);
1269 final double[][] tolKep2 = FieldNumericalPropagator.tolerances(dP, dV, orbit, OrbitType.KEPLERIAN);
1270 for (int i = 0; i < tolCart1.length; i++) {
1271 Assert.assertArrayEquals(tolKep1[i], tolKep2[i], Double.MIN_VALUE);
1272 }
1273
1274 }
1275
1276 @Test
1277 public void testShiftKeplerianEllipticTrueWithoutDerivatives() {
1278 doTestShiftKeplerianEllipticTrueWithoutDerivatives(Decimal64Field.getInstance());
1279 }
1280
1281 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticTrueWithoutDerivatives(final Field<T> field) {
1282 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.TRUE, false,
1283 18.1, 72.0, 437.3, 1601.1, 3141.8);
1284 }
1285
1286 @Test
1287 public void testShiftKeplerianEllipticTrueWithDerivatives() {
1288 doTestShiftKeplerianEllipticTrueWithDerivatives(Decimal64Field.getInstance());
1289 }
1290
1291 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticTrueWithDerivatives(final Field<T> field) {
1292 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
1293 1.14, 9.1, 140.3, 1066.7, 3306.9);
1294 }
1295
1296 @Test
1297 public void testShiftKeplerianEllipticEccentricWithoutDerivatives() {
1298 doTestShiftKeplerianEllipticEccentricWithoutDerivatives(Decimal64Field.getInstance());
1299 }
1300
1301 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticEccentricWithoutDerivatives(final Field<T> field) {
1302 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, false,
1303 18.1, 72.0, 437.3, 1601.1, 3141.8);
1304 }
1305
1306 @Test
1307 public void testShiftKeplerianEllipticEcentricWithDerivatives() {
1308 doTestShiftKeplerianEllipticEcentricWithDerivatives(Decimal64Field.getInstance());
1309 }
1310
1311 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticEcentricWithDerivatives(final Field<T> field) {
1312 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, true,
1313 1.14, 9.1, 140.3, 1066.7, 3306.9);
1314 }
1315
1316 @Test
1317 public void testShiftKeplerianEllipticMeanWithoutDerivatives() {
1318 doTestShiftKeplerianEllipticMeanWithoutDerivatives(Decimal64Field.getInstance());
1319 }
1320
1321 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticMeanWithoutDerivatives(final Field<T> field) {
1322 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.MEAN, false,
1323 18.1, 72.0, 437.3, 1601.1, 3141.8);
1324 }
1325
1326 @Test
1327 public void testShiftKeplerianEllipticMeanWithDerivatives() {
1328 doTestShiftKeplerianEllipticMeanWithDerivatives(Decimal64Field.getInstance());
1329 }
1330
1331 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianEllipticMeanWithDerivatives(final Field<T> field) {
1332 doTestShift(createEllipticOrbit(field), OrbitType.KEPLERIAN, PositionAngle.MEAN, true,
1333 1.14, 9.1, 140.3, 1066.7, 3306.9);
1334 }
1335
1336 @Test
1337 public void testShiftKeplerianHyperbolicTrueWithoutDerivatives() {
1338 doTestShiftKeplerianHyperbolicTrueWithoutDerivatives(Decimal64Field.getInstance());
1339 }
1340
1341 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicTrueWithoutDerivatives(final Field<T> field) {
1342 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.TRUE, false,
1343 0.484, 1.94, 12.1, 48.3, 108.5);
1344 }
1345
1346 @Test
1347 public void testShiftKeplerianHyperbolicTrueWithDerivatives() {
1348 doTestShiftKeplerianHyperbolicTrueWithDerivatives(Decimal64Field.getInstance());
1349 }
1350
1351 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicTrueWithDerivatives(final Field<T> field) {
1352 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
1353 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1354 }
1355
1356 @Test
1357 public void testShiftKeplerianHyperbolicEccentricWithoutDerivatives() {
1358 doTestShiftKeplerianHyperbolicEccentricWithoutDerivatives(Decimal64Field.getInstance());
1359 }
1360
1361 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicEccentricWithoutDerivatives(final Field<T> field) {
1362 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, false,
1363 0.484, 1.94, 12.1, 48.3, 108.5);
1364 }
1365
1366 @Test
1367 public void testShiftKeplerianHyperbolicEcentricWithDerivatives() {
1368 doTestShiftKeplerianHyperbolicEcentricWithDerivatives(Decimal64Field.getInstance());
1369 }
1370
1371 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicEcentricWithDerivatives(final Field<T> field) {
1372 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, true,
1373 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1374 }
1375
1376 @Test
1377 public void testShiftKeplerianHyperbolicMeanWithoutDerivatives() {
1378 doTestShiftKeplerianHyperbolicMeanWithoutDerivatives(Decimal64Field.getInstance());
1379 }
1380
1381 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicMeanWithoutDerivatives(final Field<T> field) {
1382 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.MEAN, false,
1383 0.484, 1.94, 12.1, 48.3, 108.5);
1384 }
1385
1386 @Test
1387 public void testShiftKeplerianHyperbolicMeanWithDerivatives() {
1388 doTestShiftKeplerianHyperbolicMeanWithDerivatives(Decimal64Field.getInstance());
1389 }
1390
1391 private <T extends CalculusFieldElement<T>> void doTestShiftKeplerianHyperbolicMeanWithDerivatives(final Field<T> field) {
1392 doTestShift(createHyperbolicOrbit(field), OrbitType.KEPLERIAN, PositionAngle.MEAN, true,
1393 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1394 }
1395
1396 @Test
1397 public void testShiftCartesianEllipticTrueWithoutDerivatives() {
1398 doTestShiftCartesianEllipticTrueWithoutDerivatives(Decimal64Field.getInstance());
1399 }
1400
1401 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticTrueWithoutDerivatives(final Field<T> field) {
1402 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.TRUE, false,
1403 18.1, 72.0, 437.3, 1601.1, 3141.8);
1404 }
1405
1406 @Test
1407 public void testShiftCartesianEllipticTrueWithDerivatives() {
1408 doTestShiftCartesianEllipticTrueWithDerivatives(Decimal64Field.getInstance());
1409 }
1410
1411 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticTrueWithDerivatives(final Field<T> field) {
1412 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.TRUE, true,
1413 1.14, 9.1, 140.3, 1066.7, 3306.9);
1414 }
1415
1416 @Test
1417 public void testShiftCartesianEllipticEccentricWithoutDerivatives() {
1418 doTestShiftCartesianEllipticEccentricWithoutDerivatives(Decimal64Field.getInstance());
1419 }
1420
1421 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticEccentricWithoutDerivatives(final Field<T> field) {
1422 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, false,
1423 18.1, 72.0, 437.3, 1601.1, 3141.8);
1424 }
1425
1426 @Test
1427 public void testShiftCartesianEllipticEcentricWithDerivatives() {
1428 doTestShiftCartesianEllipticEcentricWithDerivatives(Decimal64Field.getInstance());
1429 }
1430
1431 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticEcentricWithDerivatives(final Field<T> field) {
1432 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, true,
1433 1.14, 9.1, 140.3, 1066.7, 3306.9);
1434 }
1435
1436 @Test
1437 public void testShiftCartesianEllipticMeanWithoutDerivatives() {
1438 doTestShiftCartesianEllipticMeanWithoutDerivatives(Decimal64Field.getInstance());
1439 }
1440
1441 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticMeanWithoutDerivatives(final Field<T> field) {
1442 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.MEAN, false,
1443 18.1, 72.0, 437.3, 1601.1, 3141.8);
1444 }
1445
1446 @Test
1447 public void testShiftCartesianEllipticMeanWithDerivatives() {
1448 doTestShiftCartesianEllipticMeanWithDerivatives(Decimal64Field.getInstance());
1449 }
1450
1451 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianEllipticMeanWithDerivatives(final Field<T> field) {
1452 doTestShift(createEllipticOrbit(field), OrbitType.CARTESIAN, PositionAngle.MEAN, true,
1453 1.14, 9.1, 140.3, 1066.7, 3306.9);
1454 }
1455
1456 @Test
1457 public void testShiftCartesianHyperbolicTrueWithoutDerivatives() {
1458 doTestShiftCartesianHyperbolicTrueWithoutDerivatives(Decimal64Field.getInstance());
1459 }
1460
1461 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicTrueWithoutDerivatives(final Field<T> field) {
1462 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.TRUE, false,
1463 0.48, 1.93, 12.1, 48.3, 108.5);
1464 }
1465
1466 @Test
1467 public void testShiftCartesianHyperbolicTrueWithDerivatives() {
1468 doTestShiftCartesianHyperbolicTrueWithDerivatives(Decimal64Field.getInstance());
1469 }
1470
1471 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicTrueWithDerivatives(final Field<T> field) {
1472 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.TRUE, true,
1473 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1474 }
1475
1476 @Test
1477 public void testShiftCartesianHyperbolicEccentricWithoutDerivatives() {
1478 doTestShiftCartesianHyperbolicEccentricWithoutDerivatives(Decimal64Field.getInstance());
1479 }
1480
1481 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicEccentricWithoutDerivatives(final Field<T> field) {
1482 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, false,
1483 0.48, 1.93, 12.1, 48.3, 108.5);
1484 }
1485
1486 @Test
1487 public void testShiftCartesianHyperbolicEcentricWithDerivatives() {
1488 doTestShiftCartesianHyperbolicEcentricWithDerivatives(Decimal64Field.getInstance());
1489 }
1490
1491 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicEcentricWithDerivatives(final Field<T> field) {
1492 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, true,
1493 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1494 }
1495
1496 @Test
1497 public void testShiftCartesianHyperbolicMeanWithoutDerivatives() {
1498 doTestShiftCartesianHyperbolicMeanWithoutDerivatives(Decimal64Field.getInstance());
1499 }
1500
1501 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicMeanWithoutDerivatives(final Field<T> field) {
1502 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.MEAN, false,
1503 0.48, 1.93, 12.1, 48.3, 108.5);
1504 }
1505
1506 @Test
1507 public void testShiftCartesianHyperbolicMeanWithDerivatives() {
1508 doTestShiftCartesianHyperbolicMeanWithDerivatives(Decimal64Field.getInstance());
1509 }
1510
1511 private <T extends CalculusFieldElement<T>> void doTestShiftCartesianHyperbolicMeanWithDerivatives(final Field<T> field) {
1512 doTestShift(createHyperbolicOrbit(field), OrbitType.CARTESIAN, PositionAngle.MEAN, true,
1513 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1514 }
1515
1516 @Test
1517 public void testShiftCircularTrueWithoutDerivatives() {
1518 doTestShiftCircularTrueWithoutDerivatives(Decimal64Field.getInstance());
1519 }
1520
1521 private <T extends CalculusFieldElement<T>> void doTestShiftCircularTrueWithoutDerivatives(final Field<T> field) {
1522 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.TRUE, false,
1523 18.1, 72.0, 437.3, 1601.1, 3141.8);
1524 }
1525
1526 @Test
1527 public void testShiftCircularTrueWithDerivatives() {
1528 doTestShiftCircularTrueWithDerivatives(Decimal64Field.getInstance());
1529 }
1530
1531 private <T extends CalculusFieldElement<T>> void doTestShiftCircularTrueWithDerivatives(final Field<T> field) {
1532 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.TRUE, true,
1533 1.14, 9.1, 140.3, 1066.7, 3306.9);
1534 }
1535
1536 @Test
1537 public void testShiftCircularEccentricWithoutDerivatives() {
1538 doTestShiftCircularEccentricWithoutDerivatives(Decimal64Field.getInstance());
1539 }
1540
1541 private <T extends CalculusFieldElement<T>> void doTestShiftCircularEccentricWithoutDerivatives(final Field<T> field) {
1542 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.ECCENTRIC, false,
1543 18.1, 72.0, 437.3, 1601.1, 3141.8);
1544 }
1545
1546 @Test
1547 public void testShiftCircularEcentricWithDerivatives() {
1548 doTestShiftCircularEcentricWithDerivatives(Decimal64Field.getInstance());
1549 }
1550
1551 private <T extends CalculusFieldElement<T>> void doTestShiftCircularEcentricWithDerivatives(final Field<T> field) {
1552 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.ECCENTRIC, true,
1553 1.14, 9.1, 140.3, 1066.7, 3306.9);
1554 }
1555
1556 @Test
1557 public void testShiftCircularMeanWithoutDerivatives() {
1558 doTestShiftCircularMeanWithoutDerivatives(Decimal64Field.getInstance());
1559 }
1560
1561 private <T extends CalculusFieldElement<T>> void doTestShiftCircularMeanWithoutDerivatives(final Field<T> field) {
1562 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.MEAN, false,
1563 18.1, 72.0, 437.3, 1601.1, 3141.8);
1564 }
1565
1566 @Test
1567 public void testShiftCircularMeanWithDerivatives() {
1568 doTestShiftCircularMeanWithDerivatives(Decimal64Field.getInstance());
1569 }
1570
1571 private <T extends CalculusFieldElement<T>> void doTestShiftCircularMeanWithDerivatives(final Field<T> field) {
1572 doTestShift(createEllipticOrbit(field), OrbitType.CIRCULAR, PositionAngle.MEAN, true,
1573 1.14, 9.1, 140.3, 1066.7, 3306.9);
1574 }
1575
1576 @Test
1577 public void testShiftEquinoctialTrueWithoutDerivatives() {
1578 doTestShiftEquinoctialTrueWithoutDerivatives(Decimal64Field.getInstance());
1579 }
1580
1581 private <T extends CalculusFieldElement<T>> void doTestShiftEquinoctialTrueWithoutDerivatives(final Field<T> field) {
1582 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false,
1583 18.1, 72.0, 437.3, 1601.1, 3141.8);
1584 }
1585
1586 @Test
1587 public void testShiftEquinoctialTrueWithDerivatives() {
1588 doTestShiftEquinoctialTrueWithDerivatives(Decimal64Field.getInstance());
1589 }
1590
1591 private <T extends CalculusFieldElement<T>> void doTestShiftEquinoctialTrueWithDerivatives(final Field<T> field) {
1592 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.TRUE, true,
1593 1.14, 9.1, 140.3, 1066.7, 3306.9);
1594 }
1595
1596 @Test
1597 public void testShiftEquinoctialEccentricWithoutDerivatives() {
1598 doTestShiftEquinoctialEccentricWithoutDerivatives(Decimal64Field.getInstance());
1599 }
1600
1601 private <T extends CalculusFieldElement<T>> void doTestShiftEquinoctialEccentricWithoutDerivatives(final Field<T> field) {
1602 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC, false,
1603 18.1, 72.0, 437.3, 1601.1, 3141.8);
1604 }
1605
1606 @Test
1607 public void testShiftEquinoctialEcentricWithDerivatives() {
1608 doTtestShiftEquinoctialEcentricWithDerivatives(Decimal64Field.getInstance());
1609 }
1610
1611 private <T extends CalculusFieldElement<T>> void doTtestShiftEquinoctialEcentricWithDerivatives(final Field<T> field) {
1612 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC, true,
1613 1.14, 9.1, 140.3, 1066.7, 3306.9);
1614 }
1615
1616 @Test
1617 public void testShiftEquinoctialMeanWithoutDerivatives() {
1618 doTestShiftEquinoctialMeanWithoutDerivatives(Decimal64Field.getInstance());
1619 }
1620
1621 private <T extends CalculusFieldElement<T>> void doTestShiftEquinoctialMeanWithoutDerivatives(final Field<T> field) {
1622 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.MEAN, false,
1623 18.1, 72.0, 437.3, 1601.1, 3141.8);
1624 }
1625
1626 @Test
1627 public void testShiftEquinoctialMeanWithDerivatives() {
1628 doTestShiftEquinoctialMeanWithDerivatives(Decimal64Field.getInstance());
1629 }
1630
1631 private <T extends CalculusFieldElement<T>> void doTestShiftEquinoctialMeanWithDerivatives(final Field<T> field) {
1632 doTestShift(createEllipticOrbit(field), OrbitType.EQUINOCTIAL, PositionAngle.MEAN, true,
1633 1.14, 9.1, 140.3, 1066.7, 3306.9);
1634 }
1635
1636 private static <T extends CalculusFieldElement<T>> void doTestShift(final FieldCartesianOrbit<T> orbit, final OrbitType orbitType,
1637 final PositionAngle angleType, final boolean withDerivatives,
1638 final double error60s, final double error120s,
1639 final double error300s, final double error600s,
1640 final double error900s)
1641 {
1642
1643 T zero = orbit.getDate().getField().getZero();
1644
1645 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1646 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1647 final FieldNumericalPropagator<T> np = createPropagator(new FieldSpacecraftState<>(orbit), orbitType, angleType);
1648
1649
1650
1651 final FieldAbsoluteDate<T> reference = orbit.getDate().shiftedBy(60.0);
1652 final ShiftChecker<T> checker = new ShiftChecker<>(withDerivatives, orbitType, angleType,
1653 error60s,
1654 error120s, error300s,
1655 error600s, error900s);
1656 @SuppressWarnings("unchecked")
1657 FieldTimeStamped<T>[] dates = (FieldTimeStamped<T>[]) Array.newInstance(FieldTimeStamped.class, 6);
1658 dates[0] = reference;
1659 dates[1] = reference.shiftedBy( 60.0);
1660 dates[2] = reference.shiftedBy(120.0);
1661 dates[3] = reference.shiftedBy(300.0);
1662 dates[4] = reference.shiftedBy(600.0);
1663 dates[5] = reference.shiftedBy(900.0);
1664 np.addEventDetector(new FieldDateDetector<T>(zero.add(30.0), zero.add(1.0e-9), (FieldTimeStamped<T>[]) dates).
1665 withHandler(checker));
1666 np.propagate(reference.shiftedBy(1000.0));
1667 }
1668
1669 private static class ShiftChecker<T extends CalculusFieldElement<T>> implements FieldEventHandler<FieldDateDetector<T>, T> {
1670
1671 private final boolean withDerivatives;
1672 private final OrbitType orbitType;
1673 private final PositionAngle angleType;
1674 private final double error60s;
1675 private final double error120s;
1676 private final double error300s;
1677 private final double error600s;
1678 private final double error900s;
1679 private FieldSpacecraftState<T> referenceState;
1680
1681 ShiftChecker(final boolean withDerivatives, final OrbitType orbitType,
1682 final PositionAngle angleType, final double error60s,
1683 final double error120s, final double error300s,
1684 final double error600s, final double error900s) {
1685 this.withDerivatives = withDerivatives;
1686 this.orbitType = orbitType;
1687 this.angleType = angleType;
1688 this.error60s = error60s;
1689 this.error120s = error120s;
1690 this.error300s = error300s;
1691 this.error600s = error600s;
1692 this.error900s = error900s;
1693 this.referenceState = null;
1694 }
1695
1696 @Override
1697 public Action eventOccurred(final FieldSpacecraftState<T> s, final FieldDateDetector<T> detector,
1698 final boolean increasing)
1699 {
1700 if (referenceState == null) {
1701
1702 if (withDerivatives) {
1703 referenceState = s;
1704 } else {
1705
1706 final T[] stateVector = MathArrays.buildArray(s.getDate().getField(), 6);
1707 final FieldOrbit<T> o = s.getOrbit();
1708 orbitType.mapOrbitToArray(o, angleType, stateVector, null);
1709 final FieldOrbit<T> fixedOrbit = orbitType.mapArrayToOrbit(stateVector, null, angleType,
1710 o.getDate(), o.getMu(), o.getFrame());
1711 referenceState = new FieldSpacecraftState<>(fixedOrbit, s.getAttitude(), s.getMass());
1712 }
1713 } else {
1714
1715 final T dt = s.getDate().durationFrom(referenceState.getDate());
1716 final FieldSpacecraftState<T> shifted = referenceState.shiftedBy(dt);
1717 final T error = FieldVector3D.distance(shifted.getPVCoordinates().getPosition(),
1718 s.getPVCoordinates().getPosition());
1719 switch ((int) FastMath.rint(dt.getReal())) {
1720 case 60 :
1721 Assert.assertEquals(error60s, error.getReal(), 0.01 * error60s);
1722 break;
1723 case 120 :
1724 Assert.assertEquals(error120s, error.getReal(), 0.01 * error120s);
1725 break;
1726 case 300 :
1727 Assert.assertEquals(error300s, error.getReal(), 0.01 * error300s);
1728 break;
1729 case 600 :
1730 Assert.assertEquals(error600s, error.getReal(), 0.01 * error600s);
1731 break;
1732 case 900 :
1733 Assert.assertEquals(error900s, error.getReal(), 0.01 * error900s);
1734 break;
1735 default :
1736
1737 Assert.fail("no error set for dt = " + dt);
1738 break;
1739 }
1740 }
1741 return Action.CONTINUE;
1742 }
1743
1744 }
1745
1746 private static <T extends CalculusFieldElement<T>> FieldNumericalPropagator<T> createPropagator(FieldSpacecraftState<T> spacecraftState,
1747 OrbitType orbitType,
1748 PositionAngle angleType)
1749 {
1750
1751 final Field<T> field = spacecraftState.getDate().getField();
1752 final T zero = field.getZero();
1753 final double minStep = 0.001;
1754 final double maxStep = 120.0;
1755 final T positionTolerance = zero.add(0.1);
1756 final int degree = 20;
1757 final int order = 20;
1758 final double spacecraftArea = 1.0;
1759 final double spacecraftDragCoefficient = 2.0;
1760 final double spacecraftReflectionCoefficient = 2.0;
1761
1762
1763 final double[][] tol = FieldNumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
1764 final FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, tol[0], tol[1]);
1765 final FieldNumericalPropagator<T> np = new FieldNumericalPropagator<>(field, integrator);
1766 np.setOrbitType(orbitType);
1767 np.setPositionAngleType(angleType);
1768 np.setInitialState(spacecraftState);
1769
1770
1771 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1772 Constants.WGS84_EARTH_FLATTENING,
1773 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1774 final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
1775 np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
1776
1777
1778 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
1779 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
1780
1781
1782 MarshallSolarActivityFutureEstimation msafe =
1783 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
1784 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
1785 DataContext.getDefault().getDataProvidersManager().feed(msafe.getSupportedNames(), msafe);
1786 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
1787 np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
1788
1789
1790 np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(),
1791 earth.getEquatorialRadius(),
1792 new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
1793
1794 return np;
1795
1796 }
1797
1798 private static <T extends CalculusFieldElement<T>> FieldCartesianOrbit<T> createEllipticOrbit(Field<T> field) {
1799 T zero = field.getZero();
1800 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1801 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(6896874.444705),
1802 zero.add(1956581.072644),
1803 zero.add(-147476.245054));
1804 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(169.816407662),
1805 zero.add(-1126.783301861),
1806 zero.add(-7332.745712770));
1807 final TimeStampedFieldPVCoordinates<T> pv = new TimeStampedFieldPVCoordinates<>(date, position, velocity, FieldVector3D.getZero(field));
1808 final Frame frame = FramesFactory.getEME2000();
1809 final double mu = Constants.EIGEN5C_EARTH_MU;
1810 return new FieldCartesianOrbit<>(pv, frame, zero.add(mu));
1811 }
1812
1813 private static <T extends CalculusFieldElement<T>> FieldCartesianOrbit<T> createHyperbolicOrbit(Field<T> field) {
1814 T zero = field.getZero();
1815 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1816 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(224267911.905821),
1817 zero.add(290251613.109399),
1818 zero.add(45534292.777492));
1819 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-1494.068165293),
1820 zero.add(1124.771027677),
1821 zero.add(526.915286134));
1822 final TimeStampedFieldPVCoordinates<T> pv = new TimeStampedFieldPVCoordinates<>(date, position, velocity, FieldVector3D.getZero(field));
1823 final Frame frame = FramesFactory.getEME2000();
1824 final double mu = Constants.EIGEN5C_EARTH_MU;
1825 return new FieldCartesianOrbit<>(pv, frame, zero.add(mu));
1826 }
1827
1828 private class CheckingHandler<D extends FieldEventDetector<T>, T extends CalculusFieldElement<T>> implements FieldEventHandler<D, T> {
1829
1830 private final Action actionOnEvent;
1831 private boolean gotHere;
1832
1833 public CheckingHandler(final Action actionOnEvent) {
1834 this.actionOnEvent = actionOnEvent;
1835 this.gotHere = false;
1836 }
1837
1838 public void assertEvent(boolean expected) {
1839 Assert.assertEquals(expected, gotHere);
1840 }
1841
1842 public Action eventOccurred(FieldSpacecraftState<T> s, D detector, boolean increasing) {
1843 gotHere = true;
1844 return actionOnEvent;
1845 }
1846
1847 }
1848
1849 private <T extends CalculusFieldElement<T>> FieldNumericalPropagator<T> createPropagator(Field<T> field)
1850 {
1851 T zero = field.getZero();
1852 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6),
1853 zero.add(1.0e6),
1854 zero.add(4.0e6));
1855 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0),
1856 zero.add(8000.0),
1857 zero.add(1000.0));
1858 FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
1859 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
1860 FramesFactory.getEME2000(), initDate, zero.add(mu));
1861 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
1862 OrbitType type = OrbitType.EQUINOCTIAL;
1863 double[][] tolerance = FieldNumericalPropagator.tolerances(zero.add(0.001), orbit, type);
1864 AdaptiveStepsizeFieldIntegrator<T> integrator =
1865 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
1866
1867 integrator.setInitialStepSize(60);
1868 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
1869 propagator.setOrbitType(type);
1870
1871
1872 propagator.setInitialState(initialState);
1873 return propagator;
1874 }
1875
1876 private <T extends CalculusFieldElement<T>> FieldTimeStamped<T>[] toArray(final FieldAbsoluteDate<T> date) {
1877 @SuppressWarnings("unchecked")
1878 final FieldTimeStamped<T>[] array = (FieldTimeStamped<T>[]) Array.newInstance(FieldTimeStamped.class, 1);
1879 array[0] = date;
1880 return array;
1881 }
1882
1883 @Before
1884 public void setUp() {
1885 Utils.setDataRoot("regular-data:potential/shm-format");
1886 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
1887 mu = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
1888 }
1889
1890 }
1891