test orbit inverse matrices as well

diff --git a/src/world/world.cpp b/src/world/world.cpp
index c9a7a326f3d94a2f4bc496e4100da541df82e794..5ba96732596c8657b176db6e6faa8a0c811c51b3 100644 (file)
--- a/src/world/world.cpp
+++ b/src/world/world.cpp
@@ -235,7 +235,7 @@ glm::dmat4 Orbit::Matrix(double t) const noexcept {
        double P = sma * (cos(E) - ecc);
        double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
 
-       return glm::translate(glm::yawPitchRoll(asc, inc, arg), glm::dvec3(P, 0.0, -Q));
+       return glm::yawPitchRoll(asc, inc, arg) * glm::translate(glm::dvec3(P, 0.0, -Q));
 }
 
 glm::dmat4 Orbit::InverseMatrix(double t) const noexcept {
@@ -243,7 +243,7 @@ glm::dmat4 Orbit::InverseMatrix(double t) const noexcept {
        double E = mean2eccentric(M, ecc);
        double P = sma * (cos(E) - ecc);
        double Q = sma * sin(E) * sqrt(1 - (ecc * ecc));
-       return glm::transpose(glm::yawPitchRoll(asc, inc, arg)) * glm::translate(glm::dvec3(-P, 0.0, Q));
+       return glm::translate(glm::dvec3(-P, 0.0, Q)) * glm::transpose(glm::yawPitchRoll(asc, inc, arg));
 }
 
 

diff --git a/tst/world/OrbitTest.cpp b/tst/world/OrbitTest.cpp
index ea13860ea3d0d6c4d904a009124d95d3b19a1e7a..175da83db81e751fa30d5c7c1f3eb41b471b0792 100644 (file)
--- a/tst/world/OrbitTest.cpp
+++ b/tst/world/OrbitTest.cpp
@@ -21,7 +21,7 @@ void OrbitTest::tearDown() {
 }
 
 
-void OrbitTest::testDefaultOrbit() {
+void OrbitTest::testDefault() {
        Orbit orbit;
        CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
                "wrong semi-major axis on default orbit",
@@ -379,6 +379,312 @@ void OrbitTest::testMnAn() {
        );
 }
 
+void OrbitTest::testInverseDefault() {
+       Orbit orbit;
+
+       // inverse matrix should project expected orbit position back to the origin
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(1.0f, 0.0f, 0.0f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 90° position should be (0,0,sma) since the zero inclination
+       // reference plane is XZ and rotates counter-clockwise
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(0.0f, 0.0f, -1.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 180° position should be (-sma,0,0)
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-1.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 270° position should be (0,0,-sma)
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.0f, 0.0f, 1.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 360° position should be (sma,0,0), the initial position
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(1.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseSMA() {
+       Orbit orbit;
+       orbit.SemiMajorAxis(2.0);
+
+       // reference direction is +X, so at t=0, the body should be
+       // at (sma,0,0) relative to its parent
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(2.0f, 0.0f, 0.0f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 90° position should be (0,0,sma) since the zero inclination
+       // reference plane is XZ and rotates counter-clockwise
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(0.0f, 0.0f, -2.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 180° position should be (-sma,0,0)
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-2.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 270° position should be (0,0,-sma)
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.0f, 0.0f, 2.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 360° position should be (sma,0,0), the initial position
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(2.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseEcc() {
+       Orbit orbit;
+       orbit.Eccentricity(0.5);
+
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(0.5f, 0.0f, 0.0f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(-0.935130834579468f, 0.0f, -0.779740869998932f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-1.5f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(-0.935130834579468f, 0.0f, 0.779740869998932f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(0.5f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseInc() {
+       Orbit orbit;
+       orbit.Inclination(PI * 0.25); // 45°
+
+       // inclination rotates counter clockwise around +X, so at t=0 should be
+       // at (sma,0,0) relative to its parent
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(1.0f, 0.0f, 0.0f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(0.0f, 0.70710676908493f, -0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-1.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.0f, -0.70710676908493f, 0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(1.0f, 0.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseLngAsc() {
+       Orbit orbit;
+       orbit.LongitudeAscending(PI * 0.25); // 45°
+       orbit.Inclination(PI * 0.5);
+
+       // inclination rotates counter clockwise around +X, while LAN rotates it
+       // around +Y, so at t=0 should be at (sma*sin(45°),0,-sma*cos(45°))
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(0.0f, 1.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-0.70710676908493f, 0.0f, 0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.0f, -1.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseArgPe() {
+       Orbit orbit;
+       orbit.ArgumentPeriapsis(PI * 0.25); // 45°
+       orbit.Inclination(PI * 0.5);
+
+       // inclination rotates counter clockwise around +X, while APe rotates it
+       // around +Y in the rotated coordinate system, so at t=0 should be at
+       // (sma*sin(45°),0,sma*cos(45°))
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(0.70710676908493f, 0.0f, 0.70710676908493f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(0.0f, 1.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.0f, -1.0f, 0.0f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(0.70710676908493f, 0.0f, 0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
+
+void OrbitTest::testInverseMnAn() {
+       Orbit orbit;
+       orbit.MeanAnomaly(PI * 0.25); // 45°
+
+       // mean anomaly just phase shifts the orbit
+       glm::vec4 pos(orbit.InverseMatrix(0.0) * glm::vec4(0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f));
+       AssertEqual(
+               "wrong position at t=0",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 90° position should be (0,0,sma) since the zero inclination
+       // reference plane is XZ and rotates counter-clockwise
+       pos = orbit.InverseMatrix(PI_0p5) * glm::vec4(-0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=90°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 180° position should be (-sma,0,0)
+       pos = orbit.InverseMatrix(PI) * glm::vec4(-0.70710676908493f, 0.0f, 0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=180°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 270° position should be (0,0,-sma)
+       pos = orbit.InverseMatrix(PI_1p5) * glm::vec4(0.70710676908493f, 0.0f, 0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=270°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+
+       // at 360° position should be (sma,0,0), the initial position
+       pos = orbit.InverseMatrix(PI_2p0) * glm::vec4(0.70710676908493f, 0.0f, -0.70710676908493f, 1.0f);
+       AssertEqual(
+               "wrong position at t=360°",
+               glm::vec3(0.0f, 0.0f, 0.0f),
+               glm::vec3(pos) / pos.w
+       );
+}
 }
 }
 }

diff --git a/tst/world/OrbitTest.hpp b/tst/world/OrbitTest.hpp
index 3f38cce7ba6f5210db64188c7ba074ebb750812b..f9c5238a19dfba6af81fb7230632ede0daa39026 100644 (file)
--- a/tst/world/OrbitTest.hpp
+++ b/tst/world/OrbitTest.hpp
@@ -13,8 +13,7 @@ class OrbitTest
 
 CPPUNIT_TEST_SUITE(OrbitTest);
 
-CPPUNIT_TEST(testDefaultOrbit);
-
+CPPUNIT_TEST(testDefault);
 CPPUNIT_TEST(testSMA);
 CPPUNIT_TEST(testEcc);
 CPPUNIT_TEST(testInc);
@@ -22,14 +21,21 @@ CPPUNIT_TEST(testLngAsc);
 CPPUNIT_TEST(testArgPe);
 CPPUNIT_TEST(testMnAn);
 
+CPPUNIT_TEST(testInverseDefault);
+CPPUNIT_TEST(testInverseSMA);
+CPPUNIT_TEST(testInverseEcc);
+CPPUNIT_TEST(testInverseInc);
+CPPUNIT_TEST(testInverseLngAsc);
+CPPUNIT_TEST(testInverseArgPe);
+CPPUNIT_TEST(testInverseMnAn);
+
 CPPUNIT_TEST_SUITE_END();
 
 public:
        void setUp();
        void tearDown();
 
-       void testDefaultOrbit();
-
+       void testDefault();
        void testSMA();
        void testEcc();
        void testInc();
@@ -37,6 +43,14 @@ public:
        void testArgPe();
        void testMnAn();
 
+       void testInverseDefault();
+       void testInverseSMA();
+       void testInverseEcc();
+       void testInverseInc();
+       void testInverseLngAsc();
+       void testInverseArgPe();
+       void testInverseMnAn();
+
 };
 
 }