1 #include "OrbitTest.hpp"
3 #include "../assert.hpp"
5 #include "math/const.hpp"
6 #include "world/Orbit.hpp"
8 CPPUNIT_TEST_SUITE_REGISTRATION(blobs::world::test::OrbitTest);
10 using blobs::test::AssertEqual;
17 void OrbitTest::setUp() {
20 void OrbitTest::tearDown() {
24 void OrbitTest::testDefault() {
25 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
28 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
29 "wrong semi-major axis on default orbit",
30 1.0, orbit.SemiMajorAxis(), epsilon
32 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
33 "wrong eccentricity on default orbit",
34 0.0, orbit.Eccentricity(), epsilon
36 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
37 "wrong inclination on default orbit",
38 0.0, orbit.Inclination(), epsilon
40 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
41 "wrong longitude of ascending node on default orbit",
42 0.0, orbit.LongitudeAscending(), epsilon
44 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
45 "wrong argument of periapsis on default orbit",
46 0.0, orbit.ArgumentPeriapsis(), epsilon
48 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
49 "wrong mean anomaly on default orbit",
50 0.0, orbit.MeanAnomaly(), epsilon
53 // reference direction is +X, so at t=0, the body should be
54 // at (sma,0,0) relative to its parent
55 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
57 "wrong position at t=0",
58 glm::dvec3(1.0, 0.0, 0.0),
59 glm::dvec3(pos) / pos.w,
63 // at 90° position should be (0,0,sma) since the zero inclination
64 // reference plane is XZ and rotates counter-clockwise
65 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
67 "wrong position at t=90°",
68 glm::dvec3(0.0, 0.0, -1.0),
69 glm::dvec3(pos) / pos.w,
73 // at 180° position should be (-sma,0,0)
74 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
76 "wrong position at t=180°",
77 glm::dvec3(-1.0, 0.0, 0.0),
78 glm::dvec3(pos) / pos.w,
82 // at 270° position should be (0,0,-sma)
83 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
85 "wrong position at t=270°",
86 glm::dvec3(0.0, 0.0, 1.0),
87 glm::dvec3(pos) / pos.w,
91 // at 360° position should be (sma,0,0), the initial position
92 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
94 "wrong position at t=360°",
95 glm::dvec3(1.0, 0.0, 0.0),
96 glm::dvec3(pos) / pos.w,
101 void OrbitTest::testSMA() {
102 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
105 orbit.SemiMajorAxis(2.0);
106 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
107 "wrong semi-major axis on orbit",
108 2.0, orbit.SemiMajorAxis(), epsilon
111 // reference direction is +X, so at t=0, the body should be
112 // at (sma,0,0) relative to its parent
113 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
115 "wrong position at t=0",
116 glm::dvec3(2.0, 0.0, 0.0),
117 glm::dvec3(pos) / pos.w,
121 // at 90° position should be (0,0,sma) since the zero inclination
122 // reference plane is XZ and rotates counter-clockwise
123 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
125 "wrong position at t=90°",
126 glm::dvec3(0.0, 0.0, -2.0),
127 glm::dvec3(pos) / pos.w,
131 // at 180° position should be (-sma,0,0)
132 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
134 "wrong position at t=180°",
135 glm::dvec3(-2.0, 0.0, 0.0),
136 glm::dvec3(pos) / pos.w,
140 // at 270° position should be (0,0,-sma)
141 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
143 "wrong position at t=270°",
144 glm::dvec3(0.0, 0.0, 2.0),
145 glm::dvec3(pos) / pos.w,
149 // at 360° position should be (sma,0,0), the initial position
150 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
152 "wrong position at t=360°",
153 glm::dvec3(2.0, 0.0, 0.0),
154 glm::dvec3(pos) / pos.w,
159 void OrbitTest::testEcc() {
160 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
163 orbit.Eccentricity(0.5);
164 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
165 "wrong eccentricity on orbit",
166 0.5, orbit.Eccentricity(), epsilon
169 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
171 "wrong position at t=0",
172 glm::dvec3(0.5, 0.0, 0.0),
173 glm::dvec3(pos) / pos.w,
177 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
179 "wrong position at t=90°",
180 glm::dvec3(-0.93513085903671, 0.0, -0.779740887497559),
181 glm::dvec3(pos) / pos.w,
185 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
187 "wrong position at t=180°",
188 glm::dvec3(-1.5, 0.0, 0.0),
189 glm::dvec3(pos) / pos.w,
193 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
195 "wrong position at t=270°",
196 glm::dvec3(-0.93513085903671, 0.0, 0.779740887497559),
197 glm::dvec3(pos) / pos.w,
201 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
203 "wrong position at t=360°",
204 glm::dvec3(0.5, 0.0, 0.0),
205 glm::dvec3(pos) / pos.w,
210 void OrbitTest::testInc() {
211 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
214 orbit.Inclination(PI * 0.25); // 45°
215 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
216 "wrong inclination on orbit",
217 PI * 0.25, orbit.Inclination(), epsilon
220 // inclination rotates counter clockwise around +X, so at t=0 should be
221 // at (sma,0,0) relative to its parent
222 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
224 "wrong position at t=0",
225 glm::dvec3(1.0, 0.0, 0.0),
226 glm::dvec3(pos) / pos.w,
230 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
232 "wrong position at t=90°",
233 glm::dvec3(0.0, 0.707106781186548, -0.707106781186548),
234 glm::dvec3(pos) / pos.w,
238 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
240 "wrong position at t=180°",
241 glm::dvec3(-1.0, 0.0, 0.0),
242 glm::dvec3(pos) / pos.w,
246 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
248 "wrong position at t=270°",
249 glm::dvec3(0.0, -0.707106781186548, 0.707106781186548),
250 glm::dvec3(pos) / pos.w,
254 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
256 "wrong position at t=360°",
257 glm::dvec3(1.0, 0.0, 0.0),
258 glm::dvec3(pos) / pos.w,
263 void OrbitTest::testLngAsc() {
264 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
267 orbit.LongitudeAscending(PI * 0.25); // 45°
268 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
269 "wrong longitude of ascending node on orbit",
270 PI * 0.25, orbit.LongitudeAscending(), epsilon
272 // using an inclination of 90° as well to make the rotation more apparent
273 orbit.Inclination(PI * 0.5);
275 // inclination rotates counter clockwise around +X, while LAN rotates it
276 // around +Y, so at t=0 should be at (sma*sin(45°),0,-sma*cos(45°))
277 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
279 "wrong position at t=0",
280 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
281 glm::dvec3(pos) / pos.w,
285 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
287 "wrong position at t=90°",
288 glm::dvec3(0.0, 1.0, 0.0),
289 glm::dvec3(pos) / pos.w,
293 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
295 "wrong position at t=180°",
296 glm::dvec3(-0.707106781186548, 0.0, 0.707106781186548),
297 glm::dvec3(pos) / pos.w,
301 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
303 "wrong position at t=270°",
304 glm::dvec3(0.0, -1.0, 0.0),
305 glm::dvec3(pos) / pos.w,
309 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
311 "wrong position at t=360°",
312 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
313 glm::dvec3(pos) / pos.w,
318 void OrbitTest::testArgPe() {
319 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
322 orbit.ArgumentPeriapsis(PI * 0.25); // 45°
323 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
324 "wrong argument of periapsis node on orbit",
325 PI * 0.25, orbit.ArgumentPeriapsis(), epsilon
327 // using an inclination of 90° as well to make the rotation more apparent
328 orbit.Inclination(PI * 0.5);
330 // inclination rotates counter clockwise around +X, while APe rotates it
331 // around +Y in the rotated coordinate system, so at t=0 should be at
332 // (sma*sin(45°),0,sma*cos(45°))
333 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
335 "wrong position at t=0",
336 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
337 glm::dvec3(pos) / pos.w,
341 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
343 "wrong position at t=90°",
344 glm::dvec3(0.0, 1.0, 0.0),
345 glm::dvec3(pos) / pos.w,
349 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
351 "wrong position at t=180°",
352 glm::dvec3(-0.707106781186548, 0.0, -0.707106781186548),
353 glm::dvec3(pos) / pos.w,
357 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
359 "wrong position at t=270°",
360 glm::dvec3(0.0, -1.0, 0.0),
361 glm::dvec3(pos) / pos.w,
365 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
367 "wrong position at t=360°",
368 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
369 glm::dvec3(pos) / pos.w,
374 void OrbitTest::testMnAn() {
375 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
378 orbit.MeanAnomaly(PI * 0.25); // 45°
379 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
380 "wrong mean anomaly on default orbit",
381 PI * 0.25, orbit.MeanAnomaly(), epsilon
384 // mean anomaly just phase shifts the orbit
385 glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
387 "wrong position at t=0",
388 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
389 glm::dvec3(pos) / pos.w,
393 // at 90° position should be (0,0,sma) since the zero inclination
394 // reference plane is XZ and rotates counter-clockwise
395 pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
397 "wrong position at t=90°",
398 glm::dvec3(-0.707106781186548, 0.0, -0.707106781186548),
399 glm::dvec3(pos) / pos.w,
403 // at 180° position should be (-sma,0,0)
404 pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
406 "wrong position at t=180°",
407 glm::dvec3(-0.707106781186548, 0.0, 0.707106781186548),
408 glm::dvec3(pos) / pos.w,
412 // at 270° position should be (0,0,-sma)
413 pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
415 "wrong position at t=270°",
416 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
417 glm::dvec3(pos) / pos.w,
421 // at 360° position should be (sma,0,0), the initial position
422 pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
424 "wrong position at t=360°",
425 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
426 glm::dvec3(pos) / pos.w,
431 void OrbitTest::testInverseDefault() {
432 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
436 // inverse matrix should project expected orbit position back to the origin
437 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(1.0, 0.0, 0.0, 1.0));
439 "wrong position at t=0",
440 glm::dvec3(0.0, 0.0, 0.0),
441 glm::dvec3(pos) / pos.w,
445 // at 90° position should be (0,0,sma) since the zero inclination
446 // reference plane is XZ and rotates counter-clockwise
447 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 0.0, -1.0, 1.0);
449 "wrong position at t=90°",
450 glm::dvec3(0.0, 0.0, 0.0),
451 glm::dvec3(pos) / pos.w,
455 // at 180° position should be (-sma,0,0)
456 pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.0, 0.0, 0.0, 1.0);
458 "wrong position at t=180°",
459 glm::dvec3(0.0, 0.0, 0.0),
460 glm::dvec3(pos) / pos.w,
464 // at 270° position should be (0,0,-sma)
465 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 1.0, 1.0);
467 "wrong position at t=270°",
468 glm::dvec3(0.0, 0.0, 0.0),
469 glm::dvec3(pos) / pos.w,
473 // at 360° position should be (sma,0,0), the initial position
474 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(1.0, 0.0, 0.0, 1.0);
476 "wrong position at t=360°",
477 glm::dvec3(0.0, 0.0, 0.0),
478 glm::dvec3(pos) / pos.w,
483 void OrbitTest::testInverseSMA() {
484 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
487 orbit.SemiMajorAxis(2.0);
489 // reference direction is +X, so at t=0, the body should be
490 // at (sma,0,0) relative to its parent
491 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(2.0, 0.0, 0.0, 1.0));
493 "wrong position at t=0",
494 glm::dvec3(0.0, 0.0, 0.0),
495 glm::dvec3(pos) / pos.w,
499 // at 90° position should be (0,0,sma) since the zero inclination
500 // reference plane is XZ and rotates counter-clockwise
501 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 0.0, -2.0, 1.0);
503 "wrong position at t=90°",
504 glm::dvec3(0.0, 0.0, 0.0),
505 glm::dvec3(pos) / pos.w,
509 // at 180° position should be (-sma,0,0)
510 pos = orbit.InverseMatrix(PI) * glm::dvec4(-2.0, 0.0, 0.0, 1.0);
512 "wrong position at t=180°",
513 glm::dvec3(0.0, 0.0, 0.0),
514 glm::dvec3(pos) / pos.w,
518 // at 270° position should be (0,0,-sma)
519 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 2.0, 1.0);
521 "wrong position at t=270°",
522 glm::dvec3(0.0, 0.0, 0.0),
523 glm::dvec3(pos) / pos.w,
527 // at 360° position should be (sma,0,0), the initial position
528 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(2.0, 0.0, 0.0, 1.0);
530 "wrong position at t=360°",
531 glm::dvec3(0.0, 0.0, 0.0),
532 glm::dvec3(pos) / pos.w,
537 void OrbitTest::testInverseEcc() {
538 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
541 orbit.Eccentricity(0.5);
543 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(0.5, 0.0, 0.0, 1.0));
545 "wrong position at t=0",
546 glm::dvec3(0.0, 0.0, 0.0),
547 glm::dvec3(pos) / pos.w,
551 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(-0.93513085903671, 0.0, -0.779740887497559, 1.0);
553 "wrong position at t=90°",
554 glm::dvec3(0.0, 0.0, 0.0),
555 glm::dvec3(pos) / pos.w,
559 pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.5, 0.0, 0.0, 1.0);
561 "wrong position at t=180°",
562 glm::dvec3(0.0, 0.0, 0.0),
563 glm::dvec3(pos) / pos.w,
567 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(-0.93513085903671, 0.0, 0.779740887497559, 1.0);
569 "wrong position at t=270°",
570 glm::dvec3(0.0, 0.0, 0.0),
571 glm::dvec3(pos) / pos.w,
575 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.5, 0.0, 0.0, 1.0);
577 "wrong position at t=360°",
578 glm::dvec3(0.0, 0.0, 0.0),
579 glm::dvec3(pos) / pos.w,
584 void OrbitTest::testInverseInc() {
585 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
588 orbit.Inclination(PI * 0.25); // 45°
590 // inclination rotates counter clockwise around +X, so at t=0 should be
591 // at (sma,0,0) relative to its parent
592 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(1.0, 0.0, 0.0, 1.0));
594 "wrong position at t=0",
595 glm::dvec3(0.0, 0.0, 0.0),
596 glm::dvec3(pos) / pos.w,
600 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 0.707106781186548, -0.707106781186548, 1.0);
602 "wrong position at t=90°",
603 glm::dvec3(0.0, 0.0, 0.0),
604 glm::dvec3(pos) / pos.w,
608 pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.0, 0.0, 0.0, 1.0);
610 "wrong position at t=180°",
611 glm::dvec3(0.0, 0.0, 0.0),
612 glm::dvec3(pos) / pos.w,
616 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -0.707106781186548, 0.707106781186548, 1.0);
618 "wrong position at t=270°",
619 glm::dvec3(0.0, 0.0, 0.0),
620 glm::dvec3(pos) / pos.w,
624 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(1.0, 0.0, 0.0, 1.0);
626 "wrong position at t=360°",
627 glm::dvec3(0.0, 0.0, 0.0),
628 glm::dvec3(pos) / pos.w,
633 void OrbitTest::testInverseLngAsc() {
634 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
637 orbit.LongitudeAscending(PI * 0.25); // 45°
638 orbit.Inclination(PI * 0.5);
640 // inclination rotates counter clockwise around +X, while LAN rotates it
641 // around +Y, so at t=0 should be at (sma*sin(45°),0,-sma*cos(45°))
642 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(0.707106781186548, 0.0, -0.707106781186548, 1.0));
644 "wrong position at t=0",
645 glm::dvec3(0.0, 0.0, 0.0),
646 glm::dvec3(pos) / pos.w,
650 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 1.0, 0.0, 1.0);
652 "wrong position at t=90°",
653 glm::dvec3(0.0, 0.0, 0.0),
654 glm::dvec3(pos) / pos.w,
658 pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, 0.707106781186548, 1.0);
660 "wrong position at t=180°",
661 glm::dvec3(0.0, 0.0, 0.0),
662 glm::dvec3(pos) / pos.w,
666 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -1.0, 0.0, 1.0);
668 "wrong position at t=270°",
669 glm::dvec3(0.0, 0.0, 0.0),
670 glm::dvec3(pos) / pos.w,
674 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.707106781186548, 0.0, -0.707106781186548, 1.0);
676 "wrong position at t=360°",
677 glm::dvec3(0.0, 0.0, 0.0),
678 glm::dvec3(pos) / pos.w,
683 void OrbitTest::testInverseArgPe() {
684 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
687 orbit.ArgumentPeriapsis(PI * 0.25); // 45°
688 orbit.Inclination(PI * 0.5);
690 // inclination rotates counter clockwise around +X, while APe rotates it
691 // around +Y in the rotated coordinate system, so at t=0 should be at
692 // (sma*sin(45°),0,sma*cos(45°))
693 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(0.707106781186548, 0.0, 0.707106781186548, 1.0));
695 "wrong position at t=0",
696 glm::dvec3(0.0, 0.0, 0.0),
697 glm::dvec3(pos) / pos.w,
701 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 1.0, 0.0, 1.0);
703 "wrong position at t=90°",
704 glm::dvec3(0.0, 0.0, 0.0),
705 glm::dvec3(pos) / pos.w,
709 pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, -0.707106781186548, 1.0);
711 "wrong position at t=180°",
712 glm::dvec3(0.0, 0.0, 0.0),
713 glm::dvec3(pos) / pos.w,
717 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -1.0, 0.0, 1.0);
719 "wrong position at t=270°",
720 glm::dvec3(0.0, 0.0, 0.0),
721 glm::dvec3(pos) / pos.w,
725 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.707106781186548, 0.0, 0.707106781186548, 1.0);
727 "wrong position at t=360°",
728 glm::dvec3(0.0, 0.0, 0.0),
729 glm::dvec3(pos) / pos.w,
734 void OrbitTest::testInverseMnAn() {
735 constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
738 orbit.MeanAnomaly(PI * 0.25); // 45°
740 // mean anomaly just phase shifts the orbit
741 glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(0.707106781186548, 0.0, -0.707106781186548, 1.0));
743 "wrong position at t=0",
744 glm::dvec3(0.0, 0.0, 0.0),
745 glm::dvec3(pos) / pos.w,
749 // at 90° position should be (0,0,sma) since the zero inclination
750 // reference plane is XZ and rotates counter-clockwise
751 pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(-0.707106781186548, 0.0, -0.707106781186548, 1.0);
753 "wrong position at t=90°",
754 glm::dvec3(0.0, 0.0, 0.0),
755 glm::dvec3(pos) / pos.w,
759 // at 180° position should be (-sma,0,0)
760 pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, 0.707106781186548, 1.0);
762 "wrong position at t=180°",
763 glm::dvec3(0.0, 0.0, 0.0),
764 glm::dvec3(pos) / pos.w,
768 // at 270° position should be (0,0,-sma)
769 pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.707106781186548, 0.0, 0.707106781186548, 1.0);
771 "wrong position at t=270°",
772 glm::dvec3(0.0, 0.0, 0.0),
773 glm::dvec3(pos) / pos.w,
777 // at 360° position should be (sma,0,0), the initial position
778 pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.707106781186548, 0.0, -0.707106781186548, 1.0);
780 "wrong position at t=360°",
781 glm::dvec3(0.0, 0.0, 0.0),
782 glm::dvec3(pos) / pos.w,