From 24d3a64fc219c8b6437edd65a4acd781b5cdcb0f Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 15 Apr 2019 17:58:56 -0300
Subject: [PATCH] added CoordinateSystem::RotateToZ and rotation with arbitrary
dim. vector
---
Documentation/Examples/geometry_example.cc | 3 +-
Framework/Geometry/CMakeLists.txt | 1 +
Framework/Geometry/CoordinateSystem.h | 39 ++++++++++++++++++++--
3 files changed, 40 insertions(+), 3 deletions(-)
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 9c589ee71..8ed0d6f65 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -32,7 +32,8 @@ int main() {
CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
// rotations are possible, too; parameters are axis vector and angle
- CoordinateSystem cs3 = root.rotate({1_m, 0_m, 0_m}, 90 * degree_angle);
+ CoordinateSystem cs3 =
+ root.rotate(QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
// now let's define some geometrical objects:
Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index f827de1f4..c12c79b01 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -40,6 +40,7 @@ set_target_properties (
target_link_libraries (
CORSIKAgeometry
CORSIKAunits
+ CORSIKAutilities
)
target_include_directories (
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index aaf226f38..195e89bb8 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -14,6 +14,7 @@
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/sgn.h>
#include <Eigen/Dense>
#include <stdexcept>
@@ -23,6 +24,8 @@ typedef Eigen::Translation<double, 3> EigenTranslation;
namespace corsika::geometry {
class RootCoordinateSystem;
+ template <typename T>
+ class Vector;
using corsika::units::si::length_d;
@@ -59,7 +62,38 @@ namespace corsika::geometry {
return CoordinateSystem(*this, translation);
}
- auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
+ template <typename TDim>
+ auto RotateToZ(Vector<TDim> vVec) const {
+ auto const a = vVec.normalized().GetComponents(*this).eVector;
+ auto const a1 = a(0), a2 = a(1);
+
+ auto const s = utl::sgn(a(2));
+ auto const c = 1 / (1 + s * a(2));
+
+ Eigen::Matrix3d A, B;
+
+ if (s > 0) {
+ 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, 0, (a1 * a1 + a2 * a2) * c; // .
+ }
+
+ return CoordinateSystem(*this, EigenTransform(A + B));
+ }
+
+ template <typename TDim>
+ auto rotate(QuantityVector<TDim> axis, double angle) const {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
@@ -69,8 +103,9 @@ namespace corsika::geometry {
return CoordinateSystem(*this, rotation);
}
+ template <typename TDim>
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
- QuantityVector<phys::units::length_d> axis, double angle) {
+ QuantityVector<TDim> axis, double angle) {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
--
GitLab