From c10b0ac93c0daa0d6a3d8881b3c5cbaeed5f9987 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 25 Nov 2020 14:39:44 +0100
Subject: [PATCH] fixed all of Remy's comments
---
CMakeLists.txt | 1 -
.../detail/framework/geometry/BaseVector.inl | 2 +-
.../framework/geometry/CoordinateSystem.inl | 160 ++++++------
.../detail/framework/geometry/FourVector.inl | 11 +-
corsika/detail/framework/geometry/Helix.inl | 10 +-
corsika/detail/framework/geometry/Line.inl | 10 +-
corsika/detail/framework/geometry/Point.inl | 62 +++--
.../framework/geometry/QuantityVector.inl | 123 +++++----
.../detail/framework/geometry/Trajectory.inl | 6 +-
corsika/detail/framework/geometry/Vector.inl | 244 ++++++++++--------
corsika/detail/framework/utility/COMBoost.inl | 6 +-
corsika/framework/geometry/BaseVector.hpp | 13 +-
.../framework/geometry/CoordinateSystem.hpp | 102 ++++++--
corsika/framework/geometry/FourVector.hpp | 28 +-
corsika/framework/geometry/Helix.hpp | 8 +-
corsika/framework/geometry/Line.hpp | 29 +--
corsika/framework/geometry/Point.hpp | 43 ++-
corsika/framework/geometry/QuantityVector.hpp | 39 +--
corsika/framework/geometry/Sphere.hpp | 13 +-
corsika/framework/geometry/Trajectory.hpp | 6 +-
corsika/framework/geometry/Vector.hpp | 56 ++--
examples/geometry_example.cpp | 12 +-
tests/framework/testGeometry.cpp | 126 +++++----
23 files changed, 652 insertions(+), 458 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 411e5faa3..a6bfab39a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -232,7 +232,6 @@ add_subdirectory (documentation)
#+++++++++++++++++++++++++++++
# modules
#
-set (NO_BOOST 1) # do not use libboost_iostream in CorsikaData
add_subdirectory (modules/data)
add_subdirectory (modules/pythia)
add_subdirectory (modules/sibyll)
diff --git a/corsika/detail/framework/geometry/BaseVector.inl b/corsika/detail/framework/geometry/BaseVector.inl
index 7fbb362e1..e5bd52b19 100644
--- a/corsika/detail/framework/geometry/BaseVector.inl
+++ b/corsika/detail/framework/geometry/BaseVector.inl
@@ -26,7 +26,7 @@ namespace corsika {
}
template <typename TDimension>
- QuantityVector<TDimension>& BaseVector<TDimension>::quantityVector() {
+ QuantityVector<TDimension>& BaseVector<TDimension>::getQuantityVector() {
return quantityVector_;
}
diff --git a/corsika/detail/framework/geometry/CoordinateSystem.inl b/corsika/detail/framework/geometry/CoordinateSystem.inl
index e8eae9b1a..ce5aa9e93 100644
--- a/corsika/detail/framework/geometry/CoordinateSystem.inl
+++ b/corsika/detail/framework/geometry/CoordinateSystem.inl
@@ -20,76 +20,7 @@
namespace corsika {
- inline CoordinateSystemPtr CoordinateSystem::translate(
- QuantityVector<length_d> vector) const {
- EigenTransform const translation{EigenTranslation(vector.eigenVector_)};
-
- return std::make_shared<CoordinateSystem const>(*(new CoordinateSystem(
- std::make_shared<CoordinateSystem const>(*this), translation)));
- }
-
- template <typename TDim>
- CoordinateSystemPtr CoordinateSystem::rotateToZ(Vector<TDim> vVec) const {
- auto const a = vVec.normalized()
- .getComponents(std::make_shared<CoordinateSystem const>(*this))
- .getEigenVector();
- auto const a1 = a(0), a2 = a(1), a3 = a(2);
-
- Eigen::Matrix3d A, B;
-
- if (a3 > 0) {
- auto const c = 1 / (1 + a3);
- 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 {
- auto const c = 1 / (1 - a3);
- 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 std::make_shared<CoordinateSystem const>(*(new CoordinateSystem(
- std::make_shared<CoordinateSystem const>(*this), EigenTransform(A + B))));
- }
-
- template <typename TDim>
- CoordinateSystemPtr CoordinateSystem::rotate(QuantityVector<TDim> axis,
- double angle) const {
- if (axis.eigenVector_.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
-
- EigenTransform const rotation{
- Eigen::AngleAxisd(angle, axis.eigenVector_.normalized())};
-
- return std::make_shared<CoordinateSystem const>(
- CoordinateSystem(std::make_shared<CoordinateSystem const>(*this), rotation));
- }
-
- template <typename TDim>
- CoordinateSystemPtr CoordinateSystem::translateAndRotate(
- QuantityVector<length_d> translation, QuantityVector<TDim> axis, double angle) {
- if (axis.eigenVector_.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
-
- EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eigenVector_.normalized()) *
- EigenTranslation(translation.eigenVector_)};
-
- return std::make_shared<CoordinateSystem const>(CoordinateSystem(*this, transf));
- }
-
- CoordinateSystemPtr CoordinateSystem::getReferenceCS() const {
- return referenceCS_; //*(referenceCS_.get());
- }
+ CoordinateSystemPtr CoordinateSystem::getReferenceCS() const { return referenceCS_; }
EigenTransform const& CoordinateSystem::getTransform() const { return transf_; }
@@ -101,30 +32,27 @@ namespace corsika {
return !(cs == *this);
}
- /**
- * returns the transformation matrix necessary to transform primitives with coordinates
- * in \a pFrom to \a pTo, e.g.
- * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
- * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
- * the indicated CoordinateSystem).
- */
- inline EigenTransform getTransformation(CoordinateSystemPtr pFrom,
- CoordinateSystemPtr pTo) {
+ /// find transformation between two CS, using most optimal common base
+ inline EigenTransform get_transformation(CoordinateSystemPtr const& pFrom,
+ CoordinateSystemPtr const& pTo) {
CoordinateSystemPtr a{pFrom};
CoordinateSystemPtr b{pTo};
CoordinateSystemPtr commonBase{nullptr};
- while (a != b && b != nullptr) {
- a = pFrom;
+ while (a != b && b) {
- while (a != b && a != nullptr) { a = a->getReferenceCS(); }
+ // traverse pFrom
+ a = pFrom;
+ while (a != b && a) {
+ a = a->getReferenceCS();
+ }
if (a == b) break;
b = b->getReferenceCS();
}
- if (a == b && a != nullptr) {
+ if (a == b && a) {
commonBase = a;
} else {
@@ -149,4 +77,70 @@ namespace corsika {
return t;
}
+ inline CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+ QuantityVector<length_d> const& vector) {
+ EigenTransform const translation{EigenTranslation(vector.getEigenVector())};
+ return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, translation));
+ }
+
+ template <typename TDim>
+ inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+ Vector<TDim> const& vVec) {
+ auto const a = vVec.normalized().getComponents(cs).getEigenVector();
+ auto const a1 = a(0), a2 = a(1), a3 = a(2);
+
+ Eigen::Matrix3d A, B;
+
+ if (a3 > 0) {
+ auto const c = 1 / (1 + a3);
+ 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 {
+ auto const c = 1 / (1 - a3);
+ 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 std::make_shared<CoordinateSystem const>(
+ CoordinateSystem(cs, EigenTransform(A + B)));
+ }
+
+ template <typename TDim>
+ inline CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+ QuantityVector<TDim> const& axis,
+ double const angle) {
+ if (axis.getEigenVector().isZero()) {
+ throw std::runtime_error("null-vector given as axis parameter");
+ }
+
+ EigenTransform const rotation{
+ Eigen::AngleAxisd(angle, axis.getEigenVector().normalized())};
+
+ return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, rotation));
+ }
+
+ template <typename TDim>
+ inline CoordinateSystemPtr make_translationAndRotation(
+ CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+ QuantityVector<TDim> const& axis, double const angle) {
+ if (axis.getEigenVector().isZero()) {
+ throw std::runtime_error("null-vector given as axis parameter");
+ }
+
+ EigenTransform const transf{
+ Eigen::AngleAxisd(angle, axis.getEigenVector().normalized()) *
+ EigenTranslation(translation.getEigenVector())};
+
+ return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, transf));
+ }
+
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index bd39af001..454854288 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -1,8 +1,6 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
- * See file AUTHORS for a list of contributors.
- *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
@@ -115,4 +113,13 @@ namespace corsika {
return timeLike_ * timeLike_;
}
+ template <typename TTimeType, typename TSpaceVecType>
+ inline std::ostream& operator<<(std::ostream& os,
+ corsika::FourVector<TTimeType, TSpaceVecType> const qv) {
+
+ os << '(' << qv.timeLike_ << ", " << qv.spaceLike_ << ") ";
+ return os;
+ }
+
+
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Helix.inl b/corsika/detail/framework/geometry/Helix.inl
index 09cdd456c..fe740cf8e 100644
--- a/corsika/detail/framework/geometry/Helix.inl
+++ b/corsika/detail/framework/geometry/Helix.inl
@@ -11,29 +11,29 @@
#pragma once
#include <corsika/framework/geometry/Point.hpp>
-#include <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <cmath>
namespace corsika {
LengthType Helix::getRadius() const { return radius_; }
- Point Helix::getPosition(TimeType t) const {
+ Point Helix::getPosition(TimeType const t) const {
return r0_ + vPar_ * t +
(vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
- Point Helix::getPositionFromArclength(LengthType l) const {
+ Point Helix::getPositionFromArclength(LengthType const l) const {
return getPosition(getTimeFromArclength(l));
}
- LengthType Helix::getArcLength(TimeType t1, TimeType t2) const {
+ LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
- TimeType Helix::getTimeFromArclength(LengthType l) const {
+ TimeType Helix::getTimeFromArclength(LengthType const l) const {
return l / (vPar_ + vPerp_).getNorm();
}
diff --git a/corsika/detail/framework/geometry/Line.inl b/corsika/detail/framework/geometry/Line.inl
index 13f0292d1..e3a5ef6ee 100644
--- a/corsika/detail/framework/geometry/Line.inl
+++ b/corsika/detail/framework/geometry/Line.inl
@@ -1,8 +1,6 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
- * See file AUTHORS for a list of contributors.
- *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
@@ -16,17 +14,17 @@
namespace corsika {
- Point Line::getPosition(TimeType t) const { return start_point_ + velocity_ * t; }
+ Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
- Point Line::getPositionFromArclength(LengthType l) const {
+ Point Line::getPositionFromArclength(LengthType const l) const {
return start_point_ + velocity_.normalized() * l;
}
- LengthType Line::getArcLength(TimeType t1, TimeType t2) const {
+ LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
return velocity_.getNorm() * (t2 - t1);
}
- TimeType Line::getTimeFromArclength(LengthType t) const {
+ TimeType Line::getTimeFromArclength(LengthType const t) const {
return t / velocity_.getNorm();
}
diff --git a/corsika/detail/framework/geometry/Point.inl b/corsika/detail/framework/geometry/Point.inl
index 9a24699f6..119de7f04 100644
--- a/corsika/detail/framework/geometry/Point.inl
+++ b/corsika/detail/framework/geometry/Point.inl
@@ -17,72 +17,82 @@
namespace corsika {
- auto Point::getCoordinates() const { return BaseVector<length_d>::getQuantityVector(); }
+ QuantityVector<length_d> const& Point::getCoordinates() const {
+ return BaseVector<length_d>::getQuantityVector();
+ }
+
+ QuantityVector<length_d>& Point::getCoordinates() {
+ return BaseVector<length_d>::getQuantityVector();
+ }
- inline LengthType Point::getX(CoordinateSystemPtr cs) const {
- if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ inline LengthType Point::getX(CoordinateSystemPtr const& pCS) const {
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+ if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getX();
} else {
return QuantityVector<length_d>(
- getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getX();
}
}
- inline LengthType Point::getY(CoordinateSystemPtr cs) const {
- if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ inline LengthType Point::getY(CoordinateSystemPtr const& pCS) const {
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+ if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getY();
} else {
return QuantityVector<length_d>(
- getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getY();
}
}
- inline LengthType Point::getZ(CoordinateSystemPtr cs) const {
- if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ inline LengthType Point::getZ(CoordinateSystemPtr const& pCS) const {
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+ if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getZ();
} else {
return QuantityVector<length_d>(
- getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getZ();
}
}
/// this always returns a QuantityVector as triple
- auto Point::getCoordinates(CoordinateSystemPtr pCS) const {
- if (*pCS == *BaseVector<length_d>::getCoordinateSystem()) {
+ QuantityVector<length_d> Point::getCoordinates(CoordinateSystemPtr const& pCS) const {
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+ if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector();
} else {
return QuantityVector<length_d>(
- getTransformation(BaseVector<length_d>::getCoordinateSystem(), pCS) *
+ get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_);
}
}
- /*!
- * transforms the Point into another CoordinateSystem by changing its
- * coordinates interally
- */
- void Point::rebase(CoordinateSystemPtr pCS) {
- BaseVector<length_d>::quantityVector() = getCoordinates(pCS);
+ /// this always returns a QuantityVector as triple
+ QuantityVector<length_d>& Point::getCoordinates(CoordinateSystemPtr const& pCS) {
+ if (*pCS != *BaseVector<length_d>::getCoordinateSystem()) { rebase(pCS); }
+ return BaseVector<length_d>::getQuantityVector();
+ }
+
+ void Point::rebase(CoordinateSystemPtr const& pCS) {
+ BaseVector<length_d>::setQuantityVector(QuantityVector<length_d>(
+ get_transformation(BaseVector<length_d>::getCoordinateSystem(), pCS) *
+ BaseVector<length_d>::getQuantityVector().eigenVector_));
BaseVector<length_d>::setCoordinateSystem(pCS);
}
Point Point::operator+(Vector<length_d> const& pVec) const {
- return Point(BaseVector<length_d>::getCoordinateSystem(),
- getCoordinates() +
- pVec.getComponents(BaseVector<length_d>::getCoordinateSystem()));
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+ return Point(cs, getCoordinates() + pVec.getComponents(cs));
}
- /*!
- * returns the distance Vector between two points
- */
Vector<length_d> Point::operator-(Point const& pB) const {
- CoordinateSystemPtr cs = BaseVector<length_d>::getCoordinateSystem();
+ CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
}
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 2b206d3e3..3b5118abe 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -18,58 +18,63 @@
namespace corsika {
- template <typename dim>
- inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::operator[](
- size_t index) const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
+ operator[](size_t const index) const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
}
- template <typename dim>
- inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getX() const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_type
+ QuantityVector<TDimension>::getX() const {
return (*this)[0];
}
- template <typename dim>
- inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getY() const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_type
+ QuantityVector<TDimension>::getY() const {
return (*this)[1];
}
- template <typename dim>
- inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getZ() const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_type
+ QuantityVector<TDimension>::getZ() const {
return (*this)[2];
}
- template <typename dim>
- inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getNorm()
- const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_type
+ QuantityVector<TDimension>::getNorm() const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_.norm());
}
- template <typename dim>
- inline auto QuantityVector<dim>::getSquaredNorm() const {
+ template <typename TDimension>
+ inline typename QuantityVector<TDimension>::quantity_square_type
+ QuantityVector<TDimension>::getSquaredNorm() const {
using QuantitySquared =
decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
return QuantitySquared(phys::units::detail::magnitude_tag,
eigenVector_.squaredNorm());
}
- template <typename dim>
- inline QuantityVector<dim> QuantityVector<dim>::operator+(
- QuantityVector<dim> const& pQVec) const {
- return QuantityVector<dim>(eigenVector_ + pQVec.eigenVector_);
+ template <typename TDimension>
+ inline QuantityVector<TDimension> QuantityVector<TDimension>::operator+(
+ QuantityVector<TDimension> const& pQVec) const {
+ return QuantityVector<TDimension>(eigenVector_ + pQVec.eigenVector_);
}
- template <typename dim>
- inline QuantityVector<dim> QuantityVector<dim>::operator-(
- QuantityVector<dim> const& pQVec) const {
- return QuantityVector<dim>(eigenVector_ - pQVec.eigenVector_);
+ template <typename TDimension>
+ inline QuantityVector<TDimension> QuantityVector<TDimension>::operator-(
+ QuantityVector<TDimension> const& pQVec) const {
+ return QuantityVector<TDimension>(eigenVector_ - pQVec.eigenVector_);
}
- template <typename dim>
- template <typename ScalarDim>
- inline auto QuantityVector<dim>::operator*(
- phys::units::quantity<ScalarDim, double> const p) const {
- using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
+ template <typename TDimension>
+ template <typename TScalarDim>
+ inline auto QuantityVector<TDimension>::operator*(
+ phys::units::quantity<TScalarDim, double> const p) const {
+ using ResQuantity =
+ phys::units::detail::Product<TScalarDim, TDimension, double, double>;
if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
// a "quantity_type" anymore
@@ -81,69 +86,73 @@ namespace corsika {
}
}
- template <typename dim>
- template <typename ScalarDim>
- inline auto QuantityVector<dim>::operator/(
- phys::units::quantity<ScalarDim, double> const p) const {
+ template <typename TDimension>
+ template <typename TScalarDim>
+ inline auto QuantityVector<TDimension>::operator/(
+ phys::units::quantity<TScalarDim, double> const p) const {
return (*this) * (1 / p);
}
- template <typename dim>
- inline auto QuantityVector<dim>::operator*(double const p) const {
- return QuantityVector<dim>(eigenVector_ * p);
+ template <typename TDimension>
+ inline auto QuantityVector<TDimension>::operator*(double const p) const {
+ return QuantityVector<TDimension>(eigenVector_ * p);
}
- template <typename dim>
- inline auto QuantityVector<dim>::operator/(double const p) const {
- return QuantityVector<dim>(eigenVector_ / p);
+ template <typename TDimension>
+ inline auto QuantityVector<TDimension>::operator/(double const p) const {
+ return QuantityVector<TDimension>(eigenVector_ / p);
}
- template <typename dim>
- inline auto& QuantityVector<dim>::operator/=(double const p) {
+ template <typename TDimension>
+ inline auto& QuantityVector<TDimension>::operator/=(double const p) {
eigenVector_ /= p;
return *this;
}
- template <typename dim>
- inline auto& QuantityVector<dim>::operator*=(double const p) {
+ template <typename TDimension>
+ inline auto& QuantityVector<TDimension>::operator*=(double const p) {
eigenVector_ *= p;
return *this;
}
- template <typename dim>
- inline auto& QuantityVector<dim>::operator+=(QuantityVector<dim> const& pQVec) {
+ template <typename TDimension>
+ inline auto& QuantityVector<TDimension>::operator+=(
+ QuantityVector<TDimension> const& pQVec) {
eigenVector_ += pQVec.eigenVector_;
return *this;
}
- template <typename dim>
- inline auto& QuantityVector<dim>::operator-=(QuantityVector<dim> const& pQVec) {
+ template <typename TDimension>
+ inline auto& QuantityVector<TDimension>::operator-=(
+ QuantityVector<TDimension> const& pQVec) {
eigenVector_ -= pQVec.eigenVector_;
return *this;
}
- template <typename dim>
- inline auto& QuantityVector<dim>::operator-() const {
- return QuantityVector<dim>(-eigenVector_);
+ template <typename TDimension>
+ inline auto& QuantityVector<TDimension>::operator-() const {
+ return QuantityVector<TDimension>(-eigenVector_);
}
- template <typename dim>
- inline auto QuantityVector<dim>::normalized() const {
- return QuantityVector<dim>(eigenVector_.normalized());
+ template <typename TDimension>
+ inline auto QuantityVector<TDimension>::normalized() const {
+ return QuantityVector<TDimension>(eigenVector_.normalized());
}
- template <typename dim>
- inline auto QuantityVector<dim>::operator==(QuantityVector<dim> const& p) const {
+ template <typename TDimension>
+ inline auto QuantityVector<TDimension>::operator==(
+ QuantityVector<TDimension> const& p) const {
return eigenVector_ == p.eigenVector_;
}
- template <typename dim>
- inline std::ostream& operator<<(std::ostream& os, corsika::QuantityVector<dim> qv) {
- using quantity_type = phys::units::quantity<dim, double>;
+ template <typename TDimension>
+ inline std::ostream& operator<<(std::ostream& os,
+ corsika::QuantityVector<TDimension> const qv) {
+ using quantity_type = phys::units::quantity<TDimension, double>;
os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
<< qv.eigenVector_(2) << ") "
- << phys::units::to_unit_symbol<dim, double>(
+ << phys::units::to_unit_symbol<TDimension, double>(
quantity_type(phys::units::detail::magnitude_tag, 1));
return os;
}
diff --git a/corsika/detail/framework/geometry/Trajectory.inl b/corsika/detail/framework/geometry/Trajectory.inl
index 89b02ad95..a43010c19 100644
--- a/corsika/detail/framework/geometry/Trajectory.inl
+++ b/corsika/detail/framework/geometry/Trajectory.inl
@@ -17,7 +17,7 @@
namespace corsika {
template <typename TType>
- Point Trajectory<TType>::getPosition(double u) const {
+ Point Trajectory<TType>::getPosition(double const u) const {
return TType::getPosition(timeLength_ * u);
}
@@ -32,14 +32,14 @@ namespace corsika {
}
template <typename TType>
- LengthType Trajectory<TType>::getDistance(TimeType t) const {
+ LengthType Trajectory<TType>::getDistance(TimeType const t) const {
assert(t <= timeLength_);
assert(t >= 0 * second);
return TType::getArcLength(0 * second, t);
}
template <typename TType>
- void Trajectory<TType>::getLimitEndTo(LengthType limit) {
+ void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
timeLength_ = TType::getTimeFromArclength(limit);
}
diff --git a/corsika/detail/framework/geometry/Vector.inl b/corsika/detail/framework/geometry/Vector.inl
index 1df45dfc6..d11250be8 100644
--- a/corsika/detail/framework/geometry/Vector.inl
+++ b/corsika/detail/framework/geometry/Vector.inl
@@ -1,8 +1,6 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
- * See file AUTHORS for a list of contributors.
- *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
@@ -16,183 +14,209 @@
namespace corsika {
- template <typename dim>
- QuantityVector<dim> Vector<dim>::getComponents() const {
- return BaseVector<dim>::getQuantityVector();
+ template <typename TDimension>
+ inline QuantityVector<TDimension> const& Vector<TDimension>::getComponents() const {
+ return BaseVector<TDimension>::getQuantityVector();
+ }
+
+ template <typename TDimension>
+ inline QuantityVector<TDimension>& Vector<TDimension>::getComponents() {
+ return BaseVector<TDimension>::getQuantityVector();
}
- template <typename dim>
- QuantityVector<dim> Vector<dim>::getComponents(CoordinateSystemPtr pCS) const {
- if (*pCS == *BaseVector<dim>::getCoordinateSystem()) { // FIXME
- return BaseVector<dim>::getQuantityVector();
+ template <typename TDimension>
+ inline QuantityVector<TDimension> Vector<TDimension>::getComponents(
+ CoordinateSystemPtr const& pCS) const {
+ if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+ return BaseVector<TDimension>::getQuantityVector();
} else {
- return QuantityVector<dim>(
- getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
- BaseVector<dim>::getQuantityVector().eigenVector_);
+ return QuantityVector<TDimension>(
+ get_transformation(BaseVector<TDimension>::getCoordinateSystem(), pCS)
+ .linear() *
+ BaseVector<TDimension>::getQuantityVector().eigenVector_);
}
}
- template <typename dim>
- inline typename Vector<dim>::quantity_type Vector<dim>::getX(
- CoordinateSystemPtr pCS) const {
- if (*pCS == *BaseVector<dim>::getCoordinateSystem()) {
- return BaseVector<dim>::getQuantityVector()[0];
+ template <typename TDimension>
+ inline QuantityVector<TDimension>& Vector<TDimension>::getComponents(
+ CoordinateSystemPtr const& pCS) {
+ if (*pCS != *BaseVector<TDimension>::getCoordinateSystem()) { rebase(pCS); }
+ return BaseVector<TDimension>::getQuantityVector();
+ }
+
+ template <typename TDimension>
+ inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getX(
+ CoordinateSystemPtr const& pCS) const {
+ if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+ return BaseVector<TDimension>::getQuantityVector()[0];
} else {
- return QuantityVector<dim>(
- getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
- BaseVector<dim>::getQuantityVector().eigenVector_)[0];
+ return QuantityVector<TDimension>(
+ get_transformation(BaseVector<TDimension>::getCoordinateSystem(), pCS)
+ .linear() *
+ BaseVector<TDimension>::getQuantityVector().eigenVector_)[0];
}
}
- template <typename dim>
- inline typename Vector<dim>::quantity_type Vector<dim>::getY(
- CoordinateSystemPtr pCS) const {
- if (*pCS == *BaseVector<dim>::getCoordinateSystem()) {
- return BaseVector<dim>::getQuantityVector()[1];
+ template <typename TDimension>
+ inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getY(
+ CoordinateSystemPtr const& pCS) const {
+ if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+ return BaseVector<TDimension>::getQuantityVector()[1];
} else {
- return QuantityVector<dim>(
- getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
- BaseVector<dim>::getQuantityVector().eigenVector_)[1];
+ return QuantityVector<TDimension>(
+ get_transformation(BaseVector<TDimension>::getCoordinateSystem(), pCS)
+ .linear() *
+ BaseVector<TDimension>::getQuantityVector().eigenVector_)[1];
}
}
- template <typename dim>
- inline typename Vector<dim>::quantity_type Vector<dim>::getZ(
- CoordinateSystemPtr pCS) const {
- if (*pCS == *BaseVector<dim>::getCoordinateSystem()) {
- return BaseVector<dim>::getQuantityVector()[2];
+ template <typename TDimension>
+ inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getZ(
+ CoordinateSystemPtr const& pCS) const {
+ if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+ return BaseVector<TDimension>::getQuantityVector()[2];
} else {
- return QuantityVector<dim>(
- getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
- BaseVector<dim>::getQuantityVector().eigenVector_)[2];
+ return QuantityVector<TDimension>(
+ get_transformation(BaseVector<TDimension>::getCoordinateSystem(), pCS)
+ .linear() *
+ BaseVector<TDimension>::getQuantityVector().eigenVector_)[2];
}
}
- template <typename dim>
- void Vector<dim>::rebase(CoordinateSystemPtr pCS) {
- BaseVector<dim>::quantityVector() = getComponents(pCS);
- BaseVector<dim>::setCoordinateSystem(pCS);
+ template <typename TDimension>
+ void Vector<TDimension>::rebase(CoordinateSystemPtr const& pCS) {
+ BaseVector<TDimension>::setQuantityVector(QuantityVector<TDimension>(
+ get_transformation(BaseVector<TDimension>::getCoordinateSystem(), pCS).linear() *
+ BaseVector<TDimension>::getQuantityVector().eigenVector_));
+ BaseVector<TDimension>::setCoordinateSystem(pCS);
}
- template <typename dim>
- inline typename Vector<dim>::quantity_type Vector<dim>::getNorm() const {
- return BaseVector<dim>::getQuantityVector().getNorm();
+ template <typename TDimension>
+ inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getNorm() const {
+ return BaseVector<TDimension>::getQuantityVector().getNorm();
}
- template <typename dim>
- auto Vector<dim>::getSquaredNorm() const {
- return BaseVector<dim>::getQuantityVector().getSquaredNorm();
+ template <typename TDimension>
+ inline typename Vector<TDimension>::quantity_square_type
+ Vector<TDimension>::getSquaredNorm() const {
+ return BaseVector<TDimension>::getQuantityVector().getSquaredNorm();
}
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::getParallelProjectionOnto(Vector<dim2> const& pVec,
- CoordinateSystemPtr pCS) const {
+ template <typename TDimension>
+ template <typename TDimension2>
+ auto Vector<TDimension>::getParallelProjectionOnto(
+ Vector<TDimension2> const& pVec, CoordinateSystemPtr const& pCS) const {
auto const ourCompVec = getComponents(pCS);
auto const otherCompVec = pVec.getComponents(pCS);
auto const& a = ourCompVec.eigenVector_;
auto const& b = otherCompVec.eigenVector_;
- return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
+ return Vector<TDimension>(
+ pCS, QuantityVector<TDimension>(b * ((a.dot(b)) / b.squaredNorm())));
}
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::getParallelProjectionOnto(Vector<dim2> const& pVec) const {
- return getParallelProjectionOnto<dim2>(pVec, BaseVector<dim>::getCoordinateSystem());
+ template <typename TDimension>
+ template <typename TDimension2>
+ auto Vector<TDimension>::getParallelProjectionOnto(
+ Vector<TDimension2> const& pVec) const {
+ return getParallelProjectionOnto<TDimension2>(
+ pVec, BaseVector<TDimension>::getCoordinateSystem());
}
- template <typename dim>
- Vector<dim> Vector<dim>::operator+(Vector<dim> const& pVec) const {
- auto const components = getComponents(BaseVector<dim>::getCoordinateSystem()) +
- pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
- return Vector<dim>(BaseVector<dim>::getCoordinateSystem(), components);
+ template <typename TDimension>
+ Vector<TDimension> Vector<TDimension>::operator+(Vector<TDimension> const& pVec) const {
+ CoordinateSystemPtr const& cs = BaseVector<TDimension>::getCoordinateSystem();
+ auto const components = getComponents(cs) + pVec.getComponents(cs);
+ return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(), components);
}
- template <typename dim>
- Vector<dim> Vector<dim>::operator-(Vector<dim> const& pVec) const {
- auto const components =
- getComponents() - pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
- return Vector<dim>(BaseVector<dim>::getCoordinateSystem(), components);
+ template <typename TDimension>
+ Vector<TDimension> Vector<TDimension>::operator-(Vector<TDimension> const& pVec) const {
+ CoordinateSystemPtr const& cs = BaseVector<TDimension>::getCoordinateSystem();
+ return Vector<TDimension>(cs, getComponents() - pVec.getComponents(cs));
}
- template <typename dim>
- auto& Vector<dim>::operator*=(double const p) {
- BaseVector<dim>::quantityVector() *= p;
+ template <typename TDimension>
+ auto& Vector<TDimension>::operator*=(double const p) {
+ BaseVector<TDimension>::getQuantityVector() *= p;
return *this;
}
- template <typename dim>
- template <typename ScalarDim>
- auto Vector<dim>::operator*(phys::units::quantity<ScalarDim, double> const p) const {
- using ProdDim = phys::units::detail::product_d<dim, ScalarDim>;
+ template <typename TDimension>
+ template <typename TScalarDim>
+ auto Vector<TDimension>::operator*(
+ phys::units::quantity<TScalarDim, double> const p) const {
+ using ProdDim = phys::units::detail::product_d<TDimension, TScalarDim>;
- return Vector<ProdDim>(BaseVector<dim>::getCoordinateSystem(),
- BaseVector<dim>::getQuantityVector() * p);
+ return Vector<ProdDim>(BaseVector<TDimension>::getCoordinateSystem(),
+ BaseVector<TDimension>::getQuantityVector() * p);
}
- template <typename dim>
- template <typename ScalarDim>
- auto Vector<dim>::operator/(phys::units::quantity<ScalarDim, double> const p) const {
+ template <typename TDimension>
+ template <typename TScalarDim>
+ auto Vector<TDimension>::operator/(
+ phys::units::quantity<TScalarDim, double> const p) const {
return (*this) * (1 / p);
}
- template <typename dim>
- auto Vector<dim>::operator*(double const p) const {
- return Vector<dim>(BaseVector<dim>::getCoordinateSystem(),
- BaseVector<dim>::getQuantityVector() * p);
+ template <typename TDimension>
+ auto Vector<TDimension>::operator*(double const p) const {
+ return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+ BaseVector<TDimension>::getQuantityVector() * p);
}
- template <typename dim>
- auto Vector<dim>::operator/(double const p) const {
- return Vector<dim>(BaseVector<dim>::getCoordinateSystem(),
- BaseVector<dim>::quantityVector() / p);
+ template <typename TDimension>
+ auto Vector<TDimension>::operator/(double const p) const {
+ return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+ BaseVector<TDimension>::quantityVector() / p);
}
- template <typename dim>
- auto& Vector<dim>::operator+=(Vector<dim> const& pVec) {
- BaseVector<dim>::quantityVector() +=
- pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
+ template <typename TDimension>
+ auto& Vector<TDimension>::operator+=(Vector<TDimension> const& pVec) {
+ BaseVector<TDimension>::getQuantityVector() +=
+ pVec.getComponents(BaseVector<TDimension>::getCoordinateSystem());
return *this;
}
- template <typename dim>
- auto& Vector<dim>::operator-=(Vector<dim> const& pVec) {
- BaseVector<dim>::quantityVector() -=
- pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
+ template <typename TDimension>
+ auto& Vector<TDimension>::operator-=(Vector<TDimension> const& pVec) {
+ BaseVector<TDimension>::getQuantityVector() -=
+ pVec.getComponents(BaseVector<TDimension>::getCoordinateSystem());
return *this;
}
- template <typename dim>
- auto& Vector<dim>::operator-() const {
- return Vector<dim>(BaseVector<dim>::getCoordinateSystem(),
- -BaseVector<dim>::quantityVector());
+ template <typename TDimension>
+ auto& Vector<TDimension>::operator-() const {
+ return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+ -BaseVector<TDimension>::getQuantityVector());
}
- template <typename dim>
- auto Vector<dim>::normalized() const {
+ template <typename TDimension>
+ auto Vector<TDimension>::normalized() const {
return (*this) * (1 / getNorm());
}
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::cross(Vector<dim2> pV) const {
+ template <typename TDimension>
+ template <typename TDimension2>
+ auto Vector<TDimension>::cross(Vector<TDimension2> const& pV) const {
auto const c1 = getComponents().eigenVector_;
- auto const c2 = pV.getComponents(BaseVector<dim>::getCoordinateSystem()).eigenVector_;
+ auto const c2 =
+ pV.getComponents(BaseVector<TDimension>::getCoordinateSystem()).eigenVector_;
auto const bareResult = c1.cross(c2);
- using ProdDim = phys::units::detail::product_d<dim, dim2>;
- return Vector<ProdDim>(BaseVector<dim>::getCoordinateSystem(), bareResult);
+ using ProdDim = phys::units::detail::product_d<TDimension, TDimension2>;
+ return Vector<ProdDim>(BaseVector<TDimension>::getCoordinateSystem(), bareResult);
}
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::dot(Vector<dim2> pV) const {
+ template <typename TDimension>
+ template <typename TDimension2>
+ auto Vector<TDimension>::dot(Vector<TDimension2> const& pV) const {
auto const c1 = getComponents().eigenVector_;
- auto const c2 = pV.getComponents(BaseVector<dim>::getCoordinateSystem()).eigenVector_;
+ auto const c2 =
+ pV.getComponents(BaseVector<TDimension>::getCoordinateSystem()).eigenVector_;
auto const bareResult = c1.dot(c2);
- using ProdDim = phys::units::detail::product_d<dim, dim2>;
+ using ProdDim = phys::units::detail::product_d<TDimension, TDimension2>;
return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
bareResult);
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index cef0ce97c..36111a2b9 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -26,7 +26,7 @@ namespace corsika {
COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
HEPMassType const massTarget)
: originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
- , rotatedCS_{originalCS_->rotateToZ(Pprojectile.getSpaceLikeComponents())} {
+ , rotatedCS_{make_rotationToZ(originalCS_, Pprojectile.getSpaceLikeComponents())} {
auto const pProjectile = Pprojectile.getSpaceLikeComponents();
auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
@@ -48,7 +48,7 @@ namespace corsika {
COMBoost::COMBoost(Vector<hepmomentum_d> const& momentum, HEPEnergyType mass)
: originalCS_{momentum.getCoordinateSystem()}
- , rotatedCS_{originalCS_->rotateToZ(momentum)} {
+ , rotatedCS_{make_rotationToZ(originalCS_, momentum)} {
auto const squaredNorm = momentum.getSquaredNorm();
auto const norm = sqrt(squaredNorm);
auto const sinhEta = -norm / mass;
@@ -89,7 +89,7 @@ namespace corsika {
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
- pCM.eigenVector()(2) = boostedZ(1) * (1_GeV).magnitude();
+ pCM.getEigenVector()[2] = boostedZ(1) * (1_GeV).magnitude();
Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
pLab.rebase(originalCS_);
diff --git a/corsika/framework/geometry/BaseVector.hpp b/corsika/framework/geometry/BaseVector.hpp
index 21819df29..5569149c6 100644
--- a/corsika/framework/geometry/BaseVector.hpp
+++ b/corsika/framework/geometry/BaseVector.hpp
@@ -18,29 +18,32 @@ namespace corsika {
/*!
* Common base class for Vector and Point.
*
- * This holds a QuantityVector and a CoordinateSystem
+ * This holds a QuantityVector and a CoordinateSystem. The
+ * BaseVector manages resources for many geometry objects.
*
*/
template <typename TDimension>
class BaseVector {
public:
- BaseVector(CoordinateSystemPtr pCS, QuantityVector<TDimension> const& pQVector)
+ BaseVector(CoordinateSystemPtr const& pCS, QuantityVector<TDimension> const& pQVector)
: quantityVector_(pQVector)
, cs_(pCS) {}
- BaseVector() = delete;
+ BaseVector() = delete; // we only want to creat initialized
+ // objects
BaseVector(BaseVector const&) = default;
BaseVector(BaseVector&& a) = default;
BaseVector& operator=(BaseVector const&) = default;
~BaseVector() = default;
CoordinateSystemPtr getCoordinateSystem() const;
- void setCoordinateSystem(CoordinateSystemPtr cs) { cs_ = cs; }
+ void setCoordinateSystem(CoordinateSystemPtr const& cs) { cs_ = cs; }
protected:
QuantityVector<TDimension> const& getQuantityVector() const;
- QuantityVector<TDimension>& quantityVector();
+ QuantityVector<TDimension>& getQuantityVector();
+ void setQuantityVector(QuantityVector<TDimension> const& v) { quantityVector_=v; }
private:
QuantityVector<TDimension> quantityVector_;
diff --git a/corsika/framework/geometry/CoordinateSystem.hpp b/corsika/framework/geometry/CoordinateSystem.hpp
index 6cd5c4813..45d34e699 100644
--- a/corsika/framework/geometry/CoordinateSystem.hpp
+++ b/corsika/framework/geometry/CoordinateSystem.hpp
@@ -10,6 +10,7 @@ n/*
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/logging/Logging.hpp>
#include <Eigen/Dense>
#include <stdexcept>
@@ -32,19 +33,48 @@ namespace corsika {
class RootCoordinateSystem; // fwd decl
static CoordinateSystemPtr get_root_CoordinateSystem(); // fwd decl
+ inline CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+ QuantityVector<length_d> const& vector);
+
+ /**
+ * creates a new CS in which vVec points in direction of the new z-axis, \a vVec
+ */
+ template <typename TDim>
+ inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+ Vector<TDim> const& vVec);
+
+ template <typename TDim>
+ inline CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+ QuantityVector<TDim> const& axis,
+ double const angle);
+
+ template <typename TDim>
+ inline CoordinateSystemPtr make_translationAndRotation(
+ CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+ QuantityVector<TDim> const& axis, double const angle);
+
/**
- * A class to store the reference for a geometric object
+ * A class to store the reference coordinate system for a geometric object
*
* A CoordinateSystem can only be created in reference and relative
- * to other CoordinateSystem sytems. Thus, the geometric
- * transformation between all CoordinateSystems is known.
+ * to other CoordinateSystems. Thus, the geometric
+ * transformation between all CoordinateSystems is always known and stored.
*
- * Only the \sa RootCoordinateSystem can be created as a singleton
- * as global main reference point.
+ * The static (sigleton) function \sa make_root_CoordinateSystem is
+ * the only way to create and access the global top-level
+ * CoordinateSystem obect. CoordinateSystem objects should be
+ * *abosulte* *only* handled in their form of CoordinateSystemPtr,
+ * which are shared_ptr that handle the lifetime of the entire
+ * CoordinateSystem hirarchy.
*
- * Thus, new CoordinateSystems can only be created using
- * transformation: \sa rotateToZ, \sa rotate, \sa translateAndRotate
+ * Thus, new CoordinateSystem are only be created (via
+ * CoordinateSystemPtr) by transforing existing CoordinateSystem
+ * using: \sa rotateToZ, \sa rotate, or \sa translateAndRotate, see
* below.
+ *
+ * Warning: As a consequence, never try to access, modify, copy, the raw
+ * CoordinateSystem objects directly, this will almost certainly result in undefined
+ * behaviour. Only access, copy, handle them via CoordinateSystemPtr.
*/
class CoordinateSystem {
@@ -52,7 +82,7 @@ namespace corsika {
/**
* Constructor only from referenceCS, given the transformation matrix transf
*/
- CoordinateSystem(CoordinateSystemPtr referenceCS, EigenTransform const& transf)
+ CoordinateSystem(CoordinateSystemPtr const& referenceCS, EigenTransform const& transf)
: referenceCS_(referenceCS)
, transf_(transf) {}
@@ -67,23 +97,14 @@ namespace corsika {
// default resource allocation
CoordinateSystem(CoordinateSystem const&) = default;
CoordinateSystem(CoordinateSystem&&) = default;
- CoordinateSystem& operator=(CoordinateSystem const& pCS) = default;
+ CoordinateSystem& operator=(CoordinateSystem const& pCS) =
+ delete; // avoid making copies
~CoordinateSystem() = default;
- inline CoordinateSystemPtr translate(QuantityVector<length_d> vector) const;
-
/**
- * creates a new CS in which vVec points in direction of the new z-axis
+ * Checks, if this is the unique ROOT CS
*/
- template <typename TDim>
- CoordinateSystemPtr rotateToZ(Vector<TDim> vVec) const;
-
- template <typename TDim>
- CoordinateSystemPtr rotate(QuantityVector<TDim> axis, double angle) const;
-
- template <typename TDim>
- CoordinateSystemPtr translateAndRotate(QuantityVector<length_d> translation,
- QuantityVector<TDim> axis, double angle);
+ inline bool isRoot() const { return !referenceCS_; }
inline CoordinateSystemPtr getReferenceCS() const;
@@ -95,15 +116,46 @@ namespace corsika {
protected:
static CoordinateSystem createCS() { return CoordinateSystem(); }
- friend CoordinateSystemPtr get_root_CoordinateSystem(); /// this is the only way to
- /// create ONE unique root CS
+ /**
+ * \defgroup manipulation and creation function
+ * \{
+ **/
+
+ /** this is the only way to create ONE unique root CS **/
+ friend CoordinateSystemPtr get_root_CoordinateSystem();
+
+ friend CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+ QuantityVector<length_d> const& vector);
+ template <typename TDim>
+ friend CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+ Vector<TDim> const& vVec);
+ template <typename TDim>
+ friend CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+ QuantityVector<TDim> const& axis,
+ double const angle);
+ template <typename TDim>
+ friend CoordinateSystemPtr make_translationAndRotation(
+ CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+ QuantityVector<TDim> const& axis, double const angle);
+
+ /** \} **/
private:
- std::shared_ptr<CoordinateSystem const> referenceCS_;
+ CoordinateSystemPtr referenceCS_;
EigenTransform transf_;
};
- EigenTransform getTransformation(CoordinateSystemPtr c1, CoordinateSystemPtr c2);
+ /**
+ * Transformation matrix from one reference system to another.
+ *
+ * returns the transformation matrix necessary to transform primitives with coordinates
+ * in \a pFrom to \a pTo, e.g.
+ * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
+ * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
+ * the indicated CoordinateSystem).
+ */
+ inline EigenTransform get_transformation(CoordinateSystemPtr const& c1,
+ CoordinateSystemPtr const& c2);
} // namespace corsika
diff --git a/corsika/framework/geometry/FourVector.hpp b/corsika/framework/geometry/FourVector.hpp
index 55fefc2fb..4f4bea4de 100644
--- a/corsika/framework/geometry/FourVector.hpp
+++ b/corsika/framework/geometry/FourVector.hpp
@@ -59,16 +59,17 @@ namespace corsika {
decltype(std::declval<norm_type>() * std::declval<norm_type>());
public:
+ // resource management
FourVector() = default;
-
+ FourVector(FourVector&&) = default;
+ FourVector(FourVector const&) = default;
+ FourVector& operator=(const FourVector&) = default;
+ ~FourVector() = default;
+
FourVector(TTimeType const& eT, TSpaceVecType const& eS)
: timeLike_(eT)
, spaceLike_(eS) {}
- /*
- * FIXME: These Getters are mis-leading and does not favor
- * locality. Adhere to Getter/Setter
- */
/**
*
* @return timeLike_
@@ -106,18 +107,13 @@ namespace corsika {
*
* RU: then you have to decide in the constructor which avoids "lazyness"
*/
- /**
- *
- * @return
- */
bool isTimelike() const;
+ bool isSpacelike() const;
/**
* @defgroup math operators (class members)
* @{
*/
- bool isSpacelike() const;
-
FourVector& operator+=(FourVector const&);
FourVector& operator-=(FourVector const&);
FourVector& operator*=(double const);
@@ -174,7 +170,7 @@ namespace corsika {
}
friend FourVector<time_type, space_vec_type> operator/(FourVector const& a,
- double b) {
+ double const b) {
return FourVector<time_type, space_vec_type>(a.timeLike_ / b, a.spaceLike_ / b);
}
@@ -188,6 +184,14 @@ namespace corsika {
norm_square_type getTimeSquared() const;
};
+ /**
+ * streaming operator
+ **/
+
+ template <typename TTimeType, typename TSpaceVecType>
+ inline std::ostream& operator<<(std::ostream& os,
+ corsika::FourVector<TTimeType, TSpaceVecType> const& qv);
+
} // namespace corsika
#include <corsika/detail/framework/geometry/FourVector.inl>
diff --git a/corsika/framework/geometry/Helix.hpp b/corsika/framework/geometry/Helix.hpp
index 05ca40be0..826b4c3fc 100644
--- a/corsika/framework/geometry/Helix.hpp
+++ b/corsika/framework/geometry/Helix.hpp
@@ -46,13 +46,13 @@ namespace corsika {
inline LengthType getRadius() const;
- inline Point getPosition(TimeType t) const;
+ inline Point getPosition(TimeType const t) const;
- inline Point getPositionFromArclength(LengthType l) const;
+ inline Point getPositionFromArclength(LengthType const l) const;
- inline LengthType getArcLength(TimeType t1, TimeType t2) const;
+ inline LengthType getArcLength(TimeType const t1, TimeType const t2) const;
- inline TimeType getTimeFromArclength(LengthType l) const;
+ inline TimeType getTimeFromArclength(LengthType const l) const;
private:
Point r0_; ///! origin of helix, but this is in the center of the
diff --git a/corsika/framework/geometry/Line.hpp b/corsika/framework/geometry/Line.hpp
index 2c27a0bc0..3e1e18564 100644
--- a/corsika/framework/geometry/Line.hpp
+++ b/corsika/framework/geometry/Line.hpp
@@ -15,15 +15,10 @@ n/*
namespace corsika {
/**
- * \class Line
+ * Describes a straight line in space
*
- * A Line describes a movement in three dimensional space. It
- * consists of a Point `$\vec{p_0}$` and and a speed-Vector
- * `$\vec{v}$`, so that it can return GetPosition as
- * `$\vec{p_0}*\vec{v}*t$` for any value of time `$t$`.
- *
- **/
-
+ */
+
class Line {
///! \toto move this to PhysicalUnits
@@ -34,21 +29,23 @@ namespace corsika {
: start_point_(pR0)
, velocity_(pV0) {}
- inline Point getPosition(TimeType t) const;
+ inline Point getPosition(TimeType const t) const;
- inline Point getPositionFromArclength(LengthType l) const;
+ inline Point getPositionFromArclength(LengthType const l) const;
- inline LengthType getArcLength(TimeType t1, TimeType t2) const;
+ inline LengthType getArcLength(TimeType const t1, TimeType const t2) const;
- inline TimeType getTimeFromArclength(LengthType t) const;
+ inline TimeType getTimeFromArclength(LengthType const t) const;
- inline const Point& getStartPoint() const;
+ inline Point const& getStartPoint() const;
+ inline Point& startPoint() { return start_point_; }
- inline const VelocityVec& getVelocity() const;
+ inline VelocityVec const& getVelocity() const;
+ inline VelocityVec& velocity() { return velocity_; }
private:
- Point const start_point_;
- VelocityVec const velocity_;
+ Point start_point_;
+ VelocityVec velocity_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Point.hpp b/corsika/framework/geometry/Point.hpp
index 06c302e1b..c302d14c8 100644
--- a/corsika/framework/geometry/Point.hpp
+++ b/corsika/framework/geometry/Point.hpp
@@ -22,29 +22,54 @@ namespace corsika {
class Point : public BaseVector<length_d> {
public:
- Point(CoordinateSystemPtr pCS, QuantityVector<length_d> const& pQVector)
+ Point(CoordinateSystemPtr const& pCS, QuantityVector<length_d> const& pQVector)
: BaseVector<length_d>(pCS, pQVector) {}
- Point(CoordinateSystemPtr cs, LengthType x, LengthType y, LengthType z)
+ Point(CoordinateSystemPtr const& cs, LengthType x, LengthType y, LengthType z)
: BaseVector<length_d>(cs, {x, y, z}) {}
/** \todo TODO: this should be private or protected, we don NOT want to expose numbers
* without reference to outside:
*/
- inline auto getCoordinates() const;
+ inline QuantityVector<length_d> const& getCoordinates() const;
+ inline QuantityVector<length_d>& getCoordinates();
- /// this always returns a QuantityVector as triple
- inline auto getCoordinates(CoordinateSystemPtr pCS) const;
+ /**
+ this always returns a QuantityVector as triple
- inline LengthType getX(CoordinateSystemPtr pCS) const;
- inline LengthType getY(CoordinateSystemPtr pCS) const;
- inline LengthType getZ(CoordinateSystemPtr pCS) const;
+ \returns A value type QuantityVector, since it may have to create a temporary
+ object to transform to pCS.
+ **/
+ inline QuantityVector<length_d> getCoordinates(CoordinateSystemPtr const& pCS) const;
+
+ /**
+ * this always returns a QuantityVector as triple
+ *
+ * \returns A reference type QuantityVector&, but be aware, the underlying class data
+ * is actually transformed to pCS, if needed. Thus, there may be an implicit call to
+ * \sa rebase.
+ **/
+ inline QuantityVector<length_d>& getCoordinates(CoordinateSystemPtr const& pCS);
+
+ /**
+ * \defgroup access coordinate components
+ * \{
+ *
+ * Note, if you access components in a different CoordinateSystem
+ * pCS than the stored data, internally a temporary object will be
+ * created and destroyed each call. This can be avoided by using
+ * \sa rebase first.
+ **/
+ inline LengthType getX(CoordinateSystemPtr const& pCS) const;
+ inline LengthType getY(CoordinateSystemPtr const& pCS) const;
+ inline LengthType getZ(CoordinateSystemPtr const& pCS) const;
+ /** \} **/
/*!
* transforms the Point into another CoordinateSystem by changing its
* coordinates interally
*/
- inline void rebase(CoordinateSystemPtr pCS);
+ inline void rebase(CoordinateSystemPtr const& pCS);
inline Point operator+(Vector<length_d> const& pVec) const;
diff --git a/corsika/framework/geometry/QuantityVector.hpp b/corsika/framework/geometry/QuantityVector.hpp
index e1f56d380..bc7e41e4e 100644
--- a/corsika/framework/geometry/QuantityVector.hpp
+++ b/corsika/framework/geometry/QuantityVector.hpp
@@ -17,9 +17,9 @@ n/*
namespace corsika {
class CoordinateSystem; // fwd decl
- class Point;
+ class Point; // fwd decl
template <typename T>
- class Vector;
+ class Vector; // fwd decl
/*!
* A QuantityVector is a three-component container based on Eigen::Vector3d
@@ -37,43 +37,45 @@ namespace corsika {
using quantity_type =
phys::units::quantity<TDimension, double>; //< the phys::units::quantity
// corresponding to the dimension
+ using quantity_square_type =
+ decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
- QuantityVector(Eigen::Vector3d pBareVector)
+ QuantityVector(Eigen::Vector3d const& pBareVector)
: eigenVector_(pBareVector) {}
public:
typedef TDimension dimension_type; //!< should be a phys::units::dimension
- QuantityVector(quantity_type a, quantity_type b, quantity_type c)
+ QuantityVector(quantity_type const a, quantity_type const b, quantity_type const c)
: eigenVector_{a.magnitude(), b.magnitude(), c.magnitude()} {}
- QuantityVector(double a, double b, double c)
+ QuantityVector(double const a, double const b, double const c)
: eigenVector_{a, b, c} {
static_assert(
std::is_same_v<TDimension, phys::units::dimensionless_d>,
- "initialization of dimensionful QuantityVector with pure numbers not allowed!");
+ "initialization of dimensionfull QuantityVector with pure numbers not allowed!");
}
- quantity_type operator[](size_t index) const;
+ quantity_type operator[](size_t const index) const;
quantity_type getX() const;
quantity_type getY() const;
quantity_type getZ() const;
Eigen::Vector3d const& getEigenVector() const { return eigenVector_; }
- Eigen::Vector3d& eigenVector() { return eigenVector_; }
+ Eigen::Vector3d& getEigenVector() { return eigenVector_; }
quantity_type getNorm() const;
- auto getSquaredNorm() const;
+ quantity_square_type getSquaredNorm() const;
QuantityVector operator+(QuantityVector<TDimension> const& pQVec) const;
QuantityVector operator-(QuantityVector<TDimension> const& pQVec) const;
- template <typename ScalarDim>
- auto operator*(phys::units::quantity<ScalarDim, double> const p) const;
+ template <typename TScalarDim>
+ auto operator*(phys::units::quantity<TScalarDim, double> const p) const;
- template <typename ScalarDim>
- auto operator/(phys::units::quantity<ScalarDim, double> const p) const;
+ template <typename TScalarDim>
+ auto operator/(phys::units::quantity<TScalarDim, double> const p) const;
auto operator*(double const p) const;
@@ -93,25 +95,26 @@ namespace corsika {
auto operator==(QuantityVector<TDimension> const& p) const;
+ // friends:
friend class CoordinateSystem;
friend class Point;
template <typename T>
friend class corsika::Vector;
- template <typename dim>
- friend std::ostream& operator<<(std::ostream& os, QuantityVector<dim> qv);
+ template <typename TDim>
+ friend std::ostream& operator<<(std::ostream& os, QuantityVector<TDim> qv);
protected:
Eigen::Vector3d
eigenVector_; //!< the actual container where the raw numbers are stored
};
- /*
+ /**
* streaming operator
- */
+ **/
template <typename TDimension>
inline std::ostream& operator<<(std::ostream& os,
- corsika::QuantityVector<TDimension> qv);
+ corsika::QuantityVector<TDimension> const qv);
} // namespace corsika
diff --git a/corsika/framework/geometry/Sphere.hpp b/corsika/framework/geometry/Sphere.hpp
index d5f4c44a8..e014dd86d 100644
--- a/corsika/framework/geometry/Sphere.hpp
+++ b/corsika/framework/geometry/Sphere.hpp
@@ -14,6 +14,11 @@ n/*
namespace corsika {
+ /**
+ * Describes a sphere in space
+ *
+ **/
+
class Sphere : public IVolume {
public:
@@ -24,13 +29,15 @@ namespace corsika {
//! returns true if the Point p is within the sphere
inline bool isInside(Point const& p) const override;
- inline const Point& getCenter() const;
+ inline Point const& getCenter() const;
+ inline Point& getCenter() { return center_; }
inline LengthType getRadius() const;
+ inline LengthType& getRadius() { return radius_; }
private:
- Point const center_;
- LengthType const radius_;
+ Point center_;
+ LengthType radius_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Trajectory.hpp b/corsika/framework/geometry/Trajectory.hpp
index 0c26f9dcd..9aa7b917d 100644
--- a/corsika/framework/geometry/Trajectory.hpp
+++ b/corsika/framework/geometry/Trajectory.hpp
@@ -25,15 +25,15 @@ namespace corsika {
: TType(theT)
, timeLength_(timeLength) {}
- Point getPosition(double u) const;
+ Point getPosition(double const u) const;
TimeType getDuration() const;
LengthType getLength() const;
- LengthType getDistance(TimeType t) const;
+ LengthType getDistance(TimeType const t) const;
- void getLimitEndTo(LengthType limit);
+ void getLimitEndTo(LengthType const limit);
auto getNormalizedDirection() const;
diff --git a/corsika/framework/geometry/Vector.hpp b/corsika/framework/geometry/Vector.hpp
index 96a154562..83ea39dcf 100644
--- a/corsika/framework/geometry/Vector.hpp
+++ b/corsika/framework/geometry/Vector.hpp
@@ -15,7 +15,9 @@ n/*
namespace corsika {
/*!
- * A Vector represents a 3-vector in Euclidean space. It is defined by components
+ * A Vector represents a 3-vector in Euclidean space.
+ *
+ * It is defined by components
* given in a specific CoordinateSystem. It has a physical dimension ("unit")
* as part of its type, so you cannot mix up e.g. electric with magnetic fields
* (but you could calculate their cross-product to get an energy flux vector).
@@ -28,38 +30,61 @@ namespace corsika {
class Vector : public BaseVector<TDimension> {
public:
using quantity_type = phys::units::quantity<TDimension, double>;
+ using quantity_square_type =
+ decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
- public:
- Vector(CoordinateSystemPtr pCS, QuantityVector<TDimension> pQVector)
+ Vector(CoordinateSystemPtr const& pCS, QuantityVector<TDimension> const& pQVector)
: BaseVector<TDimension>(pCS, pQVector) {}
- Vector(CoordinateSystemPtr cs, quantity_type x, quantity_type y, quantity_type z)
+ Vector(CoordinateSystemPtr const& cs, quantity_type const x, quantity_type const y,
+ quantity_type const z)
: BaseVector<TDimension>(cs, QuantityVector<TDimension>(x, y, z)) {}
/*!
- * returns a QuantityVector with the components given in the "home"
+ * \returns a QuantityVector with the components given in the "home"
* CoordinateSystem of the Vector
*
* \todo this should best be protected, we don't want users to use
* bare coordinates without reference frame
*/
- QuantityVector<TDimension> getComponents() const;
+ inline QuantityVector<TDimension> const& getComponents() const;
+ inline QuantityVector<TDimension>& getComponents();
/*!
* returns a QuantityVector with the components given in an arbitrary
* CoordinateSystem
*/
- QuantityVector<TDimension> getComponents(CoordinateSystemPtr pCS) const;
+ inline QuantityVector<TDimension> getComponents(CoordinateSystemPtr const& pCS) const;
+
+ /**
+ * this always returns a QuantityVector as triple
+ *
+ * \returns A reference type QuantityVector&, but be aware, the underlying class data
+ * is actually transformed to pCS, if needed. Thus, there may be an implicit call to
+ * \sa rebase.
+ **/
+ inline QuantityVector<TDimension>& getComponents(CoordinateSystemPtr const& pCS);
+
+ /**
+ * \defgroup access coordinate components
+ * \{
+ *
+ * Note, if you access components in a different CoordinateSystem
+ * pCS than the stored data, internally a temporary object will be
+ * created and destroyed each call. This can be avoided by using
+ * \sa rebase first.
+ **/
- inline quantity_type getX(CoordinateSystemPtr pCS) const;
- inline quantity_type getY(CoordinateSystemPtr pCS) const;
- inline quantity_type getZ(CoordinateSystemPtr pCS) const;
+ inline quantity_type getX(CoordinateSystemPtr const& pCS) const;
+ inline quantity_type getY(CoordinateSystemPtr const& pCS) const;
+ inline quantity_type getZ(CoordinateSystemPtr const& pCS) const;
+ /** \} **/
/*!
* transforms the Vector into another CoordinateSystem by changing
* its components internally
*/
- inline void rebase(CoordinateSystemPtr pCS);
+ inline void rebase(CoordinateSystemPtr const& pCS);
/*!
* returns the norm/length of the Vector. Before using this method,
@@ -71,7 +96,8 @@ namespace corsika {
* returns the squared norm of the Vector. Before using this method,
* think about whether norm() might be cheaper for your computation.
*/
- inline auto getSquaredNorm() const;
+ inline quantity_square_type getSquaredNorm() const;
+
/*!
* returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection
* of this vector \f$ \vec{v}_1 \f$ along another Vector \f$ \vec{v}_2 \f$ given by
@@ -81,7 +107,7 @@ namespace corsika {
*/
template <typename TDimension2>
auto getParallelProjectionOnto(Vector<TDimension2> const& pVec,
- CoordinateSystemPtr pCS) const;
+ CoordinateSystemPtr const& pCS) const;
template <typename TDimension2>
auto getParallelProjectionOnto(Vector<TDimension2> const& pVec) const;
@@ -109,10 +135,10 @@ namespace corsika {
auto normalized() const;
template <typename TDimension2>
- auto cross(Vector<TDimension2> pV) const;
+ auto cross(Vector<TDimension2> const& pV) const;
template <typename TDimension2>
- auto dot(Vector<TDimension2> pV) const;
+ auto dot(Vector<TDimension2> const& pV) const;
};
} // namespace corsika
diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp
index 0458c257f..5cdf6043c 100644
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -11,28 +11,26 @@ n/*
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/logging/Logging.hpp>
#include <cstdlib>
#include <iostream>
#include <typeinfo>
using namespace corsika;
-using namespace corsika::units::si;
int main() {
- std::cout << "geometry_example" << std::endl;
-
// define the root coordinate system
- corsika::CoordinateSystemPtr root = get_root_CoordinateSystem();
+ CoordinateSystemPtr root = get_root_CoordinateSystem();
// another CS defined by a translation relative to the root CS
- CoordinateSystemPtr cs2 = root->translate({0_m, 0_m, 1_m});
+ CoordinateSystemPtr cs2 = make_translation(root, {0_m, 0_m, 1_m});
// rotations are possible, too; parameters are axis vector and angle
CoordinateSystemPtr cs3 =
- root->rotate(QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
-
+ make_rotation(root, 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
Point const p2(cs2, {0_m, 0_m, 0_m}); // the origin of cs2
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index e303e6419..1c03de501 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -24,76 +24,68 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
- REQUIRE(getTransformation(rootCS, rootCS).isApprox(EigenTransform::Identity()));
-
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
Point p1(rootCS, coordinates);
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
- REQUIRE((p1.getCoordinates() - coordinates).getNorm().magnitude() ==
- Approx(0).margin(absMargin));
- REQUIRE((p1.getCoordinates(rootCS) - coordinates).getNorm().magnitude() ==
- Approx(0).margin(absMargin));
-
- /*
- SECTION("unconnected CoordinateSystems") {
- CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
- REQUIRE_THROWS(CoordinateSystem::getTransformation(rootCS, rootCS2));
- }*/
+ CHECK((p1.getCoordinates() - coordinates).getNorm().magnitude() ==
+ Approx(0).margin(absMargin));
+ CHECK((p1.getCoordinates(rootCS) - coordinates).getNorm().magnitude() ==
+ Approx(0).margin(absMargin));
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
- CoordinateSystemPtr translatedCS = rootCS->translate(translationVector);
+ CoordinateSystemPtr translatedCS = make_translation(rootCS, translationVector);
- REQUIRE(*translatedCS->getReferenceCS() == *rootCS);
+ CHECK(*translatedCS->getReferenceCS() == *rootCS);
- REQUIRE((p1.getCoordinates(translatedCS) + translationVector).getNorm().magnitude() ==
- Approx(0).margin(absMargin));
+ CHECK((p1.getCoordinates(translatedCS) + translationVector).getNorm().magnitude() ==
+ Approx(0).margin(absMargin));
// Vectors are not subject to translations
- REQUIRE((v1.getComponents(rootCS) - v1.getComponents(translatedCS))
- .getNorm()
- .magnitude() == Approx(0).margin(absMargin));
+ CHECK((v1.getComponents(rootCS) - v1.getComponents(translatedCS))
+ .getNorm()
+ .magnitude() == Approx(0).margin(absMargin));
Point p2(translatedCS, {0_m, 0_m, 0_m});
- REQUIRE(((p2 - p1).getComponents() - translationVector).getNorm().magnitude() ==
- Approx(0).margin(absMargin));
+ CHECK(((p2 - p1).getComponents() - translationVector).getNorm().magnitude() ==
+ Approx(0).margin(absMargin));
}
SECTION("multiple translations") {
QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
- CoordinateSystemPtr cs2 = rootCS->translate(tv1);
+ CoordinateSystemPtr cs2 = make_translation(rootCS, tv1);
QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
- CoordinateSystemPtr cs3 = rootCS->translate(tv2);
+ CoordinateSystemPtr cs3 = make_translation(rootCS, tv2);
QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
- CoordinateSystemPtr cs4 = cs3->translate(tv3);
+ CoordinateSystemPtr cs4 = make_translation(cs3, tv3);
- REQUIRE(*cs4->getReferenceCS()->getReferenceCS() == *rootCS);
+ CHECK(*cs4->getReferenceCS()->getReferenceCS() == *rootCS);
- REQUIRE(getTransformation(cs3, cs2).isApprox(
- rootCS->translate({3_m, -5_m, 0_m})->getTransform()));
- REQUIRE(getTransformation(cs2, cs3).isApprox(
- rootCS->translate({-3_m, +5_m, 0_m})->getTransform()));
+ CHECK(get_transformation(cs3, cs2).isApprox(
+ make_translation(rootCS, {3_m, -5_m, 0_m})->getTransform()));
+ CHECK(get_transformation(cs2, cs3).isApprox(
+ make_translation(rootCS, {-3_m, +5_m, 0_m})->getTransform()));
}
SECTION("rotations") {
QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
double const angle = 90. / 180. * M_PI;
- CoordinateSystemPtr rotatedCS = rootCS->rotate(axis, angle);
- REQUIRE(*rotatedCS->getReferenceCS() == *rootCS);
+ CoordinateSystemPtr rotatedCS = make_rotation(rootCS, axis, angle);
+ CHECK(*rotatedCS->getReferenceCS() == *rootCS);
- REQUIRE(v1.getComponents(rotatedCS)[1].magnitude() ==
- Approx((-1. * tesla).magnitude()));
+ CHECK(v1.getComponents(rotatedCS)[1].magnitude() ==
+ Approx((-1. * tesla).magnitude()));
// vector norm invariant under rotation
- REQUIRE(v1.getComponents(rotatedCS).getNorm().magnitude() ==
- Approx(v1.getComponents(rootCS).getNorm().magnitude()));
+ CHECK(v1.getComponents(rotatedCS).getNorm().magnitude() ==
+ Approx(v1.getComponents(rootCS).getNorm().magnitude()));
}
SECTION("multiple rotations") {
@@ -107,20 +99,20 @@ TEST_CASE("transformations between CoordinateSystems") {
double const angle = 90. / 180. * M_PI;
- CoordinateSystemPtr rotated1 = rootCS->rotate(zAxis, angle);
- CoordinateSystemPtr rotated2 = rotated1->rotate(yAxis, angle);
- CoordinateSystemPtr rotated3 = rotated2->rotate(zAxis, -angle);
+ CoordinateSystemPtr rotated1 = make_rotation(rootCS, zAxis, angle);
+ CoordinateSystemPtr rotated2 = make_rotation(rotated1, yAxis, angle);
+ CoordinateSystemPtr rotated3 = make_rotation(rotated2, zAxis, -angle);
- CoordinateSystemPtr combined = rootCS->rotate(xAxis, -angle);
+ CoordinateSystemPtr combined = make_rotation(rootCS, xAxis, -angle);
auto comp1 = v1.getComponents(rotated3);
auto comp3 = v1.getComponents(combined);
- REQUIRE((comp1 - comp3).getNorm().magnitude() == Approx(0).margin(absMargin));
+ CHECK((comp1 - comp3).getNorm().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 csPrime = make_rotationToZ(rootCS, v);
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};
@@ -153,7 +145,7 @@ TEST_CASE("transformations between CoordinateSystems") {
SECTION("RotateToZ negative") {
Vector const v{rootCS, 0_m, 0_m, -1_m};
- auto const csPrime = rootCS->rotateToZ(v);
+ auto const csPrime = make_rotationToZ(rootCS, v);
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};
@@ -241,6 +233,52 @@ TEST_CASE("transformations between CoordinateSystems") {
}
}
+TEST_CASE("CoordinateSystem hirarchy") {
+
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
+
+ CHECK(get_transformation(rootCS, rootCS).isApprox(EigenTransform::Identity()));
+
+ // define the root coordinate system
+ CoordinateSystemPtr root = get_root_CoordinateSystem();
+ Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS
+
+ // root -> cs2
+ CoordinateSystemPtr cs2 = make_translation(root, {0_m, 0_m, 1_m});
+ Point const p2(cs2, {0_m, 0_m, -1_m});
+
+ // root -> cs2 -> cs3
+ CoordinateSystemPtr cs3 = make_translation(cs2, {0_m, 0_m, -1_m});
+ Point const p3(cs3, {0_m, 0_m, 0_m});
+
+ // root -> cs2 -> cs4
+ CoordinateSystemPtr cs4 = make_translation(cs2, {0_m, 0_m, -1_m});
+ Point const p4(cs4, {0_m, 0_m, 0_m});
+
+ // root -> cs2 -> cs4 -> cs5
+ CoordinateSystemPtr cs5 =
+ make_rotation(cs4, QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
+ Point const p5(cs5, {0_m, 0_m, 0_m});
+
+ // root -> cs6
+ CoordinateSystemPtr cs6 =
+ make_rotation(root, QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
+ Point const p6(cs6, {0_m, 0_m, 0_m}); // the origin of the root CS
+
+ // all points should be on top of each other
+
+ CHECK_FALSE(get_transformation(root, cs2).isApprox(EigenTransform::Identity()));
+ CHECK(get_transformation(root, cs3).isApprox(EigenTransform::Identity()));
+ CHECK(get_transformation(root, cs4).isApprox(EigenTransform::Identity()));
+ CHECK(get_transformation(cs5, cs6).isApprox(EigenTransform::Identity()));
+
+ CHECK((p1 - p2).getNorm().magnitude() == Approx(0).margin(absMargin));
+ CHECK((p1 - p3).getNorm().magnitude() == Approx(0).margin(absMargin));
+ CHECK((p1 - p4).getNorm().magnitude() == Approx(0).margin(absMargin));
+ CHECK((p1 - p5).getNorm().magnitude() == Approx(0).margin(absMargin));
+ CHECK((p1 - p6).getNorm().magnitude() == Approx(0).margin(absMargin));
+}
+
TEST_CASE("Sphere") {
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
@@ -255,8 +293,8 @@ TEST_CASE("Sphere") {
}
SECTION("isInside") {
- REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
- REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
+ CHECK_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
+ CHECK(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
}
}
--
GitLab