1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.orbits;
18
19 import org.hamcrest.MatcherAssert;
20 import org.hipparchus.analysis.UnivariateFunction;
21 import org.hipparchus.analysis.differentiation.DSFactory;
22 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
23 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.linear.MatrixUtils;
26 import org.hipparchus.linear.RealMatrixPreservingVisitor;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathUtils;
29 import org.junit.jupiter.api.AfterEach;
30 import org.junit.jupiter.api.Assertions;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.Utils;
34 import org.orekit.errors.OrekitIllegalArgumentException;
35 import org.orekit.errors.OrekitMessages;
36 import org.orekit.frames.Frame;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.frames.Transform;
39 import org.orekit.time.AbsoluteDate;
40 import org.orekit.time.TimeScalesFactory;
41 import org.orekit.utils.Constants;
42 import org.orekit.utils.PVCoordinates;
43 import org.orekit.utils.TimeStampedPVCoordinates;
44
45 import java.io.ByteArrayInputStream;
46 import java.io.ByteArrayOutputStream;
47 import java.io.IOException;
48 import java.io.ObjectInputStream;
49 import java.io.ObjectOutputStream;
50 import java.util.function.Function;
51
52 import static org.orekit.OrekitMatchers.relativelyCloseTo;
53
54
55 public class EquinoctialOrbitTest {
56
57
58 private AbsoluteDate date;
59
60
61 private double mu;
62
63 @Test
64 void testEquinoctialToEquinoctialEll() {
65
66 double ix = 1.200e-04;
67 double iy = -1.16e-04;
68 double inc = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4.));
69 double hx = FastMath.tan(inc / 2.) * ix / (2 * FastMath.sin(inc / 2.));
70 double hy = FastMath.tan(inc / 2.) * iy / (2 * FastMath.sin(inc / 2.));
71
72
73 EquinoctialOrbit equi =
74 new EquinoctialOrbit(42166712.0, 0.5, -0.5, hx, hy,
75 5.300, PositionAngleType.MEAN,
76 FramesFactory.getEME2000(), date, mu);
77 Vector3D pos = equi.getPosition();
78 Vector3D vit = equi.getPVCoordinates().getVelocity();
79
80 PVCoordinates pvCoordinates = new PVCoordinates(pos, vit);
81
82 EquinoctialOrbit param = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
83 Assertions.assertEquals(param.getA(), equi.getA(), Utils.epsilonTest * equi.getA());
84 Assertions.assertEquals(param.getEquinoctialEx(), equi.getEquinoctialEx(),
85 Utils.epsilonE * FastMath.abs(equi.getE()));
86 Assertions.assertEquals(param.getEquinoctialEy(), equi.getEquinoctialEy(),
87 Utils.epsilonE * FastMath.abs(equi.getE()));
88 Assertions.assertEquals(param.getHx(), equi.getHx(), Utils.epsilonAngle
89 * FastMath.abs(equi.getI()));
90 Assertions.assertEquals(param.getHy(), equi.getHy(), Utils.epsilonAngle
91 * FastMath.abs(equi.getI()));
92 Assertions.assertEquals(MathUtils.normalizeAngle(param.getLv(), equi.getLv()), equi.getLv(),
93 Utils.epsilonAngle * FastMath.abs(equi.getLv()));
94
95 }
96
97 @Test
98 void testEquinoctialToEquinoctialCirc() {
99
100 double ix = 1.200e-04;
101 double iy = -1.16e-04;
102 double inc = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4.));
103 double hx = FastMath.tan(inc / 2.) * ix / (2 * FastMath.sin(inc / 2.));
104 double hy = FastMath.tan(inc / 2.) * iy / (2 * FastMath.sin(inc / 2.));
105
106
107 EquinoctialOrbit equiCir =
108 new EquinoctialOrbit(42166712.0, 0.1e-10, -0.1e-10, hx, hy,
109 5.300, PositionAngleType.MEAN,
110 FramesFactory.getEME2000(), date, mu);
111 Vector3D posCir = equiCir.getPosition();
112 Vector3D vitCir = equiCir.getPVCoordinates().getVelocity();
113
114 PVCoordinates pvCoordinates = new PVCoordinates(posCir, vitCir);
115
116 EquinoctialOrbit paramCir = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(),
117 date, mu);
118 Assertions.assertEquals(paramCir.getA(), equiCir.getA(), Utils.epsilonTest
119 * equiCir.getA());
120 Assertions.assertEquals(paramCir.getEquinoctialEx(), equiCir.getEquinoctialEx(),
121 Utils.epsilonEcir * FastMath.abs(equiCir.getE()));
122 Assertions.assertEquals(paramCir.getEquinoctialEy(), equiCir.getEquinoctialEy(),
123 Utils.epsilonEcir * FastMath.abs(equiCir.getE()));
124 Assertions.assertEquals(paramCir.getHx(), equiCir.getHx(), Utils.epsilonAngle
125 * FastMath.abs(equiCir.getI()));
126 Assertions.assertEquals(paramCir.getHy(), equiCir.getHy(), Utils.epsilonAngle
127 * FastMath.abs(equiCir.getI()));
128 Assertions.assertEquals(MathUtils.normalizeAngle(paramCir.getLv(), equiCir.getLv()), equiCir
129 .getLv(), Utils.epsilonAngle * FastMath.abs(equiCir.getLv()));
130
131 }
132
133 @Test
134 void testEquinoctialToCartesian() {
135
136 double ix = 1.200e-04;
137 double iy = -1.16e-04;
138 double inc = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4.));
139 double hx = FastMath.tan(inc / 2.) * ix / (2 * FastMath.sin(inc / 2.));
140 double hy = FastMath.tan(inc / 2.) * iy / (2 * FastMath.sin(inc / 2.));
141
142 EquinoctialOrbit equi =
143 new EquinoctialOrbit(42166712.0, -7.900e-06, 1.100e-04, hx, hy,
144 5.300, PositionAngleType.MEAN,
145 FramesFactory.getEME2000(), date, mu);
146 Vector3D pos = equi.getPosition();
147 Vector3D vit = equi.getPVCoordinates().getVelocity();
148
149
150 double oneovera = (2. / pos.getNorm()) - vit.getNorm() * vit.getNorm() / mu;
151 Assertions.assertEquals(oneovera, 1. / equi.getA(), 1.0e-7);
152
153 Assertions.assertEquals( 0.233745668678733e+08, pos.getX(), Utils.epsilonTest * FastMath.abs(pos.getX()));
154 Assertions.assertEquals(-0.350998914352669e+08, pos.getY(), Utils.epsilonTest * FastMath.abs(pos.getY()));
155 Assertions.assertEquals(-0.150053723123334e+04, pos.getZ(), Utils.epsilonTest * FastMath.abs(pos.getZ()));
156
157 Assertions.assertEquals(2558.7096558809967, vit.getX(), Utils.epsilonTest * FastMath.abs(vit.getX()));
158 Assertions.assertEquals(1704.1586039092576, vit.getY(), Utils.epsilonTest * FastMath.abs(vit.getY()));
159 Assertions.assertEquals( 0.5013093577879, vit.getZ(), Utils.epsilonTest * FastMath.abs(vit.getZ()));
160
161 }
162
163 @Test
164 void testEquinoctialToKeplerian() {
165
166 double ix = 1.20e-4;
167 double iy = -1.16e-4;
168 double i = 2 * FastMath.asin(FastMath.sqrt((ix * ix + iy * iy) / 4));
169 double hx = FastMath.tan(i / 2) * ix / (2 * FastMath.sin(i / 2));
170 double hy = FastMath.tan(i / 2) * iy / (2 * FastMath.sin(i / 2));
171
172 EquinoctialOrbit equi =
173 new EquinoctialOrbit(42166712.0, -7.900e-6, 1.100e-4, hx, hy,
174 5.300, PositionAngleType.MEAN,
175 FramesFactory.getEME2000(), date, mu);
176 KeplerianOrbit kep = new KeplerianOrbit(equi);
177
178 Assertions.assertEquals(42166712.000, equi.getA(), Utils.epsilonTest * kep.getA());
179 Assertions.assertEquals(0.110283316961361e-03, kep.getE(), Utils.epsilonE
180 * FastMath.abs(kep.getE()));
181 Assertions.assertEquals(0.166901168553917e-03, kep.getI(), Utils.epsilonAngle
182 * FastMath.abs(kep.getI()));
183 Assertions.assertEquals(MathUtils.normalizeAngle(-3.87224326008837, kep.getPerigeeArgument()),
184 kep.getPerigeeArgument(), Utils.epsilonTest
185 * FastMath.abs(kep.getPerigeeArgument()));
186 Assertions.assertEquals(MathUtils.normalizeAngle(5.51473467358854, kep
187 .getRightAscensionOfAscendingNode()), kep
188 .getRightAscensionOfAscendingNode(), Utils.epsilonTest
189 * FastMath.abs(kep.getRightAscensionOfAscendingNode()));
190 Assertions.assertEquals(MathUtils.normalizeAngle(3.65750858649982, kep.getMeanAnomaly()), kep
191 .getMeanAnomaly(), Utils.epsilonTest * FastMath.abs(kep.getMeanAnomaly()));
192
193 }
194
195 @Test
196 void testHyperbolic1() {
197 try {
198 new EquinoctialOrbit(-42166712.0, 1.9, 0.5, 0.01, -0.02, 5.300,
199 PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
200 Assertions.fail("an exception should have been thrown");
201 } catch (OrekitIllegalArgumentException oe) {
202 Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
203 }
204 }
205
206
207 @Test
208 void testHyperbolic2() {
209 Orbit orbit = new KeplerianOrbit(-42166712.0, 1.9, 0.5, 0.01, -0.02, 5.300,
210 PositionAngleType.MEAN, FramesFactory.getEME2000(), date, mu);
211 try {
212 new EquinoctialOrbit(orbit.getPVCoordinates(), orbit.getFrame(), orbit.getMu());
213 Assertions.fail("an exception should have been thrown");
214 } catch (OrekitIllegalArgumentException oe) {
215 Assertions.assertEquals(OrekitMessages.HYPERBOLIC_ORBIT_NOT_HANDLED_AS, oe.getSpecifier());
216 }
217 }
218
219 @Test
220 void testNumericalIssue25() {
221 Vector3D position = new Vector3D(3782116.14107698, 416663.11924914, 5875541.62103057);
222 Vector3D velocity = new Vector3D(-6349.7848910501, 288.4061811651, 4066.9366759691);
223 EquinoctialOrbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
224 FramesFactory.getEME2000(),
225 new AbsoluteDate("2004-01-01T23:00:00.000",
226 TimeScalesFactory.getUTC()),
227 3.986004415E14);
228 Assertions.assertEquals(0.0, orbit.getE(), 2.0e-14);
229 }
230
231
232 @Test
233 void testAnomaly() {
234
235
236 Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
237 Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
238
239 EquinoctialOrbit p = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, mu);
240 KeplerianOrbit kep = new KeplerianOrbit(p);
241
242 double e = p.getE();
243 double eRatio = FastMath.sqrt((1 - e) / (1 + e));
244 double paPraan = kep.getPerigeeArgument()
245 + kep.getRightAscensionOfAscendingNode();
246
247 double lv = 1.1;
248
249 double lE = 2 * FastMath.atan(eRatio * FastMath.tan((lv - paPraan) / 2)) + paPraan;
250 double lM = lE - e * FastMath.sin(lE - paPraan);
251
252 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
253 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
254 p.getFrame(), p.getDate(), p.getMu());
255 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
256 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
257 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
258 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
259 p.getEquinoctialEy() , p.getHx(), p.getHy() , 0 , PositionAngleType.TRUE,
260 p.getFrame(), p.getDate(), p.getMu());
261
262 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
263 p.getEquinoctialEy() , p.getHx(), p.getHy() , lE , PositionAngleType.ECCENTRIC,
264 p.getFrame(), p.getDate(), p.getMu());
265 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
266 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
267 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
268 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
269 p.getEquinoctialEy() , p.getHx(), p.getHy() , 0 , PositionAngleType.TRUE,
270 p.getFrame(), p.getDate(), p.getMu());
271
272 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
273 p.getEquinoctialEy() , p.getHx(), p.getHy() , lM , PositionAngleType.MEAN,
274 p.getFrame(), p.getDate(), p.getMu());
275 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
276 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
277 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
278
279
280 p = new EquinoctialOrbit(p.getA(), 0 ,
281 0, p.getHx(), p.getHy() , p.getLv() , PositionAngleType.TRUE,
282 p.getFrame(), p.getDate(), p.getMu());
283
284 lE = lv;
285 lM = lE;
286
287 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
288 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
289 p.getFrame(), p.getDate(), p.getMu());
290 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
291 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
292 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
293 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
294 p.getEquinoctialEy() , p.getHx(), p.getHy() , 0 , PositionAngleType.TRUE,
295 p.getFrame(), p.getDate(), p.getMu());
296
297 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
298 p.getEquinoctialEy() , p.getHx(), p.getHy() , lE , PositionAngleType.ECCENTRIC,
299 p.getFrame(), p.getDate(), p.getMu());
300 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
301 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
302 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
303 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
304 p.getEquinoctialEy() , p.getHx(), p.getHy() , 0 , PositionAngleType.TRUE,
305 p.getFrame(), p.getDate(), p.getMu());
306
307 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
308 p.getEquinoctialEy() , p.getHx(), p.getHy() , lM , PositionAngleType.MEAN, p.getFrame(), p.getDate(), p.getMu());
309 Assertions.assertEquals(p.getLv(), lv, Utils.epsilonAngle * FastMath.abs(lv));
310 Assertions.assertEquals(p.getLE(), lE, Utils.epsilonAngle * FastMath.abs(lE));
311 Assertions.assertEquals(p.getLM(), lM, Utils.epsilonAngle * FastMath.abs(lM));
312 }
313
314 @Test
315 void testPositionVelocityNorms() {
316
317
318 EquinoctialOrbit p =
319 new EquinoctialOrbit(42166712.0, 0.5, -0.5, 1.200, 2.1,
320 0.67, PositionAngleType.TRUE,
321 FramesFactory.getEME2000(), date, mu);
322
323 double ex = p.getEquinoctialEx();
324 double ey = p.getEquinoctialEy();
325 double lv = p.getLv();
326 double ksi = 1 + ex * FastMath.cos(lv) + ey * FastMath.sin(lv);
327 double nu = ex * FastMath.sin(lv) - ey * FastMath.cos(lv);
328 double epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
329
330 double a = p.getA();
331 double na = FastMath.sqrt(p.getMu() / a);
332
333 Assertions.assertEquals(a * epsilon * epsilon / ksi, p.getPosition().getNorm(),
334 Utils.epsilonTest * FastMath.abs(p.getPosition().getNorm()));
335 Assertions.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon, p
336 .getPVCoordinates().getVelocity().getNorm(), Utils.epsilonTest
337 * FastMath.abs(p.getPVCoordinates().getVelocity().getNorm()));
338
339
340 EquinoctialOrbit pCirEqua =
341 new EquinoctialOrbit(42166712.0, 0.1e-8, 0.1e-8, 0.1e-8, 0.1e-8,
342 0.67, PositionAngleType.TRUE,
343 FramesFactory.getEME2000(), date, mu);
344
345 ex = pCirEqua.getEquinoctialEx();
346 ey = pCirEqua.getEquinoctialEy();
347 lv = pCirEqua.getLv();
348 ksi = 1 + ex * FastMath.cos(lv) + ey * FastMath.sin(lv);
349 nu = ex * FastMath.sin(lv) - ey * FastMath.cos(lv);
350 epsilon = FastMath.sqrt(1 - ex * ex - ey * ey);
351
352 a = pCirEqua.getA();
353 na = FastMath.sqrt(pCirEqua.getMu() / a);
354
355 Assertions.assertEquals(a * epsilon * epsilon / ksi, pCirEqua.getPosition()
356 .getNorm(), Utils.epsilonTest
357 * FastMath.abs(pCirEqua.getPosition().getNorm()));
358 Assertions.assertEquals(na * FastMath.sqrt(ksi * ksi + nu * nu) / epsilon, pCirEqua
359 .getPVCoordinates().getVelocity().getNorm(), Utils.epsilonTest
360 * FastMath.abs(pCirEqua.getPVCoordinates().getVelocity().getNorm()));
361 }
362
363 @Test
364 void testGeometry() {
365
366
367 EquinoctialOrbit p =
368 new EquinoctialOrbit(42166712.0, 0.5, -0.5, 1.200, 2.1,
369 0.67, PositionAngleType.TRUE,
370 FramesFactory.getEME2000(), date, mu);
371
372 Vector3D position = p.getPosition();
373 Vector3D velocity = p.getPVCoordinates().getVelocity();
374 Vector3D momentum = p.getPVCoordinates().getMomentum().normalize();
375
376 double apogeeRadius = p.getA() * (1 + p.getE());
377 double perigeeRadius = p.getA() * (1 - p.getE());
378
379 for (double lv = 0; lv <= 2 * FastMath.PI; lv += 2 * FastMath.PI / 100.) {
380 p = new EquinoctialOrbit(p.getA(), p.getEquinoctialEx(),
381 p.getEquinoctialEy() , p.getHx(), p.getHy() , lv , PositionAngleType.TRUE,
382 p.getFrame(), p.getDate(), p.getMu());
383 position = p.getPosition();
384
385
386
387
388 Assertions.assertTrue((position.getNorm() - apogeeRadius) <= (apogeeRadius * Utils.epsilonTest));
389 Assertions.assertTrue((position.getNorm() - perigeeRadius) >= (-perigeeRadius * Utils.epsilonTest));
390
391 position = position.normalize();
392 velocity = p.getPVCoordinates().getVelocity();
393 velocity = velocity.normalize();
394
395
396
397
398
399 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
400
401 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
402 }
403
404
405 EquinoctialOrbit pCirEqua =
406 new EquinoctialOrbit(42166712.0, 0.1e-8, 0.1e-8, 0.1e-8, 0.1e-8,
407 0.67, PositionAngleType.TRUE,
408 FramesFactory.getEME2000(), date, mu);
409
410 position = pCirEqua.getPosition();
411 velocity = pCirEqua.getPVCoordinates().getVelocity();
412
413 momentum = Vector3D.crossProduct(position, velocity).normalize();
414
415 apogeeRadius = pCirEqua.getA() * (1 + pCirEqua.getE());
416 perigeeRadius = pCirEqua.getA() * (1 - pCirEqua.getE());
417
418 Assertions.assertEquals(perigeeRadius, apogeeRadius, 1.e+4 * Utils.epsilonTest
419 * apogeeRadius);
420
421 for (double lv = 0; lv <= 2 * FastMath.PI; lv += 2 * FastMath.PI / 100.) {
422 pCirEqua = new EquinoctialOrbit(pCirEqua.getA(), pCirEqua.getEquinoctialEx(),
423 pCirEqua.getEquinoctialEy() , pCirEqua.getHx(), pCirEqua.getHy() , lv , PositionAngleType.TRUE,
424 pCirEqua.getFrame(), p.getDate(), p.getMu());
425 position = pCirEqua.getPosition();
426
427
428
429 Assertions.assertTrue((position.getNorm() - apogeeRadius) <= (apogeeRadius * Utils.epsilonTest));
430 Assertions.assertTrue((position.getNorm() - perigeeRadius) >= (-perigeeRadius * Utils.epsilonTest));
431
432 position = position.normalize();
433 velocity = pCirEqua.getPVCoordinates().getVelocity();
434 velocity = velocity.normalize();
435
436
437
438
439
440 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(position, momentum)) < Utils.epsilonTest);
441
442 Assertions.assertTrue(FastMath.abs(Vector3D.dotProduct(velocity, momentum)) < Utils.epsilonTest);
443 }
444 }
445
446
447 @Test
448 void testRadiusOfCurvature() {
449
450
451 EquinoctialOrbit p =
452 new EquinoctialOrbit(42166712.0, 0.5, -0.5, 1.200, 2.1,
453 0.67, PositionAngleType.TRUE,
454 FramesFactory.getEME2000(), date, mu);
455
456
457 Vector3D u = p.getPVCoordinates().getMomentum().orthogonal();
458 Vector3D v = Vector3D.crossProduct(p.getPVCoordinates().getMomentum(), u).normalize();
459
460
461 double xDot = Vector3D.dotProduct(p.getPVCoordinates().getVelocity(), u);
462 double yDot = Vector3D.dotProduct(p.getPVCoordinates().getVelocity(), v);
463 double xDotDot = Vector3D.dotProduct(p.getPVCoordinates().getAcceleration(), u);
464 double yDotDot = Vector3D.dotProduct(p.getPVCoordinates().getAcceleration(), v);
465 double dot2 = xDot * xDot + yDot * yDot;
466 double rCart = dot2 * FastMath.sqrt(dot2) /
467 FastMath.abs(xDot * yDotDot - yDot * xDotDot);
468
469
470 double ex = p.getEquinoctialEx();
471 double ey = p.getEquinoctialEy();
472 double f = ex * FastMath.cos(p.getLE()) + ey * FastMath.sin(p.getLE());
473 double oMf2 = 1 - f * f;
474 double rOrb = p.getA() * oMf2 * FastMath.sqrt(oMf2 / (1 - (ex * ex + ey * ey)));
475
476
477 Assertions.assertEquals(rCart, rOrb, 1.0e-15 * p.getA());
478
479
480
481 Assertions.assertEquals(0.8477 * p.getA(), rCart, 1.0e-4 * p.getA());
482
483 }
484
485 @Test
486 void testSymmetry() {
487
488
489 Vector3D position = new Vector3D(4512.9, 18260., -5127.);
490 Vector3D velocity = new Vector3D(134664.6, 90066.8, 72047.6);
491
492 EquinoctialOrbit p = new EquinoctialOrbit(new PVCoordinates(position, velocity),
493 FramesFactory.getEME2000(), date, mu);
494
495 Vector3D positionOffset = p.getPosition().subtract(position);
496 Vector3D velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
497
498 Assertions.assertEquals(0, positionOffset.getNorm(), 7.5e-12);
499 Assertions.assertEquals(0, velocityOffset.getNorm(), 1.0e-15);
500
501
502 position = new Vector3D(33051.2, 26184.9, -1.3E-5);
503 velocity = new Vector3D(-60376.2, 76208., 2.7E-4);
504
505 p = new EquinoctialOrbit(new PVCoordinates(position, velocity),
506 FramesFactory.getEME2000(), date, mu);
507
508 positionOffset = p.getPosition().subtract(position);
509 velocityOffset = p.getPVCoordinates().getVelocity().subtract(velocity);
510
511 Assertions.assertEquals(0, positionOffset.getNorm(), 1.1e-11);
512 Assertions.assertEquals(0, velocityOffset.getNorm(), 1.0e-15);
513 }
514
515 @Test
516 void testNonInertialFrame() throws IllegalArgumentException {
517 Assertions.assertThrows(IllegalArgumentException.class, () -> {
518 Vector3D position = new Vector3D(4512.9, 18260., -5127.);
519 Vector3D velocity = new Vector3D(134664.6, 90066.8, 72047.6);
520 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
521 new EquinoctialOrbit(pvCoordinates,
522 new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false),
523 date, mu);
524 });
525 }
526
527 @Test
528 void testJacobianReference() {
529
530 AbsoluteDate dateTca = new AbsoluteDate(2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
531 double mu = 3.986004415e+14;
532 EquinoctialOrbit orbEqu = new EquinoctialOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
533 FastMath.toRadians(40.), PositionAngleType.MEAN,
534 FramesFactory.getEME2000(), dateTca, mu);
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605 Vector3D pRef = new Vector3D(2004367.298657628707588, 6575317.978060320019722, -1518024.843913963763043);
606 Vector3D vRef = new Vector3D(5574.048661495634406, -368.839015744295409, 5009.529487849066754);
607 double[][] jRef = {
608 { 0.56305379787310628, 1.8470954710993663, -0.42643364527246025, 1370.4369387322224, -90.682848736736688 , 1231.6441195141242 },
609 { 9.52434720041122055E-008, 9.49704503778007296E-008, 4.46607520107935678E-008, 1.69704446323098610E-004, 7.05603505855828105E-005, 1.14825140460141970E-004 },
610 { -5.41784097802642701E-008, 9.54903765833015538E-008, -8.95815777332234450E-008, 1.01864980963344096E-004, -1.03194262242761416E-004, 1.40668700715197768E-004 },
611 { 1.96680305426455816E-007, -1.12388745957974467E-007, -2.27118924123407353E-007, 2.06472886488132167E-004, -1.17984506564646906E-004, -2.38427023682723818E-004 },
612 { -2.24382495052235118E-007, 1.28218568601277626E-007, 2.59108357381747656E-007, 1.89034327703662092E-004, -1.08019615830663994E-004, -2.18289640324466583E-004 },
613 { -3.04001022071876804E-007, 1.22214683774559989E-007, 1.35141804810132761E-007, -1.34034616931480536E-004, -2.14283975204169379E-004, 1.29018773893081404E-004 }
614 };
615
616 PVCoordinates pv = orbEqu.getPVCoordinates();
617 Assertions.assertEquals(0, pv.getPosition().subtract(pRef).getNorm(), 2.0e-16 * pRef.getNorm());
618 Assertions.assertEquals(0, pv.getVelocity().subtract(vRef).getNorm(), 2.0e-16 * vRef.getNorm());
619
620 double[][] jacobian = new double[6][6];
621 orbEqu.getJacobianWrtCartesian(PositionAngleType.MEAN, jacobian);
622
623 for (int i = 0; i < jacobian.length; i++) {
624 double[] row = jacobian[i];
625 double[] rowRef = jRef[i];
626 for (int j = 0; j < row.length; j++) {
627 Assertions.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 4.0e-15);
628 }
629 }
630
631 }
632
633 @Test
634 void testJacobianFinitedifferences() {
635
636 AbsoluteDate dateTca = new AbsoluteDate(2000, 4, 1, 0, 0, 0.000, TimeScalesFactory.getUTC());
637 double mu = 3.986004415e+14;
638 EquinoctialOrbit orbEqu = new EquinoctialOrbit(7000000.0, 0.01, -0.02, 1.2, 2.1,
639 FastMath.toRadians(40.), PositionAngleType.MEAN,
640 FramesFactory.getEME2000(), dateTca, mu);
641
642 for (PositionAngleType type : PositionAngleType.values()) {
643 double hP = 2.0;
644 double[][] finiteDiffJacobian = finiteDifferencesJacobian(type, orbEqu, hP);
645 double[][] jacobian = new double[6][6];
646 orbEqu.getJacobianWrtCartesian(type, jacobian);
647
648 for (int i = 0; i < jacobian.length; i++) {
649 double[] row = jacobian[i];
650 double[] rowRef = finiteDiffJacobian[i];
651 for (int j = 0; j < row.length; j++) {
652 Assertions.assertEquals(0, (row[j] - rowRef[j]) / rowRef[j], 4.0e-9);
653 }
654 }
655
656 double[][] invJacobian = new double[6][6];
657 orbEqu.getJacobianWrtParameters(type, invJacobian);
658 MatrixUtils.createRealMatrix(jacobian).
659 multiply(MatrixUtils.createRealMatrix(invJacobian)).
660 walkInRowOrder(new RealMatrixPreservingVisitor() {
661 public void start(int rows, int columns,
662 int startRow, int endRow, int startColumn, int endColumn) {
663 }
664
665 public void visit(int row, int column, double value) {
666 Assertions.assertEquals(row == column ? 1.0 : 0.0, value, 7.0e-9);
667 }
668
669 public double end() {
670 return Double.NaN;
671 }
672 });
673
674 }
675
676 }
677
678 private double[][] finiteDifferencesJacobian(PositionAngleType type, EquinoctialOrbit orbit, double hP)
679 {
680 double[][] jacobian = new double[6][6];
681 for (int i = 0; i < 6; ++i) {
682 fillColumn(type, i, orbit, hP, jacobian);
683 }
684 return jacobian;
685 }
686
687 private void fillColumn(PositionAngleType type, int i, EquinoctialOrbit orbit, double hP, double[][] jacobian) {
688
689
690
691 Vector3D p = orbit.getPosition();
692 Vector3D v = orbit.getPVCoordinates().getVelocity();
693 double hV = orbit.getMu() * hP / (v.getNorm() * p.getNormSq());
694
695 double h;
696 Vector3D dP = Vector3D.ZERO;
697 Vector3D dV = Vector3D.ZERO;
698 switch (i) {
699 case 0:
700 h = hP;
701 dP = new Vector3D(hP, 0, 0);
702 break;
703 case 1:
704 h = hP;
705 dP = new Vector3D(0, hP, 0);
706 break;
707 case 2:
708 h = hP;
709 dP = new Vector3D(0, 0, hP);
710 break;
711 case 3:
712 h = hV;
713 dV = new Vector3D(hV, 0, 0);
714 break;
715 case 4:
716 h = hV;
717 dV = new Vector3D(0, hV, 0);
718 break;
719 default:
720 h = hV;
721 dV = new Vector3D(0, 0, hV);
722 break;
723 }
724
725 EquinoctialOrbit oM4h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, -4, dP), new Vector3D(1, v, -4, dV)),
726 orbit.getFrame(), orbit.getDate(), orbit.getMu());
727 EquinoctialOrbit oM3h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, -3, dP), new Vector3D(1, v, -3, dV)),
728 orbit.getFrame(), orbit.getDate(), orbit.getMu());
729 EquinoctialOrbit oM2h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, -2, dP), new Vector3D(1, v, -2, dV)),
730 orbit.getFrame(), orbit.getDate(), orbit.getMu());
731 EquinoctialOrbit oM1h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, -1, dP), new Vector3D(1, v, -1, dV)),
732 orbit.getFrame(), orbit.getDate(), orbit.getMu());
733 EquinoctialOrbit oP1h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, +1, dP), new Vector3D(1, v, +1, dV)),
734 orbit.getFrame(), orbit.getDate(), orbit.getMu());
735 EquinoctialOrbit oP2h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, +2, dP), new Vector3D(1, v, +2, dV)),
736 orbit.getFrame(), orbit.getDate(), orbit.getMu());
737 EquinoctialOrbit oP3h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, +3, dP), new Vector3D(1, v, +3, dV)),
738 orbit.getFrame(), orbit.getDate(), orbit.getMu());
739 EquinoctialOrbit oP4h = new EquinoctialOrbit(new PVCoordinates(new Vector3D(1, p, +4, dP), new Vector3D(1, v, +4, dV)),
740 orbit.getFrame(), orbit.getDate(), orbit.getMu());
741
742 jacobian[0][i] = (-3 * (oP4h.getA() - oM4h.getA()) +
743 32 * (oP3h.getA() - oM3h.getA()) -
744 168 * (oP2h.getA() - oM2h.getA()) +
745 672 * (oP1h.getA() - oM1h.getA())) / (840 * h);
746 jacobian[1][i] = (-3 * (oP4h.getEquinoctialEx() - oM4h.getEquinoctialEx()) +
747 32 * (oP3h.getEquinoctialEx() - oM3h.getEquinoctialEx()) -
748 168 * (oP2h.getEquinoctialEx() - oM2h.getEquinoctialEx()) +
749 672 * (oP1h.getEquinoctialEx() - oM1h.getEquinoctialEx())) / (840 * h);
750 jacobian[2][i] = (-3 * (oP4h.getEquinoctialEy() - oM4h.getEquinoctialEy()) +
751 32 * (oP3h.getEquinoctialEy() - oM3h.getEquinoctialEy()) -
752 168 * (oP2h.getEquinoctialEy() - oM2h.getEquinoctialEy()) +
753 672 * (oP1h.getEquinoctialEy() - oM1h.getEquinoctialEy())) / (840 * h);
754 jacobian[3][i] = (-3 * (oP4h.getHx() - oM4h.getHx()) +
755 32 * (oP3h.getHx() - oM3h.getHx()) -
756 168 * (oP2h.getHx() - oM2h.getHx()) +
757 672 * (oP1h.getHx() - oM1h.getHx())) / (840 * h);
758 jacobian[4][i] = (-3 * (oP4h.getHy() - oM4h.getHy()) +
759 32 * (oP3h.getHy() - oM3h.getHy()) -
760 168 * (oP2h.getHy() - oM2h.getHy()) +
761 672 * (oP1h.getHy() - oM1h.getHy())) / (840 * h);
762 jacobian[5][i] = (-3 * (oP4h.getL(type) - oM4h.getL(type)) +
763 32 * (oP3h.getL(type) - oM3h.getL(type)) -
764 168 * (oP2h.getL(type) - oM2h.getL(type)) +
765 672 * (oP1h.getL(type) - oM1h.getL(type))) / (840 * h);
766
767 }
768
769 @Test
770 void testSerialization()
771 throws IOException, ClassNotFoundException {
772 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
773 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
774 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity);
775 EquinoctialOrbit orbit = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
776 Assertions.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
777
778 ByteArrayOutputStream bos = new ByteArrayOutputStream();
779 ObjectOutputStream oos = new ObjectOutputStream(bos);
780 oos.writeObject(orbit);
781
782 Assertions.assertEquals(bos.size(), 458);
783
784 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
785 ObjectInputStream ois = new ObjectInputStream(bis);
786 EquinoctialOrbit deserialized = (EquinoctialOrbit) ois.readObject();
787 Assertions.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
788 Assertions.assertEquals(orbit.getEquinoctialEx(), deserialized.getEquinoctialEx(), 1.0e-10);
789 Assertions.assertEquals(orbit.getEquinoctialEy(), deserialized.getEquinoctialEy(), 1.0e-10);
790 Assertions.assertEquals(orbit.getHx(), deserialized.getHx(), 1.0e-10);
791 Assertions.assertEquals(orbit.getHy(), deserialized.getHy(), 1.0e-10);
792 Assertions.assertEquals(orbit.getLv(), deserialized.getLv(), 1.0e-10);
793 Assertions.assertTrue(Double.isNaN(orbit.getADot()) && Double.isNaN(deserialized.getADot()));
794 Assertions.assertTrue(Double.isNaN(orbit.getEquinoctialExDot()) && Double.isNaN(deserialized.getEquinoctialExDot()));
795 Assertions.assertTrue(Double.isNaN(orbit.getEquinoctialEyDot()) && Double.isNaN(deserialized.getEquinoctialEyDot()));
796 Assertions.assertTrue(Double.isNaN(orbit.getHxDot()) && Double.isNaN(deserialized.getHxDot()));
797 Assertions.assertTrue(Double.isNaN(orbit.getHyDot()) && Double.isNaN(deserialized.getHyDot()));
798 Assertions.assertTrue(Double.isNaN(orbit.getLvDot()) && Double.isNaN(deserialized.getLvDot()));
799 Assertions.assertEquals(orbit.getDate(), deserialized.getDate());
800 Assertions.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
801 Assertions.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
802
803 }
804
805 @Test
806 void testSerializationWithDerivatives()
807 throws IOException, ClassNotFoundException {
808 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
809 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
810 double r2 = position.getNormSq();
811 double r = FastMath.sqrt(r2);
812 Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
813 1, new Vector3D(-0.1, 0.2, 0.3));
814 PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
815 EquinoctialOrbit orbit = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
816 Assertions.assertEquals(42255170.003, orbit.getA(), 1.0e-3);
817
818 ByteArrayOutputStream bos = new ByteArrayOutputStream();
819 ObjectOutputStream oos = new ObjectOutputStream(bos);
820 oos.writeObject(orbit);
821
822 Assertions.assertEquals(bos.size(), 506);
823
824 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
825 ObjectInputStream ois = new ObjectInputStream(bis);
826 EquinoctialOrbit deserialized = (EquinoctialOrbit) ois.readObject();
827 Assertions.assertEquals(orbit.getA(), deserialized.getA(), 1.0e-10);
828 Assertions.assertEquals(orbit.getEquinoctialEx(), deserialized.getEquinoctialEx(), 1.0e-10);
829 Assertions.assertEquals(orbit.getEquinoctialEy(), deserialized.getEquinoctialEy(), 1.0e-10);
830 Assertions.assertEquals(orbit.getHx(), deserialized.getHx(), 1.0e-10);
831 Assertions.assertEquals(orbit.getHy(), deserialized.getHy(), 1.0e-10);
832 Assertions.assertEquals(orbit.getLv(), deserialized.getLv(), 1.0e-10);
833 Assertions.assertEquals(orbit.getADot(), deserialized.getADot(), 1.0e-10);
834 Assertions.assertEquals(orbit.getEquinoctialExDot(), deserialized.getEquinoctialExDot(), 1.0e-10);
835 Assertions.assertEquals(orbit.getEquinoctialEyDot(), deserialized.getEquinoctialEyDot(), 1.0e-10);
836 Assertions.assertEquals(orbit.getHxDot(), deserialized.getHxDot(), 1.0e-10);
837 Assertions.assertEquals(orbit.getHyDot(), deserialized.getHyDot(), 1.0e-10);
838 Assertions.assertEquals(orbit.getLvDot(), deserialized.getLvDot(), 1.0e-10);
839 Assertions.assertEquals(orbit.getDate(), deserialized.getDate());
840 Assertions.assertEquals(orbit.getMu(), deserialized.getMu(), 1.0e-10);
841 Assertions.assertEquals(orbit.getFrame().getName(), deserialized.getFrame().getName());
842
843 }
844
845 @Test
846 void testNonKeplerianDerivatives() {
847 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
848 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
849 final Vector3D velocity = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
850 final Vector3D acceleration = new Vector3D(-7.466182457944, -2.118153357345, 0.160004048437);
851 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
852 final Frame frame = FramesFactory.getEME2000();
853 final double mu = Constants.EIGEN5C_EARTH_MU;
854 final EquinoctialOrbit orbit = new EquinoctialOrbit(pv, frame, mu);
855
856 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getA()),
857 orbit.getADot(),
858 4.3e-8);
859 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEx()),
860 orbit.getEquinoctialExDot(),
861 2.1e-15);
862 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getEquinoctialEy()),
863 orbit.getEquinoctialEyDot(),
864 5.3e-16);
865 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHx()),
866 orbit.getHxDot(),
867 4.4e-15);
868 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getHy()),
869 orbit.getHyDot(),
870 1.5e-15);
871 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLv()),
872 orbit.getLvDot(),
873 1.2e-15);
874 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLE()),
875 orbit.getLEDot(),
876 7.7e-16);
877 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getLM()),
878 orbit.getLMDot(),
879 8.8e-16);
880 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getE()),
881 orbit.getEDot(),
882 6.9e-16);
883 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getI()),
884 orbit.getIDot(),
885 3.5e-15);
886 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.TRUE)),
887 orbit.getLDot(PositionAngleType.TRUE),
888 1.2e-15);
889 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.ECCENTRIC)),
890 orbit.getLDot(PositionAngleType.ECCENTRIC),
891 7.7e-16);
892 Assertions.assertEquals(differentiate(pv, frame, mu, shifted -> shifted.getL(PositionAngleType.MEAN)),
893 orbit.getLDot(PositionAngleType.MEAN),
894 8.8e-16);
895
896 }
897
898 private <S extends Function<EquinoctialOrbit, Double>>
899 double differentiate(TimeStampedPVCoordinates pv, Frame frame, double mu, S picker) {
900 final DSFactory factory = new DSFactory(1, 1);
901 FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 0.1);
902 UnivariateDifferentiableFunction diff = differentiator.differentiate((UnivariateFunction) dt ->
903 picker.apply(new EquinoctialOrbit(pv.shiftedBy(dt), frame, mu)));
904 return diff.value(factory.variable(0, 0.0)).getPartialDerivative(1);
905 }
906
907 @Test
908 void testPositionAngleDerivatives() {
909 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
910 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
911 final Vector3D velocity = new Vector3D(166.816407662, -1106.783301861, -7372.745712770);
912 final Vector3D acceleration = new Vector3D(-7.466182457944, -2.118153357345, 0.160004048437);
913 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity, acceleration);
914 final Frame frame = FramesFactory.getEME2000();
915 final double mu = Constants.EIGEN5C_EARTH_MU;
916 final EquinoctialOrbit orbit = new EquinoctialOrbit(pv, frame, mu);
917
918 for (PositionAngleType type : PositionAngleType.values()) {
919 final EquinoctialOrbit rebuilt = new EquinoctialOrbit(orbit.getA(),
920 orbit.getEquinoctialEx(),
921 orbit.getEquinoctialEy(),
922 orbit.getHx(),
923 orbit.getHy(),
924 orbit.getL(type),
925 orbit.getADot(),
926 orbit.getEquinoctialExDot(),
927 orbit.getEquinoctialEyDot(),
928 orbit.getHxDot(),
929 orbit.getHyDot(),
930 orbit.getLDot(type),
931 type, orbit.getFrame(), orbit.getDate(), orbit.getMu());
932 MatcherAssert.assertThat(rebuilt.getA(), relativelyCloseTo(orbit.getA(), 1));
933 MatcherAssert.assertThat(rebuilt.getEquinoctialEx(), relativelyCloseTo(orbit.getEquinoctialEx(), 1));
934 MatcherAssert.assertThat(rebuilt.getEquinoctialEy(), relativelyCloseTo(orbit.getEquinoctialEy(), 1));
935 MatcherAssert.assertThat(rebuilt.getHx(), relativelyCloseTo(orbit.getHx(), 1));
936 MatcherAssert.assertThat(rebuilt.getHy(), relativelyCloseTo(orbit.getHy(), 1));
937 MatcherAssert.assertThat(rebuilt.getADot(), relativelyCloseTo(orbit.getADot(), 1));
938 MatcherAssert.assertThat(rebuilt.getEquinoctialExDot(), relativelyCloseTo(orbit.getEquinoctialExDot(), 1));
939 MatcherAssert.assertThat(rebuilt.getEquinoctialEyDot(), relativelyCloseTo(orbit.getEquinoctialEyDot(), 1));
940 MatcherAssert.assertThat(rebuilt.getHxDot(), relativelyCloseTo(orbit.getHxDot(), 1));
941 MatcherAssert.assertThat(rebuilt.getHyDot(), relativelyCloseTo(orbit.getHyDot(), 1));
942 for (PositionAngleType type2 : PositionAngleType.values()) {
943 MatcherAssert.assertThat(rebuilt.getL(type2), relativelyCloseTo(orbit.getL(type2), 1));
944 MatcherAssert.assertThat(rebuilt.getLDot(type2), relativelyCloseTo(orbit.getLDot(type2), 1));
945 }
946 }
947
948 }
949
950 @Test
951 void testEquatorialRetrograde() {
952 Vector3D position = new Vector3D(10000000.0, 0.0, 0.0);
953 Vector3D velocity = new Vector3D(0.0, -6500.0, 0.0);
954 double r2 = position.getNormSq();
955 double r = FastMath.sqrt(r2);
956 Vector3D acceleration = new Vector3D(-mu / (r * r2), position,
957 1, new Vector3D(-0.1, 0.2, 0.3));
958 PVCoordinates pvCoordinates = new PVCoordinates(position, velocity, acceleration);
959
960
961
962 EquinoctialOrbit orbit = new EquinoctialOrbit(new KeplerianOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu));
963 Assertions.assertEquals(10637829.465, orbit.getA(), 1.0e-3);
964 Assertions.assertEquals(-738.145, orbit.getADot(), 1.0e-3);
965 Assertions.assertEquals(0.05995861, orbit.getE(), 1.0e-8);
966 Assertions.assertEquals(-6.523e-5, orbit.getEDot(), 1.0e-8);
967 Assertions.assertTrue(Double.isNaN(orbit.getI()));
968 Assertions.assertTrue(Double.isNaN(orbit.getIDot()));
969 Assertions.assertTrue(Double.isNaN(orbit.getHx()));
970 Assertions.assertTrue(Double.isNaN(orbit.getHxDot()));
971 Assertions.assertTrue(Double.isNaN(orbit.getHy()));
972 Assertions.assertTrue(Double.isNaN(orbit.getHyDot()));
973 }
974
975 @Test
976 void testDerivativesConversionSymmetry() {
977 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:01:20.000", TimeScalesFactory.getUTC());
978 Vector3D position = new Vector3D(6893443.400234382, 1886406.1073757345, -589265.1150359757);
979 Vector3D velocity = new Vector3D(-281.1261461082365, -1231.6165642450928, -7348.756363469432);
980 Vector3D acceleration = new Vector3D(-7.460341170581685, -2.0415957334584527, 0.6393322823627762);
981 PVCoordinates pvCoordinates = new PVCoordinates( position, velocity, acceleration);
982 EquinoctialOrbit orbit = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(),
983 date, Constants.EIGEN5C_EARTH_MU);
984 Assertions.assertTrue(orbit.hasDerivatives());
985 double r2 = position.getNormSq();
986 double r = FastMath.sqrt(r2);
987 Vector3D keplerianAcceleration = new Vector3D(-orbit.getMu() / (r2 * r), position);
988 Assertions.assertEquals(0.0101, Vector3D.distance(keplerianAcceleration, acceleration), 1.0e-4);
989
990 for (OrbitType type : OrbitType.values()) {
991 Orbit converted = type.convertType(orbit);
992 Assertions.assertTrue(converted.hasDerivatives());
993 EquinoctialOrbit rebuilt = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(converted);
994 Assertions.assertTrue(rebuilt.hasDerivatives());
995 Assertions.assertEquals(orbit.getADot(), rebuilt.getADot(), 3.0e-13);
996 Assertions.assertEquals(orbit.getEquinoctialExDot(), rebuilt.getEquinoctialExDot(), 1.0e-15);
997 Assertions.assertEquals(orbit.getEquinoctialEyDot(), rebuilt.getEquinoctialEyDot(), 1.0e-15);
998 Assertions.assertEquals(orbit.getHxDot(), rebuilt.getHxDot(), 1.0e-15);
999 Assertions.assertEquals(orbit.getHyDot(), rebuilt.getHyDot(), 1.0e-15);
1000 Assertions.assertEquals(orbit.getLvDot(), rebuilt.getLvDot(), 1.0e-15);
1001 }
1002
1003 }
1004
1005 @Test
1006 void testToString() {
1007 Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1008 Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1009 PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1010 EquinoctialOrbit orbit = new EquinoctialOrbit(pvCoordinates, FramesFactory.getEME2000(), date, mu);
1011 Assertions.assertEquals("equinoctial parameters: {a: 4.225517000282565E7; ex: 5.927324978565528E-4; ey: -0.002062743969643666; hx: 6.401103130239252E-5; hy: -0.0017606836670756732; lv: 134.24111947709974;}",
1012 orbit.toString());
1013 }
1014
1015 @Test
1016 void testRemoveRates() {
1017
1018 final Vector3D position = new Vector3D(-29536113.0, 30329259.0, -100125.0);
1019 final Vector3D velocity = new Vector3D(-2194.0, -2141.0, -8.0);
1020 final PVCoordinates pvCoordinates = new PVCoordinates(position, velocity);
1021 final EquinoctialOrbit orbit = new EquinoctialOrbit(pvCoordinates, FramesFactory.getGCRF(), date, mu);
1022
1023 final EquinoctialOrbit orbitWithoutRates = orbit.removeRates();
1024
1025 Assertions.assertFalse(orbitWithoutRates.hasRates());
1026 Assertions.assertTrue(Double.isNaN(orbitWithoutRates.getADot()));
1027 Assertions.assertEquals(orbit.getMu(), orbitWithoutRates.getMu());
1028 Assertions.assertEquals(orbit.getDate(), orbitWithoutRates.getDate());
1029 Assertions.assertEquals(orbit.getFrame(), orbitWithoutRates.getFrame());
1030 Assertions.assertEquals(orbit.getA(), orbitWithoutRates.getA());
1031 Assertions.assertEquals(orbit.getEquinoctialEx(), orbitWithoutRates.getEquinoctialEx());
1032 Assertions.assertEquals(orbit.getEquinoctialEy(), orbitWithoutRates.getEquinoctialEy());
1033 Assertions.assertEquals(orbit.getHx(), orbitWithoutRates.getHx());
1034 Assertions.assertEquals(orbit.getHy(), orbitWithoutRates.getHy());
1035 Assertions.assertEquals(orbit.getLv(), orbitWithoutRates.getLv());
1036 }
1037
1038 @Test
1039 void testCopyNonKeplerianAcceleration() {
1040
1041 final Frame eme2000 = FramesFactory.getEME2000();
1042
1043
1044 final Vector3D position = new Vector3D(42164140, 0, 0);
1045
1046 final PVCoordinates pv = new PVCoordinates(position,
1047 new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
1048
1049 final Orbit orbit = new EquinoctialOrbit(pv, eme2000, date, mu);
1050
1051
1052 final Orbit orbitCopy = new EquinoctialOrbit(orbit);
1053
1054
1055 final Orbit shiftedOrbit = orbit.shiftedBy(10);
1056 final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10);
1057
1058 Assertions.assertEquals(0.0,
1059 Vector3D.distance(shiftedOrbit.getPosition(),
1060 shiftedOrbitCopy.getPosition()),
1061 1.0e-10);
1062 Assertions.assertEquals(0.0,
1063 Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
1064 shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1065 1.0e-10);
1066
1067 }
1068
1069 @Test
1070 void testNormalize() {
1071 EquinoctialOrbit withoutDerivatives =
1072 new EquinoctialOrbit(42166712.0, 0.005, -0.025, 0.17, 0.34,
1073 0.4, PositionAngleType.MEAN,
1074 FramesFactory.getEME2000(), date, mu);
1075 EquinoctialOrbit ref =
1076 new EquinoctialOrbit(24000000.0, -0.012, 0.01, 0.2, 0.1,
1077 -6.28, PositionAngleType.MEAN,
1078 FramesFactory.getEME2000(), date, mu);
1079
1080 EquinoctialOrbit normalized1 = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.normalize(withoutDerivatives, ref);
1081 Assertions.assertFalse(normalized1.hasDerivatives());
1082 Assertions.assertEquals(0.0, normalized1.getA() - withoutDerivatives.getA(), 1.0e-6);
1083 Assertions.assertEquals(0.0, normalized1.getEquinoctialEx() - withoutDerivatives.getEquinoctialEx(), 1.0e-10);
1084 Assertions.assertEquals(0.0, normalized1.getEquinoctialEy() - withoutDerivatives.getEquinoctialEy(), 1.0e-10);
1085 Assertions.assertEquals(0.0, normalized1.getHx() - withoutDerivatives.getHx(), 1.0e-10);
1086 Assertions.assertEquals(0.0, normalized1.getHy() - withoutDerivatives.getHy(), 1.0e-10);
1087 Assertions.assertEquals(-MathUtils.TWO_PI, normalized1.getLv() - withoutDerivatives.getLv(), 1.0e-10);
1088 Assertions.assertTrue(Double.isNaN(normalized1.getADot()));
1089 Assertions.assertTrue(Double.isNaN(normalized1.getEquinoctialExDot()));
1090 Assertions.assertTrue(Double.isNaN(normalized1.getEquinoctialEyDot()));
1091 Assertions.assertTrue(Double.isNaN(normalized1.getHxDot()));
1092 Assertions.assertTrue(Double.isNaN(normalized1.getHyDot()));
1093 Assertions.assertTrue(Double.isNaN(normalized1.getLvDot()));
1094
1095 double[] p = new double[6];
1096 OrbitType.EQUINOCTIAL.mapOrbitToArray(withoutDerivatives, PositionAngleType.TRUE, p, null);
1097 EquinoctialOrbit withDerivatives = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.mapArrayToOrbit(p,
1098 new double[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 },
1099 PositionAngleType.TRUE,
1100 withoutDerivatives.getDate(),
1101 withoutDerivatives.getMu(),
1102 withoutDerivatives.getFrame());
1103 EquinoctialOrbit normalized2 = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.normalize(withDerivatives, ref);
1104 Assertions.assertTrue(normalized2.hasDerivatives());
1105 Assertions.assertFalse(normalized1.hasDerivatives());
1106 Assertions.assertEquals(0.0, normalized1.getA() - withoutDerivatives.getA(), 1.0e-6);
1107 Assertions.assertEquals(0.0, normalized1.getEquinoctialEx() - withoutDerivatives.getEquinoctialEx(), 1.0e-10);
1108 Assertions.assertEquals(0.0, normalized1.getEquinoctialEy() - withoutDerivatives.getEquinoctialEy(), 1.0e-10);
1109 Assertions.assertEquals(0.0, normalized1.getHx() - withoutDerivatives.getHx(), 1.0e-10);
1110 Assertions.assertEquals(0.0, normalized1.getHy() - withoutDerivatives.getHy(), 1.0e-10);
1111 Assertions.assertEquals(-MathUtils.TWO_PI, normalized1.getLv() - withoutDerivatives.getLv(), 1.0e-10);
1112 Assertions.assertEquals(0.0, normalized2.getADot() - withDerivatives.getADot(), 1.0e-10);
1113 Assertions.assertEquals(0.0, normalized2.getEquinoctialExDot() - withDerivatives.getEquinoctialExDot(), 1.0e-10);
1114 Assertions.assertEquals(0.0, normalized2.getEquinoctialEyDot() - withDerivatives.getEquinoctialEyDot(), 1.0e-10);
1115 Assertions.assertEquals(0.0, normalized2.getHxDot() - withDerivatives.getHxDot(), 1.0e-10);
1116 Assertions.assertEquals(0.0, normalized2.getHyDot() - withDerivatives.getHyDot(), 1.0e-10);
1117 Assertions.assertEquals(0.0, normalized2.getLvDot() - withDerivatives.getLvDot(), 1.0e-10);
1118
1119 }
1120
1121 @Test
1122 @Deprecated
1123 void positionAngleNonRegressionOnDeprecated() {
1124
1125
1126 final double ex = 0.2;
1127 final double ey = 0.3;
1128 final double originalPositionAngle = 1.;
1129
1130 final double actualEccentricToMean = EquinoctialOrbit.eccentricToMean(originalPositionAngle, ex, ey);
1131 final double actualEccentricToTrue = EquinoctialOrbit.eccentricToTrue(originalPositionAngle, ex, ey);
1132 final double actualMeanToEccentric = EquinoctialOrbit.meanToEccentric(originalPositionAngle, ex, ey);
1133 final double actualTrueToEccentric = EquinoctialOrbit.trueToEccentric(originalPositionAngle, ex, ey);
1134
1135 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.eccentricToMean(ex, ey, originalPositionAngle),
1136 actualEccentricToMean);
1137 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.eccentricToTrue(ex, ey, originalPositionAngle),
1138 actualEccentricToTrue);
1139 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.meanToEccentric(ex, ey, originalPositionAngle),
1140 actualMeanToEccentric);
1141 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.trueToEccentric(ex, ey, originalPositionAngle),
1142 actualTrueToEccentric);
1143 }
1144
1145 @Test
1146 void testCoverageCachedPositionAngleTypeWithRates() {
1147
1148 final double semiMajorAxis = 1e4;
1149 final double ex = 0.;
1150 final double ey = 0.;
1151 final double expectedL = 0.;
1152 final double expectedLDot = 0.;
1153
1154 for (final PositionAngleType inputPositionAngleType: PositionAngleType.values()) {
1155 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1156 final EquinoctialOrbit equinoctialOrbit = new EquinoctialOrbit(semiMajorAxis, ex, ey, 0., 0.,
1157 expectedL, 0., 0., 0., 0., 0., expectedLDot,
1158 inputPositionAngleType, cachedPositionAngleType, FramesFactory.getGCRF(), date, mu);
1159 Assertions.assertEquals(expectedL, equinoctialOrbit.getLv());
1160 Assertions.assertEquals(expectedL, equinoctialOrbit.getLM());
1161 Assertions.assertEquals(expectedL, equinoctialOrbit.getLE());
1162 Assertions.assertEquals(expectedLDot, equinoctialOrbit.getLvDot());
1163 Assertions.assertEquals(expectedLDot, equinoctialOrbit.getLMDot());
1164 Assertions.assertEquals(expectedLDot, equinoctialOrbit.getLEDot());
1165 }
1166 }
1167 }
1168
1169 @Test
1170 void testCachedPositionAngleTypeTrue() {
1171
1172 final double semiMajorAxis = 1e4;
1173 final double ex = 1e-2;
1174 final double ey = -1e-3;
1175 final double expectedLv = 2.;
1176 final PositionAngleType inputPositionAngleType = PositionAngleType.TRUE;
1177
1178 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1179 final EquinoctialOrbit equinoctialOrbit = new EquinoctialOrbit(semiMajorAxis, ex, ey, 0., 0.,
1180 expectedLv, inputPositionAngleType, cachedPositionAngleType, FramesFactory.getGCRF(), date, mu);
1181 Assertions.assertEquals(expectedLv, equinoctialOrbit.getLv(), 1e-15);
1182 final double actualL = equinoctialOrbit.getL(cachedPositionAngleType);
1183 switch (cachedPositionAngleType) {
1184
1185 case MEAN:
1186 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.trueToMean(ex, ey, expectedLv),
1187 actualL);
1188 break;
1189
1190 case ECCENTRIC:
1191 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.trueToEccentric(ex, ey, expectedLv),
1192 actualL);
1193 break;
1194
1195 case TRUE:
1196 Assertions.assertEquals(expectedLv, actualL);
1197 break;
1198 }
1199 }
1200 }
1201
1202 @Test
1203 void testCachedPositionAngleTypeMean() {
1204
1205 final double semiMajorAxis = 1e4;
1206 final double ex = 1e-2;
1207 final double ey = -1e-3;
1208 final double expectedLM = 2.;
1209 final PositionAngleType inputPositionAngleType = PositionAngleType.MEAN;
1210
1211 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1212 final EquinoctialOrbit equinoctialOrbit = new EquinoctialOrbit(semiMajorAxis, ex, ey, 0., 0.,
1213 expectedLM, inputPositionAngleType, cachedPositionAngleType, FramesFactory.getGCRF(), date, mu);
1214 Assertions.assertEquals(expectedLM, equinoctialOrbit.getLM(), 1e-15);
1215 final double actualL = equinoctialOrbit.getL(cachedPositionAngleType);
1216 switch (cachedPositionAngleType) {
1217
1218 case TRUE:
1219 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.meanToTrue(ex, ey, expectedLM),
1220 actualL);
1221 break;
1222
1223 case ECCENTRIC:
1224 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.meanToEccentric(ex, ey, expectedLM),
1225 actualL);
1226 break;
1227
1228 case MEAN:
1229 Assertions.assertEquals(expectedLM, actualL);
1230 break;
1231 }
1232 }
1233 }
1234
1235 @Test
1236 void testCachedPositionAngleTypeEccentric() {
1237
1238 final double semiMajorAxis = 1e4;
1239 final double ex = 1e-2;
1240 final double ey = -1e-3;
1241 final double expectedLE = 2.;
1242 final PositionAngleType inputPositionAngleType = PositionAngleType.ECCENTRIC;
1243
1244 for (final PositionAngleType cachedPositionAngleType: PositionAngleType.values()) {
1245 final EquinoctialOrbit equinoctialOrbit = new EquinoctialOrbit(semiMajorAxis, ex, ey, 0., 0.,
1246 expectedLE, inputPositionAngleType, cachedPositionAngleType, FramesFactory.getGCRF(), date, mu);
1247 Assertions.assertEquals(expectedLE, equinoctialOrbit.getLE(), 1e-15);
1248 final double actualL = equinoctialOrbit.getL(cachedPositionAngleType);
1249 switch (cachedPositionAngleType) {
1250
1251 case TRUE:
1252 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.eccentricToTrue(ex, ey, expectedLE),
1253 actualL);
1254 break;
1255
1256 case MEAN:
1257 Assertions.assertEquals(EquinoctialLongitudeArgumentUtility.eccentricToMean(ex, ey, expectedLE),
1258 actualL);
1259 break;
1260
1261 case ECCENTRIC:
1262 Assertions.assertEquals(expectedLE, actualL);
1263 break;
1264 }
1265 }
1266 }
1267
1268 @BeforeEach
1269 public void setUp() {
1270
1271 Utils.setDataRoot("regular-data");
1272
1273
1274 date = AbsoluteDate.J2000_EPOCH;
1275
1276
1277 mu = 3.9860047e14;
1278 }
1279
1280 @AfterEach
1281 public void tearDown() {
1282 date = null;
1283 }
1284
1285 }