From 698d8120797383886d1386d3a8ac4b5fc1ca7f16 Mon Sep 17 00:00:00 2001 From: Daniel Karbach Date: Sun, 12 Nov 2017 23:38:09 +0100 Subject: [PATCH] test orbit inverse matrices as well --- src/world/world.cpp | 4 +- tst/world/OrbitTest.cpp | 308 +++++++++++++++++++++++++++++++++++++++- tst/world/OrbitTest.hpp | 22 ++- 3 files changed, 327 insertions(+), 7 deletions(-) diff --git a/src/world/world.cpp b/src/world/world.cpp index c9a7a32..5ba9673 100644 --- 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 ea13860..175da83 100644 --- 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 3f38cce..f9c5238 100644 --- 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(); + }; } -- 2.39.2