Cylindrical transformations with well-defined phi angle

src/core/observables/CylindricalLBProfileObservable.hpp

@@ -47,15 +47,14 @@ class CylindricalLBProfileObservable : public CylindricalProfileObservable {
         limits[0], limits[1], limits[2], n_bins[0], n_bins[1], n_bins[2],
         sampling_density);
     for (auto &p : sampling_positions) {
-      double theta;
-      Utils::Vector3d rotation_axis;
       auto p_cart = Utils::transform_coordinate_cylinder_to_cartesian(
           p, Utils::Vector3d{{0.0, 0.0, 1.0}});
       // We have to rotate the coordinates since the utils function assumes
       // z-axis symmetry.
-      std::tie(theta, rotation_axis) =
-          Utils::rotation_params(Utils::Vector3d{{0.0, 0.0, 1.0}}, axis);
-      p_cart = Utils::vec_rotate(rotation_axis, theta, p_cart);
+      constexpr auto z_axis = Utils::Vector3d{{0.0, 0.0, 1.0}};
+      auto const theta = angle_between(z_axis, axis);
+      auto const rot_axis = Utils::vector_product(z_axis, axis).normalize();
+      p_cart = Utils::vec_rotate(rot_axis, theta, p_cart);
       p = p_cart + center;
     }
   }

src/utils/include/utils/math/coordinate_transformation.hpp

@@ -19,57 +19,168 @@
 #ifndef UTILS_COORDINATE_TRANSFORMATION_HPP
 #define UTILS_COORDINATE_TRANSFORMATION_HPP
 
+/**
+ * @file
+ * Convert coordinates from the Cartesian system to the cylindrical system.
+ * The transformation functions are provided with three overloads:
+ * - one function for the trivial Cartesian <-> cylindrical transformation
+ * - one function to transform from/to a cylindrical system with custom axis
+ *   (extra @p axis argument, keep in mind the angle phi is under-defined)
+ * - one function to transform from/to an oriented cylindrical system with
+ *   custom axis (extra @p orientation argument, the angle phi is well-defined)
+ */
+
 #include "utils/Vector.hpp"
 #include "utils/constants.hpp"
+#include "utils/math/interval.hpp"
 #include "utils/math/vec_rotate.hpp"
 
+#include <cmath>
+
 namespace Utils {
 
-/** \brief Transform the given 3D position to cylinder coordinates with
- * longitudinal axis aligned with axis parameter.
+/**
+ * @brief Coordinate transformation from Cartesian to cylindrical coordinates.
+ * The origins and z-axis of the coordinate systems co-incide.
+ * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the
+ * original coordinate system.
+ * @param pos    %Vector to transform
  */
 inline Vector3d
-transform_coordinate_cartesian_to_cylinder(const Vector3d &pos,
-                                           const Vector3d &axis) {
+transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) {
+  auto const r = std::sqrt(pos[0] * pos[0] + pos[1] * pos[1]);
+  auto const phi = std::atan2(pos[1], pos[0]);
+  return {r, phi, pos[2]};
+}
+
+/**
+ * @brief Coordinate transformation from Cartesian to cylindrical coordinates
+ * with change of basis. The origins of the coordinate systems co-incide.
+ *
+ * If the parameter @p axis is not equal to <tt>[0, 0, 1]</tt>, the value
+ * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined.
+ * To fully define it, it is necessary to provide an orientation vector
+ * in Cartesian coordinates that will be used as the reference point
+ * (i.e. such that @f$ \phi = 0 @f$).
+ *
+ * If you don't need performance, you can call the overloaded function
+ * @ref transform_coordinate_cartesian_to_cylinder(Vector3d const &, <!--
+ * --> Vector3d const &, Vector3d const &) directly. Otherwise, call this
+ * function once on the orientation vector and cache its @f$ \phi @f$ angle;
+ * this angle can be later passed as argument @p phi0 for a 50-80% speed-up.
+ *
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Angle offset to get @f$ \phi = 0 @f$ (obtained by calling
+ *               this function on the orientation vector instead of @p pos)
+ */
+inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos,
+                                                           Vector3d const &axis,
+                                                           double phi0 = 0.) {
   static auto const z_axis = Vector3d{{0, 0, 1}};
-  double theta;
-  Vector3d rotation_axis;
-  std::tie(theta, rotation_axis) = rotation_params(axis, z_axis);
-  auto const rotated_pos = vec_rotate(rotation_axis, theta, pos);
-  auto const r = std::sqrt(rotated_pos[0] * rotated_pos[0] +
-                           rotated_pos[1] * rotated_pos[1]);
-  auto const phi = std::atan2(rotated_pos[1], rotated_pos[0]);
-  return Vector3d{r, phi, rotated_pos[2]};
+  auto const angle = angle_between(axis, z_axis);
+  auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize();
+  auto const pos_rotated = vec_rotate(rotation_axis, angle, pos);
+  auto pos_cylinder = transform_coordinate_cartesian_to_cylinder(pos_rotated);
+  if (phi0 != 0.) {
+    pos_cylinder[1] = interval(pos_cylinder[1] - phi0, -pi(), pi());
+  }
+  return pos_cylinder;
+}
+
+/**
+ * @brief Coordinate transformation from Cartesian to cylindrical coordinates
+ * with change of basis. The origins of the coordinate systems co-incide.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param orientation   Reference point (in untransformed coordinates) for
+ *                      which @f$ \phi = 0 @f$
+ */
+inline Vector3d transform_coordinate_cartesian_to_cylinder(
+    Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) {
+  auto const phi0 =
+      transform_coordinate_cartesian_to_cylinder(orientation, axis)[1];
+  return transform_coordinate_cartesian_to_cylinder(pos, axis, phi0);
 }
 
 /**
- * @brief Coordinate transformation from cylinder to cartesian coordinates.
+ * @brief Coordinate transformation from cylindrical to Cartesian coordinates.
+ * The origins and z-axis of the coordinate systems co-incide.
+ * The @f$ \phi = 0 @f$ direction corresponds to the x-axis in the
+ * transformed coordinate system.
+ * @param pos    %Vector to transform
  */
 inline Vector3d
-transform_coordinate_cylinder_to_cartesian(Vector3d const &pos,
-                                           Vector3d const &axis) {
-  Vector3d const transformed{
-      {pos[0] * std::cos(pos[1]), pos[0] * std::sin(pos[1]), pos[2]}};
+transform_coordinate_cylinder_to_cartesian(Vector3d const &pos) {
+  auto const &rho = pos[0];
+  auto const &phi = pos[1];
+  auto const &z = pos[2];
+  return {rho * std::cos(phi), rho * std::sin(phi), z};
+}
+
+/**
+ * @brief Coordinate transformation from cylindrical to Cartesian coordinates
+ * with change of basis. The origins of the coordinate systems co-incide.
+ *
+ * If the parameter @p axis is not equal to <tt>[0, 0, 1]</tt>, the value
+ * of the angle @f$ \phi @f$ in cylindrical coordinates is under-defined.
+ * To fully define it, it is necessary to provide an orientation vector
+ * in Cartesian coordinates that will be used as the reference point
+ * (i.e. such that @f$ \phi = 0 @f$).
+ *
+ * If you don't need performance, you can call the overloaded function
+ * @ref transform_coordinate_cylinder_to_cartesian(Vector3d const &, <!--
+ * --> Vector3d const &, Vector3d const &) directly. Otherwise, call function
+ * @ref transform_coordinate_cartesian_to_cylinder(Vector3d const &, <!--
+ * --> Vector3d const &, double) once on the orientation vector and cache
+ * its @f$ \phi @f$ angle; this angle can be later passed as argument
+ * @p phi0 for an 80-100% speed-up.
+ *
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Angle offset to get @f$ \phi = 0 @f$
+ */
+inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos,
+                                                           Vector3d const &axis,
+                                                           double phi0 = 0.) {
   static auto const z_axis = Vector3d{{0, 0, 1}};
-  double theta;
-  Vector3d rotation_axis;
-  std::tie(theta, rotation_axis) = rotation_params(z_axis, axis);
-  auto const rotated_pos = vec_rotate(rotation_axis, theta, transformed);
-  return rotated_pos;
+  auto const angle = angle_between(z_axis, axis);
+  auto const rotation_axis = Utils::vector_product(z_axis, axis).normalize();
+  auto const pos_cyl = Vector3d{{pos[0], pos[1] + phi0, pos[2]}};
+  auto const pos_cart = transform_coordinate_cylinder_to_cartesian(pos_cyl);
+  auto const pos_rot = vec_rotate(rotation_axis, angle, pos_cart);
+  return pos_rot;
 }
 
-/** \brief Transform the given 3D vector to cylinder coordinates with
- * symmetry axis aligned with axis parameter.
+/**
+ * @brief Coordinate transformation from cylindrical to Cartesian coordinates
+ * with change of basis. The origins of the coordinate systems co-incide.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param orientation   Reference point (in Cartesian coordinates) for
+ *                      which @f$ \phi = 0 @f$
+ */
+inline Vector3d transform_coordinate_cylinder_to_cartesian(
+    Vector3d const &pos, Vector3d const &axis, Vector3d const &orientation) {
+  auto const phi0 =
+      transform_coordinate_cartesian_to_cylinder(orientation, axis)[1];
+  return transform_coordinate_cylinder_to_cartesian(pos, axis, phi0);
+}
+
+/**
+ * @brief Vector transformation from Cartesian to cylindrical coordinates.
+ * @param vec    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param pos    Origin of the vector
  */
 inline Vector3d transform_vector_cartesian_to_cylinder(Vector3d const &vec,
                                                        Vector3d const &axis,
                                                        Vector3d const &pos) {
   static auto const z_axis = Vector3d{{0, 0, 1}};
-  double theta;
-  Vector3d rotation_axis;
-  std::tie(theta, rotation_axis) = rotation_params(axis, z_axis);
-  auto const rotated_pos = vec_rotate(rotation_axis, theta, pos);
-  auto const rotated_vec = vec_rotate(rotation_axis, theta, vec);
+  auto const angle = angle_between(axis, z_axis);
+  auto const rotation_axis = Utils::vector_product(axis, z_axis).normalize();
+  auto const rotated_pos = vec_rotate(rotation_axis, angle, pos);
+  auto const rotated_vec = vec_rotate(rotation_axis, angle, vec);
   // v_r = (x * v_x + y * v_y) / sqrt(x^2 + y^2)
   auto const v_r =
       (rotated_pos[0] * rotated_vec[0] + rotated_pos[1] * rotated_vec[1]) /

src/utils/include/utils/math/interval.hpp

@@ -0,0 +1,36 @@
+/*
+ * Copyright (C) 2020 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef ESPRESSO_INTERVAL_HPP
+#define ESPRESSO_INTERVAL_HPP
+
+#include <cmath>
+
+namespace Utils {
+/**
+ * @brief Wrap around the value of @p val in the interval <tt>[low, high]</tt>.
+ */
+template <typename T> T interval(T val, T low, T high) {
+  if (val == low or val == high)
+    return val;
+  auto const span = high - low;
+  return std::fmod(span + std::fmod(val - low, span), span) + low;
+}
+} // namespace Utils
+
+#endif // ESPRESSO_INTERVAL_HPP

src/utils/include/utils/math/vec_rotate.hpp

@@ -23,7 +23,6 @@
 #include "utils/math/sqr.hpp"
 
 #include <cmath>
-#include <tuple>
 
 namespace Utils {
 /**
@@ -50,17 +49,10 @@ inline Vector3d vec_rotate(const Vector3d &axis, double alpha,
 }
 
 /**
- * @brief Determine rotation angle and axis for rotating vec onto target_vec.
- * @param vec Vector to be rotated
- * @param target_vec Target vector
- * @return rotation angle and rotation axis
+ * @brief Determine the angle between two vectors.
  */
-inline std::tuple<double, Vector3d>
-rotation_params(Vector3d const &vec, Vector3d const &target_vec) {
-  auto const theta =
-      std::acos(vec * target_vec / (vec.norm() * target_vec.norm()));
-  auto const rotation_axis = Utils::vector_product(vec, target_vec).normalize();
-  return std::make_tuple(theta, rotation_axis);
+inline double angle_between(Vector3d const &v1, Vector3d const &v2) {
+  return std::acos(v1 * v2 / (v1.norm() * v2.norm()));
 }
 
 } // namespace Utils

src/utils/tests/CMakeLists.txt

@@ -1,6 +1,7 @@
 include(unit_test)
 
 unit_test(NAME abs_test SRC abs_test.cpp DEPENDS EspressoUtils)
+unit_test(NAME interval_test SRC interval_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME Vector_test SRC Vector_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME Factory_test SRC Factory_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME NumeratedContainer_test SRC NumeratedContainer_test.cpp DEPENDS