1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.orekit.estimation.iod;
19
20 import java.util.List;
21
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.util.FastMath;
24 import org.junit.Assert;
25 import org.junit.Before;
26 import org.junit.Test;
27 import org.orekit.Utils;
28 import org.orekit.errors.OrekitException;
29 import org.orekit.errors.OrekitMessages;
30 import org.orekit.estimation.Context;
31 import org.orekit.estimation.EstimationTestUtils;
32 import org.orekit.estimation.measurements.ObservableSatellite;
33 import org.orekit.estimation.measurements.ObservedMeasurement;
34 import org.orekit.estimation.measurements.PVMeasurementCreator;
35 import org.orekit.estimation.measurements.Position;
36 import org.orekit.frames.Frame;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.orbits.KeplerianOrbit;
39 import org.orekit.orbits.OrbitType;
40 import org.orekit.orbits.PositionAngle;
41 import org.orekit.propagation.Propagator;
42 import org.orekit.propagation.SpacecraftState;
43 import org.orekit.propagation.analytical.KeplerianPropagator;
44 import org.orekit.propagation.analytical.tle.TLE;
45 import org.orekit.propagation.analytical.tle.TLEPropagator;
46 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
47 import org.orekit.time.AbsoluteDate;
48 import org.orekit.time.TimeScale;
49 import org.orekit.time.TimeScalesFactory;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.PVCoordinates;
52
53
54
55
56
57
58
59
60
61 public class IodLambertTest {
62
63 @Test
64 public void testLambert() {
65 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
66
67 final double mu = context.initialOrbit.getMu();
68 final Frame frame = context.initialOrbit.getFrame();
69
70 final NumericalPropagatorBuilder propagatorBuilder =
71 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
72 1.0e-6, 60.0, 0.001);
73
74
75 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
76 propagatorBuilder);
77
78 final List<ObservedMeasurement<?>> measurements =
79 EstimationTestUtils.createMeasurements(propagator,
80 new PVMeasurementCreator(),
81 0.0, 1.0, 60.0);
82
83
84 final int idMeasure1 = 0;
85 final AbsoluteDate date1 = measurements.get(idMeasure1).getDate();
86
87
88
89
90 final Vector3D position1 = new Vector3D(
91 measurements.get(idMeasure1).getObservedValue()[0],
92 measurements.get(idMeasure1).getObservedValue()[1],
93 measurements.get(idMeasure1).getObservedValue()[2]);
94
95
96 final int idMeasure2 = 10;
97 final AbsoluteDate date2 = measurements.get(idMeasure2).getDate();
98
99
100
101
102 final Vector3D position2 = new Vector3D(
103 measurements.get(idMeasure2).getObservedValue()[0],
104 measurements.get(idMeasure2).getObservedValue()[1],
105 measurements.get(idMeasure2).getObservedValue()[2]);
106
107 final int nRev = 0;
108
109
110 final IodLambert iod = new IodLambert(mu);
111
112 final KeplerianOrbit orbit = iod.estimate(frame,
113 true,
114 nRev,
115 (position1), date1,
116 (position2), date2);
117
118 Assert.assertEquals(orbit.getA(), context.initialOrbit.getA(), 1.0e-9 * context.initialOrbit.getA());
119 Assert.assertEquals(orbit.getE(), context.initialOrbit.getE(), 1.0e-9 * context.initialOrbit.getE());
120 Assert.assertEquals(orbit.getI(), context.initialOrbit.getI(), 1.0e-9 * context.initialOrbit.getI());
121 }
122
123
124
125
126 @Test
127 public void testLambert2() {
128
129
130 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
131 final double mu = context.initialOrbit.getMu();
132 final Frame frame = context.initialOrbit.getFrame();
133
134
135 final NumericalPropagatorBuilder propagatorBuilder =
136 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
137 1.0e-6, 60.0, 0.001);
138
139
140 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
141 propagatorBuilder);
142
143 final SpacecraftState initialState = propagator.getInitialState();
144 final KeplerianOrbit refOrbit = new KeplerianOrbit(initialState.getOrbit());
145 double[] refOrbitArray = new double[6];
146 OrbitType.KEPLERIAN.mapOrbitToArray(refOrbit, PositionAngle.TRUE, refOrbitArray, null);
147 final Vector3D position1 = refOrbit.getPVCoordinates().getPosition();
148 final AbsoluteDate date1 = refOrbit.getDate();
149
150
151 final double T = context.initialOrbit.getKeplerianPeriod();
152
153
154
155
156
157
158
159
160 final double[] dts = new double[] {
161 T/4, T/2, 3*T/4, 2*T + T/4, 3*T + 3*T/4};
162 final int[] nRevs = new int[] {0, 0, 0, 2, 3};
163 final boolean[] posigrades = new boolean[] {true, true, false, true, false};
164
165 for (int i = 0; i < dts.length; i++) {
166
167 propagator.resetInitialState(initialState);
168
169
170 final AbsoluteDate date2 = date1.shiftedBy(dts[i]);
171 final Vector3D position2 = propagator.propagate(date2).getPVCoordinates().getPosition();
172
173
174 final IodLambert iod = new IodLambert(mu);
175
176
177 final KeplerianOrbit orbit = iod.estimate(frame, posigrades[i], nRevs[i], position1, date1, position2, date2);
178
179
180 final double relTol = 1e-12;
181 Assert.assertEquals(refOrbit.getA(), orbit.getA(), relTol * refOrbit.getA());
182 Assert.assertEquals(refOrbit.getE(), orbit.getE(), relTol * refOrbit.getE());
183 Assert.assertEquals(refOrbit.getI(), orbit.getI(), relTol * refOrbit.getI());
184 Assert.assertEquals(refOrbit.getPerigeeArgument(), orbit.getPerigeeArgument(), relTol * refOrbit.getPerigeeArgument());
185 Assert.assertEquals(refOrbit.getRightAscensionOfAscendingNode(), orbit.getRightAscensionOfAscendingNode(), relTol * refOrbit.getRightAscensionOfAscendingNode());
186 Assert.assertEquals(refOrbit.getTrueAnomaly(), orbit.getTrueAnomaly(), relTol * refOrbit.getTrueAnomaly());
187 }
188 }
189
190
191 @Test
192 public void testMultiRevolutions() {
193
194 Utils.setDataRoot("regular-data");
195 TLE aussatB1 = new TLE("1 22087U 92054A 17084.21270512 -.00000243 +00000-0 +00000-0 0 9999",
196 "2 22087 008.5687 046.5717 0005960 022.3650 173.1619 00.99207999101265");
197 final Propagator propagator = TLEPropagator.selectExtrapolator(aussatB1);
198 final Frame teme = FramesFactory.getTEME();
199
200 final AbsoluteDate t1 = new AbsoluteDate("2017-03-25T23:48:31.282", TimeScalesFactory.getUTC());
201 final Vector3D p1 = propagator.propagate(t1).getPVCoordinates(teme).getPosition();
202 final AbsoluteDate t2 = t1.shiftedBy(40000);
203 final Vector3D p2 = propagator.propagate(t2).getPVCoordinates(teme).getPosition();
204 final AbsoluteDate t3 = t1.shiftedBy(115200.0);
205 final Vector3D p3 = propagator.propagate(t3).getPVCoordinates(teme).getPosition();
206
207 IodLambert lambert = new IodLambert(Constants.WGS84_EARTH_MU);
208 KeplerianOrbit k0 = lambert.estimate(teme, true, 0, p1, t1, p2, t2);
209 Assert.assertEquals(6.08e-4, k0.getE(), 1.0e-6);
210 Assert.assertEquals(8.55, FastMath.toDegrees(k0.getI()), 0.01);
211 Assert.assertEquals(0.0, Vector3D.distance(p1, k0.getPVCoordinates(t1, teme).getPosition()), 2.0e-8);
212 Assert.assertEquals(0.0, Vector3D.distance(p2, k0.getPVCoordinates(t2, teme).getPosition()), 2.0e-7);
213
214 KeplerianOrbit k1 = lambert.estimate(teme, true, 1, p1, t1, p3, t3);
215 Assert.assertEquals(5.97e-4, k1.getE(), 1.0e-6);
216 Assert.assertEquals(8.55, FastMath.toDegrees(k1.getI()), 0.01);
217 Assert.assertEquals(0.0, Vector3D.distance(p1, k1.getPVCoordinates(t1, teme).getPosition()), 1.4e-8);
218 Assert.assertEquals(0.0, Vector3D.distance(p3, k1.getPVCoordinates(t3, teme).getPosition()), 3.0e-7);
219
220 }
221
222 @Test
223 public void testNonChronologicalObservations() {
224
225
226 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
227 final double mu = context.initialOrbit.getMu();
228 final Frame frame = context.initialOrbit.getFrame();
229
230
231 final NumericalPropagatorBuilder propagatorBuilder =
232 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
233 1.0e-6, 60.0, 0.001);
234
235
236 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
237 propagatorBuilder);
238
239 final SpacecraftState initialState = propagator.getInitialState();
240 final KeplerianOrbit refOrbit = new KeplerianOrbit(initialState.getOrbit());
241 double[] refOrbitArray = new double[6];
242 OrbitType.KEPLERIAN.mapOrbitToArray(refOrbit, PositionAngle.TRUE, refOrbitArray, null);
243
244 final Vector3D position1 = refOrbit.getPVCoordinates().getPosition();
245 final AbsoluteDate date1 = refOrbit.getDate();
246
247
248 final double T = context.initialOrbit.getKeplerianPeriod();
249
250 final double dts = -T/4;
251 final int nRevs = 0;
252 final boolean posigrades = true;
253
254
255 propagator.resetInitialState(initialState);
256
257
258 final AbsoluteDate date2 = date1.shiftedBy(dts);
259 final Vector3D position2 = propagator.propagate(date2).getPVCoordinates().getPosition();
260
261
262 final IodLambert iod = new IodLambert(mu);
263
264
265 try {
266 iod.estimate(frame, posigrades, nRevs, position1, date1, position2, date2);
267 Assert.fail("An exception should have been thrown");
268
269 } catch (OrekitException oe) {
270 Assert.assertEquals(OrekitMessages.NON_CHRONOLOGICAL_DATES_FOR_OBSERVATIONS, oe.getSpecifier());
271 }
272 }
273
274
275
276
277
278
279
280
281 @Test
282 public void testIssue533() {
283
284
285
286 double mu = Constants.EGM96_EARTH_MU;
287 Frame j2000 = FramesFactory.getEME2000();
288
289
290
291
292
293
294
295
296
297
298 for (int testCase = 1; testCase < 3; testCase++) {
299
300 double trueAnomalyDifference;
301 boolean posigrade;
302 int nRev;
303 if (testCase == 1) {
304 posigrade = true;
305 nRev = 0;
306 trueAnomalyDifference = 0.5*FastMath.PI;
307 } else {
308 posigrade = false;
309 nRev = 0;
310 trueAnomalyDifference = 1.5*FastMath.PI;
311 }
312
313
314 double a = 24000*1e3;
315 double e = 0.72;
316 double i = FastMath.toRadians(0.);
317 double aop = 0.;
318 double raan = 0.;
319 double nu0 = 0.;
320
321
322 AbsoluteDate t1 = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
323 KeplerianOrbit kep1 = new KeplerianOrbit(a, e, i, aop, raan, nu0, PositionAngle.TRUE, j2000, t1, mu);
324 Vector3D p1 = kep1.getPVCoordinates().getPosition();
325 Vector3D v1 = kep1.getPVCoordinates().getVelocity();
326
327
328 KeplerianOrbit kep2 = new KeplerianOrbit(a, e, i, aop, raan, trueAnomalyDifference, PositionAngle.TRUE, j2000, t1, mu);
329 double n = kep2.getKeplerianMeanMotion();
330 double M2 = kep2.getMeanAnomaly();
331 double delta_t = M2 / n;
332 AbsoluteDate t2 = t1.shiftedBy(delta_t);
333
334
335 PVCoordinates pv2 = kep1.getPVCoordinates(t2, j2000);
336 Vector3D p2 = pv2.getPosition();
337 Vector3D v2 = pv2.getVelocity();
338
339
340 IodLambert iodLambert = new IodLambert(mu);
341 KeplerianOrbit resultOrbit1 = iodLambert.estimate(j2000, posigrade, nRev, p1, t1, p2, t2);
342
343
344 PVCoordinates resultPv1 = resultOrbit1.getPVCoordinates(t1, j2000);
345 PVCoordinates resultPv2 = resultOrbit1.getPVCoordinates(t2, j2000);
346
347 Vector3D resultP1 = resultPv1.getPosition();
348 Vector3D resultV1 = resultPv1.getVelocity();
349 Vector3D resultP2 = resultPv2.getPosition();
350 Vector3D resultV2 = resultPv2.getVelocity();
351
352
353 double dP1 = Vector3D.distance(p1, resultP1);
354 double dV1 = Vector3D.distance(v1, resultV1);
355 double dP2 = Vector3D.distance(p2, resultP2);
356 double dV2 = Vector3D.distance(v2, resultV2);
357
358
359 double dP1Tol, dP2Tol, dV1Tol, dV2Tol;
360 if (testCase == 1) {
361 dP1Tol = 4.28e-25;
362 dV1Tol = 5.97e-12;
363 dP2Tol = 7.57e-9;
364 dV2Tol = 7.50e-12;
365 } else {
366 dP1Tol = 2.81e-25;
367 dV1Tol = 3.03e-12;
368 dP2Tol = 9.86e-7;
369 dV2Tol = 4.01e-10;
370 }
371
372
373 Assert.assertEquals(0., dP1, dP1Tol);
374 Assert.assertEquals(0., dV1, dV1Tol);
375 Assert.assertEquals(0., dP2, dP2Tol);
376 Assert.assertEquals(0., dV2, dV2Tol);
377 }
378 }
379
380 @Test
381 public void testIssue752() {
382
383
384
385
386 final Frame inertialFrame = FramesFactory.getEME2000();
387 final TimeScale utc = TimeScalesFactory.getUTC();
388
389
390 final IodLambert lambert = new IodLambert(Constants.EGM96_EARTH_MU);
391
392
393 final ObservableSatellite satellite = new ObservableSatellite(0);
394
395
396 final Vector3D posR1 = new Vector3D(22592145.603, -1599915.239, -19783950.506);
397 final Vector3D posR2 = new Vector3D(1922067.697, 4054157.051, -8925727.465);
398
399
400 AbsoluteDate dateRef = new AbsoluteDate(2018, 11, 1, 0, 0, 0.0, utc);
401 AbsoluteDate date2 = dateRef.shiftedBy(36000.0);
402
403
404 final Vector3D velR1 = new Vector3D(2000.652697, 387.688615, -2666.947760);
405 final Vector3D velR2 = new Vector3D(-3792.46619, -1777.07641, 6856.81495);
406
407
408 final KeplerianOrbit orbit = lambert.estimate(inertialFrame, true, 0,
409 new Position(dateRef, posR1, 1.0, 1.0, satellite),
410 new Position(date2, posR2, 1.0, 1.0, satellite));
411
412
413 Assert.assertEquals(0.0, orbit.getPVCoordinates().getVelocity().getNorm() - velR1.getNorm(), 1e-3);
414
415
416 final KeplerianPropagator kepler = new KeplerianPropagator(orbit);
417 Assert.assertEquals(0.0, kepler.propagate(date2).getPVCoordinates().getVelocity().getNorm() - velR2.getNorm(), 1e-3);
418
419
420 }
421
422 @Before
423 public void setUp() {
424 Utils.setDataRoot("regular-data");
425 }
426
427 }