1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.analytical;
18
19
20 import java.io.IOException;
21 import java.text.ParseException;
22
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
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.AttitudeProvider;
32 import org.orekit.attitudes.LofOffset;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
35 import org.orekit.forces.gravity.ThirdBodyAttraction;
36 import org.orekit.forces.gravity.potential.GravityFieldFactory;
37 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
38 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
39 import org.orekit.forces.maneuvers.ConstantThrustManeuver;
40 import org.orekit.forces.maneuvers.SmallManeuverAnalyticalModel;
41 import org.orekit.frames.FramesFactory;
42 import org.orekit.frames.LOFType;
43 import org.orekit.orbits.CircularOrbit;
44 import org.orekit.orbits.KeplerianOrbit;
45 import org.orekit.orbits.Orbit;
46 import org.orekit.orbits.PositionAngle;
47 import org.orekit.propagation.AdditionalStateProvider;
48 import org.orekit.propagation.BoundedPropagator;
49 import org.orekit.propagation.EphemerisGenerator;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.propagation.numerical.NumericalPropagator;
52 import org.orekit.time.AbsoluteDate;
53 import org.orekit.time.DateComponents;
54 import org.orekit.time.TimeComponents;
55 import org.orekit.time.TimeScalesFactory;
56 import org.orekit.utils.Constants;
57 import org.orekit.utils.PVCoordinates;
58
59 public class AdapterPropagatorTest {
60
61 @Test
62 public void testLowEarthOrbit() throws ParseException, IOException {
63
64 Orbit leo = new CircularOrbit(7200000.0, -1.0e-5, 2.0e-4,
65 FastMath.toRadians(98.0),
66 FastMath.toRadians(123.456),
67 0.0, PositionAngle.MEAN,
68 FramesFactory.getEME2000(),
69 new AbsoluteDate(new DateComponents(2004, 01, 01),
70 new TimeComponents(23, 30, 00.000),
71 TimeScalesFactory.getUTC()),
72 Constants.EIGEN5C_EARTH_MU);
73 double mass = 5600.0;
74 AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
75 Vector3D dV = new Vector3D(-0.1, 0.2, 0.3);
76 double f = 20.0;
77 double isp = 315.0;
78 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
79 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
80 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, 5,
81 new LofOffset(leo.getFrame(), LOFType.LVLH),
82 t0, Vector3D.ZERO, f, isp,
83 false, false, null);
84 BoundedPropagator withManeuver = getEphemeris(leo, mass, 5,
85 new LofOffset(leo.getFrame(), LOFType.LVLH),
86 t0, dV, f, isp,
87 false, false, null);
88
89
90 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
91 AdapterPropagator.DifferentialEffect effect =
92 new SmallManeuverAnalyticalModel(adapterPropagator.propagate(t0), dV.negate(), isp);
93 adapterPropagator.addEffect(effect);
94 adapterPropagator.addAdditionalStateProvider(new AdditionalStateProvider() {
95 public String getName() {
96 return "dummy 3";
97 }
98 public double[] getAdditionalState(SpacecraftState state) {
99 return new double[3];
100 }
101 });
102
103
104
105 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 1"));
106 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 2"));
107 Assert.assertTrue(adapterPropagator.isAdditionalStateManaged("dummy 3"));
108
109 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
110 t.compareTo(withoutManeuver.getMaxDate()) < 0;
111 t = t.shiftedBy(60.0)) {
112 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
113 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, leo.getFrame());
114 double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm();
115 Assert.assertEquals(0, revertError, 0.45);
116 Assert.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
117 Assert.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
118 Assert.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
119 }
120
121 }
122
123 @Test
124 public void testEccentricOrbit() throws ParseException, IOException {
125
126 Orbit heo = new KeplerianOrbit(90000000.0, 0.92, FastMath.toRadians(98.0),
127 FastMath.toRadians(12.3456),
128 FastMath.toRadians(123.456),
129 FastMath.toRadians(1.23456), PositionAngle.MEAN,
130 FramesFactory.getEME2000(),
131 new AbsoluteDate(new DateComponents(2004, 01, 01),
132 new TimeComponents(23, 30, 00.000),
133 TimeScalesFactory.getUTC()),
134 Constants.EIGEN5C_EARTH_MU);
135 double mass = 5600.0;
136 AbsoluteDate t0 = heo.getDate().shiftedBy(1000.0);
137 Vector3D dV = new Vector3D(-0.01, 0.02, 0.03);
138 double f = 20.0;
139 double isp = 315.0;
140 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
141 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
142 BoundedPropagator withoutManeuver = getEphemeris(heo, mass, 5,
143 new LofOffset(heo.getFrame(), LOFType.LVLH),
144 t0, Vector3D.ZERO, f, isp,
145 false, false, null);
146 BoundedPropagator withManeuver = getEphemeris(heo, mass, 5,
147 new LofOffset(heo.getFrame(), LOFType.LVLH),
148 t0, dV, f, isp,
149 false, false, null);
150
151
152 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
153 AdapterPropagator.DifferentialEffect effect =
154 new SmallManeuverAnalyticalModel(adapterPropagator.propagate(t0), dV.negate(), isp);
155 adapterPropagator.addEffect(effect);
156 adapterPropagator.addAdditionalStateProvider(new AdditionalStateProvider() {
157 public String getName() {
158 return "dummy 3";
159 }
160 public double[] getAdditionalState(SpacecraftState state) {
161 return new double[3];
162 }
163 });
164
165
166
167 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 1"));
168 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 2"));
169 Assert.assertTrue(adapterPropagator.isAdditionalStateManaged("dummy 3"));
170
171 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
172 t.compareTo(withoutManeuver.getMaxDate()) < 0;
173 t = t.shiftedBy(300.0)) {
174 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, heo.getFrame());
175 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, heo.getFrame());
176 double revertError = Vector3D.distance(pvWithout.getPosition(), pvReverted.getPosition());
177 Assert.assertEquals(0, revertError, 2.5e-5 * heo.getA());
178 Assert.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
179 Assert.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
180 Assert.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
181 }
182
183 }
184
185 @Test
186 public void testNonKeplerian() throws ParseException, IOException {
187
188 Orbit leo = new CircularOrbit(7204319.233600575, 4.434564637450575E-4, 0.0011736728299091088,
189 1.7211611441767323, 5.5552084166959474,
190 24950.321259193086, PositionAngle.TRUE,
191 FramesFactory.getEME2000(),
192 new AbsoluteDate(new DateComponents(2003, 9, 16),
193 new TimeComponents(23, 11, 20.264),
194 TimeScalesFactory.getUTC()),
195 Constants.EIGEN5C_EARTH_MU);
196 double mass = 4093.0;
197 AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2003, 9, 16),
198 new TimeComponents(23, 14, 40.264),
199 TimeScalesFactory.getUTC());
200 Vector3D dV = new Vector3D(0.0, 3.0, 0.0);
201 double f = 40.0;
202 double isp = 300.0;
203 double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
204 double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
205
206
207 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("g007_eigen_05c_coef", false));
208 NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getNormalizedProvider(8, 8);
209 BoundedPropagator withoutManeuver = getEphemeris(leo, mass, 10,
210 new LofOffset(leo.getFrame(), LOFType.VNC),
211 t0, Vector3D.ZERO, f, isp,
212 true, true, gravityField);
213 BoundedPropagator withManeuver = getEphemeris(leo, mass, 10,
214 new LofOffset(leo.getFrame(), LOFType.VNC),
215 t0, dV, f, isp,
216 true, true, gravityField);
217
218
219 AdapterPropagator adapterPropagator = new AdapterPropagator(withManeuver);
220 SpacecraftState state0 = adapterPropagator.propagate(t0);
221 AdapterPropagator.DifferentialEffect directEffect =
222 new SmallManeuverAnalyticalModel(state0, dV.negate(), isp);
223 AdapterPropagator.DifferentialEffect derivedEffect =
224 new J2DifferentialEffect(state0, directEffect, false,
225 GravityFieldFactory.getUnnormalizedProvider(gravityField));
226 adapterPropagator.addEffect(directEffect);
227 adapterPropagator.addEffect(derivedEffect);
228 adapterPropagator.addAdditionalStateProvider(new AdditionalStateProvider() {
229 public String getName() {
230 return "dummy 3";
231 }
232 public double[] getAdditionalState(SpacecraftState state) {
233 return new double[3];
234 }
235 });
236
237
238
239 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 1"));
240 Assert.assertFalse(adapterPropagator.isAdditionalStateManaged("dummy 2"));
241 Assert.assertTrue(adapterPropagator.isAdditionalStateManaged("dummy 3"));
242
243 double maxDelta = 0;
244 double maxNominal = 0;
245 for (AbsoluteDate t = t0.shiftedBy(0.5 * dt);
246 t.compareTo(withoutManeuver.getMaxDate()) < 0;
247 t = t.shiftedBy(600.0)) {
248 PVCoordinates pvWithout = withoutManeuver.getPVCoordinates(t, leo.getFrame());
249 PVCoordinates pvWith = withManeuver.getPVCoordinates(t, leo.getFrame());
250 PVCoordinates pvReverted = adapterPropagator.getPVCoordinates(t, leo.getFrame());
251 double nominal = new PVCoordinates(pvWithout, pvWith).getPosition().getNorm();
252 double revertError = new PVCoordinates(pvWithout, pvReverted).getPosition().getNorm();
253 maxDelta = FastMath.max(maxDelta, revertError);
254 maxNominal = FastMath.max(maxNominal, nominal);
255 Assert.assertEquals(2, adapterPropagator.propagate(t).getAdditionalState("dummy 1").length);
256 Assert.assertEquals(1, adapterPropagator.propagate(t).getAdditionalState("dummy 2").length);
257 Assert.assertEquals(3, adapterPropagator.propagate(t).getAdditionalState("dummy 3").length);
258 }
259 Assert.assertTrue(maxDelta < 120);
260 Assert.assertTrue(maxNominal > 2800);
261
262 }
263
264 private BoundedPropagator getEphemeris(final Orbit orbit, final double mass, final int nbOrbits,
265 final AttitudeProvider law,
266 final AbsoluteDate t0, final Vector3D dV,
267 final double f, final double isp,
268 final boolean sunAttraction, final boolean moonAttraction,
269 final NormalizedSphericalHarmonicsProvider gravityField)
270 throws ParseException, IOException {
271
272 SpacecraftState initialState =
273 new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
274
275
276 initialState = initialState.addAdditionalState("dummy 1", 1.25, 2.5);
277 initialState = initialState.addAdditionalState("dummy 2", 5.0);
278
279
280 final double dP = 1.0;
281 double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
282 AdaptiveStepsizeIntegrator integrator =
283 new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
284 integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
285 final NumericalPropagator propagator = new NumericalPropagator(integrator);
286 propagator.addAdditionalStateProvider(new AdditionalStateProvider() {
287 public String getName() {
288 return "dummy 2";
289 }
290 public double[] getAdditionalState(SpacecraftState state) {
291 return new double[] { 5.0 };
292 }
293 });
294 propagator.setInitialState(initialState);
295 propagator.setAttitudeProvider(law);
296
297 if (dV.getNorm() > 1.0e-6) {
298
299 final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
300 final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
301 final ConstantThrustManeuver maneuver =
302 new ConstantThrustManeuver(t0.shiftedBy(-0.5 * dt), dt , f, isp, dV.normalize());
303 propagator.addForceModel(maneuver);
304 }
305
306 if (sunAttraction) {
307 propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
308 }
309
310 if (moonAttraction) {
311 propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
312 }
313
314 if (gravityField != null) {
315 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getGTOD(false),
316 gravityField));
317 }
318
319 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
320 propagator.propagate(t0.shiftedBy(nbOrbits * orbit.getKeplerianPeriod()));
321
322 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
323
324
325
326 Assert.assertFalse(propagator.isAdditionalStateManaged("dummy 1"));
327 Assert.assertTrue(propagator.isAdditionalStateManaged("dummy 2"));
328 Assert.assertFalse(ephemeris.isAdditionalStateManaged("dummy 1"));
329 Assert.assertTrue(ephemeris.isAdditionalStateManaged("dummy 2"));
330 Assert.assertEquals(2, ephemeris.getInitialState().getAdditionalState("dummy 1").length);
331 Assert.assertEquals(1, ephemeris.getInitialState().getAdditionalState("dummy 2").length);
332
333 return ephemeris;
334
335 }
336
337 @Before
338 public void setUp() {
339 Utils.setDataRoot("regular-data:potential/icgem-format");
340 }
341
342 }