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.Vector3D;
21 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
22 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
23 import org.hipparchus.util.FastMath;
24 import org.junit.Assert;
25 import org.junit.Before;
26 import org.junit.Test;
27 import org.orekit.Utils;
28 import org.orekit.attitudes.AttitudeProvider;
29 import org.orekit.attitudes.LofOffset;
30 import org.orekit.frames.Frame;
31 import org.orekit.frames.FramesFactory;
32 import org.orekit.frames.LOFType;
33 import org.orekit.orbits.CircularOrbit;
34 import org.orekit.orbits.KeplerianOrbit;
35 import org.orekit.orbits.Orbit;
36 import org.orekit.orbits.OrbitType;
37 import org.orekit.orbits.PositionAngle;
38 import org.orekit.propagation.BoundedPropagator;
39 import org.orekit.propagation.EphemerisGenerator;
40 import org.orekit.propagation.SpacecraftState;
41 import org.orekit.propagation.numerical.NumericalPropagator;
42 import org.orekit.time.AbsoluteDate;
43 import org.orekit.time.DateComponents;
44 import org.orekit.time.TimeComponents;
45 import org.orekit.time.TimeScalesFactory;
46 import org.orekit.utils.Constants;
47 import org.orekit.utils.PVCoordinates;
48
49 public class SmallManeuverAnalyticalModelTest {
50
51 @Test
52 public void testLowEarthOrbit1() {
53
54 Orbit leo = new CircularOrbit(7200000.0, -1.0e-5, 2.0e-4,
55 FastMath.toRadians(98.0),
56 FastMath.toRadians(123.456),
57 0.0, PositionAngle.MEAN,
58 FramesFactory.getEME2000(),
59 new AbsoluteDate(new DateComponents(2004, 01, 01),
60 new TimeComponents(23, 30, 00.000),
61 TimeScalesFactory.getUTC()),
62 Constants.EIGEN5C_EARTH_MU);
63 double mass = 5600.0;
64 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
65 Vector3D dV = new Vector3D(-0.01, 0.02, 0.03);
66 double f = 20.0;
67 double isp = 315.0;
68 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, t0, Vector3D.ZERO, f, isp);
69 BoundedPropagator withManeuver = getEphemeris(leo, mass, t0, dV, f, isp);
70 SmallManeuverAnalyticalModel model =
71 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp);
72 Assert.assertEquals(t0, model.getDate());
73
74 for (AbsoluteDate t = withoutManeuver.getMinDate();
75 t.compareTo(withoutManeuver.getMaxDate()) < 0;
76 t = t.shiftedBy(60.0)) {
77 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
78 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, leo.getFrame());
79 PVCoordinates pvModel = model.apply(withoutManeuver.propagate(t)).getPVCoordinates(leo.getFrame());
80 double nominalDeltaP = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
81 double modelError = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
82 if (t.compareTo(t0) < 0) {
83
84 Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
85 Assert.assertEquals(0, modelError, 1.0e-10);
86 } else {
87
88
89 if (t.durationFrom(t0) > 0.1 * leo.getKeplerianPeriod()) {
90 Assert.assertTrue(modelError < 0.009 * nominalDeltaP);
91 }
92 Assert.assertTrue(modelError < 0.8);
93 }
94 }
95
96 }
97
98 @Test
99 public void testLowEarthOrbit2() {
100
101 Orbit leo = new CircularOrbit(7200000.0, -1.0e-5, 2.0e-4,
102 FastMath.toRadians(98.0),
103 FastMath.toRadians(123.456),
104 0.0, PositionAngle.MEAN,
105 FramesFactory.getEME2000(),
106 new AbsoluteDate(new DateComponents(2004, 01, 01),
107 new TimeComponents(23, 30, 00.000),
108 TimeScalesFactory.getUTC()),
109 Constants.EIGEN5C_EARTH_MU);
110 double mass = 5600.0;
111 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
112 Vector3D dV = new Vector3D(-0.01, 0.02, 0.03);
113 double f = 20.0;
114 double isp = 315.0;
115 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, t0, Vector3D.ZERO, f, isp);
116 BoundedPropagator withManeuver = getEphemeris(leo, mass, t0, dV, f, isp);
117 SmallManeuverAnalyticalModel model =
118 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp);
119 Assert.assertEquals(t0, model.getDate());
120
121 for (AbsoluteDate t = withoutManeuver.getMinDate();
122 t.compareTo(withoutManeuver.getMaxDate()) < 0;
123 t = t.shiftedBy(60.0)) {
124 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
125 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, leo.getFrame());
126 PVCoordinates pvModel = model.apply(withoutManeuver.propagate(t).getOrbit()).getPVCoordinates(leo.getFrame());
127 double nominalDeltaP = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
128 double modelError = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
129 if (t.compareTo(t0) < 0) {
130
131 Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
132 Assert.assertEquals(0, modelError, 1.0e-10);
133 } else {
134
135
136 if (t.durationFrom(t0) > 0.1 * leo.getKeplerianPeriod()) {
137 Assert.assertTrue(modelError < 0.009 * nominalDeltaP);
138 }
139 Assert.assertTrue(modelError < 0.8);
140 }
141 }
142
143 }
144
145 @Test
146 public void testEccentricOrbit() {
147
148 Orbit heo = new KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0),
149 FastMath.toRadians(12.3456),
150 FastMath.toRadians(123.456),
151 FastMath.toRadians(1.23456), PositionAngle.MEAN,
152 FramesFactory.getEME2000(),
153 new AbsoluteDate(new DateComponents(2004, 01, 01),
154 new TimeComponents(23, 30, 00.000),
155 TimeScalesFactory.getUTC()),
156 Constants.EIGEN5C_EARTH_MU);
157 double mass = 5600.0;
158 AbsoluteDate t0 = heo.getDate().shiftedBy(1000.0);
159 Vector3D dV = new Vector3D(-0.01, 0.02, 0.03);
160 double f = 20.0;
161 double isp = 315.0;
162 BoundedPropagator withoutManeuver = getEphemeris(heo, mass, t0, Vector3D.ZERO, f, isp);
163 BoundedPropagator withManeuver = getEphemeris(heo, mass, t0, dV, f, isp);
164 SpacecraftState s0 = withoutManeuver.propagate(t0);
165 SmallManeuverAnalyticalModel model = new SmallManeuverAnalyticalModel(s0, dV, isp);
166 Assert.assertEquals(t0, model.getDate());
167 Assert.assertEquals(0.0,
168 Vector3D.distance(dV, s0.getAttitude().getRotation().applyTo(model.getInertialDV())),
169 1.0e-15);
170 Assert.assertSame(FramesFactory.getEME2000(), model.getInertialFrame());
171
172 for (AbsoluteDate t = withoutManeuver.getMinDate();
173 t.compareTo(withoutManeuver.getMaxDate()) < 0;
174 t = t.shiftedBy(600.0)) {
175 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame());
176 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, heo.getFrame());
177 PVCoordinates pvModel = model.apply(withoutManeuver.propagate(t)).getPVCoordinates(heo.getFrame());
178 double nominalDeltaP = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
179 double modelError = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
180 if (t.compareTo(t0) < 0) {
181
182 Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
183 Assert.assertEquals(0, modelError, 1.0e-10);
184 } else {
185
186
187 if (t.durationFrom(t0) > 0.01 * heo.getKeplerianPeriod()) {
188 Assert.assertTrue(modelError < 0.005 * nominalDeltaP);
189 }
190 Assert.assertTrue(modelError < 1700);
191 }
192 }
193
194 }
195
196 @Test
197 public void testJacobian() {
198
199 Frame eme2000 = FramesFactory.getEME2000();
200 Orbit leo = new CircularOrbit(7200000.0, -1.0e-2, 2.0e-3,
201 FastMath.toRadians(98.0),
202 FastMath.toRadians(123.456),
203 0.3, PositionAngle.MEAN,
204 eme2000,
205 new AbsoluteDate(new DateComponents(2004, 01, 01),
206 new TimeComponents(23, 30, 00.000),
207 TimeScalesFactory.getUTC()),
208 Constants.EIGEN5C_EARTH_MU);
209 double mass = 5600.0;
210 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
211 Vector3D dV0 = new Vector3D(-0.1, 0.2, 0.3);
212 double f = 400.0;
213 double isp = 315.0;
214
215 for (OrbitType orbitType : OrbitType.values()) {
216 for (PositionAngle positionAngle : PositionAngle.values()) {
217 BoundedPropagator withoutManeuver = getEphemeris(orbitType.convertType(leo), mass, t0, Vector3D.ZERO, f, isp);
218
219 SpacecraftState state0 = withoutManeuver.propagate(t0);
220 SmallManeuverAnalyticalModel model = new SmallManeuverAnalyticalModel(state0, eme2000, dV0, isp);
221 Assert.assertEquals(t0, model.getDate());
222
223 Vector3D[] velDirs = new Vector3D[] { Vector3D.PLUS_I, Vector3D.PLUS_J, Vector3D.PLUS_K, Vector3D.ZERO };
224 double[] timeDirs = new double[] { 0, 0, 0, 1 };
225 double h = 1.0;
226 AbsoluteDate t1 = t0.shiftedBy(20.0);
227 for (int i = 0; i < 4; ++i) {
228
229 SmallManeuverAnalyticalModel[] models = new SmallManeuverAnalyticalModel[] {
230 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-4 * h * timeDirs[i])),
231 eme2000, new Vector3D(1, dV0, -4 * h, velDirs[i]), isp),
232 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-3 * h * timeDirs[i])),
233 eme2000, new Vector3D(1, dV0, -3 * h, velDirs[i]), isp),
234 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-2 * h * timeDirs[i])),
235 eme2000, new Vector3D(1, dV0, -2 * h, velDirs[i]), isp),
236 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-1 * h * timeDirs[i])),
237 eme2000, new Vector3D(1, dV0, -1 * h, velDirs[i]), isp),
238 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+1 * h * timeDirs[i])),
239 eme2000, new Vector3D(1, dV0, +1 * h, velDirs[i]), isp),
240 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+2 * h * timeDirs[i])),
241 eme2000, new Vector3D(1, dV0, +2 * h, velDirs[i]), isp),
242 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+3 * h * timeDirs[i])),
243 eme2000, new Vector3D(1, dV0, +3 * h, velDirs[i]), isp),
244 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+4 * h * timeDirs[i])),
245 eme2000, new Vector3D(1, dV0, +4 * h, velDirs[i]), isp),
246 };
247 double[][] array = new double[models.length][6];
248
249 Orbit orbitWithout = withoutManeuver.propagate(t1).getOrbit();
250
251
252 double c = 1.0 / (840 * h);
253 for (int j = 0; j < models.length; ++j) {
254 orbitType.mapOrbitToArray(models[j].apply(orbitWithout), positionAngle, array[j], null);
255 }
256 double[] orbitGradient = new double[6];
257 for (int k = 0; k < orbitGradient.length; ++k) {
258 double d4 = array[7][k] - array[0][k];
259 double d3 = array[6][k] - array[1][k];
260 double d2 = array[5][k] - array[2][k];
261 double d1 = array[4][k] - array[3][k];
262 orbitGradient[k] = (-3 * d4 + 32 * d3 - 168 * d2 + 672 * d1) * c;
263 }
264
265
266 double[][] jacobian = new double[6][4];
267 model.getJacobian(orbitWithout, positionAngle, jacobian);
268
269 for (int j = 0; j < orbitGradient.length; ++j) {
270 Assert.assertEquals(orbitGradient[j], jacobian[j][i], 1.6e-4 * FastMath.abs(orbitGradient[j]));
271 }
272
273 }
274
275 }
276
277 }
278
279 }
280
281 private BoundedPropagator getEphemeris(final Orbit orbit, final double mass,
282 final AbsoluteDate t0, final Vector3D dV,
283 final double f, final double isp)
284 {
285
286 AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH);
287 final SpacecraftState initialState =
288 new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
289
290
291
292 final double dP = 1.0;
293 double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
294 AdaptiveStepsizeIntegrator integrator =
295 new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
296 integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
297 final NumericalPropagator propagator = new NumericalPropagator(integrator);
298 propagator.setOrbitType(orbit.getType());
299 propagator.setInitialState(initialState);
300 propagator.setAttitudeProvider(law);
301
302 if (dV.getNorm() > 1.0e-6) {
303
304 final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
305 final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
306 final ConstantThrustManeuver maneuver =
307 new ConstantThrustManeuver(t0, dt , f, isp, dV.normalize());
308 propagator.addForceModel(maneuver);
309 }
310
311 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
312 propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod()));
313 return generator.getGeneratedEphemeris();
314
315 }
316
317 @Before
318 public void setUp() {
319 Utils.setDataRoot("regular-data");
320 }
321
322 }