IAP GITLAB

Commit e3a74d98 authored by Ralf Ulrich's avatar Ralf Ulrich

never expose CoordinateSystem pointer, only shared_ptr

parent b98680ef
......@@ -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
......@@ -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
......@@ -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();
......
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -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>;