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 org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.differentiation.Gradient;
22 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
23 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Rotation;
25 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
26 import org.hipparchus.geometry.euclidean.threed.Vector3D;
27 import org.hipparchus.util.Binary64Field;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.MathArrays;
30 import org.junit.jupiter.api.Assertions;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.Attitude;
35 import org.orekit.attitudes.AttitudeProvider;
36 import org.orekit.attitudes.FieldAttitude;
37 import org.orekit.attitudes.FrameAlignedProvider;
38 import org.orekit.attitudes.LofOffset;
39 import org.orekit.bodies.CelestialBody;
40 import org.orekit.bodies.CelestialBodyFactory;
41 import org.orekit.bodies.OneAxisEllipsoid;
42 import org.orekit.forces.BoxAndSolarArraySpacecraft;
43 import org.orekit.forces.gravity.potential.GravityFieldFactory;
44 import org.orekit.forces.gravity.potential.SHMFormatReader;
45 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
46 import org.orekit.forces.radiation.RadiationSensitive;
47 import org.orekit.frames.Frame;
48 import org.orekit.frames.FramesFactory;
49 import org.orekit.frames.LOFType;
50 import org.orekit.orbits.EquinoctialOrbit;
51 import org.orekit.orbits.FieldEquinoctialOrbit;
52 import org.orekit.orbits.FieldOrbit;
53 import org.orekit.orbits.Orbit;
54 import org.orekit.orbits.OrbitType;
55 import org.orekit.orbits.PositionAngleType;
56 import org.orekit.propagation.FieldSpacecraftState;
57 import org.orekit.propagation.PropagationType;
58 import org.orekit.propagation.SpacecraftState;
59 import org.orekit.propagation.ToleranceProvider;
60 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
61 import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
62 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
63 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
64 import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
65 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
66 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
67 import org.orekit.time.AbsoluteDate;
68 import org.orekit.time.DateComponents;
69 import org.orekit.time.FieldAbsoluteDate;
70 import org.orekit.time.TimeComponents;
71 import org.orekit.time.TimeScalesFactory;
72 import org.orekit.utils.Constants;
73 import org.orekit.utils.IERSConventions;
74 import org.orekit.utils.ParameterDriver;
75 import org.orekit.utils.ParameterDriversList;
76 import org.orekit.utils.TimeStampedFieldAngularCoordinates;
77
78 import java.io.IOException;
79 import java.text.ParseException;
80 import java.util.ArrayList;
81 import java.util.Arrays;
82 import java.util.List;
83
84 class FieldDSSTSolarRadiationPressureTest {
85
86 @Test
87 void testGetMeanElementRate() {
88 doTestGetMeanElementRate(Binary64Field.getInstance());
89 }
90
91 private <T extends CalculusFieldElement<T>> void doTestGetMeanElementRate(final Field<T> field) {
92
93 final T zero = field.getZero();
94
95 final Frame earthFrame = FramesFactory.getGCRF();
96 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2003, 9, 16, 0, 0, 0, TimeScalesFactory.getUTC());
97 final double mu = 3.986004415E14;
98
99
100
101
102
103
104 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(4.2166258E7),
105 zero.add(6.532127416888538E-6),
106 zero.add(9.978642849310487E-5),
107 zero.add(-5.69711879850274E-6),
108 zero.add(6.61038518895005E-6),
109 zero.add(8.56084687583949),
110 PositionAngleType.TRUE,
111 earthFrame,
112 initDate,
113 zero.add(mu));
114
115
116 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
117 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
118 Constants.WGS84_EARTH_FLATTENING,
119 FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
120 mu);
121
122
123 Rotation rotation = new Rotation(0.9999999999999984,
124 1.653020584550675E-8,
125 -4.028108631990782E-8,
126 -3.539139805514139E-8,
127 false);
128 AttitudeProvider attitudeProvider = new FrameAlignedProvider(rotation);
129 srp.registerAttitudeProvider(attitudeProvider);
130
131
132 FieldRotation<T> fieldRotation = new FieldRotation<>(field, rotation);
133 FieldVector3D<T> rotationRate = new FieldVector3D<>(zero, zero, zero);
134 FieldVector3D<T> rotationAcceleration = new FieldVector3D<>(zero, zero, zero);
135 TimeStampedFieldAngularCoordinates<T> orientation = new TimeStampedFieldAngularCoordinates<>(initDate,
136 fieldRotation,
137 rotationRate,
138 rotationAcceleration);
139 final FieldAttitude<T> att = new FieldAttitude<>(earthFrame, orientation);
140
141
142 final T mass = zero.add(1000.0);
143 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit, att, mass);
144 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), 1);
145
146
147 final T[] parameters = srp.getParameters(field);
148
149 srp.initializeShortPeriodTerms(auxiliaryElements,
150 PropagationType.MEAN, parameters);
151
152
153 final T[] elements = MathArrays.buildArray(field, 7);
154 Arrays.fill(elements, zero);
155 final T[] daidt = srp.getMeanElementRate(state, auxiliaryElements, parameters);
156 for (int i = 0; i < daidt.length; i++) {
157 elements[i] = daidt[i];
158 }
159
160 Assertions.assertEquals( 6.839644083279315E-8, elements[0].getReal(), 1.0e-23);
161 Assertions.assertEquals(-2.990943651374143E-11, elements[1].getReal(), 1.0e-27);
162 Assertions.assertEquals(-2.5384000825777143E-10, elements[2].getReal(), 1.0e-26);
163 Assertions.assertEquals( 2.0378397296268607E-13, elements[3].getReal(), 1.0e-29);
164 Assertions.assertEquals(-2.33389107008235E-14, elements[4].getReal(), 1.0e-30);
165 Assertions.assertEquals( 1.6082477070964215E-11, elements[5].getReal(), 1.0e-27);
166
167 }
168
169 @Test
170 void testShortPeriodTerms() {
171 doTestShortPeriodTerms(Binary64Field.getInstance());
172 }
173
174 @SuppressWarnings("unchecked")
175 private <T extends CalculusFieldElement<T>> void doTestShortPeriodTerms(final Field<T> field) {
176
177 final T zero = field.getZero();
178 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
179
180 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(7069219.9806427825),
181 zero.add(-4.5941811292223825E-4),
182 zero.add(1.309932339472599E-4),
183 zero.add(-1.002996107003202),
184 zero.add(0.570979900577994),
185 zero.add(2.62038786211518),
186 PositionAngleType.TRUE,
187 FramesFactory.getEME2000(),
188 initDate,
189 zero.add(3.986004415E14));
190
191 final FieldSpacecraftState<T> meanState = new FieldSpacecraftState<>(orbit);
192
193 final CelestialBody sun = CelestialBodyFactory.getSun();
194
195 final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
196 sun,
197 50.0, Vector3D.PLUS_J,
198 2.0, 0.1,
199 0.2, 0.6);
200
201 final AttitudeProvider attitudeProvider = new LofOffset(meanState.getFrame(),
202 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
203 0.0, 0.0, 0.0);
204
205 final DSSTForceModel srp = new DSSTSolarRadiationPressure(sun,
206 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
207 Constants.WGS84_EARTH_FLATTENING,
208 FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
209 boxAndWing,
210 meanState.getOrbit().getMu().getReal());
211
212
213 final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(meanState.getOrbit(), 1);
214
215
216 final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
217
218 srp.registerAttitudeProvider(attitudeProvider);
219 shortPeriodTerms.addAll(srp.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, srp.getParameters(field)));
220 srp.updateShortPeriodTerms(srp.getParametersAllValues(field), meanState);
221
222 T[] y = MathArrays.buildArray(field, 6);
223 Arrays.fill(y, zero);
224
225 for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
226 final T[] shortPeriodic = spt.value(meanState.getOrbit());
227 for (int i = 0; i < shortPeriodic.length; i++) {
228 y[i] = y[i].add(shortPeriodic[i]);
229 }
230 }
231
232 Assertions.assertEquals(0.3668654523023674, y[0].getReal(), 1.0e-15);
233 Assertions.assertEquals(-2.5673332283029E-10, y[1].getReal(), 1.0e-23);
234 Assertions.assertEquals(-3.84959877691874E-9, y[2].getReal(), 1.0e-23);
235 Assertions.assertEquals(-3.069285299519465E-9, y[3].getReal(), 1.0e-23);
236 Assertions.assertEquals(-4.90887054227722E-9, y[4].getReal(), 1.0e-23);
237 Assertions.assertEquals(-2.38549338428378E-9, y[5].getReal(), 1.0e-20);
238 }
239
240 @Test
241 @SuppressWarnings("unchecked")
242 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 AttitudeProvider attitudeProvider = new LofOffset(meanState.getFrame(),
264 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
265 0.0, 0.0, 0.0);
266
267
268 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
269 final DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
270 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
271 Constants.WGS84_EARTH_FLATTENING,
272 FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
273 provider.getMu());
274 srp.registerAttitudeProvider(attitudeProvider);
275
276
277 final DSSTGradientConverter converter = new DSSTGradientConverter(meanState, Utils.defaultLaw());
278
279
280 final FieldSpacecraftState<Gradient> dsState = converter.getState(srp);
281
282 final FieldAuxiliaryElements<Gradient> fieldAuxiliaryElements = new FieldAuxiliaryElements<>(dsState.getOrbit(), 1);
283
284
285 final Gradient zero = dsState.getDate().getField().getZero();
286
287
288 final List<FieldShortPeriodTerms<Gradient>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<Gradient>>();
289 shortPeriodTerms.addAll(srp.initializeShortPeriodTerms(fieldAuxiliaryElements, PropagationType.OSCULATING,
290 converter.getParametersAtStateDate(dsState, srp)));
291 srp.updateShortPeriodTerms(converter.getParameters(dsState, srp), dsState);
292 final Gradient[] shortPeriod = new Gradient[6];
293 Arrays.fill(shortPeriod, zero);
294 for (final FieldShortPeriodTerms<Gradient> spt : shortPeriodTerms) {
295 final Gradient[] spVariation = spt.value(dsState.getOrbit());
296 for (int i = 0; i < spVariation .length; i++) {
297 shortPeriod[i] = shortPeriod[i].add(spVariation[i]);
298 }
299 }
300
301 final double[][] shortPeriodJacobian = new double[6][6];
302
303 final double[] derivativesASP = shortPeriod[0].getGradient();
304 final double[] derivativesExSP = shortPeriod[1].getGradient();
305 final double[] derivativesEySP = shortPeriod[2].getGradient();
306 final double[] derivativesHxSP = shortPeriod[3].getGradient();
307 final double[] derivativesHySP = shortPeriod[4].getGradient();
308 final double[] derivativesLSP = shortPeriod[5].getGradient();
309
310
311 addToRow(derivativesASP, 0, shortPeriodJacobian);
312 addToRow(derivativesExSP, 1, shortPeriodJacobian);
313 addToRow(derivativesEySP, 2, shortPeriodJacobian);
314 addToRow(derivativesHxSP, 3, shortPeriodJacobian);
315 addToRow(derivativesHySP, 4, shortPeriodJacobian);
316 addToRow(derivativesLSP, 5, shortPeriodJacobian);
317
318
319 double[][] shortPeriodJacobianRef = new double[6][6];
320 double dP = 0.001;
321 double[] steps = ToleranceProvider.getDefaultToleranceProvider(1000000 * dP).getTolerances(orbit, orbitType)[0];
322 for (int i = 0; i < 6; i++) {
323
324 SpacecraftState stateM4 = shiftState(meanState, orbitType, -4 * steps[i], i);
325 double[] shortPeriodM4 = computeShortPeriodTerms(stateM4, srp);
326
327 SpacecraftState stateM3 = shiftState(meanState, orbitType, -3 * steps[i], i);
328 double[] shortPeriodM3 = computeShortPeriodTerms(stateM3, srp);
329
330 SpacecraftState stateM2 = shiftState(meanState, orbitType, -2 * steps[i], i);
331 double[] shortPeriodM2 = computeShortPeriodTerms(stateM2, srp);
332
333 SpacecraftState stateM1 = shiftState(meanState, orbitType, -1 * steps[i], i);
334 double[] shortPeriodM1 = computeShortPeriodTerms(stateM1, srp);
335
336 SpacecraftState stateP1 = shiftState(meanState, orbitType, 1 * steps[i], i);
337 double[] shortPeriodP1 = computeShortPeriodTerms(stateP1, srp);
338
339 SpacecraftState stateP2 = shiftState(meanState, orbitType, 2 * steps[i], i);
340 double[] shortPeriodP2 = computeShortPeriodTerms(stateP2, srp);
341
342 SpacecraftState stateP3 = shiftState(meanState, orbitType, 3 * steps[i], i);
343 double[] shortPeriodP3 = computeShortPeriodTerms(stateP3, srp);
344
345 SpacecraftState stateP4 = shiftState(meanState, orbitType, 4 * steps[i], i);
346 double[] shortPeriodP4 = computeShortPeriodTerms(stateP4, srp);
347
348 fillJacobianColumn(shortPeriodJacobianRef, i, orbitType, steps[i],
349 shortPeriodM4, shortPeriodM3, shortPeriodM2, shortPeriodM1,
350 shortPeriodP1, shortPeriodP2, shortPeriodP3, shortPeriodP4);
351
352 }
353
354 for (int m = 0; m < 6; ++m) {
355 for (int n = 0; n < 6; ++n) {
356 double error = FastMath.abs((shortPeriodJacobian[m][n] - shortPeriodJacobianRef[m][n]) / shortPeriodJacobianRef[m][n]);
357 Assertions.assertEquals(0, error, 8.3e-10);
358 }
359 }
360
361 }
362
363 @Test
364 void testSRPParametersDerivatives() throws ParseException, IOException {
365 doTestShortPeriodTermsParametersDerivatives(RadiationSensitive.REFLECTION_COEFFICIENT, 9.e-15);
366 }
367
368 @Test
369 void testMuParametersDerivatives() throws ParseException, IOException {
370 doTestShortPeriodTermsParametersDerivatives(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT, 5.e-11);
371 }
372
373 @SuppressWarnings("unchecked")
374 private void doTestShortPeriodTermsParametersDerivatives(String parameterName, double tolerance) {
375
376
377 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
378 TimeScalesFactory.getUTC());
379
380 final Orbit orbit = new EquinoctialOrbit(42164000,
381 10e-3,
382 10e-3,
383 FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
384 FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
385 PositionAngleType.TRUE,
386 FramesFactory.getEME2000(),
387 initDate,
388 3.986004415E14);
389
390 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
391
392 final SpacecraftState meanState = new SpacecraftState(orbit);
393
394
395 final AttitudeProvider attitudeProvider = new LofOffset(meanState.getFrame(),
396 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
397 0.0, 0.0, 0.0);
398
399
400 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
401 final DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
402 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
403 Constants.WGS84_EARTH_FLATTENING,
404 FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
405 provider.getMu());
406 srp.registerAttitudeProvider(attitudeProvider);
407
408 for (final ParameterDriver driver : srp.getParametersDrivers()) {
409 driver.setValue(driver.getReferenceValue());
410 driver.setSelected(driver.getName().equals(parameterName));
411 }
412
413
414 final DSSTGradientConverter converter = new DSSTGradientConverter(meanState, Utils.defaultLaw());
415
416
417 final FieldSpacecraftState<Gradient> dsState = converter.getState(srp);
418
419 final FieldAuxiliaryElements<Gradient> fieldAuxiliaryElements = new FieldAuxiliaryElements<>(dsState.getOrbit(), 1);
420
421
422 final Gradient zero = dsState.getDate().getField().getZero();
423
424
425 final List<FieldShortPeriodTerms<Gradient>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<Gradient>>();
426 shortPeriodTerms.addAll(srp.initializeShortPeriodTerms(fieldAuxiliaryElements, PropagationType.OSCULATING, converter.getParametersAtStateDate(dsState, srp)));
427 srp.updateShortPeriodTerms(converter.getParameters(dsState, srp), dsState);
428 final Gradient[] shortPeriod = new Gradient[6];
429 Arrays.fill(shortPeriod, zero);
430 for (final FieldShortPeriodTerms<Gradient> spt : shortPeriodTerms) {
431 final Gradient[] spVariation = spt.value(dsState.getOrbit());
432 for (int i = 0; i < spVariation .length; i++) {
433 shortPeriod[i] = shortPeriod[i].add(spVariation[i]);
434 }
435 }
436
437 final double[][] shortPeriodJacobian = new double[6][1];
438
439 final double[] derivativesASP = shortPeriod[0].getGradient();
440 final double[] derivativesExSP = shortPeriod[1].getGradient();
441 final double[] derivativesEySP = shortPeriod[2].getGradient();
442 final double[] derivativesHxSP = shortPeriod[3].getGradient();
443 final double[] derivativesHySP = shortPeriod[4].getGradient();
444 final double[] derivativesLSP = shortPeriod[5].getGradient();
445
446 int index = converter.getFreeStateParameters();
447 for (ParameterDriver driver : srp.getParametersDrivers()) {
448 if (driver.isSelected()) {
449 shortPeriodJacobian[0][0] += derivativesASP[index];
450 shortPeriodJacobian[1][0] += derivativesExSP[index];
451 shortPeriodJacobian[2][0] += derivativesEySP[index];
452 shortPeriodJacobian[3][0] += derivativesHxSP[index];
453 shortPeriodJacobian[4][0] += derivativesHySP[index];
454 shortPeriodJacobian[5][0] += derivativesLSP[index];
455 ++index;
456 }
457 }
458
459
460 double[][] shortPeriodJacobianRef = new double[6][1];
461 ParameterDriversList bound = new ParameterDriversList();
462 for (final ParameterDriver driver : srp.getParametersDrivers()) {
463 if (driver.getName().equals(parameterName)) {
464 driver.setSelected(true);
465 bound.add(driver);
466 } else {
467 driver.setSelected(false);
468 }
469 }
470
471 ParameterDriver selected = bound.getDrivers().get(0);
472 double p0 = selected.getReferenceValue();
473 double h = selected.getScale();
474
475 selected.setValue(p0 - 4 * h);
476 final double[] shortPeriodM4 = computeShortPeriodTerms(meanState, srp);
477
478 selected.setValue(p0 - 3 * h);
479 final double[] shortPeriodM3 = computeShortPeriodTerms(meanState, srp);
480
481 selected.setValue(p0 - 2 * h);
482 final double[] shortPeriodM2 = computeShortPeriodTerms(meanState, srp);
483
484 selected.setValue(p0 - 1 * h);
485 final double[] shortPeriodM1 = computeShortPeriodTerms(meanState, srp);
486
487 selected.setValue(p0 + 1 * h);
488 final double[] shortPeriodP1 = computeShortPeriodTerms(meanState, srp);
489
490 selected.setValue(p0 + 2 * h);
491 final double[] shortPeriodP2 = computeShortPeriodTerms(meanState, srp);
492
493 selected.setValue(p0 + 3 * h);
494 final double[] shortPeriodP3 = computeShortPeriodTerms(meanState, srp);
495
496 selected.setValue(p0 + 4 * h);
497 final double[] shortPeriodP4 = computeShortPeriodTerms(meanState, srp);
498
499 fillJacobianColumn(shortPeriodJacobianRef, 0, orbitType, h,
500 shortPeriodM4, shortPeriodM3, shortPeriodM2, shortPeriodM1,
501 shortPeriodP1, shortPeriodP2, shortPeriodP3, shortPeriodP4);
502
503 for (int i = 0; i < 6; ++i) {
504 Assertions.assertEquals(shortPeriodJacobianRef[i][0],
505 shortPeriodJacobian[i][0],
506 FastMath.abs(shortPeriodJacobianRef[i][0] * tolerance));
507 }
508
509 }
510
511 private double[] computeShortPeriodTerms(SpacecraftState state,
512 DSSTForceModel force) {
513
514 AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), 1);
515
516 List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
517 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(auxiliaryElements, PropagationType.OSCULATING, force.getParameters(state.getDate())));
518 force.updateShortPeriodTerms(force.getParametersAllValues(), state);
519
520 double[] shortPeriod = new double[6];
521 for (ShortPeriodTerms spt : shortPeriodTerms) {
522 double[] spVariation = spt.value(state.getOrbit());
523 for (int i = 0; i < spVariation.length; i++) {
524 shortPeriod[i] += spVariation[i];
525 }
526 }
527
528 return shortPeriod;
529
530 }
531
532 private void fillJacobianColumn(double[][] jacobian, int column,
533 OrbitType orbitType, double h,
534 double[] M4h, double[] M3h,
535 double[] M2h, double[] M1h,
536 double[] P1h, double[] P2h,
537 double[] P3h, double[] P4h) {
538 for (int i = 0; i < jacobian.length; ++i) {
539 jacobian[i][column] = ( -3 * (P4h[i] - M4h[i]) +
540 32 * (P3h[i] - M3h[i]) -
541 168 * (P2h[i] - M2h[i]) +
542 672 * (P1h[i] - M1h[i])) / (840 * h);
543 }
544 }
545
546 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
547 double delta, int column) {
548
549 double[][] array = stateToArray(state, orbitType);
550 array[0][column] += delta;
551
552 return arrayToState(array, orbitType, state.getFrame(), state.getDate(),
553 state.getOrbit().getMu(), state.getAttitude());
554
555 }
556
557 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
558 double[][] array = new double[2][6];
559
560 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
561 return array;
562 }
563
564 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
565 Frame frame, AbsoluteDate date, double mu,
566 Attitude attitude) {
567 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngleType.MEAN, date, mu, frame);
568 return new SpacecraftState(orbit, attitude);
569 }
570
571
572
573
574
575
576 private void addToRow(final double[] derivatives, final int index,
577 final double[][] jacobian) {
578
579 for (int i = 0; i < 6; i++) {
580 jacobian[index][i] += derivatives[i];
581 }
582
583 }
584
585 @BeforeEach
586 public void setUp() throws IOException, ParseException {
587 Utils.setDataRoot("regular-data:potential/shm-format");
588 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
589 }
590
591 }