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.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.FastMath;
26 import org.junit.Assert;
27 import org.junit.Before;
28 import org.junit.Test;
29 import org.orekit.Utils;
30 import org.orekit.attitudes.Attitude;
31 import org.orekit.attitudes.AttitudeProvider;
32 import org.orekit.attitudes.InertialProvider;
33 import org.orekit.attitudes.LofOffset;
34 import org.orekit.bodies.CelestialBodyFactory;
35 import org.orekit.frames.FramesFactory;
36 import org.orekit.frames.LOFType;
37 import org.orekit.orbits.CartesianOrbit;
38 import org.orekit.orbits.KeplerianOrbit;
39 import org.orekit.orbits.Orbit;
40 import org.orekit.orbits.PositionAngle;
41 import org.orekit.propagation.SpacecraftState;
42 import org.orekit.propagation.analytical.KeplerianPropagator;
43 import org.orekit.propagation.events.AbstractDetector;
44 import org.orekit.propagation.events.DateDetector;
45 import org.orekit.propagation.events.NodeDetector;
46 import org.orekit.propagation.events.handlers.EventHandler;
47 import org.orekit.propagation.events.handlers.StopOnIncreasing;
48 import org.orekit.propagation.numerical.NumericalPropagator;
49 import org.orekit.propagation.numerical.PartialDerivativesEquations;
50 import org.orekit.time.AbsoluteDate;
51 import org.orekit.time.DateComponents;
52 import org.orekit.time.TimeComponents;
53 import org.orekit.time.TimeScalesFactory;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.TimeStampedPVCoordinates;
56
57 public class ImpulseManeuverTest {
58 @Test
59 public void testInclinationManeuver() {
60 final Orbit initialOrbit =
61 new KeplerianOrbit(24532000.0, 0.72, 0.3, FastMath.PI, 0.4, 2.0,
62 PositionAngle.MEAN, FramesFactory.getEME2000(),
63 new AbsoluteDate(new DateComponents(2008, 06, 23),
64 new TimeComponents(14, 18, 37),
65 TimeScalesFactory.getUTC()),
66 3.986004415e14);
67 final double a = initialOrbit.getA();
68 final double e = initialOrbit.getE();
69 final double i = initialOrbit.getI();
70 final double mu = initialOrbit.getMu();
71 final double vApo = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
72 double dv = 0.99 * FastMath.tan(i) * vApo;
73 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit,
74 new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS));
75 propagator.addEventDetector(new ImpulseManeuver<NodeDetector>(new NodeDetector(initialOrbit, FramesFactory.getEME2000()),
76 new Vector3D(dv, Vector3D.PLUS_J), 400.0));
77 SpacecraftState propagated = propagator.propagate(initialOrbit.getDate().shiftedBy(8000));
78 Assert.assertEquals(0.0028257, propagated.getI(), 1.0e-6);
79 }
80
81 @Test
82 public void testInertialManeuver() {
83 final double mu = CelestialBodyFactory.getEarth().getGM();
84
85 final double initialX = 7100e3;
86 final double initialY = 0.0;
87 final double initialZ = 1300e3;
88 final double initialVx = 0;
89 final double initialVy = 8000;
90 final double initialVz = 1000;
91
92 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
93 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
94 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
95 final TimeStampedPVCoordinates state = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
96 final Orbit initialOrbit = new CartesianOrbit(state, FramesFactory.getEME2000(), mu);
97
98 final double totalPropagationTime = 0.00001;
99 final double driftTimeInSec = totalPropagationTime / 2.0;
100 final double deltaX = 0.01;
101 final double deltaY = 0.02;
102 final double deltaZ = 0.03;
103 final double isp = 300;
104
105 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
106
107 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit, new LofOffset(initialOrbit.getFrame(), LOFType.VNC));
108 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(driftTimeInSec));
109 InertialProvider attitudeOverride = new InertialProvider(new Rotation(RotationOrder.XYX,
110 RotationConvention.VECTOR_OPERATOR,
111 0, 0, 0));
112 ImpulseManeuver<DateDetector> burnAtEpoch = new ImpulseManeuver<DateDetector>(dateDetector, attitudeOverride, deltaV, isp).withThreshold(driftTimeInSec/4);
113 Assert.assertEquals(0.0, Vector3D.distance(deltaV, burnAtEpoch.getDeltaVSat()), 1.0e-15);
114 Assert.assertEquals(isp, burnAtEpoch.getIsp(), 1.0e-15);
115 Assert.assertSame(dateDetector, burnAtEpoch.getTrigger());
116 propagator.addEventDetector(burnAtEpoch);
117
118 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
119
120 final double finalVxExpected = initialVx + deltaX;
121 final double finalVyExpected = initialVy + deltaY;
122 final double finalVzExpected = initialVz + deltaZ;
123 final double maneuverTolerance = 1e-4;
124
125 final Vector3D finalVelocity = finalState.getPVCoordinates().getVelocity();
126 Assert.assertEquals(finalVxExpected, finalVelocity.getX(), maneuverTolerance);
127 Assert.assertEquals(finalVyExpected, finalVelocity.getY(), maneuverTolerance);
128 Assert.assertEquals(finalVzExpected, finalVelocity.getZ(), maneuverTolerance);
129
130 }
131
132 @Test
133 public void testBackward() {
134
135 final AbsoluteDate iniDate = new AbsoluteDate(2003, 5, 1, 17, 30, 0.0, TimeScalesFactory.getUTC());
136 final Orbit initialOrbit = new KeplerianOrbit(7e6, 1.0e-4, FastMath.toRadians(98.5),
137 FastMath.toRadians(87.0), FastMath.toRadians(216.1807),
138 FastMath.toRadians(319.779), PositionAngle.MEAN,
139 FramesFactory.getEME2000(), iniDate,
140 Constants.EIGEN5C_EARTH_MU);
141 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit,
142 new LofOffset(initialOrbit.getFrame(),
143 LOFType.VNC));
144 DateDetector dateDetector = new DateDetector(iniDate.shiftedBy(-300));
145 Vector3D deltaV = new Vector3D(12.0, 1.0, -4.0);
146 final double isp = 300;
147 ImpulseManeuver<DateDetector> maneuver =
148 new ImpulseManeuver<DateDetector>(dateDetector, deltaV, isp).
149 withMaxCheck(3600.0).
150 withThreshold(1.0e-6);
151 propagator.addEventDetector(maneuver);
152
153 SpacecraftState finalState = propagator.propagate(initialOrbit.getDate().shiftedBy(-900));
154
155 Assert.assertTrue(finalState.getMass() > propagator.getInitialState().getMass());
156 Assert.assertTrue(finalState.getDate().compareTo(propagator.getInitialState().getDate()) < 0);
157
158 }
159
160 @Test
161 public void testBackAndForth() {
162
163 final AttitudeProvider lof = new LofOffset(FramesFactory.getEME2000(), LOFType.VNC);
164 final double mu = Constants.EIGEN5C_EARTH_MU;
165 final AbsoluteDate iniDate = new AbsoluteDate(2003, 5, 1, 17, 30, 0.0, TimeScalesFactory.getUTC());
166 final Orbit pastOrbit = new KeplerianOrbit(7e6, 1.0e-4, FastMath.toRadians(98.5),
167 FastMath.toRadians(87.0), FastMath.toRadians(216.1807),
168 FastMath.toRadians(319.779), PositionAngle.MEAN,
169 FramesFactory.getEME2000(), iniDate, mu);
170 final double pastMass = 2500.0;
171 DateDetector dateDetector = new DateDetector(iniDate.shiftedBy(600));
172 Vector3D deltaV = new Vector3D(12.0, 1.0, -4.0);
173 final double isp = 300;
174 ImpulseManeuver<DateDetector> maneuver =
175 new ImpulseManeuver<DateDetector>(dateDetector,
176 new InertialProvider(Rotation.IDENTITY),
177 deltaV, isp).
178 withMaxCheck(3600.0).
179 withThreshold(1.0e-6);
180
181 double span = 900.0;
182 KeplerianPropagator forwardPropagator = new KeplerianPropagator(pastOrbit, lof, mu, pastMass);
183 forwardPropagator.addEventDetector(maneuver);
184 SpacecraftState futureState = forwardPropagator.propagate(pastOrbit.getDate().shiftedBy(span));
185
186 KeplerianPropagator backwardPropagator = new KeplerianPropagator(futureState.getOrbit(), lof,
187 mu, futureState.getMass());
188 backwardPropagator.addEventDetector(maneuver);
189 SpacecraftState rebuiltPast = backwardPropagator.propagate(pastOrbit.getDate());
190 Assert.assertEquals(0.0,
191 Vector3D.distance(pastOrbit.getPVCoordinates().getPosition(),
192 rebuiltPast.getPVCoordinates().getPosition()),
193 2.0e-8);
194 Assert.assertEquals(0.0,
195 Vector3D.distance(pastOrbit.getPVCoordinates().getVelocity(),
196 rebuiltPast.getPVCoordinates().getVelocity()),
197 2.0e-11);
198 Assert.assertEquals(pastMass, rebuiltPast.getMass(), 5.0e-13);
199
200 }
201
202 @Test
203 public void testAdditionalStateKeplerian() {
204 final double mu = CelestialBodyFactory.getEarth().getGM();
205
206 final double initialX = 7100e3;
207 final double initialY = 0.0;
208 final double initialZ = 1300e3;
209 final double initialVx = 0;
210 final double initialVy = 8000;
211 final double initialVz = 1000;
212
213 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
214 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
215 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
216 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
217 final Orbit initialOrbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), mu);
218
219 final double totalPropagationTime = 10;
220 final double deltaX = 0.01;
221 final double deltaY = 0.02;
222 final double deltaZ = 0.03;
223 final double isp = 300;
224
225 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
226
227 final AttitudeProvider attitudeProvider = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
228 final Attitude initialAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
229 final SpacecraftState initialState = new SpacecraftState(initialOrbit, initialAttitude);
230 KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit);
231 propagator.resetInitialState(initialState.addAdditionalState("testOnly", -1.0));
232 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(0.5 * totalPropagationTime));
233 InertialProvider attitudeOverride = new InertialProvider(new Rotation(RotationOrder.XYX,
234 RotationConvention.VECTOR_OPERATOR,
235 0, 0, 0));
236 ImpulseManeuver<DateDetector> burnAtEpoch = new ImpulseManeuver<DateDetector>(dateDetector, attitudeOverride, deltaV, isp).withThreshold(1.0e-3);
237 propagator.addEventDetector(burnAtEpoch);
238
239 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
240 Assert.assertEquals(1, finalState.getAdditionalStates().size());
241 Assert.assertEquals(-1.0, finalState.getAdditionalState("testOnly")[0], 1.0e-15);
242
243 }
244
245 @Test
246 public void testAdditionalStateNumerical() {
247 final double mu = CelestialBodyFactory.getEarth().getGM();
248
249 final double initialX = 7100e3;
250 final double initialY = 0.0;
251 final double initialZ = 1300e3;
252 final double initialVx = 0;
253 final double initialVy = 8000;
254 final double initialVz = 1000;
255
256 final Vector3D position = new Vector3D(initialX, initialY, initialZ);
257 final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
258 final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
259 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
260 final Orbit initialOrbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), mu);
261
262 final double totalPropagationTime = 10.0;
263 final double deltaX = 0.01;
264 final double deltaY = 0.02;
265 final double deltaZ = 0.03;
266 final double isp = 300;
267
268 final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
269
270 final AttitudeProvider attitudeProvider = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
271 final Attitude initialAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
272
273 double[][] tolerances = NumericalPropagator.tolerances(10.0, initialOrbit, initialOrbit.getType());
274 DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-3, 60, tolerances[0], tolerances[1]);
275 NumericalPropagator propagator = new NumericalPropagator(integrator);
276 propagator.setOrbitType(initialOrbit.getType());
277 PartialDerivativesEquations pde = new PartialDerivativesEquations("derivatives", propagator);
278 final SpacecraftState initialState = pde.setInitialJacobians(new SpacecraftState(initialOrbit, initialAttitude));
279 propagator.resetInitialState(initialState);
280 DateDetector dateDetector = new DateDetector(epoch.shiftedBy(0.5 * totalPropagationTime));
281 InertialProvider attitudeOverride = new InertialProvider(new Rotation(RotationOrder.XYX,
282 RotationConvention.VECTOR_OPERATOR,
283 0, 0, 0));
284 ImpulseManeuver<DateDetector> burnAtEpoch = new ImpulseManeuver<DateDetector>(dateDetector, attitudeOverride, deltaV, isp).withThreshold(1.0e-3);
285 propagator.addEventDetector(burnAtEpoch);
286
287 SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
288 Assert.assertEquals(1, finalState.getAdditionalStates().size());
289 Assert.assertEquals(36, finalState.getAdditionalState("derivatives").length);
290
291 double[][] stateTransitionMatrix = new double[6][6];
292 pde.getMapper().getStateJacobian(finalState, stateTransitionMatrix);
293 for (int i = 0; i < 6; ++i) {
294 for (int j = 0; j < 6; ++j) {
295 double sIJ = stateTransitionMatrix[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 }