1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.conversion.osc2mean;
18
19 import java.util.concurrent.TimeUnit;
20
21 import org.hipparchus.exception.MathIllegalStateException;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer;
26 import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
27 import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
28 import org.hipparchus.stat.descriptive.rank.Max;
29 import org.hipparchus.stat.descriptive.rank.Min;
30 import org.hipparchus.util.Binary64Field;
31 import org.hipparchus.util.FastMath;
32 import org.junit.jupiter.api.AfterEach;
33 import org.junit.jupiter.api.Assertions;
34 import org.junit.jupiter.api.BeforeEach;
35 import org.junit.jupiter.api.Test;
36 import org.orekit.Utils;
37 import org.orekit.errors.OrekitException;
38 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
39 import org.orekit.forces.gravity.potential.GravityFieldFactory;
40 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
41 import org.orekit.frames.Frame;
42 import org.orekit.frames.FramesFactory;
43 import org.orekit.orbits.FieldOrbit;
44 import org.orekit.orbits.KeplerianOrbit;
45 import org.orekit.orbits.Orbit;
46 import org.orekit.orbits.OrbitType;
47 import org.orekit.orbits.PositionAngleType;
48 import org.orekit.propagation.PropagationType;
49 import org.orekit.propagation.Propagator;
50 import org.orekit.propagation.SpacecraftState;
51 import org.orekit.propagation.ToleranceProvider;
52 import org.orekit.propagation.analytical.BrouwerLyddanePropagator;
53 import org.orekit.propagation.analytical.tle.TLE;
54 import org.orekit.propagation.analytical.tle.TLEPropagator;
55 import org.orekit.propagation.numerical.NumericalPropagator;
56 import org.orekit.time.AbsoluteDate;
57 import org.orekit.utils.Constants;
58 import org.orekit.utils.IERSConventions;
59
60 public class BrouwerLyddaneMeanConversionTest {
61
62 private UnnormalizedSphericalHarmonicsProvider provider;
63
64 private LeastSquaresOptimizer optimizer;
65
66 private boolean doPrint;
67
68 @BeforeEach
69 public void setUp() {
70 Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
71 provider = GravityFieldFactory.getUnnormalizedProvider(5, 0);
72 optimizer = new LevenbergMarquardtOptimizer();
73 doPrint = false;
74 }
75
76 @AfterEach
77 public void tearDown() {
78 provider = null;
79 optimizer = null;
80 }
81
82 @Test
83 void testIssue947() {
84
85 final TLE tleOrbit = new TLE("1 43196U 18015E 21055.59816856 .00000894 00000-0 38966-4 0 9996",
86 "2 43196 97.4662 188.8169 0016935 299.6845 60.2706 15.24746686170319");
87 final Propagator propagator = TLEPropagator.selectExtrapolator(tleOrbit);
88
89
90 SpacecraftState tleState = propagator.getInitialState();
91 final KeplerianOrbit osculating0 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(tleState.getOrbit());
92 SpacecraftState tleStateAtDate = propagator.propagate(propagator.getInitialState().getDate().shiftedBy(3, TimeUnit.DAYS));
93 final KeplerianOrbit osculating1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(tleStateAtDate.getOrbit());
94
95
96 final BrouwerLyddaneTheory theory = new BrouwerLyddaneTheory(provider,
97 BrouwerLyddanePropagator.M2);
98
99 final FixedPointConverter fpConvert = new FixedPointConverter(theory);
100 final LeastSquaresConverter lsConvert = new LeastSquaresConverter(theory, optimizer);
101
102
103 Assertions.assertDoesNotThrow(() -> {
104 final Orbit mean0fp = fpConvert.convertToMean(osculating0);
105 final Orbit mean0ls = lsConvert.convertToMean(osculating0);
106 if (doPrint) {
107 System.out.println("osculating : " + osculating0);
108 System.out.println("mean (fp) : " + mean0fp);
109 System.out.println("mean (ls) : " + mean0ls);
110 System.out.println("dP (m) : " + Vector3D.distance(mean0ls.getPosition(mean0ls.getFrame()),
111 mean0fp.getPosition(mean0ls.getFrame())));
112 System.out.println("iter : " + lsConvert.getIterationsNb() + " ; rms : " + lsConvert.getRMS());
113 System.out.println();
114 }
115 });
116
117 Assertions.assertDoesNotThrow(() -> {
118 final Orbit mean1fp1 = fpConvert.convertToMean(osculating1);
119 final Orbit mean1ls1 = lsConvert.convertToMean(osculating1);
120 if (doPrint) {
121 System.out.println("osculating : " + osculating1);
122 System.out.println("mean (fp) : " + mean1fp1);
123 System.out.println("mean (ls) : " + mean1ls1);
124 System.out.println("dP (m) : " + Vector3D.distance(mean1ls1.getPosition(mean1ls1.getFrame()),
125 mean1fp1.getPosition(mean1ls1.getFrame())));
126 System.out.println("iter : " + lsConvert.getIterationsNb() + " ; rms : " + lsConvert.getRMS());
127 System.out.println();
128 }
129 });
130 }
131
132 @Test
133 void testIssue1558() {
134
135 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH;
136 Orbit initialOrbit = new KeplerianOrbit(67679244.0, 0.0, 1.85850, 2.1, 2.9, 3.2,
137 PositionAngleType.TRUE, FramesFactory.getEME2000(),
138 initDate, Constants.EGM96_EARTH_MU);
139 if (doPrint) {
140 System.out.println(initialOrbit);
141 System.out.println();
142 }
143
144
145 BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(initialOrbit, provider,
146 PropagationType.MEAN,
147 BrouwerLyddanePropagator.M2);
148
149
150 final double shift = 3 * Constants.JULIAN_DAY;
151 final Orbit initialState = propagator.propagate(initDate).getOrbit();
152 final Orbit futureState = propagator.propagate(initDate.shiftedBy(shift)).getOrbit();
153 final Binary64Field field = Binary64Field.getInstance();
154 final FieldOrbit<?> initialField = initialState.getType().convertToFieldOrbit(field, initialState);
155 final FieldOrbit<?> futureField = futureState.getType().convertToFieldOrbit(field, futureState);
156
157
158 final BrouwerLyddaneTheory theory = new BrouwerLyddaneTheory(provider,
159 BrouwerLyddanePropagator.M2);
160
161 final FixedPointConverter fpConvert = new FixedPointConverter(theory);
162 final LeastSquaresConverter lsConvert = new LeastSquaresConverter(theory, optimizer);
163
164
165 OrekitException oei = Assertions.assertThrows(OrekitException.class, () -> {
166 fpConvert.convertToMean(initialState);
167 });
168 Assertions.assertTrue(oei.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
169 oei = Assertions.assertThrows(OrekitException.class, () -> {
170 fpConvert.convertToMean(initialField);
171 });
172 Assertions.assertTrue(oei.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
173 OrekitException oef = Assertions.assertThrows(OrekitException.class, () -> {
174 fpConvert.convertToMean(futureState);
175 });
176 Assertions.assertTrue(oef.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
177 oef = Assertions.assertThrows(OrekitException.class, () -> {
178 fpConvert.convertToMean(futureField);
179 });
180 Assertions.assertTrue(oef.getMessage().contains("unable to compute Brouwer-Lyddane mean parameters after"));
181
182
183 Assertions.assertDoesNotThrow(() -> {
184 final Orbit mean0ls = lsConvert.convertToMean(initialState);
185 if (doPrint) {
186 System.out.println(initialState);
187 System.out.println(mean0ls);
188 System.out.println("iter : " + lsConvert.getIterationsNb() + " ; rms : " + lsConvert.getRMS());
189 System.out.println();
190 }
191 });
192
193 Assertions.assertDoesNotThrow(() -> {
194 final Orbit mean1ls = lsConvert.convertToMean(futureState);
195 if (doPrint) {
196 System.out.println(futureState);
197 System.out.println(mean1ls);
198 System.out.println("iter : " + lsConvert.getIterationsNb() + " ; rms : " + lsConvert.getRMS());
199 System.out.println();
200 }
201 });
202 }
203
204 @Test
205 void testMeanOrbit() {
206 final KeplerianOrbit initialOsculating =
207 new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngleType.TRUE,
208 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
209 provider.getMu());
210
211
212 final BrouwerLyddaneTheory theory = new BrouwerLyddaneTheory(provider,
213 BrouwerLyddanePropagator.M2);
214
215 final FixedPointConverter fpConvert = new FixedPointConverter(theory);
216 final LeastSquaresConverter lsConvert = new LeastSquaresConverter(theory, optimizer);
217
218
219
220 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(0.1)
221 .getTolerances(initialOsculating, OrbitType.KEPLERIAN);
222 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
223 integrator.setInitialStepSize(60);
224 NumericalPropagator num = new NumericalPropagator(integrator);
225 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
226 num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider)));
227 num.setInitialState(new SpacecraftState(initialOsculating));
228 num.setOrbitType(OrbitType.KEPLERIAN);
229 num.setPositionAngleType(initialOsculating.getCachedPositionAngleType());
230 final StorelessUnivariateStatistic oscMin = new Min();
231 final StorelessUnivariateStatistic oscMax = new Max();
232 final StorelessUnivariateStatistic meanFpMin = new Min();
233 final StorelessUnivariateStatistic meanFpMax = new Max();
234 final StorelessUnivariateStatistic meanLsMin = new Min();
235 final StorelessUnivariateStatistic meanLsMax = new Max();
236 num.getMultiplexer().add(60, state -> {
237 final Orbit osc = state.getOrbit();
238 oscMin.increment(osc.getA());
239 oscMax.increment(osc.getA());
240
241 final Orbit mean0 = fpConvert.convertToMean(state.getOrbit());
242 final Orbit mean1 = lsConvert.convertToMean(state.getOrbit());
243 meanFpMin.increment(mean0.getA());
244 meanFpMax.increment(mean0.getA());
245 meanLsMin.increment(mean1.getA());
246 meanLsMax.increment(mean1.getA());
247 if (doPrint) {
248 System.out.println("dP (m) : " + Vector3D.distance(mean0.getPosition(mean0.getFrame()),
249 mean1.getPosition(mean0.getFrame())));
250 System.out.println("iter : " + lsConvert.getIterationsNb() + " ; rms : " + lsConvert.getRMS());
251 System.out.println();
252 }
253 });
254 num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
255
256
257 Assertions.assertEquals(3188.347, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
258 Assertions.assertEquals( 25.794, meanFpMax.getResult() - meanFpMin.getResult(), 1.0e-3);
259 Assertions.assertEquals( 25.794, meanLsMax.getResult() - meanLsMin.getResult(), 1.0e-3);
260
261 }
262
263 @Test
264 void testGeostationaryOrbit() {
265
266 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
267 final Frame eci = FramesFactory.getEME2000();
268
269
270 final double sma = FastMath.cbrt(Constants.IERS2010_EARTH_MU /
271 Constants.IERS2010_EARTH_ANGULAR_VELOCITY /
272 Constants.IERS2010_EARTH_ANGULAR_VELOCITY);
273 final double ecc = 0.0;
274 final double inc = 0.0;
275 final double pa = 0.0;
276 final double raan = 0.0;
277 final double lV = 0.0;
278 final Orbit orbit = new KeplerianOrbit(sma, ecc, inc, pa, raan, lV,
279 PositionAngleType.TRUE,
280 eci, date, provider.getMu());
281
282
283 final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(orbit, provider, PropagationType.MEAN,
284 BrouwerLyddanePropagator.M2);
285
286
287 final SpacecraftState state = bl.propagate(date);
288 final KeplerianOrbit orbOsc = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state.getOrbit());
289
290 if (doPrint) {
291 System.out.println(orbit);
292 System.out.println();
293 System.out.println(orbOsc);
294 }
295
296
297 final BrouwerLyddaneTheory theory = new BrouwerLyddaneTheory(provider,
298 BrouwerLyddanePropagator.M2);
299
300 final FixedPointConverter fpConvert = new FixedPointConverter(theory,
301 FixedPointConverter.DEFAULT_THRESHOLD,
302 FixedPointConverter.DEFAULT_MAX_ITERATIONS,
303 0.75);
304 final LeastSquaresConverter lsConvert = new LeastSquaresConverter(theory, optimizer);
305
306
307 Assertions.assertDoesNotThrow(() -> {
308 Orbit orbMean = fpConvert.convertToMean(orbOsc);
309 if (doPrint) {
310 System.out.println(orbMean);
311 System.out.println();
312 }
313 });
314
315 MathIllegalStateException mise = Assertions.assertThrows(MathIllegalStateException.class, () -> {
316 lsConvert.convertToMean(orbOsc);
317 });
318 Assertions.assertTrue(mise.getMessage().contains("unable to perform Q.R decomposition"));
319
320
321 final BrouwerLyddanePropagator bl2 = new BrouwerLyddanePropagator(orbit, provider,
322 PropagationType.OSCULATING,
323 BrouwerLyddanePropagator.M2);
324
325
326 final SpacecraftState state2 = bl2.propagate(date.shiftedBy(Constants.JULIAN_DAY));
327 final KeplerianOrbit orbOsc2 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(state2.getOrbit());
328 if (doPrint) {
329 System.out.println(orbOsc2);
330 }
331 Assertions.assertDoesNotThrow(() -> {
332 Orbit orbMean2 = fpConvert.convertToMean(orbOsc2);
333 if (doPrint) {
334 System.out.println(orbMean2);
335 System.out.println();
336 }
337 });
338
339 mise = Assertions.assertThrows(MathIllegalStateException.class, () -> {
340 lsConvert.convertToMean(orbOsc2);
341 });
342 Assertions.assertTrue(mise.getMessage().contains("unable to perform Q.R decomposition"));
343 }
344 }