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