fixed issue with multiplication with dimensionless result, solution

requires C++17

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)

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

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

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: