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.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.List;
23
24 import org.hipparchus.exception.MathIllegalArgumentException;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.util.FastMath;
27 import org.junit.After;
28 import org.junit.Assert;
29 import org.junit.Before;
30 import org.junit.Test;
31 import org.orekit.Utils;
32 import org.orekit.bodies.CelestialBodyFactory;
33 import org.orekit.errors.OrekitException;
34 import org.orekit.errors.OrekitMessages;
35 import org.orekit.frames.Frame;
36 import org.orekit.frames.FramesFactory;
37 import org.orekit.orbits.CartesianOrbit;
38 import org.orekit.orbits.EquinoctialOrbit;
39 import org.orekit.orbits.KeplerianOrbit;
40 import org.orekit.orbits.Orbit;
41 import org.orekit.orbits.PositionAngle;
42 import org.orekit.propagation.AdditionalStateProvider;
43 import org.orekit.propagation.SpacecraftState;
44 import org.orekit.time.AbsoluteDate;
45 import org.orekit.time.DateComponents;
46 import org.orekit.time.TimeComponents;
47 import org.orekit.time.TimeScale;
48 import org.orekit.time.TimeScalesFactory;
49 import org.orekit.utils.PVCoordinates;
50 import org.orekit.utils.TimeStampedPVCoordinates;
51
52
53 public class TabulatedEphemerisTest {
54
55 @Test
56 public void testInterpolationFullEcksteinHechlerOrbit() {
57
58
59 checkInterpolation(new StateFilter() {
60 public SpacecraftState filter(SpacecraftState state) {
61 return state;
62 }
63 }, 6.62e-4, 1.89e-5);
64 }
65
66 @Test
67 public void testInterpolationKeplerianAcceleration() {
68
69 checkInterpolation(new StateFilter() {
70 public SpacecraftState filter(SpacecraftState state) {
71 CartesianOrbit c = (CartesianOrbit) state.getOrbit();
72 Vector3D p = c.getPVCoordinates().getPosition();
73 Vector3D v = c.getPVCoordinates().getVelocity();
74 double r2 = p.getNormSq();
75 double r = FastMath.sqrt(r2);
76 Vector3D kepA = new Vector3D(-c.getMu() / (r * r2), c.getPVCoordinates().getPosition());
77 return new SpacecraftState(new CartesianOrbit(new TimeStampedPVCoordinates(c.getDate(),
78 p, v, kepA),
79 c.getFrame(), c.getMu()),
80 state.getAttitude(),
81 state.getMass(),
82 state.getAdditionalStatesValues());
83 }
84 }, 8.5, 0.22);
85 }
86
87 private void checkInterpolation(StateFilter f, double expectedDP, double expectedDV)
88 {
89
90 double mass = 2500;
91 double a = 7187990.1979844316;
92 double e = 0.5e-4;
93 double i = 1.7105407051081795;
94 double omega = 1.9674147913622104;
95 double OMEGA = FastMath.toRadians(261);
96 double lv = 0;
97
98 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
99 TimeComponents.H00,
100 TimeScalesFactory.getUTC());
101 final AbsoluteDate finalDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
102 TimeComponents.H00,
103 TimeScalesFactory.getUTC());
104 double deltaT = finalDate.durationFrom(initDate);
105
106 Orbit transPar = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE,
107 FramesFactory.getEME2000(), initDate, mu);
108
109 int nbIntervals = 720;
110 EcksteinHechlerPropagator eck =
111 new EcksteinHechlerPropagator(transPar, mass,
112 ae, mu, c20, c30, c40, c50, c60);
113 AdditionalStateProvider provider = new AdditionalStateProvider() {
114
115 public String getName() {
116 return "dt";
117 }
118
119 public double[] getAdditionalState(SpacecraftState state) {
120 return new double[] { state.getDate().durationFrom(initDate) };
121 }
122 };
123 eck.addAdditionalStateProvider(provider);
124 try {
125 eck.addAdditionalStateProvider(provider);
126 Assert.fail("an exception should have been thrown");
127 } catch (OrekitException oe) {
128 Assert.assertEquals(OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE,
129 oe.getSpecifier());
130 }
131 List<SpacecraftState> tab = new ArrayList<SpacecraftState>(nbIntervals + 1);
132 for (int j = 0; j<= nbIntervals; j++) {
133 AbsoluteDate current = initDate.shiftedBy((j * deltaT) / nbIntervals);
134 tab.add(f.filter(eck.propagate(current)));
135 }
136
137 try {
138 new Ephemeris(tab, nbIntervals + 2);
139 Assert.fail("an exception should have been thrown");
140 } catch (MathIllegalArgumentException miae) {
141
142 }
143 Ephemeris te = new Ephemeris(tab, 2);
144
145 Assert.assertEquals(0.0, te.getMaxDate().durationFrom(finalDate), 1.0e-9);
146 Assert.assertEquals(0.0, te.getMinDate().durationFrom(initDate), 1.0e-9);
147
148 double maxP = 0;
149 double maxV = 0;
150 for (double dt = 0; dt < 3600; dt += 1) {
151 AbsoluteDate date = initDate.shiftedBy(dt);
152 CartesianOrbit c1 = (CartesianOrbit) eck.propagate(date).getOrbit();
153 CartesianOrbit c2 = (CartesianOrbit) te.propagate(date).getOrbit();
154 maxP = FastMath.max(maxP,
155 Vector3D.distance(c1.getPVCoordinates().getPosition(),
156 c2.getPVCoordinates().getPosition()));
157 maxV = FastMath.max(maxV,
158 Vector3D.distance(c1.getPVCoordinates().getVelocity(),
159 c2.getPVCoordinates().getVelocity()));
160 }
161
162 Assert.assertEquals(expectedDP, maxP, 0.1 * expectedDP);
163 Assert.assertEquals(expectedDV, maxV, 0.1 * expectedDV);
164
165 }
166
167 private interface StateFilter {
168 SpacecraftState filter(final SpacecraftState state);
169 }
170
171 @Test
172 public void testPiWraping() {
173
174 TimeScale utc= TimeScalesFactory.getUTC();
175 Frame frame = FramesFactory.getEME2000();
176 double mu=CelestialBodyFactory.getEarth().getGM();
177 AbsoluteDate t0 = new AbsoluteDate(2009, 10, 29, 0, 0, 0, utc);
178
179 AbsoluteDate t1 = new AbsoluteDate(t0, 1320.0);
180 Vector3D p1 = new Vector3D(-0.17831296727974E+08, 0.67919502669856E+06, -0.16591008368477E+07);
181 Vector3D v1 = new Vector3D(-0.38699705630724E+04, -0.36209408682762E+04, -0.16255053872347E+03);
182 SpacecraftState s1 = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(p1, v1), frame, t1, mu));
183
184 AbsoluteDate t2 = new AbsoluteDate(t0, 1440.0);
185 Vector3D p2 = new Vector3D(-0.18286942572033E+08, 0.24442124296930E+06, -0.16777961761695E+07);
186 Vector3D v2 = new Vector3D(-0.37252897467918E+04, -0.36246628128896E+04, -0.14917724596280E+03);
187 SpacecraftState s2 = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(p2, v2), frame, t2, mu));
188
189 AbsoluteDate t3 = new AbsoluteDate(t0, 1560.0);
190 Vector3D p3 = new Vector3D(-0.18725635245837E+08, -0.19058407701834E+06, -0.16949352249614E+07);
191 Vector3D v3 = new Vector3D(-0.35873348682393E+04, -0.36248828501784E+04, -0.13660045394149E+03);
192 SpacecraftState s3 = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(p3, v3), frame, t3, mu));
193
194 Ephemeris ephem= new Ephemeris(Arrays.asList(s1, s2, s3), 2);
195
196 AbsoluteDate tA = new AbsoluteDate(t0, 24 * 60);
197 Vector3D pA = ephem.propagate(tA).getPVCoordinates(frame).getPosition();
198 Assert.assertEquals(1.766,
199 Vector3D.distance(pA, s1.shiftedBy(tA.durationFrom(s1.getDate())).getPVCoordinates(frame).getPosition()),
200 1.0e-3);
201 Assert.assertEquals(0.000,
202 Vector3D.distance(pA, s2.shiftedBy(tA.durationFrom(s2.getDate())).getPVCoordinates(frame).getPosition()),
203 1.0e-3);
204 Assert.assertEquals(1.556,
205 Vector3D.distance(pA, s3.shiftedBy(tA.durationFrom(s3.getDate())).getPVCoordinates(frame).getPosition()),
206 1.0e-3);
207
208 AbsoluteDate tB = new AbsoluteDate(t0, 25 * 60);
209 Vector3D pB = ephem.propagate(tB).getPVCoordinates(frame).getPosition();
210 Assert.assertEquals(2.646,
211 Vector3D.distance(pB, s1.shiftedBy(tB.durationFrom(s1.getDate())).getPVCoordinates(frame).getPosition()),
212 1.0e-3);
213 Assert.assertEquals(2.619,
214 Vector3D.distance(pB, s2.shiftedBy(tB.durationFrom(s2.getDate())).getPVCoordinates(frame).getPosition()),
215 1.0e-3);
216 Assert.assertEquals(2.632,
217 Vector3D.distance(pB, s3.shiftedBy(tB.durationFrom(s3.getDate())).getPVCoordinates(frame).getPosition()),
218 1.0e-3);
219
220 AbsoluteDate tC = new AbsoluteDate(t0, 26 * 60);
221 Vector3D pC = ephem.propagate(tC).getPVCoordinates(frame).getPosition();
222 Assert.assertEquals(6.851,
223 Vector3D.distance(pC, s1.shiftedBy(tC.durationFrom(s1.getDate())).getPVCoordinates(frame).getPosition()),
224 1.0e-3);
225 Assert.assertEquals(1.605,
226 Vector3D.distance(pC, s2.shiftedBy(tC.durationFrom(s2.getDate())).getPVCoordinates(frame).getPosition()),
227 1.0e-3);
228 Assert.assertEquals(0.000,
229 Vector3D.distance(pC, s3.shiftedBy(tC.durationFrom(s3.getDate())).getPVCoordinates(frame).getPosition()),
230 1.0e-3);
231
232 }
233
234
235 @Test
236 public void testGetFrame() throws MathIllegalArgumentException, IllegalArgumentException, OrekitException {
237
238 Frame frame = FramesFactory.getICRF();
239 AbsoluteDate date = AbsoluteDate.JULIAN_EPOCH;
240
241 SpacecraftState state = new SpacecraftState(
242 new KeplerianOrbit(1e9, 0.01, 1, 1, 1, 1, PositionAngle.TRUE, frame, date, mu));
243 Ephemeris ephem = new Ephemeris(Arrays.asList(state, state.shiftedBy(1)), 2);
244
245
246 Assert.assertSame(ephem.getFrame(), frame);
247 }
248
249 @Before
250 public void setUp() {
251 Utils.setDataRoot("regular-data");
252 mu = 3.9860047e14;
253 ae = 6.378137e6;
254 c20 = -1.08263e-3;
255 c30 = 2.54e-6;
256 c40 = 1.62e-6;
257 c50 = 2.3e-7;
258 c60 = -5.5e-7;
259 }
260
261 @After
262 public void tearDown() {
263 mu = Double.NaN;
264 ae = Double.NaN;
265 c20 = Double.NaN;
266 c30 = Double.NaN;
267 c40 = Double.NaN;
268 c50 = Double.NaN;
269 c60 = Double.NaN;
270 }
271
272 private double mu;
273 private double ae;
274 private double c20;
275 private double c30;
276 private double c40;
277 private double c50;
278 private double c60;
279
280 }