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