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.Collection;
23 import java.util.Collections;
24 import java.util.HashSet;
25 import java.util.List;
26 import java.util.Set;
27 import java.util.stream.Stream;
28
29 import org.hamcrest.MatcherAssert;
30 import org.hipparchus.CalculusFieldElement;
31 import org.hipparchus.Field;
32 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
33 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
34 import org.hipparchus.geometry.euclidean.threed.Vector3D;
35 import org.hipparchus.ode.FieldODEIntegrator;
36 import org.hipparchus.ode.events.Action;
37 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
38 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
39 import org.hipparchus.ode.nonstiff.DormandPrince54FieldIntegrator;
40 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
41 import org.hipparchus.util.Binary64Field;
42 import org.hipparchus.util.FastMath;
43 import org.hipparchus.util.MathArrays;
44 import org.hipparchus.util.MathUtils;
45 import org.hipparchus.util.Precision;
46 import org.junit.jupiter.api.Assertions;
47 import org.junit.jupiter.api.BeforeEach;
48 import org.junit.jupiter.api.Test;
49 import org.orekit.OrekitMatchers;
50 import org.orekit.Utils;
51 import org.orekit.attitudes.AttitudeProvider;
52 import org.orekit.attitudes.FrameAlignedProvider;
53 import org.orekit.attitudes.LofOffset;
54 import org.orekit.bodies.CelestialBody;
55 import org.orekit.bodies.CelestialBodyFactory;
56 import org.orekit.bodies.OneAxisEllipsoid;
57 import org.orekit.errors.OrekitException;
58 import org.orekit.forces.BoxAndSolarArraySpacecraft;
59 import org.orekit.forces.ForceModel;
60 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
61 import org.orekit.forces.gravity.potential.GRGSFormatReader;
62 import org.orekit.forces.gravity.potential.GravityFieldFactory;
63 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
64 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
65 import org.orekit.forces.gravity.potential.TideSystem;
66 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
67 import org.orekit.frames.Frame;
68 import org.orekit.frames.FramesFactory;
69 import org.orekit.frames.LOFType;
70 import org.orekit.models.earth.atmosphere.Atmosphere;
71 import org.orekit.models.earth.atmosphere.HarrisPriester;
72 import org.orekit.orbits.FieldCartesianOrbit;
73 import org.orekit.orbits.FieldCircularOrbit;
74 import org.orekit.orbits.FieldEquinoctialOrbit;
75 import org.orekit.orbits.FieldKeplerianOrbit;
76 import org.orekit.orbits.FieldOrbit;
77 import org.orekit.orbits.KeplerianOrbit;
78 import org.orekit.orbits.OrbitType;
79 import org.orekit.orbits.PositionAngleType;
80 import org.orekit.propagation.FieldBoundedPropagator;
81 import org.orekit.propagation.FieldEphemerisGenerator;
82 import org.orekit.propagation.FieldPropagator;
83 import org.orekit.propagation.FieldSpacecraftState;
84 import org.orekit.propagation.PropagationType;
85 import org.orekit.propagation.SpacecraftState;
86 import org.orekit.propagation.ToleranceProvider;
87 import org.orekit.propagation.events.EventDetector;
88 import org.orekit.propagation.events.FieldAltitudeDetector;
89 import org.orekit.propagation.events.FieldDateDetector;
90 import org.orekit.propagation.events.FieldEventDetector;
91 import org.orekit.propagation.events.FieldLatitudeCrossingDetector;
92 import org.orekit.propagation.events.handlers.FieldEventHandler;
93 import org.orekit.propagation.numerical.FieldNumericalPropagator;
94 import org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution;
95 import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
96 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
97 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
98 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
99 import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
100 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
101 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
102 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
103 import org.orekit.time.AbsoluteDate;
104 import org.orekit.time.DateComponents;
105 import org.orekit.time.FieldAbsoluteDate;
106 import org.orekit.time.TimeComponents;
107 import org.orekit.time.TimeScale;
108 import org.orekit.time.TimeScalesFactory;
109 import org.orekit.utils.Constants;
110 import org.orekit.utils.FieldPVCoordinates;
111 import org.orekit.utils.IERSConventions;
112 import org.orekit.utils.ParameterDriver;
113 import org.orekit.utils.TimeStampedFieldPVCoordinates;
114
115 public class FieldDSSTPropagatorTest {
116
117
118
119
120
121
122 @Test
123 public void testIssue1029() {
124 doTestIssue1029(Binary64Field.getInstance());
125 }
126
127 private <T extends CalculusFieldElement<T>> void doTestIssue1029(Field<T> field) {
128
129 Utils.setDataRoot("regular-data:potential/grgs-format");
130 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
131
132
133 final T zero = field.getZero();
134
135
136 final FieldAbsoluteDate<T> orbitEpoch = new FieldAbsoluteDate<>(field, 2023, 2, 18, TimeScalesFactory.getUTC());
137 final Frame inertial = FramesFactory.getCIRF(IERSConventions.IERS_2010, true);
138 final FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(42166000.0), zero.add(0.00028), zero.add(FastMath.toRadians(0.05)),
139 zero.add(FastMath.toRadians(66.0)), zero.add(FastMath.toRadians(270.0)),
140 zero.add(FastMath.toRadians(11.94)), PositionAngleType.MEAN,
141 inertial, orbitEpoch, zero.add(Constants.WGS84_EARTH_MU));
142 final FieldEquinoctialOrbit<T> equinoctial = new FieldEquinoctialOrbit<>(orbit);
143
144
145 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-3).getTolerances(equinoctial, OrbitType.EQUINOCTIAL);
146 final AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 3600.0, 86400.0, tol[0], tol[1]);
147 final FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.OSCULATING);
148
149
150 final Frame ecefFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
151 final UnnormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getUnnormalizedProvider(8, 8);
152 propagator.addForceModel(new DSSTZonal(gravityProvider));
153 propagator.addForceModel(new DSSTTesseral(ecefFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, gravityProvider));
154 propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun(), gravityProvider.getMu()));
155 propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon(), gravityProvider.getMu()));
156
157
158 propagator.setInitialState(new FieldSpacecraftState<>(equinoctial, zero.add(6000.0)), PropagationType.MEAN);
159 FieldSpacecraftState<T> propagated = propagator.propagate(orbitEpoch.shiftedBy(20.0 * Constants.JULIAN_DAY));
160
161
162 Assertions.assertEquals(0.0, propagated.getDate().durationFrom(orbitEpoch.shiftedBy(20.0 * Constants.JULIAN_DAY)).getReal(), Double.MIN_VALUE);
163 Assertions.assertEquals(4.216464862956647E7, propagated.getOrbit().getA().getReal(), Double.MIN_VALUE);
164
165 }
166
167 @Test
168 public void testIssue363() {
169 doTestIssue363(Binary64Field.getInstance());
170 }
171
172 private <T extends CalculusFieldElement<T>> void doTestIssue363(Field<T> field) {
173 Utils.setDataRoot("regular-data");
174 final T zero = field.getZero();
175 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
176 FieldCircularOrbit<T> orbit = new FieldCircularOrbit<>(zero.add(7389068.5), zero.add(1.0e-15), zero.add(1.0e-15), zero.add(1.709573), zero.add(1.308398), zero.add(0), PositionAngleType.MEAN,
177 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
178 date, zero.add(Constants.WGS84_EARTH_MU));
179 FieldSpacecraftState<T> osculatingState = new FieldSpacecraftState<>(orbit, zero.add(1116.2829));
180
181 List<DSSTForceModel> dsstForceModels = new ArrayList<>();
182
183 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu().getReal()));
184 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu().getReal()));
185
186 FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
187 Assertions.assertEquals( 0.421, osculatingState.getOrbit().getA().subtract(meanState.getOrbit().getA()).getReal(), 1.0e-3);
188 Assertions.assertEquals(-5.23e-8, osculatingState.getOrbit().getEquinoctialEx().subtract(meanState.getOrbit().getEquinoctialEx()).getReal(), 1.0e-10);
189 Assertions.assertEquals(15.22e-8, osculatingState.getOrbit().getEquinoctialEy().subtract(meanState.getOrbit().getEquinoctialEy()).getReal(), 1.0e-10);
190 Assertions.assertEquals(-3.15e-8, osculatingState.getOrbit().getHx().subtract(meanState.getOrbit().getHx()).getReal(), 1.0e-10);
191 Assertions.assertEquals( 2.83e-8, osculatingState.getOrbit().getHy().subtract(meanState.getOrbit().getHy()).getReal(), 1.0e-10);
192 Assertions.assertEquals(15.96e-8, osculatingState.getOrbit().getLM().subtract(meanState.getOrbit().getLM()).getReal(), 1.0e-10);
193
194 }
195
196 @Test
197 public void testIssue364() {
198 doTestIssue364(Binary64Field.getInstance());
199 }
200
201 private <T extends CalculusFieldElement<T>> void doTestIssue364(Field<T> field) {
202 Utils.setDataRoot("regular-data");
203 final T zero = field.getZero();
204 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, "2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
205 FieldCircularOrbit<T> orbit = new FieldCircularOrbit<>(zero.add(7389068.5), zero.add(0.0), zero.add(0.0), zero.add(1.709573), zero.add(1.308398), zero.add(0), PositionAngleType.MEAN,
206 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
207 date, zero.add(Constants.WGS84_EARTH_MU));
208 FieldSpacecraftState<T> osculatingState = new FieldSpacecraftState<>(orbit, zero.add(1116.2829));
209
210 List<DSSTForceModel> dsstForceModels = new ArrayList<>();
211
212 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu().getReal()));
213 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu().getReal()));
214
215 FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
216 Assertions.assertEquals( 0.421, osculatingState.getOrbit().getA().subtract(meanState.getOrbit().getA()).getReal(), 1.0e-3);
217 Assertions.assertEquals(-5.23e-8, osculatingState.getOrbit().getEquinoctialEx().subtract(meanState.getOrbit().getEquinoctialEx()).getReal(), 1.0e-10);
218 Assertions.assertEquals(15.22e-8, osculatingState.getOrbit().getEquinoctialEy().subtract(meanState.getOrbit().getEquinoctialEy()).getReal(), 1.0e-10);
219 Assertions.assertEquals(-3.15e-8, osculatingState.getOrbit().getHx().subtract(meanState.getOrbit().getHx()).getReal(), 1.0e-10);
220 Assertions.assertEquals( 2.83e-8, osculatingState.getOrbit().getHy().subtract(meanState.getOrbit().getHy()).getReal(), 1.0e-10);
221 Assertions.assertEquals(15.96e-8, osculatingState.getOrbit().getLM().subtract(meanState.getOrbit().getLM()).getReal(), 1.0e-10);
222
223 }
224
225 @Test
226 public void testHighDegreesSetting() {
227 doTestHighDegreesSetting(Binary64Field.getInstance());
228 }
229
230 private <T extends CalculusFieldElement<T>> void doTestHighDegreesSetting(Field<T> field) {
231
232 final T zero = field.getZero();
233 Utils.setDataRoot("regular-data:potential/grgs-format");
234 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
235 int earthDegree = 36;
236 int earthOrder = 36;
237 int eccPower = 4;
238 final UnnormalizedSphericalHarmonicsProvider provider =
239 GravityFieldFactory.getUnnormalizedProvider(earthDegree, earthOrder);
240 final org.orekit.frames.Frame earthFrame =
241 FramesFactory.getITRF(IERSConventions.IERS_2010, true);
242 final DSSTForceModel force =
243 new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
244 earthDegree, earthOrder, eccPower, earthDegree + eccPower,
245 earthDegree, earthOrder, eccPower);
246 final Collection<DSSTForceModel> forces = new ArrayList<>();
247 forces.add(force);
248 TimeScale tai = TimeScalesFactory.getTAI();
249 FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, "2015-07-01", tai);
250 Frame eci = FramesFactory.getGCRF();
251 FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(
252 zero.add(7120000.0), zero.add(1.0e-3), zero.add(FastMath.toRadians(60.0)),
253 zero.add(FastMath.toRadians(120.0)), zero.add(FastMath.toRadians(47.0)),
254 zero.add(FastMath.toRadians(12.0)),
255 PositionAngleType.TRUE, eci, initialDate, zero.add(Constants.EIGEN5C_EARTH_MU));
256 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
257 FieldSpacecraftState<T> oscuState = FieldDSSTPropagator.computeOsculatingState(state, null, forces);
258 Assertions.assertEquals(7119927.148, oscuState.getOrbit().getA().getReal(), 0.001);
259 }
260
261 @Test
262 public void testEphemerisDates() {
263 doTestEphemerisDates(Binary64Field.getInstance());
264 }
265
266 private <T extends CalculusFieldElement<T>> void doTestEphemerisDates(Field<T> field) {
267 final T zero = field.getZero();
268
269 TimeScale tai = TimeScalesFactory.getTAI();
270 FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, "2015-07-01", tai);
271 FieldAbsoluteDate<T> startDate = new FieldAbsoluteDate<>(field, "2015-07-03", tai).shiftedBy(-0.1);
272 FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(field, "2015-07-04", tai);
273 Frame eci = FramesFactory.getGCRF();
274 FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(
275 zero.add(600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS), zero, zero, zero, zero, zero,
276 PositionAngleType.TRUE, eci, initialDate, zero.add(Constants.EIGEN5C_EARTH_MU));
277 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(orbit, OrbitType.EQUINOCTIAL);
278 FieldPropagator<T> prop = new FieldDSSTPropagator<>(field,
279 new DormandPrince853FieldIntegrator<>(field, 0.1, 500, tol[0], tol[1]));
280 prop.resetInitialState(new FieldSpacecraftState<>(new FieldCartesianOrbit<>(orbit)));
281
282
283 final FieldEphemerisGenerator<T> generator = prop.getEphemerisGenerator();
284 prop.propagate(startDate, endDate);
285 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
286
287
288 TimeStampedFieldPVCoordinates<T> actualPV = ephemeris.getPVCoordinates(startDate, eci);
289 TimeStampedFieldPVCoordinates<T> expectedPV = orbit.getPVCoordinates(startDate, eci);
290 MatcherAssert.assertThat(actualPV.getPosition().toVector3D(),
291 OrekitMatchers.vectorCloseTo(expectedPV.getPosition().toVector3D(), 1.0));
292 MatcherAssert.assertThat(actualPV.getVelocity().toVector3D(),
293 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity().toVector3D(), 1.0));
294 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate).getReal(),
295 OrekitMatchers.closeTo(0, 0));
296 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate).getReal(),
297 OrekitMatchers.closeTo(0, 0));
298
299 FieldAbsoluteDate<T> date = endDate.shiftedBy(-0.11);
300 Assertions.assertEquals(
301 ephemeris.propagate(date).getDate().durationFrom(date).getReal(), 0, 0);
302 }
303
304 @Test
305 public void testNoExtrapolation() {
306 doTestNoExtrapolation(Binary64Field.getInstance());
307 }
308
309 private <T extends CalculusFieldElement<T>> void doTestNoExtrapolation(Field<T> field) {
310 FieldSpacecraftState<T> state = getLEOState(field);
311 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, state);
312
313
314 final FieldSpacecraftState<T> finalState = dsstPropagator.propagate(state.getDate());
315
316
317 final FieldVector3D<T> initialPosition = state.getPosition();
318 final FieldVector3D<T> initialVelocity = state.getPVCoordinates().getVelocity();
319
320
321 final FieldVector3D<T> finalPosition = finalState.getPosition();
322 final FieldVector3D<T> finalVelocity = finalState.getPVCoordinates().getVelocity();
323
324
325 Assertions.assertEquals(initialPosition.getX().getReal(), finalPosition.getX().getReal(), 0.0);
326 Assertions.assertEquals(initialPosition.getY().getReal(), finalPosition.getY().getReal(), 0.0);
327 Assertions.assertEquals(initialPosition.getZ().getReal(), finalPosition.getZ().getReal(), 0.0);
328 Assertions.assertEquals(initialVelocity.getX().getReal(), finalVelocity.getX().getReal(), 0.0);
329 Assertions.assertEquals(initialVelocity.getY().getReal(), finalVelocity.getY().getReal(), 0.0);
330 Assertions.assertEquals(initialVelocity.getZ().getReal(), finalVelocity.getZ().getReal(), 0.0);
331 }
332
333 @Test
334 public void testKepler() {
335 doTestKepler(Binary64Field.getInstance());
336 }
337
338 private <T extends CalculusFieldElement<T>> void doTestKepler(Field<T> field) {
339 FieldSpacecraftState<T> state = getLEOState(field);
340 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, state);
341
342
343 final double dt = 3200.;
344 final FieldSpacecraftState<T> finalState = dsstPropagator.propagate(state.getDate().shiftedBy(dt));
345
346
347 final T n = FastMath.sqrt(state.getOrbit().getA().reciprocal().multiply(state.getOrbit().getMu())).divide(state.getOrbit().getA());
348 Assertions.assertEquals(state.getOrbit().getA().getReal(), finalState.getOrbit().getA().getReal(), 0.);
349 Assertions.assertEquals(state.getOrbit().getEquinoctialEx().getReal(), finalState.getOrbit().getEquinoctialEx().getReal(), 0.);
350 Assertions.assertEquals(state.getOrbit().getEquinoctialEy().getReal(), finalState.getOrbit().getEquinoctialEy().getReal(), 0.);
351 Assertions.assertEquals(state.getOrbit().getHx().getReal(), finalState.getOrbit().getHx().getReal(), 0.);
352 Assertions.assertEquals(state.getOrbit().getHy().getReal(), finalState.getOrbit().getHy().getReal(), 0.);
353 Assertions.assertEquals(state.getOrbit().getLM().add(n.multiply(dt)).getReal(), finalState.getOrbit().getLM().getReal(), 1.e-14);
354
355 }
356
357 @Test
358 public void testEphemeris() {
359 doTestEphemeris(Binary64Field.getInstance());
360 }
361
362 private <T extends CalculusFieldElement<T>> void doTestEphemeris(Field<T> field) {
363 FieldSpacecraftState<T> state = getGEOState(field);
364 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, state);
365
366
367 final FieldEphemerisGenerator<T> generator = dsstPropagator.getEphemerisGenerator();
368
369
370 final double dt = 2. * Constants.JULIAN_DAY;
371 dsstPropagator.propagate(state.getDate().shiftedBy(5. * dt));
372
373
374 FieldBoundedPropagator<T> ephem = generator.getGeneratedEphemeris();
375
376
377 final FieldSpacecraftState<T> s = ephem.propagate(state.getDate().shiftedBy(dt));
378
379
380 final T n = FastMath.sqrt(state.getOrbit().getA().reciprocal().multiply(state.getOrbit().getMu())).divide(state.getOrbit().getA());
381 Assertions.assertEquals(state.getOrbit().getA().getReal(), s.getOrbit().getA().getReal(), 0.);
382 Assertions.assertEquals(state.getOrbit().getEquinoctialEx().getReal(), s.getOrbit().getEquinoctialEx().getReal(), 0.);
383 Assertions.assertEquals(state.getOrbit().getEquinoctialEy().getReal(), s.getOrbit().getEquinoctialEy().getReal(), 0.);
384 Assertions.assertEquals(state.getOrbit().getHx().getReal(), s.getOrbit().getHx().getReal(), 0.);
385 Assertions.assertEquals(state.getOrbit().getHy().getReal(), s.getOrbit().getHy().getReal(), 0.);
386 Assertions.assertEquals(state.getOrbit().getLM().add(n.multiply(dt)).getReal(), s.getOrbit().getLM().getReal(), 1.5e-14);
387
388 }
389
390 @Test
391 public void testPropagationWithCentralBody() {
392 doTestPropagationWithCentralBody(Binary64Field.getInstance());
393 }
394
395 private <T extends CalculusFieldElement<T>> void doTestPropagationWithCentralBody(Field<T> field) {
396
397 final T zero = field.getZero();
398
399
400 final UnnormalizedSphericalHarmonicsProvider provider =
401 GravityFieldFactory.getUnnormalizedProvider(4, 4);
402 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
403
404
405 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2007, 4, 16, 0, 46, 42.400,
406 TimeScalesFactory.getUTC());
407 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(26559890.),
408 zero.add(0.0041632),
409 zero.add(FastMath.toRadians(55.2)),
410 zero.add(FastMath.toRadians(315.4985)),
411 zero.add(FastMath.toRadians(130.7562)),
412 zero.add(FastMath.toRadians(44.2377)),
413 PositionAngleType.MEAN,
414 FramesFactory.getEME2000(),
415 initDate,
416 zero.add(provider.getMu()));
417
418
419 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, new FieldSpacecraftState<>(orbit));
420 dsstPropagator.addForceModel(new DSSTZonal(provider, 4, 3, 9));
421 dsstPropagator.addForceModel(new DSSTTesseral(earthFrame,
422 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
423 4, 4, 4, 8, 4, 4, 2));
424
425
426 final FieldSpacecraftState<T> state = dsstPropagator.propagate(initDate.shiftedBy(5. * 86400.));
427
428
429
430
431
432
433
434
435 Assertions.assertEquals(2.655992081E7, state.getOrbit().getA().getReal(), 1.e2);
436 Assertions.assertEquals(0.2731622444E-03, state.getOrbit().getEquinoctialEx().getReal(), 2.e-8);
437 Assertions.assertEquals(0.4164167597E-02, state.getOrbit().getEquinoctialEy().getReal(), 2.e-8);
438 Assertions.assertEquals(-0.3399607878, state.getOrbit().getHx().getReal(), 5.e-8);
439 Assertions.assertEquals(0.3971568634, state.getOrbit().getHy().getReal(), 2.e-6);
440 Assertions.assertEquals(140.6375352,
441 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), zero.add(FastMath.PI)).getReal()),
442 5.e-3);
443 }
444
445 @Test
446 public void testPropagationWithThirdBody() throws IOException {
447 doTestPropagationWithThirdBody(Binary64Field.getInstance());
448 }
449
450 private <T extends CalculusFieldElement<T>> void doTestPropagationWithThirdBody(Field<T> field) throws IOException {
451
452 final T zero = field.getZero();
453
454
455 final UnnormalizedSphericalHarmonicsProvider provider =
456 GravityFieldFactory.getUnnormalizedProvider(2, 0);
457 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
458 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
459 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
460 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
461 provider, 2, 0, 0, 2, 2, 0, 0);
462
463
464 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
465 DSSTForceModel sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu());
466
467
468 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2003, 7, 1, 0, 0, 00.000,
469 TimeScalesFactory.getUTC());
470 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(42163393.),
471 zero.add(0.2684),
472 zero.add(FastMath.toRadians(63.435)),
473 zero.add(FastMath.toRadians(270.0)),
474 zero.add(FastMath.toRadians(285.0)),
475 zero.add(FastMath.toRadians(344.0)),
476 PositionAngleType.MEAN,
477 FramesFactory.getEME2000(),
478 initDate,
479 zero.add(provider.getMu()));
480
481
482 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, new FieldSpacecraftState<>(orbit));
483 dsstPropagator.addForceModel(zonal);
484 dsstPropagator.addForceModel(tesseral);
485 dsstPropagator.addForceModel(moon);
486 dsstPropagator.addForceModel(sun);
487
488
489 final FieldSpacecraftState<T> state = dsstPropagator.propagate(initDate.shiftedBy(5. * 86400.));
490
491
492
493
494
495
496
497
498 Assertions.assertEquals(42163393.0, state.getOrbit().getA().getReal(), 1.e-1);
499 Assertions.assertEquals(-0.2592789733084587, state.getOrbit().getEquinoctialEx().getReal(), 5.e-7);
500 Assertions.assertEquals(-0.06893353670734315, state.getOrbit().getEquinoctialEy().getReal(), 2.e-7);
501 Assertions.assertEquals( 0.1595005111738418, state.getOrbit().getHx().getReal(), 2.e-7);
502 Assertions.assertEquals(-0.5968524904937771, state.getOrbit().getHy().getReal(), 5.e-8);
503 Assertions.assertEquals(183.9386620425922,
504 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), zero.add(FastMath.PI)).getReal()),
505 3.e-2);
506 }
507
508 @Test
509 public void testTooSmallMaxDegree() {
510 Assertions.assertThrows(OrekitException.class, () -> {
511 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 1, 0, 3);
512 });
513 }
514
515 @Test
516 public void testTooLargeMaxDegree() {
517 Assertions.assertThrows(OrekitException.class, () -> {
518 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 8, 0, 8);
519 });
520 }
521
522 @Test
523 public void testWrongMaxPower() {
524 Assertions.assertThrows(OrekitException.class, () -> {
525 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(8, 8), 4, 4, 4);
526 });
527 }
528
529 @Test
530 public void testPropagationWithDrag() {
531 doTestPropagationWithDrag(Binary64Field.getInstance());
532 }
533
534 private <T extends CalculusFieldElement<T>> void doTestPropagationWithDrag(Field<T> field) {
535
536 final T zero = field.getZero();
537
538 final UnnormalizedSphericalHarmonicsProvider provider =
539 GravityFieldFactory.getUnnormalizedProvider(2, 0);
540 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
541 DSSTForceModel zonal = new DSSTZonal(provider, 2, 0, 5);
542 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
543 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
544 provider, 2, 0, 0, 2, 2, 0, 0);
545
546
547 final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
548 Constants.WGS84_EARTH_FLATTENING,
549 earthFrame);
550 final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
551
552 final double cd = 2.0;
553 final double area = 25.0;
554 DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area, provider.getMu());
555
556
557 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2003, 7, 1, 0, 0, 00.000,
558 TimeScalesFactory.getUTC());
559 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(7204535.848109440),
560 zero.add(0.0012402238462686),
561 zero.add(FastMath.toRadians(98.74341600466740)),
562 zero.add(FastMath.toRadians(111.1990175076630)),
563 zero.add(FastMath.toRadians(43.32990110790340)),
564 zero.add(FastMath.toRadians(68.66852509725620)),
565 PositionAngleType.MEAN,
566 FramesFactory.getEME2000(),
567 initDate,
568 zero.add(provider.getMu()));
569
570
571 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, new FieldSpacecraftState<>(orbit));
572 dsstPropagator.addForceModel(zonal);
573 dsstPropagator.addForceModel(tesseral);
574 dsstPropagator.addForceModel(drag);
575
576
577 final FieldSpacecraftState<T> state = dsstPropagator.propagate(initDate.shiftedBy(5. * 86400.));
578
579
580
581
582
583
584
585
586 Assertions.assertEquals(7204521.657141485, state.getOrbit().getA().getReal(), 6.e-1);
587 Assertions.assertEquals(-0.001016800430994036, state.getOrbit().getEquinoctialEx().getReal(), 5.e-8);
588 Assertions.assertEquals(0.0007093755541595772, state.getOrbit().getEquinoctialEy().getReal(), 2.e-8);
589 Assertions.assertEquals(0.7757573478894775, state.getOrbit().getHx().getReal(), 5.e-8);
590 Assertions.assertEquals(0.8698955648709271, state.getOrbit().getHy().getReal(), 5.e-8);
591 Assertions.assertEquals(193.0939742953394,
592 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), zero.add(FastMath.PI)).getReal()),
593 2.e-3);
594
595
596 Assertions.assertEquals(((DSSTAtmosphericDrag)drag).getAtmosphere(), atm);
597
598 final double atmosphericMaxConstant = 1000000.0;
599 Assertions.assertEquals(((DSSTAtmosphericDrag)drag).getRbar(), atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 1e-9);
600 }
601
602 @Test
603 public void testPropagationWithSolarRadiationPressure() {
604 doTestPropagationWithSolarRadiationPressure(Binary64Field.getInstance());
605 }
606
607 private <T extends CalculusFieldElement<T>> void doTestPropagationWithSolarRadiationPressure(Field<T> field) {
608
609 final T zero = field.getZero();
610
611 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
612 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
613 DSSTForceModel tesseral = new DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(),
614 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
615 provider, 2, 0, 0, 2, 2, 0, 0);
616
617
618 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
619 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
620 Constants.WGS84_EARTH_FLATTENING,
621 FramesFactory.getITRF(IERSConventions.IERS_2010, false)),
622 provider.getMu());
623
624
625 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, 2003, 9, 16, 0, 0, 00.000,
626 TimeScalesFactory.getUTC());
627 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(42166258.),
628 zero.add(0.0001),
629 zero.add(FastMath.toRadians(0.001)),
630 zero.add(FastMath.toRadians(315.4985)),
631 zero.add(FastMath.toRadians(130.7562)),
632 zero.add(FastMath.toRadians(44.2377)),
633 PositionAngleType.MEAN,
634 FramesFactory.getGCRF(),
635 initDate,
636 zero.add(provider.getMu()));
637
638
639 final FieldDSSTPropagator<T> dsstPropagatorp = new FieldDSSTPropagator<>(field, new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(86400.)));
640 dsstPropagatorp.setInitialState(new FieldSpacecraftState<>(orbit), PropagationType.MEAN);
641 dsstPropagatorp.addForceModel(zonal);
642 dsstPropagatorp.addForceModel(tesseral);
643 dsstPropagatorp.addForceModel(srp);
644
645
646 final FieldSpacecraftState<T> state = dsstPropagatorp.propagate(initDate.shiftedBy(10. * 86400.));
647
648
649
650
651
652
653
654
655 Assertions.assertEquals(42166257.99807995, state.getOrbit().getA().getReal(), 1.2);
656 Assertions.assertEquals(-0.1781865038201885e-05, state.getOrbit().getEquinoctialEx().getReal(), 3.e-7);
657 Assertions.assertEquals(-0.1191876027555493e-03, state.getOrbit().getEquinoctialEy().getReal(), 4.e-6);
658 Assertions.assertEquals(-0.5624363171289686e-05, state.getOrbit().getHx().getReal(), 4.e-9);
659 Assertions.assertEquals( 0.6618387121369373e-05, state.getOrbit().getHy().getReal(), 3.e-10);
660 Assertions.assertEquals(140.3496229467104,
661 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), zero.add(FastMath.PI)).getReal()),
662 2.e-4);
663 }
664
665 @Test
666 public void testStopEvent() {
667 doTestStopEvent(Binary64Field.getInstance());
668 }
669
670 private <T extends CalculusFieldElement<T>> void doTestStopEvent(Field<T> field) {
671 FieldSpacecraftState<T> state = getLEOState(field);
672 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, state);
673
674 final FieldAbsoluteDate<T> stopDate = state.getDate().shiftedBy(1000);
675 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<>(Action.STOP);
676 FieldDateDetector<T> detector = new FieldDateDetector<>(field, stopDate).withHandler(checking);
677 dsstPropagator.addEventDetector(detector);
678 checking.assertEvent(false);
679 final FieldSpacecraftState<T> finalState = dsstPropagator.propagate(state.getDate().shiftedBy(3200));
680 checking.assertEvent(true);
681 Assertions.assertEquals(0, finalState.getDate().durationFrom(stopDate).getReal(), 1.0e-10);
682 }
683
684 @Test
685 public void testContinueEvent() {
686 doTestContinueEvent(Binary64Field.getInstance());
687 }
688
689 private <T extends CalculusFieldElement<T>> void doTestContinueEvent(Field<T> field) {
690 FieldSpacecraftState<T> state = getLEOState(field);
691 final FieldDSSTPropagator<T> dsstPropagator = setDSSTProp(field, state);
692
693 final FieldAbsoluteDate<T> resetDate = state.getDate().shiftedBy(1000);
694 CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<>(Action.CONTINUE);
695 FieldDateDetector<T> detector = new FieldDateDetector<>(field, resetDate).withHandler(checking);
696 dsstPropagator.addEventDetector(detector);
697 final double dt = 3200;
698 checking.assertEvent(false);
699 final FieldSpacecraftState<T> finalState = dsstPropagator.propagate(state.getDate().shiftedBy(dt));
700 checking.assertEvent(true);
701 final T n = FastMath.sqrt(state.getOrbit().getA().reciprocal().multiply(state.getOrbit().getMu())).divide(state.getOrbit().getA());
702 Assertions.assertEquals(state.getOrbit().getA().getReal(), finalState.getOrbit().getA().getReal(), 1.0e-10);
703 Assertions.assertEquals(state.getOrbit().getEquinoctialEx().getReal(), finalState.getOrbit().getEquinoctialEx().getReal(), 1.0e-10);
704 Assertions.assertEquals(state.getOrbit().getEquinoctialEy().getReal(), finalState.getOrbit().getEquinoctialEy().getReal(), 1.0e-10);
705 Assertions.assertEquals(state.getOrbit().getHx().getReal(), finalState.getOrbit().getHx().getReal(), 1.0e-10);
706 Assertions.assertEquals(state.getOrbit().getHy().getReal(), finalState.getOrbit().getHy().getReal(), 1.0e-10);
707 Assertions.assertEquals(state.getOrbit().getLM().add(n.multiply(dt)).getReal(), finalState.getOrbit().getLM().getReal(), 6.0e-10);
708 }
709
710 @Test
711 public void testIssue157() {
712 doTestIssue157(Binary64Field.getInstance());
713 }
714
715 private <T extends CalculusFieldElement<T>> void doTestIssue157(Field<T> field) {
716 Utils.setDataRoot("regular-data:potential/icgem-format");
717 final T zero = field.getZero();
718 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
719 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
720 FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(13378000), zero.add(0.05), zero.add(0), zero.add(0), zero.add(FastMath.PI), zero.add(0), PositionAngleType.MEAN,
721 FramesFactory.getTOD(false),
722 new FieldAbsoluteDate<>(field, 2003, 5, 6, TimeScalesFactory.getUTC()),
723 zero.add(nshp.getMu()));
724 T period = orbit.getKeplerianPeriod();
725 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(orbit, OrbitType.EQUINOCTIAL);
726 AdaptiveStepsizeFieldIntegrator<T> integrator =
727 new DormandPrince853FieldIntegrator<>(field, period.getReal() / 100, period.getReal() * 100, tolerance[0], tolerance[1]);
728 integrator.setInitialStepSize(period.multiply(10.).getReal());
729 FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.MEAN);
730 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
731 Constants.WGS84_EARTH_FLATTENING,
732 FramesFactory.getGTOD(false));
733 CelestialBody sun = CelestialBodyFactory.getSun();
734 CelestialBody moon = CelestialBodyFactory.getMoon();
735 propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
736 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
737 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
738 nshp, 8, 8, 4, 12, 8, 8, 4));
739 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
740 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
741 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
742 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
743
744
745 propagator.setInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)), PropagationType.OSCULATING);
746 FieldSpacecraftState<T> finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
747
748
749
750
751 Assertions.assertEquals(2187.2, orbit.getA().subtract(finalState.getOrbit().getA()).getReal(), 1.0);
752
753 propagator.setInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)), PropagationType.MEAN);
754 finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
755
756
757 Assertions.assertEquals(1475.90, orbit.getA().subtract(finalState.getOrbit().getA()).getReal(), 1.0);
758
759 }
760
761
762
763
764 @Test
765 public void testEphemerisGeneration() {
766 doTestEphemerisGeneration(Binary64Field.getInstance());
767 }
768
769 private <T extends CalculusFieldElement<T>> void doTestEphemerisGeneration(Field<T> field){
770
771
772
773
774 Utils.setDataRoot("regular-data:potential/icgem-format");
775 final T zero = field.getZero();
776 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
777 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
778 FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(13378000), zero.add(0.05), zero.add(0), zero.add(0), zero.add(FastMath.PI), zero.add(0), PositionAngleType.MEAN,
779 FramesFactory.getTOD(false),
780 new FieldAbsoluteDate<>(field, 2003, 5, 6, TimeScalesFactory.getUTC()),
781 zero.add(nshp.getMu()));
782 T period = orbit.getKeplerianPeriod();
783 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(orbit, OrbitType.EQUINOCTIAL);
784 AdaptiveStepsizeFieldIntegrator<T> integrator =
785 new DormandPrince853FieldIntegrator<>(field, period.getReal() / 100, period.getReal() * 100, tolerance[0], tolerance[1]);
786 integrator.setInitialStepSize(period.multiply(10.).getReal());
787 FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.OSCULATING);
788 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
789 Constants.WGS84_EARTH_FLATTENING,
790 FramesFactory.getGTOD(false));
791 CelestialBody sun = CelestialBodyFactory.getSun();
792 CelestialBody moon = CelestialBodyFactory.getMoon();
793 propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
794 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
795 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
796 nshp, 8, 8, 4, 12, 8, 8, 4));
797 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
798 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
799 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
800 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
801 propagator.setInterpolationGridToMaxTimeGap(zero.add(0.5 * Constants.JULIAN_DAY));
802
803
804
805
806
807
808 final double nDays = 5.;
809
810
811 propagator.setInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)), PropagationType.MEAN);
812 final List<FieldSpacecraftState<T>> states = new ArrayList<>();
813 propagator.setStepHandler(zero.add(600), states::add);
814 propagator.propagate(orbit.getDate().shiftedBy(nDays * Constants.JULIAN_DAY));
815
816
817 propagator.setInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)), PropagationType.MEAN);
818 final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
819 propagator.propagate(orbit.getDate().shiftedBy(nDays * Constants.JULIAN_DAY));
820 FieldBoundedPropagator<T> ephemeris = generator.getGeneratedEphemeris();
821
822 T maxError = zero;
823 for (final FieldSpacecraftState<T> state : states) {
824 final FieldSpacecraftState<T> fromEphemeris = ephemeris.propagate(state.getDate());
825 final T error = FieldVector3D.distance(state.getPosition(), fromEphemeris.getPosition());
826 maxError = FastMath.max(maxError, error);
827 }
828
829
830
831
832
833 Assertions.assertEquals(0.0, maxError.getReal(), Precision.SAFE_MIN);
834 }
835
836 @Test
837 public void testGetInitialOsculatingState() {
838 doTestGetInitialOsculatingState(Binary64Field.getInstance());
839 }
840
841 private <T extends CalculusFieldElement<T>> void doTestGetInitialOsculatingState(Field<T> field) {
842 final FieldSpacecraftState<T> initialState = getGEOState(field);
843
844
845 final T minStep = initialState.getOrbit().getKeplerianPeriod().multiply(0.1);
846 final T maxStep = initialState.getOrbit().getKeplerianPeriod().multiply(10.0);
847 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-1).getTolerances(initialState.getOrbit(), OrbitType.EQUINOCTIAL);
848 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep.getReal(), maxStep.getReal(), tol[0], tol[1]);
849
850
851 FieldDSSTPropagator<T> prop = new FieldDSSTPropagator<>(field, integrator, PropagationType.MEAN);
852
853 final UnnormalizedSphericalHarmonicsProvider provider =
854 GravityFieldFactory.getUnnormalizedProvider(4, 0);
855 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
856 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
857 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
858 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
859 provider, 4, 0, 4, 8, 4, 0, 2);
860 prop.addForceModel(zonal);
861 prop.addForceModel(tesseral);
862
863
864 prop.setInitialState(initialState, PropagationType.MEAN);
865
866 Assertions.assertEquals(initialState, prop.getInitialState());
867
868
869 Assertions.assertEquals(initialState, prop.propagate(initialState.getDate()));
870 }
871
872 @Test
873 public void testMeanToOsculatingState() {
874 doTestMeanToOsculatingState(Binary64Field.getInstance());
875 }
876 private <T extends CalculusFieldElement<T>> void doTestMeanToOsculatingState(Field<T> field) {
877 final FieldSpacecraftState<T> meanState = getGEOState(field);
878
879 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
880 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
881 final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
882 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
883 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
884 provider, 2, 0, 0, 2, 2, 0, 0);
885
886 final Collection<DSSTForceModel> forces = new ArrayList<>();
887 forces.add(zonal);
888 forces.add(tesseral);
889
890 final FieldSpacecraftState<T> osculatingState = FieldDSSTPropagator.computeOsculatingState(meanState, null, forces);
891 Assertions.assertEquals(1559.1,
892 FieldVector3D.distance(meanState.getPosition(),
893 osculatingState.getPosition()).getReal(),
894 1.0);
895 }
896
897 @Test
898 public void testOsculatingToMeanState() {
899 doTestOsculatingToMeanState(Binary64Field.getInstance());
900 }
901
902 private <T extends CalculusFieldElement<T>> void doTestOsculatingToMeanState(Field<T> field) {
903 final FieldSpacecraftState<T> meanState = getGEOState(field);
904
905 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
906 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
907
908 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
909 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
910 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
911 provider, 2, 0, 0, 2, 2, 0, 0);
912
913 final Collection<DSSTForceModel> forces = new ArrayList<>();
914 forces.add(zonal);
915 forces.add(tesseral);
916
917 final FieldSpacecraftState<T> osculatingState = FieldDSSTPropagator.computeOsculatingState(meanState, null, forces);
918
919
920 final FieldSpacecraftState<T> computedMeanState = FieldDSSTPropagator.computeMeanState(osculatingState, null, forces);
921
922 Assertions.assertEquals(meanState.getOrbit().getA().getReal(), computedMeanState.getOrbit().getA().getReal(), 2.0e-8);
923 Assertions.assertEquals(0.0,
924 FieldVector3D.distance(meanState.getPosition(),
925 computedMeanState.getPosition()).getReal(),
926 2.0e-8);
927 }
928
929 @Test
930 public void testShortPeriodCoefficients() {
931 doTestShortPeriodCoefficients(Binary64Field.getInstance());
932 }
933
934 private <T extends CalculusFieldElement<T>> void doTestShortPeriodCoefficients(Field<T> field) {
935 Utils.setDataRoot("regular-data:potential/icgem-format");
936 final T zero = field.getZero();
937 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
938 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(4, 4);
939 FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(13378000), zero.add(0.05), zero.add(0), zero.add(0), zero.add(FastMath.PI), zero.add(0), PositionAngleType.MEAN,
940 FramesFactory.getTOD(false),
941 new FieldAbsoluteDate<>(field, 2003, 5, 6, TimeScalesFactory.getUTC()),
942 zero.add(nshp.getMu()));
943 T period = orbit.getKeplerianPeriod();
944 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(orbit, OrbitType.EQUINOCTIAL);
945 AdaptiveStepsizeFieldIntegrator<T> integrator =
946 new DormandPrince853FieldIntegrator<>(field, period.getReal() / 100, period.getReal() * 100, tolerance[0], tolerance[1]);
947 integrator.setInitialStepSize(period.multiply(10).getReal());
948 FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.OSCULATING);
949 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
950 Constants.WGS84_EARTH_FLATTENING,
951 FramesFactory.getGTOD(false));
952 CelestialBody sun = CelestialBodyFactory.getSun();
953 CelestialBody moon = CelestialBodyFactory.getMoon();
954 propagator.addForceModel(new DSSTZonal(nshp, 4, 3, 9));
955 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
956 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
957 nshp, 4, 4, 4, 8, 4, 4, 2));
958 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
959 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
960 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
961 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
962
963 final FieldAbsoluteDate<T> finalDate = orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY);
964 propagator.resetInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)));
965 final FieldSpacecraftState<T> stateNoConfig = propagator.propagate(finalDate);
966 Assertions.assertEquals(0, stateNoConfig.getAdditionalDataValues().size());
967
968 propagator.setSelectedCoefficients(new HashSet<>());
969 propagator.resetInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)));
970 final FieldSpacecraftState<T> stateConfigEmpty = propagator.propagate(finalDate);
971 Assertions.assertEquals(234, stateConfigEmpty.getAdditionalDataValues().size());
972
973 final Set<String> selected = new HashSet<>();
974 selected.add("DSST-3rd-body-Moon-s[7]");
975 selected.add("DSST-central-body-tesseral-c[-2][3]");
976 propagator.setSelectedCoefficients(selected);
977 propagator.resetInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)));
978 final FieldSpacecraftState<T> stateConfigeSelected = propagator.propagate(finalDate);
979 Assertions.assertEquals(selected.size(), stateConfigeSelected.getAdditionalDataValues().size());
980
981 propagator.setSelectedCoefficients(null);
982 propagator.resetInitialState(new FieldSpacecraftState<>(orbit, zero.add(45.0)));
983 final FieldSpacecraftState<T> stateConfigNull = propagator.propagate(finalDate);
984 Assertions.assertEquals(0, stateConfigNull.getAdditionalDataValues().size());
985
986 }
987
988 @Test
989 public void testIssueMeanInclination() {
990 doTestIssueMeanInclination(Binary64Field.getInstance());
991 }
992
993 private <T extends CalculusFieldElement<T>> void doTestIssueMeanInclination(Field<T> field) {
994
995 final T zero = field.getZero();
996 final double earthAe = 6378137.0;
997 final double earthMu = 3.9860044E14;
998 final double earthJ2 = 0.0010826;
999
1000
1001 FieldOrbit<T> orb = new FieldKeplerianOrbit<>(new TimeStampedFieldPVCoordinates<>(new FieldAbsoluteDate<>(field, "1992-10-08T15:20:38.821",
1002 TimeScalesFactory.getUTC()),
1003 new FieldVector3D<>(zero.add(5392808.809823), zero.add(-4187618.3357927715), zero.add(-44206.638015847195)),
1004 new FieldVector3D<>(zero.add(2337.4472786270794), zero.add(2474.0146611860464), zero.add(6778.507766114648)),
1005 FieldVector3D.getZero(field)),
1006 FramesFactory.getTOD(false), zero.add(earthMu));
1007 final FieldSpacecraftState<T> ss = new FieldSpacecraftState<>(orb);
1008 final UnnormalizedSphericalHarmonicsProvider provider =
1009 GravityFieldFactory.getUnnormalizedProvider(earthAe, earthMu, TideSystem.UNKNOWN,
1010 new double[][] { { 0.0 }, { 0.0 }, { -earthJ2 } },
1011 new double[][] { { 0.0 }, { 0.0 }, { 0.0 } });
1012 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
1013 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
1014 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
1015 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
1016 provider, 2, 0, 0, 2, 2, 0, 0);
1017 final Collection<DSSTForceModel> forces = new ArrayList<>();
1018 forces.add(zonal);
1019 forces.add(tesseral);
1020
1021 final FieldOrbit<T> meanOrb = FieldDSSTPropagator.computeMeanState(ss, null, forces).getOrbit();
1022 Assertions.assertEquals(0.0164196, FastMath.toDegrees(orb.getI().subtract(meanOrb.getI()).getReal()), 1.0e-7);
1023 }
1024
1025 @Test
1026 public void testIssue257() {
1027 doTestIssue257(Binary64Field.getInstance());
1028 }
1029
1030 private <T extends CalculusFieldElement<T>> void doTestIssue257(Field<T> field) {
1031 final FieldSpacecraftState<T> meanState = getGEOState(field);
1032
1033
1034 final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), meanState.getOrbit().getMu().getReal());
1035 final DSSTForceModel sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), meanState.getOrbit().getMu().getReal());
1036
1037 final Collection<DSSTForceModel> forces = new ArrayList<>();
1038 forces.add(moon);
1039 forces.add(sun);
1040
1041 final FieldSpacecraftState<T> osculatingState = FieldDSSTPropagator.computeOsculatingState(meanState, null, forces);
1042 Assertions.assertEquals(734.3,
1043 FieldVector3D.distance(meanState.getPosition(),
1044 osculatingState.getPosition()).getReal(),
1045 1.0);
1046
1047 final FieldSpacecraftState<T> computedMeanState = FieldDSSTPropagator.computeMeanState(osculatingState, null, forces);
1048 Assertions.assertEquals(734.3,
1049 FieldVector3D.distance(osculatingState.getPosition(),
1050 computedMeanState.getPosition()).getReal(),
1051 1.0);
1052
1053 Assertions.assertEquals(0.0,
1054 FieldVector3D.distance(computedMeanState.getPosition(),
1055 meanState.getPosition()).getReal(),
1056 5.0e-6);
1057
1058 }
1059
1060 @Test
1061 public void testIssue339() {
1062 doTestIssue339(Binary64Field.getInstance());
1063 }
1064
1065 private <T extends CalculusFieldElement<T>> void doTestIssue339(Field<T> field) {
1066
1067 final FieldSpacecraftState<T> osculatingState = getLEOState(field);
1068
1069 final CelestialBody sun = CelestialBodyFactory.getSun();
1070 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1071 Constants.WGS84_EARTH_FLATTENING,
1072 FramesFactory.getITRF(IERSConventions.IERS_2010,
1073 true));
1074 final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
1075 sun,
1076 50.0, Vector3D.PLUS_J,
1077 2.0, 0.1,
1078 0.2, 0.6);
1079 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
1080 final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(),
1081 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
1082 0.0, 0.0, 0.0);
1083
1084
1085 final Collection<DSSTForceModel> forces = new ArrayList<>();
1086 forces.add(new DSSTSolarRadiationPressure(sun, earth, boxAndWing, osculatingState.getOrbit().getMu().getReal()));
1087 forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getOrbit().getMu().getReal()));
1088
1089 final FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1090 Assertions.assertEquals(0.522,
1091 FieldVector3D.distance(osculatingState.getPosition(),
1092 meanState.getPosition()).getReal(),
1093 0.001);
1094
1095 final FieldSpacecraftState<T> computedOsculatingState = FieldDSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1096 Assertions.assertEquals(0.0,
1097 FieldVector3D.distance(osculatingState.getPosition(),
1098 computedOsculatingState.getPosition()).getReal(),
1099 5.0e-6);
1100
1101 }
1102
1103 @Test
1104 public void testIssue339WithAccelerations() {
1105 doTestIssue339WithAccelerations(Binary64Field.getInstance());
1106 }
1107
1108 private <T extends CalculusFieldElement<T>> void doTestIssue339WithAccelerations(Field<T> field) {
1109 final FieldSpacecraftState<T> osculatingState = getLEOStatePropagatedBy30Minutes(field);
1110 final CelestialBody sun = CelestialBodyFactory.getSun();
1111 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
1112 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1113 final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0, sun, 50.0, Vector3D.PLUS_J, 2.0, 0.1, 0.2, 0.6);
1114 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
1115 final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(), LOFType.LVLH_CCSDS, RotationOrder.XYZ, 0.0, 0.0, 0.0);
1116
1117 final Collection<DSSTForceModel> forces = new ArrayList<>();
1118 forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getOrbit().getMu().getReal()));
1119 final FieldSpacecraftState<T> meanState = FieldDSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1120 final FieldSpacecraftState<T> computedOsculatingState = FieldDSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1121 Assertions.assertEquals(0.0,
1122 FieldVector3D.distance(osculatingState.getPosition(), computedOsculatingState.getPosition()).getReal(),
1123 5.0e-6);
1124 }
1125
1126 @Test
1127 public void testIssue613() {
1128 doTestIssue613(Binary64Field.getInstance());
1129 }
1130
1131 private <T extends CalculusFieldElement<T>> void doTestIssue613(final Field<T> field) {
1132 final T zero = field.getZero();
1133
1134 final FieldSpacecraftState<T> state = getLEOState(field);
1135
1136
1137 final Frame itrf = FramesFactory .getITRF(IERSConventions.IERS_2010, true);
1138
1139
1140 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 4);
1141 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1142 Constants.WGS84_EARTH_FLATTENING, itrf);
1143
1144 final List<FieldEventDetector<T>> events = new ArrayList<>();
1145 events.add(new FieldAltitudeDetector<>(zero.add(85.5), earth));
1146 events.add(new FieldLatitudeCrossingDetector<>(field, earth, 0.0));
1147
1148
1149 final List<DSSTForceModel> forceModels = new ArrayList<>();
1150 forceModels.add(new DSSTZonal(provider));
1151 forceModels.add(new DSSTTesseral(itrf, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
1152 forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu()));
1153 forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu()));
1154
1155
1156 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(10).getTolerances(state.getOrbit(), OrbitType.EQUINOCTIAL);
1157 final FieldODEIntegrator<T> integrator = new DormandPrince54FieldIntegrator<>(field, 60.0, 3600.0, tol[0], tol[1]);
1158 final FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.OSCULATING);
1159 for (DSSTForceModel force : forceModels) {
1160 propagator.addForceModel(force);
1161 }
1162 for (FieldEventDetector<T> event : events) {
1163 propagator.addEventDetector(event);
1164 }
1165 propagator.setInitialState(state);
1166
1167
1168 final FieldSpacecraftState<T> finalState = propagator.propagate(state.getDate().shiftedBy(86400.0));
1169
1170
1171 Assertions.assertEquals(finalState.getOrbit().getMu().getReal(), 3.986004415E14, Double.MIN_VALUE);
1172 }
1173
1174 @Test
1175 public void testIssue704() {
1176 doTestIssue704(Binary64Field.getInstance());
1177 }
1178
1179 private <T extends CalculusFieldElement<T>> void doTestIssue704(final Field<T> field) {
1180
1181
1182 final FieldOrbit<T> orbit = getLEOState(field).getOrbit();
1183 final FieldPVCoordinates<T> pv = orbit.getPVCoordinates();
1184
1185
1186 final T dP = field.getZero().add(10.0);
1187
1188
1189 final T r2 = pv.getPosition().getNormSq();
1190 final T v = pv.getVelocity().getNorm();
1191 final T dV = dP.multiply(orbit.getMu()).divide(v.multiply(r2));
1192
1193
1194 final double[][] tol1 = FieldDSSTPropagator.tolerances(dP, orbit);
1195 final double[][] tol2 = FieldDSSTPropagator.tolerances(dP, dV, orbit);
1196 for (int i = 0; i < tol1.length; i++) {
1197 Assertions.assertArrayEquals(tol1[i], tol2[i], Double.MIN_VALUE);
1198 }
1199
1200 }
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239 @Test
1240 public void testIssue672() {
1241 doTestIssue672(Binary64Field.getInstance());
1242 }
1243
1244 private <T extends CalculusFieldElement<T>> void doTestIssue672(final Field<T> field) {
1245
1246
1247
1248
1249
1250 final UnnormalizedSphericalHarmonicsProvider provider =
1251 GravityFieldFactory.getUnnormalizedProvider(3, 3);
1252 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
1253
1254
1255 FieldSpacecraftState<T> initialState = getLEOState(field);
1256
1257
1258 final T zero = field.getZero();
1259 initialState.getDate();
1260 final T minStep = initialState.getOrbit().getKeplerianPeriod();
1261 final T maxStep = minStep.multiply(100.);
1262 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(initialState.getOrbit(), OrbitType.EQUINOCTIAL);
1263 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep.getReal(), maxStep.getReal(), tol[0], tol[1]);
1264 final FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.OSCULATING);
1265 propagator.setInitialState(initialState, PropagationType.OSCULATING);
1266 propagator.addForceModel(new DSSTZonal(provider));
1267 propagator.addForceModel(new DSSTTesseral(earthFrame,
1268 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
1269
1270
1271
1272
1273
1274 final FieldSpacecraftState<T> state = propagator.propagate(initialState.getDate().shiftedBy(3600.));
1275
1276
1277
1278
1279
1280 Assertions.assertNotNull(state);
1281 }
1282
1283
1284
1285
1286
1287 @Test
1288 public void testIssue1461() {
1289 doTestIssue1461(Binary64Field.getInstance());
1290 }
1291
1292
1293 private <T extends CalculusFieldElement<T>> void doTestIssue1461(final Field<T> field) {
1294
1295
1296
1297
1298 final T zero = field.getZero();
1299
1300 final double step = 60.;
1301 final int nStep = 10;
1302 final Frame gcrf = FramesFactory.getGCRF();
1303 final AttitudeProvider attitude = new FrameAlignedProvider(gcrf);
1304
1305
1306 final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(field,
1307 new KeplerianOrbit(7000e3, 1.e-3, FastMath.toRadians(98.), 0., 0., 0.,
1308 PositionAngleType.ECCENTRIC, gcrf,
1309 AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU));
1310 final FieldAbsoluteDate<T> t0 = initialOrbit.getDate();
1311 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(initialOrbit);
1312 final ClassicalRungeKuttaFieldIntegrator<T> integrator = new ClassicalRungeKuttaFieldIntegrator<>(field, zero.newInstance(step));
1313
1314
1315 final FieldDSSTPropagator<T> propagator = new FieldDSSTPropagator<>(field, integrator, PropagationType.MEAN);
1316
1317
1318 propagator.setInitialState(initialState, PropagationType.OSCULATING);
1319
1320
1321 propagator.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0)));
1322 propagator.setAttitudeProvider(attitude);
1323
1324
1325 final double initialMeanSma = FieldDSSTPropagator.computeMeanState(initialState,
1326 attitude,
1327 propagator.getAllForceModels()).getOrbit().getA().getReal();
1328
1329
1330
1331
1332 final List<Double> propagatedMeanSma = new ArrayList<>();
1333 for (int i = 1; i <= nStep; i++) {
1334 propagatedMeanSma.add(propagator.propagate(t0.shiftedBy(i * step)).getOrbit().getA().getReal());
1335 }
1336
1337
1338
1339
1340
1341 for (Double meanSma : propagatedMeanSma) {
1342 Assertions.assertEquals(initialMeanSma, meanSma, 0.);
1343 }
1344 }
1345
1346 private <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> getGEOState(final Field<T> field) {
1347 final T zero = field.getZero();
1348
1349 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
1350 TimeScalesFactory.getUTC());
1351 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(zero.add(42164000),
1352 zero.add(10e-3),
1353 zero.add(10e-3),
1354 zero.add(FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3)),
1355 zero.add(FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3)),
1356 zero.add(0.1),
1357 PositionAngleType.TRUE,
1358 FramesFactory.getEME2000(),
1359 initDate,
1360 zero.add(3.986004415E14));
1361 return new FieldSpacecraftState<>(orbit);
1362 }
1363
1364 private <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> getLEOState(final Field<T> field) {
1365 final T zero = field.getZero();
1366 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
1367 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
1368
1369 final FieldAbsoluteDate<T> initDate = new FieldAbsoluteDate<>(field, new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1370 return new FieldSpacecraftState<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
1371 FramesFactory.getEME2000(),
1372 initDate,
1373 zero.add(3.986004415E14)));
1374 }
1375
1376 private <T extends CalculusFieldElement<T>> FieldDSSTPropagator<T> setDSSTProp(Field<T> field,
1377 FieldSpacecraftState<T> initialState) {
1378 final T zero = field.getZero();
1379 initialState.getDate();
1380 final T minStep = initialState.getOrbit().getKeplerianPeriod();
1381 final T maxStep = minStep.multiply(100.);
1382 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(initialState.getOrbit(), OrbitType.EQUINOCTIAL);
1383 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep.getReal(), maxStep.getReal(), tol[0], tol[1]);
1384 final FieldDSSTPropagator<T> dsstProp = new FieldDSSTPropagator<>(field, integrator);
1385 dsstProp.setInitialState(initialState, PropagationType.MEAN);
1386
1387 return dsstProp;
1388 }
1389
1390 private class CheckingHandler<D extends FieldEventDetector<T>, T extends CalculusFieldElement<T>> implements FieldEventHandler<T> {
1391
1392 private final Action actionOnEvent;
1393 private boolean gotHere;
1394
1395 public CheckingHandler(final Action actionOnEvent) {
1396 this.actionOnEvent = actionOnEvent;
1397 this.gotHere = false;
1398 }
1399
1400 public void assertEvent(boolean expected) {
1401 Assertions.assertEquals(expected, gotHere);
1402 }
1403
1404 @Override
1405 public Action eventOccurred(FieldSpacecraftState<T> s, FieldEventDetector<T> detector, boolean increasing) {
1406 gotHere = true;
1407 return actionOnEvent;
1408 }
1409
1410 }
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420 @Override
1421 protected List<ParameterDriver> getParametersDriversWithoutMu() {
1422 return Collections.emptyList();
1423 }
1424
1425
1426 @Override
1427 protected double[] getLLimits(SpacecraftState state,
1428 AuxiliaryElements auxiliaryElements) {
1429 return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0),
1430 FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0) };
1431 }
1432
1433
1434 @Override
1435 protected <T extends CalculusFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
1436 FieldAuxiliaryElements<T> auxiliaryElements) {
1437 final Field<T> field = state.getDate().getField();
1438 final T zero = field.getZero();
1439 final T[] tab = MathArrays.buildArray(field, 2);
1440 tab[0] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).subtract(FastMath.PI);
1441 tab[1] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).add(FastMath.PI);
1442 return tab;
1443 }
1444
1445 }
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458 @Override
1459 public void init(final SpacecraftState s0, final AbsoluteDate t) {
1460 this.initialized = true;
1461 this.accComputed = false;
1462 }
1463
1464
1465 @Override
1466 public boolean dependsOnPositionOnly() {
1467 return false;
1468 }
1469
1470
1471 @Override
1472 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1473 return Vector3D.ZERO;
1474 }
1475
1476
1477 @Override
1478 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1479 this.accComputed = true;
1480 return FieldVector3D.getZero(s.getDate().getField());
1481 }
1482
1483
1484 @Override
1485 public Stream<EventDetector> getEventDetectors() {
1486 return Stream.empty();
1487 }
1488
1489
1490 @Override
1491 public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
1492 return Stream.empty();
1493 }
1494
1495
1496 @Override
1497 public List<ParameterDriver> getParametersDrivers() {
1498 return Collections.emptyList();
1499 }
1500
1501 }
1502
1503 @BeforeEach
1504 public void setUp() throws IOException, ParseException {
1505 Utils.setDataRoot("regular-data:potential/shm-format");
1506 }
1507
1508 private <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> getLEOStatePropagatedBy30Minutes(Field<T> field) {
1509
1510 final T zero = field.getZero();
1511 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
1512 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
1513
1514 final FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1515 final FieldCartesianOrbit<T> osculatingOrbit = new FieldCartesianOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getTOD(IERSConventions.IERS_1996, false),
1516 initialDate, zero.add(Constants.WGS84_EARTH_MU));
1517
1518
1519 double minStep = 0.001;
1520 double maxstep = 1000.0;
1521 T positionTolerance = zero.add(10.0);
1522 OrbitType propagationType = OrbitType.EQUINOCTIAL;
1523 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(positionTolerance.getReal()).getTolerances(osculatingOrbit, propagationType);
1524 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxstep, tolerances[0], tolerances[1]);
1525 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
1526 propagator.setOrbitType(propagationType);
1527
1528 NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
1529 ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
1530 propagator.addForceModel(holmesFeatherstone);
1531 propagator.setInitialState(new FieldSpacecraftState<>(osculatingOrbit));
1532
1533 return propagator.propagate(new FieldAbsoluteDate<>(initialDate, 1800.));
1534 }
1535
1536 }