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