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 static org.hamcrest.Matchers.contains;
20 import static org.hamcrest.Matchers.is;
21 import static org.hamcrest.Matchers.nullValue;
22
23 import java.io.IOException;
24 import java.text.ParseException;
25 import java.util.ArrayList;
26 import java.util.Collection;
27 import java.util.Collections;
28 import java.util.HashSet;
29 import java.util.List;
30 import java.util.Set;
31 import java.util.stream.Stream;
32
33 import org.hamcrest.Matcher;
34 import org.hamcrest.MatcherAssert;
35 import org.hamcrest.Matchers;
36 import org.hipparchus.CalculusFieldElement;
37 import org.hipparchus.Field;
38 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
39 import org.hipparchus.geometry.euclidean.threed.RotationOrder;
40 import org.hipparchus.geometry.euclidean.threed.Vector3D;
41 import org.hipparchus.linear.MatrixUtils;
42 import org.hipparchus.linear.RealMatrix;
43 import org.hipparchus.ode.ODEIntegrator;
44 import org.hipparchus.ode.events.Action;
45 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
46 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
47 import org.hipparchus.ode.nonstiff.DormandPrince54Integrator;
48 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
49 import org.hipparchus.util.FastMath;
50 import org.hipparchus.util.MathArrays;
51 import org.hipparchus.util.MathUtils;
52 import org.hipparchus.util.Precision;
53 import org.junit.jupiter.api.AfterEach;
54 import org.junit.jupiter.api.Assertions;
55 import org.junit.jupiter.api.BeforeEach;
56 import org.junit.jupiter.api.Test;
57 import org.orekit.OrekitMatchers;
58 import org.orekit.Utils;
59 import org.orekit.attitudes.AttitudeProvider;
60 import org.orekit.attitudes.AttitudeProviderModifier;
61 import org.orekit.attitudes.FrameAlignedProvider;
62 import org.orekit.attitudes.LofOffset;
63 import org.orekit.bodies.CelestialBody;
64 import org.orekit.bodies.CelestialBodyFactory;
65 import org.orekit.bodies.OneAxisEllipsoid;
66 import org.orekit.errors.OrekitException;
67 import org.orekit.forces.BoxAndSolarArraySpacecraft;
68 import org.orekit.forces.ForceModel;
69 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
70 import org.orekit.forces.gravity.potential.GRGSFormatReader;
71 import org.orekit.forces.gravity.potential.GravityFieldFactory;
72 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
73 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
74 import org.orekit.forces.gravity.potential.TideSystem;
75 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
76 import org.orekit.forces.maneuvers.ImpulseManeuver;
77 import org.orekit.frames.Frame;
78 import org.orekit.frames.FramesFactory;
79 import org.orekit.frames.LOFType;
80 import org.orekit.models.earth.atmosphere.Atmosphere;
81 import org.orekit.models.earth.atmosphere.HarrisPriester;
82 import org.orekit.orbits.CartesianOrbit;
83 import org.orekit.orbits.CircularOrbit;
84 import org.orekit.orbits.EquinoctialOrbit;
85 import org.orekit.orbits.KeplerianOrbit;
86 import org.orekit.orbits.Orbit;
87 import org.orekit.orbits.OrbitType;
88 import org.orekit.orbits.PositionAngleType;
89 import org.orekit.propagation.*;
90 import org.orekit.propagation.events.AltitudeDetector;
91 import org.orekit.propagation.events.ApsideDetector;
92 import org.orekit.propagation.events.DateDetector;
93 import org.orekit.propagation.events.EventDetector;
94 import org.orekit.propagation.events.FieldEventDetector;
95 import org.orekit.propagation.events.LatitudeCrossingDetector;
96 import org.orekit.propagation.events.NodeDetector;
97 import org.orekit.propagation.events.handlers.CountAndContinue;
98 import org.orekit.propagation.events.handlers.EventHandler;
99 import org.orekit.propagation.numerical.NumericalPropagator;
100 import org.orekit.propagation.semianalytical.dsst.forces.AbstractGaussianContribution;
101 import org.orekit.propagation.semianalytical.dsst.forces.DSSTAtmosphericDrag;
102 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
103 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
104 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
105 import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
106 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
107 import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
108 import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
109 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
110 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
111 import org.orekit.time.AbsoluteDate;
112 import org.orekit.time.DateComponents;
113 import org.orekit.time.TimeComponents;
114 import org.orekit.time.TimeScale;
115 import org.orekit.time.TimeScalesFactory;
116 import org.orekit.utils.Constants;
117 import org.orekit.utils.IERSConventions;
118 import org.orekit.utils.PVCoordinates;
119 import org.orekit.utils.ParameterDriver;
120 import org.orekit.utils.TimeSpanMap;
121 import org.orekit.utils.TimeStampedPVCoordinates;
122
123 public class DSSTPropagatorTest {
124
125 private DSSTPropagator dsstProp;
126
127 @Test
128 public void testIssue1703() {
129
130 final AbsoluteDate orbitEpoch = new AbsoluteDate(2023, 2, 18, TimeScalesFactory.getUTC());
131 final Frame inertial = FramesFactory.getCIRF(IERSConventions.IERS_2010, true);
132 final KeplerianOrbit orbit = new KeplerianOrbit(42166000.0, 0.00028, FastMath.toRadians(0.05), FastMath.toRadians(66.0),
133 FastMath.toRadians(270.0), FastMath.toRadians(11.94), PositionAngleType.MEAN,
134 inertial, orbitEpoch, Constants.WGS84_EARTH_MU);
135
136
137 DSSTPropagator propagator = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(60.0), PropagationType.OSCULATING);
138 propagator.setInitialState(new SpacecraftState(orbit), PropagationType.MEAN);
139
140
141 Assertions.assertFalse(propagator.initialIsOsculating());
142
143
144 propagator.setupMatricesComputation("stm", null, null);
145 propagator.setUpStmAndJacobianGenerators();
146
147
148 Assertions.assertFalse(propagator.initialIsOsculating());
149 }
150
151
152
153
154
155
156 @Test
157 public void testIssue1029() {
158
159 Utils.setDataRoot("regular-data:potential/grgs-format");
160 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
161
162
163 final AbsoluteDate orbitEpoch = new AbsoluteDate(2023, 2, 18, TimeScalesFactory.getUTC());
164 final Frame inertial = FramesFactory.getCIRF(IERSConventions.IERS_2010, true);
165 final KeplerianOrbit orbit = new KeplerianOrbit(42166000.0, 0.00028, FastMath.toRadians(0.05), FastMath.toRadians(66.0),
166 FastMath.toRadians(270.0), FastMath.toRadians(11.94), PositionAngleType.MEAN,
167 inertial, orbitEpoch, Constants.WGS84_EARTH_MU);
168 final EquinoctialOrbit equinoctial = new EquinoctialOrbit(orbit);
169
170
171 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.e-3).getTolerances(equinoctial, OrbitType.EQUINOCTIAL);
172 final AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(3600.0, 86400.0, tol[0], tol[1]);
173 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
174
175
176 final Frame ecefFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
177 final UnnormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getUnnormalizedProvider(8, 8);
178 propagator.addForceModel(new DSSTZonal(gravityProvider));
179 propagator.addForceModel(new DSSTTesseral(ecefFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, gravityProvider));
180 propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun(), gravityProvider.getMu()));
181 propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon(), gravityProvider.getMu()));
182
183
184 propagator.setInitialState(new SpacecraftState(equinoctial).withMass(6000.0), PropagationType.MEAN);
185 SpacecraftState propagated = propagator.propagate(orbitEpoch.shiftedBy(20.0 * Constants.JULIAN_DAY));
186
187
188 Assertions.assertEquals(0.0, propagated.getDate().durationFrom(orbitEpoch.shiftedBy(20.0 * Constants.JULIAN_DAY)), Double.MIN_VALUE);
189 MatcherAssert.assertThat( propagated.getOrbit().getA(),
190 Matchers.closeTo(4.216464862956647E7, 5e-7));
191
192 }
193
194 @Test
195 public void testIssue363() {
196 Utils.setDataRoot("regular-data");
197 AbsoluteDate date = new AbsoluteDate("2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
198 CircularOrbit orbit = new CircularOrbit(7389068.5, 1.0e-15, 1.0e-15, 1.709573, 1.308398, 0, PositionAngleType.MEAN,
199 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
200 date, Constants.WGS84_EARTH_MU);
201 SpacecraftState osculatingState = new SpacecraftState(orbit).withMass(1116.2829);
202
203 List<DSSTForceModel> dsstForceModels = new ArrayList<DSSTForceModel>();
204
205 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu()));
206 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu()));
207
208 SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
209 Assertions.assertEquals( 0.421, osculatingState.getOrbit().getA() - meanState.getOrbit().getA(), 1.0e-3);
210 Assertions.assertEquals(-5.23e-8, osculatingState.getOrbit().getEquinoctialEx() - meanState.getOrbit().getEquinoctialEx(), 1.0e-10);
211 Assertions.assertEquals(15.22e-8, osculatingState.getOrbit().getEquinoctialEy() - meanState.getOrbit().getEquinoctialEy(), 1.0e-10);
212 Assertions.assertEquals(-3.15e-8, osculatingState.getOrbit().getHx() - meanState.getOrbit().getHx(), 1.0e-10);
213 Assertions.assertEquals( 2.83e-8, osculatingState.getOrbit().getHy() - meanState.getOrbit().getHy(), 1.0e-10);
214 Assertions.assertEquals(15.96e-8, osculatingState.getOrbit().getLM() - meanState.getOrbit().getLM(), 1.0e-10);
215
216 }
217
218 @Test
219 public void testIssue364() {
220 Utils.setDataRoot("regular-data");
221 AbsoluteDate date = new AbsoluteDate("2003-06-18T00:00:00.000", TimeScalesFactory.getUTC());
222 CircularOrbit orbit = new CircularOrbit(7389068.5, 0.0, 0.0, 1.709573, 1.308398, 0, PositionAngleType.MEAN,
223 FramesFactory.getTOD(IERSConventions.IERS_2010, false),
224 date, Constants.WGS84_EARTH_MU);
225 SpacecraftState osculatingState = new SpacecraftState(orbit).withMass(1116.2829);
226
227 List<DSSTForceModel> dsstForceModels = new ArrayList<DSSTForceModel>();
228
229 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), orbit.getMu()));
230 dsstForceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), orbit.getMu()));
231
232 SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, null, dsstForceModels);
233 Assertions.assertEquals( 0.421, osculatingState.getOrbit().getA() - meanState.getOrbit().getA(), 1.0e-3);
234 Assertions.assertEquals(-5.23e-8, osculatingState.getOrbit().getEquinoctialEx() - meanState.getOrbit().getEquinoctialEx(), 1.0e-10);
235 Assertions.assertEquals(15.22e-8, osculatingState.getOrbit().getEquinoctialEy() - meanState.getOrbit().getEquinoctialEy(), 1.0e-10);
236 Assertions.assertEquals(-3.15e-8, osculatingState.getOrbit().getHx() - meanState.getOrbit().getHx(), 1.0e-10);
237 Assertions.assertEquals( 2.83e-8, osculatingState.getOrbit().getHy() - meanState.getOrbit().getHy(), 1.0e-10);
238 Assertions.assertEquals(15.96e-8, osculatingState.getOrbit().getLM() - meanState.getOrbit().getLM(), 1.0e-10);
239
240 }
241
242 @Test
243 public void testHighDegreesSetting() {
244
245 Utils.setDataRoot("regular-data:potential/grgs-format");
246 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
247 int earthDegree = 36;
248 int earthOrder = 36;
249 int eccPower = 4;
250 final UnnormalizedSphericalHarmonicsProvider provider =
251 GravityFieldFactory.getUnnormalizedProvider(earthDegree, earthOrder);
252 final org.orekit.frames.Frame earthFrame =
253 FramesFactory.getITRF(IERSConventions.IERS_2010, true);
254 final DSSTForceModel force =
255 new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
256 earthDegree, earthOrder, eccPower, earthDegree + eccPower,
257 earthDegree, earthOrder, eccPower);
258 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
259 forces.add(force);
260 TimeScale tai = TimeScalesFactory.getTAI();
261 AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
262 Frame eci = FramesFactory.getGCRF();
263 KeplerianOrbit orbit = new KeplerianOrbit(
264 7120000.0, 1.0e-3, FastMath.toRadians(60.0),
265 FastMath.toRadians(120.0), FastMath.toRadians(47.0),
266 FastMath.toRadians(12.0),
267 PositionAngleType.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
268 SpacecraftState oscuState = DSSTPropagator.computeOsculatingState(new SpacecraftState(orbit), null, forces);
269 Assertions.assertEquals(7119927.147704, oscuState.getOrbit().getA(), 1.e-6);
270 }
271
272 @Test
273 public void testEphemerisDates() {
274
275 TimeScale tai = TimeScalesFactory.getTAI();
276 AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
277 AbsoluteDate startDate = new AbsoluteDate("2015-07-03", tai).shiftedBy(-0.1);
278 AbsoluteDate endDate = new AbsoluteDate("2015-07-04", tai);
279 Frame eci = FramesFactory.getGCRF();
280 KeplerianOrbit orbit = new KeplerianOrbit(
281 600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0, 0, 0, 0,
282 PositionAngleType.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
283 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1).getTolerances(orbit, OrbitType.EQUINOCTIAL);
284 Propagator prop = new DSSTPropagator(
285 new DormandPrince853Integrator(0.1, 500, tol[0], tol[1]));
286 prop.resetInitialState(new SpacecraftState(new CartesianOrbit(orbit)));
287
288
289 EphemerisGenerator generator = prop.getEphemerisGenerator();
290 prop.propagate(startDate, endDate);
291 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
292
293
294 TimeStampedPVCoordinates actualPV = ephemeris.getPVCoordinates(startDate, eci);
295 TimeStampedPVCoordinates expectedPV = orbit.getPVCoordinates(startDate, eci);
296 MatcherAssert.assertThat(actualPV.getPosition(),
297 OrekitMatchers.vectorCloseTo(expectedPV.getPosition(), 1.0));
298 MatcherAssert.assertThat(actualPV.getVelocity(),
299 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity(), 1.0));
300 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate),
301 OrekitMatchers.closeTo(0, 0));
302 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate),
303 OrekitMatchers.closeTo(0, 0));
304
305 AbsoluteDate date = endDate.shiftedBy(-0.11);
306 Assertions.assertEquals(
307 ephemeris.propagate(date).getDate().durationFrom(date), 0, 0);
308 }
309
310 @Test
311 public void testNoExtrapolation() {
312 SpacecraftState state = getLEOState();
313 setDSSTProp(state);
314
315
316 final SpacecraftState finalState = dsstProp.propagate(state.getDate());
317
318
319 final Vector3D initialPosition = state.getPosition();
320 final Vector3D initialVelocity = state.getVelocity();
321
322
323 final Vector3D finalPosition = finalState.getPosition();
324 final Vector3D finalVelocity = finalState.getVelocity();
325
326
327 Assertions.assertEquals(initialPosition.getX(), finalPosition.getX(), 0.0);
328 Assertions.assertEquals(initialPosition.getY(), finalPosition.getY(), 0.0);
329 Assertions.assertEquals(initialPosition.getZ(), finalPosition.getZ(), 0.0);
330 Assertions.assertEquals(initialVelocity.getX(), finalVelocity.getX(), 0.0);
331 Assertions.assertEquals(initialVelocity.getY(), finalVelocity.getY(), 0.0);
332 Assertions.assertEquals(initialVelocity.getZ(), finalVelocity.getZ(), 0.0);
333 }
334
335 @Test
336 public void testKepler() {
337 SpacecraftState state = getLEOState();
338 setDSSTProp(state);
339
340 Assertions.assertEquals(2, dsstProp.getSatelliteRevolution());
341 dsstProp.setSatelliteRevolution(17);
342 Assertions.assertEquals(17, dsstProp.getSatelliteRevolution());
343
344
345 final double dt = 3200.;
346 final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
347
348
349 final double n = FastMath.sqrt(state.getOrbit().getMu() / state.getOrbit().getA()) / state.getOrbit().getA();
350 Assertions.assertEquals(state.getOrbit().getA(), finalState.getOrbit().getA(), 0.);
351 Assertions.assertEquals(state.getOrbit().getEquinoctialEx(), finalState.getOrbit().getEquinoctialEx(), 0.);
352 Assertions.assertEquals(state.getOrbit().getEquinoctialEy(), finalState.getOrbit().getEquinoctialEy(), 0.);
353 Assertions.assertEquals(state.getOrbit().getHx(), finalState.getOrbit().getHx(), 0.);
354 Assertions.assertEquals(state.getOrbit().getHy(), finalState.getOrbit().getHy(), 0.);
355 Assertions.assertEquals(state.getOrbit().getLM() + n * dt, finalState.getOrbit().getLM(), 1.e-14);
356
357 }
358
359 @Test
360 public void testEphemeris() {
361 SpacecraftState state = getGEOState();
362 setDSSTProp(state);
363
364
365 EphemerisGenerator generator = dsstProp.getEphemerisGenerator();
366
367
368 final double dt = 2. * Constants.JULIAN_DAY;
369 dsstProp.propagate(state.getDate().shiftedBy(5. * dt));
370
371
372 BoundedPropagator ephem = generator.getGeneratedEphemeris();
373
374
375 final SpacecraftState s = ephem.propagate(state.getDate().shiftedBy(dt));
376
377
378 final double n = FastMath.sqrt(state.getOrbit().getMu() / state.getOrbit().getA()) / state.getOrbit().getA();
379 Assertions.assertEquals(state.getOrbit().getA(), s.getOrbit().getA(), 0.);
380 Assertions.assertEquals(state.getOrbit().getEquinoctialEx(), s.getOrbit().getEquinoctialEx(), 0.);
381 Assertions.assertEquals(state.getOrbit().getEquinoctialEy(), s.getOrbit().getEquinoctialEy(), 0.);
382 Assertions.assertEquals(state.getOrbit().getHx(), s.getOrbit().getHx(), 0.);
383 Assertions.assertEquals(state.getOrbit().getHy(), s.getOrbit().getHy(), 0.);
384 Assertions.assertEquals(state.getOrbit().getLM() + n * dt, s.getOrbit().getLM(), 1.5e-14);
385
386 }
387
388 @Test
389 public void testImpulseManeuver() {
390 final Orbit initialOrbit = new KeplerianOrbit(24532000.0, 0.72, 0.3, FastMath.PI, 0.4, 2.0, PositionAngleType.MEAN, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2008, 06, 23), new TimeComponents(14, 18, 37), TimeScalesFactory.getUTC()), 3.986004415e14);
391 final double a = initialOrbit.getA();
392 final double e = initialOrbit.getE();
393 final double i = initialOrbit.getI();
394 final double mu = initialOrbit.getMu();
395 final double vApo = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
396 double dv = 0.99 * FastMath.tan(i) * vApo;
397
398
399 setDSSTProp(new SpacecraftState(initialOrbit));
400
401
402 dsstProp.setAttitudeProvider(new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS));
403 dsstProp.addEventDetector(new ImpulseManeuver(new NodeDetector(initialOrbit, FramesFactory.getEME2000()), new Vector3D(dv, Vector3D.PLUS_J), 400.0));
404 SpacecraftState propagated = dsstProp.propagate(initialOrbit.getDate().shiftedBy(8000));
405
406 Assertions.assertEquals(0.0028257, propagated.getOrbit().getI(), 1.0e-6);
407 }
408
409 @Test
410 public void testPropagationWithCentralBody() throws Exception {
411
412
413 final UnnormalizedSphericalHarmonicsProvider provider =
414 GravityFieldFactory.getUnnormalizedProvider(4, 4);
415 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
416
417
418 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400,
419 TimeScalesFactory.getUTC());
420 final Orbit orbit = new KeplerianOrbit(26559890.,
421 0.0041632,
422 FastMath.toRadians(55.2),
423 FastMath.toRadians(315.4985),
424 FastMath.toRadians(130.7562),
425 FastMath.toRadians(44.2377),
426 PositionAngleType.MEAN,
427 FramesFactory.getEME2000(),
428 initDate,
429 provider.getMu());
430
431
432 setDSSTProp(new SpacecraftState(orbit));
433 dsstProp.addForceModel(new DSSTZonal(provider, 4, 3, 9));
434 dsstProp.addForceModel(new DSSTTesseral(earthFrame,
435 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
436 4, 4, 4, 8, 4, 4, 2));
437
438
439 final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
440
441
442
443
444
445
446
447
448 Assertions.assertEquals(26559920.81, state.getOrbit().getA(), 1.e-1);
449 Assertions.assertEquals(0.2731622444E-03, state.getOrbit().getEquinoctialEx(), 2.e-8);
450 Assertions.assertEquals(0.4164167597E-02, state.getOrbit().getEquinoctialEy(), 2.e-8);
451 Assertions.assertEquals(-0.3399607878, state.getOrbit().getHx(), 5.e-8);
452 Assertions.assertEquals(0.3971568634, state.getOrbit().getHy(), 2.e-6);
453 Assertions.assertEquals(140.6375352,
454 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), FastMath.PI)),
455 5.e-3);
456 }
457
458 @Test
459 public void testPropagationWithThirdBody() throws IOException {
460
461
462 final UnnormalizedSphericalHarmonicsProvider provider =
463 GravityFieldFactory.getUnnormalizedProvider(2, 0);
464 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
465 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
466 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
467 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
468 provider, 2, 0, 0, 2, 2, 0, 0);
469
470
471 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
472 DSSTForceModel sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu());
473
474
475 final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
476 TimeScalesFactory.getUTC());
477 final Orbit orbit = new KeplerianOrbit(42163393.,
478 0.2684,
479 FastMath.toRadians(63.435),
480 FastMath.toRadians(270.0),
481 FastMath.toRadians(285.0),
482 FastMath.toRadians(344.0),
483 PositionAngleType.MEAN,
484 FramesFactory.getEME2000(),
485 initDate,
486 provider.getMu());
487
488
489 setDSSTProp(new SpacecraftState(orbit));
490 dsstProp.addForceModel(zonal);
491 dsstProp.addForceModel(tesseral);
492 dsstProp.addForceModel(moon);
493 dsstProp.addForceModel(sun);
494
495
496 final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
497
498
499
500
501
502
503
504
505 Assertions.assertEquals(42163393.0, state.getOrbit().getA(), 1.e-1);
506 Assertions.assertEquals(-0.2592789733084587, state.getOrbit().getEquinoctialEx(), 5.e-7);
507 Assertions.assertEquals(-0.06893353670734315, state.getOrbit().getEquinoctialEy(), 2.e-7);
508 Assertions.assertEquals( 0.1595005111738418, state.getOrbit().getHx(), 2.e-7);
509 Assertions.assertEquals(-0.5968524904937771, state.getOrbit().getHy(), 5.e-8);
510 Assertions.assertEquals(183.9386620425922,
511 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), FastMath.PI)),
512 3.e-2);
513 }
514
515 @Test
516 public void testTooSmallMaxDegree() {
517 Assertions.assertThrows(OrekitException.class, () -> {
518 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 1, 0, 3);
519 });
520 }
521
522 @Test
523 public void testTooLargeMaxDegree() {
524 Assertions.assertThrows(OrekitException.class, () -> {
525 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0), 8, 0, 8);
526 });
527 }
528
529 @Test
530 public void testWrongMaxPower() {
531 Assertions.assertThrows(OrekitException.class, () -> {
532 new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(8, 8), 4, 4, 4);
533 });
534 }
535
536 @Test
537 public void testPropagationWithDrag() {
538
539
540 final UnnormalizedSphericalHarmonicsProvider provider =
541 GravityFieldFactory.getUnnormalizedProvider(2, 0);
542 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
543 DSSTForceModel zonal = new DSSTZonal(provider, 2, 0, 5);
544 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
545 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
546 provider, 2, 0, 0, 2, 2, 0, 0);
547
548
549 final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
550 Constants.WGS84_EARTH_FLATTENING,
551 earthFrame);
552 final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
553
554 final double cd = 2.0;
555 final double area = 25.0;
556 DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area, provider.getMu());
557
558
559 final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000,
560 TimeScalesFactory.getUTC());
561 final Orbit orbit = new KeplerianOrbit(7204535.848109440,
562 0.0012402238462686,
563 FastMath.toRadians(98.74341600466740),
564 FastMath.toRadians(111.1990175076630),
565 FastMath.toRadians(43.32990110790340),
566 FastMath.toRadians(68.66852509725620),
567 PositionAngleType.MEAN,
568 FramesFactory.getEME2000(),
569 initDate,
570 provider.getMu());
571
572
573 setDSSTProp(new SpacecraftState(orbit));
574 dsstProp.addForceModel(zonal);
575 dsstProp.addForceModel(tesseral);
576 dsstProp.addForceModel(drag);
577
578
579 final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
580
581
582
583
584
585
586
587
588 Assertions.assertEquals(7204521.657141485, state.getOrbit().getA(), 6.e-1);
589 Assertions.assertEquals(-0.001016800430994036, state.getOrbit().getEquinoctialEx(), 5.e-8);
590 Assertions.assertEquals(0.0007093755541595772, state.getOrbit().getEquinoctialEy(), 2.e-8);
591 Assertions.assertEquals(0.7757573478894775, state.getOrbit().getHx(), 5.e-8);
592 Assertions.assertEquals(0.8698955648709271, state.getOrbit().getHy(), 5.e-8);
593 Assertions.assertEquals(193.0939742953394,
594 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), FastMath.PI)),
595 2.e-3);
596
597
598 Assertions.assertEquals(((DSSTAtmosphericDrag)drag).getAtmosphere(), atm);
599
600 final double atmosphericMaxConstant = 1000000.0;
601 Assertions.assertEquals(((DSSTAtmosphericDrag)drag).getRbar(), atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 1e-9);
602 }
603
604 @Test
605 public void testPropagationWithSolarRadiationPressure() {
606
607
608 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
609 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
610 DSSTForceModel tesseral = new DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(),
611 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
612 provider, 2, 0, 0, 2, 2, 0, 0);
613
614
615 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
616 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
617 Constants.WGS84_EARTH_FLATTENING,
618 CelestialBodyFactory.getEarth().getBodyOrientedFrame()),
619 provider.getMu());
620
621
622 final AbsoluteDate initDate = new AbsoluteDate(2003, 9, 16, 0, 0, 00.000,
623 TimeScalesFactory.getUTC());
624 final Orbit orbit = new KeplerianOrbit(42166258.,
625 0.0001,
626 FastMath.toRadians(0.001),
627 FastMath.toRadians(315.4985),
628 FastMath.toRadians(130.7562),
629 FastMath.toRadians(44.2377),
630 PositionAngleType.MEAN,
631 FramesFactory.getGCRF(),
632 initDate,
633 provider.getMu());
634
635
636 dsstProp = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(86400.));
637 dsstProp.setInitialState(new SpacecraftState(orbit), PropagationType.MEAN);
638 dsstProp.addForceModel(zonal);
639 dsstProp.addForceModel(tesseral);
640 dsstProp.addForceModel(srp);
641
642
643 final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(10. * 86400.));
644
645
646
647
648
649
650
651
652 Assertions.assertEquals(42166257.99807995, state.getOrbit().getA(), 0.9);
653 Assertions.assertEquals(-0.1781865038201885e-05, state.getOrbit().getEquinoctialEx(), 3.e-7);
654 Assertions.assertEquals(-0.1191876027555493e-03, state.getOrbit().getEquinoctialEy(), 4.e-6);
655 Assertions.assertEquals(-0.5624363171289686e-05, state.getOrbit().getHx(), 4.e-9);
656 Assertions.assertEquals( 0.6618387121369373e-05, state.getOrbit().getHy(), 3.e-10);
657 Assertions.assertEquals(140.3496229467104,
658 FastMath.toDegrees(MathUtils.normalizeAngle(state.getOrbit().getLM(), FastMath.PI)),
659 2.e-4);
660 }
661
662 @Test
663 public void testStopEvent() {
664 SpacecraftState state = getLEOState();
665 setDSSTProp(state);
666
667 final AbsoluteDate stopDate = state.getDate().shiftedBy(1000);
668 CheckingHandler checking = new CheckingHandler(Action.STOP);
669 dsstProp.addEventDetector(new DateDetector(stopDate).withHandler(checking));
670 checking.assertEvent(false);
671 final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(3200));
672 checking.assertEvent(true);
673 Assertions.assertEquals(0, finalState.getDate().durationFrom(stopDate), 1.0e-10);
674 }
675
676 @Test
677 public void testContinueEvent() {
678 SpacecraftState state = getLEOState();
679 setDSSTProp(state);
680
681 final AbsoluteDate resetDate = state.getDate().shiftedBy(1000);
682 CheckingHandler checking = new CheckingHandler(Action.CONTINUE);
683 dsstProp.addEventDetector(new DateDetector(resetDate).withHandler(checking));
684 final double dt = 3200;
685 checking.assertEvent(false);
686 final SpacecraftState finalState = dsstProp.propagate(state.getDate().shiftedBy(dt));
687 checking.assertEvent(true);
688 final double n = FastMath.sqrt(state.getOrbit().getMu() / state.getOrbit().getA()) / state.getOrbit().getA();
689 Assertions.assertEquals(state.getOrbit().getA(), finalState.getOrbit().getA(), 1.0e-10);
690 Assertions.assertEquals(state.getOrbit().getEquinoctialEx(), finalState.getOrbit().getEquinoctialEx(), 1.0e-10);
691 Assertions.assertEquals(state.getOrbit().getEquinoctialEy(), finalState.getOrbit().getEquinoctialEy(), 1.0e-10);
692 Assertions.assertEquals(state.getOrbit().getHx(), finalState.getOrbit().getHx(), 1.0e-10);
693 Assertions.assertEquals(state.getOrbit().getHy(), finalState.getOrbit().getHy(), 1.0e-10);
694 Assertions.assertEquals(state.getOrbit().getLM() + n * dt, finalState.getOrbit().getLM(), 6.0e-10);
695 }
696
697 @Test
698 public void testIssue157() {
699 Utils.setDataRoot("regular-data:potential/icgem-format");
700 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
701 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
702 Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngleType.MEAN,
703 FramesFactory.getTOD(false),
704 new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
705 nshp.getMu());
706 double period = orbit.getKeplerianPeriod();
707 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.EQUINOCTIAL);
708 AdaptiveStepsizeIntegrator integrator =
709 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
710 integrator.setInitialStepSize(10 * period);
711 DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.MEAN);
712 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
713 Constants.WGS84_EARTH_FLATTENING,
714 FramesFactory.getGTOD(false));
715 CelestialBody sun = CelestialBodyFactory.getSun();
716 CelestialBody moon = CelestialBodyFactory.getMoon();
717 propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
718 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
719 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
720 nshp, 8, 8, 4, 12, 8, 8, 4));
721 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
722 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
723 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
724 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
725
726
727 propagator.setInitialState(new SpacecraftState(orbit).withMass(45.0), PropagationType.OSCULATING);
728 SpacecraftState finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
729
730
731
732
733 Assertions.assertEquals(2187.2, orbit.getA() - finalState.getOrbit().getA(), 1.0);
734
735 propagator.setInitialState(new SpacecraftState(orbit).withMass(45.0), PropagationType.MEAN);
736 finalState = propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
737
738
739 Assertions.assertEquals(1475.9, orbit.getA() - finalState.getOrbit().getA(), 1.0);
740
741 }
742
743
744
745
746 @Test
747 public void testEphemerisGeneration() {
748
749
750
751
752 Utils.setDataRoot("regular-data:potential/icgem-format");
753 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
754 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
755 Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngleType.MEAN,
756 FramesFactory.getTOD(false),
757 new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
758 nshp.getMu());
759 double period = orbit.getKeplerianPeriod();
760 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.EQUINOCTIAL);
761 AdaptiveStepsizeIntegrator integrator =
762 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
763 integrator.setInitialStepSize(10 * period);
764 DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
765 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
766 Constants.WGS84_EARTH_FLATTENING,
767 FramesFactory.getGTOD(false));
768 CelestialBody sun = CelestialBodyFactory.getSun();
769 CelestialBody moon = CelestialBodyFactory.getMoon();
770 propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
771 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
772 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
773 nshp, 8, 8, 4, 12, 8, 8, 4));
774 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
775 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
776 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
777 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
778 propagator.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
779
780
781
782
783
784
785
786 final double nDays = 5.;
787
788
789 propagator.setInitialState(new SpacecraftState(orbit).withMass( 45.0), PropagationType.MEAN);
790 final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
791 propagator.setStepHandler(600, currentState -> states.add(currentState));
792 propagator.propagate(orbit.getDate().shiftedBy(nDays * Constants.JULIAN_DAY));
793
794
795 propagator.setInitialState(new SpacecraftState(orbit).withMass(45.0), PropagationType.MEAN);
796 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
797 propagator.propagate(orbit.getDate().shiftedBy(nDays * Constants.JULIAN_DAY));
798 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
799
800 double maxError = 0;
801 for (final SpacecraftState state : states) {
802 final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
803 final double error = Vector3D.distance(state.getPosition(),
804 fromEphemeris.getPosition());
805 maxError = FastMath.max(maxError, error);
806 }
807
808
809
810
811
812 Assertions.assertEquals(0.0, maxError, Precision.SAFE_MIN);
813 }
814
815 @Test
816 public void testGetInitialOsculatingState() throws IllegalArgumentException, OrekitException {
817 final SpacecraftState initialState = getGEOState();
818
819
820 final double minStep = initialState.getOrbit().getKeplerianPeriod() * 0.1;
821 final double maxStep = initialState.getOrbit().getKeplerianPeriod() * 10.0;
822 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.e-1).getTolerances(initialState.getOrbit(), OrbitType.CARTESIAN);
823 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
824
825
826 DSSTPropagator prop = new DSSTPropagator(integrator, PropagationType.MEAN);
827
828 final UnnormalizedSphericalHarmonicsProvider provider =
829 GravityFieldFactory.getUnnormalizedProvider(4, 0);
830 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
831 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
832 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
833 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
834 provider, 4, 0, 4, 8, 4, 0, 2);
835 prop.addForceModel(zonal);
836 prop.addForceModel(tesseral);
837
838
839 prop.setInitialState(initialState, PropagationType.MEAN);
840
841 Assertions.assertEquals(initialState, prop.getInitialState());
842
843
844 Assertions.assertEquals(initialState, prop.propagate(initialState.getDate()));
845 }
846
847 @Test
848 public void testMeanToOsculatingState() throws IllegalArgumentException, OrekitException {
849 final SpacecraftState meanState = getGEOState();
850
851 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
852 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
853 final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
854 final DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
855 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
856 provider, 2, 0, 0, 2, 2, 0, 0);
857
858 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
859 forces.add(zonal);
860 forces.add(tesseral);
861
862 final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
863 Assertions.assertEquals(1559.1,
864 Vector3D.distance(meanState.getPosition(),
865 osculatingState.getPosition()),
866 1.0);
867 }
868
869 @Test
870 public void testOsculatingToMeanState() throws IllegalArgumentException, OrekitException {
871 final SpacecraftState meanState = getGEOState();
872
873 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
874 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
875
876 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
877 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
878 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
879 provider, 2, 0, 0, 2, 2, 0, 0);
880
881 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
882 forces.add(zonal);
883 forces.add(tesseral);
884
885 final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
886
887
888 final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
889
890 Assertions.assertEquals(meanState.getOrbit().getA(), computedMeanState.getOrbit().getA(), 2.0e-8);
891 Assertions.assertEquals(0.0,
892 Vector3D.distance(meanState.getPosition(),
893 computedMeanState.getPosition()),
894 2.0e-8);
895 }
896
897 @Test
898 public void testShortPeriodCoefficients() {
899 Utils.setDataRoot("regular-data:potential/icgem-format");
900 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
901 UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(4, 4);
902 Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngleType.MEAN,
903 FramesFactory.getTOD(false),
904 new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
905 nshp.getMu());
906 double period = orbit.getKeplerianPeriod();
907 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.EQUINOCTIAL);
908 AdaptiveStepsizeIntegrator integrator =
909 new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
910 integrator.setInitialStepSize(10 * period);
911 DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
912 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
913 Constants.WGS84_EARTH_FLATTENING,
914 FramesFactory.getGTOD(false));
915 CelestialBody sun = CelestialBodyFactory.getSun();
916 CelestialBody moon = CelestialBodyFactory.getMoon();
917 propagator.addForceModel(new DSSTZonal(nshp, 4, 3, 9));
918 propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(),
919 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
920 nshp, 4, 4, 4, 8, 4, 4, 2));
921 propagator.addForceModel(new DSSTThirdBody(sun, nshp.getMu()));
922 propagator.addForceModel(new DSSTThirdBody(moon, nshp.getMu()));
923 propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180, nshp.getMu()));
924 propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth, nshp.getMu()));
925
926 final AbsoluteDate finalDate = orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY);
927 propagator.resetInitialState(new SpacecraftState(orbit).withMass( 45.0));
928 final SpacecraftState stateNoConfig = propagator.propagate(finalDate);
929 Assertions.assertEquals(0, stateNoConfig.getAdditionalDataValues().size());
930
931 Assertions.assertNull(propagator.getSelectedCoefficients());
932 propagator.setSelectedCoefficients(new HashSet<String>());
933 Assertions.assertNotNull(propagator.getSelectedCoefficients());
934 Assertions.assertTrue(propagator.getSelectedCoefficients().isEmpty());
935 propagator.resetInitialState(new SpacecraftState(orbit).withMass( 45.0));
936 final SpacecraftState stateConfigEmpty = propagator.propagate(finalDate);
937 Assertions.assertEquals(234, stateConfigEmpty.getAdditionalDataValues().size());
938
939 final Set<String> selected = new HashSet<String>();
940 selected.add("DSST-3rd-body-Moon-s[7]");
941 selected.add("DSST-central-body-tesseral-c[-2][3]");
942 propagator.setSelectedCoefficients(selected);
943 Assertions.assertEquals(2, propagator.getSelectedCoefficients().size());
944 propagator.resetInitialState(new SpacecraftState(orbit).withMass(45.0));
945 final SpacecraftState stateConfigeSelected = propagator.propagate(finalDate);
946 Assertions.assertEquals(selected.size(), stateConfigeSelected.getAdditionalDataValues().size());
947
948 propagator.setSelectedCoefficients(null);
949 propagator.resetInitialState(new SpacecraftState(orbit).withMass(45.0));
950 final SpacecraftState stateConfigNull = propagator.propagate(finalDate);
951 Assertions.assertEquals(0, stateConfigNull.getAdditionalDataValues().size());
952
953 }
954
955 @Test
956 public void testIssueMeanInclination() {
957
958 final double earthAe = 6378137.0;
959 final double earthMu = 3.9860044E14;
960 final double earthJ2 = 0.0010826;
961
962
963 Orbit orb = new KeplerianOrbit(new TimeStampedPVCoordinates(new AbsoluteDate("1992-10-08T15:20:38.821",
964 TimeScalesFactory.getUTC()),
965 new Vector3D(5392808.809823, -4187618.3357927715, -44206.638015847195),
966 new Vector3D(2337.4472786270794, 2474.0146611860464, 6778.507766114648)),
967 FramesFactory.getTOD(false), earthMu);
968 final SpacecraftState ss = new SpacecraftState(orb);
969 final UnnormalizedSphericalHarmonicsProvider provider =
970 GravityFieldFactory.getUnnormalizedProvider(earthAe, earthMu, TideSystem.UNKNOWN,
971 new double[][] { { 0.0 }, { 0.0 }, { -earthJ2 } },
972 new double[][] { { 0.0 }, { 0.0 }, { 0.0 } });
973 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
974 DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
975 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
976 Constants.WGS84_EARTH_ANGULAR_VELOCITY,
977 provider, 2, 0, 0, 2, 2, 0, 0);
978 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
979 forces.add(zonal);
980 forces.add(tesseral);
981
982 final Orbit meanOrb = DSSTPropagator.computeMeanState(ss, null, forces).getOrbit();
983 Assertions.assertEquals(0.0164196, FastMath.toDegrees(orb.getI() - meanOrb.getI()), 1.0e-7);
984 }
985
986 @Test
987 public void testIssue257() {
988 final SpacecraftState meanState = getGEOState();
989
990
991 final DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), meanState.getOrbit().getMu());
992 final DSSTForceModel sun = new DSSTThirdBody(CelestialBodyFactory.getSun(), meanState.getOrbit().getMu());
993
994 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
995 forces.add(moon);
996 forces.add(sun);
997
998 final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
999 Assertions.assertEquals(734.3,
1000 Vector3D.distance(meanState.getPosition(),
1001 osculatingState.getPosition()),
1002 1.0);
1003
1004 final SpacecraftState computedMeanState = DSSTPropagator.computeMeanState(osculatingState, null, forces);
1005 Assertions.assertEquals(734.3,
1006 Vector3D.distance(osculatingState.getPosition(),
1007 computedMeanState.getPosition()),
1008 1.0);
1009
1010 Assertions.assertEquals(0.0,
1011 Vector3D.distance(computedMeanState.getPosition(),
1012 meanState.getPosition()),
1013 5.0e-6);
1014
1015 }
1016
1017 @Test
1018 public void testIssue339() {
1019
1020 final SpacecraftState osculatingState = getLEOState();
1021
1022 final CelestialBody sun = CelestialBodyFactory.getSun();
1023 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1024 Constants.WGS84_EARTH_FLATTENING,
1025 FramesFactory.getITRF(IERSConventions.IERS_2010,
1026 true));
1027 final BoxAndSolarArraySpacecraft boxAndWing = new BoxAndSolarArraySpacecraft(5.0, 2.0, 2.0,
1028 sun,
1029 50.0, Vector3D.PLUS_J,
1030 2.0, 0.1,
1031 0.2, 0.6);
1032 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
1033 final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(),
1034 LOFType.LVLH_CCSDS, RotationOrder.XYZ,
1035 0.0, 0.0, 0.0);
1036
1037
1038 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
1039 forces.add(new DSSTSolarRadiationPressure(sun, earth, boxAndWing, osculatingState.getOrbit().getMu()));
1040 forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getOrbit().getMu()));
1041
1042 final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1043 Assertions.assertEquals(0.522,
1044 Vector3D.distance(osculatingState.getPosition(),
1045 meanState.getPosition()),
1046 0.001);
1047
1048 final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1049 Assertions.assertEquals(0.0,
1050 Vector3D.distance(osculatingState.getPosition(),
1051 computedOsculatingState.getPosition()),
1052 5.0e-6);
1053
1054 }
1055
1056
1057 @Test
1058 public void testIssue1907() {
1059
1060 final SpacecraftState state = getLEOState();
1061
1062
1063 final Frame itrf = FramesFactory .getITRF(IERSConventions.IERS_2010, true);
1064
1065
1066 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 4);
1067
1068
1069 final List<EventDetector> events = new ArrayList<>();
1070 events.add(new ApsideDetector(state.getOrbit()).withHandler(new ApsideHandlerWithResetState()));
1071
1072
1073 final List<DSSTForceModel> forceModels = new ArrayList<>();
1074 forceModels.add(new DSSTZonal(provider));
1075 forceModels.add(new DSSTTesseral(itrf, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
1076
1077
1078 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(state.getOrbit(), OrbitType.EQUINOCTIAL);
1079 final ODEIntegrator integrator = new DormandPrince54Integrator(60.0, 3600.0, tol[0], tol[1]);
1080 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1081 for (DSSTForceModel force : forceModels) {
1082 propagator.addForceModel(force);
1083 }
1084 for (EventDetector event : events) {
1085 propagator.addEventDetector(event);
1086 }
1087 propagator.setInitialState(state);
1088
1089
1090 Assertions.assertDoesNotThrow(() -> propagator.propagate(state.getDate().shiftedBy(3600)));
1091 }
1092
1093 public class ApsideHandlerWithResetState implements EventHandler {
1094 @Override
1095 public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
1096 return Action.RESET_STATE;
1097 }
1098 }
1099
1100 @Test
1101 public void testIssue613() {
1102
1103 final SpacecraftState state = getLEOState();
1104
1105
1106 final Frame itrf = FramesFactory .getITRF(IERSConventions.IERS_2010, true);
1107
1108
1109 final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 4);
1110 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1111 Constants.WGS84_EARTH_FLATTENING, itrf);
1112
1113 final List<EventDetector> events = new ArrayList<>();
1114 events.add(new AltitudeDetector(85.5, earth));
1115 events.add(new LatitudeCrossingDetector(earth, 0.0));
1116
1117
1118 final List<DSSTForceModel> forceModels = new ArrayList<>();
1119 forceModels.add(new DSSTZonal(provider));
1120 forceModels.add(new DSSTTesseral(itrf, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
1121 forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu()));
1122 forceModels.add(new DSSTThirdBody(CelestialBodyFactory.getSun(), provider.getMu()));
1123
1124
1125 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(state.getOrbit(), OrbitType.EQUINOCTIAL);
1126 final ODEIntegrator integrator = new DormandPrince54Integrator(60.0, 3600.0, tol[0], tol[1]);
1127 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1128 for (DSSTForceModel force : forceModels) {
1129 propagator.addForceModel(force);
1130 }
1131 for (EventDetector event : events) {
1132 propagator.addEventDetector(event);
1133 }
1134 propagator.setInitialState(state);
1135
1136
1137 final SpacecraftState finalState = propagator.propagate(state.getDate().shiftedBy(86400.0));
1138
1139
1140 Assertions.assertEquals(finalState.getOrbit().getMu(), 3.986004415E14, Double.MIN_VALUE);
1141 }
1142
1143 @Test
1144 public void testIssue339WithAccelerations() {
1145 final SpacecraftState osculatingState = getLEOStatePropagatedBy30Minutes();
1146 final CelestialBody sun = CelestialBodyFactory.getSun();
1147 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
1148 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1149 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);
1150 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
1151 final AttitudeProvider attitudeProvider = new LofOffset(osculatingState.getFrame(), LOFType.LVLH_CCSDS, RotationOrder.XYZ, 0.0, 0.0, 0.0);
1152
1153 final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
1154 forces.add(new DSSTAtmosphericDrag(atmosphere, boxAndWing, osculatingState.getOrbit().getMu()));
1155 final SpacecraftState meanState = DSSTPropagator.computeMeanState(osculatingState, attitudeProvider, forces);
1156 final SpacecraftState computedOsculatingState = DSSTPropagator.computeOsculatingState(meanState, attitudeProvider, forces);
1157 Assertions.assertEquals(0.0,
1158 Vector3D.distance(osculatingState.getPosition(), computedOsculatingState.getPosition()),
1159 5.0e-6);
1160 }
1161
1162 @Test
1163 @Deprecated
1164 public void testIssue704() {
1165
1166
1167 final Orbit orbit = getLEOState().getOrbit();
1168 final PVCoordinates pv = orbit.getPVCoordinates();
1169
1170
1171 final double dP = 10.0;
1172
1173
1174 final double r2 = pv.getPosition().getNormSq();
1175 final double v = pv.getVelocity().getNorm();
1176 final double dV = orbit.getMu() * dP / (v * r2);
1177
1178
1179 final double[][] tol1 = DSSTPropagator.tolerances(dP, orbit);
1180 final double[][] tol2 = DSSTPropagator.tolerances(dP, dV, orbit);
1181 for (int i = 0; i < tol1.length; i++) {
1182 Assertions.assertArrayEquals(tol1[i], tol2[i], Double.MIN_VALUE);
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 @Test
1219 public void testIssue672() {
1220
1221
1222
1223
1224
1225 final UnnormalizedSphericalHarmonicsProvider provider =
1226 GravityFieldFactory.getUnnormalizedProvider(3, 3);
1227 final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
1228
1229
1230 final AbsoluteDate initDate = new AbsoluteDate(2007, 4, 16, 0, 46, 42.400,
1231 TimeScalesFactory.getUTC());
1232 final Orbit orbit = new KeplerianOrbit(26559890.,
1233 0.0041632,
1234 FastMath.toRadians(55.2),
1235 FastMath.toRadians(315.4985),
1236 FastMath.toRadians(130.7562),
1237 FastMath.toRadians(44.2377),
1238 PositionAngleType.MEAN,
1239 FramesFactory.getEME2000(),
1240 initDate,
1241 provider.getMu());
1242
1243
1244 final SpacecraftState initialState = new SpacecraftState(orbit);
1245 final double minStep = initialState.getOrbit().getKeplerianPeriod();
1246 final double maxStep = 100. * minStep;
1247 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(initialState.getOrbit(), OrbitType.EQUINOCTIAL);
1248 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1249 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1250 propagator.setInitialState(initialState, PropagationType.MEAN);
1251 propagator.addForceModel(new DSSTZonal(provider));
1252 propagator.addForceModel(new DSSTTesseral(earthFrame,
1253 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));
1254
1255
1256
1257
1258
1259 final SpacecraftState state = propagator.propagate(initDate.shiftedBy(1. * 86400.));
1260
1261
1262
1263
1264
1265 Assertions.assertNotNull(state);
1266 }
1267
1268
1269
1270
1271
1272 @Test
1273 public void testIssue1461() {
1274
1275
1276
1277
1278 final double step = 60.;
1279 final int nStep = 10;
1280 final Frame gcrf = FramesFactory.getGCRF();
1281 final AttitudeProvider attitude = new FrameAlignedProvider(gcrf);
1282
1283
1284 final Orbit initialOrbit = new KeplerianOrbit(7000e3, 1.e-3, FastMath.toRadians(98.), 0., 0., 0.,
1285 PositionAngleType.ECCENTRIC, gcrf,
1286 AbsoluteDate.ARBITRARY_EPOCH, Constants.EGM96_EARTH_MU);
1287 final AbsoluteDate t0 = initialOrbit.getDate();
1288 final SpacecraftState initialState = new SpacecraftState(initialOrbit);
1289 final ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(step);
1290
1291
1292 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.MEAN);
1293
1294
1295 propagator.setInitialState(initialState, PropagationType.OSCULATING);
1296
1297
1298 propagator.addForceModel(new DSSTZonal(GravityFieldFactory.getUnnormalizedProvider(2, 0)));
1299 propagator.setAttitudeProvider(attitude);
1300
1301
1302 final double initialMeanSma = DSSTPropagator.computeMeanState(initialState,
1303 attitude,
1304 propagator.getAllForceModels()).getOrbit().getA();
1305
1306
1307
1308
1309 final List<Double> propagatedMeanSma = new ArrayList<>();
1310 for (int i = 1; i <= nStep; i++) {
1311 propagatedMeanSma.add(propagator.propagate(t0.shiftedBy(i * step)).getOrbit().getA());
1312 }
1313
1314
1315
1316
1317
1318 for (Double meanSma : propagatedMeanSma) {
1319 Assertions.assertEquals(initialMeanSma, meanSma, 0.);
1320 }
1321 }
1322
1323
1324
1325
1326
1327
1328 @Test
1329 public void testJacobianWrtParametersIssue986() {
1330
1331 Frame eci = FramesFactory.getGCRF();
1332 AbsoluteDate epoch = AbsoluteDate.ARBITRARY_EPOCH;
1333 UnnormalizedSphericalHarmonicsProvider harmonics =
1334 GravityFieldFactory.getUnnormalizedProvider(4, 4);
1335 double mu = harmonics.getMu();
1336
1337 double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 500e3;
1338
1339
1340
1341
1342
1343 double e = 1e-17;
1344 Orbit orbit = new EquinoctialOrbit(new KeplerianOrbit(
1345 a, e, 1, 0, 0, 0,
1346 PositionAngleType.MEAN, eci, epoch, mu));
1347 final SpacecraftState initialState = new SpacecraftState(orbit);
1348 final double[][] tols = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(orbit, OrbitType.EQUINOCTIAL);
1349 ODEIntegrator integrator = new DormandPrince853Integrator(
1350 100 * 60,
1351 Constants.JULIAN_DAY,
1352 tols[0],
1353 tols[1]);
1354 final double dt = 2 * orbit.getKeplerianPeriod();
1355 AbsoluteDate end = epoch.shiftedBy(dt);
1356 final double n = orbit.getKeplerianMeanMotion();
1357 RealMatrix expectedStm = MatrixUtils.createRealIdentityMatrix(6);
1358
1359 final double dmda = -3.0 * n * dt / (2 * a);
1360 expectedStm.setEntry(5, 0, dmda);
1361 int ulps = 5;
1362 double absTol = 1e-18;
1363 final double[] expectedStmMdot = {0, 0, 0, 0, 0, dt};
1364 final Matcher<double[]> stmMdotMatcherExact = OrekitMatchers
1365 .doubleArrayContaining(expectedStmMdot, ulps);
1366 final Matcher<double[]> stmMdotMatcherClose = OrekitMatchers
1367 .doubleArrayContaining(expectedStmMdot, 1e-30, 15);
1368
1369 final double dmdmu = dt / (2 * Math.sqrt(mu * a * a * a));
1370 final Matcher<double[]> stmMuMatcherClose = OrekitMatchers
1371 .doubleArrayContaining(new double[]{0, 0, 0, 0, 0, dmdmu}, 1e-27, 15);
1372
1373 double twoParameterAbsTol = 1e-12;
1374
1375
1376
1377 double osculatingAbsTol = 1e-12;
1378
1379
1380 DSSTPropagator dsst = new DSSTPropagator(integrator, PropagationType.MEAN);
1381 dsst.setInitialState(initialState, PropagationType.MEAN);
1382
1383 dsst.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu));
1384
1385
1386 DSSTHarvester harvester = dsst.setupMatricesComputation("stm", null, null);
1387 SpacecraftState state = dsst.propagate(end);
1388 RealMatrix stm = harvester.getStateTransitionMatrix(state);
1389 RealMatrix stmParameters = harvester.getParametersJacobian(state);
1390
1391
1392 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1393 MatcherAssert.assertThat(
1394 harvester.getPositionAngleType(),
1395 is(PositionAngleType.MEAN));
1396 MatcherAssert.assertThat(stmParameters, nullValue());
1397 MatcherAssert.assertThat(stm, OrekitMatchers.matrixCloseTo(expectedStm, absTol));
1398
1399
1400 dsst = new DSSTPropagator(integrator, PropagationType.MEAN);
1401 dsst.setInitialState(initialState, PropagationType.MEAN);
1402 MDot mDot = new MDot();
1403 mDot.getParametersDrivers().get(0).setSelected(true);
1404 dsst.addForceModel(mDot);
1405
1406 dsst.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu));
1407
1408
1409 harvester = dsst.setupMatricesComputation("stm", null, null);
1410 state = dsst.propagate(end);
1411 stm = harvester.getStateTransitionMatrix(state);
1412 stmParameters = harvester.getParametersJacobian(state);
1413
1414
1415 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1416 MatcherAssert.assertThat(
1417 harvester.getPositionAngleType(),
1418 is(PositionAngleType.MEAN));
1419 MatcherAssert.assertThat(stm, OrekitMatchers.matrixCloseTo(expectedStm, absTol));
1420 MatcherAssert.assertThat(harvester.getJacobiansColumnsNames(), contains("MDot"));
1421 MatcherAssert.assertThat(stmParameters.getColumn(0), stmMdotMatcherExact);
1422 MatcherAssert.assertThat(stmParameters.getColumnDimension(), is(1));
1423 MatcherAssert.assertThat(stmParameters.getRowDimension(), is(6));
1424
1425
1426 dsst = new DSSTPropagator(integrator, PropagationType.MEAN);
1427 dsst.setInitialState(initialState, PropagationType.MEAN);
1428 mDot = new MDot();
1429 mDot.getParametersDrivers().get(0).setSelected(true);
1430 dsst.addForceModel(mDot);
1431
1432 DSSTThirdBody third = new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu);
1433 third.getParametersDrivers().get(1).setSelected(true);
1434 dsst.addForceModel(third);
1435
1436
1437 harvester = dsst.setupMatricesComputation("stm", null, null);
1438 state = dsst.propagate(end);
1439 stm = harvester.getStateTransitionMatrix(state);
1440 stmParameters = harvester.getParametersJacobian(state);
1441
1442
1443 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1444 MatcherAssert.assertThat(
1445 harvester.getPositionAngleType(),
1446 is(PositionAngleType.MEAN));
1447 MatcherAssert.assertThat(stm,
1448 OrekitMatchers.matrixCloseTo(expectedStm, twoParameterAbsTol));
1449
1450 MatcherAssert.assertThat(harvester.getJacobiansColumnsNames(),
1451 contains("MDot", "Spancentral attraction coefficient0"));
1452 MatcherAssert.assertThat(stmParameters.getColumn(0), stmMdotMatcherClose);
1453 MatcherAssert.assertThat(stmParameters.getColumn(1), stmMuMatcherClose);
1454 MatcherAssert.assertThat(stmParameters.getColumnDimension(), is(2));
1455 MatcherAssert.assertThat(stmParameters.getRowDimension(), is(6));
1456
1457
1458 dsst = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1459 dsst.setInitialState(initialState, PropagationType.MEAN);
1460
1461 dsst.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu));
1462
1463
1464 harvester = dsst.setupMatricesComputation("stm", null, null);
1465 harvester.initializeFieldShortPeriodTerms(initialState);
1466 harvester.updateFieldShortPeriodTerms(initialState);
1467 harvester.setReferenceState(initialState);
1468 state = dsst.propagate(end);
1469 stm = harvester.getStateTransitionMatrix(state);
1470 stmParameters = harvester.getParametersJacobian(state);
1471
1472
1473 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1474 MatcherAssert.assertThat(harvester.getPositionAngleType(), is(PositionAngleType.MEAN));
1475 MatcherAssert.assertThat(stmParameters, nullValue());
1476 MatcherAssert.assertThat(stm,
1477 OrekitMatchers.matrixCloseTo(expectedStm, osculatingAbsTol));
1478
1479
1480 dsst = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1481 dsst.setInitialState(initialState, PropagationType.MEAN);
1482 mDot = new MDot();
1483 mDot.getParametersDrivers().get(0).setSelected(true);
1484 dsst.addForceModel(mDot);
1485
1486 dsst.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu));
1487
1488
1489 harvester = dsst.setupMatricesComputation("stm", null, null);
1490 harvester.initializeFieldShortPeriodTerms(initialState);
1491 harvester.updateFieldShortPeriodTerms(initialState);
1492 harvester.setReferenceState(initialState);
1493 state = dsst.propagate(end);
1494 stm = harvester.getStateTransitionMatrix(state);
1495 stmParameters = harvester.getParametersJacobian(state);
1496
1497
1498 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1499 MatcherAssert.assertThat(harvester.getPositionAngleType(), is(PositionAngleType.MEAN));
1500 MatcherAssert.assertThat(stm,
1501 OrekitMatchers.matrixCloseTo(expectedStm, osculatingAbsTol));
1502 MatcherAssert.assertThat(harvester.getJacobiansColumnsNames(), contains("MDot"));
1503 MatcherAssert.assertThat(stmParameters.getColumn(0), stmMdotMatcherExact);
1504 MatcherAssert.assertThat(stmParameters.getColumnDimension(), is(1));
1505 MatcherAssert.assertThat(stmParameters.getRowDimension(), is(6));
1506
1507
1508 dsst = new DSSTPropagator(integrator, PropagationType.OSCULATING);
1509 dsst.setInitialState(initialState, PropagationType.MEAN);
1510 mDot = new MDot();
1511 mDot.getParametersDrivers().get(0).setSelected(true);
1512 dsst.addForceModel(mDot);
1513
1514 third = new DSSTThirdBody(CelestialBodyFactory.getPluto(), mu);
1515 third.getParametersDrivers().get(1).setSelected(true);
1516 dsst.addForceModel(third);
1517
1518
1519 harvester = dsst.setupMatricesComputation("stm", null, null);
1520 harvester.initializeFieldShortPeriodTerms(initialState);
1521 harvester.updateFieldShortPeriodTerms(initialState);
1522 harvester.setReferenceState(initialState);
1523 state = dsst.propagate(end);
1524 stm = harvester.getStateTransitionMatrix(state);
1525 stmParameters = harvester.getParametersJacobian(state);
1526
1527
1528 MatcherAssert.assertThat(harvester.getOrbitType(), is(OrbitType.EQUINOCTIAL));
1529 MatcherAssert.assertThat(
1530 harvester.getPositionAngleType(),
1531 is(PositionAngleType.MEAN));
1532 MatcherAssert.assertThat(stm,
1533 OrekitMatchers.matrixCloseTo(expectedStm, twoParameterAbsTol));
1534
1535 MatcherAssert.assertThat(harvester.getJacobiansColumnsNames(),
1536 contains("MDot", "Spancentral attraction coefficient0"));
1537 MatcherAssert.assertThat(stmParameters.getColumn(0), stmMdotMatcherClose);
1538 MatcherAssert.assertThat(stmParameters.getColumn(1), stmMuMatcherClose);
1539 MatcherAssert.assertThat(stmParameters.getColumnDimension(), is(2));
1540 MatcherAssert.assertThat(stmParameters.getRowDimension(), is(6));
1541 }
1542
1543
1544 private static class MDot implements DSSTForceModel {
1545
1546
1547 private final ParameterDriver mDot = new ParameterDriver(
1548 "MDot",
1549
1550 new TimeSpanMap<>("MDot"),
1551 new TimeSpanMap<>(0.0),
1552 0,
1553 1,
1554 Double.NEGATIVE_INFINITY,
1555 Double.POSITIVE_INFINITY);
1556
1557 @Override
1558 public List<ParameterDriver> getParametersDrivers() {
1559 return Collections.singletonList(mDot);
1560 }
1561
1562 @Override
1563 public double[] getMeanElementRate(SpacecraftState state,
1564 AuxiliaryElements auxiliaryElements,
1565 double[] parameters) {
1566 return new double[]{0, 0, 0, 0, 0, parameters[0]};
1567 }
1568
1569 @Override
1570 public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(
1571 FieldSpacecraftState<T> state,
1572 FieldAuxiliaryElements<T> auxiliaryElements,
1573 T[] parameters) {
1574
1575 final Field<T> field = state.getOrbit().getA().getField();
1576 final T zero = field.getZero();
1577 final T[] rates = MathArrays.buildArray(field, 6);
1578 rates[0] = zero;
1579 rates[1] = zero;
1580 rates[2] = zero;
1581 rates[3] = zero;
1582 rates[4] = zero;
1583 rates[5] = parameters[0];
1584 return rates;
1585 }
1586
1587
1588
1589 @Override
1590 public List<ShortPeriodTerms> initializeShortPeriodTerms(
1591 AuxiliaryElements auxiliaryElements,
1592 PropagationType type,
1593 double[] parameters) {
1594 return Collections.emptyList();
1595 }
1596
1597 @Override
1598 public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>>
1599 initializeShortPeriodTerms(FieldAuxiliaryElements<T> auxiliaryElements,
1600 PropagationType type,
1601 T[] parameters) {
1602 return Collections.emptyList();
1603 }
1604
1605 @Override
1606 public void registerAttitudeProvider(AttitudeProvider provider) {
1607 }
1608
1609 @Override
1610 public void updateShortPeriodTerms(double[] parameters,
1611 SpacecraftState... meanStates) {
1612 }
1613
1614 @SuppressWarnings("unchecked")
1615 @Override
1616 public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(
1617 T[] parameters, FieldSpacecraftState<T>... meanStates) {
1618 }
1619
1620 }
1621
1622 private SpacecraftState getGEOState() throws IllegalArgumentException, OrekitException {
1623
1624 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 05, 21), new TimeComponents(1, 0, 0.),
1625 TimeScalesFactory.getUTC());
1626 final Orbit orbit = new EquinoctialOrbit(42164000,
1627 10e-3,
1628 10e-3,
1629 FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3),
1630 FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1,
1631 PositionAngleType.TRUE,
1632 FramesFactory.getEME2000(),
1633 initDate,
1634 3.986004415E14);
1635 return new SpacecraftState(orbit);
1636 }
1637
1638 private SpacecraftState getLEOState() throws IllegalArgumentException, OrekitException {
1639 final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1640 final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1641
1642 final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1643 return new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
1644 FramesFactory.getEME2000(),
1645 initDate,
1646 3.986004415E14));
1647 }
1648
1649 private void setDSSTProp(SpacecraftState initialState) {
1650 initialState.getDate();
1651 final double minStep = initialState.getOrbit().getKeplerianPeriod();
1652 final double maxStep = 100. * minStep;
1653 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1.).getTolerances(initialState.getOrbit(), OrbitType.EQUINOCTIAL);
1654 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1655 dsstProp = new DSSTPropagator(integrator);
1656 dsstProp.setInitialState(initialState, PropagationType.MEAN);
1657 }
1658
1659 private static class CheckingHandler implements EventHandler {
1660
1661 private final Action actionOnEvent;
1662 private boolean gotHere;
1663
1664 public CheckingHandler(final Action actionOnEvent) {
1665 this.actionOnEvent = actionOnEvent;
1666 this.gotHere = false;
1667 }
1668
1669 public void assertEvent(boolean expected) {
1670 Assertions.assertEquals(expected, gotHere);
1671 }
1672
1673 @Override
1674 public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
1675 gotHere = true;
1676 return actionOnEvent;
1677 }
1678
1679 }
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689 @Override
1690 protected List<ParameterDriver> getParametersDriversWithoutMu() {
1691 return Collections.emptyList();
1692 }
1693
1694
1695 @Override
1696 protected double[] getLLimits(SpacecraftState state,
1697 AuxiliaryElements auxiliaryElements) {
1698 return new double[] { -FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0),
1699 FastMath.PI + MathUtils.normalizeAngle(state.getOrbit().getLv(), 0) };
1700 }
1701
1702
1703 @Override
1704 protected <T extends CalculusFieldElement<T>> T[] getLLimits(FieldSpacecraftState<T> state,
1705 FieldAuxiliaryElements<T> auxiliaryElements) {
1706 final Field<T> field = state.getDate().getField();
1707 final T zero = field.getZero();
1708 final T[] tab = MathArrays.buildArray(field, 2);
1709 tab[0] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).subtract(FastMath.PI);
1710 tab[1] = MathUtils.normalizeAngle(state.getOrbit().getLv(), zero).add(FastMath.PI);
1711 return tab;
1712 }
1713
1714 }
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727 @Override
1728 public void init(final SpacecraftState s0, final AbsoluteDate t) {
1729 this.initialized = true;
1730 this.accComputed = false;
1731 }
1732
1733
1734 @Override
1735 public boolean dependsOnPositionOnly() {
1736 return false;
1737 }
1738
1739
1740 @Override
1741 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1742 this.accComputed = true;
1743 return Vector3D.ZERO;
1744 }
1745
1746
1747 @Override
1748 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1749 return FieldVector3D.getZero(s.getDate().getField());
1750 }
1751
1752
1753 @Override
1754 public Stream<EventDetector> getEventDetectors() {
1755 return Stream.empty();
1756 }
1757
1758
1759 @Override
1760 public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
1761 return Stream.empty();
1762 }
1763
1764
1765 @Override
1766 public List<ParameterDriver> getParametersDrivers() {
1767 return Collections.emptyList();
1768 }
1769
1770 }
1771
1772 @BeforeEach
1773 public void setUp() throws IOException, ParseException {
1774 Utils.setDataRoot("regular-data:potential/shm-format");
1775 }
1776
1777 @AfterEach
1778 public void tearDown() {
1779 dsstProp = null;
1780 }
1781
1782 private SpacecraftState getLEOStatePropagatedBy30Minutes() throws IllegalArgumentException, OrekitException {
1783
1784 final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
1785 final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
1786
1787 final AbsoluteDate initialDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
1788 final CartesianOrbit osculatingOrbit = new CartesianOrbit(new PVCoordinates(position, velocity), FramesFactory.getTOD(IERSConventions.IERS_1996, false),
1789 initialDate, Constants.WGS84_EARTH_MU);
1790
1791
1792 double minStep = 0.001;
1793 double maxstep = 1000.0;
1794 double positionTolerance = 10.0;
1795 OrbitType propagationType = OrbitType.EQUINOCTIAL;
1796 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(osculatingOrbit, propagationType);
1797 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
1798 NumericalPropagator propagator = new NumericalPropagator(integrator);
1799 propagator.setOrbitType(propagationType);
1800
1801 NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
1802 ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
1803 propagator.addForceModel(holmesFeatherstone);
1804 propagator.setInitialState(new SpacecraftState(osculatingOrbit));
1805
1806 return propagator.propagate(new AbsoluteDate(initialDate, 1800.));
1807 }
1808
1809 @Test
1810 void testIssue1788() {
1811
1812 final ClassicalRungeKuttaIntegrator integrator = new ClassicalRungeKuttaIntegrator(1000);
1813 final CountAndContinue countAndContinue = new CountAndContinue();
1814 final KeplerianOrbit orbit = new KeplerianOrbit(42166000.0, 0.00028, FastMath.toRadians(0.05), FastMath.toRadians(66.0),
1815 FastMath.toRadians(270.0), FastMath.toRadians(11.94), PositionAngleType.MEAN,
1816 FramesFactory.getGCRF(), AbsoluteDate.ARBITRARY_EPOCH, Constants.WGS84_EARTH_MU);
1817 final SpacecraftState initialState = new SpacecraftState(orbit);
1818 final double timeOfFlight = 1e4;
1819 final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING,
1820 new TestAttitudeProvider(countAndContinue, initialState.getDate().shiftedBy(timeOfFlight/2)));
1821 propagator.setInitialState(initialState);
1822
1823 propagator.propagate(initialState.getDate().shiftedBy(timeOfFlight));
1824
1825 Assertions.assertEquals(1, countAndContinue.getCount());
1826 }
1827
1828 private static class TestAttitudeProvider implements AttitudeProviderModifier {
1829
1830 private final DateDetector detector;
1831
1832 TestAttitudeProvider(final CountAndContinue countingHandler,
1833 final AbsoluteDate date) {
1834 this.detector = new DateDetector(date).withHandler(countingHandler);
1835 }
1836
1837 @Override
1838 public Stream<EventDetector> getEventDetectors() {
1839 return Stream.of(detector);
1840 }
1841
1842 @Override
1843 public AttitudeProvider getUnderlyingAttitudeProvider() {
1844 return new FrameAlignedProvider(FramesFactory.getEME2000());
1845 }
1846 }
1847 }