test and fix of CoordinateSystem::RotateToZ()

Framework/Geometry/CoordinateSystem.h

@@ -72,19 +72,19 @@ namespace corsika::geometry {
       Eigen::Matrix3d A, B;
 
       if (s > 0) {
-        A << 1, 0, -a1,                     // comment to prevent clang-format
-            0, 1, -a2,                      // .
-            a1, a2, 1;                      // .
+        A << 1, 0, a1,                      // comment to prevent clang-format
+            0, 1, a2,                       // .
+            -a1, -a2, 1;                    // .
         B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
             -a1 * a2 * c, -a2 * a2 * c, 0,  // .
             0, 0, -(a1 * a1 + a2 * a2) * c; // .
 
       } else {
         A << 1, 0, a1,                      // .
-            0, -1, -a2,                     // .
-            a1, a2, -1;                     // .
-        B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-            +a1 * a2 * c, +a2 * a2 * c, 0,  // .
+            0, -1, a2,                      // .
+            a1, -a2, -1;                    // .
+        B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
+            -a1 * a2 * c, +a2 * a2 * c, 0,  // .
             0, 0, (a1 * a1 + a2 * a2) * c;  // .
       }
 

Framework/Geometry/testGeometry.cc

@@ -123,6 +123,63 @@ TEST_CASE("transformations between CoordinateSystems") {
     auto comp3 = v1.GetComponents(combined);
     REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
   }
+
+  SECTION("RotateToZ positive") {
+    Vector const v{rootCS, 0_m, 1_m, 1_m};
+    auto const csPrime = rootCS.RotateToZ(v);
+    auto const transform = csPrime.GetTransform().matrix();
+    Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
+    Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
+    Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
+
+    CHECK(zPrime.dot(v).magnitude() > 0);
+    CHECK(zPrime.GetComponents(rootCS)[1].magnitude() ==
+          Approx(zPrime.GetComponents(rootCS)[2].magnitude()));
+    CHECK(zPrime.GetComponents(rootCS)[0].magnitude() == Approx(0));
+
+    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+              yPrime.GetComponents(rootCS).eVector) == Approx(0));
+    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+              xPrime.GetComponents(rootCS).eVector) == Approx(0));
+    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+              zPrime.GetComponents(rootCS).eVector) == Approx(0));
+
+    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+              yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+              xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+              zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+  }
+
+  SECTION("RotateToZ negative") {
+    Vector const v{rootCS, 0_m, 0_m, -1_m};
+    auto const csPrime = rootCS.RotateToZ(v);
+    auto const transform = csPrime.GetTransform().matrix();
+    Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
+    Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
+    Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
+
+    CHECK(zPrime.dot(v).magnitude() > 0);
+    CHECK(xPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
+          Approx(0));
+    CHECK(yPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
+          Approx(0));
+
+    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+              yPrime.GetComponents(rootCS).eVector) == Approx(0));
+    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+              xPrime.GetComponents(rootCS).eVector) == Approx(0));
+    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+              zPrime.GetComponents(rootCS).eVector) == Approx(0));
+
+    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+              yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+              xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+              zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+  }
 }
 
 TEST_CASE("Sphere") {