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.errors.OrekitException;
31 import org.orekit.frames.Frame;
32 import org.orekit.frames.FramesFactory;
33 import org.orekit.frames.LOFType;
34 import org.orekit.orbits.CircularOrbit;
35 import org.orekit.orbits.KeplerianOrbit;
36 import org.orekit.orbits.Orbit;
37 import org.orekit.orbits.OrbitType;
38 import org.orekit.orbits.PositionAngle;
39 import org.orekit.propagation.BoundedPropagator;
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() throws OrekitException {
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() throws OrekitException {
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() throws OrekitException {
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 SmallManeuverAnalyticalModel model =
165 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0), dV, isp);
166 Assert.assertEquals(t0, model.getDate());
167
168 for (AbsoluteDate t = withoutManeuver.getMinDate();
169 t.compareTo(withoutManeuver.getMaxDate()) < 0;
170 t = t.shiftedBy(600.0)) {
171 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame());
172 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, heo.getFrame());
173 PVCoordinates pvModel = model.apply(withoutManeuver.propagate(t)).getPVCoordinates(heo.getFrame());
174 double nominalDeltaP = new PVCoordinates(pvWith, pvWithout).getPosition().getNorm();
175 double modelError = new PVCoordinates(pvWith, pvModel).getPosition().getNorm();
176 if (t.compareTo(t0) < 0) {
177
178 Assert.assertEquals(0, nominalDeltaP, 1.0e-10);
179 Assert.assertEquals(0, modelError, 1.0e-10);
180 } else {
181
182
183 if (t.durationFrom(t0) > 0.01 * heo.getKeplerianPeriod()) {
184 Assert.assertTrue(modelError < 0.005 * nominalDeltaP);
185 }
186 Assert.assertTrue(modelError < 1700);
187 }
188 }
189
190 }
191
192 @Test
193 public void testJacobian() throws OrekitException {
194
195 Frame eme2000 = FramesFactory.getEME2000();
196 Orbit leo = new CircularOrbit(7200000.0, -1.0e-2, 2.0e-3,
197 FastMath.toRadians(98.0),
198 FastMath.toRadians(123.456),
199 0.3, PositionAngle.MEAN,
200 eme2000,
201 new AbsoluteDate(new DateComponents(2004, 01, 01),
202 new TimeComponents(23, 30, 00.000),
203 TimeScalesFactory.getUTC()),
204 Constants.EIGEN5C_EARTH_MU);
205 double mass = 5600.0;
206 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
207 Vector3D dV0 = new Vector3D(-0.1, 0.2, 0.3);
208 double f = 400.0;
209 double isp = 315.0;
210
211 for (OrbitType orbitType : OrbitType.values()) {
212 for (PositionAngle positionAngle : PositionAngle.values()) {
213 BoundedPropagator withoutManeuver = getEphemeris(orbitType.convertType(leo), mass, t0, Vector3D.ZERO, f, isp);
214
215 SpacecraftState state0 = withoutManeuver.propagate(t0);
216 SmallManeuverAnalyticalModel model = new SmallManeuverAnalyticalModel(state0, eme2000, dV0, isp);
217 Assert.assertEquals(t0, model.getDate());
218
219 Vector3D[] velDirs = new Vector3D[] { Vector3D.PLUS_I, Vector3D.PLUS_J, Vector3D.PLUS_K, Vector3D.ZERO };
220 double[] timeDirs = new double[] { 0, 0, 0, 1 };
221 double h = 1.0;
222 AbsoluteDate t1 = t0.shiftedBy(20.0);
223 for (int i = 0; i < 4; ++i) {
224
225 SmallManeuverAnalyticalModel[] models = new SmallManeuverAnalyticalModel[] {
226 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-4 * h * timeDirs[i])),
227 eme2000, new Vector3D(1, dV0, -4 * h, velDirs[i]), isp),
228 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-3 * h * timeDirs[i])),
229 eme2000, new Vector3D(1, dV0, -3 * h, velDirs[i]), isp),
230 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-2 * h * timeDirs[i])),
231 eme2000, new Vector3D(1, dV0, -2 * h, velDirs[i]), isp),
232 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-1 * h * timeDirs[i])),
233 eme2000, new Vector3D(1, dV0, -1 * h, velDirs[i]), isp),
234 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+1 * h * timeDirs[i])),
235 eme2000, new Vector3D(1, dV0, +1 * h, velDirs[i]), isp),
236 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+2 * h * timeDirs[i])),
237 eme2000, new Vector3D(1, dV0, +2 * h, velDirs[i]), isp),
238 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+3 * h * timeDirs[i])),
239 eme2000, new Vector3D(1, dV0, +3 * h, velDirs[i]), isp),
240 new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+4 * h * timeDirs[i])),
241 eme2000, new Vector3D(1, dV0, +4 * h, velDirs[i]), isp),
242 };
243 double[][] array = new double[models.length][6];
244
245 Orbit orbitWithout = withoutManeuver.propagate(t1).getOrbit();
246
247
248 double c = 1.0 / (840 * h);
249 for (int j = 0; j < models.length; ++j) {
250 orbitType.mapOrbitToArray(models[j].apply(orbitWithout), positionAngle, array[j], null);
251 }
252 double[] orbitGradient = new double[6];
253 for (int k = 0; k < orbitGradient.length; ++k) {
254 double d4 = array[7][k] - array[0][k];
255 double d3 = array[6][k] - array[1][k];
256 double d2 = array[5][k] - array[2][k];
257 double d1 = array[4][k] - array[3][k];
258 orbitGradient[k] = (-3 * d4 + 32 * d3 - 168 * d2 + 672 * d1) * c;
259 }
260
261
262 double[][] jacobian = new double[6][4];
263 model.getJacobian(orbitWithout, positionAngle, jacobian);
264
265 for (int j = 0; j < orbitGradient.length; ++j) {
266 Assert.assertEquals(orbitGradient[j], jacobian[j][i], 1.6e-4 * FastMath.abs(orbitGradient[j]));
267 }
268
269 }
270
271 }
272
273 }
274
275 }
276
277 private BoundedPropagator getEphemeris(final Orbit orbit, final double mass,
278 final AbsoluteDate t0, final Vector3D dV,
279 final double f, final double isp)
280 throws OrekitException {
281
282 AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH);
283 final SpacecraftState initialState =
284 new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
285
286
287
288 final double dP = 1.0;
289 double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
290 AdaptiveStepsizeIntegrator integrator =
291 new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
292 integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
293 final NumericalPropagator propagator = new NumericalPropagator(integrator);
294 propagator.setOrbitType(orbit.getType());
295 propagator.setInitialState(initialState);
296 propagator.setAttitudeProvider(law);
297
298 if (dV.getNorm() > 1.0e-6) {
299
300 final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
301 final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
302 final ConstantThrustManeuver maneuver =
303 new ConstantThrustManeuver(t0, dt , f, isp, dV.normalize());
304 propagator.addForceModel(maneuver);
305 }
306
307 propagator.setEphemerisMode();
308 propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod()));
309 return propagator.getGeneratedEphemeris();
310
311 }
312
313 @Before
314 public void setUp() {
315 Utils.setDataRoot("regular-data");
316 }
317
318 }