From ada19b6a826968c8996a38fac3e1f1fecacb0794 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 26 Nov 2020 18:47:43 +0100 Subject: [PATCH 1/9] utils: Add modulo function for radians --- src/utils/include/utils/math/interval.hpp | 36 ++++++++++++++++++++ src/utils/tests/CMakeLists.txt | 1 + src/utils/tests/interval_test.cpp | 40 +++++++++++++++++++++++ 3 files changed, 77 insertions(+) create mode 100644 src/utils/include/utils/math/interval.hpp create mode 100644 src/utils/tests/interval_test.cpp diff --git a/src/utils/include/utils/math/interval.hpp b/src/utils/include/utils/math/interval.hpp new file mode 100644 index 00000000000..b2ddf33ee58 --- /dev/null +++ b/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 . + */ +#ifndef ESPRESSO_INTERVAL_HPP +#define ESPRESSO_INTERVAL_HPP + +#include + +namespace Utils { +/** + * @brief Wrap around the value of @p val in the interval [low, high]. + */ +template 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 diff --git a/src/utils/tests/CMakeLists.txt b/src/utils/tests/CMakeLists.txt index 1bb9474a482..87caee85934 100644 --- a/src/utils/tests/CMakeLists.txt +++ b/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 diff --git a/src/utils/tests/interval_test.cpp b/src/utils/tests/interval_test.cpp new file mode 100644 index 00000000000..c1433e6f66e --- /dev/null +++ b/src/utils/tests/interval_test.cpp @@ -0,0 +1,40 @@ +/* + * Copyright (C) 2020 The ESPResSo project + * + * 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 . + */ + +#define BOOST_TEST_MODULE interval test +#define BOOST_TEST_DYN_LINK +#include + +#include +#include + +#include + +using Utils::Vector3d; + +BOOST_AUTO_TEST_CASE(interval_test) { + constexpr auto eps = 1e-14; + constexpr auto pi = Utils::pi(); + BOOST_CHECK_EQUAL(Utils::interval(-pi, -pi, pi), -pi); + BOOST_CHECK_EQUAL(Utils::interval(pi, -pi, pi), pi); + for (int i = -24; i < 24; ++i) { + double const theta = i * pi / 6; + double const wrapped = Utils::interval(theta, -pi, pi); + double const ref_val = std::atan2(std::sin(theta), std::cos(theta)); + BOOST_CHECK_SMALL(abs(fmod(wrapped - ref_val, 2 * pi)), eps); + } +} From 9c55cab973e1f7df7ceaa95ef0b0d8421f8e1f27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 26 Nov 2020 18:53:26 +0100 Subject: [PATCH 2/9] tests: Add tolerance for floating point operations --- src/utils/tests/coordinate_transformation.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index 714a11700e1..9121302d5d5 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -29,6 +29,7 @@ using Utils::Vector3d; BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_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 +52,14 @@ 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(cylinder_to_cartesian_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 +82,8 @@ 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); } } From 809009d73b6e0a4a6f940b722150355951796533 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 26 Nov 2020 20:04:41 +0100 Subject: [PATCH 3/9] utils: Refactor coordinate transformation code Provides 3 overloads to transform between Cartesian and: - the cylindrical system - a cylindrical system with custom axis (+axis parameter) - an oriented cylindrical system with custom axis (+phi0 parameter) --- .../utils/math/coordinate_transformation.hpp | 126 ++++++++++++++---- src/utils/tests/coordinate_transformation.cpp | 125 +++++++++++++++++ 2 files changed, 222 insertions(+), 29 deletions(-) diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index 5316016e2de..d80e261c525 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/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 +#include + 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) { + 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]) / diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index 9121302d5d5..ba9b320c5b0 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -29,6 +29,16 @@ 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( @@ -58,7 +68,69 @@ BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_test) { } } +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( @@ -87,3 +159,56 @@ BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) { 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); + } + } +} From b849f904a4f2202e531adafdddf55568a2287d3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 27 Nov 2020 15:49:13 +0100 Subject: [PATCH 4/9] utils: Fix regressions --- src/utils/include/utils/math/coordinate_transformation.hpp | 2 +- src/utils/tests/interval_test.cpp | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index d80e261c525..7b5b5c8c50c 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -65,7 +65,7 @@ inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos, 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) { + if (phi0 != 0.) { pos_cylinder[1] = interval(pos_cylinder[1] - phi0, -pi(), pi()); } return pos_cylinder; diff --git a/src/utils/tests/interval_test.cpp b/src/utils/tests/interval_test.cpp index c1433e6f66e..8f619d3cd10 100644 --- a/src/utils/tests/interval_test.cpp +++ b/src/utils/tests/interval_test.cpp @@ -24,8 +24,6 @@ #include -using Utils::Vector3d; - BOOST_AUTO_TEST_CASE(interval_test) { constexpr auto eps = 1e-14; constexpr auto pi = Utils::pi(); @@ -35,6 +33,6 @@ BOOST_AUTO_TEST_CASE(interval_test) { double const theta = i * pi / 6; double const wrapped = Utils::interval(theta, -pi, pi); double const ref_val = std::atan2(std::sin(theta), std::cos(theta)); - BOOST_CHECK_SMALL(abs(fmod(wrapped - ref_val, 2 * pi)), eps); + BOOST_CHECK_SMALL(std::abs(std::fmod(wrapped - ref_val, 2 * pi)), eps); } } From 91a7d702533c93f74d1c464a09e112f964ed0ce4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 30 Nov 2020 17:02:08 +0100 Subject: [PATCH 5/9] utils: Simplify transformation code Remove the Utils::rotation_params() function in favor of a more basic utility function: Utils::angle_between(). This makes clear at the call sites what is actually calculated. --- .../CylindricalLBProfileObservable.hpp | 9 ++-- .../utils/math/coordinate_transformation.hpp | 44 +++++++++++-------- src/utils/include/utils/math/vec_rotate.hpp | 14 ++---- src/utils/tests/vec_rotate_test.cpp | 16 +++---- 4 files changed, 38 insertions(+), 45 deletions(-) diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp index 5edd4b7ca7c..e2c27d1cf33 100644 --- a/src/core/observables/CylindricalLBProfileObservable.hpp +++ b/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; } } diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index 7b5b5c8c50c..9cb4cccd1ff 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -36,7 +36,6 @@ #include "utils/math/vec_rotate.hpp" #include -#include namespace Utils { @@ -56,14 +55,16 @@ transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) { * 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$ + * @param phi0 Value for which @f$ \phi = 0 @f$, is obtained by the + * overloaded function that takes an orientation vector. */ 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}}; - auto const rot = rotation_params(axis, z_axis); - auto const pos_rotated = vec_rotate(std::get<1>(rot), std::get<0>(rot), pos); + 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()); @@ -76,12 +77,14 @@ inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos, * 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$ + * @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 &phi0) { - auto const phi = transform_coordinate_cartesian_to_cylinder(phi0, axis)[1]; - return transform_coordinate_cartesian_to_cylinder(pos, axis, phi); + 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); } /** @@ -100,16 +103,18 @@ transform_coordinate_cylinder_to_cartesian(Vector3d const &pos) { * @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$ + * @param phi0 Value for which @f$ \phi = 0 @f$, is obtained by the + * overloaded function that takes an orientation vector. */ 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}}; - auto const rot = rotation_params(z_axis, axis); + 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(std::get<1>(rot), std::get<0>(rot), pos_cart); + auto const pos_rot = vec_rotate(rotation_axis, angle, pos_cart); return pos_rot; } @@ -117,12 +122,14 @@ inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, * @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$ + * @param orientation Reference point (in untransformed coordinates) for + * 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); + 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); } /** @@ -135,9 +142,10 @@ 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}}; - 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); + 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]) / diff --git a/src/utils/include/utils/math/vec_rotate.hpp b/src/utils/include/utils/math/vec_rotate.hpp index 4dbc624f168..559f297766f 100644 --- a/src/utils/include/utils/math/vec_rotate.hpp +++ b/src/utils/include/utils/math/vec_rotate.hpp @@ -23,7 +23,6 @@ #include "utils/math/sqr.hpp" #include -#include 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 -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 diff --git a/src/utils/tests/vec_rotate_test.cpp b/src/utils/tests/vec_rotate_test.cpp index 18b9b14ef45..f1d5727a2cd 100644 --- a/src/utils/tests/vec_rotate_test.cpp +++ b/src/utils/tests/vec_rotate_test.cpp @@ -22,11 +22,9 @@ #include #include #include -using Utils::vec_rotate; #include #include -#include BOOST_AUTO_TEST_CASE(rotation) { using std::cos; @@ -43,20 +41,16 @@ BOOST_AUTO_TEST_CASE(rotation) { auto const expected = cos(t) * v + sin(t) * vector_product(k, v) + (1. - cos(t)) * (k * v) * k; - auto const is = vec_rotate(k, t, v); + auto const is = Utils::vec_rotate(k, t, v); auto const rel_diff = (expected - is).norm() / expected.norm(); BOOST_CHECK(rel_diff < std::numeric_limits::epsilon()); } -BOOST_AUTO_TEST_CASE(rotation_params) { - Utils::Vector3d v1 = {1.0, 0.0, 0.0}; - Utils::Vector3d v2 = {1.0, 1.0, 0.0}; +BOOST_AUTO_TEST_CASE(angle_between) { + Utils::Vector3d const v1 = {1.0, 0.0, 0.0}; + Utils::Vector3d const v2 = {1.0, 1.0, 0.0}; - double angle; - Utils::Vector3d rotation_axis; - std::tie(angle, rotation_axis) = Utils::rotation_params(v1, v2); + auto const angle = Utils::angle_between(v1, v2); BOOST_CHECK_CLOSE(angle, Utils::pi() / 4.0, 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v1), 1e-7); - BOOST_CHECK_SMALL((rotation_axis * v2), 1e-7); } From 692a5d0cd8b6e3c8cc9e0fb6522f896a70b0c95d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 2 Dec 2020 12:46:28 +0100 Subject: [PATCH 6/9] utils: Fix doxygen documentation --- .../include/utils/math/coordinate_transformation.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index 9cb4cccd1ff..5fcd30632ee 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -40,7 +40,7 @@ namespace Utils { /** - * @brief Coordinate transformation from cylindrical to Cartesian coordinates. + * @brief Coordinate transformation from Cartesian to cylindrical coordinates. * @param pos %Vector to transform */ inline Vector3d @@ -51,7 +51,7 @@ transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) { } /** - * @brief Coordinate transformation from cylindrical to Cartesian coordinates + * @brief Coordinate transformation from Cartesian to cylindrical coordinates * with change of basis. * @param pos %Vector to transform * @param axis Longitudinal axis of the cylindrical coordinates @@ -73,7 +73,7 @@ inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos, } /** - * @brief Coordinate transformation from cylindrical to Cartesian coordinates + * @brief Coordinate transformation from Cartesian to cylindrical coordinates * with change of basis. * @param pos %Vector to transform * @param axis Longitudinal axis of the cylindrical coordinates @@ -122,7 +122,7 @@ inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, * @brief Coordinate transformation from cylindrical to Cartesian coordinates. * @param pos %Vector to transform * @param axis Longitudinal axis of the cylindrical coordinates - * @param orientation Reference point (in untransformed coordinates) for + * @param orientation Reference point (in Cartesian coordinates) for * which @f$ \phi = 0 @f$ */ inline Vector3d transform_coordinate_cylinder_to_cartesian( @@ -133,7 +133,7 @@ inline Vector3d transform_coordinate_cylinder_to_cartesian( } /** - * @brief Vector transformation from cylindrical to Cartesian coordinates. + * @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 From 09d6dc738ec6efa2d668e9232089d54fdb634ec9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 2 Dec 2020 12:48:03 +0100 Subject: [PATCH 7/9] tests: Remove duplicated checks --- src/utils/tests/coordinate_transformation.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index ba9b320c5b0..815b5f32fbc 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -109,14 +109,11 @@ BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_with_phi_test) { 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); } } } @@ -201,14 +198,11 @@ BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_with_axis_with_phi_test) { 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); } } } From dfcd1f4455661a7d4f6c041e75b26d6f91d58ef0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 7 Dec 2020 17:51:40 +0100 Subject: [PATCH 8/9] utils: Improve documentation --- .../utils/math/coordinate_transformation.hpp | 53 +++++++++++++++---- 1 file changed, 44 insertions(+), 9 deletions(-) diff --git a/src/utils/include/utils/math/coordinate_transformation.hpp b/src/utils/include/utils/math/coordinate_transformation.hpp index 5fcd30632ee..09d17ff97c7 100644 --- a/src/utils/include/utils/math/coordinate_transformation.hpp +++ b/src/utils/include/utils/math/coordinate_transformation.hpp @@ -27,7 +27,7 @@ * - 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) + * custom axis (extra @p orientation argument, the angle phi is well-defined) */ #include "utils/Vector.hpp" @@ -41,6 +41,9 @@ namespace Utils { /** * @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 @@ -52,11 +55,24 @@ transform_coordinate_cartesian_to_cylinder(Vector3d const &pos) { /** * @brief Coordinate transformation from Cartesian to cylindrical coordinates - * with change of basis. + * with change of basis. The origins of the coordinate systems co-incide. + * + * If the parameter @p axis is not equal to [0, 0, 1], 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 Value for which @f$ \phi = 0 @f$, is obtained by the - * overloaded function that takes an orientation vector. + * @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, @@ -74,7 +90,7 @@ inline Vector3d transform_coordinate_cartesian_to_cylinder(Vector3d const &pos, /** * @brief Coordinate transformation from Cartesian to cylindrical coordinates - * with change of basis. + * 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 @@ -89,6 +105,9 @@ inline Vector3d transform_coordinate_cartesian_to_cylinder( /** * @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 @@ -100,11 +119,26 @@ transform_coordinate_cylinder_to_cartesian(Vector3d const &pos) { } /** - * @brief Coordinate transformation from cylindrical to Cartesian coordinates. + * @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 [0, 0, 1], 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 Value for which @f$ \phi = 0 @f$, is obtained by the - * overloaded function that takes an orientation vector. + * @param phi0 Angle offset to get @f$ \phi = 0 @f$ */ inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, Vector3d const &axis, @@ -119,7 +153,8 @@ inline Vector3d transform_coordinate_cylinder_to_cartesian(Vector3d const &pos, } /** - * @brief Coordinate transformation from cylindrical to Cartesian coordinates. + * @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 From 51830528fecb630810e5b17b097bb3c2384d9553 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 7 Dec 2020 20:11:15 +0100 Subject: [PATCH 9/9] tests: Add cylindrical transformation test case --- src/utils/tests/coordinate_transformation.cpp | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/utils/tests/coordinate_transformation.cpp b/src/utils/tests/coordinate_transformation.cpp index 815b5f32fbc..34327232123 100644 --- a/src/utils/tests/coordinate_transformation.cpp +++ b/src/utils/tests/coordinate_transformation.cpp @@ -25,6 +25,7 @@ #include #include +#include using Utils::Vector3d; @@ -116,6 +117,29 @@ BOOST_AUTO_TEST_CASE(cartesian_to_cylinder_with_axis_with_phi_test) { BOOST_CHECK_SMALL(abs(v_cyl[i] - v_ref[i]), eps); } } + // check transformation of random vectors + { + std::subtract_with_carry_engine rng(2); + auto const r_uniform = [&rng]() { + return static_cast(rng() - rng.min()) / (rng.max() - rng.min()); + }; + for (int trial = 0; trial < 100; ++trial) { + Vector3d const v1{r_uniform(), r_uniform(), r_uniform()}; + Vector3d const v2{r_uniform(), r_uniform(), r_uniform()}; + auto const a = Utils::vector_product(v1, v2) / v1.norm() / v2.norm(); + auto const v1_v1 = transform_coordinate_cartesian_to_cylinder(v1, a, v1); + auto const v2_v1 = transform_coordinate_cartesian_to_cylinder(v2, a, v1); + auto const v1_v2 = transform_coordinate_cartesian_to_cylinder(v1, a, v2); + Vector3d const v1_v1_ref{v1.norm(), 0.0, 0.0}; + Vector3d const v2_v1_ref{v2.norm(), Utils::angle_between(v1, v2), 0.0}; + Vector3d const v1_v2_ref{v1.norm(), -Utils::angle_between(v1, v2), 0.0}; + for (int i = 0; i < 3; ++i) { + BOOST_CHECK_SMALL(abs(v1_v1[i] - v1_v1_ref[i]), eps); + BOOST_CHECK_SMALL(abs(v2_v1[i] - v2_v1_ref[i]), eps); + BOOST_CHECK_SMALL(abs(v1_v2[i] - v1_v2_ref[i]), eps); + } + } + } } BOOST_AUTO_TEST_CASE(cylinder_to_cartesian_test) {