1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.analytical;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
21 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
22 import org.hipparchus.util.FastMath;
23 import org.junit.jupiter.api.Assertions;
24 import org.junit.jupiter.api.BeforeEach;
25 import org.junit.jupiter.api.Test;
26 import org.orekit.TestUtils;
27 import org.orekit.Utils;
28 import org.orekit.attitudes.AttitudeProvider;
29 import org.orekit.attitudes.LofOffset;
30 import org.orekit.bodies.CelestialBodyFactory;
31 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
32 import org.orekit.forces.gravity.ThirdBodyAttraction;
33 import org.orekit.forces.gravity.potential.GravityFieldFactory;
34 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
35 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
36 import org.orekit.forces.maneuvers.ConstantThrustManeuver;
37 import org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel;
38 import org.orekit.frames.FramesFactory;
39 import org.orekit.frames.LOFType;
40 import org.orekit.orbits.CircularOrbit;
41 import org.orekit.orbits.KeplerianOrbit;
42 import org.orekit.orbits.Orbit;
43 import org.orekit.orbits.PositionAngleType;
44 import org.orekit.propagation.*;
45 import org.orekit.propagation.numerical.NumericalPropagator;
46 import org.orekit.time.AbsoluteDate;
47 import org.orekit.time.DateComponents;
48 import org.orekit.time.TimeComponents;
49 import org.orekit.time.TimeScalesFactory;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.PVCoordinates;
52
53 import java.io.IOException;
54 import java.text.ParseException;
55 import java.util.ArrayList;
56 import java.util.List;
57
58 public class AdapterPropagatorTest {
59
60 @Test
61 public void testLowEarthOrbit() throws ParseException, IOException {
62
63 Orbit leo = new CircularOrbit(7200000.0, -1.0e-5, 2.0e-4,
64 FastMath.toRadians(98.0),
65 FastMath.toRadians(123.456),
66 0.0, PositionAngleType.MEAN,
67 FramesFactory.getEME2000(),
68 new AbsoluteDate(new DateComponents(2004, 01, 01),
69 new TimeComponents(23, 30, 00.000),
70 TimeScalesFactory.getUTC()),
71 Constants.EIGEN5C_EARTH_MU);
72 double mass = 5600.0;
73 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
74 Vector3D dV = new Vector3D(-0.1, 0.2, 0.3);
75 double f = 20.0;
76 double isp = 315.0;
77 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
78 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
79 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, 5,
80 new LofOffset(leo.getFrame(), LOFType.LVLH),
81 t0, Vector3D.ZERO, f, isp,
82 false, false, null);
83 BoundedPropagator withManeuver = getEphemeris(leo, mass, 5,
84 new LofOffset(leo.getFrame(), LOFType.LVLH),
85 t0, dV, f, isp,
86 false, false, null);
87
88
89 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
90 AdapterPropagator.DifferentialEffect effect =
91 new SmallManeuverAnalyticalModel(adapterPropagator.propagate(t0), dV.negate(), isp);
92 adapterPropagator.addEffect(effect);
93 adapterPropagator.addAdditionalDataProvider(new AdditionalDataProvider<double[]>() {
94 public String getName() {
95 return "dummy 3";
96 }
97 public double[] getAdditionalData(SpacecraftState state) {
98 return new double[3];
99 }
100 });
101
102
103
104 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 1"));
105 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 2"));
106 Assertions.assertTrue(adapterPropagator.isAdditionalDataManaged("dummy 3"));
107
108 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
109 t.compareTo(withoutManeuver.getMaxDate()) < 0;
110 t = t.shiftedBy(60.0)) {
111 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
112 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, leo.getFrame());
113 double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm();
114 Assertions.assertEquals(0, revertError, 0.45);
115 Assertions.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
116 Assertions.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
117 Assertions.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
118 }
119
120 }
121
122 @Test
123 public void testEccentricOrbit() throws ParseException, IOException {
124
125 Orbit heo = new KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0),
126 FastMath.toRadians(12.3456),
127 FastMath.toRadians(123.456),
128 FastMath.toRadians(1.23456), PositionAngleType.MEAN,
129 FramesFactory.getEME2000(),
130 new AbsoluteDate(new DateComponents(2004, 01, 01),
131 new TimeComponents(23, 30, 00.000),
132 TimeScalesFactory.getUTC()),
133 Constants.EIGEN5C_EARTH_MU);
134 double mass = 5600.0;
135 AbsoluteDate t0 = heo.getDate().shiftedBy(1000.0);
136 Vector3D dV = new Vector3D(-0.01, 0.02, 0.03);
137 double f = 20.0;
138 double isp = 315.0;
139 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
140 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
141 BoundedPropagator withoutManeuver = getEphemeris(heo, mass, 5,
142 new LofOffset(heo.getFrame(), LOFType.LVLH),
143 t0, Vector3D.ZERO, f, isp,
144 false, false, null);
145 BoundedPropagator withManeuver = getEphemeris(heo, mass, 5,
146 new LofOffset(heo.getFrame(), LOFType.LVLH),
147 t0, dV, f, isp,
148 false, false, null);
149
150
151 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
152 AdapterPropagator.DifferentialEffect effect =
153 new SmallManeuverAnalyticalModel(adapterPropagator.propagate(t0), dV.negate(), isp);
154 adapterPropagator.addEffect(effect);
155 adapterPropagator.addAdditionalDataProvider(new AdditionalDataProvider<double[]>() {
156 public String getName() {
157 return "dummy 3";
158 }
159 public double[] getAdditionalData(SpacecraftState state) {
160 return new double[3];
161 }
162 });
163
164
165
166 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 1"));
167 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 2"));
168 Assertions.assertTrue(adapterPropagator.isAdditionalDataManaged("dummy 3"));
169
170 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
171 t.compareTo(withoutManeuver.getMaxDate()) < 0;
172 t = t.shiftedBy(300.0)) {
173 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame());
174 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, heo.getFrame());
175 double revertError = Vector3D.distance(pvWithout.getPosition(), pvReverted.getPosition());
176 Assertions.assertEquals(0, revertError, 2.5e-5 * heo.getA());
177 Assertions.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
178 Assertions.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
179 Assertions.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
180 }
181
182 }
183
184 @Test
185 public void testNonKeplerian() throws ParseException, IOException {
186
187 Orbit leo = new CircularOrbit(7204319.233600575, 4.434564637450575E-4, 0.0011736728299091088,
188 1.7211611441767323, 5.5552084166959474,
189 24950.321259193086, PositionAngleType.TRUE,
190 FramesFactory.getEME2000(),
191 new AbsoluteDate(new DateComponents(2003, 9, 16),
192 new TimeComponents(23, 11, 20.264),
193 TimeScalesFactory.getUTC()),
194 Constants.EIGEN5C_EARTH_MU);
195 double mass = 4093.0;
196 AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2003, 9, 16),
197 new TimeComponents(23, 14, 40.264),
198 TimeScalesFactory.getUTC());
199 Vector3D dV = new Vector3D(0.0, 3.0, 0.0);
200 double f = 40.0;
201 double isp = 300.0;
202 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
203 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
204
205
206 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("g007_eigen_05c_coef", false));
207 NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getNormalizedProvider(8, 8);
208 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, 10,
209 new LofOffset(leo.getFrame(), LOFType.VNC),
210 t0, Vector3D.ZERO, f, isp,
211 true, true, gravityField);
212 BoundedPropagator withManeuver = getEphemeris(leo, mass, 10,
213 new LofOffset(leo.getFrame(), LOFType.VNC),
214 t0, dV, f, isp,
215 true, true, gravityField);
216
217
218 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
219 SpacecraftState state0 = adapterPropagator.propagate(t0);
220 AdapterPropagator.DifferentialEffect directEffect =
221 new SmallManeuverAnalyticalModel(state0, dV.negate(), isp);
222 AdapterPropagator.DifferentialEffect derivedEffect =
223 new J2DifferentialEffect(state0, directEffect, false,
224 GravityFieldFactory.getUnnormalizedProvider(gravityField));
225 adapterPropagator.addEffect(directEffect);
226 adapterPropagator.addEffect(derivedEffect);
227 adapterPropagator.addAdditionalDataProvider(new AdditionalDataProvider<double[]>() {
228 public String getName() {
229 return "dummy 3";
230 }
231 public double[] getAdditionalData(SpacecraftState state) {
232 return new double[3];
233 }
234 });
235
236
237
238 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 1"));
239 Assertions.assertFalse(adapterPropagator.isAdditionalDataManaged("dummy 2"));
240 Assertions.assertTrue(adapterPropagator.isAdditionalDataManaged("dummy 3"));
241
242 double maxDelta = 0;
243 double maxNominal = 0;
244 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
245 t.compareTo(withoutManeuver.getMaxDate()) < 0;
246 t = t.shiftedBy(600.0)) {
247 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
248 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, leo.getFrame());
249 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, leo.getFrame());
250 double nominal = new PVCoordinates(pvWithout, pvWith).getPosition().getNorm();
251 double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm();
252 maxDelta = FastMath.max(maxDelta, revertError);
253 maxNominal = FastMath.max(maxNominal, nominal);
254 Assertions.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
255 Assertions.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
256 Assertions.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
257 }
258 Assertions.assertTrue(maxDelta < 120);
259 Assertions.assertTrue(maxNominal > 2800);
260
261 }
262
263
264 @Test
265 public void testIssue949() {
266
267 final AbsoluteDate initialDate = new AbsoluteDate();
268 final Orbit initialOrbit = TestUtils.getDefaultOrbit(initialDate);
269
270
271 final Orbit finalOrbit = initialOrbit.shiftedBy(1);
272 final List<SpacecraftState> states = new ArrayList<>();
273 states.add(new SpacecraftState(initialOrbit));
274 states.add(new SpacecraftState(finalOrbit));
275
276 final Ephemeris ephemeris = new Ephemeris(states, 2);
277 final AdapterPropagator adapterPropagator = new AdapterPropagator(ephemeris);
278
279
280 final AdditionalDataProvider<double[]> additionalDataProvider = TestUtils.getAdditionalProviderWithInit();
281 adapterPropagator.addAdditionalDataProvider(additionalDataProvider);
282
283
284 Assertions.assertDoesNotThrow(() -> adapterPropagator.propagate(finalOrbit.getDate()), "No error should have been thrown");
285
286 }
287
288 private BoundedPropagator getEphemeris(final Orbit orbit, final double mass, final int nbOrbits,
289 final AttitudeProvider law,
290 final AbsoluteDate t0, final Vector3D dV,
291 final double f, final double isp,
292 final boolean sunAttraction, final boolean moonAttraction,
293 final NormalizedSphericalHarmonicsProvider gravityField)
294 throws ParseException, IOException {
295
296 SpacecraftState initialState =
297 new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
298
299
300 initialState = initialState.addAdditionalData("dummy 1", new double[]{1.25, 2.5});
301 initialState = initialState.addAdditionalData("dummy 2", 5.0);
302
303
304 final double dP = 1.0;
305 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, orbit.getType());
306 AdaptiveStepsizeIntegrator integrator =
307 new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
308 integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
309 final NumericalPropagator propagator = new NumericalPropagator(integrator);
310 propagator.addAdditionalDataProvider(new AdditionalDataProvider<Double>() {
311 public String getName() {
312 return "dummy 2";
313 }
314 public Double getAdditionalData(SpacecraftState state) {
315 return 5.0;
316 }
317 });
318 propagator.setInitialState(initialState);
319 propagator.setAttitudeProvider(law);
320
321 if (dV.getNorm() > 1.0e-6) {
322
323 final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
324 final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
325 final ConstantThrustManeuver maneuver =
326 new ConstantThrustManeuver(t0.shiftedBy(-0.5 * dt), dt , f, isp, dV.normalize());
327 propagator.addForceModel(maneuver);
328 }
329
330 if (sunAttraction) {
331 propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
332 }
333
334 if (moonAttraction) {
335 propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
336 }
337
338 if (gravityField != null) {
339 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getGTOD(false),
340 gravityField));
341 }
342
343 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
344 propagator.propagate(t0.shiftedBy(nbOrbits * orbit.getKeplerianPeriod()));
345
346 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
347
348
349
350 Assertions.assertFalse(propagator.isAdditionalDataManaged("dummy 1"));
351 Assertions.assertTrue(propagator.isAdditionalDataManaged("dummy 2"));
352 Assertions.assertFalse(ephemeris.isAdditionalDataManaged("dummy 1"));
353 Assertions.assertTrue(ephemeris.isAdditionalDataManaged("dummy 2"));
354 Assertions.assertEquals(2, ephemeris.getInitialState().getAdditionalState("dummy 1").length);
355 Assertions.assertEquals(1, ephemeris.getInitialState().getAdditionalState("dummy 2").length);
356
357 return ephemeris;
358
359 }
360
361 @BeforeEach
362 public void setUp() {
363 Utils.setDataRoot("regular-data:potential/icgem-format");
364 }
365
366 }