From 3c2ca22a9bf9e6ca74d5ba17de3a0e42c09d09ac Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Fri, 10 Aug 2018 18:34:40 +0200
Subject: [PATCH] fixed issue with multiplication with dimensionless result,
 solution requires C++17

---
 CMakeLists.txt                      |  2 +-
 Framework/Geometry/QuantityVector.h | 14 ++++++++++----
 Framework/Geometry/Vector.h         | 20 ++++++++++++++------
 ThirdParty/phys/units/quantity.hpp  |  2 +-
 4 files changed, 26 insertions(+), 12 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index ab861c930..2437c992d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,7 +6,7 @@ project (corsika VERSION 8.0.0 DESCRIPTION "CORSIKA C++ PROJECT " LANGUAGES CXX)
 set (CMAKE_INSTALL_MESSAGE LAZY)
 
 # --std=c++14
-set (CMAKE_CXX_STANDARD 14)
+set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
 
 #add_custom_target (corsika_pre_build)
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 36cf0d567..025469620 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -58,10 +58,16 @@ public:
     template <typename ScalarDim>
     auto operator*(phys::units::quantity<ScalarDim, double> const p) const
     {
-        return QuantityVector<typename phys::units::detail::Product<ScalarDim, dim, double, double>::dimension_type>(eVector * p.magnitude());
-        // TODO: this function does not work if the result is dimensionless, as
-        // dimensionless quantities are "cast" back to plain old double in PhysUnits.
-        // Either change PhysUnits, or cover this case with a template specialization?
+        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());
+        }
     }
     
     auto operator*(double const p) const
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 47b91f7cc..5670aaaf1 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -55,18 +55,18 @@ public:
     }
         
     template <typename dim2>
-    auto parallelProjectionOnto(BaseVector<dim2> const& pVec, CoordinateSystem const& pCS) const
+    auto parallelProjectionOnto(Vector<dim2> const& pVec, CoordinateSystem const& pCS) const
     {
         auto const ourCompVec = getComponents(pCS);
-        auto const otherCompVec = pVec.getComponents(pVec);
+        auto const otherCompVec = pVec.getComponents(pCS);
         auto const& a = ourCompVec.eVector;
         auto const& b = otherCompVec.eVector;
         
-        return Vector<dim>(pCS, (a * b) / b.squaredNorm() * b);
+        return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
     }
     
     template <typename dim2>
-    auto parallelProjectionOnto(BaseVector<dim2> const& pVec)
+    auto parallelProjectionOnto(Vector<dim2> const& pVec) const
     {
         return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
     }
@@ -92,8 +92,16 @@ public:
     template <typename ScalarDim>
     auto operator*(phys::units::quantity<ScalarDim, double> const p) const
     {
-        using res_dim = typename decltype(BaseVector<dim>::qVector * p)::dimension;
-        return Vector<res_dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);        
+        using ResQuantity = decltype(BaseVector<dim>::qVector * p);
+        if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
+        {
+            return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
+        }
+        else
+        {
+            using res_dim = typename decltype(BaseVector<dim>::qVector * p)::dimension;
+            return Vector<res_dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);        
+        }
     }
     
     auto operator*(double const p) const
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index d64e56bca..29ab27f65 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -340,7 +340,7 @@ private:
 
     enum { has_dimension = ! Dims::is_all_zero };
 
-    static_assert( has_dimension, "quantity dimensions must not all be zero" );
+    // static_assert( has_dimension, "quantity dimensions must not all be zero" );
 
 private:
     // friends:
-- 
GitLab