Cylindrical transformations with well-defined phi angle

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

@@ -19,57 +19,125 @@
 #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 phi0 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>
+#include <tuple>
+
 namespace Utils {
 
-/** \brief Transform the given 3D position to cylinder coordinates with
- * longitudinal axis aligned with axis parameter.
+/**
+ * @brief Coordinate transformation from cylindrical to Cartesian coordinates.
+ * @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 cylindrical to Cartesian coordinates
+ * with change of basis.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Reference angle for @f$ \phi = 0 @f$
+ */
+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 rot = rotation_params(axis, z_axis);
+  auto const pos_rotated = vec_rotate(std::get<1>(rot), std::get<0>(rot), 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 cylindrical to Cartesian coordinates
+ * with change of basis.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Reference point at which @f$ \phi = 0 @f$
+ */
+inline Vector3d transform_coordinate_cartesian_to_cylinder(
+    Vector3d const &pos, Vector3d const &axis, Vector3d const &phi0) {
+  auto const phi = transform_coordinate_cartesian_to_cylinder(phi0, axis)[1];
+  return transform_coordinate_cartesian_to_cylinder(pos, axis, phi);
 }
 
 /**
- * @brief Coordinate transformation from cylinder to cartesian coordinates.
+ * @brief Coordinate transformation from cylindrical to Cartesian coordinates.
+ * @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.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Reference angle at which @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 rot = rotation_params(z_axis, axis);
+  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(std::get<1>(rot), std::get<0>(rot), 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.
+ * @param pos    %Vector to transform
+ * @param axis   Longitudinal axis of the cylindrical coordinates
+ * @param phi0   Reference point at which @f$ \phi = 0 @f$
+ */
+inline Vector3d transform_coordinate_cylinder_to_cartesian(
+    Vector3d const &pos, Vector3d const &axis, Vector3d const &phi0) {
+  auto const phi = transform_coordinate_cartesian_to_cylinder(phi0, axis)[1];
+  return transform_coordinate_cylinder_to_cartesian(pos, axis, phi);
+}
+
+/**
+ * @brief Vector transformation from cylindrical to Cartesian 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 rot = rotation_params(axis, z_axis);
+  auto const rotated_pos = vec_rotate(std::get<1>(rot), std::get<0>(rot), pos);
+  auto const rotated_vec = vec_rotate(std::get<1>(rot), std::get<0>(rot), 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/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

src/utils/tests/coordinate_transformation.cpp

@@ -29,6 +29,17 @@
 using Utils::Vector3d;
 
 BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) {
+  constexpr auto eps = 1e-14;
+  auto const pos = Vector3d{{1.0, 3.3, 2.0}};
+  auto const cyl = transform_coordinate_cartesian_to_cylinder(pos);
+  BOOST_CHECK_SMALL(abs(cyl[0] - std::sqrt(pos[0] * pos[0] + pos[1] * pos[1])),
+                    eps);
+  BOOST_CHECK_SMALL(abs(cyl[1] - std::atan2(pos[1], pos[0])), eps);
+  BOOST_CHECK_SMALL(abs(cyl[2] - pos[2]), eps);
+}
+
+BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_test) {
+  constexpr auto eps = 1e-14;
   Vector3d const cart_coord{{1.0, 3.3, 2.0}};
   auto const transformed_x = transform_coordinate_cartesian_to_cylinder(
       cart_coord, Vector3d{{1, 0, 0}});
@@ -51,13 +62,76 @@ BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) {
        std::atan2(cart_coord[1], cart_coord[0]), cart_coord[2]}};
 
   for (int i = 0; i < 3; ++i) {
-    BOOST_CHECK(transformed_x[i] == expected_x[i]);
-    BOOST_CHECK(transformed_y[i] == expected_y[i]);
-    BOOST_CHECK(transformed_z[i] == expected_z[i]);
+    BOOST_CHECK_SMALL(abs(transformed_x[i] - expected_x[i]), eps);
+    BOOST_CHECK_SMALL(abs(transformed_y[i] - expected_y[i]), eps);
+    BOOST_CHECK_SMALL(abs(transformed_z[i] - expected_z[i]), eps);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_with_phi_test) {
+  constexpr auto eps = 1e-14;
+  // tilted orthogonal basis
+  auto const x =
+      (Vector3d{{1, 0, 0}} - (1. / 3.) * Vector3d{{1, 1, 1}}).normalize();
+  auto const y = (Vector3d{{0, 1, -1}}).normalize();
+  auto const z = (Vector3d{{1, 1, 1}}).normalize();
+  // check simple transformation without orientation (phi is random)
+  {
+    auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z);
+    auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z);
+    auto const x_ref = Vector3d{{1.0, x_cyl[1], 0.0}};
+    auto const y_ref = Vector3d{{1.0, y_cyl[1], 0.0}};
+    auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}};
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(x_cyl[i] - x_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(y_cyl[i] - y_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cyl[i] - z_ref[i]), eps);
+    }
+  }
+  // check transformation with orientation (phi is only random for r=0)
+  {
+    auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y);
+    auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y);
+    auto const x_ref = Vector3d{{1.0, -Utils::pi() / 2.0, 0.0}};
+    auto const y_ref = Vector3d{{1.0, 0.0, 0.0}};
+    auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}};
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(x_cyl[i] - x_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(y_cyl[i] - y_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cyl[i] - z_ref[i]), eps);
+    }
+  }
+  // check transformation with orientation for another angle
+  {
+    auto const u = vec_rotate(z, Utils::pi() / 3.0, x);
+    auto const v = vec_rotate(z, Utils::pi() / 3.0, y);
+    auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y);
+    auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y);
+    auto const u_ref = Vector3d{{1.0, Utils::pi() * (1. / 3. - 1. / 2.), 0.0}};
+    auto const v_ref = Vector3d{{1.0, Utils::pi() / 3.0, 0.0}};
+    auto const z_ref = Vector3d{{0.0, z_cyl[1], 1.0}};
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(u_cyl[i] - u_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(v_cyl[i] - v_ref[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cyl[i] - z_ref[i]), eps);
+    }
   }
 }
 
 BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) {
+  constexpr auto eps = 1e-14;
+  auto const cyl = Vector3d{{1.0, Utils::pi() / 4, 2.0}};
+  auto const pos = transform_coordinate_cylinder_to_cartesian(cyl);
+  BOOST_CHECK_SMALL(abs(pos[0] - std::sqrt(2) / 2), eps);
+  BOOST_CHECK_SMALL(abs(pos[1] - std::sqrt(2) / 2), eps);
+  BOOST_CHECK_SMALL(abs(pos[2] - cyl[2]), eps);
+}
+
+BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_test) {
+  constexpr auto eps = 1e-14;
   Vector3d const cylinder_coord{{1.2, 3.123, 42.0}};
   auto const transformed_x = transform_coordinate_cylinder_to_cartesian(
       cylinder_coord, Vector3d{{1, 0, 0}});
@@ -80,8 +154,61 @@ BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) {
       {cylinder_coord[0] * std::cos(cylinder_coord[1]),
        cylinder_coord[0] * std::sin(cylinder_coord[1]), cylinder_coord[2]}};
   for (int i = 0; i < 3; ++i) {
-    BOOST_CHECK(transformed_x[i] == expected_x[i]);
-    BOOST_CHECK(transformed_y[i] == expected_y[i]);
-    BOOST_CHECK(transformed_z[i] == expected_z[i]);
+    BOOST_CHECK_SMALL(abs(transformed_x[i] - expected_x[i]), eps);
+    BOOST_CHECK_SMALL(abs(transformed_y[i] - expected_y[i]), eps);
+    BOOST_CHECK_SMALL(abs(transformed_z[i] - expected_z[i]), eps);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_with_phi_test) {
+  constexpr auto eps = 1e-14;
+  // tilted orthogonal basis
+  auto const x =
+      (Vector3d{{1, 0, 0}} - (1. / 3.) * Vector3d{{1, 1, 1}}).normalize();
+  auto const y = (Vector3d{{0, 1, -1}}).normalize();
+  auto const z = (Vector3d{{1, 1, 1}}).normalize();
+  // check simple transformation without orientation
+  {
+    auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z);
+    auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z);
+    auto const x_cart = transform_coordinate_cylinder_to_cartesian(x_cyl, z);
+    auto const y_cart = transform_coordinate_cylinder_to_cartesian(y_cyl, z);
+    auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z);
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(x_cart[i] - x[i]), eps);
+      BOOST_CHECK_SMALL(abs(y_cart[i] - y[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cart[i] - z[i]), eps);
+    }
+  }
+  // check transformation with orientation
+  {
+    auto const x_cyl = transform_coordinate_cartesian_to_cylinder(x, z, y);
+    auto const y_cyl = transform_coordinate_cartesian_to_cylinder(y, z, y);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y);
+    auto const x_cart = transform_coordinate_cylinder_to_cartesian(x_cyl, z, y);
+    auto const y_cart = transform_coordinate_cylinder_to_cartesian(y_cyl, z, y);
+    auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z, y);
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(x_cart[i] - x[i]), eps);
+      BOOST_CHECK_SMALL(abs(y_cart[i] - y[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cart[i] - z[i]), eps);
+    }
+  }
+  // check transformation with orientation for another angle
+  {
+    auto const u = vec_rotate(z, Utils::pi() / 3.0, x);
+    auto const v = vec_rotate(z, Utils::pi() / 3.0, y);
+    auto const u_cyl = transform_coordinate_cartesian_to_cylinder(u, z, y);
+    auto const v_cyl = transform_coordinate_cartesian_to_cylinder(v, z, y);
+    auto const z_cyl = transform_coordinate_cartesian_to_cylinder(z, z, y);
+    auto const u_cart = transform_coordinate_cylinder_to_cartesian(u_cyl, z, y);
+    auto const v_cart = transform_coordinate_cylinder_to_cartesian(v_cyl, z, y);
+    auto const z_cart = transform_coordinate_cylinder_to_cartesian(z_cyl, z, y);
+    for (int i = 0; i < 3; ++i) {
+      BOOST_CHECK_SMALL(abs(u_cart[i] - u[i]), eps);
+      BOOST_CHECK_SMALL(abs(v_cart[i] - v[i]), eps);
+      BOOST_CHECK_SMALL(abs(z_cart[i] - z[i]), eps);
+    }
   }
 }