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