1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.semianalytical.dsst.forces;
18
19 import java.io.IOException;
20 import java.lang.reflect.InvocationTargetException;
21 import java.lang.reflect.Method;
22 import java.text.ParseException;
23 import java.util.ArrayList;
24 import java.util.Arrays;
25 import java.util.List;
26
27 import org.hipparchus.exception.LocalizedCoreFormats;
28 import org.hipparchus.util.FastMath;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.BeforeEach;
31 import org.junit.jupiter.api.Test;
32 import org.orekit.Utils;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.bodies.OneAxisEllipsoid;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.forces.gravity.potential.GravityFieldFactory;
37 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
38 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
39 import org.orekit.frames.Frame;
40 import org.orekit.frames.FramesFactory;
41 import org.orekit.orbits.EquinoctialOrbit;
42 import org.orekit.orbits.KeplerianOrbit;
43 import org.orekit.orbits.Orbit;
44 import org.orekit.orbits.PositionAngleType;
45 import org.orekit.propagation.PropagationType;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.TimeScalesFactory;
50 import org.orekit.utils.Constants;
51
52 class DSSTTesseralTest {
53
54 @Test
55 void testGetMeanElementRate() {
56
57
58 final UnnormalizedSphericalHarmonicsProvider provider =
59 GravityFieldFactory.getUnnormalizedProvider(4, 4);
60
61 final Frame frame = FramesFactory.getEME2000();
62 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
63 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
64
65
66
67
68
69
70
71 final Orbit orbit = new EquinoctialOrbit(2.655989E7,
72 2.719455286199036E-4,
73 0.0041543085910249414,
74 -0.3412974060023717,
75 0.3960084733107685,
76 8.566537840341699,
77 PositionAngleType.TRUE,
78 frame,
79 initDate,
80 provider.getMu());
81
82 final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
83
84 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
85
86 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
87 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
88 4, 4, 4, 8, 4, 4, 2);
89
90 final double[] parameters = tesseral.getParameters(orbit.getDate());
91
92
93 tesseral.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, parameters);
94
95 final double[] elements = new double[7];
96 Arrays.fill(elements, 0.0);
97
98 final double[] daidt = tesseral.getMeanElementRate(state, auxiliaryElements, parameters);
99 for (int i = 0; i < daidt.length; i++) {
100 elements[i] = daidt[i];
101 }
102 Assertions.assertEquals(7.125576870652436E-5 , elements[0], 1.e-20);
103 Assertions.assertEquals(-1.1135134574790914E-11, elements[1], 1.e-26);
104 Assertions.assertEquals(2.302319084099073E-11 , elements[2], 1.e-26);
105 Assertions.assertEquals(2.4994484564991748E-12 , elements[3], 1.e-27);
106 Assertions.assertEquals(1.381385271417345E-13 , elements[4], 1.e-28);
107 Assertions.assertEquals(5.815883045595586E-12 , elements[5], 1.e-27);
108 }
109
110 @Test
111 void testShortPeriodTerms() throws IllegalArgumentException {
112
113 Utils.setDataRoot("regular-data:potential/icgem-format");
114 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
115 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
116 Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngleType.MEAN,
117 FramesFactory.getTOD(false),
118 new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
119 nshp.getMu());
120
121 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
122 Constants.WGS84_EARTH_FLATTENING,
123 FramesFactory.getGTOD(false));
124
125
126 final DSSTForceModel force = new DSSTTesseral(earth.getBodyFrame(),
127 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
128 nshp, 8, 8, 4, 12, 8, 8, 4);
129
130
131 final SpacecraftState meanState = new SpacecraftState(orbit, 45.0);
132
133
134 final AuxiliaryElements aux = new AuxiliaryElements(orbit, 1);
135
136 final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
137
138 force.registerAttitudeProvider(null);
139
140 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters(meanState.getDate())));
141 force.updateShortPeriodTerms(force.getParametersAllValues(), meanState);
142
143 double[] y = new double[6];
144 for (final ShortPeriodTerms spt : shortPeriodTerms) {
145 final double[] shortPeriodic = spt.value(meanState.getOrbit());
146 for (int i = 0; i < shortPeriodic.length; i++) {
147 y[i] += shortPeriodic[i];
148 }
149 }
150 Assertions.assertEquals(5.192409957353241 , y[0], 1.e-15);
151 Assertions.assertEquals(9.660364749662038E-7 , y[1], 1.e-22);
152 Assertions.assertEquals(1.5420089871620561E-6, y[2], 1.e-21);
153 Assertions.assertEquals(-4.99441460131262E-8 , y[3], 1.e-22);
154 Assertions.assertEquals(-4.500974242661193E-8, y[4], 1.e-22);
155 Assertions.assertEquals(-2.785213556107623E-7, y[5], 1.e-21);
156 }
157
158 @Test
159 void testIssue625() {
160
161
162 final UnnormalizedSphericalHarmonicsProvider provider =
163 GravityFieldFactory.getUnnormalizedProvider(4, 4);
164
165 final Frame frame = FramesFactory.getEME2000();
166 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
167 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
168
169
170
171
172
173
174
175 final Orbit orbit = new EquinoctialOrbit(2.655989E7,
176 2.719455286199036E-4,
177 0.0041543085910249414,
178 -0.3412974060023717,
179 0.3960084733107685,
180 8.566537840341699,
181 PositionAngleType.TRUE,
182 frame,
183 initDate,
184 provider.getMu());
185
186 final SpacecraftState state = new SpacecraftState(orbit, 1000.0);
187
188 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
189
190
191 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
192 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
193 4, 4, 4, 8, 4, 4, 2);
194 tesseral.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, tesseral.getParameters(orbit.getDate()));
195
196
197 final DSSTForceModel tesseralDefault = new DSSTTesseral(earthFrame,
198 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
199 provider);
200 tesseralDefault.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, tesseralDefault.getParameters(orbit.getDate()));
201
202
203 final double[] elements = tesseral.getMeanElementRate(state, auxiliaryElements, tesseral.getParameters(orbit.getDate()));
204
205
206 final double[] elementsDefault = tesseralDefault.getMeanElementRate(state, auxiliaryElements, tesseralDefault.getParameters(orbit.getDate()));
207
208
209 for (int i = 0; i < 6; i++) {
210 Assertions.assertEquals(elements[i], elementsDefault[i], Double.MIN_VALUE);
211 }
212
213 }
214
215 @Test
216 void testIssue736() {
217
218
219 final UnnormalizedSphericalHarmonicsProvider provider =
220 GravityFieldFactory.getUnnormalizedProvider(4, 4);
221
222
223 final Frame frame = FramesFactory.getEME2000();
224 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
225 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
226
227
228 final Orbit orbit = new EquinoctialOrbit(2.655989E7, 2.719455286199036E-4, 0.0041543085910249414,
229 -0.3412974060023717, 0.3960084733107685,
230 8.566537840341699, PositionAngleType.TRUE,
231 frame, initDate, provider.getMu());
232
233
234 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider);
235 final double[] parameters = tesseral.getParameters(orbit.getDate());
236
237
238 tesseral.initializeShortPeriodTerms(new AuxiliaryElements(orbit, 1), PropagationType.MEAN, parameters);
239
240
241 final Orbit shfitedOrbit = new EquinoctialOrbit(2.655989E7, 0.02, 0.0041543085910249414,
242 -0.3412974060023717, 0.3960084733107685,
243 8.566537840341699, PositionAngleType.TRUE,
244 frame, initDate, provider.getMu());
245
246 final double[] elements = tesseral.getMeanElementRate(new SpacecraftState(shfitedOrbit), new AuxiliaryElements(shfitedOrbit, 1), parameters);
247
248
249
250
251 for (int i = 0; i < elements.length; i++) {
252 Assertions.assertTrue(elements[i] != 0);
253 }
254
255 }
256
257
258
259
260
261 @Test
262 void testIssue672() {
263
264
265
266
267
268 final UnnormalizedSphericalHarmonicsProvider provider =
269 GravityFieldFactory.getUnnormalizedProvider(2, 2);
270 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
271
272
273 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400,
274 TimeScalesFactory.getUTC());
275 final Orbit orbit = new KeplerianOrbit(26559890.,
276 0.0041632,
277 FastMath.toRadians(55.2),
278 FastMath.toRadians(315.4985),
279 FastMath.toRadians(130.7562),
280 FastMath.toRadians(44.2377),
281 PositionAngleType.MEAN,
282 FramesFactory.getEME2000(),
283 initDate,
284 provider.getMu());
285
286
287 final SpacecraftState initialState = new SpacecraftState(orbit);
288
289
290 final DSSTTesseral dsstTesseral =
291 new DSSTTesseral(earthFrame,
292 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider);
293
294
295 final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
296 final AuxiliaryElements aux = new AuxiliaryElements(orbit, 1);
297 shortPeriodTerms.addAll(dsstTesseral.initializeShortPeriodTerms(aux,
298 PropagationType.OSCULATING,
299 dsstTesseral.getParameters(orbit.getDate())));
300
301
302
303
304 dsstTesseral.updateShortPeriodTerms(dsstTesseral.getParametersAllValues(), initialState);
305
306
307
308
309
310 Assertions.assertNotNull(shortPeriodTerms);
311
312 }
313
314
315
316
317
318
319 @Test
320 void testIssue672OutOfRangeException() {
321
322
323
324
325
326
327 int degree = 3;
328 int order = 2;
329 UnnormalizedSphericalHarmonicsProvider provider =
330 GravityFieldFactory.getUnnormalizedProvider(degree, order);
331 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
332
333 try {
334
335
336
337 new DSSTTesseral(earthFrame,
338 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
339 degree, order, 3, FastMath.min(12, degree + 4),
340 degree, order, FastMath.min(4, degree - 2));
341 Assertions.fail("An exception should have been thrown");
342 } catch (OrekitException oe) {
343
344
345 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
346 }
347
348
349
350
351
352
353
354 degree = 2;
355 order = 0;
356 provider = GravityFieldFactory.getUnnormalizedProvider(degree, order);
357
358
359
360 final DSSTTesseral dsstTesseral = new DSSTTesseral(earthFrame,
361 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
362 degree, order, 4, FastMath.min(12, degree + 4),
363 degree, order, FastMath.min(4, degree - 2));
364
365 Assertions.assertNotNull(dsstTesseral);
366
367 }
368
369 @Test
370 void testOutOfRangeException() {
371
372 final UnnormalizedSphericalHarmonicsProvider provider =
373 GravityFieldFactory.getUnnormalizedProvider(1, 0);
374
375 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
376 Constants.WGS84_EARTH_FLATTENING,
377 FramesFactory.getGTOD(false));
378 try {
379 @SuppressWarnings("unused")
380 final DSSTForceModel tesseral = new DSSTTesseral(earth.getBodyFrame(),
381 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
382 provider);
383 Assertions.fail("An exception should have been thrown");
384 } catch (OrekitException oe) {
385 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
386 }
387 }
388
389 @Test
390 void testGetMaxEccPow()
391 throws NoSuchMethodException, SecurityException, IllegalAccessException, IllegalArgumentException, InvocationTargetException {
392 final UnnormalizedSphericalHarmonicsProvider provider =
393 GravityFieldFactory.getUnnormalizedProvider(4, 4);;
394 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
395 final DSSTTesseral force = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider);
396 Method getMaxEccPow = DSSTTesseral.class.getDeclaredMethod("getMaxEccPow", Double.TYPE);
397 getMaxEccPow.setAccessible(true);
398 Assertions.assertEquals(3, getMaxEccPow.invoke(force, 0.0));
399 Assertions.assertEquals(4, getMaxEccPow.invoke(force, 0.01));
400 Assertions.assertEquals(7, getMaxEccPow.invoke(force, 0.08));
401 Assertions.assertEquals(10, getMaxEccPow.invoke(force, 0.15));
402 Assertions.assertEquals(12, getMaxEccPow.invoke(force, 0.25));
403 Assertions.assertEquals(15, getMaxEccPow.invoke(force, 0.35));
404 Assertions.assertEquals(20, getMaxEccPow.invoke(force, 1.0));
405 }
406
407 @BeforeEach
408 public void setUp() throws IOException, ParseException {
409 Utils.setDataRoot("regular-data:potential/shm-format");
410 }
411 }