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