tst/world/OrbitTest.cpp

   1 #include "OrbitTest.hpp"
   2
   3 #include "../assert.hpp"
   4
   5 #include "math/const.hpp"
   6 #include "world/Orbit.hpp"
   7
   8 CPPUNIT_TEST_SUITE_REGISTRATION(blobs::world::test::OrbitTest);
   9
  10 using blobs::test::AssertEqual;
  11
  12
  13 namespace blobs {
  14 namespace world {
  15 namespace test {
  16
  17 void OrbitTest::setUp() {
  18 }
  19
  20 void OrbitTest::tearDown() {
  21 }
  22
  23
  24 void OrbitTest::testDefault() {
  25         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
  26
  27         Orbit orbit;
  28         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  29                 "wrong semi-major axis on default orbit",
  30                 1.0, orbit.SemiMajorAxis(), epsilon
  31         );
  32         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  33                 "wrong eccentricity on default orbit",
  34                 0.0, orbit.Eccentricity(), epsilon
  35         );
  36         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  37                 "wrong inclination on default orbit",
  38                 0.0, orbit.Inclination(), epsilon
  39         );
  40         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  41                 "wrong longitude of ascending node on default orbit",
  42                 0.0, orbit.LongitudeAscending(), epsilon
  43         );
  44         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  45                 "wrong argument of periapsis on default orbit",
  46                 0.0, orbit.ArgumentPeriapsis(), epsilon
  47         );
  48         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
  49                 "wrong mean anomaly on default orbit",
  50                 0.0, orbit.MeanAnomaly(), epsilon
  51         );
  52
  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));
  56         AssertEqual(
  57                 "wrong position at t=0",
  58                 glm::dvec3(1.0, 0.0, 0.0),
  59                 glm::dvec3(pos) / pos.w,
  60                 epsilon
  61         );
  62
  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);
  66         AssertEqual(
  67                 "wrong position at t=90°",
  68                 glm::dvec3(0.0, 0.0, -1.0),
  69                 glm::dvec3(pos) / pos.w,
  70                 epsilon
  71         );
  72
  73         // at 180° position should be (-sma,0,0)
  74         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
  75         AssertEqual(
  76                 "wrong position at t=180°",
  77                 glm::dvec3(-1.0, 0.0, 0.0),
  78                 glm::dvec3(pos) / pos.w,
  79                 epsilon
  80         );
  81
  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);
  84         AssertEqual(
  85                 "wrong position at t=270°",
  86                 glm::dvec3(0.0, 0.0, 1.0),
  87                 glm::dvec3(pos) / pos.w,
  88                 epsilon
  89         );
  90
  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);
  93         AssertEqual(
  94                 "wrong position at t=360°",
  95                 glm::dvec3(1.0, 0.0, 0.0),
  96                 glm::dvec3(pos) / pos.w,
  97                 epsilon
  98         );
  99 }
 100
 101 void OrbitTest::testSMA() {
 102         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 103
 104         Orbit orbit;
 105         orbit.SemiMajorAxis(2.0);
 106         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
 107                 "wrong semi-major axis on orbit",
 108                 2.0, orbit.SemiMajorAxis(), epsilon
 109         );
 110
 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));
 114         AssertEqual(
 115                 "wrong position at t=0",
 116                 glm::dvec3(2.0, 0.0, 0.0),
 117                 glm::dvec3(pos) / pos.w,
 118                 epsilon
 119         );
 120
 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);
 124         AssertEqual(
 125                 "wrong position at t=90°",
 126                 glm::dvec3(0.0, 0.0, -2.0),
 127                 glm::dvec3(pos) / pos.w,
 128                 epsilon
 129         );
 130
 131         // at 180° position should be (-sma,0,0)
 132         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 133         AssertEqual(
 134                 "wrong position at t=180°",
 135                 glm::dvec3(-2.0, 0.0, 0.0),
 136                 glm::dvec3(pos) / pos.w,
 137                 epsilon
 138         );
 139
 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);
 142         AssertEqual(
 143                 "wrong position at t=270°",
 144                 glm::dvec3(0.0, 0.0, 2.0),
 145                 glm::dvec3(pos) / pos.w,
 146                 epsilon
 147         );
 148
 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);
 151         AssertEqual(
 152                 "wrong position at t=360°",
 153                 glm::dvec3(2.0, 0.0, 0.0),
 154                 glm::dvec3(pos) / pos.w,
 155                 epsilon
 156         );
 157 }
 158
 159 void OrbitTest::testEcc() {
 160         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 161
 162         Orbit orbit;
 163         orbit.Eccentricity(0.5);
 164         CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(
 165                 "wrong eccentricity on orbit",
 166                 0.5, orbit.Eccentricity(), epsilon
 167         );
 168
 169         glm::dvec4 pos(orbit.Matrix(0.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0));
 170         AssertEqual(
 171                 "wrong position at t=0",
 172                 glm::dvec3(0.5, 0.0, 0.0),
 173                 glm::dvec3(pos) / pos.w,
 174                 epsilon
 175         );
 176
 177         pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 178         AssertEqual(
 179                 "wrong position at t=90°",
 180                 glm::dvec3(-0.93513085903671, 0.0, -0.779740887497559),
 181                 glm::dvec3(pos) / pos.w,
 182                 epsilon
 183         );
 184
 185         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 186         AssertEqual(
 187                 "wrong position at t=180°",
 188                 glm::dvec3(-1.5, 0.0, 0.0),
 189                 glm::dvec3(pos) / pos.w,
 190                 epsilon
 191         );
 192
 193         pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 194         AssertEqual(
 195                 "wrong position at t=270°",
 196                 glm::dvec3(-0.93513085903671, 0.0, 0.779740887497559),
 197                 glm::dvec3(pos) / pos.w,
 198                 epsilon
 199         );
 200
 201         pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 202         AssertEqual(
 203                 "wrong position at t=360°",
 204                 glm::dvec3(0.5, 0.0, 0.0),
 205                 glm::dvec3(pos) / pos.w,
 206                 epsilon
 207         );
 208 }
 209
 210 void OrbitTest::testInc() {
 211         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 212
 213         Orbit orbit;
 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
 218         );
 219
 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));
 223         AssertEqual(
 224                 "wrong position at t=0",
 225                 glm::dvec3(1.0, 0.0, 0.0),
 226                 glm::dvec3(pos) / pos.w,
 227                 epsilon
 228         );
 229
 230         pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 231         AssertEqual(
 232                 "wrong position at t=90°",
 233                 glm::dvec3(0.0, 0.707106781186548, -0.707106781186548),
 234                 glm::dvec3(pos) / pos.w,
 235                 epsilon
 236         );
 237
 238         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 239         AssertEqual(
 240                 "wrong position at t=180°",
 241                 glm::dvec3(-1.0, 0.0, 0.0),
 242                 glm::dvec3(pos) / pos.w,
 243                 epsilon
 244         );
 245
 246         pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 247         AssertEqual(
 248                 "wrong position at t=270°",
 249                 glm::dvec3(0.0, -0.707106781186548, 0.707106781186548),
 250                 glm::dvec3(pos) / pos.w,
 251                 epsilon
 252         );
 253
 254         pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 255         AssertEqual(
 256                 "wrong position at t=360°",
 257                 glm::dvec3(1.0, 0.0, 0.0),
 258                 glm::dvec3(pos) / pos.w,
 259                 epsilon
 260         );
 261 }
 262
 263 void OrbitTest::testLngAsc() {
 264         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 265
 266         Orbit orbit;
 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
 271         );
 272         // using an inclination of 90° as well to make the rotation more apparent
 273         orbit.Inclination(PI * 0.5);
 274
 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));
 278         AssertEqual(
 279                 "wrong position at t=0",
 280                 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
 281                 glm::dvec3(pos) / pos.w,
 282                 epsilon
 283         );
 284
 285         pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 286         AssertEqual(
 287                 "wrong position at t=90°",
 288                 glm::dvec3(0.0, 1.0, 0.0),
 289                 glm::dvec3(pos) / pos.w,
 290                 epsilon
 291         );
 292
 293         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 294         AssertEqual(
 295                 "wrong position at t=180°",
 296                 glm::dvec3(-0.707106781186548, 0.0, 0.707106781186548),
 297                 glm::dvec3(pos) / pos.w,
 298                 epsilon
 299         );
 300
 301         pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 302         AssertEqual(
 303                 "wrong position at t=270°",
 304                 glm::dvec3(0.0, -1.0, 0.0),
 305                 glm::dvec3(pos) / pos.w,
 306                 epsilon
 307         );
 308
 309         pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 310         AssertEqual(
 311                 "wrong position at t=360°",
 312                 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
 313                 glm::dvec3(pos) / pos.w,
 314                 epsilon
 315         );
 316 }
 317
 318 void OrbitTest::testArgPe() {
 319         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 320
 321         Orbit orbit;
 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
 326         );
 327         // using an inclination of 90° as well to make the rotation more apparent
 328         orbit.Inclination(PI * 0.5);
 329
 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));
 334         AssertEqual(
 335                 "wrong position at t=0",
 336                 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
 337                 glm::dvec3(pos) / pos.w,
 338                 epsilon
 339         );
 340
 341         pos = orbit.Matrix(PI * 0.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 342         AssertEqual(
 343                 "wrong position at t=90°",
 344                 glm::dvec3(0.0, 1.0, 0.0),
 345                 glm::dvec3(pos) / pos.w,
 346                 epsilon
 347         );
 348
 349         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 350         AssertEqual(
 351                 "wrong position at t=180°",
 352                 glm::dvec3(-0.707106781186548, 0.0, -0.707106781186548),
 353                 glm::dvec3(pos) / pos.w,
 354                 epsilon
 355         );
 356
 357         pos = orbit.Matrix(PI * 1.5) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 358         AssertEqual(
 359                 "wrong position at t=270°",
 360                 glm::dvec3(0.0, -1.0, 0.0),
 361                 glm::dvec3(pos) / pos.w,
 362                 epsilon
 363         );
 364
 365         pos = orbit.Matrix(PI * 2.0) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 366         AssertEqual(
 367                 "wrong position at t=360°",
 368                 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
 369                 glm::dvec3(pos) / pos.w,
 370                 epsilon
 371         );
 372 }
 373
 374 void OrbitTest::testMnAn() {
 375         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 376
 377         Orbit orbit;
 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
 382         );
 383
 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));
 386         AssertEqual(
 387                 "wrong position at t=0",
 388                 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
 389                 glm::dvec3(pos) / pos.w,
 390                 epsilon
 391         );
 392
 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);
 396         AssertEqual(
 397                 "wrong position at t=90°",
 398                 glm::dvec3(-0.707106781186548, 0.0, -0.707106781186548),
 399                 glm::dvec3(pos) / pos.w,
 400                 epsilon
 401         );
 402
 403         // at 180° position should be (-sma,0,0)
 404         pos = orbit.Matrix(PI) * glm::dvec4(0.0, 0.0, 0.0, 1.0);
 405         AssertEqual(
 406                 "wrong position at t=180°",
 407                 glm::dvec3(-0.707106781186548, 0.0, 0.707106781186548),
 408                 glm::dvec3(pos) / pos.w,
 409                 epsilon
 410         );
 411
 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);
 414         AssertEqual(
 415                 "wrong position at t=270°",
 416                 glm::dvec3(0.707106781186548, 0.0, 0.707106781186548),
 417                 glm::dvec3(pos) / pos.w,
 418                 epsilon
 419         );
 420
 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);
 423         AssertEqual(
 424                 "wrong position at t=360°",
 425                 glm::dvec3(0.707106781186548, 0.0, -0.707106781186548),
 426                 glm::dvec3(pos) / pos.w,
 427                 epsilon
 428         );
 429 }
 430
 431 void OrbitTest::testInverseDefault() {
 432         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 433
 434         Orbit orbit;
 435
 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));
 438         AssertEqual(
 439                 "wrong position at t=0",
 440                 glm::dvec3(0.0, 0.0, 0.0),
 441                 glm::dvec3(pos) / pos.w,
 442                 epsilon
 443         );
 444
 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);
 448         AssertEqual(
 449                 "wrong position at t=90°",
 450                 glm::dvec3(0.0, 0.0, 0.0),
 451                 glm::dvec3(pos) / pos.w,
 452                 epsilon
 453         );
 454
 455         // at 180° position should be (-sma,0,0)
 456         pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.0, 0.0, 0.0, 1.0);
 457         AssertEqual(
 458                 "wrong position at t=180°",
 459                 glm::dvec3(0.0, 0.0, 0.0),
 460                 glm::dvec3(pos) / pos.w,
 461                 epsilon
 462         );
 463
 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);
 466         AssertEqual(
 467                 "wrong position at t=270°",
 468                 glm::dvec3(0.0, 0.0, 0.0),
 469                 glm::dvec3(pos) / pos.w,
 470                 epsilon
 471         );
 472
 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);
 475         AssertEqual(
 476                 "wrong position at t=360°",
 477                 glm::dvec3(0.0, 0.0, 0.0),
 478                 glm::dvec3(pos) / pos.w,
 479                 epsilon
 480         );
 481 }
 482
 483 void OrbitTest::testInverseSMA() {
 484         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 485
 486         Orbit orbit;
 487         orbit.SemiMajorAxis(2.0);
 488
 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));
 492         AssertEqual(
 493                 "wrong position at t=0",
 494                 glm::dvec3(0.0, 0.0, 0.0),
 495                 glm::dvec3(pos) / pos.w,
 496                 epsilon
 497         );
 498
 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);
 502         AssertEqual(
 503                 "wrong position at t=90°",
 504                 glm::dvec3(0.0, 0.0, 0.0),
 505                 glm::dvec3(pos) / pos.w,
 506                 epsilon
 507         );
 508
 509         // at 180° position should be (-sma,0,0)
 510         pos = orbit.InverseMatrix(PI) * glm::dvec4(-2.0, 0.0, 0.0, 1.0);
 511         AssertEqual(
 512                 "wrong position at t=180°",
 513                 glm::dvec3(0.0, 0.0, 0.0),
 514                 glm::dvec3(pos) / pos.w,
 515                 epsilon
 516         );
 517
 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);
 520         AssertEqual(
 521                 "wrong position at t=270°",
 522                 glm::dvec3(0.0, 0.0, 0.0),
 523                 glm::dvec3(pos) / pos.w,
 524                 epsilon
 525         );
 526
 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);
 529         AssertEqual(
 530                 "wrong position at t=360°",
 531                 glm::dvec3(0.0, 0.0, 0.0),
 532                 glm::dvec3(pos) / pos.w,
 533                 epsilon
 534         );
 535 }
 536
 537 void OrbitTest::testInverseEcc() {
 538         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 539
 540         Orbit orbit;
 541         orbit.Eccentricity(0.5);
 542
 543         glm::dvec4 pos(orbit.InverseMatrix(0.0) * glm::dvec4(0.5, 0.0, 0.0, 1.0));
 544         AssertEqual(
 545                 "wrong position at t=0",
 546                 glm::dvec3(0.0, 0.0, 0.0),
 547                 glm::dvec3(pos) / pos.w,
 548                 epsilon
 549         );
 550
 551         pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(-0.93513085903671, 0.0, -0.779740887497559, 1.0);
 552         AssertEqual(
 553                 "wrong position at t=90°",
 554                 glm::dvec3(0.0, 0.0, 0.0),
 555                 glm::dvec3(pos) / pos.w,
 556                 epsilon
 557         );
 558
 559         pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.5, 0.0, 0.0, 1.0);
 560         AssertEqual(
 561                 "wrong position at t=180°",
 562                 glm::dvec3(0.0, 0.0, 0.0),
 563                 glm::dvec3(pos) / pos.w,
 564                 epsilon
 565         );
 566
 567         pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(-0.93513085903671, 0.0, 0.779740887497559, 1.0);
 568         AssertEqual(
 569                 "wrong position at t=270°",
 570                 glm::dvec3(0.0, 0.0, 0.0),
 571                 glm::dvec3(pos) / pos.w,
 572                 epsilon
 573         );
 574
 575         pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.5, 0.0, 0.0, 1.0);
 576         AssertEqual(
 577                 "wrong position at t=360°",
 578                 glm::dvec3(0.0, 0.0, 0.0),
 579                 glm::dvec3(pos) / pos.w,
 580                 epsilon
 581         );
 582 }
 583
 584 void OrbitTest::testInverseInc() {
 585         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 586
 587         Orbit orbit;
 588         orbit.Inclination(PI * 0.25); // 45°
 589
 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));
 593         AssertEqual(
 594                 "wrong position at t=0",
 595                 glm::dvec3(0.0, 0.0, 0.0),
 596                 glm::dvec3(pos) / pos.w,
 597                 epsilon
 598         );
 599
 600         pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 0.707106781186548, -0.707106781186548, 1.0);
 601         AssertEqual(
 602                 "wrong position at t=90°",
 603                 glm::dvec3(0.0, 0.0, 0.0),
 604                 glm::dvec3(pos) / pos.w,
 605                 epsilon
 606         );
 607
 608         pos = orbit.InverseMatrix(PI) * glm::dvec4(-1.0, 0.0, 0.0, 1.0);
 609         AssertEqual(
 610                 "wrong position at t=180°",
 611                 glm::dvec3(0.0, 0.0, 0.0),
 612                 glm::dvec3(pos) / pos.w,
 613                 epsilon
 614         );
 615
 616         pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -0.707106781186548, 0.707106781186548, 1.0);
 617         AssertEqual(
 618                 "wrong position at t=270°",
 619                 glm::dvec3(0.0, 0.0, 0.0),
 620                 glm::dvec3(pos) / pos.w,
 621                 epsilon
 622         );
 623
 624         pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(1.0, 0.0, 0.0, 1.0);
 625         AssertEqual(
 626                 "wrong position at t=360°",
 627                 glm::dvec3(0.0, 0.0, 0.0),
 628                 glm::dvec3(pos) / pos.w,
 629                 epsilon
 630         );
 631 }
 632
 633 void OrbitTest::testInverseLngAsc() {
 634         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 635
 636         Orbit orbit;
 637         orbit.LongitudeAscending(PI * 0.25); // 45°
 638         orbit.Inclination(PI * 0.5);
 639
 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));
 643         AssertEqual(
 644                 "wrong position at t=0",
 645                 glm::dvec3(0.0, 0.0, 0.0),
 646                 glm::dvec3(pos) / pos.w,
 647                 epsilon
 648         );
 649
 650         pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 1.0, 0.0, 1.0);
 651         AssertEqual(
 652                 "wrong position at t=90°",
 653                 glm::dvec3(0.0, 0.0, 0.0),
 654                 glm::dvec3(pos) / pos.w,
 655                 epsilon
 656         );
 657
 658         pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, 0.707106781186548, 1.0);
 659         AssertEqual(
 660                 "wrong position at t=180°",
 661                 glm::dvec3(0.0, 0.0, 0.0),
 662                 glm::dvec3(pos) / pos.w,
 663                 epsilon
 664         );
 665
 666         pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -1.0, 0.0, 1.0);
 667         AssertEqual(
 668                 "wrong position at t=270°",
 669                 glm::dvec3(0.0, 0.0, 0.0),
 670                 glm::dvec3(pos) / pos.w,
 671                 epsilon
 672         );
 673
 674         pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.707106781186548, 0.0, -0.707106781186548, 1.0);
 675         AssertEqual(
 676                 "wrong position at t=360°",
 677                 glm::dvec3(0.0, 0.0, 0.0),
 678                 glm::dvec3(pos) / pos.w,
 679                 epsilon
 680         );
 681 }
 682
 683 void OrbitTest::testInverseArgPe() {
 684         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 685
 686         Orbit orbit;
 687         orbit.ArgumentPeriapsis(PI * 0.25); // 45°
 688         orbit.Inclination(PI * 0.5);
 689
 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));
 694         AssertEqual(
 695                 "wrong position at t=0",
 696                 glm::dvec3(0.0, 0.0, 0.0),
 697                 glm::dvec3(pos) / pos.w,
 698                 epsilon
 699         );
 700
 701         pos = orbit.InverseMatrix(PI * 0.5) * glm::dvec4(0.0, 1.0, 0.0, 1.0);
 702         AssertEqual(
 703                 "wrong position at t=90°",
 704                 glm::dvec3(0.0, 0.0, 0.0),
 705                 glm::dvec3(pos) / pos.w,
 706                 epsilon
 707         );
 708
 709         pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, -0.707106781186548, 1.0);
 710         AssertEqual(
 711                 "wrong position at t=180°",
 712                 glm::dvec3(0.0, 0.0, 0.0),
 713                 glm::dvec3(pos) / pos.w,
 714                 epsilon
 715         );
 716
 717         pos = orbit.InverseMatrix(PI * 1.5) * glm::dvec4(0.0, -1.0, 0.0, 1.0);
 718         AssertEqual(
 719                 "wrong position at t=270°",
 720                 glm::dvec3(0.0, 0.0, 0.0),
 721                 glm::dvec3(pos) / pos.w,
 722                 epsilon
 723         );
 724
 725         pos = orbit.InverseMatrix(PI * 2.0) * glm::dvec4(0.707106781186548, 0.0, 0.707106781186548, 1.0);
 726         AssertEqual(
 727                 "wrong position at t=360°",
 728                 glm::dvec3(0.0, 0.0, 0.0),
 729                 glm::dvec3(pos) / pos.w,
 730                 epsilon
 731         );
 732 }
 733
 734 void OrbitTest::testInverseMnAn() {
 735         constexpr double epsilon = 10.0 * std::numeric_limits<double>::epsilon();
 736
 737         Orbit orbit;
 738         orbit.MeanAnomaly(PI * 0.25); // 45°
 739
 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));
 742         AssertEqual(
 743                 "wrong position at t=0",
 744                 glm::dvec3(0.0, 0.0, 0.0),
 745                 glm::dvec3(pos) / pos.w,
 746                 epsilon
 747         );
 748
 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);
 752         AssertEqual(
 753                 "wrong position at t=90°",
 754                 glm::dvec3(0.0, 0.0, 0.0),
 755                 glm::dvec3(pos) / pos.w,
 756                 epsilon
 757         );
 758
 759         // at 180° position should be (-sma,0,0)
 760         pos = orbit.InverseMatrix(PI) * glm::dvec4(-0.707106781186548, 0.0, 0.707106781186548, 1.0);
 761         AssertEqual(
 762                 "wrong position at t=180°",
 763                 glm::dvec3(0.0, 0.0, 0.0),
 764                 glm::dvec3(pos) / pos.w,
 765                 epsilon
 766         );
 767
 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);
 770         AssertEqual(
 771                 "wrong position at t=270°",
 772                 glm::dvec3(0.0, 0.0, 0.0),
 773                 glm::dvec3(pos) / pos.w,
 774                 epsilon
 775         );
 776
 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);
 779         AssertEqual(
 780                 "wrong position at t=360°",
 781                 glm::dvec3(0.0, 0.0, 0.0),
 782                 glm::dvec3(pos) / pos.w,
 783                 epsilon
 784         );
 785 }
 786 }
 787 }
 788 }