1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.maneuvers;
18
19 import org.hipparchus.geometry.euclidean.threed.Rotation;
20 import org.hipparchus.geometry.euclidean.threed.RotationConvention;
21 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.linear.RealMatrix;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.FastMath;
26 import org.junit.jupiter.api.Assertions;
27 import org.junit.jupiter.api.BeforeEach;
28 import org.junit.jupiter.api.Test;
29 import org.orekit.Utils;
30 import org.orekit.attitudes.Attitude;
31 import org.orekit.attitudes.AttitudeProvider;
32 import org.orekit.attitudes.FrameAlignedProvider;
33 import org.orekit.attitudes.LofOffset;
34 import org.orekit.bodies.CelestialBodyFactory;
35 import org.orekit.frames.Frame;
36 import org.orekit.frames.FramesFactory;
37 import org.orekit.frames.LOFType;
38 import org.orekit.orbits.CartesianOrbit;
39 import org.orekit.orbits.KeplerianOrbit;
40 import org.orekit.orbits.Orbit;
41 import org.orekit.orbits.PositionAngleType;
42 import org.orekit.propagation.MatricesHarvester;
43 import org.orekit.propagation.SpacecraftState;
44 import org.orekit.propagation.analytical.KeplerianPropagator;
45 import org.orekit.propagation.events.AbstractDetector;
46 import org.orekit.propagation.events.AdaptableInterval;
47 import org.orekit.propagation.events.DateDetector;
48 import org.orekit.propagation.events.NodeDetector;
49 import org.orekit.propagation.events.handlers.EventHandler;
50 import org.orekit.propagation.events.handlers.StopOnIncreasing;
51 import org.orekit.propagation.numerical.NumericalPropagator;
52 import org.orekit.time.AbsoluteDate;
53 import org.orekit.time.DateComponents;
54 import org.orekit.time.TimeComponents;
55 import org.orekit.time.TimeScalesFactory;
56 import org.orekit.utils.Constants;
57 import org.orekit.utils.TimeStampedPVCoordinates;
58
59 public class ImpulseManeuverTest {
60 @Test
61 public void testInclinationManeuver() {
62 final Orbit initialOrbit =
63 new KeplerianOrbit(24532000.0, 0.72, 0.3, FastMath.PI, 0.4, 2.0 + 4 * FastMath.PI,
64 PositionAngleType.MEAN, FramesFactory.getEME2000(),
65 new AbsoluteDate(new DateComponents(2008, 06, 23),
66 new TimeComponents(14, 18, 37),
67 TimeScalesFactory.getUTC()),
68 3.986004415e14);
69 final double a = initialOrbit.getA();
70 final double e = initialOrbit.getE();
71 final double i = initialOrbit.getI();
72 final double mu = initialOrbit.getMu();
73 final double vApo = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
74 double dv = 0.99 * FastMath.tan(i) * vApo;
75 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit,
76 new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS));
77 propagator.addEventDetector(new ImpulseManeuver(new NodeDetector(initialOrbit, FramesFactory.getEME2000()),
78 new Vector3D(dv, Vector3D.PLUS_J), 400.0));
79 SpacecraftState propagated = propagator.propagate(initialOrbit.getDate().shiftedBy(8000));
80 Assertions.assertEquals(0.0028257, propagated.getI(), 1.0e-6);
81 Assertions.assertEquals(0.442476 + 6 * FastMath.PI, ((KeplerianOrbit) propagated.getOrbit()).getLv(), 1.0e-6);
82 }
83
84 @Test
85 public void testInertialManeuver() {
86 final double mu = CelestialBodyFactory.getEarth().getGM();
87
88 final double initialX = 7100e3;
89 final double initialY = 0.0;
90 final double initialZ = 1300e3;
91 final double initialVx = 0;
92 final double initialVy = 8000;
93 final double initialVz = 1000;
94
95 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
96 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
97 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
98 final TimeStampedPVCoordinates state = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
99 final Orbit initialOrbit = new CartesianOrbit(state, FramesFactory.getEME2000(), mu);
100
101 final double totalPropagationTime = 0.00001;
102 final double driftTimeInSec = totalPropagationTime / 2.0;
103 final double deltaX = 0.01;
104 final double deltaY = 0.02;
105 final double deltaZ = 0.03;
106 final double isp = 300;
107
108 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
109
110 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit, new LofOffset(initialOrbit.getFrame(), LOFType.VNC));
111 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(driftTimeInSec));
112 FrameAlignedProvider attitudeOverride = new FrameAlignedProvider(new Rotation(RotationOrder.XYX,
113 RotationConvention.VECTOR_OPERATOR,
114 0, 0, 0));
115 ImpulseManeuver burnAtEpoch = new ImpulseManeuver(dateDetector, attitudeOverride, deltaV, isp).withThreshold(driftTimeInSec/4);
116 Assertions.assertEquals(0.0, Vector3D.distance(deltaV, burnAtEpoch.getDeltaVSat()), 1.0e-15);
117 Assertions.assertEquals(isp, burnAtEpoch.getIsp(), 1.0e-15);
118 Assertions.assertSame(dateDetector, burnAtEpoch.getTrigger());
119 propagator.addEventDetector(burnAtEpoch);
120
121 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
122
123 final double finalVxExpected = initialVx + deltaX;
124 final double finalVyExpected = initialVy + deltaY;
125 final double finalVzExpected = initialVz + deltaZ;
126 final double maneuverTolerance = 1e-4;
127
128 final Vector3D finalVelocity = finalState.getPVCoordinates().getVelocity();
129 Assertions.assertEquals(finalVxExpected, finalVelocity.getX(), maneuverTolerance);
130 Assertions.assertEquals(finalVyExpected, finalVelocity.getY(), maneuverTolerance);
131 Assertions.assertEquals(finalVzExpected, finalVelocity.getZ(), maneuverTolerance);
132
133 }
134
135 @Test
136 public void testBackward() {
137
138 final AbsoluteDate iniDate = new AbsoluteDate(2003, 5, 1, 17, 30, 0.0, TimeScalesFactory.getUTC());
139 final Orbit initialOrbit = new KeplerianOrbit(7e6, 1.0e-4, FastMath.toRadians(98.5),
140 FastMath.toRadians(87.0), FastMath.toRadians(216.1807),
141 FastMath.toRadians(319.779), PositionAngleType.MEAN,
142 FramesFactory.getEME2000(), iniDate,
143 Constants.EIGEN5C_EARTH_MU);
144 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit,
145 new LofOffset(initialOrbit.getFrame(),
146 LOFType.VNC));
147 DateDetector dateDetector = new DateDetector(iniDate.shiftedBy(-300));
148 Vector3D deltaV = new Vector3D(12.0, 1.0, -4.0);
149 final double isp = 300;
150 ImpulseManeuver maneuver =
151 new ImpulseManeuver(dateDetector, deltaV, isp).
152 withMaxCheck(3600.0).
153 withThreshold(1.0e-6);
154 propagator.addEventDetector(maneuver);
155
156 SpacecraftState finalState = propagator.propagate(initialOrbit.getDate().shiftedBy(-900));
157
158 Assertions.assertTrue(finalState.getMass() > propagator.getInitialState().getMass());
159 Assertions.assertTrue(finalState.getDate().compareTo(propagator.getInitialState().getDate()) < 0);
160
161 }
162
163 @Test
164 public void testBackAndForth() {
165
166 final AttitudeProvider lof = new LofOffset(FramesFactory.getEME2000(), LOFType.VNC);
167 final double mu = Constants.EIGEN5C_EARTH_MU;
168 final AbsoluteDate iniDate = new AbsoluteDate(2003, 5, 1, 17, 30, 0.0, TimeScalesFactory.getUTC());
169 final Orbit pastOrbit = new KeplerianOrbit(7e6, 1.0e-4, FastMath.toRadians(98.5),
170 FastMath.toRadians(87.0), FastMath.toRadians(216.1807),
171 FastMath.toRadians(319.779), PositionAngleType.MEAN,
172 FramesFactory.getEME2000(), iniDate, mu);
173 final double pastMass = 2500.0;
174 DateDetector dateDetector = new DateDetector(iniDate.shiftedBy(600));
175 Vector3D deltaV = new Vector3D(12.0, 1.0, -4.0);
176 final double isp = 300;
177 ImpulseManeuver maneuver =
178 new ImpulseManeuver(dateDetector,
179 new FrameAlignedProvider(Rotation.IDENTITY),
180 deltaV, isp).
181 withMaxCheck(3600.0).
182 withThreshold(1.0e-6);
183
184 double span = 900.0;
185 KeplerianPropagator forwardPropagator = new KeplerianPropagator(pastOrbit, lof, mu, pastMass);
186 forwardPropagator.addEventDetector(maneuver);
187 SpacecraftState futureState = forwardPropagator.propagate(pastOrbit.getDate().shiftedBy(span));
188
189 KeplerianPropagator backwardPropagator = new KeplerianPropagator(futureState.getOrbit(), lof,
190 mu, futureState.getMass());
191 backwardPropagator.addEventDetector(maneuver);
192 SpacecraftState rebuiltPast = backwardPropagator.propagate(pastOrbit.getDate());
193 Assertions.assertEquals(0.0,
194 Vector3D.distance(pastOrbit.getPosition(),
195 rebuiltPast.getPosition()),
196 2.0e-8);
197 Assertions.assertEquals(0.0,
198 Vector3D.distance(pastOrbit.getPVCoordinates().getVelocity(),
199 rebuiltPast.getPVCoordinates().getVelocity()),
200 2.0e-11);
201 Assertions.assertEquals(pastMass, rebuiltPast.getMass(), 5.0e-13);
202
203 }
204
205 @Test
206 public void testAdditionalStateKeplerian() {
207 final double mu = CelestialBodyFactory.getEarth().getGM();
208
209 final double initialX = 7100e3;
210 final double initialY = 0.0;
211 final double initialZ = 1300e3;
212 final double initialVx = 0;
213 final double initialVy = 8000;
214 final double initialVz = 1000;
215
216 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
217 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
218 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
219 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
220 final Orbit initialOrbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), mu);
221
222 final double totalPropagationTime = 10;
223 final double deltaX = 0.01;
224 final double deltaY = 0.02;
225 final double deltaZ = 0.03;
226 final double isp = 300;
227
228 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
229
230 final AttitudeProvider attitudeProvider = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
231 final Attitude initialAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
232 final SpacecraftState initialState = new SpacecraftState(initialOrbit, initialAttitude);
233 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);
234 propagator.resetInitialState(initialState.addAdditionalState("testOnly", -1.0));
235 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(0.5 * totalPropagationTime));
236 FrameAlignedProvider attitudeOverride = new FrameAlignedProvider(new Rotation(RotationOrder.XYX,
237 RotationConvention.VECTOR_OPERATOR,
238 0, 0, 0));
239 ImpulseManeuver burnAtEpoch = new ImpulseManeuver(dateDetector, attitudeOverride, deltaV, isp).withThreshold(1.0e-3);
240 propagator.addEventDetector(burnAtEpoch);
241
242 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
243 Assertions.assertEquals(1, finalState.getAdditionalStatesValues().size());
244 Assertions.assertEquals(-1.0, finalState.getAdditionalState("testOnly")[0], 1.0e-15);
245
246 }
247
248 @Test
249 public void testAdditionalStateNumerical() {
250 final double mu = CelestialBodyFactory.getEarth().getGM();
251
252 final double initialX = 7100e3;
253 final double initialY = 0.0;
254 final double initialZ = 1300e3;
255 final double initialVx = 0;
256 final double initialVy = 8000;
257 final double initialVz = 1000;
258
259 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
260 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
261 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
262 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
263 final Orbit initialOrbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), mu);
264
265 final double totalPropagationTime = 10.0;
266 final double deltaX = 0.01;
267 final double deltaY = 0.02;
268 final double deltaZ = 0.03;
269 final double isp = 300;
270
271 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
272
273 final AttitudeProvider attitudeProvider = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
274 final Attitude initialAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
275
276 double[][] tolerances = NumericalPropagator.tolerances(10.0, initialOrbit, initialOrbit.getType());
277 DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-3, 60, tolerances[0], tolerances[1]);
278 NumericalPropagator propagator = new NumericalPropagator(integrator);
279 propagator.setOrbitType(initialOrbit.getType());
280 MatricesHarvester harvester = propagator.setupMatricesComputation("derivatives", null, null);
281 propagator.resetInitialState(new SpacecraftState(initialOrbit, initialAttitude));
282 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(0.5 * totalPropagationTime));
283 FrameAlignedProvider attitudeOverride = new FrameAlignedProvider(new Rotation(RotationOrder.XYX,
284 RotationConvention.VECTOR_OPERATOR,
285 0, 0, 0));
286 ImpulseManeuver burnAtEpoch = new ImpulseManeuver(dateDetector, attitudeOverride, deltaV, isp).withThreshold(1.0e-3);
287 propagator.addEventDetector(burnAtEpoch);
288
289 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
290 Assertions.assertEquals(1, finalState.getAdditionalStatesValues().size());
291 Assertions.assertEquals(36, finalState.getAdditionalState("derivatives").length);
292
293 RealMatrix stateTransitionMatrix = harvester.getStateTransitionMatrix(finalState);
294 for (int i = 0; i < 6; ++i) {
295 for (int j = 0; j < 6; ++j) {
296 double sIJ = stateTransitionMatrix.getEntry(i, j);
297 if (j == i) {
298
299 Assertions.assertEquals(1.0, sIJ, 2.0e-4);
300 } else if (j == i + 3) {
301
302 Assertions.assertEquals(totalPropagationTime, sIJ, 4.0e-5 * totalPropagationTime);
303 } else {
304
305 Assertions.assertEquals(0, sIJ, 1.0e-4);
306 }
307 }
308 }
309
310 }
311
312
313
314
315
316 @Test
317 public void testIssue663() {
318
319
320 final Orbit initialOrbit =
321 new KeplerianOrbit(24532000.0, 0.72, 0.3, FastMath.PI, 0.4, 2.0,
322 PositionAngleType.MEAN, FramesFactory.getEME2000(),
323 new AbsoluteDate(new DateComponents(2008, 06, 23),
324 new TimeComponents(14, 18, 37),
325 TimeScalesFactory.getUTC()),
326 3.986004415e14);
327
328 final InitializationDetector trigger = new InitializationDetector();
329
330
331 final AttitudeProvider attitudeOverride = new LofOffset(FramesFactory.getEME2000(), LOFType.TNW);
332 final Vector3D deltaVSat = Vector3D.PLUS_I;
333 final double isp = 1500.0;
334 final ImpulseManeuver maneuver = new ImpulseManeuver(trigger, attitudeOverride, deltaVSat, isp);
335
336
337 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit, attitudeOverride);
338 propagator.addEventDetector(maneuver);
339
340
341 Assertions.assertFalse(trigger.initialized);
342 propagator.propagate(initialOrbit.getDate().shiftedBy(3600.0));
343 Assertions.assertTrue(trigger.initialized);
344
345 }
346
347 @Test
348 public void testControl3DVectorCostType() {
349
350
351
352 final double a = Constants.EGM96_EARTH_EQUATORIAL_RADIUS + 2000.e3;
353 final double e = 1e-4;
354 final double i = 0.5;
355 final double pa = 0.;
356 final double raan = 6.;
357 final double anomaly = 1.;
358 final PositionAngleType positionAngleType = PositionAngleType.MEAN;
359 final Frame gcrf = FramesFactory.getGCRF();
360 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
361 final double mu = Constants.EGM96_EARTH_MU;
362
363 final KeplerianOrbit orbit = new KeplerianOrbit(a, e, i, pa, raan, anomaly, positionAngleType, gcrf, date, mu);
364
365
366 final DateDetector dateDetector = new DateDetector(date.shiftedBy(1000.));
367 final AbsoluteDate endPropagationDate = dateDetector.getDate().shiftedBy(10000.);
368 final AttitudeProvider provider = new FrameAlignedProvider(gcrf);
369 final Vector3D deltaV = new Vector3D(2., -1., 0.5);
370 final double initialMass = 1000.;
371 final double isp = 100.;
372
373
374 final KeplerianPropagator propagatorNone = new KeplerianPropagator(orbit, provider, orbit.getMu(), initialMass);
375 final KeplerianPropagator propagatorNorm1 = new KeplerianPropagator(orbit, provider, orbit.getMu(), initialMass);
376 final KeplerianPropagator propagatorNorm2 = new KeplerianPropagator(orbit, provider, orbit.getMu(), initialMass);
377 final KeplerianPropagator propagatorNormInf = new KeplerianPropagator(orbit, provider, orbit.getMu(), initialMass);
378
379
380 propagatorNone.addEventDetector(new ImpulseManeuver(dateDetector, provider, deltaV, isp,
381 Control3DVectorCostType.NONE));
382 propagatorNorm1.addEventDetector(new ImpulseManeuver(dateDetector, provider, deltaV, isp,
383 Control3DVectorCostType.ONE_NORM));
384 propagatorNorm2.addEventDetector(new ImpulseManeuver(dateDetector, provider, deltaV, isp,
385 Control3DVectorCostType.TWO_NORM));
386 propagatorNormInf.addEventDetector(new ImpulseManeuver(dateDetector, provider, deltaV, isp,
387 Control3DVectorCostType.INF_NORM));
388
389
390 final double finalMassWithNone = propagatorNone.propagate(endPropagationDate).getMass();
391 final double finalMassWithNorm1 = propagatorNorm1.propagate(endPropagationDate).getMass();
392 final double finalMassWithNorm2 = propagatorNorm2.propagate(endPropagationDate).getMass();
393 final double finalMassWithNormInf = propagatorNormInf.propagate(endPropagationDate).getMass();
394
395
396
397 Assertions.assertNotEquals(finalMassWithNorm1, finalMassWithNorm2);
398 Assertions.assertNotEquals(finalMassWithNorm1, finalMassWithNormInf);
399 Assertions.assertNotEquals(finalMassWithNormInf, finalMassWithNorm2);
400
401
402 Assertions.assertEquals(initialMass, finalMassWithNone);
403 final double factorExponential = -1. / (isp * Constants.G0_STANDARD_GRAVITY);
404 Assertions.assertEquals(initialMass * FastMath.exp(deltaV.getNorm1() * factorExponential), finalMassWithNorm1);
405 Assertions.assertEquals(initialMass * FastMath.exp(deltaV.getNorm() * factorExponential), finalMassWithNorm2);
406 Assertions.assertEquals(initialMass * FastMath.exp(deltaV.getNormInf() * factorExponential), finalMassWithNormInf);
407 }
408
409
410 @BeforeEach
411 public void setUp() {
412 Utils.setDataRoot("regular-data");
413 }
414
415
416 private class InitializationDetector extends AbstractDetector<InitializationDetector> {
417
418
419 private boolean initialized;
420
421
422
423
424 InitializationDetector() {
425 super(AbstractDetector.DEFAULT_MAXCHECK, AbstractDetector.DEFAULT_THRESHOLD,
426 AbstractDetector.DEFAULT_MAX_ITER, new StopOnIncreasing());
427 this.initialized = false;
428 }
429
430
431 @Override
432 public void init(final SpacecraftState s0, final AbsoluteDate t) {
433 super.init(s0, t);
434 this.initialized = true;
435 }
436
437
438 @Override
439 public double g(SpacecraftState s) {
440 return 1;
441 }
442
443
444 @Override
445 protected InitializationDetector create(AdaptableInterval newMaxCheck, double newThreshold, int newMaxIter, EventHandler newHandler) {
446 return new InitializationDetector();
447 }
448
449 }
450
451 }