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