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.text.ParseException;
21 import java.util.ArrayList;
22 import java.util.Arrays;
23 import java.util.List;
24
25 import org.hipparchus.CalculusFieldElement;
26 import org.hipparchus.Field;
27 import org.hipparchus.analysis.differentiation.Gradient;
28 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
29 import org.hipparchus.stat.descriptive.StreamingStatistics;
30 import org.hipparchus.util.Binary64;
31 import org.hipparchus.util.Binary64Field;
32 import org.hipparchus.util.FastMath;
33 import org.hipparchus.util.MathArrays;
34 import org.junit.jupiter.api.Assertions;
35 import org.junit.jupiter.api.BeforeEach;
36 import org.junit.jupiter.api.Test;
37 import org.orekit.Utils;
38 import org.orekit.attitudes.Attitude;
39 import org.orekit.attitudes.AttitudeProvider;
40 import org.orekit.attitudes.FrameAlignedProvider;
41 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
42 import org.orekit.forces.gravity.potential.GRGSFormatReader;
43 import org.orekit.forces.gravity.potential.GravityFieldFactory;
44 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
45 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
46 import org.orekit.frames.Frame;
47 import org.orekit.frames.FramesFactory;
48 import org.orekit.orbits.EquinoctialOrbit;
49 import org.orekit.orbits.FieldCircularOrbit;
50 import org.orekit.orbits.FieldEquinoctialOrbit;
51 import org.orekit.orbits.FieldOrbit;
52 import org.orekit.orbits.Orbit;
53 import org.orekit.orbits.OrbitType;
54 import org.orekit.orbits.PositionAngleType;
55 import org.orekit.propagation.*;
56 import org.orekit.propagation.numerical.FieldNumericalPropagator;
57 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
58 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
59 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
60 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
61 import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
62 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
63 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
64 import org.orekit.time.AbsoluteDate;
65 import org.orekit.time.DateComponents;
66 import org.orekit.time.FieldAbsoluteDate;
67 import org.orekit.time.TimeComponents;
68 import org.orekit.time.TimeScalesFactory;
69 import org.orekit.utils.IERSConventions;
70 import org.orekit.utils.ParameterDriver;
71 import org.orekit.utils.ParameterDriversList;
72
73 public class FieldDSSTZonalTest {
74
75 @Test
76 public void testGetMeanElementRate() {
77 doTestGetMeanElementRate(Binary64Field.getInstance());
78 }
79
80 private <T extends CalculusFieldElement<T>> void doTestGetMeanElementRate(final Field<T> field) {
81
82 final T zero = field.getZero();
83
84
85 final UnnormalizedSphericalHarmonicsProvider provider =
86 GravityFieldFactory.getUnnormalizedProvider(4, 4);
87
88 final Frame earthFrame = FramesFactory.getEME2000();
89 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2007, 04, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
90
91
92
93
94
95
96
97 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(2.655989E7),
98 zero.add(2.719455286199036E-4),
99 zero.add(0.0041543085910249414),
100 zero.add(-0.3412974060023717),
101 zero.add(0.3960084733107685),
102 zero.add(8.566537840341699),
103 PositionAngleType.TRUE,
104 earthFrame,
105 initDate,
106 zero.add(3.986004415E14));
107
108 final T mass = zero.add(1000.0);
109 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit, mass);
110
111 final DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
112
113 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), 1);
114
115
116 final T[] parameters = zonal.getParameters(field, state.getDate());
117
118 zonal.initializeShortPeriodTerms(auxiliaryElements,
119 PropagationType.MEAN, parameters);
120
121 final T[] elements = MathArrays.buildArray(field, 7);
122 Arrays.fill(elements, zero);
123
124 final T[] daidt = zonal.getMeanElementRate(state, auxiliaryElements, parameters);
125 for (int i = 0; i < daidt.length; i++) {
126 elements[i] = daidt[i];
127 }
128
129 Assertions.assertEquals(0.0, elements[0].getReal(), 1.e-25);
130 Assertions.assertEquals(1.3909396722346468E-11, elements[1].getReal(), 3.e-26);
131 Assertions.assertEquals(-2.0275977261372793E-13, elements[2].getReal(), 3.e-27);
132 Assertions.assertEquals(3.087141512018238E-9, elements[3].getReal(), 1.e-24);
133 Assertions.assertEquals(2.6606317310148797E-9, elements[4].getReal(), 4.e-24);
134 Assertions.assertEquals(-3.659904725206694E-9, elements[5].getReal(), 1.e-24);
135
136 }
137
138 @Test
139 public void testShortPeriodTerms() {
140 doTestShortPeriodTerms(Binary64Field.getInstance());
141 }
142
143 @SuppressWarnings("unchecked")
144 private <T extends CalculusFieldElement<T>> void doTestShortPeriodTerms(final Field<T> field) {
145 final T zero = field.getZero();
146
147 final FieldSpacecraftState<T> meanState = getGEOState(field);
148
149 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
150 final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
151
152
153 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(meanState.getOrbit(), 1);
154
155
156 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
157
158 zonal.registerAttitudeProvider(null);
159 shortPeriodTerms.addAll(zonal.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, zonal.getParameters(field)));
160 zonal.updateShortPeriodTerms(zonal.getParametersAllValues(field), meanState);
161
162 T[] y = MathArrays.buildArray(field, 6);
163 Arrays.fill(y, zero);
164 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
165 final T[] shortPeriodic = spt.value(meanState.getOrbit());
166 for (int i = 0; i < shortPeriodic.length; i++) {
167 y[i] = y[i].add(shortPeriodic[i]);
168 }
169 }
170
171 Assertions.assertEquals(35.005618980090276, y[0].getReal(), 1.e-15);
172 Assertions.assertEquals(3.75891551882889E-5, y[1].getReal(), 1.e-20);
173 Assertions.assertEquals(3.929119925563796E-6, y[2].getReal(), 1.e-21);
174 Assertions.assertEquals(-1.1781951949124315E-8, y[3].getReal(), 1.e-24);
175 Assertions.assertEquals(-3.2134924513679615E-8, y[4].getReal(), 1.e-24);
176 Assertions.assertEquals(-1.1607392915997098E-6, y[5].getReal(), 1.e-21);
177 }
178
179 @Test
180 public void testIssue625() {
181 doTestIssue625(Binary64Field.getInstance());
182 }
183
184 private <T extends CalculusFieldElement<T>> void doTestIssue625(final Field<T> field) {
185
186 Utils.setDataRoot("regular-data:potential/grgs-format");
187 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
188
189 final T zero = field.getZero();
190
191 final Frame earthFrame = FramesFactory.getEME2000();
192 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2007, 04, 16, 0, 46, 42.400, TimeScalesFactory.getUTC());
193
194
195
196
197
198
199
200 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(2.655989E6),
201 zero.add(2.719455286199036E-4),
202 zero.add(0.0041543085910249414),
203 zero.add(-0.3412974060023717),
204 zero.add(0.3960084733107685),
205 zero.add(8.566537840341699),
206 PositionAngleType.TRUE,
207 earthFrame,
208 initDate,
209 zero.add(3.986004415E14));
210
211 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit, zero.add(1000.0));
212
213 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), 1);
214
215
216 final UnnormalizedSphericalHarmonicsProvider provider =
217 GravityFieldFactory.getUnnormalizedProvider(32, 32);
218
219
220 final DSSTZonal zonal = new DSSTZonal(provider, 32, 4, 65);
221 zonal.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, zonal.getParameters(field, state.getDate()));
222
223
224 final DSSTZonal zonalDefault = new DSSTZonal(provider);
225 zonalDefault.initializeShortPeriodTerms(auxiliaryElements, PropagationType.MEAN, zonalDefault.getParameters(field, state.getDate()));
226
227
228 final T[] elements = zonal.getMeanElementRate(state, auxiliaryElements, zonal.getParameters(field, state.getDate()));
229
230
231 final T[] elementsDefault = zonalDefault.getMeanElementRate(state, auxiliaryElements, zonalDefault.getParameters(field, state.getDate()));
232
233
234 for (int i = 0; i < 6; i++) {
235 Assertions.assertEquals(elements[i].getReal(), elementsDefault[i].getReal(), Double.MIN_VALUE);
236 }
237
238 }
239
240 @Test
241 @SuppressWarnings("unchecked")
242 public void testShortPeriodTermsStateDerivatives() {
243
244
245 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
246 TimeScalesFactory.getUTC());
247
248 final Orbit orbit = new EquinoctialOrbit(42164000,
249 10e-3,
250 10e-3,
251 FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
252 FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
253 PositionAngleType.TRUE,
254 FramesFactory.getEME2000(),
255 initDate,
256 3.986004415E14);
257
258 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
259
260 final SpacecraftState meanState = new SpacecraftState(orbit);
261
262
263 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
264 final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
265
266
267 final DSSTGradientConverter converter = new DSSTGradientConverter(meanState, Utils.defaultLaw());
268
269
270 final FieldSpacecraftState<Gradient> dsState = converter.getState(zonal);
271
272 final FieldAuxiliaryElements<Gradient> fieldAuxiliaryElements = new FieldAuxiliaryElements<>(dsState.getOrbit(), 1);
273
274
275 final Gradient zero = dsState.getDate().getField().getZero();
276
277
278 final List<FieldShortPeriodTerms<Gradient>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<Gradient>>();
279 shortPeriodTerms.addAll(zonal.initializeShortPeriodTerms(fieldAuxiliaryElements, PropagationType.OSCULATING,
280 converter.getParametersAtStateDate(dsState, zonal)));
281 zonal.updateShortPeriodTerms(converter.getParameters(dsState, zonal), dsState);
282 final Gradient[] shortPeriod = new Gradient[6];
283 Arrays.fill(shortPeriod, zero);
284 for (final FieldShortPeriodTerms<Gradient> spt : shortPeriodTerms) {
285 final Gradient[] spVariation = spt.value(dsState.getOrbit());
286 for (int i = 0; i < spVariation .length; i++) {
287 shortPeriod[i] = shortPeriod[i].add(spVariation[i]);
288 }
289 }
290
291 final double[][] shortPeriodJacobian = new double[6][6];
292
293 final double[] derivativesASP = shortPeriod[0].getGradient();
294 final double[] derivativesExSP = shortPeriod[1].getGradient();
295 final double[] derivativesEySP = shortPeriod[2].getGradient();
296 final double[] derivativesHxSP = shortPeriod[3].getGradient();
297 final double[] derivativesHySP = shortPeriod[4].getGradient();
298 final double[] derivativesLSP = shortPeriod[5].getGradient();
299
300
301 addToRow(derivativesASP, 0, shortPeriodJacobian);
302 addToRow(derivativesExSP, 1, shortPeriodJacobian);
303 addToRow(derivativesEySP, 2, shortPeriodJacobian);
304 addToRow(derivativesHxSP, 3, shortPeriodJacobian);
305 addToRow(derivativesHySP, 4, shortPeriodJacobian);
306 addToRow(derivativesLSP, 5, shortPeriodJacobian);
307
308
309 double[][] shortPeriodJacobianRef = new double[6][6];
310 double dP = 0.001;
311 double[] steps = ToleranceProvider.getDefaultToleranceProvider(1000000 * dP).getTolerances(orbit, orbitType)[0];
312 for (int i = 0; i < 6; i++) {
313
314 SpacecraftState stateM4 = shiftState(meanState, orbitType, -4 * steps[i], i);
315 double[] shortPeriodM4 = computeShortPeriodTerms(stateM4, zonal);
316
317 SpacecraftState stateM3 = shiftState(meanState, orbitType, -3 * steps[i], i);
318 double[] shortPeriodM3 = computeShortPeriodTerms(stateM3, zonal);
319
320 SpacecraftState stateM2 = shiftState(meanState, orbitType, -2 * steps[i], i);
321 double[] shortPeriodM2 = computeShortPeriodTerms(stateM2, zonal);
322
323 SpacecraftState stateM1 = shiftState(meanState, orbitType, -1 * steps[i], i);
324 double[] shortPeriodM1 = computeShortPeriodTerms(stateM1, zonal);
325
326 SpacecraftState stateP1 = shiftState(meanState, orbitType, 1 * steps[i], i);
327 double[] shortPeriodP1 = computeShortPeriodTerms(stateP1, zonal);
328
329 SpacecraftState stateP2 = shiftState(meanState, orbitType, 2 * steps[i], i);
330 double[] shortPeriodP2 = computeShortPeriodTerms(stateP2, zonal);
331
332 SpacecraftState stateP3 = shiftState(meanState, orbitType, 3 * steps[i], i);
333 double[] shortPeriodP3 = computeShortPeriodTerms(stateP3, zonal);
334
335 SpacecraftState stateP4 = shiftState(meanState, orbitType, 4 * steps[i], i);
336 double[] shortPeriodP4 = computeShortPeriodTerms(stateP4, zonal);
337
338 fillJacobianColumn(shortPeriodJacobianRef, i, orbitType, steps[i],
339 shortPeriodM4, shortPeriodM3, shortPeriodM2, shortPeriodM1,
340 shortPeriodP1, shortPeriodP2, shortPeriodP3, shortPeriodP4);
341
342 }
343
344 for (int m = 0; m < 6; ++m) {
345 for (int n = 0; n < 6; ++n) {
346 double error = FastMath.abs((shortPeriodJacobian[m][n] - shortPeriodJacobianRef[m][n]) / shortPeriodJacobianRef[m][n]);
347 Assertions.assertEquals(0, error, 9.6e-10);
348 }
349 }
350
351 }
352
353 @Test
354 @SuppressWarnings("unchecked")
355 public void testShortPeriodTermsMuParametersDerivatives() {
356
357
358 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
359 TimeScalesFactory.getUTC());
360
361 final Orbit orbit = new EquinoctialOrbit(42164000,
362 10e-3,
363 10e-3,
364 FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
365 FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
366 PositionAngleType.TRUE,
367 FramesFactory.getEME2000(),
368 initDate,
369 3.986004415E14);
370
371 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
372
373 final SpacecraftState meanState = new SpacecraftState(orbit);
374
375 final double[] stateVector = new double[6];
376 OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngleType.MEAN, stateVector, null);
377
378
379 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
380 final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
381
382 for (final ParameterDriver driver : zonal.getParametersDrivers()) {
383 driver.setValue(driver.getReferenceValue());
384 driver.setSelected(driver.getName().equals(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT));
385 }
386
387
388 final DSSTGradientConverter converter = new DSSTGradientConverter(meanState, Utils.defaultLaw());
389
390
391 final FieldSpacecraftState<Gradient> dsState = converter.getState(zonal);
392
393 final FieldAuxiliaryElements<Gradient> fieldAuxiliaryElements = new FieldAuxiliaryElements<>(dsState.getOrbit(), 1);
394
395
396 final Gradient zero = dsState.getDate().getField().getZero();
397
398
399 final List<FieldShortPeriodTerms<Gradient>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<Gradient>>();
400 shortPeriodTerms.addAll(zonal.initializeShortPeriodTerms(fieldAuxiliaryElements, PropagationType.OSCULATING,
401 converter.getParametersAtStateDate(dsState, zonal)));
402 zonal.updateShortPeriodTerms(converter.getParameters(dsState, zonal), dsState);
403 final Gradient[] shortPeriod = new Gradient[6];
404 Arrays.fill(shortPeriod, zero);
405 for (final FieldShortPeriodTerms<Gradient> spt : shortPeriodTerms) {
406 final Gradient[] spVariation = spt.value(dsState.getOrbit());
407 for (int i = 0; i < spVariation .length; i++) {
408 shortPeriod[i] = shortPeriod[i].add(spVariation[i]);
409 }
410 }
411
412 final double[][] shortPeriodJacobian = new double[6][1];
413
414 final double[] derivativesASP = shortPeriod[0].getGradient();
415 final double[] derivativesExSP = shortPeriod[1].getGradient();
416 final double[] derivativesEySP = shortPeriod[2].getGradient();
417 final double[] derivativesHxSP = shortPeriod[3].getGradient();
418 final double[] derivativesHySP = shortPeriod[4].getGradient();
419 final double[] derivativesLSP = shortPeriod[5].getGradient();
420
421 int index = converter.getFreeStateParameters();
422 for (ParameterDriver driver : zonal.getParametersDrivers()) {
423 if (driver.isSelected()) {
424 shortPeriodJacobian[0][0] += derivativesASP[index];
425 shortPeriodJacobian[1][0] += derivativesExSP[index];
426 shortPeriodJacobian[2][0] += derivativesEySP[index];
427 shortPeriodJacobian[3][0] += derivativesHxSP[index];
428 shortPeriodJacobian[4][0] += derivativesHySP[index];
429 shortPeriodJacobian[5][0] += derivativesLSP[index];
430 ++index;
431 }
432 }
433
434
435 double[][] shortPeriodJacobianRef = new double[6][1];
436 ParameterDriversList bound = new ParameterDriversList();
437 for (final ParameterDriver driver : zonal.getParametersDrivers()) {
438 if (driver.getName().equals(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)) {
439 driver.setSelected(true);
440 bound.add(driver);
441 } else {
442 driver.setSelected(false);
443 }
444 }
445
446 ParameterDriver selected = bound.getDrivers().get(0);
447 double p0 = selected.getReferenceValue();
448 double h = selected.getScale();
449
450 selected.setValue(p0 - 4 * h);
451 final double[] shortPeriodM4 = computeShortPeriodTerms(meanState, zonal);
452
453 selected.setValue(p0 - 3 * h);
454 final double[] shortPeriodM3 = computeShortPeriodTerms(meanState, zonal);
455
456 selected.setValue(p0 - 2 * h);
457 final double[] shortPeriodM2 = computeShortPeriodTerms(meanState, zonal);
458
459 selected.setValue(p0 - 1 * h);
460 final double[] shortPeriodM1 = computeShortPeriodTerms(meanState, zonal);
461
462 selected.setValue(p0 + 1 * h);
463 final double[] shortPeriodP1 = computeShortPeriodTerms(meanState, zonal);
464
465 selected.setValue(p0 + 2 * h);
466 final double[] shortPeriodP2 = computeShortPeriodTerms(meanState, zonal);
467
468 selected.setValue(p0 + 3 * h);
469 final double[] shortPeriodP3 = computeShortPeriodTerms(meanState, zonal);
470
471 selected.setValue(p0 + 4 * h);
472 final double[] shortPeriodP4 = computeShortPeriodTerms(meanState, zonal);
473
474 fillJacobianColumn(shortPeriodJacobianRef, 0, orbitType, h,
475 shortPeriodM4, shortPeriodM3, shortPeriodM2, shortPeriodM1,
476 shortPeriodP1, shortPeriodP2, shortPeriodP3, shortPeriodP4);
477
478 for (int i = 0; i < 6; ++i) {
479 double error = FastMath.abs((shortPeriodJacobian[i][0] - shortPeriodJacobianRef[i][0]) / stateVector[i]) * h;
480 Assertions.assertEquals(0, error, 1.3e-18);
481 }
482
483 }
484
485 private <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> getGEOState(final Field<T> field) {
486
487 final T zero = field.getZero();
488
489 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
490 TimeScalesFactory.getUTC());
491 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(42164000),
492 zero.add(10e-3),
493 zero.add(10e-3),
494 zero.add(FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3)),
495 zero.add(FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3)),
496 zero.add(0.1),
497 PositionAngleType.TRUE,
498 FramesFactory.getEME2000(),
499 initDate,
500 zero.add(3.986004415E14));
501 return new FieldSpacecraftState<>(orbit);
502 }
503
504
505
506
507
508 @Test
509 void testIssue1104() {
510
511 final boolean printResults = false;
512
513 final Field<Binary64> field = Binary64Field.getInstance();
514
515
516 final Frame gcrf = FramesFactory.getGCRF();
517 final Frame tod = FramesFactory.getTOD(IERSConventions.IERS_2010, true);
518 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
519
520
521
522
523
524 double diMax = 9.615e-5;
525 double dOmMax = 3.374e-3;
526 double dLmMax = 1.128e-2;
527 doTestIssue1104(gcrf, gcrf, field, printResults, diMax, dOmMax, dLmMax);
528
529
530
531
532
533
534
535 diMax = 1.059e-5;
536 dOmMax = 2.789e-3;
537 dLmMax = 1.040e-2;
538 doTestIssue1104(tod, tod, field, printResults, diMax, dOmMax, dLmMax);
539
540
541
542
543
544
545 diMax = 1.067e-5;
546 dOmMax = 2.789e-3;
547 dLmMax = 1.040e-2;
548 doTestIssue1104(gcrf, itrf, field, printResults, diMax, dOmMax, dLmMax);
549
550
551
552
553
554
555 diMax = 1.059e-5;
556 dOmMax = 2.789e-3;
557 dLmMax = 1.040e-2;
558 doTestIssue1104(tod, itrf, field, printResults, diMax, dOmMax, dLmMax);
559
560
561
562 }
563
564
565 private <T extends CalculusFieldElement<T>> void doTestIssue1104(final Frame inertialFrame,
566 final Frame bodyFixedFrame,
567 final Field<T> field,
568 final boolean printResults,
569 final double diMax,
570 final double dOmMax,
571 final double dLmMax) {
572
573
574
575
576
577 final T zero = field.getZero();
578 final double step = 60.;
579 final double nOrb = 50.;
580
581 final FieldAbsoluteDate<T> t0 = new FieldAbsoluteDate<>(field);
582
583
584 final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
585
586
587 final int degree = 2;
588 final int order = 0;
589 final UnnormalizedSphericalHarmonicsProvider unnormalized =
590 GravityFieldFactory.getConstantUnnormalizedProvider(degree, order, t0.toAbsoluteDate());
591 final NormalizedSphericalHarmonicsProvider normalized =
592 GravityFieldFactory.getConstantNormalizedProvider(degree, order, t0.toAbsoluteDate());
593
594
595 final double mass = 150.;
596 final double a = 6906780.35;
597 final double ex = 5.09E-4;
598 final double ey = 1.24e-3;
599 final double i = FastMath.toRadians(97.49);
600 final double raan = FastMath.toRadians(-94.607);
601 final double alphaM = FastMath.toRadians(0.);
602 final FieldCircularOrbit<T> oscCircOrbit0 = new FieldCircularOrbit<>(
603 zero.newInstance(a),
604 zero.newInstance(ex),
605 zero.newInstance(ey),
606 zero.newInstance(i),
607 zero.newInstance(raan),
608 zero.newInstance(alphaM),
609 PositionAngleType.MEAN,
610 inertialFrame,
611 t0,
612 zero.newInstance(unnormalized.getMu()));
613
614 final FieldOrbit<T> oscOrbit0 = new FieldEquinoctialOrbit<>(oscCircOrbit0);
615 final FieldSpacecraftState<T> oscState0 = new FieldSpacecraftState<>(oscOrbit0, zero.newInstance(mass));
616 final AttitudeProvider attProvider = new FrameAlignedProvider(inertialFrame);
617
618
619 final double duration = nOrb * oscOrbit0.getKeplerianPeriod().getReal();
620 final FieldAbsoluteDate<T> tf = t0.shiftedBy(duration);
621
622
623 final ClassicalRungeKuttaFieldIntegrator<T> integrator =
624 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.newInstance(step));
625
626 final FieldNumericalPropagator<T> numProp = new FieldNumericalPropagator<>(field, integrator);
627 numProp.setOrbitType(oscOrbit0.getType());
628 numProp.setInitialState(oscState0);
629 numProp.setAttitudeProvider(attProvider);
630 numProp.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, normalized));
631 final FieldEphemerisGenerator<T> numEphemGen = numProp.getEphemerisGenerator();
632
633
634 final ClassicalRungeKuttaFieldIntegrator<T> integratorDsst =
635 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.newInstance(step));
636 final FieldDSSTPropagator<T> dsstProp = new FieldDSSTPropagator<T>(field, integratorDsst, PropagationType.OSCULATING);
637 dsstProp.setInitialState(oscState0, PropagationType.OSCULATING);
638 dsstProp.setAttitudeProvider(attProvider);
639 final DSSTForceModel zonal = new DSSTZonal(bodyFixedFrame, unnormalized);
640 dsstProp.addForceModel(zonal);
641 final FieldEphemerisGenerator<T> dsstEphemGen = dsstProp.getEphemerisGenerator();
642
643
644
645
646
647
648 final StreamingStatistics dI = new StreamingStatistics();
649 final StreamingStatistics dOm = new StreamingStatistics();
650 final StreamingStatistics dLM = new StreamingStatistics();
651
652
653 numProp.propagate(t0, tf);
654 dsstProp.propagate(t0, tf);
655
656 final FieldBoundedPropagator<T> numEphem = numEphemGen.getGeneratedEphemeris();
657 final FieldBoundedPropagator<T> dsstEphem = dsstEphemGen.getGeneratedEphemeris();
658
659
660 for (double dt = 0; dt < duration; dt += step) {
661
662
663 final FieldAbsoluteDate<T> t = t0.shiftedBy(dt);
664
665
666 final FieldCircularOrbit<T> num = new FieldCircularOrbit<>(numEphem.propagate(t).getOrbit());
667 final FieldCircularOrbit<T> dsst = new FieldCircularOrbit<>(dsstEphem.propagate(t).getOrbit());
668 dI.addValue(FastMath.toDegrees(dsst.getI().getReal() - num.getI().getReal()));
669 dOm.addValue(FastMath.toDegrees(dsst.getRightAscensionOfAscendingNode().getReal() -
670 num.getRightAscensionOfAscendingNode().getReal()));
671 dLM.addValue(FastMath.toDegrees(dsst.getLM().getReal() - num.getLM().getReal()));
672 }
673
674
675
676
677
678 if (printResults) {
679 System.out.println("Inertial frame : " + inertialFrame.toString());
680 System.out.println("Body-Fixed frame: " + bodyFixedFrame.toString());
681 System.out.println("\ndi\n" + dI.toString());
682 System.out.println("\ndĪ©\n" + dOm.toString());
683 System.out.println("\ndLM\n" + dLM.toString());
684 }
685
686
687 Assertions.assertEquals(diMax, FastMath.max(FastMath.abs(dI.getMax()), FastMath.abs(dI.getMin())), 1.e-8);
688 Assertions.assertEquals(dOmMax, FastMath.max(FastMath.abs(dOm.getMax()), FastMath.abs(dOm.getMin())), 1.e-6);
689 Assertions.assertEquals(dLmMax, FastMath.max(FastMath.abs(dLM.getMax()), FastMath.abs(dLM.getMin())), 1.e-5);
690 }
691
692 private double[] computeShortPeriodTerms(SpacecraftState state,
693 DSSTForceModel force) {
694
695 AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
696
697 List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
698 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(auxiliaryElements, PropagationType.OSCULATING, force.getParameters(state.getDate())));
699 force.updateShortPeriodTerms(force.getParametersAllValues(), state);
700
701 double[] shortPeriod = new double[6];
702 for (ShortPeriodTerms spt : shortPeriodTerms) {
703 double[] spVariation = spt.value(state.getOrbit());
704 for (int i = 0; i < spVariation.length; i++) {
705 shortPeriod[i] += spVariation[i];
706 }
707 }
708
709 return shortPeriod;
710
711 }
712
713 private void fillJacobianColumn(double[][] jacobian, int column,
714 OrbitType orbitType, double h,
715 double[] M4h, double[] M3h,
716 double[] M2h, double[] M1h,
717 double[] P1h, double[] P2h,
718 double[] P3h, double[] P4h) {
719 for (int i = 0; i < jacobian.length; ++i) {
720 jacobian[i][column] = ( -3 * (P4h[i] - M4h[i]) +
721 32 * (P3h[i] - M3h[i]) -
722 168 * (P2h[i] - M2h[i]) +
723 672 * (P1h[i] - M1h[i])) / (840 * h);
724 }
725 }
726
727 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
728 double delta, int column) {
729
730 double[][] array = stateToArray(state, orbitType);
731 array[0][column] += delta;
732
733 return arrayToState(array, orbitType, state.getFrame(), state.getDate(),
734 state.getOrbit().getMu(), state.getAttitude());
735
736 }
737
738 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
739 double[][] array = new double[2][6];
740
741 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
742 return array;
743 }
744
745 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
746 Frame frame, AbsoluteDate date, double mu,
747 Attitude attitude) {
748 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngleType.MEAN, date, mu, frame);
749 return new SpacecraftState(orbit, attitude);
750 }
751
752
753
754
755
756
757 private void addToRow(final double[] derivatives, final int index,
758 final double[][] jacobian) {
759
760 for (int i = 0; i < 6; i++) {
761 jacobian[index][i] += derivatives[i];
762 }
763
764 }
765
766 @BeforeEach
767 public void setUp() throws IOException, ParseException {
768 Utils.setDataRoot("regular-data:potential/shm-format");
769 }
770
771 }