diff --git a/corsika/detail/framework/geometry/BaseVector.inl b/corsika/detail/framework/geometry/BaseVector.inl
index b1b7361ceb8943e80441986e4962acb20bd9302f..7fbb362e198dce7b0db3577c16f31f3974b73449 100644
--- a/corsika/detail/framework/geometry/BaseVector.inl
+++ b/corsika/detail/framework/geometry/BaseVector.inl
@@ -15,12 +15,19 @@
namespace corsika {
- template <typename dim>
- auto const& BaseVector<dim>::GetCoordinateSystem() const
- {
- return *cs;
- }
+ template <typename TDimension>
+ CoordinateSystemPtr BaseVector<TDimension>::getCoordinateSystem() const {
+ return cs_;
+ }
+ template <typename TDimension>
+ QuantityVector<TDimension> const& BaseVector<TDimension>::getQuantityVector() const {
+ return quantityVector_;
+ }
-} // namespace corsika
+ template <typename TDimension>
+ QuantityVector<TDimension>& BaseVector<TDimension>::quantityVector() {
+ return quantityVector_;
+ }
+} // namespace corsika
diff --git a/corsika/detail/framework/geometry/CoordinateSystem.inl b/corsika/detail/framework/geometry/CoordinateSystem.inl
index 4ec390bd2900887257b0c2a52caabfcb2c85bc4f..e8eae9b1a82d88c888701bc1894e1299e49afd90 100644
--- a/corsika/detail/framework/geometry/CoordinateSystem.inl
+++ b/corsika/detail/framework/geometry/CoordinateSystem.inl
@@ -11,140 +11,142 @@
#pragma once
#include <corsika/framework/geometry/QuantityVector.hpp>
-#include <Eigen/Dense>
-#include <stdexcept>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <Eigen/Dense>
-namespace corsika {
+#include <memory>
+#include <stdexcept>
+namespace corsika {
- CoordinateSystem& CoordinateSystem::operator=(const corsika::CoordinateSystem& pCS)
- {
- reference = pCS.reference;
- transf = pCS.transf;
- return *this;
+ 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; // .
}
- inline CoordinateSystem CoordinateSystem::translate(QuantityVector<length_d> vector) const
- {
- EigenTransform const translation{EigenTranslation(vector.eVector)};
-
- return CoordinateSystem(*this, translation);
- }
+ return std::make_shared<CoordinateSystem const>(*(new CoordinateSystem(
+ std::make_shared<CoordinateSystem const>(*this), EigenTransform(A + B))));
+ }
- template <typename TDim>
- auto CoordinateSystem::RotateToZ(Vector<TDim> vVec) const
- {
- auto const a = vVec.normalized().GetComponents(*this).eVector;
- 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 CoordinateSystem(*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");
}
- template <typename TDim>
- auto CoordinateSystem::rotate(QuantityVector<TDim> axis, double angle) const
- {
- if (axis.eVector.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
+ EigenTransform const rotation{
+ Eigen::AngleAxisd(angle, axis.eigenVector_.normalized())};
- EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
+ return std::make_shared<CoordinateSystem const>(
+ CoordinateSystem(std::make_shared<CoordinateSystem const>(*this), rotation));
+ }
- return CoordinateSystem(*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");
}
- template <typename TDim>
- auto CoordinateSystem::translateAndRotate(QuantityVector<length_d> translation, QuantityVector<TDim> axis, double angle)
- {
- if (axis.eVector.isZero()) {
- throw std::runtime_error("null-vector given as axis parameter");
- }
+ EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eigenVector_.normalized()) *
+ EigenTranslation(translation.eigenVector_)};
- EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
- EigenTranslation(translation.eVector)};
+ return std::make_shared<CoordinateSystem const>(CoordinateSystem(*this, transf));
+ }
- return CoordinateSystem(*this, transf);
- }
+ CoordinateSystemPtr CoordinateSystem::getReferenceCS() const {
+ return referenceCS_; //*(referenceCS_.get());
+ }
- CoordinateSystem const* CoordinateSystem::GetReference() const
- {
- return reference;
- }
+ EigenTransform const& CoordinateSystem::getTransform() const { return transf_; }
- const EigenTransform& CoordinateSystem::GetTransform() const
- {
- return transf;
- }
+ inline bool CoordinateSystem::operator==(CoordinateSystem const& cs) const {
+ return referenceCS_ == cs.referenceCS_ && transf_.matrix() == cs.transf_.matrix();
+ }
- /**
- * 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(CoordinateSystem const& pFrom,
- CoordinateSystem const& pTo) {
- CoordinateSystem const* a{&pFrom};
- CoordinateSystem const* b{&pTo};
- CoordinateSystem const* commonBase{nullptr};
+ inline bool CoordinateSystem::operator!=(CoordinateSystem const& cs) const {
+ return !(cs == *this);
+ }
- while (a != b && b != nullptr) {
- a = &pFrom;
+ /**
+ * 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) {
+ CoordinateSystemPtr a{pFrom};
+ CoordinateSystemPtr b{pTo};
+ CoordinateSystemPtr commonBase{nullptr};
- while (a != b && a != nullptr) { a = a->GetReference(); }
+ while (a != b && b != nullptr) {
+ a = pFrom;
- if (a == b) break;
+ while (a != b && a != nullptr) { a = a->getReferenceCS(); }
- b = b->GetReference();
- }
+ if (a == b) break;
- if (a == b && a != nullptr) {
- commonBase = a;
+ b = b->getReferenceCS();
+ }
- } else {
- throw std::runtime_error("no connection between coordinate systems found!");
- }
+ if (a == b && a != nullptr) {
+ commonBase = a;
- EigenTransform t = EigenTransform::Identity();
- auto* p = &pFrom;
+ } else {
+ throw std::runtime_error("no connection between coordinate systems found!");
+ }
- while (p != commonBase) {
- t = p->GetTransform() * t;
- p = p->GetReference();
- }
+ EigenTransform t = EigenTransform::Identity();
+ CoordinateSystemPtr p = pFrom;
- p = &pTo;
+ while ((*p) != (*commonBase)) {
+ t = p->getTransform() * t;
+ p = p->getReferenceCS();
+ }
- while (p != commonBase) {
- t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
- p = p->GetReference();
- }
+ p = pTo;
- return t;
+ while (*p != *commonBase) {
+ t = t * p->getTransform().inverse(Eigen::TransformTraits::Isometry);
+ p = p->getReferenceCS();
}
-} // namespace corsika
+ return t;
+ }
+} // namespace corsika
diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index f3de79d9e81518f034ead3b309c4f7186997a9e3..bd39af0014b891de70b566f270db523bb2b22286 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -27,7 +27,7 @@ namespace corsika {
}
template <typename TTimeType, typename TSpaceVecType>
- const TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents()
+ TSpaceVecType const& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents()
const {
return spaceLike_;
}
@@ -35,7 +35,7 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
typename FourVector<TTimeType, TSpaceVecType>::norm_square_type
FourVector<TTimeType, TSpaceVecType>::getNormSqr() const {
- return getTimeSquared() - spaceLike_.squaredNorm();
+ return getTimeSquared() - spaceLike_.getSquaredNorm();
}
template <typename TTimeType, typename TSpaceVecType>
@@ -47,17 +47,17 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
bool FourVector<TTimeType, TSpaceVecType>::isTimelike() const {
- return getTimeSquared() < spaceLike_.squaredNorm();
+ return getTimeSquared() < spaceLike_.getSquaredNorm();
}
template <typename TTimeType, typename TSpaceVecType>
bool FourVector<TTimeType, TSpaceVecType>::isSpacelike() const {
- return getTimeSquared() > spaceLike_.squaredNorm();
+ return getTimeSquared() > spaceLike_.getSquaredNorm();
}
template <typename TTimeType, typename TSpaceVecType>
FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator+=(
- const FourVector& b) {
+ FourVector const& b) {
timeLike_ += b.timeLike_;
spaceLike_ += b.spaceLike_;
@@ -66,7 +66,7 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator-=(
- const FourVector& b) {
+ FourVector const& b) {
timeLike_ -= b.timeLike_;
spaceLike_ -= b.spaceLike_;
return *this;
@@ -74,7 +74,7 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator*=(
- const double b) {
+ double const b) {
timeLike_ *= b;
spaceLike_ *= b;
return *this;
@@ -82,22 +82,22 @@ namespace corsika {
template <typename TTimeType, typename TSpaceVecType>
FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator/=(
- const double b) {
+ double const b) {
timeLike_ /= b;
- spaceLike_.GetComponents() /= b;
+ spaceLike_.getComponents() /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::operator/(
- const double b) {
+ double const b) {
*this /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
typename FourVector<TTimeType, TSpaceVecType>::norm_type
- FourVector<TTimeType, TSpaceVecType>::operator*(const FourVector& b) {
+ FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/geometry/Helix.inl b/corsika/detail/framework/geometry/Helix.inl
index 1f9d3f94f1b815a3f2fe0420951fd4a7160de430..09cdd456cdf6548a24e702c99eef15dceb0807f2 100644
--- a/corsika/detail/framework/geometry/Helix.inl
+++ b/corsika/detail/framework/geometry/Helix.inl
@@ -30,11 +30,11 @@ namespace corsika {
}
LengthType Helix::getArcLength(TimeType t1, TimeType t2) const {
- return (vPar_ + vPerp_).norm() * (t2 - t1);
+ return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
TimeType Helix::getTimeFromArclength(LengthType l) const {
- return l / (vPar_ + vPerp_).norm();
+ return l / (vPar_ + vPerp_).getNorm();
}
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Line.inl b/corsika/detail/framework/geometry/Line.inl
index 60e78b622a43c37cf69a75c5b1841029cd472d9b..13f0292d1a812c36c8b39ec49069fcef1f011978 100644
--- a/corsika/detail/framework/geometry/Line.inl
+++ b/corsika/detail/framework/geometry/Line.inl
@@ -16,20 +16,22 @@
namespace corsika {
- Point Line::GetPosition(TimeType t) const { return r0 + v0 * t; }
+ Point Line::getPosition(TimeType t) const { return start_point_ + velocity_ * t; }
- Point Line::PositionFromArclength(LengthType l) const {
- return r0 + v0.normalized() * l;
+ Point Line::getPositionFromArclength(LengthType l) const {
+ return start_point_ + velocity_.normalized() * l;
}
- LengthType Line::ArcLength(TimeType t1, TimeType t2) const {
- return v0.norm() * (t2 - t1);
+ LengthType Line::getArcLength(TimeType t1, TimeType t2) const {
+ return velocity_.getNorm() * (t2 - t1);
}
- TimeType Line::TimeFromArclength(LengthType t) const { return t / v0.norm(); }
+ TimeType Line::getTimeFromArclength(LengthType t) const {
+ return t / velocity_.getNorm();
+ }
- const Point& Line::GetR0() const { return r0; }
+ Point const& Line::getStartPoint() const { return start_point_; }
- const Line::VelocityVec& Line::GetV0() const { return v0; }
+ Line::VelocityVec const& Line::getVelocity() const { return velocity_; }
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Plane.inl b/corsika/detail/framework/geometry/Plane.inl
index 470a598114ba21d888fe83762c22c1bb7555b20a..af0160fd6f50b0d2989802f56d1113b95b76ce19 100644
--- a/corsika/detail/framework/geometry/Plane.inl
+++ b/corsika/detail/framework/geometry/Plane.inl
@@ -16,16 +16,16 @@
namespace corsika {
- inline bool Plane::IsAbove(Point const& vP) const {
- return fNormal.dot(vP - fCenter) > LengthType::zero();
+ inline bool Plane::isAbove(Point const& vP) const {
+ return normal_.dot(vP - center_) > LengthType::zero();
}
- inline LengthType Plane::DistanceTo(corsika::Point const& vP) const {
- return (fNormal * (vP - fCenter).dot(fNormal)).norm();
+ inline LengthType Plane::getDistanceTo(corsika::Point const& vP) const {
+ return (normal_ * (vP - center_).dot(normal_)).norm();
}
- inline Point const& Plane::GetCenter() const { return fCenter; }
+ inline Point const& Plane::getCenter() const { return center_; }
- inline Plane::DimLessVec const& Plane::GetNormal() const { return fNormal; }
+ inline Plane::DimLessVec const& Plane::getNormal() const { return normal_; }
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Point.inl b/corsika/detail/framework/geometry/Point.inl
index b2c01105b2c152943f18ed698529c71697079d4b..9a24699f6de981cf9d30a13961bcb92a781902e9 100644
--- a/corsika/detail/framework/geometry/Point.inl
+++ b/corsika/detail/framework/geometry/Point.inl
@@ -17,66 +17,73 @@
namespace corsika {
+ auto Point::getCoordinates() const { return BaseVector<length_d>::getQuantityVector(); }
- // TODO: this should be private or protected, we don NOT want to expose numbers
- // without reference to outside:
- auto Point::GetCoordinates() const
- {
- return BaseVector<length_d>::qVector;
+ inline LengthType Point::getX(CoordinateSystemPtr cs) const {
+ if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ return BaseVector<length_d>::getQuantityVector().getX();
+ } else {
+ return QuantityVector<length_d>(
+ getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ BaseVector<length_d>::getQuantityVector().eigenVector_)
+ .getX();
}
+ }
- auto Point::GetX() const
- {
- return BaseVector<length_d>::qVector.GetX();
+ inline LengthType Point::getY(CoordinateSystemPtr cs) const {
+ if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ return BaseVector<length_d>::getQuantityVector().getY();
+ } else {
+ return QuantityVector<length_d>(
+ getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ BaseVector<length_d>::getQuantityVector().eigenVector_)
+ .getY();
}
+ }
- auto Point::GetY() const
- {
- return BaseVector<length_d>::qVector.GetY();
+ inline LengthType Point::getZ(CoordinateSystemPtr cs) const {
+ if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
+ return BaseVector<length_d>::getQuantityVector().getZ();
+ } else {
+ return QuantityVector<length_d>(
+ getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
+ BaseVector<length_d>::getQuantityVector().eigenVector_)
+ .getZ();
}
+ }
- auto Point::GetZ() const
- {
- return BaseVector<length_d>::qVector.GetZ();
+ /// this always returns a QuantityVector as triple
+ auto Point::getCoordinates(CoordinateSystemPtr pCS) const {
+ if (*pCS == *BaseVector<length_d>::getCoordinateSystem()) {
+ return BaseVector<length_d>::getQuantityVector();
+ } else {
+ return QuantityVector<length_d>(
+ getTransformation(BaseVector<length_d>::getCoordinateSystem(), pCS) *
+ BaseVector<length_d>::getQuantityVector().eigenVector_);
}
+ }
- /// this always returns a QuantityVector as triple
- auto Point::GetCoordinates( CoordinateSystem const& pCS) const
- {
- if (&pCS == BaseVector<length_d>::cs) {
- return BaseVector<length_d>::qVector;
- } else {
- return QuantityVector<length_d>(
- getTransformation(*BaseVector<length_d>::cs, pCS) *
- BaseVector<length_d>::qVector.eVector);
- }
- }
+ /*!
+ * transforms the Point into another CoordinateSystem by changing its
+ * coordinates interally
+ */
+ void Point::rebase(CoordinateSystemPtr pCS) {
+ BaseVector<length_d>::quantityVector() = getCoordinates(pCS);
+ BaseVector<length_d>::setCoordinateSystem(pCS);
+ }
- /*!
- * transforms the Point into another CoordinateSystem by changing its
- * coordinates interally
- */
- void Point::rebase(CoordinateSystem const& pCS)
- {
- BaseVector<length_d>::qVector = GetCoordinates(pCS);
- BaseVector<length_d>::cs = &pCS;
- }
-
- Point Point::operator+(Vector<length_d> const& pVec) const
- {
- return Point(*BaseVector<length_d>::cs,
- GetCoordinates() + pVec.GetComponents(*BaseVector<length_d>::cs));
- }
-
- /*!
- * returns the distance Vector between two points
- */
- Vector<length_d> Point::operator-( Point const& pB) const
- {
- auto& cs = *BaseVector<length_d>::cs;
- return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
- }
+ Point Point::operator+(Vector<length_d> const& pVec) const {
+ return Point(BaseVector<length_d>::getCoordinateSystem(),
+ getCoordinates() +
+ pVec.getComponents(BaseVector<length_d>::getCoordinateSystem()));
+ }
+ /*!
+ * returns the distance Vector between two points
+ */
+ Vector<length_d> Point::operator-(Point const& pB) const {
+ CoordinateSystemPtr cs = BaseVector<length_d>::getCoordinateSystem();
+ return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
+ }
} // namespace corsika
-
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index 9ac8513a93f7334ecbbe2140bf77cf89dcc630b9..2b206d3e33daf773428fce9b6fb51980a385fdda 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -18,135 +18,134 @@
namespace corsika {
- template <typename dim>
- inline auto QuantityVector<dim>::operator[](size_t index) const {
- return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
- }
+ template <typename dim>
+ inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::operator[](
+ size_t index) const {
+ return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::GetX() const
- {
- return (*this)[0];
- }
+ template <typename dim>
+ inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getX() const {
+ return (*this)[0];
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::GetY() const
- {
- return (*this)[1];
- }
+ template <typename dim>
+ inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getY() const {
+ return (*this)[1];
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::GetZ() const
- {
- return (*this)[2];
- }
+ template <typename dim>
+ inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getZ() const {
+ return (*this)[2];
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::norm() const
- {
- return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
- }
+ template <typename dim>
+ inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getNorm()
+ const {
+ return quantity_type(phys::units::detail::magnitude_tag, eigenVector_.norm());
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::squaredNorm() const
- {
- using QuantitySquared =
- decltype(std::declval<Quantity>() * std::declval<Quantity>());
- return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
- }
+ template <typename dim>
+ inline auto QuantityVector<dim>::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 auto QuantityVector<dim>::operator+(QuantityVector<dim> const& pQVec) const
- {
- return QuantityVector<dim>(eVector + pQVec.eVector);
- }
+ template <typename dim>
+ inline QuantityVector<dim> QuantityVector<dim>::operator+(
+ QuantityVector<dim> const& pQVec) const {
+ return QuantityVector<dim>(eigenVector_ + pQVec.eigenVector_);
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::operator-(QuantityVector<dim> const& pQVec) const
- {
- return QuantityVector<dim>(eVector - pQVec.eVector);
- }
+ template <typename dim>
+ inline QuantityVector<dim> QuantityVector<dim>::operator-(
+ QuantityVector<dim> const& pQVec) const {
+ return QuantityVector<dim>(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>;
-
- if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
- // a "Quantity" anymore
- {
- return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
- } else {
- return QuantityVector<typename ResQuantity::dimension_type>(eVector *
- p.magnitude());
- }
- }
+ 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 dim>
- template <typename ScalarDim>
- inline auto QuantityVector<dim>::operator/(phys::units::quantity<ScalarDim, double> const p) const
+ if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
+ // a "quantity_type" anymore
{
- return (*this) * (1 / p);
+ return QuantityVector<phys::units::dimensionless_d>(eigenVector_ * p.magnitude());
+ } else {
+ return QuantityVector<typename ResQuantity::dimension_type>(eigenVector_ *
+ p.magnitude());
}
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::operator*(double const p) const
- {
- return QuantityVector<dim>(eVector * p);
- }
+ template <typename dim>
+ template <typename ScalarDim>
+ inline auto QuantityVector<dim>::operator/(
+ phys::units::quantity<ScalarDim, double> const p) const {
+ return (*this) * (1 / p);
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::operator/(double const p) const
- {
- return QuantityVector<dim>(eVector / p);
- }
+ template <typename dim>
+ inline auto QuantityVector<dim>::operator*(double const p) const {
+ return QuantityVector<dim>(eigenVector_ * p);
+ }
- template <typename dim>
- inline auto& QuantityVector<dim>::operator/=(double const p) {
- eVector /= p;
- return *this;
- }
+ template <typename dim>
+ inline auto QuantityVector<dim>::operator/(double const p) const {
+ return QuantityVector<dim>(eigenVector_ / p);
+ }
- template <typename dim>
- inline auto& QuantityVector<dim>::operator*=(double const p) {
- eVector *= p;
- return *this;
- }
+ template <typename dim>
+ inline auto& QuantityVector<dim>::operator/=(double const p) {
+ eigenVector_ /= p;
+ return *this;
+ }
- template <typename dim>
- inline auto& QuantityVector<dim>::operator+=(QuantityVector<dim> const& pQVec) {
- eVector += pQVec.eVector;
- return *this;
- }
+ template <typename dim>
+ inline auto& QuantityVector<dim>::operator*=(double const p) {
+ eigenVector_ *= p;
+ return *this;
+ }
- template <typename dim>
- inline auto& QuantityVector<dim>::operator-=(QuantityVector<dim> const& pQVec) {
- eVector -= pQVec.eVector;
- return *this;
- }
+ template <typename dim>
+ inline auto& QuantityVector<dim>::operator+=(QuantityVector<dim> const& pQVec) {
+ eigenVector_ += pQVec.eigenVector_;
+ return *this;
+ }
- template <typename dim>
- inline auto& QuantityVector<dim>::operator-() const {
- return QuantityVector<dim>(-eVector);
- }
+ template <typename dim>
+ inline auto& QuantityVector<dim>::operator-=(QuantityVector<dim> const& pQVec) {
+ eigenVector_ -= pQVec.eigenVector_;
+ return *this;
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::normalized() const { return QuantityVector<dim>(eVector.normalized()); }
+ template <typename dim>
+ inline auto& QuantityVector<dim>::operator-() const {
+ return QuantityVector<dim>(-eigenVector_);
+ }
- template <typename dim>
- inline auto QuantityVector<dim>::operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
+ template <typename dim>
+ inline auto QuantityVector<dim>::normalized() const {
+ return QuantityVector<dim>(eigenVector_.normalized());
+ }
+ template <typename dim>
+ inline auto QuantityVector<dim>::operator==(QuantityVector<dim> const& p) const {
+ return eigenVector_ == p.eigenVector_;
+ }
template <typename dim>
- inline auto& operator<<(std::ostream& os, corsika::QuantityVector<dim> qv) {
- using Quantity = phys::units::quantity<dim, double>;
+ inline std::ostream& operator<<(std::ostream& os, corsika::QuantityVector<dim> qv) {
+ using quantity_type = phys::units::quantity<dim, double>;
- os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
+ os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
+ << qv.eigenVector_(2) << ") "
<< phys::units::to_unit_symbol<dim, double>(
- Quantity(phys::units::detail::magnitude_tag, 1));
+ quantity_type(phys::units::detail::magnitude_tag, 1));
return os;
}
} // namespace corsika
-
diff --git a/corsika/detail/framework/geometry/Sphere.inl b/corsika/detail/framework/geometry/Sphere.inl
index 1f3735655bb68f8c34d50cd7a37701c972479eb6..9a6688eb506fab5159bdc009d2b72512544888f4 100644
--- a/corsika/detail/framework/geometry/Sphere.inl
+++ b/corsika/detail/framework/geometry/Sphere.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.
@@ -12,17 +10,15 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/geometry/Volume.hpp>
namespace corsika {
- //! returns true if the Point p is within the sphere
- bool Sphere::Contains(Point const& p) const {
- return fRadius * fRadius > (fCenter - p).squaredNorm();
+ bool Sphere::isInside(Point const& p) const {
+ return radius_ * radius_ > (center_ - p).getSquaredNorm();
}
- const Point& Sphere::GetCenter() const { return fCenter; }
+ Point const& Sphere::getCenter() const { return center_; }
- LengthType Sphere::GetRadius() const { return fRadius; }
+ LengthType Sphere::getRadius() const { return radius_; }
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Trajectory.inl b/corsika/detail/framework/geometry/Trajectory.inl
index 6937ab8e30520299f54a0362d178385e94526224..89b02ad95a7a0425d14379d4b81285cae94f9afc 100644
--- a/corsika/detail/framework/geometry/Trajectory.inl
+++ b/corsika/detail/framework/geometry/Trajectory.inl
@@ -16,37 +16,37 @@
namespace corsika {
- template <typename T>
- Point Trajectory<T>::GetPosition(double u) const {
- return T::GetPosition(fTimeLength * u);
+ template <typename TType>
+ Point Trajectory<TType>::getPosition(double u) const {
+ return TType::getPosition(timeLength_ * u);
}
- template <typename T>
- TimeType Trajectory<T>::GetDuration() const {
- return fTimeLength;
+ template <typename TType>
+ TimeType Trajectory<TType>::getDuration() const {
+ return timeLength_;
}
- template <typename T>
- LengthType Trajectory<T>::GetLength() const {
- return GetDistance(fTimeLength);
+ template <typename TType>
+ LengthType Trajectory<TType>::getLength() const {
+ return getDistance(timeLength_);
}
- template <typename T>
- LengthType Trajectory<T>::GetDistance(TimeType t) const {
- assert(t <= fTimeLength);
+ template <typename TType>
+ LengthType Trajectory<TType>::getDistance(TimeType t) const {
+ assert(t <= timeLength_);
assert(t >= 0 * second);
- return T::ArcLength(0 * second, t);
+ return TType::getArcLength(0 * second, t);
}
- template <typename T>
- void Trajectory<T>::LimitEndTo(LengthType limit) {
- fTimeLength = T::TimeFromArclength(limit);
+ template <typename TType>
+ void Trajectory<TType>::getLimitEndTo(LengthType limit) {
+ timeLength_ = TType::getTimeFromArclength(limit);
}
- template <typename T>
- auto Trajectory<T>::NormalizedDirection() const {
- static_assert(std::is_same_v<T, corsika::Line>);
- return T::GetV0().normalized();
+ template <typename TType>
+ auto Trajectory<TType>::getNormalizedDirection() const {
+ static_assert(std::is_same_v<TType, corsika::Line>);
+ return TType::getVelocity().normalized();
}
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Vector.inl b/corsika/detail/framework/geometry/Vector.inl
index 419bce29833971ded68824ba242763690392ba18..1df45dfc6f640f45bc39560072e290f76102b67e 100644
--- a/corsika/detail/framework/geometry/Vector.inl
+++ b/corsika/detail/framework/geometry/Vector.inl
@@ -14,182 +14,188 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-
namespace corsika {
-
- template <typename dim>
- auto Vector<dim>::GetComponents() const
- {
- return BaseVector<dim>::qVector;
- }
-
-
- template <typename dim>
- auto Vector<dim>::GetComponents(CoordinateSystem const& pCS) const
- {
- if (&pCS == BaseVector<dim>::cs) {
- return BaseVector<dim>::qVector;
- } else {
- return QuantityVector<dim>(
- getTransformation(*BaseVector<dim>::cs, pCS).linear() *
- BaseVector<dim>::qVector.eVector);
- }
- }
-
- template <typename dim>
- void Vector<dim>::rebase(CoordinateSystem const& pCS)
- {
- BaseVector<dim>::qVector = GetComponents(pCS);
- BaseVector<dim>::cs = &pCS;
- }
-
- template <typename dim>
- auto Vector<dim>::norm() const
- {
- return BaseVector<dim>::qVector.norm();
- }
-
- template <typename dim>
- auto Vector<dim>::GetNorm() const
- {
- return BaseVector<dim>::qVector.norm();
- }
-
- template <typename dim>
- auto Vector<dim>::squaredNorm() const
- {
- return BaseVector<dim>::qVector.squaredNorm();
- }
-
- template <typename dim>
- auto Vector<dim>::GetSquaredNorm() const
- {
- return BaseVector<dim>::qVector.squaredNorm();
- }
-
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::parallelProjectionOnto(Vector<dim2> const& pVec,
- CoordinateSystem const& pCS) const
- {
- auto const ourCompVec = GetComponents(pCS);
- auto const otherCompVec = pVec.GetComponents(pCS);
- auto const& a = ourCompVec.eVector;
- auto const& b = otherCompVec.eVector;
-
- return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
- }
-
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::parallelProjectionOnto(Vector<dim2> const& pVec) const
- {
- return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
- }
-
- template <typename dim>
- auto Vector<dim>::operator+(Vector<dim> const& pVec) const
- {
- auto const components =
- GetComponents(*BaseVector<dim>::cs) + pVec.GetComponents(*BaseVector<dim>::cs);
- return Vector<dim>(*BaseVector<dim>::cs, components);
- }
-
- template <typename dim>
- auto Vector<dim>::operator-(Vector<dim> const& pVec) const
- {
- auto const components = GetComponents() - pVec.GetComponents(*BaseVector<dim>::cs);
- return Vector<dim>(*BaseVector<dim>::cs, components);
- }
-
- template <typename dim>
- auto& Vector<dim>::operator*=(double const p)
- {
- BaseVector<dim>::qVector *= 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>;
-
- return Vector<ProdDim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
- }
-
- template <typename dim>
- template <typename ScalarDim>
- auto Vector<dim>::operator/(phys::units::quantity<ScalarDim, double> const p) const
- {
- return (*this) * (1 / p);
- }
-
- template <typename dim>
- auto Vector<dim>::operator*(double const p) const
- {
- return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
- }
-
- template <typename dim>
- auto Vector<dim>::operator/(double const p) const
- {
- return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
- }
-
- template <typename dim>
- auto& Vector<dim>::operator+=(Vector<dim> const& pVec)
- {
- BaseVector<dim>::qVector += pVec.GetComponents(*BaseVector<dim>::cs);
- return *this;
- }
-
- template <typename dim>
- auto& Vector<dim>::operator-=(Vector<dim> const& pVec)
- {
- BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
- return *this;
- }
-
- template <typename dim>
- auto& Vector<dim>::operator-() const
- {
- return Vector<dim>(*BaseVector<dim>::cs, -BaseVector<dim>::qVector);
- }
-
- template <typename dim>
- auto Vector<dim>::normalized() const
- {
- return (*this) * (1 / norm());
- }
-
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::cross(Vector<dim2> pV) const
- {
- auto const c1 = GetComponents().eVector;
- auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
- auto const bareResult = c1.cross(c2);
-
- using ProdDim = phys::units::detail::product_d<dim, dim2>;
- return Vector<ProdDim>(*BaseVector<dim>::cs, bareResult);
- }
-
- template <typename dim>
- template <typename dim2>
- auto Vector<dim>::dot(Vector<dim2> pV) const
- {
- auto const c1 = GetComponents().eVector;
- auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
- auto const bareResult = c1.dot(c2);
-
- using ProdDim = phys::units::detail::product_d<dim, dim2>;
-
- return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
- bareResult);
- }
-
+ template <typename dim>
+ QuantityVector<dim> Vector<dim>::getComponents() const {
+ return BaseVector<dim>::getQuantityVector();
+ }
+
+ template <typename dim>
+ QuantityVector<dim> Vector<dim>::getComponents(CoordinateSystemPtr pCS) const {
+ if (*pCS == *BaseVector<dim>::getCoordinateSystem()) { // FIXME
+ return BaseVector<dim>::getQuantityVector();
+ } else {
+ return QuantityVector<dim>(
+ getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
+ BaseVector<dim>::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];
+ } else {
+ return QuantityVector<dim>(
+ getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
+ BaseVector<dim>::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];
+ } else {
+ return QuantityVector<dim>(
+ getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
+ BaseVector<dim>::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];
+ } else {
+ return QuantityVector<dim>(
+ getTransformation(BaseVector<dim>::getCoordinateSystem(), pCS).linear() *
+ BaseVector<dim>::getQuantityVector().eigenVector_)[2];
+ }
+ }
+
+ template <typename dim>
+ void Vector<dim>::rebase(CoordinateSystemPtr pCS) {
+ BaseVector<dim>::quantityVector() = getComponents(pCS);
+ BaseVector<dim>::setCoordinateSystem(pCS);
+ }
+
+ template <typename dim>
+ inline typename Vector<dim>::quantity_type Vector<dim>::getNorm() const {
+ return BaseVector<dim>::getQuantityVector().getNorm();
+ }
+
+ template <typename dim>
+ auto Vector<dim>::getSquaredNorm() const {
+ return BaseVector<dim>::getQuantityVector().getSquaredNorm();
+ }
+
+ template <typename dim>
+ template <typename dim2>
+ auto Vector<dim>::getParallelProjectionOnto(Vector<dim2> const& pVec,
+ CoordinateSystemPtr 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())));
+ }
+
+ template <typename dim>
+ template <typename dim2>
+ auto Vector<dim>::getParallelProjectionOnto(Vector<dim2> const& pVec) const {
+ return getParallelProjectionOnto<dim2>(pVec, BaseVector<dim>::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 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 dim>
+ auto& Vector<dim>::operator*=(double const p) {
+ BaseVector<dim>::quantityVector() *= 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>;
+
+ return Vector<ProdDim>(BaseVector<dim>::getCoordinateSystem(),
+ BaseVector<dim>::getQuantityVector() * p);
+ }
+
+ template <typename dim>
+ template <typename ScalarDim>
+ auto Vector<dim>::operator/(phys::units::quantity<ScalarDim, 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 dim>
+ auto Vector<dim>::operator/(double const p) const {
+ return Vector<dim>(BaseVector<dim>::getCoordinateSystem(),
+ BaseVector<dim>::quantityVector() / p);
+ }
+
+ template <typename dim>
+ auto& Vector<dim>::operator+=(Vector<dim> const& pVec) {
+ BaseVector<dim>::quantityVector() +=
+ pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
+ return *this;
+ }
+
+ template <typename dim>
+ auto& Vector<dim>::operator-=(Vector<dim> const& pVec) {
+ BaseVector<dim>::quantityVector() -=
+ pVec.getComponents(BaseVector<dim>::getCoordinateSystem());
+ return *this;
+ }
+
+ template <typename dim>
+ auto& Vector<dim>::operator-() const {
+ return Vector<dim>(BaseVector<dim>::getCoordinateSystem(),
+ -BaseVector<dim>::quantityVector());
+ }
+
+ template <typename dim>
+ auto Vector<dim>::normalized() const {
+ return (*this) * (1 / getNorm());
+ }
+
+ template <typename dim>
+ template <typename dim2>
+ auto Vector<dim>::cross(Vector<dim2> pV) const {
+ auto const c1 = getComponents().eigenVector_;
+ auto const c2 = pV.getComponents(BaseVector<dim>::getCoordinateSystem()).eigenVector_;
+ auto const bareResult = c1.cross(c2);
+
+ using ProdDim = phys::units::detail::product_d<dim, dim2>;
+ return Vector<ProdDim>(BaseVector<dim>::getCoordinateSystem(), bareResult);
+ }
+
+ template <typename dim>
+ template <typename dim2>
+ auto Vector<dim>::dot(Vector<dim2> pV) const {
+ auto const c1 = getComponents().eigenVector_;
+ auto const c2 = pV.getComponents(BaseVector<dim>::getCoordinateSystem()).eigenVector_;
+ auto const bareResult = c1.dot(c2);
+
+ using ProdDim = phys::units::detail::product_d<dim, dim2>;
+
+ return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
+ bareResult);
+ }
} // namespace corsika
-
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index 338114b9ff260230cee8bdc786d739d48ddf6dac..cef0ce97cfba55aaa6a8e5f6f108b937122e8883 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -17,21 +17,21 @@
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/logging/Logging.h>
+#include <corsika/framework/logging/Logging.hpp>
// using namespace corsika::units::si;
namespace corsika {
COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
- const HEPMassType massTarget)
- : originalCS_{Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()}
- , rotatedCS_{originalCS_.RotateToZ(Pprojectile.GetSpaceLikeComponents())} {
- auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
- auto const pProjNormSquared = pProjectile.squaredNorm();
+ HEPMassType const massTarget)
+ : originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
+ , rotatedCS_{originalCS_->rotateToZ(Pprojectile.getSpaceLikeComponents())} {
+ auto const pProjectile = Pprojectile.getSpaceLikeComponents();
+ auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
- auto const eProjectile = Pprojectile.GetTimeLikeComponent();
+ auto const eProjectile = Pprojectile.getTimeLikeComponent();
auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
auto const s =
massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
@@ -42,15 +42,14 @@ namespace corsika {
setBoost(coshEta, sinhEta);
- C8LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta, coshEta,
- boost_.determinant() - 1);
+ CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+ coshEta, boost_.determinant() - 1);
}
- COMBoost::COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
- units::si::HEPEnergyType mass)
- : originalCS_{momentum.GetCoordinateSystem()}
- , rotatedCS_{originalCS_.RotateToZ(momentum)} {
- auto const squaredNorm = momentum.squaredNorm();
+ COMBoost::COMBoost(Vector<hepmomentum_d> const& momentum, HEPEnergyType mass)
+ : originalCS_{momentum.getCoordinateSystem()}
+ , rotatedCS_{originalCS_->rotateToZ(momentum)} {
+ auto const squaredNorm = momentum.getSquaredNorm();
auto const norm = sqrt(squaredNorm);
auto const sinhEta = -norm / mass;
auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
@@ -58,13 +57,12 @@ namespace corsika {
}
template <typename FourVector>
- FourVector COMBoost::toCoM(const FourVector& p) const {
- using namespace corsika::units::si;
- auto pComponents = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
- Eigen::Vector3d eVecRotated = pComponents.eVector;
+ FourVector COMBoost::toCoM(FourVector const& p) const {
+ auto pComponents = p.getSpaceLikeComponents().getComponents(rotatedCS_);
+ Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
Eigen::Vector2d lab;
- lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ lab << (p.getTimeLikeComponent() * (1 / 1_GeV)),
(eVecRotated(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = boost_ * lab;
@@ -72,39 +70,37 @@ namespace corsika {
eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
- return FourVector(E_CoM,
- corsika::geometry::Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
+ return FourVector(E_CoM, Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
}
template <typename FourVector>
- FourVector COMBoost::fromCoM(const FourVector& p) const {
- using namespace corsika::units::si;
- auto pCM = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
- auto const Ecm = p.GetTimeLikeComponent();
+ FourVector COMBoost::fromCoM(FourVector const& p) const {
+ auto pCM = p.getSpaceLikeComponents().getComponents(rotatedCS_);
+ auto const Ecm = p.getTimeLikeComponent();
Eigen::Vector2d com;
- com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
+ com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
- C8LOG_TRACE(
+ CORSIKA_LOG_TRACE(
"COMBoost::fromCoM Ecm={} GeV"
" pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
- Ecm / 1_GeV, pCM / 1_GeV, pCM.norm() / 1_GeV, p.GetNorm() / 1_GeV);
+ Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV, p.getNorm() / 1_GeV);
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
- pCM.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
+ pCM.eigenVector()(2) = boostedZ(1) * (1_GeV).magnitude();
- geometry::Vector<typename decltype(pCM)::dimension> pLab{rotatedCS_, pCM};
+ Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
pLab.rebase(originalCS_);
FourVector f(E_lab, pLab);
- C8LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
- " plab={} GeV (norm={} GeV) "
- " GeV), invariant mass = {}",
- E_lab / 1_GeV, f.GetNorm() / 1_GeV, pLab.GetComponents(),
- pLab.norm() / 1_GeV);
+ CORSIKA_LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
+ " plab={} GeV (norm={} GeV) "
+ " GeV), invariant mass = {}",
+ E_lab / 1_GeV, f.getNorm() / 1_GeV, pLab.getComponents(),
+ pLab.getNorm() / 1_GeV);
return f;
}
@@ -114,6 +110,6 @@ namespace corsika {
inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
}
- CoordinateSystem const& COMBoost::GetRotatedCS() const { return rotatedCS_; }
+ CoordinateSystemPtr COMBoost::getRotatedCS() const { return rotatedCS_; }
} // namespace corsika
diff --git a/corsika/framework/core/PhysicalUnits.hpp b/corsika/framework/core/PhysicalUnits.hpp
index 6de5bd4d0352de65c6859c7009a550b9d8d6e330..6388b2188815f2d679aac2289eef803ccca0c9d3 100644
--- a/corsika/framework/core/PhysicalUnits.hpp
+++ b/corsika/framework/core/PhysicalUnits.hpp
@@ -74,6 +74,7 @@ namespace corsika::units::si {
using sigma_d = phys::units::area_d;
/// add the unit-types
+ using DimensionlessType = phys::units::quantity<phys::units::dimensionless_d, double>;
using LengthType = phys::units::quantity<phys::units::length_d, double>;
using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
using SpeedType = phys::units::quantity<phys::units::speed_d, double>;
@@ -95,7 +96,7 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
template <typename DimFrom, typename DimTo>
- auto constexpr ConversionFactorHEPToSI() {
+ auto constexpr conversion_factor_HEP_to_SI() {
static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
DimFrom::dim7 == 0,
@@ -120,7 +121,7 @@ namespace corsika::units::si {
}
template <typename DimFrom>
- auto constexpr ConversionFactorSIToHEP() {
+ auto constexpr conversion_factor_SI_to_HEP() {
static_assert(DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
DimFrom::dim7 == 0 && DimFrom::dim8 == 0,
"must be pure L, M, T type");
@@ -138,14 +139,18 @@ namespace corsika::units::si {
}
template <typename DimTo, typename DimFrom>
- auto constexpr ConvertHEPToSI(quantity<DimFrom> q) {
- return ConversionFactorHEPToSI<DimFrom, DimTo>() * q;
+ auto constexpr convert_HEP_to_SI(quantity<DimFrom> q) {
+ return conversion_factor_HEP_to_SI<DimFrom, DimTo>() * q;
}
template <typename DimFrom>
- auto constexpr ConvertSIToHEP(quantity<DimFrom> q) {
- return ConversionFactorSIToHEP<DimFrom>() * q;
+ auto constexpr convert_SI_to_HEP(quantity<DimFrom> q) {
+ return conversion_factor_SI_to_HEP<DimFrom>() * q;
}
+
+ template <typename T>
+ bool operator==(DimensionlessType a, T b){ return a.magnitude() == b; }
+
} // end namespace corsika::units::si
/**
diff --git a/corsika/framework/geometry/BaseVector.hpp b/corsika/framework/geometry/BaseVector.hpp
index 3545ec7e14f93203df544b626b6795c734d451b0..21819df294d18cc6f52f6feadb3acda470ab32a4 100644
--- a/corsika/framework/geometry/BaseVector.hpp
+++ b/corsika/framework/geometry/BaseVector.hpp
@@ -11,33 +11,40 @@ n/*
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <memory>
+
namespace corsika {
/*!
- * Common base class for Vector and Point. Currently it does basically nothing.
- */
- /*
- * FIXME Many potential issues:
- * 1. does this class really need to be templated ?
- * 2. copy constructor, assignment operator not implemented
- * 3. this member pointer is quite scary...
+ * Common base class for Vector and Point.
+ *
+ * This holds a QuantityVector and a CoordinateSystem
+ *
*/
- template <typename dim>
+ template <typename TDimension>
class BaseVector {
public:
- /*
- * FIXME Why to copy pQVector twice?
- */
- BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
- : qVector(pQVector)
- , cs(&pCS) {}
+ BaseVector(CoordinateSystemPtr pCS, QuantityVector<TDimension> const& pQVector)
+ : quantityVector_(pQVector)
+ , cs_(pCS) {}
- auto const& GetCoordinateSystem() const;
+ BaseVector() = delete;
+ BaseVector(BaseVector const&) = default;
+ BaseVector(BaseVector&& a) = default;
+ BaseVector& operator=(BaseVector const&) = default;
+ ~BaseVector() = default;
+
+ CoordinateSystemPtr getCoordinateSystem() const;
+ void setCoordinateSystem(CoordinateSystemPtr cs) { cs_ = cs; }
protected:
- QuantityVector<dim> qVector;
- CoordinateSystem const* cs;
+ QuantityVector<TDimension> const& getQuantityVector() const;
+ QuantityVector<TDimension>& quantityVector();
+
+ private:
+ QuantityVector<TDimension> quantityVector_;
+ CoordinateSystemPtr cs_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/CoordinateSystem.hpp b/corsika/framework/geometry/CoordinateSystem.hpp
index 9b6b6ac0297f95f55fb30f44da76c8419a1c439f..6cd5c4813941698766ebcf870ae70d5a2a300969 100644
--- a/corsika/framework/geometry/CoordinateSystem.hpp
+++ b/corsika/framework/geometry/CoordinateSystem.hpp
@@ -13,65 +13,97 @@ n/*
#include <Eigen/Dense>
#include <stdexcept>
-
-/*
- * FIXME Review this global typedef.
- */
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
+#include <memory>
namespace corsika {
- class RootCoordinateSystem;
+ typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
+ typedef Eigen::Translation<double, 3> EigenTranslation;
template <typename T>
- class Vector;
+ class Vector; // fwd decl
+
+ class CoordinateSystem; // fwd decl
+ /**
+ * To refer to CoordinateSystems, only the CoordinateSystemPtr must be used.
+ */
+ using CoordinateSystemPtr = std::shared_ptr<CoordinateSystem const>;
+
+ class RootCoordinateSystem; // fwd decl
+ static CoordinateSystemPtr get_root_CoordinateSystem(); // fwd decl
+
+ /**
+ * A class to store the reference 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.
+ *
+ * Only the \sa RootCoordinateSystem can be created as a singleton
+ * as global main reference point.
+ *
+ * Thus, new CoordinateSystems can only be created using
+ * transformation: \sa rotateToZ, \sa rotate, \sa translateAndRotate
+ * below.
+ */
class CoordinateSystem {
- CoordinateSystem const* reference = nullptr;
- EigenTransform transf;
-
- CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
- : reference(&reference)
- , transf(transf) {}
+ /**
+ * Constructor only from referenceCS, given the transformation matrix transf
+ */
+ CoordinateSystem(CoordinateSystemPtr referenceCS, EigenTransform const& transf)
+ : referenceCS_(referenceCS)
+ , transf_(transf) {}
+ /**
+ * for creating the root CS
+ */
CoordinateSystem()
- : // for creating the root CS
- transf(EigenTransform::Identity()) {}
+ : referenceCS_(nullptr)
+ , transf_(EigenTransform::Identity()) {}
public:
- // FIXME missing test for self assignment
- inline CoordinateSystem& operator=(const CoordinateSystem& pCS);
+ // default resource allocation
+ CoordinateSystem(CoordinateSystem const&) = default;
+ CoordinateSystem(CoordinateSystem&&) = default;
+ CoordinateSystem& operator=(CoordinateSystem const& pCS) = default;
+ ~CoordinateSystem() = default;
- inline CoordinateSystem translate(QuantityVector<length_d> vector) const;
+ inline CoordinateSystemPtr translate(QuantityVector<length_d> vector) const;
/**
* creates a new CS in which vVec points in direction of the new z-axis
*/
template <typename TDim>
- auto RotateToZ(Vector<TDim> vVec) const;
+ CoordinateSystemPtr rotateToZ(Vector<TDim> vVec) const;
template <typename TDim>
- auto rotate(QuantityVector<TDim> axis, double angle) const;
+ CoordinateSystemPtr rotate(QuantityVector<TDim> axis, double angle) const;
template <typename TDim>
- auto translateAndRotate(QuantityVector<length_d> translation,
- QuantityVector<TDim> axis, double angle);
+ CoordinateSystemPtr translateAndRotate(QuantityVector<length_d> translation,
+ QuantityVector<TDim> axis, double angle);
- inline CoordinateSystem const* GetReference() const;
+ inline CoordinateSystemPtr getReferenceCS() const;
- inline const EigenTransform& GetTransform() const;
+ inline EigenTransform const& getTransform() const;
+
+ inline bool operator==(CoordinateSystem const&) const;
+ inline bool operator!=(CoordinateSystem const&) const;
protected:
- static CoordinateSystem CreateCS() { return CoordinateSystem(); }
+ static CoordinateSystem createCS() { return CoordinateSystem(); }
+
+ friend CoordinateSystemPtr get_root_CoordinateSystem(); /// this is the only way to
+ /// create ONE unique root CS
- friend corsika::RootCoordinateSystem; /// this is the only class that can
- /// create ONE unique root CS
+ private:
+ std::shared_ptr<CoordinateSystem const> referenceCS_;
+ EigenTransform transf_;
};
- EigenTransform getTransformation(CoordinateSystem const& c1,
- CoordinateSystem const& c2);
+ EigenTransform getTransformation(CoordinateSystemPtr c1, CoordinateSystemPtr c2);
} // namespace corsika
diff --git a/corsika/framework/geometry/FourVector.hpp b/corsika/framework/geometry/FourVector.hpp
index 8290493050e7e971b3556c442b35c4980429392b..55fefc2fb5a30e74772648864de84cabc1607796 100644
--- a/corsika/framework/geometry/FourVector.hpp
+++ b/corsika/framework/geometry/FourVector.hpp
@@ -44,7 +44,7 @@ namespace corsika {
public:
using space_vec_type = typename std::decay<TSpaceVecType>::type;
- using space_type = typename space_vec_type::Quantity;
+ using space_type = typename space_vec_type::quantity_type;
using time_type = typename std::decay<TTimeType>::type;
//! check the types and the physical units here:
@@ -119,11 +119,11 @@ namespace corsika {
bool isSpacelike() const;
FourVector& operator+=(FourVector const&);
- FourVector& operator-=(const FourVector&);
- FourVector& operator*=(const double);
- FourVector& operator/=(const double);
- FourVector& operator/(const double);
-
+ FourVector& operator-=(FourVector const&);
+ FourVector& operator*=(double const);
+ FourVector& operator/=(double const);
+ FourVector& operator/(double const);
+
/**
Scalar product of two FourVectors
@@ -132,10 +132,10 @@ namespace corsika {
for this. You cannot mix different conventions due to
unit-checking.
*/
- norm_type operator*(const FourVector& b);
-
+ norm_type operator*(FourVector const& b);
+
/** @} */
-
+
protected:
//! the data members
TTimeType timeLike_;
@@ -169,7 +169,7 @@ namespace corsika {
}
friend FourVector<time_type, space_vec_type> operator*(FourVector const& a,
- const double b) {
+ double const b) {
return FourVector<time_type, space_vec_type>(a.timeLike_ * b, a.spaceLike_ * b);
}
diff --git a/corsika/framework/geometry/Helix.hpp b/corsika/framework/geometry/Helix.hpp
index 55b071586bd8c7d62d4bb5819143d519074eacde..05ca40be07ac3bbd1a6ad4b7515c1717b5565c72 100644
--- a/corsika/framework/geometry/Helix.hpp
+++ b/corsika/framework/geometry/Helix.hpp
@@ -42,12 +42,12 @@ namespace corsika {
, vPar_(pvPar)
, vPerp_(pvPerp)
, uPerp_(vPerp_.cross(vPar_.normalized()))
- , radius_(pvPar.norm() / abs(pOmegaC)) {}
+ , radius_(pvPar.getNorm() / abs(pOmegaC)) {}
inline LengthType getRadius() const;
inline Point getPosition(TimeType t) const;
-
+
inline Point getPositionFromArclength(LengthType l) const;
inline LengthType getArcLength(TimeType t1, TimeType t2) const;
@@ -59,7 +59,7 @@ namespace corsika {
///! "cylinder" on which the helix rotates
FrequencyType omegaC_; ///! speed of angular rotation
VelocityVec vPar_; ///! speed along direction of "cylinder"
- VelocityVec vPerp_, uPerp_;
+ VelocityVec vPerp_, uPerp_;
LengthType radius_;
};
diff --git a/corsika/framework/geometry/Volume.hpp b/corsika/framework/geometry/IVolume.hpp
similarity index 80%
rename from corsika/framework/geometry/Volume.hpp
rename to corsika/framework/geometry/IVolume.hpp
index 0b21df8b1749dd2e09eb00c0f0464f5c3e49939c..4983b3de16272bc3c8b70f3f8a566ce1df6b6828 100644
--- a/corsika/framework/geometry/Volume.hpp
+++ b/corsika/framework/geometry/IVolume.hpp
@@ -12,13 +12,13 @@ n/*
namespace corsika {
- class Volume {
+ class IVolume {
public:
//! returns true if the Point p is within the volume
- virtual bool Contains(Point const& p) const = 0;
+ virtual bool isInside(Point const& p) const = 0;
- virtual ~Volume() = default;
+ virtual ~IVolume() = default;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Line.hpp b/corsika/framework/geometry/Line.hpp
index 6392b4609f82cf74abcb0e9df183666aa9b9085b..2c27a0bc06fc394aa156558af70bc298027eb3d1 100644
--- a/corsika/framework/geometry/Line.hpp
+++ b/corsika/framework/geometry/Line.hpp
@@ -26,27 +26,29 @@ namespace corsika {
class Line {
+ ///! \toto move this to PhysicalUnits
using VelocityVec = Vector<SpeedType::dimension_type>;
- Point const r0;
- VelocityVec const v0;
-
public:
Line(Point const& pR0, VelocityVec const& pV0)
- : r0(pR0)
- , v0(pV0) {}
+ : start_point_(pR0)
+ , velocity_(pV0) {}
+
+ inline Point getPosition(TimeType t) const;
- inline Point GetPosition(TimeType t) const;
+ inline Point getPositionFromArclength(LengthType l) const;
- inline Point PositionFromArclength(LengthType l) const;
+ inline LengthType getArcLength(TimeType t1, TimeType t2) const;
- inline LengthType ArcLength(TimeType t1, TimeType t2) const;
+ inline TimeType getTimeFromArclength(LengthType t) const;
- inline TimeType TimeFromArclength(LengthType t) const;
+ inline const Point& getStartPoint() const;
- inline const Point& GetR0() const;
+ inline const VelocityVec& getVelocity() const;
- inline const VelocityVec& GetV0() const;
+ private:
+ Point const start_point_;
+ VelocityVec const velocity_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Plane.hpp b/corsika/framework/geometry/Plane.hpp
index d66ad195f59589f67fefe060970d30b794785ad7..928b04ebeda61f7f6f98dc73a1e94022d9ad3df6 100644
--- a/corsika/framework/geometry/Plane.hpp
+++ b/corsika/framework/geometry/Plane.hpp
@@ -16,23 +16,25 @@ namespace corsika {
class Plane {
+ ///! \todo move to PhysicalUnits
using DimLessVec = Vector<dimensionless_d>;
- Point const fCenter;
- DimLessVec const fNormal;
-
public:
Plane(Point const& vCenter, DimLessVec const& vNormal)
- : fCenter(vCenter)
- , fNormal(vNormal.normalized()) {}
+ : center_(vCenter)
+ , normal_(vNormal.normalized()) {}
+
+ bool isAbove(Point const& vP) const;
- bool IsAbove(Point const& vP) const;
+ LengthType getDistanceTo(corsika::Point const& vP) const;
- LengthType DistanceTo(corsika::Point const& vP) const;
+ Point const& getCenter() const;
- Point const& GetCenter() const;
+ DimLessVec const& getNormal() const;
- DimLessVec const& GetNormal() const;
+ public:
+ Point const center_;
+ DimLessVec const normal_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Point.hpp b/corsika/framework/geometry/Point.hpp
index 0fd1303cd09d48f4413f2f06921aec5fc29af461..06c302e1b5de46172f9f73e0a11d205fa335e2a4 100644
--- a/corsika/framework/geometry/Point.hpp
+++ b/corsika/framework/geometry/Point.hpp
@@ -22,30 +22,29 @@ namespace corsika {
class Point : public BaseVector<length_d> {
public:
- Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
+ Point(CoordinateSystemPtr pCS, QuantityVector<length_d> const& pQVector)
: BaseVector<length_d>(pCS, pQVector) {}
- Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
+ Point(CoordinateSystemPtr cs, LengthType x, LengthType y, LengthType z)
: BaseVector<length_d>(cs, {x, y, z}) {}
- // TODO: this should be private or protected, we don NOT want to expose numbers
- // without reference to outside:
- inline auto GetCoordinates() const;
-
- inline auto GetX() const;
-
- inline auto GetY() const;
-
- inline auto GetZ() const;
+ /** \todo TODO: this should be private or protected, we don NOT want to expose numbers
+ * without reference to outside:
+ */
+ inline auto getCoordinates() const;
/// this always returns a QuantityVector as triple
- inline auto GetCoordinates(CoordinateSystem const& pCS) const;
+ inline auto getCoordinates(CoordinateSystemPtr pCS) const;
+
+ inline LengthType getX(CoordinateSystemPtr pCS) const;
+ inline LengthType getY(CoordinateSystemPtr pCS) const;
+ inline LengthType getZ(CoordinateSystemPtr pCS) const;
/*!
* transforms the Point into another CoordinateSystem by changing its
* coordinates interally
*/
- inline void rebase(CoordinateSystem const& pCS);
+ inline void rebase(CoordinateSystemPtr 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 b363104484446178a7ea1f6c60f74b1adf3fd99e..e1f56d380f9663b36ab0b04a6f740646a0668957 100644
--- a/corsika/framework/geometry/QuantityVector.hpp
+++ b/corsika/framework/geometry/QuantityVector.hpp
@@ -11,56 +11,63 @@ n/*
#include <Eigen/Dense>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <iostream>
+#include <ostream>
#include <utility>
namespace corsika {
+ class CoordinateSystem; // fwd decl
+ class Point;
+ template <typename T>
+ class Vector;
+
/*!
* A QuantityVector is a three-component container based on Eigen::Vector3d
* with a phys::units::si::dimension. Arithmethic operators are defined that
* propagate the dimensions by dimensional analysis.
+ *
+ * \todo review QuantityVector: can this be a protected-only
+ * inheritance to other objects? We don't want to expose access to
+ * low-level Eigen objects anywhere.
*/
- template <typename dim>
+ template <typename TDimension>
class QuantityVector {
public:
- using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
- // corresponding to the dimension
+ using quantity_type =
+ phys::units::quantity<TDimension, double>; //< the phys::units::quantity
+ // corresponding to the dimension
- public:
- Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
+ QuantityVector(Eigen::Vector3d pBareVector)
+ : eigenVector_(pBareVector) {}
- typedef dim dimension; //!< should be a phys::units::dimension
+ public:
+ typedef TDimension dimension_type; //!< should be a phys::units::dimension
- QuantityVector(Quantity a, Quantity b, Quantity c)
- : eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
+ QuantityVector(quantity_type a, quantity_type b, quantity_type c)
+ : eigenVector_{a.magnitude(), b.magnitude(), c.magnitude()} {}
QuantityVector(double a, double b, double c)
- : eVector{a, b, c} {
+ : eigenVector_{a, b, c} {
static_assert(
- std::is_same_v<dim, phys::units::dimensionless_d>,
+ std::is_same_v<TDimension, phys::units::dimensionless_d>,
"initialization of dimensionful QuantityVector with pure numbers not allowed!");
}
- QuantityVector(Eigen::Vector3d pBareVector)
- : eVector(pBareVector) {}
-
- auto operator[](size_t index) const;
+ quantity_type operator[](size_t 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_; }
- auto GetX() const;
+ quantity_type getNorm() const;
- auto GetY() const;
+ auto getSquaredNorm() const;
- auto GetZ() const;
+ QuantityVector operator+(QuantityVector<TDimension> const& pQVec) const;
- auto norm() const;
-
- auto squaredNorm() const;
-
- auto operator+(QuantityVector<dim> const& pQVec) const;
-
- auto operator-(QuantityVector<dim> const& pQVec) const;
+ QuantityVector operator-(QuantityVector<TDimension> const& pQVec) const;
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const;
@@ -76,23 +83,35 @@ namespace corsika {
auto& operator*=(double const p);
- auto& operator+=(QuantityVector<dim> const& pQVec);
+ auto& operator+=(QuantityVector<TDimension> const& pQVec);
- auto& operator-=(QuantityVector<dim> const& pQVec);
+ auto& operator-=(QuantityVector<TDimension> const& pQVec);
auto& operator-() const;
auto normalized() const;
- auto operator==(QuantityVector<dim> const& p) const;
+ auto operator==(QuantityVector<TDimension> const& p) const;
+
+ 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);
+
+ protected:
+ Eigen::Vector3d
+ eigenVector_; //!< the actual container where the raw numbers are stored
};
/*
- * FIXME free function operators not implemented.
+ * streaming operator
*/
- template <typename dim>
- auto& operator<<(std::ostream& os, corsika::QuantityVector<dim> qv);
+ template <typename TDimension>
+ inline std::ostream& operator<<(std::ostream& os,
+ corsika::QuantityVector<TDimension> qv);
} // namespace corsika
diff --git a/corsika/framework/geometry/RootCoordinateSystem.hpp b/corsika/framework/geometry/RootCoordinateSystem.hpp
index 64422fa6f5186ce46e60efcc85fb30e67aa521dc..cc64e5ff9e76a9deffe9303031155827eb928433 100644
--- a/corsika/framework/geometry/RootCoordinateSystem.hpp
+++ b/corsika/framework/geometry/RootCoordinateSystem.hpp
@@ -9,29 +9,23 @@ n/*
#pragma once
#include <corsika/framework/geometry/CoordinateSystem.hpp>
-#include <corsika/framework/utility/Singleton.hpp>
-/*!
- * This is the only way to get a root-coordinate system, and it is a
- * singleton. All other CoordinateSystems must be relative to the
- * RootCoordinateSystem
- */
+#include <memory>
namespace corsika {
- class RootCoordinateSystem : public corsika::Singleton<RootCoordinateSystem> {
-
- friend class corsika::Singleton<RootCoordinateSystem>;
-
- protected:
- RootCoordinateSystem() {}
-
- public:
- corsika::CoordinateSystem& GetRootCoordinateSystem() { return fRootCS; }
- const corsika::CoordinateSystem& GetRootCoordinateSystem() const { return fRootCS; }
-
- private:
- corsika::CoordinateSystem fRootCS; // THIS IS IT
- };
+ /*!
+ * Singleton factory function to produce the root CoordinateSystem
+ *
+ * This is the only way to get a root-coordinate system, and it is a
+ * singleton. All other CoordinateSystems must be relative to the
+ * RootCoordinateSystem
+ */
+
+ static std::shared_ptr<CoordinateSystem const> get_root_CoordinateSystem() {
+ static std::shared_ptr<CoordinateSystem const> rootCS(
+ new CoordinateSystem); // THIS IS IT
+ return rootCS;
+ }
} // namespace corsika
diff --git a/corsika/framework/geometry/Sphere.hpp b/corsika/framework/geometry/Sphere.hpp
index 187f2816cb036b7eade98d8558604b1a3a8e20cf..d5f4c44a863f3d495dec35e404ea77669f9780cb 100644
--- a/corsika/framework/geometry/Sphere.hpp
+++ b/corsika/framework/geometry/Sphere.hpp
@@ -10,25 +10,27 @@ n/*
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Volume.hpp>
+#include <corsika/framework/geometry/IVolume.hpp>
namespace corsika {
- class Sphere : public Volume {
- Point const fCenter;
- LengthType const fRadius;
+ class Sphere : public IVolume {
public:
Sphere(Point const& pCenter, LengthType const pRadius)
- : fCenter(pCenter)
- , fRadius(pRadius) {}
+ : center_(pCenter)
+ , radius_(pRadius) {}
//! returns true if the Point p is within the sphere
- inline bool Contains(Point const& p) const override;
+ inline bool isInside(Point const& p) const override;
- inline const Point& GetCenter() const;
+ inline const Point& getCenter() const;
- inline LengthType GetRadius() const;
+ inline LengthType getRadius() const;
+
+ private:
+ Point const center_;
+ LengthType const radius_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Trajectory.hpp b/corsika/framework/geometry/Trajectory.hpp
index a8869a952fdf40da466debe302a738ecc58e96c6..0c26f9dcd4666269cf2bbcd6865336de4667aebe 100644
--- a/corsika/framework/geometry/Trajectory.hpp
+++ b/corsika/framework/geometry/Trajectory.hpp
@@ -14,34 +14,32 @@
namespace corsika {
- template <typename T>
- class Trajectory : public T {
-
- TimeType fTimeLength;
+ template <typename TType>
+ class Trajectory : public TType {
public:
- using T::ArcLength;
- using T::GetPosition;
+ using TType::getArcLength;
+ using TType::getPosition;
+
+ Trajectory(TType const& theT, TimeType timeLength)
+ : TType(theT)
+ , timeLength_(timeLength) {}
- Trajectory(T const& theT, TimeType timeLength)
- : T(theT)
- , fTimeLength(timeLength) {}
+ Point getPosition(double u) const;
- /*Point GetPosition(TimeType t) const {
- return fTraj.GetPosition(t + fTStart);
- }*/
+ TimeType getDuration() const;
- Point GetPosition(double u) const;
+ LengthType getLength() const;
- TimeType GetDuration() const;
+ LengthType getDistance(TimeType t) const;
- LengthType GetLength() const;
+ void getLimitEndTo(LengthType limit);
- LengthType GetDistance(TimeType t) const;
+ auto getNormalizedDirection() const;
- void LimitEndTo(LengthType limit);
+ private:
+ TimeType timeLength_;
- auto NormalizedDirection() const;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Vector.hpp b/corsika/framework/geometry/Vector.hpp
index 7ca7d905e8cbf91d7316c58e035e4b54ae4f3093..96a154562c560e0d9372bdea6c70cad145acd762 100644
--- a/corsika/framework/geometry/Vector.hpp
+++ b/corsika/framework/geometry/Vector.hpp
@@ -12,63 +12,66 @@ n/*
#include <corsika/framework/geometry/BaseVector.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
-/*!
- * 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).
- *
- * When transforming coordinate systems, a Vector is subject to the rotational
- * part only and invariant under translations.
- */
-
namespace corsika {
- template <typename dim>
- class Vector : public BaseVector<dim> {
+ /*!
+ * 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).
+ *
+ * When transforming coordinate systems, a Vector is subject to the rotational
+ * part only and invariant under translations.
+ */
+
+ template <typename TDimension>
+ class Vector : public BaseVector<TDimension> {
public:
- using Quantity = phys::units::quantity<dim, double>;
+ using quantity_type = phys::units::quantity<TDimension, double>;
public:
- Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
- : BaseVector<dim>(pCS, pQVector) {}
+ Vector(CoordinateSystemPtr pCS, QuantityVector<TDimension> pQVector)
+ : BaseVector<TDimension>(pCS, pQVector) {}
- Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z)
- : BaseVector<dim>(cs, QuantityVector<dim>(x, y, z)) {}
+ Vector(CoordinateSystemPtr cs, quantity_type x, quantity_type y, quantity_type z)
+ : BaseVector<TDimension>(cs, QuantityVector<TDimension>(x, y, z)) {}
/*!
* 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
*/
- auto GetComponents() const;
+ QuantityVector<TDimension> getComponents() const;
/*!
* returns a QuantityVector with the components given in an arbitrary
* CoordinateSystem
*/
- auto GetComponents(CoordinateSystem const& pCS) const;
+ QuantityVector<TDimension> getComponents(CoordinateSystemPtr pCS) const;
+
+ inline quantity_type getX(CoordinateSystemPtr pCS) const;
+ inline quantity_type getY(CoordinateSystemPtr pCS) const;
+ inline quantity_type getZ(CoordinateSystemPtr pCS) const;
/*!
* transforms the Vector into another CoordinateSystem by changing
* its components internally
*/
- void rebase(CoordinateSystem const& pCS);
+ inline void rebase(CoordinateSystemPtr pCS);
/*!
* returns the norm/length of the Vector. Before using this method,
* think about whether squaredNorm() might be cheaper for your computation.
*/
- auto norm() const;
-
- auto GetNorm() const;
+ inline quantity_type getNorm() const;
/*!
* returns the squared norm of the Vector. Before using this method,
* think about whether norm() might be cheaper for your computation.
*/
- auto squaredNorm() const;
-
- auto GetSquaredNorm() const;
+ inline auto 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
@@ -76,40 +79,40 @@ namespace corsika {
* \vec{v}_{\parallel} = \frac{\vec{v}_1 \cdot \vec{v}_2}{\vec{v}_2^2} \vec{v}_2
* \f]
*/
- template <typename dim2>
- auto parallelProjectionOnto(Vector<dim2> const& pVec,
- CoordinateSystem const& pCS) const;
- template <typename dim2>
- auto parallelProjectionOnto(Vector<dim2> const& pVec) const;
+ template <typename TDimension2>
+ auto getParallelProjectionOnto(Vector<TDimension2> const& pVec,
+ CoordinateSystemPtr pCS) const;
+ template <typename TDimension2>
+ auto getParallelProjectionOnto(Vector<TDimension2> const& pVec) const;
- auto operator+(Vector<dim> const& pVec) const;
+ Vector operator+(Vector<TDimension> const& pVec) const;
- auto operator-(Vector<dim> const& pVec) const;
+ Vector operator-(Vector<TDimension> const& pVec) const;
auto& operator*=(double const p);
- template <typename ScalarDim>
- auto operator*(phys::units::quantity<ScalarDim, 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;
+ template <typename TScalarDim>
+ auto operator/(phys::units::quantity<TScalarDim, double> const p) const;
auto operator*(double const p) const;
auto operator/(double const p) const;
- auto& operator+=(Vector<dim> const& pVec);
+ auto& operator+=(Vector<TDimension> const& pVec);
- auto& operator-=(Vector<dim> const& pVec);
+ auto& operator-=(Vector<TDimension> const& pVec);
auto& operator-() const;
auto normalized() const;
- template <typename dim2>
- auto cross(Vector<dim2> pV) const;
+ template <typename TDimension2>
+ auto cross(Vector<TDimension2> pV) const;
- template <typename dim2>
- auto dot(Vector<dim2> pV) const;
+ template <typename TDimension2>
+ auto dot(Vector<TDimension2> pV) const;
};
} // namespace corsika
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index 08932465b1f54792126179d26564dfd04ccff8ef..85f11d53ddce40434d2bd0493ed3b1b7da5ec329 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -8,7 +8,6 @@
#pragma once
-
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
@@ -24,38 +23,33 @@ namespace corsika {
*/
class COMBoost {
- Eigen::Matrix2d boost_, inverseBoost_;
- corsika::CoordinateSystem const &originalCS_, rotatedCS_;
-
- void setBoost(double coshEta, double sinhEta);
public:
//! construct a COMBoost given four-vector of projectile and mass of target
- COMBoost(
- const corsika::geometry::FourVector<
- corsika::units::si::HEPEnergyType,
- corsika::geometry::Vector<corsika::units::si::hepmomentum_d>>& Pprojectile,
- const corsika::units::si::HEPEnergyType massTarget);
+ COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
+ HEPEnergyType const massTarget);
//! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
- COMBoost(Vector<units::si::hepmomentum_d> const& momentum,
- units::si::HEPEnergyType mass);
+ COMBoost(Vector<hepmomentum_d> const& momentum, HEPEnergyType mass);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
template <typename FourVector>
- FourVector toCoM(const FourVector& p) const;
-
+ FourVector toCoM(FourVector const& p) const;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
template <typename FourVector>
- FourVector fromCoM(const FourVector& p) const;
+ FourVector fromCoM(FourVector const& p) const;
+
+ CoordinateSystemPtr getRotatedCS() const;
+ protected:
+ void setBoost(double coshEta, double sinhEta);
+
+ private:
+ Eigen::Matrix2d boost_, inverseBoost_;
+ CoordinateSystemPtr originalCS_, rotatedCS_;
- CoordinateSystem const& GetRotatedCS() const;
};
} // namespace corsika
#include <corsika/detail/framework/utility/COMBoost.inl>
-
-
-
diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp
index 5eb68de94b954456eb4e173fd8956492e3308bb7..0458c257f3befd14cf24945f4966afe8d19b1e3d 100644
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -24,15 +24,14 @@ int main() {
std::cout << "geometry_example" << std::endl;
// define the root coordinate system
- corsika::CoordinateSystem& root =
- corsika::RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ corsika::CoordinateSystemPtr root = get_root_CoordinateSystem();
// another CS defined by a translation relative to the root CS
- CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
+ CoordinateSystemPtr cs2 = root->translate({0_m, 0_m, 1_m});
// rotations are possible, too; parameters are axis vector and angle
- CoordinateSystem cs3 =
- root.rotate(QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
+ CoordinateSystemPtr 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
@@ -41,30 +40,30 @@ int main() {
Vector<length_d> const diff =
p2 -
p1; // the distance between the points, basically the translation vector given above
- auto const norm = diff.squaredNorm(); // squared length with the right dimension
+ auto const norm = diff.getSquaredNorm(); // squared length with the right dimension
// print the components of the vector as given in the different CS
- std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl;
- std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2)
+ std::cout << "p2-p1 components in root: " << diff.getComponents(root) << std::endl;
+ std::cout << "p2-p1 components in cs2: " << diff.getComponents(cs2)
<< std::endl; // by definition invariant under translations
- std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3)
+ std::cout << "p2-p1 components in cs3: " << diff.getComponents(cs3)
<< std::endl; // but not under rotations
std::cout << "p2-p1 norm^2: " << norm << std::endl;
assert(norm == 1 * meter * meter);
Sphere s(p1, 10_m); // define a sphere around a point with a radius
- std::cout << "p1 inside s: " << s.Contains(p2) << std::endl;
- assert(s.Contains(p2) == 1);
+ std::cout << "p1 inside s: " << s.isInside(p2) << std::endl;
+ assert(s.isInside(p2) == 1);
Sphere s2(p1, 3_um); // another sphere
- std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
- assert(s2.Contains(p2) == 0);
+ std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl;
+ assert(s2.isInside(p2) == 0);
// let's try parallel projections:
auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});
auto const v2 = Vector<length_d>(root, {1_m, 0_m, 0_m});
- auto const v3 = v1.parallelProjectionOnto(v2);
+ auto const v3 = v1.getParallelProjectionOnto(v2);
// cross product
auto const cross =
@@ -72,10 +71,10 @@ int main() {
// if a CS is not given as parameter for getComponents(), the components
// in the "home" CS are returned
- std::cout << "v1: " << v1.GetComponents() << std::endl;
- std::cout << "v2: " << v2.GetComponents() << std::endl;
- std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl;
- std::cout << "normalized cross product of v1 x v2" << cross.GetComponents()
+ std::cout << "v1: " << v1.getComponents() << std::endl;
+ std::cout << "v2: " << v2.getComponents() << std::endl;
+ std::cout << "parallel projection of v1 onto v2: " << v3.getComponents() << std::endl;
+ std::cout << "normalized cross product of v1 x v2" << cross.getComponents()
<< std::endl;
return EXIT_SUCCESS;
diff --git a/examples/helix_example.cpp b/examples/helix_example.cpp
index 364b7f6039f948da8143cdd5180d469139e17695..74994dc95df59dc115a83a1575793a5a1d7edec2 100644
--- a/examples/helix_example.cpp
+++ b/examples/helix_example.cpp
@@ -20,8 +20,7 @@ using namespace corsika;
using namespace corsika::units::si;
int main() {
- corsika::CoordinateSystem& root =
- corsika::RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ corsika::CoordinateSystemPtr root = get_root_CoordinateSystem();
Point const r0(root, {0_m, 0_m, 0_m});
auto const omegaC = 2 * M_PI * 1_Hz;
@@ -39,7 +38,7 @@ int main() {
auto& positions = *arr;
for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) {
- auto const r = h.GetPosition(t).GetCoordinates();
+ auto const r = h.getPosition(t).getCoordinates();
positions[i][0] = t / 1_s;
positions[i][1] = r[0] / 1_m;
diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt
index 81bab8f2273995fa69c6586b9d5a2e78b5e90889..5b345ba5855f8ee14393eca4720e8716e9e780ee 100644
--- a/tests/framework/CMakeLists.txt
+++ b/tests/framework/CMakeLists.txt
@@ -3,7 +3,7 @@ set (test_framework_sources
# testCascade.cpp this is most important, but whole content of former Processes folder missing yet
testClassTimer.cpp
#testCombinedStack.cpp #FIXME: reenable this
- #testCOMBoost.cpp #FIXME: reenable this
+ testCOMBoost.cpp #FIXME: reenable this
#testCorsikaFenv.cpp # does not work because of use of exceptions in catch2
testFourVector.cpp
testHelix.cpp
diff --git a/tests/framework/testCOMBoost.cpp b/tests/framework/testCOMBoost.cpp
index d26222a7be0624abd0e41793486fdf81f8857e37..e0d22810a74f6876e2ccb73cfed45de4d0bfda6c 100644
--- a/tests/framework/testCOMBoost.cpp
+++ b/tests/framework/testCOMBoost.cpp
@@ -14,28 +14,141 @@
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
-#include <iostream>
-
using namespace corsika;
double constexpr absMargin = 1e-6;
-CoordinateSystem const& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
// helper function for energy-momentum
// relativistic energy
auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
- return sqrt(m * m + p.squaredNorm());
+ return sqrt(m * m + p.getSquaredNorm());
};
auto const momentum = [](HEPEnergyType E, HEPMassType m) { return sqrt(E * E - m * m); };
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
- return E * E - p.squaredNorm();
+ return E * E - p.getSquaredNorm();
};
+TEST_CASE("rotation") {
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1_GeV;
+ HEPMassType const targetMass = 1.0e300_eV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ FourVector const PprojLab(eProjectileLab, pProjectileLab);
+
+ Vector<hepenergy_d> e1(rootCS, {1_GeV, 0_GeV, 0_GeV});
+ Vector<hepenergy_d> e2(rootCS, {0_GeV, 1_GeV, 0_GeV});
+ Vector<hepenergy_d> e3(rootCS, {0_GeV, 0_GeV, 1_GeV});
+
+ // define boost to com frame
+ SECTION("pos. z-axis") {
+ COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, 1_GeV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ }
+
+ SECTION("y-axis in upper half") {
+ COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, 1_meV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ }
+
+ SECTION("x-axis in upper half") {
+ COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ }
+
+ SECTION("neg. z-axis") {
+ COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, -1_GeV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ }
+
+ SECTION("x-axis lower half") {
+ COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ }
+
+ SECTION("y-axis lower half") {
+ COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, -1_meV}}}, targetMass);
+ CoordinateSystemPtr rotCS = boost.getRotatedCS();
+
+ CHECK(e1.getX(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+ CHECK(e1.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e1.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+
+ CHECK(e2.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getY(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e2.getZ(rotCS) / 1_GeV == Approx(1).margin(absMargin));
+
+ CHECK(e3.getX(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ CHECK(e3.getY(rotCS) / 1_GeV == Approx(-1).margin(absMargin));
+ CHECK(e3.getZ(rotCS) / 1_GeV == Approx(0).margin(absMargin));
+ }
+}
+
TEST_CASE("boosts") {
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
@@ -52,7 +165,7 @@ TEST_CASE("boosts") {
HEPMassType const projectileMass = 1._GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 20_GeV, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ FourVector const PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
@@ -65,28 +178,28 @@ TEST_CASE("boosts") {
// sum of momenta in CoM, should be 0
auto const sumPCoM =
- PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
- CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+ PprojCoM.getSpaceLikeComponents() + PtargCoM.getSpaceLikeComponents();
+ CHECK(sumPCoM.getNorm() / 1_GeV == Approx(0).margin(absMargin));
// mandelstam-s should be invariant under transformation
CHECK(s(eProjectileLab + eTargetLab,
- pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ pProjectileLab.getComponents() + pTargetLab.getComponents()) /
1_GeV / 1_GeV ==
- Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
- PprojCoM.GetSpaceLikeComponents().GetComponents() +
- PtargCoM.GetSpaceLikeComponents().GetComponents()) /
+ Approx(s(PprojCoM.getTimeLikeComponent() + PtargCoM.getTimeLikeComponent(),
+ PprojCoM.getSpaceLikeComponents().getComponents() +
+ PtargCoM.getSpaceLikeComponents().getComponents()) /
1_GeV / 1_GeV));
// boost back...
auto const PprojBack = boost.fromCoM(PprojCoM);
// ...should yield original values before the boosts
- CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() ==
+ CHECK(PprojBack.getTimeLikeComponent() / PprojLab.getTimeLikeComponent() ==
Approx(1));
- CHECK(
- (PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
- PprojLab.GetSpaceLikeComponents().norm() ==
- Approx(0).margin(absMargin));
+ CHECK((PprojBack.getSpaceLikeComponents() - PprojLab.getSpaceLikeComponents())
+ .getNorm() /
+ PprojLab.getSpaceLikeComponents().getNorm() ==
+ Approx(0).margin(absMargin));
}
/*
@@ -99,7 +212,7 @@ TEST_CASE("boosts") {
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_GeV, -20_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ FourVector const PprojLab(eProjectileLab, pProjectileLab);
auto const sqrt_s_lab =
sqrt(s(eProjectileLab + targetMass, pProjectileLab.GetComponents(rootCS)));
@@ -120,8 +233,8 @@ TEST_CASE("boosts") {
// sum of momenta in CoM, should be 0
auto const sumPCoM =
- PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
- CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+ PprojCoM.getSpaceLikeComponents() + PtargCoM.getSpaceLikeComponents();
+ CHECK(sumPCoM.getNorm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
@@ -130,7 +243,7 @@ TEST_CASE("boosts") {
SECTION("Test boost along tilted axis") {
- const HEPMomentumType P0 = 1_PeV;
+ HEPMomentumType const P0 = 1_PeV;
double theta = 33.;
double phi = -10.;
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
@@ -144,7 +257,7 @@ TEST_CASE("boosts") {
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab(rootCS, {px, py, pz});
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ FourVector const PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
@@ -157,8 +270,8 @@ TEST_CASE("boosts") {
// sum of momenta in CoM, should be 0
auto const sumPCoM =
- PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
- CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
+ PprojCoM.getSpaceLikeComponents() + PtargCoM.getSpaceLikeComponents();
+ CHECK(sumPCoM.getNorm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
@@ -171,7 +284,7 @@ TEST_CASE("boosts") {
HEPMomentumType P0 = 1_ZeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -P0}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ FourVector const PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
@@ -184,8 +297,8 @@ TEST_CASE("boosts") {
// sum of momenta in CoM, should be 0
auto const sumPCoM =
- PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
- CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
+ PprojCoM.getSpaceLikeComponents() + PtargCoM.getSpaceLikeComponents();
+ CHECK(sumPCoM.getNorm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
}
diff --git a/tests/framework/testFourVector.cpp b/tests/framework/testFourVector.cpp
index b4cda441b9c9d9328d4812ad73c6a95ec3d9f571..a5982da826ce1c7d72b2a21659a67bf9e9b676eb 100644
--- a/tests/framework/testFourVector.cpp
+++ b/tests/framework/testFourVector.cpp
@@ -8,20 +8,20 @@
#include <catch2/catch.hpp>
-#include <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <cmath>
+
using namespace corsika;
TEST_CASE("four vectors") {
// this is just needed as a baseline
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
/*
Test: P2 = E2 - p2 all in [GeV]
@@ -73,7 +73,7 @@ TEST_CASE("four vectors") {
FourVector p1(E1, P1);
- const double check = 100 * 100 - 10 * 10 - 5 * 5 - 15 * 15; // for dummies...
+ double const check = 100 * 100 - 10 * 10 - 5 * 5 - 15 * 15; // for dummies...
CHECK(p1.getNormSqr() / 1_GeV / 1_GeV == Approx(check));
CHECK(p1.getNorm() / 1_GeV == Approx(sqrt(check)));
@@ -88,7 +88,7 @@ TEST_CASE("four vectors") {
TimeType T2 = 10_m / constants::c;
Vector<length_d> P2(rootCS, {10_m, 5_m, 5_m});
- const double check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
+ double const check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
FourVector p2(T2, P2);
@@ -108,8 +108,8 @@ TEST_CASE("four vectors") {
HEPEnergyType E2 = 0_GeV;
Vector<hepmomentum_d> P2(rootCS, {10_GeV, 0_GeV, 0_GeV});
- FourVector p1(E1, P1);
- const FourVector p2(E2, P2);
+ FourVector const p1(E1, P1);
+ FourVector const p2(E2, P2);
CHECK(p1.getNorm() / 1_GeV == Approx(100.));
CHECK(p2.getNorm() / 1_GeV == Approx(10.));
@@ -156,22 +156,22 @@ TEST_CASE("four vectors") {
TimeType T = 10_m / constants::c;
Vector<length_d> P(rootCS, {10_m, 5_m, 5_m});
- const TimeType T_c = 10_m / constants::c;
- const Vector<length_d> P_c(rootCS, {10_m, 5_m, 5_m});
+ TimeType const T_c = 10_m / constants::c;
+ Vector<length_d> const P_c(rootCS, {10_m, 5_m, 5_m});
/*
this does not compile, and it shoudn't!
FourVector<TimeType&, Vector<length_d>&> p0(T_c, P_c);
*/
FourVector<TimeType&, Vector<length_d>&> p1(T, P);
- FourVector<const TimeType&, const Vector<length_d>&> p2(T, P);
- FourVector<const TimeType&, const Vector<length_d>&> p3(T_c, P_c);
+ FourVector<TimeType const&, Vector<length_d> const&> p2(T, P);
+ FourVector<TimeType const&, Vector<length_d> const&> p3(T_c, P_c);
p1 *= 10;
// p2 *= 10; // this does not compile, and it shoudn't !
// p3 *= 10; // this does not compile, and it shoudn't !!
- const double check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
+ double const check = 10 * 10 - 10 * 10 - 5 * 5 - 5 * 5; // for dummies...
CHECK(p1.getNormSqr() / (1_m * 1_m) == Approx(10. * 10. * check));
CHECK(p2.getNorm() / 1_m == Approx(10 * sqrt(abs(check))));
CHECK(p3.getNorm() / 1_m == Approx(sqrt(abs(check))));
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index b3eb21e070229726e97bd8f2aeb62e76274e681b..e303e641976593afe5cf355ce4e236cace23d31c 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -22,8 +22,7 @@ using namespace corsika;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
REQUIRE(getTransformation(rootCS, rootCS).isApprox(EigenTransform::Identity()));
@@ -33,68 +32,68 @@ TEST_CASE("transformations between CoordinateSystems") {
QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
Vector<magnetic_flux_density_d> v1(rootCS, components);
- REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() ==
+ REQUIRE((p1.getCoordinates() - coordinates).getNorm().magnitude() ==
Approx(0).margin(absMargin));
- REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
+ REQUIRE((p1.getCoordinates(rootCS) - coordinates).getNorm().magnitude() ==
Approx(0).margin(absMargin));
/*
SECTION("unconnected CoordinateSystems") {
CoordinateSystem rootCS2 = CoordinateSystem::CreateRootCS();
- REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
+ REQUIRE_THROWS(CoordinateSystem::getTransformation(rootCS, rootCS2));
}*/
SECTION("translations") {
QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
- CoordinateSystem translatedCS = rootCS.translate(translationVector);
+ CoordinateSystemPtr translatedCS = rootCS->translate(translationVector);
- REQUIRE(translatedCS.GetReference() == &rootCS);
+ REQUIRE(*translatedCS->getReferenceCS() == *rootCS);
- REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() ==
+ REQUIRE((p1.getCoordinates(translatedCS) + translationVector).getNorm().magnitude() ==
Approx(0).margin(absMargin));
// Vectors are not subject to translations
- REQUIRE(
- (v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() ==
- Approx(0).margin(absMargin));
+ REQUIRE((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).norm().magnitude() ==
+ REQUIRE(((p2 - p1).getComponents() - translationVector).getNorm().magnitude() ==
Approx(0).margin(absMargin));
}
SECTION("multiple translations") {
QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
- CoordinateSystem cs2 = rootCS.translate(tv1);
+ CoordinateSystemPtr cs2 = rootCS->translate(tv1);
QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
- CoordinateSystem cs3 = rootCS.translate(tv2);
+ CoordinateSystemPtr cs3 = rootCS->translate(tv2);
QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
- CoordinateSystem cs4 = cs3.translate(tv3);
+ CoordinateSystemPtr cs4 = cs3->translate(tv3);
- REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
+ REQUIRE(*cs4->getReferenceCS()->getReferenceCS() == *rootCS);
REQUIRE(getTransformation(cs3, cs2).isApprox(
- rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
+ rootCS->translate({3_m, -5_m, 0_m})->getTransform()));
REQUIRE(getTransformation(cs2, cs3).isApprox(
- rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()));
+ rootCS->translate({-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;
- CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
- REQUIRE(rotatedCS.GetReference() == &rootCS);
+ CoordinateSystemPtr rotatedCS = rootCS->rotate(axis, angle);
+ REQUIRE(*rotatedCS->getReferenceCS() == *rootCS);
- REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() ==
+ REQUIRE(v1.getComponents(rotatedCS)[1].magnitude() ==
Approx((-1. * tesla).magnitude()));
// vector norm invariant under rotation
- REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() ==
- Approx(v1.GetComponents(rootCS).norm().magnitude()));
+ REQUIRE(v1.getComponents(rotatedCS).getNorm().magnitude() ==
+ Approx(v1.getComponents(rootCS).getNorm().magnitude()));
}
SECTION("multiple rotations") {
@@ -108,15 +107,79 @@ TEST_CASE("transformations between CoordinateSystems") {
double const angle = 90. / 180. * M_PI;
- CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle);
- CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle);
- CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle);
+ CoordinateSystemPtr rotated1 = rootCS->rotate(zAxis, angle);
+ CoordinateSystemPtr rotated2 = rotated1->rotate(yAxis, angle);
+ CoordinateSystemPtr rotated3 = rotated2->rotate(zAxis, -angle);
- CoordinateSystem combined = rootCS.rotate(xAxis, -angle);
+ CoordinateSystemPtr combined = rootCS->rotate(xAxis, -angle);
- auto comp1 = v1.GetComponents(rotated3);
- auto comp3 = v1.GetComponents(combined);
- REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
+ auto comp1 = v1.getComponents(rotated3);
+ auto comp3 = v1.getComponents(combined);
+ REQUIRE((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);
+ 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(xPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
+ CHECK(yPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
+ CHECK((zPrime.dot(v) / 1_m).magnitude() == Approx(5 * sqrt(2)));
+
+ CHECK(zPrime.getComponents(rootCS)[1].magnitude() ==
+ Approx(zPrime.getComponents(rootCS)[2].magnitude()));
+ CHECK(zPrime.getComponents(rootCS)[0].magnitude() == Approx(0));
+
+ CHECK(xPrime.getComponents(rootCS).getEigenVector().dot(
+ yPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+ CHECK(zPrime.getComponents(rootCS).getEigenVector().dot(
+ xPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+ CHECK(yPrime.getComponents(rootCS).getEigenVector().dot(
+ zPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+
+ CHECK(yPrime.getComponents(rootCS).getEigenVector().dot(
+ yPrime.getComponents(rootCS).getEigenVector()) ==
+ Approx((5_m * 5_m).magnitude()));
+ CHECK(xPrime.getComponents(rootCS).getEigenVector().dot(
+ xPrime.getComponents(rootCS).getEigenVector()) ==
+ Approx((5_m * 5_m).magnitude()));
+ CHECK(zPrime.getComponents(rootCS).getEigenVector().dot(
+ zPrime.getComponents(rootCS).getEigenVector()) ==
+ 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);
+ 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).getEigenVector().dot(
+ v.getComponents().getEigenVector()) == Approx(0));
+ CHECK(yPrime.getComponents(rootCS).getEigenVector().dot(
+ v.getComponents().getEigenVector()) == Approx(0));
+
+ CHECK(xPrime.getComponents(rootCS).getEigenVector().dot(
+ yPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+ CHECK(zPrime.getComponents(rootCS).getEigenVector().dot(
+ xPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+ CHECK(yPrime.getComponents(rootCS).getEigenVector().dot(
+ zPrime.getComponents(rootCS).getEigenVector()) == Approx(0));
+
+ CHECK(yPrime.getComponents(rootCS).getEigenVector().dot(
+ yPrime.getComponents(rootCS).getEigenVector()) ==
+ Approx((5_m * 5_m).magnitude()));
+ CHECK(xPrime.getComponents(rootCS).getEigenVector().dot(
+ xPrime.getComponents(rootCS).getEigenVector()) ==
+ Approx((5_m * 5_m).magnitude()));
+ CHECK(zPrime.getComponents(rootCS).getEigenVector().dot(
+ zPrime.getComponents(rootCS).getEigenVector()) ==
+ Approx((5_m * 5_m).magnitude()));
}
SECTION("RotateToZ positive") {
@@ -179,28 +242,26 @@ TEST_CASE("transformations between CoordinateSystems") {
}
TEST_CASE("Sphere") {
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
- SECTION("GetCenter") {
- CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
+ SECTION("getCenter") {
+ CHECK((sphere.getCenter().getCoordinates(rootCS) -
QuantityVector<length_d>(0_m, 3_m, 4_m))
- .norm()
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
- CHECK(sphere.GetRadius() / 5_m == Approx(1));
+ CHECK(sphere.getRadius() / 5_m == Approx(1));
}
- SECTION("Contains") {
- REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
- REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
+ 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})));
}
}
TEST_CASE("Trajectories") {
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
@@ -209,26 +270,28 @@ TEST_CASE("Trajectories") {
Line const line(r0, v0);
CHECK(
- (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
- .norm()
+ (line.getPosition(2_s).getCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
- CHECK((line.PositionFromArclength(4_m).GetCoordinates() -
+ CHECK((line.getPositionFromArclength(4_m).getCoordinates() -
QuantityVector<length_d>(4_m, 0_m, 0_m))
- .norm()
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
- CHECK((line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
- .norm()
+ CHECK((line.getPosition(7_s) -
+ line.getPositionFromArclength(line.getArcLength(0_s, 7_s)))
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
- LineTrajectory base(line, t);
- CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
+ Trajectory<Line> base(line, t);
+ CHECK(line.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates());
- CHECK((base.GetVelocity(0).normalized().GetComponents(rootCS) -
+ CHECK(base.getArcLength(1_s, 2_s) / 1_m == Approx(3));
+
+ CHECK((base.getNormalizedDirection().getComponents(rootCS) -
QuantityVector<dimensionless_d>{1, 0, 0})
- .eVector.norm() == Approx(0).margin(absMargin));
+ .getNorm() == Approx(0).margin(absMargin));
}
-
}
diff --git a/tests/framework/testHelix.cpp b/tests/framework/testHelix.cpp
index 74338e8ab27f2a80c44afd868f912a777f904fd1..5acb065ab96e7599bf1f2f7cb6110e992bd9c77f 100644
--- a/tests/framework/testHelix.cpp
+++ b/tests/framework/testHelix.cpp
@@ -21,8 +21,7 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("Helix class") {
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ const CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Helix") {
@@ -37,20 +36,20 @@ TEST_CASE("Helix class") {
Helix const helix(r0, omegaC, vPar, vPerp);
- CHECK((helix.getPosition(1_s).GetCoordinates() -
+ CHECK((helix.getPosition(1_s).getCoordinates() -
QuantityVector<length_d>(0_m, 0_m, 4_m))
- .norm()
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
- CHECK((helix.getPosition(0.25_s).GetCoordinates() -
+ CHECK((helix.getPosition(0.25_s).getCoordinates() -
QuantityVector<length_d>(-3_m / (2 * M_PI), -3_m / (2 * M_PI), 1_m))
- .norm()
+ .getNorm()
.magnitude() == Approx(0).margin(absMargin));
- CHECK(
- (helix.getPosition(7_s) - helix.getPositionFromArclength(helix.getArcLength(0_s, 7_s)))
- .norm()
- .magnitude() == Approx(0).margin(absMargin));
+ CHECK((helix.getPosition(7_s) -
+ helix.getPositionFromArclength(helix.getArcLength(0_s, 7_s)))
+ .getNorm()
+ .magnitude() == Approx(0).margin(absMargin));
/*
// we have to consider this, if we need it
@@ -61,6 +60,4 @@ TEST_CASE("Helix class") {
CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
*/
}
-
}
-
diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp
index 26f4f918ffe969aeced6af35b196e693bf1e72b9..f87af1c2a2ee56da960d9ea530bcb9b328cc4826 100644
--- a/tests/framework/testUnits.cpp
+++ b/tests/framework/testUnits.cpp
@@ -19,17 +19,17 @@ using namespace corsika::units::si;
TEST_CASE("PhysicalUnits", "[Units]") {
SECTION("Consistency") {
- REQUIRE(1_m / 1_m == Approx(1));
- // REQUIRE_FALSE( 1_m/1_s == 1 ); // static assert
+ CHECK(1_m / 1_m == Approx(1));
+ // CHECK_FALSE( 1_m/1_s == 1 ); // static assert
}
SECTION("Constructors") {
[[maybe_unused]] auto E1 = 10_GeV;
- REQUIRE(E1 == 10_GeV);
+ CHECK(E1 == 10_GeV);
LengthType l1 = 10_nm;
[[maybe_unused]] auto l2 = l1;
- REQUIRE(l2 == l1);
+ CHECK(l2 == l1);
LengthType arr0[5];
arr0[0] = 5_m;
@@ -41,54 +41,54 @@ TEST_CASE("PhysicalUnits", "[Units]") {
[[maybe_unused]] std::array<HEPEnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
auto p1 = 10_s * newton;
- REQUIRE(p1 == 10_s * newton);
+ CHECK(p1 == 10_s * newton);
}
SECTION("Powers in literal units") {
- REQUIRE(1_s / 1_ds == Approx(1e1));
- REQUIRE(1_m / 1_cm == Approx(1e2));
- REQUIRE(1_m / 1_mm == Approx(1e3));
- REQUIRE(1_V / 1_uV == Approx(1e6));
- REQUIRE(1_s / 1_ns == Approx(1e9));
- REQUIRE(1_eV / 1_peV == Approx(1e12));
- REQUIRE(1_A / 1_fA == Approx(1e15));
- REQUIRE(1_mol / 1_amol == Approx(1e18));
- REQUIRE(1_K / 1_zK == Approx(1e21));
- REQUIRE(1_K / 1_yK == Approx(1e24));
- REQUIRE(1_b / 1_mb == Approx(1e3));
-
- REQUIRE(1_A / 1_hA == Approx(1e-2));
- REQUIRE(1_m / 1_km == Approx(1e-3));
- REQUIRE(1_m / 1_Mm == Approx(1e-6));
- REQUIRE(1_V / 1_GV == Approx(1e-9));
- REQUIRE(1_s / 1_Ts == Approx(1e-12));
- REQUIRE(1_eV / 1_PeV == Approx(1e-15));
- REQUIRE(1_A / 1_EA == Approx(1e-18));
- REQUIRE(1_K / 1_ZK == Approx(1e-21));
- REQUIRE(1_mol / 1_Ymol == Approx(1e-24));
-
- REQUIRE(std::min(1_A, 2_A) == 1_A);
+ CHECK(1_s / 1_ds == Approx(1e1));
+ CHECK(1_m / 1_cm == Approx(1e2));
+ CHECK(1_m / 1_mm == Approx(1e3));
+ CHECK(1_V / 1_uV == Approx(1e6));
+ CHECK(1_s / 1_ns == Approx(1e9));
+ CHECK(1_eV / 1_peV == Approx(1e12));
+ CHECK(1_A / 1_fA == Approx(1e15));
+ CHECK(1_mol / 1_amol == Approx(1e18));
+ CHECK(1_K / 1_zK == Approx(1e21));
+ CHECK(1_K / 1_yK == Approx(1e24));
+ CHECK(1_b / 1_mb == Approx(1e3));
+
+ CHECK(1_A / 1_hA == Approx(1e-2));
+ CHECK(1_m / 1_km == Approx(1e-3));
+ CHECK(1_m / 1_Mm == Approx(1e-6));
+ CHECK(1_V / 1_GV == Approx(1e-9));
+ CHECK(1_s / 1_Ts == Approx(1e-12));
+ CHECK(1_eV / 1_PeV == Approx(1e-15));
+ CHECK(1_A / 1_EA == Approx(1e-18));
+ CHECK(1_K / 1_ZK == Approx(1e-21));
+ CHECK(1_mol / 1_Ymol == Approx(1e-24));
+
+ CHECK(std::min(1_A, 2_A) == 1_A);
}
SECTION("Powers and units") {
- REQUIRE(1 * ampere / 1_A == Approx(1e0));
- REQUIRE(mega * bar / bar == Approx(1e6));
+ CHECK(1 * ampere / 1_A == Approx(1e0));
+ CHECK(mega * bar / bar == Approx(1e6));
}
SECTION("Formulas") {
const HEPEnergyType E2 = 20_GeV * 2;
- REQUIRE(E2 == 40_GeV);
- REQUIRE(E2 / 1_GeV == Approx(40));
+ CHECK(E2 == 40_GeV);
+ CHECK(E2 / 1_GeV == Approx(40));
const MassType m = 1_kg;
const SpeedType v = 1_m / 1_s;
- REQUIRE(m * v == 1_s * newton);
+ CHECK(m * v == 1_s * newton);
const double lgE = log10(E2 / 1_GeV);
- REQUIRE(lgE == Approx(log10(40.)));
+ CHECK(lgE == Approx(log10(40.)));
const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV;
- REQUIRE(E3 == 180_GeV);
+ CHECK(E3 == 180_GeV);
}
SECTION("Output") {
@@ -96,61 +96,61 @@ TEST_CASE("PhysicalUnits", "[Units]") {
const HEPEnergyType E = 5_eV;
std::stringstream stream;
stream << E;
- REQUIRE(stream.str() == std::string("5 eV"));
+ CHECK(stream.str() == std::string("5 eV"));
}
{
const HEPEnergyType E = 5_EeV;
std::stringstream stream;
stream << E;
- REQUIRE(stream.str() == std::string("5e+18 eV"));
+ CHECK(stream.str() == std::string("5e+18 eV"));
}
}
SECTION("Special") {
const LengthType farAway = std::numeric_limits<double>::infinity() * meter;
- REQUIRE(farAway > 100000_m);
- REQUIRE_FALSE(farAway < 1e19 * meter);
+ CHECK(farAway > 100000_m);
+ CHECK_FALSE(farAway < 1e19 * meter);
}
SECTION("static_pow") {
using namespace corsika::units;
double x = 235.7913;
- REQUIRE(1 == static_pow<0, double>(x));
- REQUIRE(x == static_pow<1, double>(x));
- REQUIRE(x * x == static_pow<2, double>(x));
- REQUIRE(1 / x == static_pow<-1, double>(x));
- REQUIRE(1 / x / x == static_pow<-2, double>(x));
+ CHECK(1 == static_pow<0, double>(x));
+ CHECK(x == static_pow<1, double>(x));
+ CHECK(x * x == static_pow<2, double>(x));
+ CHECK(1 / x == static_pow<-1, double>(x));
+ CHECK(1 / x / x == static_pow<-2, double>(x));
}
SECTION("HEP/SI conversion") {
auto const invEnergy = 1 / 197.326978_MeV; // should be convertible to length or time
- LengthType const length = ConvertHEPToSI<LengthType::dimension_type>(invEnergy);
- REQUIRE((length / 1_fm) == Approx(1));
+ LengthType const length = convert_HEP_to_SI<LengthType::dimension_type>(invEnergy);
+ CHECK((length / 1_fm) == Approx(1));
- TimeType const time = ConvertHEPToSI<TimeType::dimension_type>(invEnergy);
- REQUIRE((time / (1_fm / constants::c)) == Approx(1));
+ TimeType const time = convert_HEP_to_SI<TimeType::dimension_type>(invEnergy);
+ CHECK((time / (1_fm / constants::c)) == Approx(1));
auto const protonMass = 938.272'081'3_MeV; // convertible to mass or SI energy
- MassType protonMassSI = ConvertHEPToSI<MassType::dimension_type>(protonMass);
- REQUIRE((protonMassSI / 1.672'621'898e-27_kg) == Approx(1));
- REQUIRE((protonMassSI / (1.007'276 * constants::u)) == Approx(1));
+ MassType protonMassSI = convert_HEP_to_SI<MassType::dimension_type>(protonMass);
+ CHECK((protonMassSI / 1.672'621'898e-27_kg) == Approx(1));
+ CHECK((protonMassSI / (1.007'276 * constants::u)) == Approx(1));
}
SECTION("SI/HEP conversion") {
- REQUIRE(ConvertSIToHEP(constants::c) == Approx(1));
- REQUIRE(ConvertSIToHEP(constants::hBar) == Approx(1));
+ CHECK(convert_SI_to_HEP(constants::c) == Approx(1));
+ CHECK(convert_SI_to_HEP(constants::hBar) == Approx(1));
{
auto const invLength = 1 / 197.326978_fm; // should be convertible to HEPEnergy
- HEPEnergyType const energy = ConvertSIToHEP(invLength);
- REQUIRE(energy / 1_MeV == Approx(1));
+ HEPEnergyType const energy = convert_SI_to_HEP(invLength);
+ CHECK(energy / 1_MeV == Approx(1));
}
- REQUIRE(ConvertSIToHEP(6.5823e-25_s) * 1_GeV == Approx(1).epsilon(1e-4));
+ CHECK(convert_SI_to_HEP(6.5823e-25_s) * 1_GeV == Approx(1).epsilon(1e-4));
- REQUIRE(ConvertSIToHEP(3.8938e-32 * meter * meter) * 1_GeV * 1_GeV ==
+ CHECK(convert_SI_to_HEP(3.8938e-32 * meter * meter) * 1_GeV * 1_GeV ==
Approx(1).epsilon(1e-4));
}
}