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