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